線性物件與三角形的相交測試
本篇文章介紹一種三維空間下線性物件與三角形的相交測試演算法和實現，線性物件包括：直接、射線和線段，它們可以用類似的形式來表示，先以直接為例來闡述演算法，射線和線段是直接的特殊情況。
文章目錄：
演算法介紹
原始碼實現
參考
演算法介紹
在三維空間上，判斷線性物件
與三角形
是否相交，若相交，求出交點，三角形所在的平面是
。
圖1。 直線與三角形的四種關係
三維空間的模型通常是由一個個三角形面片構成的，線性物件與三角形的相交檢測也顯得特別的重要。如圖1所示，一條直線與三角形可能有四種關係：1）直線
與三角形相交；2）直線
與三角形所在的平面相交，但與三角形不相交；3）直線
與在三角形所在的平面上；4）直線
與三角形所在的平面平行。 一種樸素的檢測演算法是計算線性物件與三角形所在的平面的交點，再判斷交點是否在三角形內。如果交點在三角形內，則線性物件與三角形相交；否則，不相交。這種演算法效率不僅效率較低，而且更容易造成浮點數誤差。因為計算機中的浮點數佔用位元組有限，能表示的精度有限，有效數較少，求直線與平面的交點，由於對同一個變數多次的乘法或者除法操作，造成誤差的疊加，最後計算出來的結果由於精度問題，可能造成較大的誤差。例如，三角形頂點的座標數值較大，計算的直線與平面的交點，會由於精度誤差導致求出的交點偏離平面，即可以檢測到交點不在平面上。
Möller&Trumbore（1997）［1］提出了一種更加高效的射線與三角形相交測試的演算法，而且乘法或者除法操作造成的誤差疊加較少，精度損失更小，因此效能上更加健壯。
三角形內的任意一點都可以用它相對於三角形的頂點的位置來定義：
其中，
，稱為重心座標（barycentric coordinate）。直線可以用引數方程表示為
，因此計算直線與三角的交點的等式為
整理，得
![\left[ \begin{matrix} -\vec{d} & {{V}_{1}}-{{V}_{0}} & {{V}_{2}}-{{V}_{0}} \\ \end{matrix} \right]\left[ \begin{matrix} t \\ u \\ v \\ \end{matrix} \right]=\left[ P-{{V}_{0}} \right] \tag{3}\\](https://img.heatask.com/upload/website_attach/I4UDk7tLuRHr6=TAV=RL+sbNdDK10FRf6bNDw_vATlw07O7PtegEffKDoZQtpUiExvngUphALXJuJd1tx=BtpYoipTHJpUrgiWzFU_tPRsoiqW-zawW_hwNgg=Y3gDaJqkIiTteDXxN0ivAPtEUFpYoRoZktaViAkOmg=nMbLXJ2qW-zawW_hwNgg=Y3gDaJqkIiTteDXxn0ivAPtEUFpYoRoZktaViAkOmg=nMbLXJ2qW-zawEBhauggb7JobRh+bt+aOr_neOoK4GEZcJNqSbEUEZ+eFBDqkZiaOJynua66hP=OWEBhaQVe+k++ViAx25kV7mHqiWrkhL1awF5qSNggb73gdn4r_D4VNy1LXa1qWhyawEmqDJ4pTh8UUiEx_Dg=yTPn5U0ivVNT6w_++Ugg=OJobRhK2A0w7ZPRZU03vAPtegEffKDoZQtpEvyqkIRaOr_+UWrkyePtEnXpYoioZk5pUiExsn2UpVtLXa2rQP=OBwRfSbAUV73gDPgiW0g=ytOS.jpeg)
這相當於是一個齊次線性方程組，
是它的解，可以採用克拉姆法則來求解，得到
![\left[ \begin{matrix} t \\ u \\ v \\ \end{matrix} \right]=\frac{1}{\left| \begin{matrix} -\vec{d} & {{E}_{1}} & {{E}_{2}} \\ \end{matrix} \right|}\left[ \begin{matrix} |T & {{E}_{1}} & {{E}_{2}}| \\ |-\vec{d} & T & {{E}_{2}}| \\ |-\vec{d} & {{E}_{1}} & T| \\ \end{matrix} \right] \tag{4}\\](https://img.heatask.com/upload/website_attach/QIMoQf+x+2jptUwYU78Xtwg9SneLa7i2Tzfo_PeovJMG=apd5g7qsvTnYCg3-uywF6C0iJmoVDD4v24wT7e7C0x1Zv9c-nBiV50pD3OJA2ZdAKHxWlb9K26HsO+c-uMiTtZfwEgLVD4HAKgxAruos6gqNI7=_JWg-10rEH995FKsu2xg-v2JM6g=O3O7t7Rn_zR02+xNVDD0FLxKVyuL=UjKx3OotuMc4CKFiPTXSu0prjqr-oOXE6gqO39c-uMakA41Rdg=VDD45fgVBvb9RMSeYCT3TBUJAF0rEMd95fBzFLxfVyFMNpgtsl9mHKVc4JC3wETM_20prqXd5pm9KxNeZvO7tuywF5Z02A49ScLG+SYbjRG4tsEHs6FcTBidSfcb_fF95fRG+SNd5g7qsvTnYCg3-uywF6C0iJmoVDD4v24wT7e7C0x1Zv9c-uykFKLfwExKW20progVA_29KxOxYCT3bAykA10rEH99IsvdFLxKVyF0MHgqO0dmHKVnkAD2wETMVDDHAKGd5gT9K26eZv9cTBUdfF0pDKj9j20prjNd5pmZC0XpZ7FcTBVxV50rDvO6AC0G+SOAVyF0R9gqO3Oo=2ykF5ZfwEgLVD4HAKG_Wl26LppHsO4nzgJc4JCiwETMW20s+oGd5pN9KxY2YCTHPuykFle02+F95fLdFLXkWqm9Kx6eZvO7tuywF5Z02A49ScLG+SYbjRG4tsEHs6FcTBidSfcb_fF95FLdFLje-nIVC0x5s3Oo=2ywF50pDBM8v.jpeg)
其中，
，依據線性代數的知識，易知
，因此等式（4）可以寫成：
![\left[ \begin{matrix} t \\ u \\ v \\ \end{matrix} \right]=\frac{1}{\left( \vec{d}\times {{E}_{2}} \right)\cdot {{E}_{1}}}\left[ \begin{matrix} \left( T\times {{E}_{1}} \right)\cdot {{E}_{2}} \\ \left( \vec{d}\times {{E}_{2}} \right)\cdot T \\ \left( T\times {{E}_{1}} \right)\cdot \vec{d} \\ \end{matrix} \right] \tag{5}\\](https://img.heatask.com/upload/website_attach/gMEfcGanm1X5_1_nSfvng2F2y3G891BT_xGfaVYszIsEj=+82C1uzZzMT-Z1JPDmtHwrZ1qs6lhOhkM+Qf92Tw_t-zDaJ2vn-NC+QPE3S179-OqIx5qtHkbvHw0aJP7n_zrXJSCF6lBF-OAID6P4zPFHtG5scGR6jFCJRBnEPEmggkjVR1BH=PFN=NVKj1PCyxerNHRV6lhATRjHwAP=h1N8_NV5NB7vHS2IZVaEPeRWTRjLwAP=_1r5zY0aIShETQRGRRE5T9WEHGT3wAQlTw_t-zVKj2R6zHsfJSRv6lBFy0loR1qtH3SvH1nGIShtPNCJQ_8EPevEHG=82JltHkH4zXdJISBE-NC+Q1SEb+psTRjHjzmGi0r0T-ImJP7n-NC+Q1AEbwvEbi5Yy1BHtu=6hXZuJPDHtzCJQxEEPesITRjHeFB1=PFH=NDLC-mCgxRkdBEgy=D3t00Nj7lpTw_ET-I1QRDH=lbrNHRg6lhATRjLx5qtHkbvHw0aJP7nTQy8bmCHfmDg0Ap82C10zZbvH1nPIShtTQy8dmqlbOZ9TRjHwAQA=UFH=AmLCHQCTQR9JSiA71D3t1-Wi6vGTwS1T-ZMX1s3_F8FKPnEPEZEHTcYx5qpTwNshHZfcBDC0F8FJSCFfcplzVcYwAQATw_E=TV7QEWvHmUBJSiA6lhA-OAID7Unz55MT-7oISBsyHEndBnEPEm4zkc82JUwTw_t-zDLC-ZvHSeXJSCFbNJWTRjHjzmGi0r0T-ImJPDmt0wVZ1426lBA-OAID7DzzhFHtQVLCHvvHSerN-aDz.jpeg)
如果得到的解
，那麼直線與平面的交點位於三角形內；否則，位於三角形外。注意，若
或
或
，則交點在三角形的邊上。
直線與三角形相交測試的演算法偽碼如下所示：
在三維空間上，判斷直線
與三角形
是否相交，若相交，求出交點，三角形所在的平面是
如果需要求射線與三角形的交，必需限定
的取值範圍，
必需滿足條件
，若計算出來的
，則射線與三角形不相交。同樣，如果需要求線段與三角形的交，線段的引數方程仍然可以表示為
，線段的兩個端點分別是
和
，則
必需滿足條件
，因此，若計算出來的
或
，則線段與三角形不相交。
注意到，在上述的演算法偽碼中，已經把一些中間變數儲存起來，避免重複計算，可以提高演算法的速度，例如中間變數
。
原始碼實現
基礎庫原始碼連結，參見
這裡
，下面是前面所描述的演算法的實現。
#include
“SrGeometricTools。h”
#include
“SrDataType。h”
namespace
{
class
SrSegment3D
{
public
：
SrPoint3D
mPoint1
；
SrPoint3D
mPoint2
；
}；
class
SrRay3D
{
public
：
SrPoint3D
mBase
；
SrPoint3D
mDirection
；
}；
class
SrLine3D
{
public
：
SrPoint3D
mBase
；
SrPoint3D
mDirection
；
}；
/**
\brief 3D triangle class。
This is a 3D triangle class with public data members。
*/
class
SrTriangle3D
{
public
：
/**
\brief Default constructor， endpoints is set to （0，0）。
*/
SrTriangle3D
（）
{
mPoint
［
0
］
=
mPoint
［
1
］
=
mPoint
［
1
］
=
SrVector3D
（
0
，
0
，
0
）；
}
/**
\brief The line is initialized by three points。
*/
SrTriangle3D
（
const
SrPoint3D
&
p0
，
const
SrPoint3D
&
p1
，
const
SrPoint3D
&
p2
）
{
mPoint
［
0
］
=
p0
；
mPoint
［
1
］
=
p1
；
mPoint
［
2
］
=
p2
；
}
/**
\brief Destructor。 Do nothing。
*/
~
SrTriangle3D
（）
{}
/**
\brief The triangle is valid if the three points are not on a line。
*/
bool
isValid
（）
const
{
SrVector3D
norm
=
normal
（）；
if
（
EQUAL
（
norm
。
x
，
0
）
&&
EQUAL
（
norm
。
y
，
0
）
&&
EQUAL
（
norm
。
z
，
0
））
return
false
；
return
true
；
}
/**
\brief Compute the normal using Newell method devised by Martin Newell。It relates to the order of the points。
*/
const
SrVector3D
normal
（）
const
{
SrVector3D
norm
=
（
mPoint
［
1
］
-
mPoint
［
0
］）。
cross
（
mPoint
［
2
］
-
mPoint
［
0
］）；
return
norm
；
}
/**
\brief Judge whether or not the line hits the triangle。
\return SR_OVERLAPPING if the line is on the plane where the triangle lies；
SR_PARALLEL if the line is parallel to the triangle。
SR_DISJOINT if the line misses the triangle；
SR_INTERSECTING if the line intersects with the triangle and intersection point is not on the edge。
\param［out］ result If the line intersects with the triangle and the intersection point information is put in ‘result’。
*/
int
lineHitTest
（
const
SrLine3D
&
line
，
SrPoint3D
&
/*［OUT］*/
result
）
const
{
SrVector3D
ret
；
int
retFlag
=
linearIntersectTriangle
（
line
。
mBase
，
line
。
mDirection
，
ret
）；
if
（
retFlag
！=
SR_INTERSECTING
）
return
retFlag
；
result
=
line
。
mBase
+
ret
。
z
*
line
。
mDirection
；
//check whether the intersection point is on the edge。
//if（ EQUAL（ret。x+ret。y-1。0，0） ||
// EQUAL（ret。x，0） ||
// EQUAL（ret。y，0） ）
// return SR_INTERSECTING_EDGE；
return
SR_INTERSECTING
；
}
/**
\brief Judge whether or not the ray hits the triangle。
\return SR_OVERLAPPING if the ray is on the plane where the triangle lies；
SR_PARALLEL if the ray is parallel to the triangle。
SR_DISJOINT if the ray misses the triangle；
SR_INTERSECTING if the ray intersects with the triangle and intersection point is not on the edge。
\param［out］ result If the ray intersects with the triangle and the intersection point information is put in ‘result’。
*/
int
rayHitTest
（
const
SrRay3D
&
ray
，
SrPoint3D
&
/*［OUT］*/
result
）
const
{
SrVector3D
ret
；
int
retFlag
=
linearIntersectTriangle
（
ray
。
mBase
，
ray
。
mDirection
，
ret
）；
if
（
retFlag
！=
SR_INTERSECTING
）
return
retFlag
；
//check t。 t>=0
if
（
LESS
（
ret
。
z
，
0
））
return
SR_DISJOINT
；
result
=
ray
。
mBase
+
ret
。
z
*
ray
。
mDirection
；
//check whether the intersection point is on the edge。
//if（ EQUAL（ret。x+ret。y-1。0，0） ||
// EQUAL（ret。x，0） ||
// EQUAL（ret。y，0） ）
// return SR_INTERSECTING_EDGE；
return
SR_INTERSECTING
；
}
/**
\brief Judge whether or not the segment hits the triangle。
\return SR_OVERLAPPING if the segment is on the plane where the triangle lies；
SR_PARALLEL if the segment is parallel to the triangle。
SR_DISJOINT if the segment misses the triangle；
SR_INTERSECTING if the segment intersects with the triangle and intersection point is not on the edge。
\param［out］ result If the segment intersects with the triangle and the intersection point information is put in ‘result’。
*/
int
segmentHitTest
（
const
SrSegment3D
&
segment
，
SrPoint3D
&
/*［OUT］*/
result
）
const
{
SrVector3D
ret
；
SrVector3D
direction
=
segment
。
mPoint2
-
segment
。
mPoint1
；
int
retFlag
=
linearIntersectTriangle
（
segment
。
mPoint1
，
direction
，
ret
）；
if
（
retFlag
！=
SR_INTERSECTING
）
return
retFlag
；
//check t。 0<= t <=1
if
（
LESS
（
ret
。
z
，
0
）
||
GREATER
（
ret
。
z
，
1
））
return
SR_DISJOINT
；
result
=
segment
。
mPoint1
+
ret
。
z
*
direction
；
//check whether the intersection point is on the edge。
//if（ EQUAL（ret。x+ret。y-1。0，0） ||
// EQUAL（ret。x，0） ||
// EQUAL（ret。y，0） ）
// return SR_INTERSECTING_EDGE；
return
SR_INTERSECTING
；
}
public
：
SrPoint3D
mPoint
［
3
］；
private
：
/**
\brief Judge whether or not the line hits the triangle。The line is base + t*direction。
It is not necessary for direction to be unit length。
\return SR_OVERLAPPING if the line is on the plane where the triangle lies；
SR_PARALLEL if the line is parallel to the triangle。
SR_DISJOINT if the line misses the triangle；
SR_INTERSECTING if the line intersects with the triangle。
\param［out］ result return the result of （u，v，t） if the line intersects with the triangle。
*/
/*
——————————————————————————————————————
Devised by Moller and Trumbore，1997。 intersection> Any point in a triangle can be defined in terms of its position relative to the triangle’s vertices： Qu，v = （1-u-v）V0 + uV1 + vV2 ，0<=u<=1，0<=v<=1，0<=u+v<=1 For the linear component–triangle intersection， P + t*^d = Qu，v which can be expanded and applied Cramer’s rule to： |t| | |P-V0 V1-V0 V2-V0| | |u| = 1/|-^d V1-V0 V2-V0| | |-^d P-V0 V2-V0| | |v| | |-^d V1-V0 P-V0| | | （（P-V0）×（V1-V0））·（V2-V0） | = 1/（^d×（V2-v0））·（V1-V0） | （^d×（V2-V0））·（P-V0） | | （（P-V0）×（V1-V0））·^d | —————————————————————————————————————— */ int linearIntersectTriangle （ const SrPoint3D & base ， const SrVector3D & direction ， SrVector3D & result ） const { SrReal u ， v ， tmp ； SrVector3D e1 ， e2 ， p ， s ， q ； e1 = mPoint ［ 1 ］ - mPoint ［ 0 ］； e2 = mPoint ［ 2 ］ - mPoint ［ 0 ］； p = direction 。 cross （ e2 ）； tmp = p 。 dot （ e1 ）； //If the line is perpendicular to the normal of triangle。 if （ EQUAL （ tmp ， 0 ）） { //The line is on the plane。 p = e1 。 cross （ e2 ）； if （ EQUAL （ p 。 dot （ base - mPoint ［ 0 ］）， 0 ）） return SR_OVERLAPPING ； //The line is parallel to the plane。 return SR_PARALLEL ； } s = base - mPoint ［ 0 ］； u = p 。 dot （ s ） / tmp ； if （ LESS （ u ， 0 ） || GREATER （ u ， 1 ）） return SR_DISJOINT ； q = s 。 cross （ e1 ）； v = q 。 dot （ direction ） / tmp ； if （ LESS （ v ， 0 ） || GREATER （ v ， 1 ） || GREATER （ u + v ， 1 ）） return SR_DISJOINT ； result 。 x = u ； result 。 y = v ； result 。 z = e2 。 dot （ q ） / tmp ； return SR_INTERSECTING ； } }； } 參考 ［1］ Tomas Möller， and Ben Trumbore。 “Fast， minimum storage ray-triangle intersection。” Journal of graphics tools， vol。2， no。1， pp。21-28， 1997。 ［2］ Philip J。 Schneider， and David H。 Eberly。 Geometric Tools for Computer Graphics。 Morgan Kaufmann， 2002。
