[Bf-blender-cvs] [9d1031b0119] blender-v2.81-release: Fixed delaunay check, was causing 'desperation' messages.

Howard Trickey noreply at git.blender.org
Tue Nov 5 19:27:11 CET 2019


Commit: 9d1031b011994f508c766593eb8b149077ebe84a
Author: Howard Trickey
Date:   Tue Nov 5 13:23:20 2019 -0500
Branches: blender-v2.81-release
https://developer.blender.org/rB9d1031b011994f508c766593eb8b149077ebe84a

Fixed delaunay check, was causing 'desperation' messages.

Check was losing precision -- adjust by translating points
before calculating circumcircle.
Also, needed to check for flippability of edges before flipping.

===================================================================

M	source/blender/blenlib/intern/delaunay_2d.c

===================================================================

diff --git a/source/blender/blenlib/intern/delaunay_2d.c b/source/blender/blenlib/intern/delaunay_2d.c
index af4fa9fa54e..f297d8cfeed 100644
--- a/source/blender/blenlib/intern/delaunay_2d.c
+++ b/source/blender/blenlib/intern/delaunay_2d.c
@@ -124,6 +124,17 @@ static int CCW_test(const double a[2], const double b[2], const double c[2], con
 
   /* This is twice the signed area of triangle abc. */
   det = (b[0] - a[0]) * (c[1] - a[1]) - (c[0] - a[0]) * (b[1] - a[1]);
+  if (eps == 0.0) {
+    if (det > 0) {
+      return 1;
+    }
+    else if (det < 0) {
+      return -1;
+    }
+    else {
+      return 0;
+    }
+  }
   ab = len_v2v2_db(a, b);
   if (ab <= eps) {
     return 0;
@@ -617,8 +628,9 @@ static bool locate_point_final(const double p[2],
   double lambda, close[2];
   bool done = false;
 #ifdef DEBUG_CDT
-  int dbglevel = 0;
-  if (dbglevel > 0) {
+  int dbg_level = 0;
+
+  if (dbg_level > 0) {
     fprintf(stderr, "locate_point_final %d\n", try_neighbors);
     dump_se(tri_se, "tri_se");
     fprintf(stderr, "\n");
@@ -628,7 +640,7 @@ static bool locate_point_final(const double p[2],
   i = 0;
   do {
 #ifdef DEBUG_CDT
-    if (dbglevel > 1) {
+    if (dbg_level > 1) {
       fprintf(stderr, "%d: ", i);
       dump_se(se, "search se");
     }
@@ -640,7 +652,7 @@ static bool locate_point_final(const double p[2],
     if (len_close_p < epsilon) {
       if (len_v2v2_db(p, a) < epsilon) {
 #ifdef DEBUG_CDT
-        if (dbglevel > 0) {
+        if (dbg_level > 0) {
           fprintf(stderr, "OnVert case a (%.2f,%.2f)\n", F2(a));
         }
 #endif
@@ -651,7 +663,7 @@ static bool locate_point_final(const double p[2],
       }
       else if (len_v2v2_db(p, b) < epsilon) {
 #ifdef DEBUG_CDT
-        if (dbglevel > 0) {
+        if (dbg_level > 0) {
           fprintf(stderr, "OnVert case b (%.2f,%.2f)\n", F2(b));
         }
 #endif
@@ -662,7 +674,7 @@ static bool locate_point_final(const double p[2],
       }
       else if (lambda > 0.0 && lambda < 1.0) {
 #ifdef DEBUG_CDT
-        if (dbglevel > 0) {
+        if (dbg_level > 0) {
           fprintf(stderr, "OnEdge case, lambda=%f\n", lambda);
           dump_se(se, "se");
         }
@@ -682,7 +694,7 @@ static bool locate_point_final(const double p[2],
   } while (se != tri_se && !done);
   if (!done) {
 #ifdef DEBUG_CDT
-    if (dbglevel > 1) {
+    if (dbg_level > 1) {
       fprintf(stderr,
               "not done, dist_inside=%f %f %f\n",
               dist_inside[0],
@@ -692,7 +704,7 @@ static bool locate_point_final(const double p[2],
 #endif
     if (dist_inside[0] >= 0.0 && dist_inside[1] >= 0.0 && dist_inside[2] >= 0.0) {
 #ifdef DEBUG_CDT
-      if (dbglevel > 0) {
+      if (dbg_level > 0) {
         fprintf(stderr, "InFace case\n");
         dump_se_cycle(tri_se, "tri", 10);
       }
@@ -740,14 +752,14 @@ static bool locate_point_final(const double p[2],
           r_lr->edge_lambda = lambda;
         }
 #ifdef DEBUG_CDT
-        if (dbglevel > 0) {
+        if (dbg_level > 0) {
           fprintf(
               stderr, "desperation case kind=%u lambda=%f\n", r_lr->loc_kind, r_lr->edge_lambda);
           dump_se(r_lr->se, "se");
           BLI_assert(0); /* While developing, catch these "should not happens" */
         }
 #endif
-        fprintf(stderr, "desperation!\n");  // TODO: remove
+        fprintf(stderr, "desperation! point=(%g,%g)\n", p[0], p[1]);  // TODO: remove
         return true;
       }
     }
@@ -769,9 +781,9 @@ static LocateResult locate_point(CDT_state *cdt, const double p[2])
   int visit = ++cdt->visit_count;
   int loop_count = 0;
 #ifdef DEBUG_CDT
-  int dbglevel = 0;
+  int dbg_level = 0;
 
-  if (dbglevel > 0) {
+  if (dbg_level > 0) {
     fprintf(stderr, "locate_point (%.2f,%.2f), visit_index=%d\n", F2(p), visit);
   }
 #endif
@@ -790,7 +802,7 @@ static LocateResult locate_point(CDT_state *cdt, const double p[2])
     v = cdt->vert_array[i];
     dist_squared = len_squared_v2v2_db(p, v->co);
 #ifdef DEBUG_CDT
-    if (dbglevel > 0) {
+    if (dbg_level > 0) {
       fprintf(stderr, "try start vert %d, dist_squared=%f\n", i, dist_squared);
       dump_v(v, "v");
     }
@@ -806,7 +818,7 @@ static LocateResult locate_point(CDT_state *cdt, const double p[2])
     BLI_assert(cur_se->face != cdt->outer_face);
   }
 #ifdef DEBUG_CDT
-  if (dbglevel > 0) {
+  if (dbg_level > 0) {
     dump_se(cur_se, "start vert edge");
   }
 #endif
@@ -815,7 +827,7 @@ static LocateResult locate_point(CDT_state *cdt, const double p[2])
     /* Find edge of cur_tri that separates p and t's centroid,
      * and where other tri over the edge is unvisited. */
 #ifdef DEBUG_CDT
-    if (dbglevel > 0) {
+    if (dbg_level > 0) {
       dump_se_cycle(cur_se, "cur search face", 5);
     }
 #endif
@@ -826,10 +838,10 @@ static LocateResult locate_point(CDT_state *cdt, const double p[2])
     a = cur_se->vert->co;
     b = cur_se->next->vert->co;
     c = cur_se->next->next->vert->co;
-    if (CCW_test(a, b, p, epsilon) >= 0 && CCW_test(b, c, p, epsilon) >= 0 &&
-        CCW_test(c, a, p, epsilon) >= 0) {
+    if (CCW_test(a, b, p, 0.0) >= 0 && CCW_test(b, c, p, 0.0) >= 0 &&
+        CCW_test(c, a, p, 0.0) >= 0) {
 #ifdef DEBUG_CDT
-      if (dbglevel > 1) {
+      if (dbg_level > 1) {
         fprintf(stderr, "p in current triangle\n");
       }
 #endif
@@ -844,28 +856,28 @@ static LocateResult locate_point(CDT_state *cdt, const double p[2])
       b = next_se->next->vert->co;
       c = next_se->next->next->vert->co;
 #ifdef DEBUG_CDT
-      if (dbglevel > 1) {
+      if (dbg_level > 1) {
         dump_se(next_se, "search edge");
-        fprintf(stderr, "tri centroid=(%.2f,%.2f)\n", F2(cur_tri->centroid));
+        fprintf(stderr, "tri centroid=(%.3f,%.3f)\n", F2(cur_tri->centroid));
         validate_face_centroid(next_se);
       }
 #endif
       next_se_sym = sym(next_se);
-      if (CCW_test(a, b, p, epsilon) <= 0 && next_se->face != cdt->outer_face) {
+      if (CCW_test(a, b, p, 0.0) <= 0 && next_se->face != cdt->outer_face) {
 #ifdef DEBUG_CDT
-        if (dbglevel > 1) {
+        if (dbg_level > 1) {
           fprintf(stderr, "CCW_test(a, b, p) <= 0\n");
         }
 #endif
 #ifdef DEBUG_CDT
-        if (dbglevel > 0) {
+        if (dbg_level > 0) {
           dump_se(next_se_sym, "next_se_sym");
           fprintf(stderr, "next_se_sym face visit=%d\n", next_se_sym->face->visit_index);
         }
 #endif
         if (next_se_sym->face->visit_index != visit) {
 #ifdef DEBUG_CDT
-          if (dbglevel > 0) {
+          if (dbg_level > 0) {
             fprintf(stderr, "found edge to cross\n");
           }
 #endif
@@ -890,32 +902,50 @@ static LocateResult locate_point(CDT_state *cdt, const double p[2])
   return lr;
 }
 
-/** return true if circumcircle(v1, v2, v3) does not contain p. */
+/** Return true if circumcircle(v1, v2, v3) does not contain p.
+ * To avoid possible infinite flip loops, we will say true even if p is inside the circle
+ * but less than epsilon from the boundary; or if v1, v2, v3, form a straight line.
+ */
 static bool delaunay_check(CDTVert *v1, CDTVert *v2, CDTVert *v3, CDTVert *p, const double epsilon)
 {
-  double a, b, c, d, z1, z2, z3;
-  const double *p1, *p2, *p3;
-  double cen[2], r, len_pc;
-  /* To do epislon test, need center and radius of circumcircle. */
-  p1 = v1->co;
-  p2 = v2->co;
-  p3 = v3->co;
-  z1 = dot_v2v2_db(p1, p1);
-  z2 = dot_v2v2_db(p2, p2);
-  z3 = dot_v2v2_db(p3, p3);
-  a = p1[0] * (p2[1] - p3[1]) - p1[1] * (p2[0] - p3[0]) + p2[0] * p3[1] - p3[0] * p2[1];
-  b = z1 * (p3[1] - p2[1]) + z2 * (p1[1] - p3[1]) + z3 * (p2[1] - p1[1]);
-  c = z1 * (p2[0] - p3[0]) + z2 * (p3[0] - p1[0]) + z3 * (p1[0] - p2[0]);
-  d = z1 * (p3[0] * p2[1] - p2[0] * p3[1]) + z2 * (p1[0] * p3[1] - p3[0] * p1[1]) +
-      z3 * (p2[0] * p1[1] - p1[0] * p2[1]);
-  if (a == 0.0) {
-    return true; /* Not really, but this shouldn't happen. */
-  }
-  cen[0] = -b / (2 * a);
-  cen[1] = -c / (2 * a);
-  r = sqrt((b * b + c * c - 4 * a * d) / (4 * a * a));
-  len_pc = len_v2v2_db(p->co, cen);
-  return (len_pc >= (r - epsilon));
+  double x1, y1, x2, y2, den, cenx, ceny, rad, pc, a, b, w, z, q, s;
+
+  /* Find center and radius of circumcircle of v1,v2,v3.
+   * Transform coords so v3 is at origin to help reduce floating point error
+   * when coordinates are far from (0,0) but close together.
+   */
+  x1 = v1->co[0] - v3->co[0];
+  y1 = v1->co[1] - v3->co[1];
+  x2 = v2->co[0] - v3->co[0];
+  y2 = v2->co[1] - v3->co[1];
+  den = 2.0 * (x1 * y2 - x2 * y1);
+  if (UNLIKELY(den == 0.0)) {
+    /* v1, v2, v3 are in a line. */
+    return true;
+  }
+  /* cen[0] = det(x1**2 + y1**2, y1, x2**2 + y2**2, y2) / den
+   * cen[1] = det(x1, x1**2 + y1**2, x2, x2**2 + y2**2) / den
+   * den = 2 * det(x1, y1, x2, y2)
+   */
+  a = x1 * x1 + y1 * y1;
+  b = x2 * x2 + y2 * y2;
+  cenx = (a * y2 - b * y1) / den;
+  ceny = (x1 * b - x2 * a) / den;
+  w = x1 - cenx;
+  z = y1 - ceny;
+  rad = sqrt(w * w + z * z);
+  q = p->co[0] - v3->co[0] - cenx;
+  s = p->co[1] - v3->co[1] - ceny;
+  pc = sqrt(q * q + s * s);
+  return (pc >= rad - epsilon);
+}
+
+/* Return true if we can flip edge v1-v3 to edge v2-v4 inside quad v1v2v3v4 (in CCW order).
+ * We can do this if angles v4-v1-v2 and v2-v3-v4 are both CCW or straight.
+ */
+static inline bool can_flip(CDTVert *v1, CDTVert *v2, CDTVert *v3, CDTVert *v4)
+{
+  return CCW_test(v4->co, v1->co, v2->co, 0.0) >= 0 && CCW_test(v2->co, v3->co, v4->co, 0.0) >= 0;
 }
 
 /** Use LinkNode linked list as stack of SymEdges, allocating from cdt->listpool. */
@@ -955,12 +985,12 @@ static void flip(SymEdge *se, CDT_state *cdt)
   CDTFace *t1, *t2;
   CDTVert *v1, *v2;
 #ifdef DEBUG_CDT
-  const int dbglevel = 0;
+  const int dbg_level = 0;
 #endif
 
   sesym = sym(se);
 #ifdef DEBUG_CDT
-  if (dbglevel > 0) {
+  if (dbg_level > 0) {
     fprintf(stderr, "flip\n");
     dump_se(se, "se");
     dump_se(sesym, "sesym");
@@ -975,7 +1005,7 @@ static void flip(SymEdge *se, CDT_state *cdt)
   csym = sym(c);
   dsym = sym(d);
 #ifdef DEBUG_CDT
-  if (dbglevel > 1) {
+  if (dbg_level > 1) {
     dump_se(a, "a");
     dump_se(b, "b");
     dump_se(c, "c");
@@ -1020,7 +1050,7 @@ static void flip(SymEdge *se, CDT_state *cdt)
   calc_face_cent

@@ Diff output truncated at 10240 characters. @@



More information about the Bf-blender-cvs mailing list