三角形平面内任意点P
可由三角形三个顶点表示:
P
=
w
A
+
u
B
+
v
C
P = wA + uB + vC
P=wA+uB+vC
根据重心坐标性质w + u + v = 1
:
P
=
(
1
−
u
−
v
)
A
+
u
B
+
v
C
P = (1 - u - v)A + uB + vC
P=(1−u−v)A+uB+vC
P
=
A
−
u
A
−
v
A
+
u
B
+
v
C
=
A
+
u
(
B
−
A
)
+
v
(
C
−
A
)
P=A - uA - vA + uB + vC = A + u(B - A) + v(C - A)
P=A−uA−vA+uB+vC=A+u(B−A)+v(C−A)
光线表示为:
P
=
O
+
t
D
P=O+tD
P=O+tD
得到:
O
+
t
D
=
A
+
u
(
B
−
A
)
+
v
(
C
−
A
)
O
−
A
=
−
t
D
+
u
(
B
−
A
)
+
v
(
C
−
A
)
\begin{array}{l} O+tD & = & A + u(B - A) + v(C - A)\\ O-A & = & -tD+u(B-A)+v(C-A) \end{array}
O+tDO−A==A+u(B−A)+v(C−A)−tD+u(B−A)+v(C−A)
写成矩阵形式:
[
−
D
(
B
−
A
)
(
C
−
A
)
]
[
t
u
v
]
=
O
−
A
\begin{bmatrix} -D & (B-A) & (C-A) \end{bmatrix} \begin{bmatrix} t\\u\\v \end{bmatrix} =O-A
[−D(B−A)(C−A)]⎣⎡tuv⎦⎤=O−A
令:
T
=
O
−
A
E
1
=
B
−
A
E
2
=
C
−
A
T = O - A\\E_1 = B - A \\E_2 = C - A
T=O−AE1=B−AE2=C−A
根据克莱姆法则:
[ t u v ] = 1 [ ∣ − D E 1 E 2 ∣ ] [ ∣ T E 1 E 2 ∣ ∣ − D T E 2 ∣ ∣ − D E 1 T ∣ ] , \left[ \begin{array}{r} t \\ u \\ v\end{array}\right] = {1 \over \left[ \left| \begin{array}{r} -D & E_1 & E_2 \end{array}\right| \right]} \left[ \begin{array}{c} \left| \begin{array}{c} T & E_1 & E_2\end{array}\right| \\ \left| \begin{array}{c} -D & T & E_2\end{array}\right| \\ \left| \begin{array}{c} -D & E_1 & T\end{array}\right| \\ \end{array}\right], ⎣⎡tuv⎦⎤=[∣∣−DE1E2∣∣]1⎣⎡∣∣TE1E2∣∣∣∣−DTE2∣∣∣∣−DE1T∣∣⎦⎤,
再令:
P
=
(
D
×
E
2
)
Q
=
(
T
×
E
1
)
P=(D \times E_2) \\Q=(T \times E_1)
P=(D×E2)Q=(T×E1)
最终得到:
[ t u v ] = 1 ( D × E 2 ) ⋅ E 1 [ ( T × E 1 ) ⋅ E 2 ( D × E 2 ) ⋅ T ( T × E 1 ) ⋅ D ] = = 1 P ⋅ E 1 [ Q ⋅ E 2 P ⋅ T Q ⋅ D ] , \left[ \begin{array}{r} t \\ u \\ v\end{array}\right] = {1 \over {(D \times E_2) \cdot E_1}} \left[ \begin{array}{c} (T \times E_1) \cdot E_2 \\ (D \times E_2) \cdot T \\ (T \times E_1) \cdot D \end{array}\right] = = {1 \over {P \cdot E_1}} \left[ \begin{array}{c} Q \cdot E_2 \\ P \cdot T \\ Q \cdot D \end{array}\right], ⎣⎡tuv⎦⎤=(D×E2)⋅E11⎣⎡(T×E1)⋅E2(D×E2)⋅T(T×E1)⋅D⎦⎤==P⋅E11⎣⎡Q⋅E2P⋅TQ⋅D⎦⎤,
代码如下
bool rayTriangleIntersect(
const Vec3f &O, const Vec3f &D,
const Vec3f &A, const Vec3f &B, const Vec3f &C,
float &t, float &u, float &v)
{
Vec3f E1 = B - A;
Vec3f E2 = C - A;
Vec3f P = D.crossProduct(E2);
float det = E1.dotProduct(P);
#ifdef CULLING
// if the determinant is negative the triangle is backfacing
// if the determinant is close to 0, the ray misses the triangle
if (det < kEpsilon) return false;
#else
// ray and triangle are parallel if det is close to 0
if (fabs(det) < kEpsilon) return false;
#endif
float invDet = 1 / det;
Vec3f T = O - A;
u = T.dotProduct(P) * invDet;
if (u < 0 || u > 1) return false;
Vec3f Q = T.crossProduct(E1);
v = D.dotProduct(Q) * invDet;
if (v < 0 || u + v > 1) return false;
t = E2.dotProduct(Q) * invDet;
return true;
}