判断点是否在三角形内

参考文章：

https://www.cnblogs.com/graphics/archive/2010/08/05/1793393.html

https://blog.codingnow.com/2018/11/float_precision_problem.html#comment-47738

https://www.zhihu.com/question/26022206

 

重心法：

三角形的三个点在同一个平面上，如果选中其中一个点，其他两个点不过是相对该点的位移而已，比如选择点A作为起点，那么点B相当于在AB方向移动一段距离得到，而点C相当于在AC方向移动一段距离得到。

所以对于平面内任意一点，都可以由如下方程来表示

P = A +  u * (C – A) + v * (B - A) // 方程1

如果系数u或v为负值，那么相当于朝相反的方向移动，即BA或CA方向。那么如果想让P位于三角形ABC内部，u和v必须满足什么条件呢？有如下三个条件

u >= 0

v >= 0

u + v <= 1

几个边界情况，当u = 0且v = 0时，就是点A，当u = 0,v = 1时，就是点B，而当u = 1, v = 0时，就是点C

整理方程1得到P – A = u(C - A) + v(B - A)

令v0 = C – A, v1 = B – A, v2 = P – A，则v2 = u * v0 + v * v1，现在是一个方程，两个未知数，无法解出u和v，将等式两边分别点乘v0和v1的到两个等式

(v2) • v0 = (u * v0 + v * v1) • v0

(v2) • v1 = (u * v0 + v * v1) • v1

注意到这里u和v是数，而v0，v1和v2是向量，所以可以将点积展开得到下面的式子。

v2 • v0 = u * (v0 • v0) + v * (v1 • v0)  // 式1

v2 • v1 = u * (v0 • v1) + v * (v1• v1)   // 式2

解这个方程得到

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))

 

代码如下：

main.cpp

#include <stdio.h>

#define float double

static float dtVdot2D(const float v0[2], const float v1[2])
{
    return v0[0] * v1[0] + v0[1] * v1[1];
}

static float *dtVsub(float p[2], const float v0[2], const float v1[2])
{
    p[0] = v0[0] - v1[0];
    p[1] = v0[1] - v1[1];
    return p;
}

static int dtClosestHeightPointTriangle(
        const float p[2], const float a[2], const float b[2], const float c[2], float *h)
{
    float v0[2], v1[2], v2[2];

    dtVsub(v0, c,a);
    dtVsub(v1, b,a);
    dtVsub(v2, p,a);

    float dot00 = dtVdot2D(v0, v0);
    float dot01 = dtVdot2D(v0, v1);
    float dot02 = dtVdot2D(v0, v2);
    float dot11 = dtVdot2D(v1, v1);
    float dot12 = dtVdot2D(v1, v2);

    // Compute barycentric coordinates
    float InvDenom = 1.0f / (dot00 * dot11 - dot01 * dot01);
    float u = (dot11 * dot02 - dot01 * dot12) * InvDenom;
    float v = (dot00 * dot12 - dot01 * dot02) * InvDenom;

    // The (sloppy) epsilon is needed to allow to get height of points which
    // are interpolated along the edges of the triangles.
    float EPS = 1e-4f;

    // If point lies inside the triangle, return interpolated ycoord.
    if (u >= -EPS && v >= -EPS && (u+v) <= 1+EPS) {
        *h = a[1] + (v0[1]*u + v1[1]*v);
        return 1;
    }
    return 0;
}


/*
static int dtClosestHeightPointTriangle(
        const float p[2], const float a[2], const float b[2],const float c[2], float *h)
{
    float v0[2], v1[2], v2[2];

    dtVsub(v0, c,a);
    dtVsub(v1, b,a);
    dtVsub(v2, p,a);

    float Denom = v0[0] * v1[2] - v0[2] * v1[0];
    float u = v1[2] * v2[0] - v1[0] * v2[2];
    float v = v0[0] * v2[2] - v0[2] * v2[0];

    if (Denom < 0) {
        Denom = -Denom;
        u = -u;
        v = -v;
    }

    float EPS = - 1e-4f * Denom;

    if (u >= EPS && v >= EPS && (u+v) <= Denom - EPS) {
        *h = a[1] + (v0[1]*u + v1[1]*v) / Denom;
        return 1;
    }
    return 0;
}
*/

int main()
{
    float a[2] = {261.137939, 8.13000488};
    float b[2] = {73.6379318, 8.13000488};
    float c[2] = {76.9379349, 10.2300053};
    float p[2] = {74.4069519, 8.61938190};
    float h;

    int r = dtClosestHeightPointTriangle(p, a, b, c, &h);

    printf("%d %f\n", r, h);

    return 0;
}

运行结果如下：

 

单精度与双精度是什么意思，有什么区别：

 

 

精度优化：

我认为问题出在 dot00 * dot11 - dot01 * dot01 这样的运算上。dot00 点乘已经是单个量的平方，在测试数据中，大约这个量会是 261 - 73 = 188 ，小数点前大约是 8bit 的信息含量，如果我们计算 dot00 * dot11 ，差不多会得到一个这个量的 4 次方的结果，也就是 28bit ~ 32bit 之间。

但是 float 本身的有效精度才 23bit ，对一个 2^32 的数字做加减法，本身的误差就可能在 2 ~ 2^9 左右，这个误差是相当巨大的。

这段程序一个明显可以改进的地方是把乘 InvDenom 从 u v 中去掉，但 Denom 看起来可能是负数，需要增加符号判断。那么代码应该写成：

    float Denom = (dot00 * dot11 - dot01 * dot01);
    float u = (dot11 * dot02 - dot01 * dot12);
    float v = (dot00 * dot12 - dot01 * dot02);

    if (Denom < 0) {
        Denom = -Denom;
        u = -u;
        v = -v;
    }
    float EPS = 1e-4f * Denom ;

    // If point lies inside the triangle, return interpolated ycoord.
    if (u >= -EPS && v >= -EPS && (u+v) <= Denom+EPS) {
        *h = a[1] + (v0[1]*u + v1[1]*v) / Denom;
        return 1;
    }      

光这样写还是不够，其实我们应该进一步把 dot00 * dot11 - dot01 * dot01 展开为 (v0[0] * v1[1] - v0[1] * v1[0]) * (v0[0] * v1[1] - v1[0] * v0[1]) 。这样，就不会在四次方的基础上再做加减法，而是在二次方的基础上先做加减，再做乘法。这样就最大化的保持了精度。

由于简化过后，发现 Denom 是个平方数，一定为正，所以可以去掉符号判断。我优化过的函数是这样的：

static int
dtClosestHeightPointTriangle(const float p[3], const float a[3], const float b[3],const float c[3], float *h) {
    float v0[3], v1[3], v2[3];

    dtVsub(v0, c,a);
    dtVsub(v1, b,a);
    dtVsub(v2, p,a);

    float Denom = (v0[0] * v1[2] - v0[2] * v1[0]) * (v0[0] * v1[2] - v1[0] * v0[2]);
    float u = (v1[0] * v2[2] - v1[2] * v2[0]) * (v1[0] * v0[2] - v0[0] * v1[2]);
    float v = (v0[0] * v2[2] - v0[2] * v2[0]) * (v0[0] * v1[2] - v1[0] * v0[2]);

    float EPS = - 1e-4f * Denom;

    if (u >= EPS && v >= EPS && (u+v) <= Denom - EPS) {
        *h = a[1] + (v0[1]*u + v1[1]*v) / Denom;
        return 1;
    }
    return 0;
}

由于 u,v,denom 都有共同的项 (v0[0] * v1[2] - v1[0] * v0[2]) 可以约掉，能进一步保留精度。但需要处理一下符号问题。所以最终版本是这样的：

static int
dtClosestHeightPointTriangle(const float p[3], const float a[3], const float b[3],const float c[3], float *h) {
    float v0[3], v1[3], v2[3];

    dtVsub(v0, c,a);
    dtVsub(v1, b,a);
    dtVsub(v2, p,a);

    float Denom = v0[0] * v1[2] - v0[2] * v1[0];
    float u = v1[2] * v2[0] - v1[0] * v2[2];
    float v = v0[0] * v2[2] - v0[2] * v2[0];

    if (Denom < 0) {
        Denom = -Denom;
        u = -u;
        v = -v;
    }

    float EPS = - 1e-4f * Denom;

    if (u >= EPS && v >= EPS && (u+v) <= Denom - EPS) {
        *h = a[1] + (v0[1]*u + v1[1]*v) / Denom;
        return 1;
    }
    return 0;
}

最后这个优化版本即为代码中的注释部分。