Möller–Trumbore射线三角面相交算法（The Möller–Trumbore ray-triangle intersection algorithm) 是一种计算机图形学中经典的算法，用来计算射线和三维空间中三角形的交点。该方法的优点是方法速度快，存储空间少。
输入输出输入： 给定3维空间中的三个点构成一个三角面片，再给定一个射线起点和方向向量；
输出： 求射线和三角平面在三维空间中的交点。
结果展示
输入： 三角形点[[0, 0, 1], [0, 1, 0], [1, 0, 0]]
方向[2,2,2]
输出： 蓝色交点（0.333，0.333，0.333）
计算权重系数
import numpy as np
import pyvista as pv
def ray_triangle_intersection(ray_start, ray_vec, triangle):
"""Moeller–Trumbore intersection algorithm.
Parameters
----------
ray_start[n, 3] : np.ndarray
Length three numpy array representing start of point.
ray_vec : np.ndarray
Direction of the ray.
triangle : np.ndarray
``3 x 3`` numpy array containing the three vertices of a
triangle.
Returns
-------
bool
``True`` when there is an intersection.
tuple
Length three tuple containing the distance ``t``, and the
intersection in unit triangle ``u``, ``v`` coordinates. When
there is no intersection, these values will be:
``[np.nan, np.nan, np.nan]``
"""
# define a null intersection
null_inter = np.array([np.nan, np.nan, np.nan])
# break down triangle into the individual points
v1, v2, v3 = triangle
eps = 0.000001
# compute edges
edge1 = v2 - v1
edge2 = v3 - v1
pvec = np.cross(ray_vec, edge2)
det = edge1.dot(pvec)
if abs(det) < eps: # no intersection
return False, null_inter
inv_det = 1.0 / det
tvec = ray_start - v1
u = tvec.dot(pvec) * inv_det
if u < 0.0 or u > 1.0: # if not intersection
return False, null_inter
qvec = np.cross(tvec, edge1)
v = ray_vec.dot(qvec) * inv_det
if v < 0.0 or u + v > 1.0: # if not intersection
return False, null_inter
t = edge2.dot(qvec) * inv_det
if t < eps:
return False, null_inter
return True, np.array([t, u, v])
主函数，计算交点
# Create a basic triangle within pyvista
points = np.array([[0, 0, 1], [0, 1, 0], [1, 0, 0]])
faces = np.array([3, 0, 1, 2])# 三个点，points的0，1，2点
tri = pv.PolyData(points, faces)# faces是空间三角形面，points是三角形顶点
# cast a ray above pointed downwards
start = np.array([0, 0, 0])
direction = np.array([2, 2, 2])
# compute if the intersection exists
inter, tuv = ray_triangle_intersection(start, direction, points)
t, u, v = tuv
print('Intersected', inter)
print('t:', t)
print('u:', u)
print('v:', v)
a, b, c = (1 - u - v), u, v
point = tri.points[0] * a + tri.points[1] * b + tri.points[2] * c
print(point)
可视化
if inter:
# reconstruct intersection point in barycentric coordinates. See
# https://en.wikipedia.org/wiki/Barycentric_coordinate_system
a, b, c = (1 - u - v), u, v
point = tri.points[0] * a + tri.points[1] * b + tri.points[2] * c
pl = pv.Plotter()
pl.add_text(f'Intersected at ({point[0]:.3}, {point[0]:.3}, {point[0]:.3})', font_size=26)
pl.add_mesh(tri)
_ = pl.add_arrows(
np.array([start]),
np.array([direction]),
show_scalar_bar=False,
color='r',
style='wireframe',
)
pl.add_points(np.array([point]), point_size=20, render_points_as_spheres=True, color='b')
pl.add_point_labels(tri, [f'a = {1 - u - v:.3}', f'b = {u:.3}', f'c = {v:.3}'], font_size=40)
pl.show_bounds()
pl.camera_position = 'xy'
pl.show()
else: # no intersection
pl = pv.Plotter()
pl.add_text('No intersection')
_ = pl.add_arrows(
np.array([start]),
np.array([direction]),
show_scalar_bar=False,
color='r',
style='wireframe',
)
pl.add_mesh(tri)
pl.show_bounds()
pl.camera_position = 'xy'
pl.show()
原理
射线的参数： o是射线的起点，D是射线的方向。
如果有一个点从起点o出发，沿着方向D移动任意长度t，到达终点R。如果t不同那么R也不同。通过无数不同的R，我们就能构成整条射线。
如下图所示，射线的起点为P0，沿着方向u前进（红线），P0+tu就构成了整条射线。
通过上述概念，我们可以推广到三角面上，即用两个方向向量|AB|、|AC|和其权重u、v来表示一个空间上的三角形（三个顶点为A，B，C）。
如上图所示，我们给定三角形的参数方程，（1-u-v）V0+uV1+vV2。上图三角形内任意一点都可以由向量和其权重值来表示。比如A点就是在两个方向上移动了0个单位距离。点C在向量|AC|方向移动了1个距离单位。
点P就是沿着AC方向移动了一段距离，又沿着AB移动了一段距离，他们的和向量|AP|就是P点走的方向。
那么三角形内的所有点都由参数u和v控制就可以得到。
于是求解射线和三角形的交点就是求解下面方程的未知数t、u、v，其他都是已知的。
表示为线性方程组：
下面解这个方程组，用到两个知识点，一是克莱姆法则，二是向量的混合积。
令E1 = V1 - V0，E2 = V2 - V0，T = O - V0，得
根据克莱姆法则，可得到t,u,v的解
将这三个解联合起来，得
根据混合积公式：
得
令
得：
代码实现：利用numpy计算和pyvista模块可视化
https://docs.pyvista.org/examples/99-advanced/ray-trace-moeller.html
