标签：
——读Computer Graphics Principles and Practice 3rd Edition第七章时遇见课文正文和代码中的错误，作记。
本文旨在阐释一种算法，用于在三维空间中寻找某一线(ray)与某一三角形的交点。此算法是计算机图形学中的基础算法之一。
为了阐释此算法，必须先引入一组预设概念，借以使用计算机语言来描述三维空间中的线与三角形。
我们首先给出这些概念的定义及数据结构。
我们使用笛卡尔坐标系 (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; } // 构造函数及运算符重载略去 }
坐标矢量(coordinate vector)的概念源于线性代数中的向量(Vector)。
在三维空间中，我们认为坐标矢量v是定义在R3上的一个向量，用于描述一个点在三维空间中发生的运动。
注意坐标矢量与点的不同：坐标矢量表述位置的变化，而非存在，因此，不同于点，它是独立于坐标系选取的。
public struct Vector3D { public Double X { get; set; } public Double Y { get; set; } public Double Z { get; set; } // 构造函数及运算符重载略去 }
向量中常用的运算，包括求取长度、向量加法、标量乘法、矢量叉乘与矢量点乘等，对于坐标矢量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.某点P可以与一个坐标矢量v相加，得到一个新点Q = P + v。Q即是点P经历了v所描述的运动之后所到达的新点。
2.任意两点P与Q的差值P - Q是一个坐标矢量。由1，我们可得P = Q + ( P - Q)
注意，点仅仅只有这两个基本运算。两点之和是没有意义的：两个描述存在的概念相加，在数值上确实可以得到一个结果，但这个结果并不独立于坐标系，因而在应用中并没有意义。
有一类点在数值上的线性组合，可以通过变形将其转换为有意义的基本运算。具体地：
对于一组点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.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的取值，我们即可得到射线或线段。
标定三维空间中平面的方式有很多种。我们选取最为简洁的点+法向量(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)); } // 构造函数等略去 }
我们通过三角形的三个顶点标定它。
注意到，我们经常需要访问三角形所处的平面。因此，我们同时在它的数据结构中存储一个平面。
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);
// ... } // .. }
有了上述定义，我们的问题可以被描述为，给定线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。
注意到我们的相交判定主体被分成了线三角形共面检测、 (共面的)线线相交检测、(三维的)线面相交检测、(共面的)点在三角形内检测这几部分。
下面分别叙述这几个过程。
三角形通过Triangle.Plane性质确定了其所在面。因此，所述问题等价于判断给定的Line line是否处于给定的Plane plane上。
如果line在plane上，则line与plane的法向量垂直，且line上任意一点与plane上任意一点的连线也与plane的法向量垂直。
本算法仅仅涉及到共面的直线与线段的相交检测。但为通用起见，首先介绍如何检验三维空间中的两条线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; }
如果线在面上，任意返回一个线上的点即可；如果线面平行，返回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.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; } }
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); } }
标签：
原文地址：http://www.cnblogs.com/luxkodev/p/4781060.html
