Python光线追踪：3D空间多边形与射线的首次相交计算

我们可以通过向量化的Möller-Trumbore射线-三角形相交算法实现需求，仅依赖numpy即可完成所有计算，不需要额外引入其他库。思路如下：

• 所有输入的4顶点多边形（四边形）先拆分为2个三角形，排除用于闭合的第5个冗余顶点
• 对每条射线（固定x、z坐标，沿y轴方向发射），批量计算和所有三角形的交点
• 过滤有效交点后取y坐标最小的即为最先命中的交点值

完整代码实现 #

import numpy as np
import matplotlib.pyplot as plt

def ray_triangle_intersect(ray_orig, ray_dir, v0, v1, v2):
    """
    向量化Möller-Trumbore算法，计算射线与三角形的交点
    参数：
        ray_orig: 射线原点，形状(N, 3) 或 (3,)
        ray_dir: 射线方向，形状(N, 3) 或 (3,)
        v0, v1, v2: 三角形三个顶点，形状(M, 3) 或 (3,)
    返回：
        t: 射线方向上的交点距离，无交点则为inf
    """
    eps = 1e-8
    # 扩展维度适配批量计算
    ray_orig = np.atleast_2d(ray_orig)
    ray_dir = np.atleast_2d(ray_dir)
    v0 = np.atleast_2d(v0)
    v1 = np.atleast_2d(v1)
    v2 = np.atleast_2d(v2)
    
    edge1 = v1 - v0
    edge2 = v2 - v0
    h = np.cross(ray_dir[:, None, :], edge2[None, :, :])
    a = np.einsum('ijk,ijk->ij', edge1[None, :, :], h)
    
    # 射线与三角形平行的情况
    mask = np.abs(a) < eps
    a[mask] = eps
    
    f = 1.0 / a
    s = ray_orig[:, None, :] - v0[None, :, :]
    u = f * np.einsum('ijk,ijk->ij', s, h)
    
    # u不在[0,1]区间，无交点
    mask_u = (u < 0.0) | (u > 1.0)
    q = np.cross(s, edge1[None, :, :])
    v = f * np.einsum('ijk,ijk->ij', ray_dir[:, None, :], q)
    
    # v不在[0,1]区间，或者u+v>1，无交点
    mask_v = (v < 0.0) | (u + v > 1.0)
    t = f * np.einsum('ijk,ijk->ij', edge2[None, :, :], q)
    
    # 交点在射线反方向的情况
    mask_t = t < eps
    
    # 所有无效情况设置为inf
    t[mask | mask_u | mask_v | mask_t] = np.inf
    return t

def compute_ray_hits(quads, x_grid, z_grid, ray_dir=np.array([0, 1, 0])):
    """
    计算x-z网格上所有沿y方向发射的射线的最近交点y值
    参数：
        quads: 输入四边形列表，每个四边形形状为(4,3)（已经去掉第5个闭合顶点）
        x_grid: x坐标网格，1维数组
        z_grid: z坐标网格，1维数组
        ray_dir: 射线方向，默认沿y轴正方向
    返回：
        hit_y: 最近交点y值网格，形状为(len(z_grid), len(x_grid))，无交点则为inf
    """
    # 把所有四边形拆成三角形
    tris = []
    for quad in quads:
        # 四边形v0,v1,v2,v3拆为v0,v1,v2和v0,v2,v3两个三角形
        tris.append(quad[[0,1,2]])
        tris.append(quad[[0,2,3]])
    tris = np.array(tris)
    v0, v1, v2 = tris[:,0], tris[:,1], tris[:,2]
    
    # 生成所有射线原点
    X, Z = np.meshgrid(x_grid, z_grid)
    ray_orig = np.stack([X.ravel(), np.full_like(X.ravel(), -100), Z.ravel()], axis=1)
    
    # 批量计算所有射线和所有三角形的交点t值
    t_vals = ray_triangle_intersect(ray_orig, ray_dir, v0, v1, v2)
    
    # 取每个射线的最小t值，计算对应y坐标
    min_t = np.min(t_vals, axis=1)
    hit_y = ray_orig[:,1] + min_t * ray_dir[1]
    hit_y = hit_y.reshape(X.shape)
    return hit_y

# ------------ 测试示例 ------------
# 示例数据预处理，去掉第5个闭合顶点
square1 = np.array([
    [0, 0, 0],
    [1, 0, 0],
    [1, 0.5, 1],
    [0, 0.5, 1],
    [0, 0, 0]])[:4]
square2 = square1 + 0.5 * np.ones_like(square1)

# 测试单条射线：x=0.75, z=0.75
test_ray_orig = np.array([0.75, -100, 0.75])
test_ray_dir = np.array([0,1,0])
# 把两个四边形拆成三角形
all_tris = np.concatenate([
    square1[[0,1,2]][None,:], square1[[0,2,3]][None,:],
    square2[[0,1,2]][None,:], square2[[0,2,3]][None,:]
])
v0, v1, v2 = all_tris[:,0], all_tris[:,1], all_tris[:,2]
t = ray_triangle_intersect(test_ray_orig, test_ray_dir, v0, v1, v2)
min_t = np.min(t)
hit_y = test_ray_orig[1] + min_t * test_ray_dir[1]
print(f"x=0.75、z=0.75处射线最先命中的y坐标为：{hit_y:.4f}")

# 批量生成x-z网格的命中结果示例
x_grid = np.linspace(-1, 2, 100)
z_grid = np.linspace(-1, 2, 100)
hit_y_grid = compute_ray_hits([square1, square2], x_grid, z_grid)

# 可视化结果
plt.figure(figsize=(5,4))
plt.pcolormesh(x_grid, z_grid, hit_y_grid, shading='auto')
plt.colorbar(label='命中y坐标')
plt.xlabel('x')
plt.ylabel('z')
plt.title('射线最近命中y坐标分布')
plt.gca().set_aspect('equal')
plt.show()

输出说明 #

示例测试中x=0.75、z=0.75处的射线会先命中square1，输出的y坐标约为0.375，符合几何计算结果。批量计算后得到的网格可以直接用于生成深度图或渲染图像。

内容的提问来源于stack exchange，提问作者John Smith