三維空間中線與三角形相交判定


——讀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.兩個不共線坐標矢量vw通過線性組合(linear combination)所構成的新矢量k= xv + ywvw所標定的平面共面

2.兩個不共線坐標矢量vw通過叉乘所得的新矢量k= v × w垂直於vw所標定的平面,這三個向量構成一個右手系。

 

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

因此,我們可以通過預計算兩個坐標矢量vw,用以與此式相乘並消元:

預計算v s.t. v 垂直於(B - A) 且(C - A)·= 1; 則 (P - A)·= γ;

預計算w s.t. w 垂直於(C - A) 且(B - A)·= 1; 則 (P - A)·= β.

最后,計算α = 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);
}
}

 


注意!

本站转载的文章为个人学习借鉴使用,本站对版权不负任何法律责任。如果侵犯了您的隐私权益,请联系我们删除。



 
粤ICP备14056181号  © 2014-2021 ITdaan.com