雲風的 BLOG: 判斷點是否在三角形內的演算法精度問題
阿新 • • 發佈:2018-12-23
我看了一下原始碼,把這個函式提出來,改寫了一點點,方便獨立測試。
#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[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.6193819 }; 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; }
最後這個優化版本可以通過測試。