Свой 3D движок Python Пересечение отрезка и прямоугольника

Не дождался ответа, поэтому несколько часов решал уравнения и все получилось. вот код для этой функции. Возможно можно как-то еще до оптимизировать numpy

import numba as nb
import numpy as np
import time
import math
import random


@nb.njit(cache=True, fastmath=True)
def _3points_to_eq_of_plane(_3v3):
    v3_1 = _3v3[0]
    v3_2 = _3v3[1]
    v3_3 = _3v3[2]
    x1 = v3_1[0]
    y1 = v3_1[1]
    z1 = v3_1[2]
    x2 = v3_2[0]
    y2 = v3_2[1]
    z2 = v3_2[2]
    x3 = v3_3[0]
    y3 = v3_3[1]
    z3 = v3_3[2]
    vector1 = np.array([x2 - x1, y2 - y1, z2 - z1])
    vector2 = np.array([x3 - x1, y3 - y1, z3 - z1])

    cross_product = [vector1[1] * vector2[2] - vector1[2] * vector2[1],
                     -1 * (vector1[0] * vector2[2] - vector1[2] * vector2[0]),
                     vector1[0] * vector2[1] - vector1[1] * vector2[0]]

    a = cross_product[0]
    b = cross_product[1]
    c = cross_product[2]
    d = - (cross_product[0] * x1 + cross_product[1] * y1 + cross_product[2] * z1)
    return np.array([a, b, c, d])


@nb.njit(cache=True, fastmath=True)
def xyz_of_intersection_plane_and_line(a, b, c, d, x1, x2, y1, y2, z1, z2):
    x21 = (x2 - x1)
    y21 = (y2 - y1)
    z21 = (z2 - z1)
    zxy21_c = (z21 + x21) * c / y21
    y = (y1 * zxy21_c - z1 * c - d - x1 * a) / (zxy21_c + b)
    x = (y - y1) * x21 / y21 + x1
    z = (-d - a * x - b * y) / c
    return np.array([x, y, z])

@nb.njit(cache=True, fastmath=True)
def point_to_pix_on_rect(view_rect, v3, size_x, size_y):
    xsx = (view_rect[1][0] - view_rect[0][0]) / (size_x - 1)
    xsy = (view_rect[1][1] - view_rect[0][1]) / (size_x - 1)
    xsz = (view_rect[1][2] - view_rect[0][2]) / (size_x - 1)
    ysx = (view_rect[3][0] - view_rect[0][0]) / (size_y - 1)
    ysy = (view_rect[3][1] - view_rect[0][1]) / (size_y - 1)
    ysz = (view_rect[3][2] - view_rect[0][2]) / (size_y - 1)
    c1x = view_rect[0][0]
    c1y = view_rect[0][1]
    c1z = view_rect[0][2]
    x = v3[0]
    y = v3[1]
    z = v3[2]
    # Здесь система из 3 линейный уравнений с 2 неизвестными
    # Выбираем лучшие - самая большая разница координат для меньшей ошибки
    dx = abs(view_rect[2][0] - view_rect[0][0])
    dy = abs(view_rect[2][1] - view_rect[0][1])
    dz = abs(view_rect[2][2] - view_rect[0][2])
    if ysx == 0 and ysz == 0:
        i = (x - c1x) / xsx
        j = (y - c1y) / ysy
    elif ysy == 0 and ysz == 0:
        i = (x - c1x) / ysx
        j = (y - c1y) / xsx
    elif ysy == 0 and ysx == 0:
        i = (x - c1x) / ysz
        j = (y - c1y) / ysz
    elif min(dx, dy, dz) == dx:
        # i * xsx + j * ysx = x - c1x
        # i * xsy + j * ysy = y - c1y
        print(ysx)
        i = (y - c1y + x * ysy / ysx + c1x * ysy / ysx) / (xsy - xsx * ysy / ysx)
        j = (x - c1x - i * xsx) / ysx

    elif min(dx, dy, dz) == dy:
        # i * xsx + j * ysx = x - c1x
        # i * xsz + j * ysz = z - c1z

        i = (z - c1z + x * ysz / ysx + c1x * ysz / ysx) / (xsz - xsx * ysz / ysx)
        j = (x - c1x - i * xsx) / ysx

    elif min(dx, dy, dz) == dz:
        # i * xsz + j * ysz = z - c1z
        # i * xsy + j * ysy = y - c1y

        i = (y - c1y + z * ysy / ysz + c1z * ysy / ysz) / (xsy - xsz * ysy / ysz)
        j = (y - c1y - i * xsy) / ysy
    return np.array([i, j])

@nb.njit(cache=True, fastmath=True)
def rect_view_triangle_size_to_pix(rect, view, triangle, size_x, size_y):
    points_from_rect = np.zeros((3, 3))
    points_from_rect[0] = rect[0]
    points_from_rect[1] = rect[1]
    points_from_rect[2] = rect[2]
    plane = _3points_to_eq_of_plane(points_from_rect)
    a = plane[0]
    b = plane[1]
    c = plane[2]
    d = plane[3]
    v31 = triangle[0]
    v32 = triangle[1]
    v33 = triangle[2]
    x = view[0]
    y = view[1]
    z = view[2]
    x1 = v31[0]
    y1 = v31[1]
    z1 = v31[2]
    x2 = v32[0]
    y2 = v32[1]
    z2 = v32[2]
    x3 = v33[0]
    y3 = v33[1]
    z3 = v33[2]
    v31 = xyz_of_intersection_plane_and_line(a, b, c, d, x, x1, y, y1, z, z1)
    v32 = xyz_of_intersection_plane_and_line(a, b, c, d, x, x2, y, y2, z, z2)
    v33 = xyz_of_intersection_plane_and_line(a, b, c, d, x, x3, y, y3, z, z3)
    v1 = point_to_pix_on_rect(rect, v31, size_x, size_y)
    v2 = point_to_pix_on_rect(rect, v32, size_x, size_y)
    v3 = point_to_pix_on_rect(rect, v33, size_x, size_y)
    ret = np.zeros((3, 2))
    ret[0] = v1
    ret[1] = v2
    ret[2] = v3
    return ret


@nb.njit(cache=True, fastmath=True, parallel=True)
def get_triangles_on_screen(triangles, rect, size_x, size_y, x, y, z):
    view = np.zeros(3)
    view[0] = x
    view[1] = y
    view[2] = z
    lt = len(triangles)
    result = np.zeros((lt, 3, 2))
    for i in nb.prange(lt):
        result[i] = rect_view_triangle_size_to_pix(rect, view, triangles[i], size_x, size_y)
    return result

Остальное не покажу. картинка рендерится пол секунды 470000 полигонов. triangles - массив треугольников rect - прямоугольник size_x - ширина окна size_y - высота окна x, y, z - положение камеры