[Bf-blender-cvs] [0bee94f920f] master: Fixed delaunay check, was causing 'desperation' messages.
Howard Trickey
noreply at git.blender.org
Tue Nov 5 19:14:24 CET 2019
Commit: 0bee94f920f956cea6bc41a884c8a9ec96c4a03b
Author: Howard Trickey
Date: Tue Nov 5 13:12:34 2019 -0500
Branches: master
https://developer.blender.org/rB0bee94f920f956cea6bc41a884c8a9ec96c4a03b
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