计算几何---点与三角形面片的关系

点是否再三角形面之内的位置关系判定

　　常用的方法是，如果点在三角形内部，则连接与三角形的三个顶点，并计算组成的组成矢量的角度和，如果为2pei， 则在三角形内部，否则不再三角形内不。理论上该方法可行，但是效率非常慢。

　　因为要计算夹角，需要多次计算平方根。

1. 判定是否同边的技术原理（Same Side Technique）

　　快速且简洁的方法如下：

　　如图，三角形，线段AB, BC, and CA，每一条都半分剩下的部分，且剩下的在三角形的外部部分，我们可以利用这个特性做如下处理：

　　如果一个点在三角形内部，则点必然在AB下面， BC左面，AC右边。如果不符合这个特性，则判定结果比在三角形外面，因此这个判定很快。

（1）计算一个点在一个线段的位置关系

　　通过叉积的方式，用[B - A ]和[p - A ]组成的叉积，即AB * Ap, 如果结果为正（叉积朝向屏幕外面），则p 在AB的左边， 否则（叉积朝向屏幕里面）在AB的右边。

　　因为在3D空间中，AB方向线段，与点p的位置关系，可以借助点C来判定，因为ABC组成一个平面，因此C可以作为参考值。

　　继而，判定任意点p， 如果[B - A ]和[p - A ]叉积方向，与[B - A ]和[C - A]的叉积方向不相同，则可直接判定没有在三角形内部； 如果与叉积方向相同，则可以继续对余下的边进行同样的判定方式。

　　判断两个叉积结果的方向是否同向，使用点乘的方式。

　　实现方面的方法非常容易，可以实现一个方法来判定是否在相同的边的位置。在判定的时候可以分别调用：

function SameSide(p1,p2, a,b)
    cp1 = CrossProduct(b-a, p1-a)
    cp2 = CrossProduct(b-a, p2-a)
    if DotProduct(cp1, cp2) >= 0 then return true
    else return false

function PointInTriangle(p, a,b,c)
    if SameSide(p,a, b,c) and SameSide(p,b, a,c)
        and SameSide(p,c, a,b) then return true
    else return false

　　这样的方法非常高效，不需要计算平方根，弧度反求，等等。

2 判定中心的技术原理（Barycentric Technique）

　　还有一种方法可以执行的更快，原理更简单的方法，与数学知识有一点点关联。

　　因为三角形的三个点确定一个平面，在这个平面中，可以任意选择一点，其他的点都是作为相对位置关系。如以点A为例，然后选定[C-A][B-A]两个向量作为该平面的坐标分量，因此三角形平面中任意一点都可以用这连个向量AB，AC表示。

　　即　　P = A + u * (C - A) + v * (B - A)

注意：如果系数u或者v都< 0, 则p点必然在三角形区域之外。同样，如果u 或者v > 1, 则 点 p则在三角形之外，最后，如果 u + v > 1, 同样越过BC边，而到三角形之外

　　因此可以计算出P在u，v系数下的表示。但是怎么获取p点在上述公式中的系数u，v呢？

　　可以使用下述的公式推导：

    // We have two unknowns (u and v) so we need two equations to solve
    // for them.  Dot both sides by v0 to get one and dot both sides by
    // v1 to get a second.
    (v2) . v0 = (u * v0 + v * v1) . v0
    (v2) . v1 = (u * v0 + v * v1) . v1

    // Distribute v0 and v1
    v2 . v0 = u * (v0 . v0) + v * (v1 . v0)
    v2 . v1 = u * (v0 . v1) + v * (v1 . v1)

    // Now we have two equations and two unknowns and can solve one
    // equation for one variable and substitute into the other.  Or
    // if you‘re lazy like me, fire up Mathematica and save yourself
    // some handwriting.
    Solve[v2.v0 == {u(v0.v0) + v(v1.v0), v2.v1 == u(v0.v1) + v(v1.v1)}, {u, v}]
    u = ((v1.v1)(v2.v0)-(v1.v0)(v2.v1)) / ((v0.v0)(v1.v1) - (v0.v1)(v1.v0))
    v = ((v0.v0)(v2.v1)-(v0.v1)(v2.v0)) / ((v0.v0)(v1.v1) - (v0.v1)(v1.v0))

　　下面是用于编程的快速实现：

// Compute vector，计算向量AC， AB， CP
v0 = C - A
v1 = B - A
v2 = P - A

// Compute dot products 计算点乘
dot00 = dot(v0, v0)
dot01 = dot(v0, v1)
dot02 = dot(v0, v2)
dot11 = dot(v1, v1)
dot12 = dot(v1, v2)

// Compute barycentric coordinates 计算系数坐标
invDenom = 1 / (dot00 * dot11 - dot01 * dot01)
u = (dot11 * dot02 - dot01 * dot12) * invDenom
v = (dot00 * dot12 - dot01 * dot02) * invDenom

// Check if point is in triangle 检测
return (u >= 0) && (v >= 0) && (u + v < 1)

endl;