« Skynet 1.2.0 | 返回首页 | 3d engine 项目招聘 »

判断点是否在三角形内的算法精度问题

今天一个同事反应,在使用 recastnavigation 库时,判断一个点是否在一个三角形内,遇到了精度问题,而且精度误差很大。

具体是 dtClosestHeightPointTriangle 这个函数。

他给出了一组测试参数,abc 三点为 {261.137939, 8.13000488} , {73.6379318, 8.13000488}, {76.9379349, 10.2300053} ,测试 p 为 {74.4069519 , 8.6093819 } 应该在这个三角形内,但是这个函数计算出来并不是。

我看了一下源代码,把这个函数提出来,改写了一点点,方便独立测试。

#include <stdio.h>

// #define float double

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

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

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 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;
}

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

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

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

    return 0;
}

如果你在前面加上 #define float double ,把所有 float 换成双精度,那么测试是可以通过的。

我认为问题出在 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;
}

最后这个优化版本可以通过测试。

Comments

@jj

我写 blog 的时候抄错了一个数字,是 8.6093819 不是 8.6193819 。另外,之前程序的 main 函数里应该是 [3] 。一并改过来了。

这个错误不影响这篇 blog 的其它主要部分(原算法依然有精度问题,修正后的算法可以解决这个问题)。

云风就是强

p 点确实不在三角形中呢,方法在 URL 中,我把 double 换成 long long ,再为你的点先乘上 10^8 再放进去算,因为没有除法,long long 也够长,算出来 p 点不在三角形中。

花式扛精 @旗舰

单精度能解决问题,挺好。为了避免单精度数字带来的误差,我的fortran程序里统统用双精度小数,但是这样会降低计算效率。

@旗舰 这个是修复BUG,不是优化,不要听风就是雨

还是有欠考虑的地方,一要从三角形任意顶点出发,且构造最小面积凸边形,如果能遍历4点一定在三角形外,只能遍历3点,如果3点都是三角形的顶点则在三角形外,否则在三角形内,应该没有bug了吧

按凸边形规则遍历

顺时针遍历,只能遍历三个点就在三角形内,否则在三角形外

好有哲理!!那么,如果写代码的时候就做这个优化算不算过度优化

Post a comment

非这个主题相关的留言请到:留言本