首页 > 代码库 > 计算几何---点与三角形面片的关系

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

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

  常用的方法是,如果点在三角形内部,则连接与三角形的三个顶点,并计算组成的组成矢量的角度和,如果为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;

  

  

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