From a0892bb690fec36edd7999d6729e0826d72dab86 Mon Sep 17 00:00:00 2001 From: Howard Trickey Date: Sat, 21 Dec 2019 12:23:02 -0500 Subject: Fix crash in delaunay triangulation due to epsilon issues. --- source/blender/blenlib/intern/delaunay_2d.c | 586 +++++++++++++++++++-------- tests/gtests/blenlib/BLI_delaunay_2d_test.cc | 272 +++++++++++-- 2 files changed, 647 insertions(+), 211 deletions(-) diff --git a/source/blender/blenlib/intern/delaunay_2d.c b/source/blender/blenlib/intern/delaunay_2d.c index 27f7a2a4430..350f99ce434 100644 --- a/source/blender/blenlib/intern/delaunay_2d.c +++ b/source/blender/blenlib/intern/delaunay_2d.c @@ -53,6 +53,7 @@ typedef struct CDTVert { SymEdge *symedge; /* Some edge attached to it. */ LinkNode *input_ids; /* List of corresponding vertex input ids. */ int index; /* Index into array that cdt keeps. */ + int visit_index; /* Which visit epoch has this been seen. */ } CDTVert; typedef struct CDTEdge { @@ -110,7 +111,9 @@ static void dump_v(const CDTVert *v, const char *lab); static void dump_se_cycle(const SymEdge *se, const char *lab, const int limit); static void dump_id_list(const LinkNode *id_list, const char *lab); static void dump_cdt(const CDT_state *cdt, const char *lab); +static void dump_cdt_vert_neighborhood(CDT_state *cdt, int v, int maxdist, const char *lab); static void cdt_draw(CDT_state *cdt, const char *lab); +static void write_cdt_input_to_file(const CDT_input *inp); static void validate_face_centroid(SymEdge *se); static void validate_cdt(CDT_state *cdt, bool check_all_tris); #endif @@ -151,14 +154,11 @@ static int CCW_test(const double a[2], const double b[2], const double c[2], con return 0; } -/** return true if a -- b -- c are in that order, assuming they are on a straight line. */ -static bool in_line(const double a[2], const double b[2], const double c[2]) +/** return true if a -- b -- c are in that order, assuming they are on a straight line according to + * CCW_test. */ +static bool in_line(const double a[2], const double b[2], const double c[2], double eps) { - double dir_ab[2], dir_ac[2]; - - sub_v2_v2v2_db(dir_ab, a, b); - sub_v2_v2v2_db(dir_ac, a, c); - return dot_v2v2_db(dir_ab, dir_ac) >= 0.0; + return fabs(len_v2v2_db(a, c) - (len_v2v2_db(a, b) + len_v2v2_db(b, c))) <= eps; } #ifndef NDEBUG @@ -208,6 +208,7 @@ static CDTVert *add_cdtvert(CDT_state *cdt, double x, double y) } BLI_assert(cdt->vert_array_len < cdt->vert_array_len_alloc); v->index = cdt->vert_array_len; + v->visit_index = 0; cdt->vert_array[cdt->vert_array_len++] = v; return v; } @@ -365,8 +366,8 @@ static CDTEdge *add_diagonal(CDT_state *cdt, SymEdge *s1, SymEdge *s2) CDTFace *fold, *fnew; SymEdge *sdiag, *sdiagsym; SymEdge *s1prev, *s1prevsym, *s2prev, *s2prevsym, *se; - BLI_assert(reachable(s1, s2, 20)); - BLI_assert(reachable(s2, s1, 20)); + BLI_assert(reachable(s1, s2, 2000)); + BLI_assert(reachable(s2, s1, 2000)); fold = s1->face; fnew = add_cdtface(cdt); s1prev = prev(s1); @@ -385,7 +386,7 @@ static CDTEdge *add_diagonal(CDT_state *cdt, SymEdge *s1, SymEdge *s2) s2->rot = sdiagsym; sdiagsym->rot = s2prevsym; #ifdef DEBUG_CDT - BLI_assert(reachable(s2, sdiag, 20)); + BLI_assert(reachable(s2, sdiag, 2000)); #endif for (se = s2; se != sdiag; se = se->next) { se->face = fnew; @@ -1505,7 +1506,8 @@ static void add_edge_constraint( SymEdge *se; CDTEdge *edge; int ccw1, ccw2, isect; - int i, search_count; + int i, search_count, visit; + double curco[2]; double lambda; const double epsilon = cdt->epsilon; bool done, state_through_vert; @@ -1553,7 +1555,7 @@ static void add_edge_constraint( * in = a * out = b * - * crossings[0] will have in = NULL, and crossings[last] will have out = NULL + * crossings[0] will have in = NULL, and crossings[last] will have out = NULL. */ if (r_edges) { *r_edges = NULL; @@ -1569,6 +1571,9 @@ static void add_edge_constraint( dump_se(vse1, " se1"); dump_se(vse2, " se2"); } + if (dbg_level > 2) { + dump_cdt(cdt, "before insert_segment"); + } } #endif if (v1 == v2) { @@ -1579,181 +1584,292 @@ static void add_edge_constraint( #endif return; } - state_through_vert = true; + /* Quick check for special case where segment is one of edges from v1. */ + tstart = t = vse1; done = false; - t = vse1; - search_count = 0; - while (!done) { - /* Invariant: crossings[0 .. BLI_array_len(crossings)] has crossing info for path up to - * but not including the crossing of edge t, which will either be through a vert - * (if state_through_vert is true) or through edge t not at either end. - * In the latter case, t->face is the face that ray v1--v2 goes through after path-so-far. - */ -#ifdef DEBUG_CDT - if (dbg_level > 1) { - fprintf( - stderr, "top of insert_segment main loop, state_through_vert=%d\n", state_through_vert); - dump_se_cycle(t, "current t ", 4); - } -#endif - if (state_through_vert) { - /* Invariant: ray v1--v2 contains t->vert. */ - cdata.in = (BLI_array_len(crossings) == 0) ? NULL : t; - cdata.out = NULL; /* To be filled in if this isn't final. */ + do { + if (t->next->vert == v2) { + cdata.in = NULL; + cdata.out = t->next; cdata.lambda = 0.0; cdata.vert = t->vert; BLI_array_append(crossings, cdata); - if (t->vert == v2) { + cdata.in = t->next->rot; + cdata.out = NULL; + cdata.lambda = 0.0; + cdata.vert = v2; + BLI_array_append(crossings, cdata); #ifdef DEBUG_CDT - if (dbg_level > 0) { - fprintf(stderr, "found v2, so done\n"); - } + if (dbg_level > 1) { + fprintf(stderr, "special one segment case\n"); + dump_se(t, " "); + } #endif - done = true; + done = true; + break; + } + t = t->rot; + } while (t != tstart); + if (!done) { + /* To prevent infinite loop in the face of epsilon tests that might lead us back to + * an already-visited (vertex, face) pair, use visit indices. + */ + visit = ++cdt->visit_count; + state_through_vert = true; + done = false; + t = vse1; + search_count = 0; + while (!done) { + /* Invariant: crossings[0 .. BLI_array_len(crossings)] has crossing info for path up to + * but not including the crossing of edge t, which will either be through a vert + * (if state_through_vert is true) or through edge t not at either end. + * In the latter case, t->face is the face that ray v1--v2 goes through after path-so-far. + * Rather than continually looking for intersection of v1--v2, we keep track of + * last vertex or intersection point in curco,because that may be slightly off the ray + * v1--v2. + */ +#ifdef DEBUG_CDT + if (dbg_level > 1) { + fprintf(stderr, + "top of insert_segment main loop, state_through_vert=%d\n", + state_through_vert); + dump_se_cycle(t, "current t ", 4); } - else { - /* Do ccw scan of triangles around t->vert to find exit triangle for ray v1--v2. */ - tstart = t; - tout = NULL; - do { - va = t->next->vert; - vb = t->next->next->vert; - ccw1 = CCW_test(t->vert->co, va->co, v2->co, epsilon); - ccw2 = CCW_test(t->vert->co, vb->co, v2->co, epsilon); +#endif + BLI_assert(t->vert->visit_index != visit || t->face->visit_index != visit); + t->vert->visit_index = visit; + t->face->visit_index = visit; + if (state_through_vert) { + /* Invariant: ray vcur--v2 contains t->vert. */ + cdata.in = (BLI_array_len(crossings) == 0) ? NULL : t; + cdata.out = NULL; /* To be filled in if this isn't final. */ + cdata.lambda = 0.0; + cdata.vert = t->vert; + BLI_array_append(crossings, cdata); + if (t->vert == v2) { #ifdef DEBUG_CDT - if (dbg_level > 1) { - fprintf(stderr, "non-final through vert case\n"); - dump_v(va, " va"); - dump_v(vb, " vb"); - fprintf(stderr, "ccw1=%d, ccw2=%d\n", ccw1, ccw2); + if (dbg_level > 0) { + fprintf(stderr, "found v2, so done\n"); } #endif - if (ccw1 == 0 && in_line(t->vert->co, va->co, v2->co)) { + done = true; + } + else { + /* Do ccw scan of triangles around t->vert to find exit triangle for ray vcur--v2. */ + tstart = t; + tout = NULL; + do { + va = t->next->vert; + vb = t->next->next->vert; + ccw1 = CCW_test(t->vert->co, va->co, v2->co, epsilon); + ccw2 = CCW_test(t->vert->co, vb->co, v2->co, epsilon); +#ifdef DEBUG_CDT + if (dbg_level > 1) { + fprintf(stderr, "non-final through vert case\n"); + dump_v(va, " va"); + dump_v(vb, " vb"); + fprintf(stderr, "ccw1=%d, ccw2=%d\n", ccw1, ccw2); + } +#endif + if (ccw1 == 0 && in_line(t->vert->co, va->co, v2->co, epsilon) && + va->visit_index != visit) { #ifdef DEBUG_CDT - if (dbg_level > 0) { - fprintf(stderr, "ray goes through va\n"); + if (dbg_level > 0) { + fprintf(stderr, "ray goes through va\n"); + } +#endif + state_through_vert = true; + tout = t; + t = t->next; + break; } + else if (ccw2 == 0 && in_line(t->vert->co, vb->co, v2->co, epsilon) && + vb->visit_index != visit) { +#ifdef DEBUG_CDT + if (dbg_level > 0) { + fprintf(stderr, "ray goes through vb\n"); + } #endif - state_through_vert = true; - tout = t; - t = t->next; - break; - } - else if (ccw2 == 0 && in_line(t->vert->co, vb->co, v2->co)) { + state_through_vert = true; + t = t->next->next; + tout = sym(t); + break; + } + else if (ccw1 > 0 && ccw2 < 0 && + (t->next->vert->visit_index != visit || + t->next->face->visit_index != visit)) { +#ifdef DEBUG_CDT + if (dbg_level > 0) { + fprintf(stderr, "segment intersects\n"); + } +#endif + state_through_vert = false; + tout = t; + t = t->next; + break; + } + t = t->rot; #ifdef DEBUG_CDT - if (dbg_level > 0) { - fprintf(stderr, "ray goes through vb\n"); + if (dbg_level > 1) { + dump_se_cycle(t, "next rot tri", 4); } #endif - state_through_vert = true; - t = t->next->next; - tout = sym(t); - break; - } - else if (ccw1 > 0 && ccw2 < 0) { + } while (t != tstart); + if (tout == NULL) { + /* With exact arithmetic this shouldn't happen, but maybe the epsilon tests made it so + * that we want to go back to a previous vertex. + * As desperation measure, pick unvisited vertex that is closest in line with + * destination. + */ #ifdef DEBUG_CDT if (dbg_level > 0) { - fprintf(stderr, "segment intersects\n"); + fprintf(stderr, "add_edge_constraint desperation search\n"); } #endif - state_through_vert = false; - tout = t; - t = t->next; - break; - } - t = t->rot; + SymEdge *bestt = NULL; + double dot, bestdot = -2.0; + double dir_tv_v2[2], dir_tvnext_v2[2]; + sub_v2_v2v2_db(dir_tv_v2, v2->co, t->vert->co); + do { + if (t->next->vert->visit_index != visit) { + sub_v2_v2v2_db(dir_tvnext_v2, v2->co, t->next->vert->co); + dot = dot_v2v2_db(dir_tv_v2, dir_tvnext_v2); + if (dot > bestdot) { + bestdot = dot; + bestt = t->next; + } + } + t = t->rot; + } while (t != tstart); + if (bestt == NULL) { + /* No unvisited place to go! Give up on adding this edge constraint. */ +#ifdef DEBUG_CDT + fprintf(stderr, "could not add edge constraint\n"); +#endif + return; + } + else { #ifdef DEBUG_CDT - if (dbg_level > 1) { - dump_se_cycle(t, "next rot tri", 4); - } + if (dbg_level > 0) { + fprintf(stderr, "add_edge_constraint desperation search chose to go through\n"); + dump_v(bestt->vert, "desperation vert"); + } #endif - } while (t != tstart); - BLI_assert(tout != NULL); /* TODO: something sensible for "this can't happen" */ - crossings[BLI_array_len(crossings) - 1].out = tout; - } - } - else { /* State is "through edge", not "through vert" */ - /* Invariant: ray v1--v2 intersects segment t->edge, not at either end. - * and t->face is the face we have just passed through. */ - va = t->vert; - vb = t->next->vert; -#ifdef DEBUG_CDT - if (dbg_level > 1) { - fprintf(stderr, "through edge case\n"); - dump_v(va, " va"); - dump_v(vb, " vb"); + tout = bestt; + t = t->next; + } + } + crossings[BLI_array_len(crossings) - 1].out = tout; + } } -#endif - isect = isect_seg_seg_v2_lambda_mu_db(va->co, vb->co, v1->co, v2->co, &lambda, NULL); - /* TODO: something sensible for "this can't happen" */ - BLI_assert(isect == ISECT_LINE_LINE_CROSS); - UNUSED_VARS_NDEBUG(isect); + else { /* State is "through edge", not "through vert" */ + /* Invariant: ray v1--v2 intersects segment t->edge, not at either end. + * and t->face is the face we have just passed through. + * Whatever we make t next should not have both vert and face visited. + */ + va = t->vert; + vb = t->next->vert; + /* Get curco; cdata should have data from last time through the loop still. */ + if (cdata.lambda == 0.0) { + copy_v3_v3_db(curco, cdata.vert->co); + } + else { + interp_v2_v2v2_db(curco, cdata.in->vert->co, cdata.in->next->vert->co, lambda); + } #ifdef DEBUG_CDT - if (dbg_level > 0) { - fprintf(stderr, "intersect point at %f along va--vb\n", lambda); - if (dbg_level == 1) { + if (dbg_level > 1) { + fprintf(stderr, "through edge case, curco=(%f,%f)\n", F2(curco)); dump_v(va, " va"); dump_v(vb, " vb"); } - } #endif - tout = sym(t)->next; - cdata.in = t; - cdata.out = tout; - cdata.lambda = lambda; - cdata.vert = NULL; /* To be filled in with edge split vertex later. */ - BLI_array_append(crossings, cdata); + isect = isect_seg_seg_v2_lambda_mu_db(va->co, vb->co, curco, v2->co, &lambda, NULL); + if (isect != ISECT_LINE_LINE_CROSS) { + /* Shouldn't happen. Just pick something. */ #ifdef DEBUG_CDT - if (dbg_level > 0) { - dump_se_cycle(tout, "next search tri", 4); - } + if (dbg_level > 1) { + fprintf(stderr, "add_edge_constraint no intersect found, using lambda = 0.5\n"); + } #endif - /* 'tout' is 'symedge' from 'vb' to third vertex, 'vc'. */ - BLI_assert(tout->vert == va); - vc = tout->next->vert; - ccw1 = CCW_test(v1->co, v2->co, vc->co, epsilon); + lambda = 0.5; + } #ifdef DEBUG_CDT - if (dbg_level > 1) { - fprintf(stderr, "now searching with third vertex "); - dump_v(vc, "vc"); - fprintf(stderr, "ccw(v1, v2, vc) = %d\n", ccw1); - } + if (dbg_level > 0) { + fprintf(stderr, "intersect point at %f along va--vb\n", lambda); + if (dbg_level == 1) { + dump_v(va, " va"); + dump_v(vb, " vb"); + } + } #endif - if (ccw1 == -1) { - /* v1--v2 should intersect vb--vc. */ + tout = sym(t)->next; + cdata.in = t; + cdata.out = tout; + cdata.lambda = lambda; + cdata.vert = NULL; /* To be filled in with edge split vertex later. */ + BLI_array_append(crossings, cdata); #ifdef DEBUG_CDT - if (dbg_level > 1) { - fprintf(stderr, "v1--v2 intersects vb--vc\n"); + if (dbg_level > 0) { + dump_se_cycle(tout, "next search tri", 4); + dump_se(tout, "tout"); } #endif - t = tout->next; - state_through_vert = false; - } - else if (ccw1 == 1) { - /* v1--v2 should intersect va--vc. */ + /* 'tout' is 'symedge' from 'va' to third vertex, 'vc'. */ + BLI_assert(tout->vert == va); + vc = tout->next->vert; + ccw1 = CCW_test(curco, v2->co, vc->co, epsilon); #ifdef DEBUG_CDT if (dbg_level > 1) { - fprintf(stderr, "v1--v2 intersects va--vc\n"); + fprintf(stderr, "now searching with third vertex "); + dump_v(vc, "vc"); + fprintf(stderr, "ccw(vcur, v2, vc) = %d\n", ccw1); } #endif - t = tout; - state_through_vert = false; - } - else { - /* ccw1 == 0. */ + if (ccw1 == -1 && (vc->visit_index != visit || tout->next->face->visit_index != visit)) { + /* vcur--v2 should intersect vb--vc. */ #ifdef DEBUG_CDT - if (dbg_level > 1) { - fprintf(stderr, "ccw==0 case, so going through or to vc\n"); + if (dbg_level > 1) { + fprintf(stderr, "curco--v2 intersects vb--vc\n"); + } +#endif + t = tout->next; + state_through_vert = false; } + else if (ccw1 == 1 && tout->face->visit_index != visit) { + /* vcur--v2 should intersect va--vc. */ +#ifdef DEBUG_CDT + if (dbg_level > 1) { + fprintf(stderr, "curco--v2 intersects va--vc\n"); + } +#endif + t = tout; + state_through_vert = false; + } + else if (ccw1 == 0 && vc->visit_index != visit) { +#ifdef DEBUG_CDT + if (dbg_level > 1) { + fprintf(stderr, "ccw==0 case, so going through or to vc\n"); + } #endif - t = tout->next; - state_through_vert = true; + t = tout->next; + state_through_vert = true; + } + else { +#ifdef DEBUG_CDT + fprintf(stderr, "add_edge_constraint desperation search 2\n"); +#endif + if (tout->face->visit_index != visit) { + /* Treat as if an intersection of va-vc. */ + t = tout; + state_through_vert = false; + } + } + } + if (++search_count > 10000) { + fprintf(stderr, "infinite loop? bailing out\n"); + BLI_assert(0); /* Catch these while developing. */ + break; } - } - if (++search_count > 10000) { - fprintf(stderr, "infinite loop? bailing out\n"); - BLI_assert(0); /* Catch these while developing. */ - break; } } #ifdef DEBUG_CDT @@ -1968,7 +2084,7 @@ static void dissolve_symedge(CDT_state *cdt, SymEdge *se) if (se->face->symedge == se) { se->face->symedge = se->next; } - if (symse->face->symedge == se) { + if (symse->face->symedge == symse) { symse->face->symedge = symse->next; } } @@ -2380,7 +2496,13 @@ CDT_result *BLI_delaunay_2d_cdt_calc(const CDT_input *input, const CDT_output_ty CDTEdge *face_edge; SymEdge *face_symedge; #ifdef DEBUG_CDT - int dbg_level = 1; + int dbg_level = 0; +#endif + +#ifdef DEBUG_CDT + if (dbg_level == -1) { + write_cdt_input_to_file(input); + } #endif if ((nv > 0 && input->vert_coords == NULL) || (ne > 0 && input->edges == NULL) || @@ -2446,6 +2568,15 @@ CDT_result *BLI_delaunay_2d_cdt_calc(const CDT_input *input, const CDT_output_ty continue; } add_edge_constraint(cdt, verts[v1], verts[v2], i, NULL); +#ifdef DEBUG_CDT + if (dbg_level > 3) { + char namebuf[60]; + sprintf(namebuf, "after edge constraint %d = (%d,%d)\n", i, v1, v2); + cdt_draw(cdt, namebuf); + dump_cdt(cdt, namebuf); + validate_cdt(cdt, true); + } +#endif } #ifdef DEBUG_CDT if (dbg_level > 2) { @@ -2579,20 +2710,37 @@ void BLI_delaunay_2d_cdt_free(CDT_result *result) #ifdef DEBUG_CDT -static void dump_se(const SymEdge *se, const char *lab) +static const char *vertname(const CDTVert *v) { - if (se->next) { - fprintf( - stderr, "%s((%.3f,%.3f)->(%.3f,%.3f))\n", lab, F2(se->vert->co), F2(se->next->vert->co)); + static char vertnamebuf[20]; + + if (v->index < 4) { + sprintf(vertnamebuf, "[%c]", "ABCD"[v->index]); } else { - fprintf(stderr, "%s((%.3f,%.3f)->NULL)\n", lab, F2(se->vert->co)); + sprintf(vertnamebuf, "[%d]", v->index - 4); } + return vertnamebuf; } static void dump_v(const CDTVert *v, const char *lab) { - fprintf(stderr, "%s(%.3f,%.3f)\n", lab, F2(v->co)); + fprintf(stderr, "%s%s(%.3f,%.3f)\n", lab, vertname(v), F2(v->co)); +} + +static void dump_se(const SymEdge *se, const char *lab) +{ + if (se->next) { + fprintf(stderr, + "%s%s((%.3f,%.3f)->(%.3f,%.3f))\n", + lab, + vertname(se->vert), + F2(se->vert->co), + F2(se->next->vert->co)); + } + else { + fprintf(stderr, "%s%s((%.3f,%.3f)->NULL)\n", lab, vertname(se->vert), F2(se->vert->co)); + } } static void dump_se_cycle(const SymEdge *se, const char *lab, const int limit) @@ -2622,21 +2770,12 @@ static void dump_id_list(const LinkNode *id_list, const char *lab) } } -static const char *vertname(CDTVert *v) -{ - static char vertnamebuf[20]; - - if (v->index < 4) { - sprintf(vertnamebuf, "[%c]", "ABCD"[v->index]); - } - else { - sprintf(vertnamebuf, "[%d]", v->index - 4); - } - return vertnamebuf; -} - +/* If filter_fn != NULL, only dump vert v its edges when filter_fn(cdt, v, filter_data) is true. */ # define PL(p) (POINTER_AS_UINT(p) & 0xFFFF) -static void dump_cdt(const CDT_state *cdt, const char *lab) +static void dump_cdt_filtered(const CDT_state *cdt, + bool (*filter_fn)(const CDT_state *, int, void *), + void *filter_data, + const char *lab) { LinkNode *ln; CDTVert *v, *vother; @@ -2648,6 +2787,9 @@ static void dump_cdt(const CDT_state *cdt, const char *lab) fprintf(stderr, "\nCDT %s\n", lab); fprintf(stderr, "\nVERTS\n"); for (i = 0; i < cdt->vert_array_len; i++) { + if (filter_fn && !filter_fn(cdt, i, filter_data)) { + continue; + } v = cdt->vert_array[i]; fprintf(stderr, "%s %x: (%f,%f) symedge=%x\n", vertname(v), PL(v), F2(v->co), PL(v->symedge)); dump_id_list(v->input_ids, " "); @@ -2664,6 +2806,9 @@ static void dump_cdt(const CDT_state *cdt, const char *lab) fprintf(stderr, "\n"); } } + if (filter_fn) { + return; + } fprintf(stderr, "\nEDGES\n"); for (ln = cdt->edges; ln; ln = ln->next) { e = (CDTEdge *)ln->link; @@ -2709,6 +2854,55 @@ static void dump_cdt(const CDT_state *cdt, const char *lab) } # undef PL +static void dump_cdt(const CDT_state *cdt, const char *lab) +{ + dump_cdt_filtered(cdt, NULL, NULL, lab); +} + +typedef struct ReachableFilterData { + int vstart_index; + int maxdist; +} ReachableFilterData; + +/* Stupid algorithm will repeat itself. Don't use for large n. */ +static bool reachable_filter(const CDT_state *cdt, int v_index, void *filter_data) +{ + CDTVert *v; + SymEdge *se; + ReachableFilterData *rfd_in = (ReachableFilterData *)filter_data; + ReachableFilterData rfd_next; + + if (v_index == rfd_in->vstart_index) { + return true; + } + if (rfd_in->maxdist <= 0 || v_index < 0 || v_index >= cdt->vert_array_len) { + return false; + } + else { + v = cdt->vert_array[v_index]; + se = v->symedge; + if (se != NULL) { + rfd_next.vstart_index = rfd_in->vstart_index; + rfd_next.maxdist = rfd_in->maxdist - 1; + do { + if (reachable_filter(cdt, se->next->vert->index, &rfd_next)) { + return true; + } + se = se->rot; + } while (se != v->symedge); + } + } + return false; +} + +static void dump_cdt_vert_neighborhood(CDT_state *cdt, int v, int maxdist, const char *lab) +{ + ReachableFilterData rfd; + rfd.vstart_index = v; + rfd.maxdist = maxdist; + dump_cdt_filtered(cdt, reachable_filter, &rfd, lab); +} + /** * Make an html file with svg in it to display the argument cdt. * Mouse-overs will reveal the coordinates of vertices and edges. @@ -2716,8 +2910,12 @@ static void dump_cdt(const CDT_state *cdt, const char *lab) * The first call creates DRAWFILE; subsequent calls append to it. */ # define DRAWFILE "/tmp/debug_draw.html" -# define MAX_DRAW_WIDTH 1000 -# define MAX_DRAW_HEIGHT 700 +# define MAX_DRAW_WIDTH 2000 +# define MAX_DRAW_HEIGHT 1400 +# define THIN_LINE 1 +# define THICK_LINE 3 +# define VERT_RADIUS 3 +# define DRAW_VERT_LABELS 1 static void cdt_draw(CDT_state *cdt, const char *lab) { static bool append = false; @@ -2727,9 +2925,7 @@ static void cdt_draw(CDT_state *cdt, const char *lab) double maxx = cdt->maxx + draw_margin; double miny = cdt->miny - draw_margin; double maxy = cdt->maxy + draw_margin; - double width = maxx - minx; - double height = maxy - miny; - double aspect = height / width; + double width, height, aspect; int view_width, view_height; double scale; LinkNode *ln; @@ -2737,6 +2933,10 @@ static void cdt_draw(CDT_state *cdt, const char *lab) CDTEdge *e; int i, strokew; + /* Note: to debug a small area: assign custom min's/max's here. */ + width = maxx - minx; + height = maxy - miny; + aspect = height / width; view_width = MAX_DRAW_WIDTH; view_height = (int)(view_width * aspect); if (view_height > MAX_DRAW_HEIGHT) { @@ -2767,28 +2967,36 @@ static void cdt_draw(CDT_state *cdt, const char *lab) } u = e->symedges[0].vert; v = e->symedges[1].vert; - strokew = is_constrained_edge(e) ? 5 : 2; + strokew = is_constrained_edge(e) ? THICK_LINE : THIN_LINE; fprintf(f, "\n", + "x1=\"%f\" y1=\"%f\" x2=\"%f\" y2=\"%f\">\n", strokew, SX(u->co[0]), SY(u->co[1]), SX(v->co[0]), SY(v->co[1])); - fprintf( - f, " (%.3f,%.3f)(%.3f,%.3f)\n", u->co[0], u->co[1], v->co[0], v->co[1]); + fprintf(f, " %s", vertname(u)); + fprintf(f, "%s\n", vertname(v)); fprintf(f, "\n"); } i = cdt->output_prepared ? NUM_BOUND_VERTS : 0; for (; i < cdt->vert_array_len; i++) { v = cdt->vert_array[i]; fprintf(f, - "\n", + "\n", SX(v->co[0]), - SY(v->co[1])); - fprintf(f, " (%.3f,%.3f)\n", v->co[0], v->co[1]); + SY(v->co[1]), + VERT_RADIUS); + fprintf(f, " %s(%.3f,%.3f)\n", vertname(v), v->co[0], v->co[1]); fprintf(f, "\n"); +# if DRAW_VERT_LABELS + fprintf(f, + "%s\n", + SX(v->co[0]) + VERT_RADIUS, + SY(v->co[1]) - VERT_RADIUS, + vertname(v)); +# endif } fprintf(f, "\n\n"); @@ -2798,6 +3006,28 @@ static void cdt_draw(CDT_state *cdt, const char *lab) # undef SY } +# define CDTFILE "/tmp/cdtinput.txt" +static void write_cdt_input_to_file(const CDT_input *inp) +{ + int i, j; + FILE *f = fopen(CDTFILE, "w"); + + fprintf(f, "%d %d %d\n", inp->verts_len, inp->edges_len, inp->faces_len); + for (i = 0; i < inp->verts_len; i++) { + fprintf(f, "%.17f %.17f\n", inp->vert_coords[i][0], inp->vert_coords[i][1]); + } + for (i = 0; i < inp->edges_len; i++) { + fprintf(f, "%d %d\n", inp->edges[i][0], inp->edges[i][1]); + } + for (i = 0; i < inp->faces_len; i++) { + for (j = 0; j < inp->faces_len_table[i]; j++) { + fprintf(f, "%d ", inp->faces[j + inp->faces_start_table[i]]); + } + fprintf(f, "\n"); + } + fclose(f); +} + # ifndef NDEBUG /* Only used in assert. */ /** * Is a visible from b: i.e., ab crosses no edge of cdt? diff --git a/tests/gtests/blenlib/BLI_delaunay_2d_test.cc b/tests/gtests/blenlib/BLI_delaunay_2d_test.cc index 8b29128eb0d..85e8ae3f141 100644 --- a/tests/gtests/blenlib/BLI_delaunay_2d_test.cc +++ b/tests/gtests/blenlib/BLI_delaunay_2d_test.cc @@ -15,6 +15,10 @@ extern "C" { #include #include +#define DO_REGULAR_TESTS 1 +#define DO_RANDOM_TESTS 0 +#define DO_FILE_TESTS 0 + static void fill_input_verts(CDT_input *r_input, float (*vcos)[2], int nverts) { r_input->verts_len = nverts; @@ -43,6 +47,117 @@ static void add_input_faces( r_input->faces_len_table = faces_len_table; } +/* The spec should have the form: + * #verts #edges #faces + * [#verts lines) + * [#edges lines] + * ... [#faces lines] + */ +static void fill_input_from_string(CDT_input *r_input, const char *spec) +{ + std::string line; + std::vector> faces; + int i, j; + int nverts, nedges, nfaces; + float(*p)[2]; + int(*e)[2]; + int *farr; + int *flen; + int *fstart; + + std::istringstream ss(spec); + getline(ss, line); + std::istringstream hdrss(line); + hdrss >> nverts >> nedges >> nfaces; + if (nverts == 0) { + return; + } + p = (float(*)[2])MEM_malloc_arrayN(nverts, 2 * sizeof(float), __func__); + if (nedges > 0) { + e = (int(*)[2])MEM_malloc_arrayN(nedges, 2 * sizeof(int), __func__); + } + if (nfaces > 0) { + flen = (int *)MEM_malloc_arrayN(nfaces, sizeof(int), __func__); + fstart = (int *)MEM_malloc_arrayN(nfaces, sizeof(int), __func__); + } + i = 0; + while (i < nverts && getline(ss, line)) { + std::istringstream iss(line); + iss >> p[i][0] >> p[i][1]; + i++; + } + i = 0; + while (i < nedges && getline(ss, line)) { + std::istringstream ess(line); + ess >> e[i][0] >> e[i][1]; + i++; + } + i = 0; + while (i < nfaces && getline(ss, line)) { + std::istringstream fss(line); + int v; + faces.push_back(std::vector()); + while (fss >> v) { + faces[i].push_back(v); + } + i++; + } + fill_input_verts(r_input, p, nverts); + if (nedges > 0) { + add_input_edges(r_input, e, nedges); + } + if (nfaces > 0) { + for (i = 0; i < nfaces; i++) { + flen[i] = (int)faces[i].size(); + if (i == 0) { + fstart[i] = 0; + } + else { + fstart[i] = fstart[i - 1] + flen[i - 1]; + } + } + farr = (int *)MEM_malloc_arrayN(fstart[nfaces - 1] + flen[nfaces - 1], sizeof(int), __func__); + for (i = 0; i < nfaces; i++) { + for (j = 0; j < (int)faces[i].size(); j++) { + farr[fstart[i] + j] = faces[i][j]; + } + } + add_input_faces(r_input, farr, fstart, flen, nfaces); + } +} + +static void fill_input_from_file(CDT_input *in, const char *filename) +{ + std::FILE *fp = std::fopen(filename, "rb"); + if (fp) { + std::string contents; + std::fseek(fp, 0, SEEK_END); + contents.resize(std::ftell(fp)); + std::rewind(fp); + std::fread(&contents[0], 1, contents.size(), fp); + std::fclose(fp); + fill_input_from_string(in, contents.c_str()); + } + else { + printf("couldn't open file %s\n", filename); + } +} + +static void free_spec_arrays(CDT_input *in) +{ + if (in->vert_coords) { + MEM_freeN(in->vert_coords); + } + if (in->edges) { + MEM_freeN(in->edges); + } + if (in->faces_len_table) { + MEM_freeN(in->faces_len_table); + MEM_freeN(in->faces_start_table); + MEM_freeN(in->faces); + } +} + /* which output vert index goes with given input vertex? -1 if not found */ static int get_output_vert_index(const CDT_result *r, int in_index) { @@ -187,6 +302,7 @@ static void dump_result(CDT_result *r) } } +#if DO_REGULAR_TESTS TEST(delaunay, Empty) { CDT_input in; @@ -692,7 +808,108 @@ TEST(delaunay, TriInTri) BLI_delaunay_2d_cdt_free(out); } -#if 0 +TEST(delaunay, DiamondInSquare) +{ + CDT_input in; + CDT_result *out; + const char *spec = R"(8 0 2 + 0.0 0.0 + 1.0 0.0 + 1.0 1.0 + 0.0 1.0 + 0.14644660940672627 0.5 + 0.5 0.14644660940672627 + 0.8535533905932737 0.5 + 0.5 0.8535533905932737 + 0 1 2 3 + 4 5 6 7 + )"; + fill_input_from_string(&in, spec); + out = BLI_delaunay_2d_cdt_calc(&in, CDT_CONSTRAINTS_VALID_BMESH); + EXPECT_EQ(out->verts_len, 8); + EXPECT_EQ(out->edges_len, 10); + EXPECT_EQ(out->faces_len, 3); + free_spec_arrays(&in); +} + +TEST(delaunay, DiamondInSquareWire) +{ + CDT_input in; + CDT_result *out; + const char *spec = R"(8 8 0 + 0.0 0.0 + 1.0 0.0 + 1.0 1.0 + 0.0 1.0 + 0.14644660940672627 0.5 + 0.5 0.14644660940672627 + 0.8535533905932737 0.5 + 0.5 0.8535533905932737 + 0 1 + 1 2 + 2 3 + 3 0 + 4 5 + 5 6 + 6 7 + 7 4 + )"; + fill_input_from_string(&in, spec); + out = BLI_delaunay_2d_cdt_calc(&in, CDT_CONSTRAINTS); + EXPECT_EQ(out->verts_len, 8); + EXPECT_EQ(out->edges_len, 8); + EXPECT_EQ(out->faces_len, 2); + free_spec_arrays(&in); +} + +TEST(delaunay, ClosePts) +{ + CDT_input in; + CDT_result *out; + const char *spec = R"(7 2 1 + 0.46876350045204163 0.06087132915854454 + 0.46865847706794739 0.03632887825369835 + 0.49176687002182007 0.03632888197898865 + 0.49166208505630493 0.06087132543325424 + 0.49171400070190430 0.04841339960694313 + 0.49171534180641174 0.04839951172471046 + 0.49045535922050476 0.06087132915854454 + 4 5 + 6 4 + 0 1 2 3 + )"; + fill_input_from_string(&in, spec); + out = BLI_delaunay_2d_cdt_calc(&in, CDT_FULL); + EXPECT_EQ(out->verts_len, 7); + EXPECT_EQ(out->edges_len, 12); + EXPECT_EQ(out->faces_len, 6); + free_spec_arrays(&in); +} + +TEST(delaunay, ClosePts2) +{ + CDT_input in; + CDT_result *out; + const char *spec = R"(6 1 1 + -0.17878936231136322 -0.44374340772628784 + -0.17871695756912231 -0.45601493120193481 + -0.17544384300708771 -0.45601493120193481 + -0.17537136375904083 -0.44374340772628784 + -0.17544738948345184 -0.45602506399154663 + -0.17872454226016998 -0.45472940802574158 + 4 5 + 0 1 2 3 + )"; + fill_input_from_string(&in, spec); + out = BLI_delaunay_2d_cdt_calc(&in, CDT_FULL); + EXPECT_EQ(out->verts_len, 6); + EXPECT_EQ(out->edges_len, 10); + EXPECT_EQ(out->faces_len, 5); + free_spec_arrays(&in); +} +#endif + +#if DO_RANDOM_TESTS enum { RANDOM_PTS, RANDOM_SEGS, @@ -780,17 +997,17 @@ static void rand_delaunay_test(int test_kind, TEST(delaunay, randompts) { - rand_delaunay_test(RANDOM_PTS, 7, 1, CDT_FULL); + rand_delaunay_test(RANDOM_PTS, 7, 100, CDT_FULL); } TEST(delaunay, randomsegs) { - rand_delaunay_test(RANDOM_SEGS, 7, 1, CDT_FULL); + rand_delaunay_test(RANDOM_SEGS, 7, 100, CDT_FULL); } TEST(delaunay, randompoly) { - rand_delaunay_test(RANDOM_POLY, 7, 1, CDT_FULL); + rand_delaunay_test(RANDOM_POLY, 7, 100, CDT_FULL); } TEST(delaunay, randompoly_inside) @@ -809,48 +1026,36 @@ TEST(delaunay, randompoly_validbmesh) } #endif -#if 0 -/* For debugging or timing large examples. - * The given file should have one point per line, as a space-separated pair of floats +#if DO_FILE_TESTS +/* For timing large examples of points only. + * The given file should have one point per line, as a space-separated pair of floats. */ static void points_from_file_test(const char *filename) { - std::ifstream f; - std::string line; - struct XY { - float x; - float y; - } xy; - std::vector pts; - int i, n; CDT_input in; CDT_result *out; double tstart; - float (*p)[2]; - f.open(filename); - while (getline(f, line)) { - std::istringstream iss(line); - iss >> xy.x >> xy.y; - pts.push_back(xy); - } - n = (int)pts.size(); - fprintf(stderr, "read %d points\n", (int)pts.size()); - p = (float (*)[2])MEM_malloc_arrayN(n, 2 * sizeof(float), "delaunay"); - for (i = 0; i < n; i++) { - xy = pts[i]; - p[i][0] = xy.x; - p[i][1] = xy.y; - } + fill_input_from_file(&in, filename); tstart = PIL_check_seconds_timer(); - fill_input_verts(&in, p, n); out = BLI_delaunay_2d_cdt_calc(&in, CDT_FULL); fprintf(stderr, "time to triangulate=%f seconds\n", PIL_check_seconds_timer() - tstart); BLI_delaunay_2d_cdt_free(out); - MEM_freeN(p); + free_spec_arrays(&in); +} + +TEST(delaunay, debug) +{ + CDT_input in; + CDT_result *out; + fill_input_from_file(&in, "/tmp/cdtinput.txt"); + out = BLI_delaunay_2d_cdt_calc(&in, CDT_FULL); + BLI_delaunay_2d_cdt_free(out); + free_spec_arrays(&in); } -# define POINTFILEROOT "/tmp/" +# if 0 +# define POINTFILEROOT "/tmp/" TEST(delaunay, terrain1) { @@ -866,4 +1071,5 @@ TEST(delaunay, terrain3) { points_from_file_test(POINTFILEROOT "points3.txt"); } +# endif #endif -- cgit v1.2.3