——讀Computer Graphics Principles and Practice 3rd Edition第七章時遇見課文正文和代碼中的錯誤，作記。
本文旨在闡釋一種算法，用於在三維空間中尋找某一線(ray)與某一三角形的交點。此算法是計算機圖形學中的基礎算法之一。
1.預設概念
為了闡釋此算法，必須先引入一組預設概念，借以使用計算機語言來描述三維空間中的線與三角形。
我們首先給出這些概念的定義及數據結構。
1.1 點(Point3D)
我們使用笛卡爾坐標系 (Cartesian coordinates)來描述三維空間。具體地，我們使用的坐標系都是右手系(right-handed coordinate system)。
我們規定，三維空間中的點總是相對於某個坐標系的。因此，我們通過其在某一坐標系中的坐標(x, y, z)來定義它。
public struct Point3D { public Double X { get; set; } public Double Y { get; set; } public Double Z { get; set; } // 構造函數及運算符重載略去 }
1.2 坐標矢量(Vector3D)
坐標矢量(coordinate vector)的概念源於線性代數中的向量(Vector)。
在三維空間中，我們認為坐標矢量v是定義在R3上的一個向量，用於描述一個點在三維空間中發生的運動。
注意坐標矢量與點的不同：坐標矢量表述位置的變化，而非存在，因此，不同於點，它是獨立於坐標系選取的。
public struct Vector3D { public Double X { get; set; } public Double Y { get; set; } public Double Z { get; set; } // 構造函數及運算符重載略去 }
1.3 點與坐標矢量的運算
1.3.1 坐標矢量的運算
向量中常用的運算，包括求取長度、向量加法、標量乘法、矢量叉乘與矢量點乘等，對於坐標矢量Vector3D來說同樣有意義。
public struct Vector3D { ... public static Double GetLength(Vector3D v) { return Math.Sqrt(v.X*v.X+v.Y*v.Y+v.Z*v.Z); } public static Vector3D GetSum(Vector3D v, Vector3D w) { return new Vector3D(v.X+w.X, v.Y+w.Y, v.Z+w.X); } public static Vector3D GetScalaProduct(Vector3D v, Double t) { return new Vector3D(v.X*t, v.Y*t, v.Z*t); } public static Vector3D GetCrossProduct(Vector3D v, Vector3D w) { return new Vector3D(v.Y*w.Z - v.Z*w.Y, v.X*w.Z - v.Z*w.X, v.X*w.Y - v.Y*w.X); } public static Doble GetDotProduct(Vector3D v, Vector3D w) { return v.X*w.X + v.Y*w.Y + v.Z*w.Z; } }
一些常用運算的幾何概念：
1.兩個不共線坐標矢量v和w通過線性組合(linear combination)所構成的新矢量k= xv + yw與v與w所標定的平面共面。
2.兩個不共線坐標矢量v和w通過叉乘所得的新矢量k= v × w垂直於v與w所標定的平面，這三個向量構成一個右手系。
1.3.2 點與坐標矢量的交互運算
點表示存在，而坐標矢量表示變化。因而，我們可以給出一組點與矢量交互運算的定義：
1.某點P可以與一個坐標矢量v相加，得到一個新點Q = P + v。Q即是點P經歷了v所描述的運動之后所到達的新點。
2.任意兩點P與Q的差值P - Q是一個坐標矢量。由1，我們可得P = Q + ( P - Q)
注意，點僅僅只有這兩個基本運算。兩點之和是沒有意義的：兩個描述存在的概念相加，在數值上確實可以得到一個結果，但這個結果並不獨立於坐標系，因而在應用中並沒有意義。
1.3.3 點的仿射組合(affine combination)
有一類點在數值上的線性組合，可以通過變形將其轉換為有意義的基本運算。具體地：
對於一組點P1, P2, ... , Pn，以及一組標量x1, x2, ... , xn; 如果(x1+ x2+ ... + xn) = 1，由於
x1P1 + x2P2 + ... + xnPn = P1 + x2( P2 - P1) + x3( P3 - P1) + ... + xn( Pn - P1), 而右式是一個合格的結果為點的基本運算表達式；
因而，我們將上述線性組合稱為點P1, P2, ... , Pn關於標量x1, x2, ... , xn;的仿射組合，並賦予其定義。
不難看出，兩個不同點的仿射組合描述了其所構成直線上的所有點；三個不共線點的仿射組合描述了其所構成平面上的所有點。
1.4 線
由1.3.3可知，我們能夠通過仿射組合αP + βQ s.t. α+β = 1的形式來描述線。
定義(Q - P)為方向矢量d，我們能夠得出更為簡潔的描述：ι : R→Point3D , t→ P + td
public struct Line { public Point3D P { get; private set; } public Vector3D V { get; private set; } public Point3D GetPointAt (Double t) { return P + t * V; } public Boolean Contain (Point3D point) { Double t; if(V.X != 0) { t = (point.X - P.X)/V.X; return point == GetPointAt(t); } else if(V.Y != 0) { t = (point.Y - P.Y)/V.Y; return point == GetPointAt(t); } else if(V.Y != 0) { t = (point.Y - P.Y)/V.Y; return point == GetPointAt(t); } else throw new InvalidLineException(); } // 細節略去 }
限制t的取值，我們即可得到射線或線段。
1.5 平面
標定三維空間中平面的方式有很多種。我們選取最為簡潔的點+法向量(normal vector)模式。此模式下，平面由平面上的某個點P，以及垂直於它的法向量n來標定。我們規定，n必須標准化(normalized)。
此模式下我們非常容易確定某給定的點或線是否處於該平面內。
public struct Plane { public Point3D P { get; private set; } public Vector3D N { get; private set; } public Boolean Contain(Point3D point) { return (point == this.P) || (Vector3D.GetDotProduct(point - this.P, this.N) == 0); } public Boolean Contain(Line line) { return (Vector3D.GetDotProduct(line.V, this.N) == 0 && this.Contain(line.P)); } // 構造函數等略去 }
1.6 三角形
我們通過三角形的三個頂點標定它。
注意到，我們經常需要訪問三角形所處的平面。因此，我們同時在它的數據結構中存儲一個平面。
public struct Triangle { public Point3D A { get; private set; } public Point3D B { get; private set; } public Point3D C { get; private set; } public Plane Plane { get; private set; } public Triangle (Point3D a, Point3D b, Point3D c) { A = a; B = b; C= c; Plane = new Plane(A, Vector3D.GetCrossProduct(B-A, C-B);
// ... } // .. }
2.相交判定
有了上述定義，我們的問題可以被描述為，給定線Line line, 三角形Triangle triangle, 尋找一個他們的交點Point3D? intersectionPoint.
此方法約定，當二者不相交時返回null；當二者相交於一個點時，返回該交點；當二者相交於無數點時，返回其中任意一個交點。
我們的實現大致分為兩步。
第一步，確定line與triangle是否共面。
第二步，如果共面，則依次測試line與triangle的三條邊線段是否相交。【Line-Line Segment Intersection Test】如果相交則返回交點，如果不相交則返回null。
如果不共面，則測試line是否與tirangle所在平面相交。【Line-Plane Intersection Test】如果不相交，返回null；如果相交，則繼續測試交點是否落在triangle內。【Point-in-Triangle Test】如果是，返回該交點；否則返回null。
注意到我們的相交判定主體被分成了線三角形共面檢測、 (共面的)線線相交檢測、(三維的)線面相交檢測、(共面的)點在三角形內檢測這幾部分。
下面分別敘述這幾個過程。
2.1 線-三角形共面(Line-Triangle coplanarity)檢測
三角形通過Triangle.Plane性質確定了其所在面。因此，所述問題等價於判斷給定的Line line是否處於給定的Plane plane上。
如果line在plane上，則line與plane的法向量垂直，且line上任意一點與plane上任意一點的連線也與plane的法向量垂直。
2.2 線線相交檢測
本算法僅僅涉及到共面的直線與線段的相交檢測。但為通用起見，首先介紹如何檢驗三維空間中的兩條線Line l, j是否可能共面。
此過程非常簡單，僅僅需要首先通過線A與線B上的某點確定一個面（實際上僅僅用到該面的法線），再確定線B是否同樣屬於該面即可：
public struct Line { // ... public static Plane? Coplane(Line l, Line j) { Vector3D n = Vector3D.GetCrossProduct(l.V, j.P - l.P); if (Vector3D.GetDotProduct(n, j.V) == 0) { return new Plane(l.V, n); } else { return null; } } }
下面，如何確定兩條共面直線是否相交、交於何處呢？精髓在於首先確定二者是否平行。如果平行則進一步判斷共線，否則通過方程計算出交點：
public struct Line { // ... public Double? IntersectAt(Line l) { Plane? p = Coplane(this, l); if(p != null) { Vector3D perpendic = Vector3D.GetCrossPoduct(p.N, this.V); Double testResult = Vector3D.GetDotProduct(perpendic, l.V); if(testResult !=0) { return Vector3D.GetDotProduct(this.P - l.P, perpendic)/testResult; } else if(l.Contain(this.P)) { return 0; } else { return null; } } else return null; } public static Point3D? Intersect(Line l, Line j) { Double? result = l.IntersectAt(j); if(result != null) { return l.GetPointAt(result.Value); } else return null; } }
那么，如何確定直線與線段的交點呢？
由點A、點B確定的線段AB可以看作一個受限制的Line實例。我們將此實例的P性質初始化為A，V性質初始化為(B - A)，則當t屬於區間[0, 1]時，方法GetPointAt(t)所返回的點即是線段上的點。依據這一敘述，我們可以在Line類中定義相應的方法。
有了此方法后，給定三角形triangle以及某線段line，如果此時已經確定了該線段在triangle所在的平面上，我們只需
public struct Line { // ... public Double? IntersectAt(Triangle triangle) { // ... // 此時，已經確定此線段與三角形共面 Line AB = new Line(triangle.A , triangle.B - triangle.A); Double? result = AB.IntersectAsLineSegmentAt(this); if (result == null) { Line AC = new Line(triangle.A , triangle.C - triangle.A); result = AC.IntersectAsLineSegmentAt(this); if(result == null) { Line BC = new Line(triangle.B , triangle.C - triangle.B); result = BC.IntersectAsLineSegmentAt(this); } } return result; // ... } public Double? IntersectAsLineSegmentAt(Line l) { Double? result = this.IntersectAt(l); if(result != null && result<=1 && result >=0) { return result; } return null; }
2.3 線面相交檢測
如果線在面上，任意返回一個線上的點即可；如果線面平行，返回null；否則，線面必交於一點，我們通過面的法向量決定交點。
public struct Line { // ... public Double? IntersectAt(Plane plane) { Double dotP = Vector3D.GetDotProduct(this.V, plane.N); if(dotP == 0) { if(plane.Contain(this)) return 0; else return null; } else { return Vector3D.GetDotProduct(plane.P - this.P, this.N)/dotP; } } public Point3D? Intersect(Plane plane) { Double? t = this.Intersect(plane); if(t == null) return null; else { return this.GetPointAt(t.Value); } } } pulic struct Plane { // ... public Point3D? Intersect(Line l) { return l.Intersect(this); } }
2.4 點在三角形內檢測
回憶一下2.2節，在定義線與線段相交的時候，我們談到如何判定某直線上的點P屬於給定的線段AB：通過AB坐標重定義該直線的表達，並在新表達下檢查點P對應的t取值是否屬於[0,1]。
事實上，我們那時所取的t值可以看作是由點A、B所確定的仿射變換 tA + (1 - t)B 中的參數t。當t屬於[0, 1]時，該仿射變換中所有的系數都大於等於0，因而最終所確定的點必定在線段AB上。
這一原理同樣適用於三角形：
給定三角形的三個頂點A、B、C，我們可以得到一個仿射變換αA + βB + γC s.t. α + β + γ = 1. 當所有系數α、β、γ均大於等於0時，由此仿射變換所確定的點P必將落在該三角形內部（等於0則對應邊界）。
由此定義可以衍生出三角形所在平面的質心坐標系(Barycentric Coordinates)：
給定三角形ABC；對於其所處平面上的任意一點P，有且僅有一組(α, β, γ)為點P相對於△ABC質心坐標，使得α + β + γ = 1，且在任意笛卡爾坐標系下，A、B、C的坐標確定之后的仿射變換αA + βB + γC定義了點P在該笛卡爾系中的坐標。
理解質心坐標系的核心在於理解仿射變換。如果仍不太明白，建議多類比一下僅僅包含兩點的仿射變換tA + (1 - t)B及其幾何意義：線段。
有了質心坐標的相關知識，回到我們的問題：想要檢測面上的一點P是否包含於面上的三角形ABC，我們只需要逆推點P關於△ABC的質心坐標，並檢驗是否其所有系數均大於0即可。那如何逆推坐標呢？
觀察仿射變換αA + βB + γC在代入了約束條件之后的形態：P = A + β(B - A) + γ(C - A).
移項之后，我們得到一個僅僅關於坐標矢量的美好式子：(P - A ) = β(B - A) + γ(C - A). 將此式展開到實數層面，我們將會得到一個關於β、γ的一組方程。由於P與ABC共面，此方程組必有解。然而由於可能存在的0系數，此解難以通過計算機語言得出。
我們的替代方案是，在Triangle類的構造函數中預計算兩個輔助量，並通過該輔助量在矢量層面上消元以確定質心坐標。
繼續考察式子(P - A ) = β(B - A) + γ(C - A)。我們注意到，如果在式子左右同時點乘坐標矢量n，同樣可能將其轉化為實數方程：
(P - A )·n = β(B - A)·n + γ(C - A)·n
因此，我們可以通過預計算兩個坐標矢量v、w，用以與此式相乘並消元：
預計算v s.t. v 垂直於(B - A) 且(C - A)·v = 1; 則 (P - A)·v = γ；
預計算w s.t. w 垂直於(C - A) 且(B - A)·w = 1; 則 (P - A)·w = β.
最后，計算α = 1 - β - γ; 即可用以確定該點是否在三角形內。
public struct Triangle { // ... private Vecter3D _v; pirvate Vector3D _w; public Triangle (Point3D a, Point3D b, Point3D c) { // ... // 預計算_v與_w _v = Vector3D.GetCrossProduct(this.Plane.N, this.B - this.A); _v /= Vector3D.GetDotProduct(this.C - this.A, _v); _w = Vector3D.GetCrossProduct(this.Plane.N, this.C - this.A); _w /= Vector3D.GetDotProduct(this.B - this.A, _w); } public Boolean Contain(Point3D point) { if(this.Plane.Contain(point)) { Vector3D AP = point - this.A; Double gamma = Vector3D.GetDotProduct(AP, _v); if(gamma >=0 && gamma <= 1) { Double beta = Vector3D.GetDotProduct(AP, _w); if(beta >=0 && beta <= 1) { Double alpha = 1 - gamma - beta; if(alpha >= 0 && alpha <=1) return true; } } } return false; } }
3.最終實現
public struct Line { // ... public Double? IntersectAt(Triangle triangle) { Double? result=null; if(triangle.Plane.Contain(this)) { Line AB = new Line(triangle.A , triangle.B - triangle.A); result = AB.IntersectAsLineSegmentAt(this); if (result == null) { Line AC = new Line(triangle.A , triangle.C - triangle.A); result = AC.IntersectAsLineSegmentAt(this); if(result == null) { Line BC = new Line(triangle.B , triangle.C - triangle.B); result = BC.IntersectAsLineSegmentAt(this); } } } else { result = this.IntersectAt(triangle.Plane); if(result != null) { Point3D point = this.GetPositionAt(result); if(!triangle.Contain(point)) result = null; } } return result; } public Point3D? Intersect(Triangle triangle) { Double? t = this.IntersectAt(triangle); if(t == null) return null; else return this.GetPointAt(t.Value); } } public class Triangle { // ... public Point3D? Intersect( Line line) { return line.Intersect(this); } }