diff options
author | Howard Trickey <howard.trickey@gmail.com> | 2019-12-21 20:23:02 +0300 |
---|---|---|
committer | Howard Trickey <howard.trickey@gmail.com> | 2019-12-21 20:23:02 +0300 |
commit | a0892bb690fec36edd7999d6729e0826d72dab86 (patch) | |
tree | da1d94d5e26a026dabeabaaa67fb2769e1c11d53 /source/blender/blenlib | |
parent | 51d8d790d75f61cc30af08f03db92c38b196d793 (diff) |
Fix crash in delaunay triangulation due to epsilon issues.
Diffstat (limited to 'source/blender/blenlib')
-rw-r--r-- | source/blender/blenlib/intern/delaunay_2d.c | 586 |
1 files changed, 408 insertions, 178 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, "<line fill=\"none\" stroke=\"black\" stroke-width=\"%d\" " - "x1=\"%.3f\" y1=\"%.3f\" x2=\"%.3f\" y2=\"%.3f\">\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, " <title>(%.3f,%.3f)(%.3f,%.3f)</title>\n", u->co[0], u->co[1], v->co[0], v->co[1]); + fprintf(f, " <title>%s", vertname(u)); + fprintf(f, "%s</title>\n", vertname(v)); fprintf(f, "</line>\n"); } i = cdt->output_prepared ? NUM_BOUND_VERTS : 0; for (; i < cdt->vert_array_len; i++) { v = cdt->vert_array[i]; fprintf(f, - "<circle fill=\"black\" cx=\"%.3f\" cy=\"%.3f\" r=\"5\">\n", + "<circle fill=\"black\" cx=\"%f\" cy=\"%f\" r=\"%d\">\n", SX(v->co[0]), - SY(v->co[1])); - fprintf(f, " <title>(%.3f,%.3f)</title>\n", v->co[0], v->co[1]); + SY(v->co[1]), + VERT_RADIUS); + fprintf(f, " <title>%s(%.3f,%.3f)</title>\n", vertname(v), v->co[0], v->co[1]); fprintf(f, "</circle>\n"); +# if DRAW_VERT_LABELS + fprintf(f, + "<text x=\"%f\" y=\"%f\" font-size=\"small\">%s</text>\n", + SX(v->co[0]) + VERT_RADIUS, + SY(v->co[1]) - VERT_RADIUS, + vertname(v)); +# endif } fprintf(f, "</svg>\n</div>\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? |