0
0
0
---
判断点是否在三角形内的算法精度问题
云风 发表于 2018年11月09日 10:36 | Hits: 1769
今天一个同事反应,在使用recastnavigation库时,判断一个点是否在一个三角形内,遇到了精度问题,而且精度误差很大。
具体是dtClosestHeightPointTriangle这个函数。
他给出了一组测试参数,abc 三点为 {261.137939, 8.13000488} , {73.6379318, 8.13000488}, {76.9379349, 10.2300053} ,测试 p 为 {74.4069519 , 8.6193819 } 应该在这个三角形内,但是这个函数计算出来并不是。
我看了一下源代码,把这个函数提出来,改写了一点点,方便独立测试。
#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) 如果你在前面加上 #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 = -EPS && v >= -EPS && (u+v) 光这样写还是不够,其实我们应该进一步把 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) 由于 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 = EPS && v >= EPS && (u+v) 最后这个优化版本可以通过测试。
原文链接: https://blog.codingnow.com/2018/11/float_precision_problem.html
我要给这篇文章打分:
可以不填写评论, 而只是打分. 如果发表评论, 你可以给的分值是-5到+5, 否则, 你只能评-1, +1两种分数. 你的评论可能需要审核.