/* * This program is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License * as published by the Free Software Foundation; either version 2 * of the License, or (at your option) any later version. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with this program; if not, write to the Free Software Foundation, * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. */ /** \file * \ingroup bli * * Constrained 2d Delaunay Triangulation. */ #include "MEM_guardedalloc.h" #include "BLI_array.h" #include "BLI_bitmap.h" #include "BLI_linklist.h" #include "BLI_math.h" #include "BLI_memarena.h" #include "BLI_mempool.h" #include "BLI_delaunay_2d.h" /* Uncomment this define to get helpful debugging functions etc. defined. */ #define DEBUG_CDT struct CDTEdge; struct CDTFace; struct CDTVert; typedef struct SymEdge { struct SymEdge *next; /* In face, doing CCW traversal of face. */ struct SymEdge *rot; /* CCW around vert. */ struct CDTVert *vert; /* Vert at origin. */ struct CDTEdge *edge; /* Undirected edge this is for. */ struct CDTFace *face; /* Face on left side. */ } SymEdge; typedef struct CDTVert { double co[2]; /* Coordinate. */ 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 merge_to_index; /* Index of a CDTVert that this has merged to. -1 if no merge. */ } CDTVert; typedef struct CDTEdge { LinkNode *input_ids; /* List of input edge ids that this is part of. */ SymEdge symedges[2]; /* The directed edges for this edge. */ bool in_queue; /* Used in flipping algorithm. */ } CDTEdge; typedef struct CDTFace { SymEdge *symedge; /* A symedge in face; only used during output. */ LinkNode *input_ids; /* List of input face ids that this is part of. */ int visit_index; /* Which visit epoch has this been seen. */ bool deleted; /* Marks this face no longer used. */ bool in_queue; /* Used in remove_small_features algorithm. */ } CDTFace; typedef struct CDT_state { LinkNode *edges; LinkNode *faces; CDTFace *outer_face; CDTVert **vert_array; int vert_array_len; int vert_array_len_alloc; int input_vert_tot; double minx; double miny; double maxx; double maxy; double margin; int visit_count; int face_edge_offset; MemArena *arena; BLI_mempool *listpool; double epsilon; double epsilon_squared; bool output_prepared; } CDT_state; #define DLNY_ARENASIZE 1 << 14 /** * This margin means that will only be a 1 degree possible * concavity on outside if remove all border touching triangles. */ #define DLNY_MARGIN_PCT 2000.0 #ifdef DEBUG_CDT # ifdef __GNUC__ # define ATTU __attribute__((unused)) # else # define ATTU # endif # define F2(p) p[0], p[1] # define F3(p) p[0], p[1], p[2] ATTU static void dump_se(const SymEdge *se, const char *lab); ATTU static void dump_v(const CDTVert *v, const char *lab); ATTU static void dump_se_cycle(const SymEdge *se, const char *lab, const int limit); ATTU static void dump_id_list(const LinkNode *id_list, const char *lab); ATTU static void dump_cdt(const CDT_state *cdt, const char *lab); ATTU static void dump_cdt_vert_neighborhood(CDT_state *cdt, int v, int maxdist, const char *lab); ATTU static void cdt_draw(CDT_state *cdt, const char *lab); ATTU static void cdt_draw_vertex_region(CDT_state *cdt, int v, double dist, const char *lab); ATTU static void write_cdt_input_to_file(const CDT_input *inp); ATTU static void validate_cdt(CDT_state *cdt, bool check_all_tris, bool check_delaunay, bool check_visibility); #endif static void exactinit(void); static double orient2d(const double *pa, const double *pb, const double *pc); static double incircle(const double *pa, const double *pb, const double *pc, const double *pd); /** Return other #SymEdge for same #CDTEdge as se. */ BLI_INLINE SymEdge *sym(const SymEdge *se) { return se->next->rot; } /** Return SymEdge whose next is se. */ BLI_INLINE SymEdge *prev(const SymEdge *se) { return se->rot->next->rot; } /** Return true if a -- b -- c are in that order, assuming they are on a straight line according to * orient2d and we know the order is either abc or bac. * This means ab . ac and bc . ac must both be non-negative. */ static bool in_line(const double a[2], const double b[2], const double c[2]) { double ab[2], bc[2], ac[2]; sub_v2_v2v2_db(ab, b, a); sub_v2_v2v2_db(bc, c, b); sub_v2_v2v2_db(ac, c, a); if (dot_v2v2_db(ab, ac) < 0.0) { return false; } return dot_v2v2_db(bc, ac) >= 0.0; } #ifndef NDEBUG /** Is s2 reachable from s1 by next pointers with < limit hops? */ static bool reachable(SymEdge *s1, SymEdge *s2, int limit) { int count = 0; for (SymEdge *s = s1; s && count < limit; s = s->next) { if (s == s2) { return true; } count++; } return false; } #endif /** Using array to store these instead of linked list so can make a random selection from them. */ static CDTVert *add_cdtvert(CDT_state *cdt, double x, double y) { CDTVert *v = BLI_memarena_alloc(cdt->arena, sizeof(*v)); v->co[0] = x; v->co[1] = y; v->input_ids = NULL; v->symedge = NULL; if (cdt->vert_array_len == cdt->vert_array_len_alloc) { CDTVert **old_array = cdt->vert_array; cdt->vert_array_len_alloc *= 4; cdt->vert_array = BLI_memarena_alloc(cdt->arena, cdt->vert_array_len_alloc * sizeof(cdt->vert_array[0])); memmove(cdt->vert_array, old_array, cdt->vert_array_len * sizeof(cdt->vert_array[0])); } BLI_assert(cdt->vert_array_len < cdt->vert_array_len_alloc); v->index = cdt->vert_array_len; v->merge_to_index = -1; cdt->vert_array[cdt->vert_array_len++] = v; return v; } static CDTEdge *add_cdtedge( CDT_state *cdt, CDTVert *v1, CDTVert *v2, CDTFace *fleft, CDTFace *fright) { CDTEdge *e = BLI_memarena_alloc(cdt->arena, sizeof(*e)); SymEdge *se = &e->symedges[0]; SymEdge *sesym = &e->symedges[1]; e->input_ids = NULL; e->in_queue = false; BLI_linklist_prepend_arena(&cdt->edges, (void *)e, cdt->arena); se->edge = sesym->edge = e; se->face = fleft; sesym->face = fright; se->vert = v1; if (v1->symedge == NULL) { v1->symedge = se; } sesym->vert = v2; if (v2->symedge == NULL) { v2->symedge = sesym; } se->next = sesym->next = se->rot = sesym->rot = NULL; return e; } static CDTFace *add_cdtface(CDT_state *cdt) { CDTFace *f = BLI_memarena_alloc(cdt->arena, sizeof(*f)); f->visit_index = 0; f->deleted = false; f->symedge = NULL; f->input_ids = NULL; f->in_queue = false; BLI_linklist_prepend_arena(&cdt->faces, (void *)f, cdt->arena); return f; } static bool id_in_list(const LinkNode *id_list, int id) { const LinkNode *ln; for (ln = id_list; ln; ln = ln->next) { if (POINTER_AS_INT(ln->link) == id) { return true; } } return false; } /** is any id in (range_start, range_start+1, ... , range_end) in id_list? */ static bool id_range_in_list(const LinkNode *id_list, int range_start, int range_end) { const LinkNode *ln; int id; for (ln = id_list; ln; ln = ln->next) { id = POINTER_AS_INT(ln->link); if (id >= range_start && id <= range_end) { return true; } } return false; } static void add_to_input_ids(LinkNode **dst, int input_id, CDT_state *cdt) { if (!id_in_list(*dst, input_id)) { BLI_linklist_prepend_arena(dst, POINTER_FROM_INT(input_id), cdt->arena); } } static void add_list_to_input_ids(LinkNode **dst, const LinkNode *src, CDT_state *cdt) { const LinkNode *ln; for (ln = src; ln; ln = ln->next) { add_to_input_ids(dst, POINTER_AS_INT(ln->link), cdt); } } BLI_INLINE bool is_border_edge(const CDTEdge *e, const CDT_state *cdt) { return e->symedges[0].face == cdt->outer_face || e->symedges[1].face == cdt->outer_face; } BLI_INLINE bool is_constrained_edge(const CDTEdge *e) { return e->input_ids != NULL; } BLI_INLINE bool is_deleted_edge(const CDTEdge *e) { return e->symedges[0].next == NULL; } BLI_INLINE bool is_original_vert(const CDTVert *v, CDT_state *cdt) { return (v->index < cdt->input_vert_tot); } /** Is there already an edge between a and b? */ static bool exists_edge(const CDTVert *a, const CDTVert *b) { SymEdge *se, *ss; se = a->symedge; if (se->next->vert == b) { return true; } for (ss = se->rot; ss != se; ss = ss->rot) { if (ss->next->vert == b) { return true; } } return false; } /** Is the vertex v incident on face f? */ static bool vert_touches_face(const CDTVert *v, const CDTFace *f) { SymEdge *se = v->symedge; do { if (se->face == f) { return true; } } while ((se = se->rot) != v->symedge); return false; } /** * Assume s1 and s2 are both SymEdges in a face with > 3 sides, * and one is not the next of the other. * Add an edge from s1->v to s2->v, splitting the face in two. * The original face will continue to be associated with the subface * that has s1, and a new face will be made for s2's new face. * Return the new diagonal's CDTEdge *. */ static CDTEdge *add_diagonal(CDT_state *cdt, SymEdge *s1, SymEdge *s2) { CDTEdge *ediag; CDTFace *fold, *fnew; SymEdge *sdiag, *sdiagsym; SymEdge *s1prev, *s1prevsym, *s2prev, *s2prevsym, *se; BLI_assert(reachable(s1, s2, 20000)); BLI_assert(reachable(s2, s1, 20000)); fold = s1->face; fnew = add_cdtface(cdt); s1prev = prev(s1); s1prevsym = sym(s1prev); s2prev = prev(s2); s2prevsym = sym(s2prev); ediag = add_cdtedge(cdt, s1->vert, s2->vert, fnew, fold); sdiag = &ediag->symedges[0]; sdiagsym = &ediag->symedges[1]; sdiag->next = s2; sdiagsym->next = s1; s2prev->next = sdiagsym; s1prev->next = sdiag; s1->rot = sdiag; sdiag->rot = s1prevsym; s2->rot = sdiagsym; sdiagsym->rot = s2prevsym; #ifdef DEBUG_CDT BLI_assert(reachable(s2, sdiag, 2000)); #endif for (se = s2; se != sdiag; se = se->next) { se->face = fnew; } add_list_to_input_ids(&fnew->input_ids, fold->input_ids, cdt); return ediag; } /** * Add a dangling edge from an isolated v to the vert at se in the same face as se->face. */ static CDTEdge *add_vert_to_symedge_edge(CDT_state *cdt, CDTVert *v, SymEdge *se) { CDTEdge *e; SymEdge *se_rot, *se_rotsym, *new_se, *new_se_sym; se_rot = se->rot; se_rotsym = sym(se_rot); e = add_cdtedge(cdt, v, se->vert, se->face, se->face); new_se = &e->symedges[0]; new_se_sym = &e->symedges[1]; new_se->next = se; new_se_sym->next = new_se; new_se->rot = new_se; new_se_sym->rot = se_rot; se->rot = new_se_sym; se_rotsym->next = new_se_sym; return e; } /* Connect the verts of se1 and se2, assuming that currently those two SymEdges are on * the outer boundary (have face == outer_face) of two components that are isolated from * each other. */ static CDTEdge *connect_separate_parts(CDT_state *cdt, SymEdge *se1, SymEdge *se2) { CDTEdge *e; SymEdge *se1_rot, *se1_rotsym, *se2_rot, *se2_rotsym, *new_se, *new_se_sym; ; BLI_assert(se1->face == cdt->outer_face && se2->face == cdt->outer_face); se1_rot = se1->rot; se1_rotsym = sym(se1_rot); se2_rot = se2->rot; se2_rotsym = sym(se2_rot); e = add_cdtedge(cdt, se1->vert, se2->vert, cdt->outer_face, cdt->outer_face); new_se = &e->symedges[0]; new_se_sym = &e->symedges[1]; new_se->next = se2; new_se_sym->next = se1; new_se->rot = se1_rot; new_se_sym->rot = se2_rot; se1->rot = new_se; se2->rot = new_se_sym; se1_rotsym->next = new_se; se2_rotsym->next = new_se_sym; return e; } /** * Split \a se at fraction \a lambda, * and return the new #CDTEdge that is the new second half. * Copy the edge input_ids into the new one. */ static CDTEdge *split_edge(CDT_state *cdt, SymEdge *se, double lambda) { const double *a, *b; double p[2]; CDTVert *v; CDTEdge *e; SymEdge *sesym, *newse, *newsesym, *senext, *sesymprev, *sesymprevsym; /* Split e at lambda. */ a = se->vert->co; b = se->next->vert->co; sesym = sym(se); sesymprev = prev(sesym); sesymprevsym = sym(sesymprev); senext = se->next; p[0] = (1.0 - lambda) * a[0] + lambda * b[0]; p[1] = (1.0 - lambda) * a[1] + lambda * b[1]; v = add_cdtvert(cdt, p[0], p[1]); e = add_cdtedge(cdt, v, se->next->vert, se->face, sesym->face); sesym->vert = v; newse = &e->symedges[0]; newsesym = &e->symedges[1]; se->next = newse; newsesym->next = sesym; newse->next = senext; newse->rot = sesym; sesym->rot = newse; senext->rot = newsesym; newsesym->rot = sesymprevsym; sesymprev->next = newsesym; if (newsesym->vert->symedge == sesym) { newsesym->vert->symedge = newsesym; } add_list_to_input_ids(&e->input_ids, se->edge->input_ids, cdt); return e; } /** * Delete an edge from the structure. The new combined face on either side of * the deleted edge will be the one that was e's face. * There will be now an unused face, marked by setting its deleted flag, * and an unused #CDTEdge, marked by setting the next and rot pointers of * its SymEdges to NULL. *
* . v2 . * / \ / \ * /f|j\ / \ * / | \ / \ * | * A | B A * \ e| / \ / * \ | / \ / * \h|i/ \ / * . v1 . ** Also handle variant cases where one or both ends * are attached only to e. */ static void delete_edge(CDT_state *cdt, SymEdge *e) { SymEdge *esym, *f, *h, *i, *j, *k, *jsym, *hsym; CDTFace *aface, *bface; CDTVert *v1, *v2; bool v1_isolated, v2_isolated; esym = sym(e); v1 = e->vert; v2 = esym->vert; aface = e->face; bface = esym->face; f = e->next; h = prev(e); i = esym->next; j = prev(esym); jsym = sym(j); hsym = sym(h); v1_isolated = (i == e); v2_isolated = (f == esym); if (!v1_isolated) { h->next = i; i->rot = hsym; } if (!v2_isolated) { j->next = f; f->rot = jsym; } if (!v1_isolated && !v2_isolated && aface != bface) { for (k = i; k != f; k = k->next) { k->face = aface; } } /* If e was representative symedge for v1 or v2, fix that. */ if (v1_isolated) { v1->symedge = NULL; } else if (v1->symedge == e) { v1->symedge = i; } if (v2_isolated) { v2->symedge = NULL; } else if (v2->symedge == esym) { v2->symedge = f; } /* Mark SymEdge as deleted by setting all its pointers to NULL. */ e->next = e->rot = NULL; esym->next = esym->rot = NULL; if (!v1_isolated && !v2_isolated && aface != bface) { bface->deleted = true; if (cdt->outer_face == bface) { cdt->outer_face = aface; } } } static CDT_state *new_cdt_init(const CDT_input *in) { int i; MemArena *arena = BLI_memarena_new(DLNY_ARENASIZE, __func__); CDT_state *cdt = BLI_memarena_calloc(arena, sizeof(CDT_state)); cdt->epsilon = (double)in->epsilon; cdt->epsilon_squared = cdt->epsilon * cdt->epsilon; cdt->arena = arena; cdt->input_vert_tot = in->verts_len; cdt->vert_array_len_alloc = 2 * in->verts_len; cdt->vert_array = BLI_memarena_alloc(arena, cdt->vert_array_len_alloc * sizeof(*cdt->vert_array)); cdt->listpool = BLI_mempool_create( sizeof(LinkNode), 128 + 4 * in->verts_len, 128 + in->verts_len, 0); for (i = 0; i < in->verts_len; i++) { add_cdtvert(cdt, (double)(in->vert_coords[i][0]), (double)(in->vert_coords[i][1])); } cdt->outer_face = add_cdtface(cdt); return cdt; } static void new_cdt_free(CDT_state *cdt) { BLI_mempool_destroy(cdt->listpool); BLI_memarena_free(cdt->arena); } typedef struct SiteInfo { CDTVert *v; int orig_index; } SiteInfo; static int site_lexicographic_cmp(const void *a, const void *b) { const SiteInfo *s1 = a; const SiteInfo *s2 = b; const double *co1 = s1->v->co; const double *co2 = s2->v->co; if (co1[0] < co2[0]) { return -1; } else if (co1[0] > co2[0]) { return 1; } else if (co1[1] < co2[1]) { return -1; } else if (co1[1] > co2[1]) { return 1; } else if (s1->orig_index < s2->orig_index) { return -1; } else if (s1->orig_index > s2->orig_index) { return 1; } return 0; } BLI_INLINE bool vert_left_of_symedge(CDTVert *v, SymEdge *se) { return orient2d(v->co, se->vert->co, se->next->vert->co) > 0.0; } BLI_INLINE bool vert_right_of_symedge(CDTVert *v, SymEdge *se) { return orient2d(v->co, se->next->vert->co, se->vert->co) > 0.0; } /* Is se above basel? */ BLI_INLINE bool dc_tri_valid(SymEdge *se, SymEdge *basel, SymEdge *basel_sym) { return orient2d(se->next->vert->co, basel_sym->vert->co, basel->vert->co) > 0.0; } /* Delaunay triangulate sites[start} to sites[end-1]. * Assume sites are lexicographically sorted by coordinate. * Return SymEdge of ccw convex hull at left-most point in *r_le * and that of right-most point of cw convex null in *r_re. */ static void dc_tri( CDT_state *cdt, SiteInfo *sites, int start, int end, SymEdge **r_le, SymEdge **r_re) { int n = end - start; int n2; CDTVert *v1, *v2, *v3; CDTEdge *ea, *eb, *ebasel; SymEdge *ldo, *ldi, *rdi, *rdo, *basel, *basel_sym, *lcand, *rcand, *t; double orient; bool valid_lcand, valid_rcand; #ifdef DEBUG_CDT char label_buf[100]; int dbg_level = 0; if (dbg_level > 0) { fprintf(stderr, "DC_TRI start=%d end=%d\n", start, end); } #endif BLI_assert(r_le != NULL && r_re != NULL); if (n <= 1) { *r_le = NULL; *r_re = NULL; return; } if (n <= 3) { v1 = sites[start].v; v2 = sites[start + 1].v; ea = add_cdtedge(cdt, v1, v2, cdt->outer_face, cdt->outer_face); ea->symedges[0].next = &ea->symedges[1]; ea->symedges[1].next = &ea->symedges[0]; ea->symedges[0].rot = &ea->symedges[0]; ea->symedges[1].rot = &ea->symedges[1]; if (n == 2) { *r_le = &ea->symedges[0]; *r_re = &ea->symedges[1]; return; } v3 = sites[start + 2].v; eb = add_vert_to_symedge_edge(cdt, v3, &ea->symedges[1]); orient = orient2d(v1->co, v2->co, v3->co); if (orient > 0.0) { add_diagonal(cdt, &eb->symedges[0], &ea->symedges[0]); *r_le = &ea->symedges[0]; *r_re = &eb->symedges[0]; } else if (orient < 0.0) { add_diagonal(cdt, &ea->symedges[0], &eb->symedges[0]); *r_le = ea->symedges[0].rot; *r_re = eb->symedges[0].rot; } else { /* Collinear points. Just return a line. */ *r_le = &ea->symedges[0]; *r_re = &eb->symedges[0]; } return; } /* Here: n >= 4. Divide and conquer. */ n2 = n / 2; BLI_assert(n2 >= 2 && end - (start + n2) >= 2); /* Delaunay triangulate two halves, L and R. */ dc_tri(cdt, sites, start, start + n2, &ldo, &ldi); dc_tri(cdt, sites, start + n2, end, &rdi, &rdo); #ifdef DEBUG_CDT if (dbg_level > 0) { fprintf(stderr, "\nDC_TRI merge step for start=%d, end=%d\n", start, end); dump_se(ldo, "ldo"); dump_se(ldi, "ldi"); dump_se(rdi, "rdi"); dump_se(rdo, "rdo"); if (dbg_level > 1) { sprintf(label_buf, "dc_tri(%d,%d)(%d,%d)", start, start + n2, start + n2, end); /* dump_cdt(cdt, label_buf); */ cdt_draw(cdt, label_buf); } } #endif /* Find lower common tangent of L and R. */ for (;;) { if (vert_left_of_symedge(rdi->vert, ldi)) { ldi = ldi->next; } else if (vert_right_of_symedge(ldi->vert, rdi)) { rdi = sym(rdi)->rot; /* Previous edge to rdi with same right face. */ } else { break; } } #ifdef DEBUG_CDT if (dbg_level > 0) { fprintf(stderr, "common lower tangent is between\n"); dump_se(rdi, "rdi"); dump_se(ldi, "ldi"); } #endif ebasel = connect_separate_parts(cdt, sym(rdi)->next, ldi); basel = &ebasel->symedges[0]; basel_sym = &ebasel->symedges[1]; #ifdef DEBUG_CDT if (dbg_level > 1) { dump_se(basel, "basel"); cdt_draw(cdt, "after basel made"); } #endif if (ldi->vert == ldo->vert) { ldo = basel_sym; } if (rdi->vert == rdo->vert) { rdo = basel; } /* Merge loop. */ for (;;) { /* Locate the first point lcand->next->vert encountered by rising bubble, * and delete L edges out of basel->next->vert that fail the circle test. */ lcand = basel_sym->rot; rcand = basel_sym->next; #ifdef DEBUG_CDT if (dbg_level > 1) { fprintf(stderr, "\ntop of merge loop\n"); dump_se(lcand, "lcand"); dump_se(rcand, "rcand"); dump_se(basel, "basel"); } #endif if (dc_tri_valid(lcand, basel, basel_sym)) { #ifdef DEBUG_CDT if (dbg_level > 1) { fprintf(stderr, "found valid lcand\n"); dump_se(lcand, " lcand"); } #endif while (incircle(basel_sym->vert->co, basel->vert->co, lcand->next->vert->co, lcand->rot->next->vert->co) > 0.0) { #ifdef DEBUG_CDT if (dbg_level > 1) { fprintf(stderr, "incircle says to remove lcand\n"); dump_se(lcand, " lcand"); } #endif t = lcand->rot; delete_edge(cdt, sym(lcand)); lcand = t; } } /* Symmetrically, locate first R point to be hit and delete R edges. */ if (dc_tri_valid(rcand, basel, basel_sym)) { #ifdef DEBUG_CDT if (dbg_level > 1) { fprintf(stderr, "found valid rcand\n"); dump_se(rcand, " rcand"); } #endif while (incircle(basel_sym->vert->co, basel->vert->co, rcand->next->vert->co, sym(rcand)->next->next->vert->co) > 0.0) { #ifdef DEBUG_CDT if (dbg_level > 0) { fprintf(stderr, "incircle says to remove rcand\n"); dump_se(lcand, " rcand"); } #endif t = sym(rcand)->next; delete_edge(cdt, rcand); rcand = t; } } /* If both lcand and rcand are invalid, then basel is the common upper tangent. */ valid_lcand = dc_tri_valid(lcand, basel, basel_sym); valid_rcand = dc_tri_valid(rcand, basel, basel_sym); #ifdef DEBUG_CDT if (dbg_level > 0) { fprintf( stderr, "after bubbling up, valid_lcand=%d, valid_rcand=%d\n", valid_lcand, valid_rcand); dump_se(lcand, "lcand"); dump_se(rcand, "rcand"); } #endif if (!valid_lcand && !valid_rcand) { break; } /* The next cross edge to be connected is to either lcand->next->vert or rcand->next->vert; * if both are valid, choose the appropriate one using the incircle test. */ if (!valid_lcand || (valid_rcand && incircle(lcand->next->vert->co, lcand->vert->co, rcand->vert->co, rcand->next->vert->co) > 0.0)) { #ifdef DEBUG_CDT if (dbg_level > 0) { fprintf(stderr, "connecting rcand\n"); dump_se(basel_sym, " se1=basel_sym"); dump_se(rcand->next, " se2=rcand->next"); } #endif ebasel = add_diagonal(cdt, rcand->next, basel_sym); } else { #ifdef DEBUG_CDT if (dbg_level > 0) { fprintf(stderr, "connecting lcand\n"); dump_se(sym(lcand), " se1=sym(lcand)"); dump_se(basel_sym->next, " se2=basel_sym->next"); } #endif ebasel = add_diagonal(cdt, basel_sym->next, sym(lcand)); } basel = &ebasel->symedges[0]; basel_sym = &ebasel->symedges[1]; BLI_assert(basel_sym->face == cdt->outer_face); #ifdef DEBUG_CDT if (dbg_level > 2) { cdt_draw(cdt, "after adding new crossedge"); // dump_cdt(cdt, "after adding new crossedge"); } #endif } *r_le = ldo; *r_re = rdo; BLI_assert(sym(ldo)->face == cdt->outer_face && rdo->face == cdt->outer_face); } /* Guibas-Stolfi Divide-and_Conquer algorithm. */ static void dc_triangulate(CDT_state *cdt, SiteInfo *sites, int nsites) { int i, j, n; SymEdge *le, *re; /* Compress sites in place to eliminated verts that merge to others. */ i = 0; j = 0; while (j < nsites) { /* Invariante: sites[0..i-1] have non-merged verts from 0..(j-1) in them. */ sites[i] = sites[j++]; if (sites[i].v->merge_to_index < 0) { i++; } } n = i; if (n == 0) { return; } dc_tri(cdt, sites, 0, n, &le, &re); } /** * Do a Delaunay Triangulation of the points in cdt->vert_array. * This is only a first step in the Constrained Delaunay triangulation, * because it doesn't yet deal with the segment constraints. * The algorithm used is the Divide & Conquer algorithm from the * Guibas-Stolfi "Primitives for the Manipulation of General Subdivision * and the Computation of Voronoi Diagrams" paper. * The data structure here is similar to but not exactly the same as * the quad-edge structure described in that paper. * The incircle and ccw tests are done using Shewchuk's exact * primitives (see below), so that this routine is robust. * * As a preprocessing step, we want to merge all vertices that are * within cdt->epsilon of each other. This is accomplished by lexicographically * sorting the coordinates first (which is needed anyway for the D&C algorithm). * The CDTVerts with merge_to_index not equal to -1 are after this regarded * as having been merged into the vertex with the corresponding index. */ static void initial_triangulation(CDT_state *cdt) { int i, j, n; SiteInfo *sites; double *ico, *jco; double xend, yend, xcur; double epsilon = cdt->epsilon; double epsilon_squared = cdt->epsilon_squared; #ifdef SJF_WAY CDTEdge *e; CDTVert *va, *vb; #endif #ifdef DEBUG_CDT int dbg_level = 0; if (dbg_level > 0) { fprintf(stderr, "\nINITIAL TRIANGULATION\n\n"); } #endif /* First sort the vertices by lexicographic order of their * coordinates, breaking ties by putting earlier original-index * vertices first. */ n = cdt->vert_array_len; if (n <= 1) { return; } sites = MEM_malloc_arrayN(n, sizeof(SiteInfo), __func__); for (i = 0; i < n; i++) { sites[i].v = cdt->vert_array[i]; sites[i].orig_index = i; } qsort(sites, n, sizeof(SiteInfo), site_lexicographic_cmp); #ifdef DEBUG_CDT if (dbg_level > 0) { fprintf(stderr, "after sorting\n"); for (i = 0; i < n; i++) { fprintf(stderr, "%d: orig index: %d, (%f,%f)\n", i, sites[i].orig_index, F2(sites[i].v->co)); } } #endif /* Now dedup according to user-defined epsilon. * We will merge a vertex into an earlier-indexed vertex * that is within epsilon (Euclidean distance). * Merges may cascade. So we may end up merging two things * that are farther than epsilon by transitive merging. Oh well. * Assume that merges are rare, so use simple searches in the * lexicographic ordering - likely we will soon hit y's with * the same x that are farther away than epsilon, and then * skipping ahead to the next biggest x, are likely to soon * find one of those farther away than epsilon. */ for (i = 0; i < n - 1; i++) { ico = sites[i].v->co; /* Start j at next place that has both x and y coords within epsilon. */ xend = ico[0] + epsilon; yend = ico[1] + epsilon; j = i + 1; while (j < n) { jco = sites[j].v->co; if (jco[0] > xend) { break; /* No more j's to process. */ } else if (jco[1] > yend) { /* Get past any string of v's with the same x and too-big y. */ xcur = jco[0]; while (++j < n) { if (sites[j].v->co[0] > xcur) { break; } } BLI_assert(j == n || sites[j].v->co[0] > xcur); if (j == n) { break; } jco = sites[j].v->co; if (jco[0] > xend || jco[1] > yend) { break; } } /* When here, vertex i and j are within epsilon by box test. * The Euclidean distance test is stricter, so need to do it too, now. */ BLI_assert(j < n && jco[0] <= xend && jco[1] <= yend); if (len_squared_v2v2_db(ico, jco) <= epsilon_squared) { sites[j].v->merge_to_index = (sites[i].v->merge_to_index == -1) ? sites[i].orig_index : sites[i].v->merge_to_index; #ifdef DEBUG_CDT if (dbg_level > 0) { fprintf(stderr, "merged orig vert %d to %d\n", sites[j].orig_index, sites[j].v->merge_to_index); } #endif } j++; } } /* Now add non-dup vertices into triangulation in lexicographic order. */ dc_triangulate(cdt, sites, n); MEM_freeN(sites); } /** Use LinkNode linked list as stack of SymEdges, allocating from cdt->listpool. */ typedef LinkNode *Stack; BLI_INLINE void push(Stack *stack, SymEdge *se, CDT_state *cdt) { BLI_linklist_prepend_pool(stack, se, cdt->listpool); } BLI_INLINE SymEdge *pop(Stack *stack, CDT_state *cdt) { return (SymEdge *)BLI_linklist_pop_pool(stack, cdt->listpool); } BLI_INLINE bool is_empty(Stack *stack) { return *stack == NULL; } /** * Re-triangulates, assuring constrained delaunay condition, * the pseudo-polygon that cycles from se. * "pseudo" because a vertex may be repeated. * See Anglada paper, "An Improved incremental algorithm * for constructing restricted Delaunay triangulations". */ static void re_delaunay_triangulate(CDT_state *cdt, SymEdge *se) { SymEdge *ss, *first, *cse; CDTVert *a, *b, *c, *v; CDTEdge *ebc, *eca; int count; #ifdef DEBUG_CDT SymEdge *last; const int dbg_level = 0; #endif if (se->face == cdt->outer_face || sym(se)->face == cdt->outer_face) { return; } #ifdef DEBUG_CDT if (dbg_level > 0) { fprintf(stderr, "retriangulate"); dump_se_cycle(se, "poly ", 1000); } #endif /* 'se' is a diagonal just added, and it is base of area to retriangulate (face on its left) */ count = 1; for (ss = se->next; ss != se; ss = ss->next) { count++; } if (count <= 3) { #ifdef DEBUG_CDT if (dbg_level > 0) { fprintf(stderr, "nothing to do\n"); } #endif return; } /* First and last are the SymEdges whose verts are first and last off of base, * continuing from 'se'. */ first = se->next->next; /* We want to make a triangle with 'se' as base and some other c as 3rd vertex. */ a = se->vert; b = se->next->vert; c = first->vert; cse = first; #ifdef DEBUG_CDT last = prev(se); if (dbg_level > 1) { dump_se(first, "first"); dump_se(last, "last"); dump_v(a, "a"); dump_v(b, "b"); dump_v(c, "c"); } #endif for (ss = first->next; ss != se; ss = ss->next) { v = ss->vert; if (incircle(a->co, b->co, c->co, v->co) > 0.0) { c = v; cse = ss; #ifdef DEBUG_CDT if (dbg_level > 1) { dump_v(c, "new c "); } #endif } } /* Add diagonals necessary to make abc a triangle. */ #ifdef DEBUG_CDT if (dbg_level > 0) { fprintf(stderr, "make triangle abc exist where\n"); dump_v(a, " a"); dump_v(b, " b"); dump_v(c, " c"); } #endif ebc = NULL; eca = NULL; if (!exists_edge(b, c)) { ebc = add_diagonal(cdt, se->next, cse); #ifdef DEBUG_CDT if (dbg_level > 1) { fprintf(stderr, "added edge ebc\n"); dump_se(&ebc->symedges[0], " ebc"); } #endif } if (!exists_edge(c, a)) { eca = add_diagonal(cdt, cse, se); #ifdef DEBUG_CDT if (dbg_level > 1) { fprintf(stderr, "added edge eca\n"); dump_se(&eca->symedges[0], " eca"); } #endif } /* Now recurse. */ if (ebc) { re_delaunay_triangulate(cdt, &ebc->symedges[1]); } if (eca) { re_delaunay_triangulate(cdt, &eca->symedges[1]); } } /** * Add a constrained edge between v1 and v2 to cdt structure. * This may result in a number of #CDTEdges created, due to intersections * and partial overlaps with existing cdt vertices and edges. * Each created #CDTEdge will have input_id added to its input_ids list. * * If \a r_edges is not NULL, the #CDTEdges generated or found that go from * v1 to v2 are put into that linked list, in order. * * Assumes that #BLI_constrained_delaunay_get_output has not been called yet. */ static void add_edge_constraint( CDT_state *cdt, CDTVert *v1, CDTVert *v2, int input_id, LinkNode **r_edges) { CDTVert *va, *vb, *vc; SymEdge *vse1; #ifdef DEBUG_CDT SymEdge *vse2; #endif SymEdge *t, *tstart, *tout, *tnext; SymEdge *se; CDTEdge *edge; int isect; double orient1, orient2; int i, search_count; double curco[2]; double lambda; bool done, state_through_vert; LinkNodePair edge_list = {NULL, NULL}; typedef struct CrossData { double lambda; CDTVert *vert; SymEdge *in; SymEdge *out; } CrossData; CrossData cdata; CrossData *crossings = NULL; CrossData *cd; BLI_array_staticdeclare(crossings, 128); #ifdef DEBUG_CDT int dbg_level = 0; #endif /* Find path through structure from v1 to v2 and record how we got there in crossings. * In crossings array, each CrossData is populated as follows: * * If ray from previous node goes through a face, not along an edge: * * _ B * / |\ * - - | \ * prev........X \ * \ d | \C * -- | / * \ a| b/ * - - | / * \ A * * lambda = fraction of way along AB where X is. * vert = NULL initially, will later get new node that splits AB * in = a (SymEdge from A->B, whose face the ray goes through) * out = b (SymEdge from A->C, whose face the ray goes through next * * If the ray from the previous node goes directly to an existing vertex, say A * in the previous diagram, maybe along an existing edge like d in that diagram * but if prev had lambda !=0 then there may be no such edge d, then: * * lambda = 0 * vert = A * in = a * out = b * * crossings[0] will have in = NULL, and crossings[last] will have out = NULL. */ if (r_edges) { *r_edges = NULL; } vse1 = v1->symedge; #ifdef DEBUG_CDT if (dbg_level > 0) { vse2 = v2->symedge; if (dbg_level > 2) { // dump_cdt(cdt, "before insert_segment"); cdt_draw(cdt, "before insert segment"); } fprintf(stderr, "\ninsert_segment %d\n", input_id); dump_v(v1, " 1"); dump_v(v2, " 2"); if (dbg_level > 1) { dump_se(vse1, " se1"); dump_se(vse2, " se2"); } } #endif if (v1 == v2) { #ifdef DEBUG_CDT if (dbg_level > 0) { fprintf(stderr, "segment between same vertices, ignored\n"); } #endif return; } /* Quick check for special case where segment is one of edges from v1. */ tstart = t = vse1; done = false; 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); 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, "special one segment case\n"); dump_se(t, " "); } #endif done = true; break; } t = t->rot; } while (t != tstart); if (!done) { 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); } #endif 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, "found v2, so done\n"); } #endif 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; orient1 = orient2d(t->vert->co, va->co, v2->co); #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, "orient1=%g\n", orient1); } #endif if (orient1 == 0.0 && in_line(t->vert->co, va->co, v2->co)) { #ifdef DEBUG_CDT if (dbg_level > 1) { fprintf(stderr, "ray goes through va\n"); } #endif state_through_vert = true; tout = t; t = t->next; break; } else if (t->face != cdt->outer_face) { orient2 = orient2d(t->vert->co, vb->co, v2->co); #ifdef DEBUG_CDT if (dbg_level > 1) { fprintf(stderr, "orient2=%g\n", orient2); } #endif if (orient2 == 0.0 && in_line(t->vert->co, vb->co, v2->co)) { #ifdef DEBUG_CDT if (dbg_level > 1) { fprintf(stderr, "ray goes through vb\n"); } #endif state_through_vert = true; t = t->next->next; tout = sym(t); break; } else if (orient1 > 0.0 && orient2 < 0.0) { #ifdef DEBUG_CDT if (dbg_level > 1) { 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 > 1) { dump_se_cycle(t, "next rot tri", 4); } #endif } while (t != tstart); BLI_assert(tout != NULL); 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; /* Get curco; cdata should have data from last time through the loop still. */ if (cdata.lambda == 0.0) { copy_v2_v2_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 > 1) { fprintf(stderr, "through edge case, curco=(%f,%f)\n", F2(curco)); dump_v(va, " va"); dump_v(vb, " vb"); } #endif isect = isect_seg_seg_v2_lambda_mu_db(va->co, vb->co, curco, v2->co, &lambda, NULL); if (isect == ISECT_LINE_LINE_NONE || isect == ISECT_LINE_LINE_EXACT) { /* The orient tests say that there is an intersection between * va and vb, but the inexect isect routine has either put the * intersection exactly on one of the endpoints or just outside * one of them. * Or this is an exact intersection at one of the curco / v2 ends. * If lambda is outside of range, move the intersection to somewhere * just inside the segment. * Could also snap to an endpoint and redo this as a "through vert" * case, but the short edge will be cleaned up later and this seems * less risky to get into "impossible" cases. */ if (lambda <= 0.0) { lambda = 4.0 * (double)FLT_EPSILON; } else if (lambda >= 1.0) { lambda = 1.0 - 4.0 * (double)FLT_EPSILON; } } if (isect == ISECT_LINE_LINE_COLINEAR) { #ifdef DEBUG_CDT if (dbg_level > 1) { fprintf(stderr, "intersect is collinear, treating as through vert\n"); } dump_v(va, "va"); #endif state_through_vert = true; continue; } #ifdef DEBUG_CDT interp_v2_v2v2_db(curco, va->co, vb->co, lambda); if (dbg_level > 0) { fprintf(stderr, "intersect point at lambda=%.17g along va--vb\n", lambda); fprintf(stderr, "which is (%g,%g)\n", F2(curco)); if (dbg_level > 1) { 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); #ifdef DEBUG_CDT if (dbg_level > 0) { dump_se_cycle(tout, "next search tri", 4); dump_se(tout, "tout"); } #endif /* 'tout' is 'symedge' from 'va' to third vertex, 'vc'. */ BLI_assert(tout->vert == va); vc = tout->next->vert; orient1 = orient2d(curco, v2->co, vc->co); #ifdef DEBUG_CDT if (dbg_level > 1) { fprintf(stderr, "now searching with third vertex "); dump_v(vc, "vc"); fprintf(stderr, "orient2d(vcur, v2, vc) = %g\n", orient1); } #endif if (orient1 < 0.0) { /* vcur--v2 should intersect vb--vc. */ #ifdef DEBUG_CDT if (dbg_level > 1) { fprintf(stderr, "curco--v2 intersects vb--vc\n"); } #endif t = tout->next; state_through_vert = false; } else if (orient1 > 0.0) { /* 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 (orient1 == 0.0) { #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; } else { #ifdef DEBUG_CDT fprintf(stderr, "add_edge_constraint desperation search needed\n"); #endif } } if (++search_count > 1000000) { fprintf(stderr, "infinite loop? bailing out\n"); BLI_assert(0); /* Catch these while developing. */ break; } } } #ifdef DEBUG_CDT if (dbg_level > 0) { fprintf(stderr, "Crossing info gathered:\n"); for (i = 0; i < BLI_array_len(crossings); i++) { cd = &crossings[i]; fprintf(stderr, "%d:\n", i); if (cd->vert != NULL) { dump_v(cd->vert, " vert: "); } else { fprintf(stderr, " lambda=%f along in\n", cd->lambda); } if (cd->in) { dump_se(cd->in, " in: "); } if (cd->out) { dump_se(cd->out, " out: "); } } } #endif if (BLI_array_len(crossings) == 2) { /* For speed, handle special case of segment must have already been there. */ se = crossings[1].in; if (se->next->vert != v1) { se = prev(se); } BLI_assert(se->vert == v1 || se->next->vert == v1); #ifdef DEBUG_CDT if (dbg_level > 0) { fprintf(stderr, "segment already there: "); dump_se(se, ""); } #endif add_to_input_ids(&se->edge->input_ids, input_id, cdt); if (r_edges != NULL) { BLI_linklist_append_pool(&edge_list, se->edge, cdt->listpool); } } else { /* Insert all intersection points. */ for (i = 0; i < BLI_array_len(crossings); i++) { cd = &crossings[i]; if (cd->lambda != 0.0 && is_constrained_edge(cd->in->edge)) { edge = split_edge(cdt, cd->in, cd->lambda); cd->vert = edge->symedges[0].vert; #ifdef DEBUG_CDT if (dbg_level > 1) { fprintf(stderr, "insert vert for crossing %d: ", i); dump_v(cd->vert, "inserted"); } #endif } } /* Remove any crossed, non-intersected edges. */ for (i = 0; i < BLI_array_len(crossings); i++) { cd = &crossings[i]; if (cd->lambda != 0.0 && !is_constrained_edge(cd->in->edge)) { delete_edge(cdt, cd->in); #ifdef DEBUG_CDT if (dbg_level > 1) { fprintf(stderr, "delete edge for crossing %d\n", i); } #endif } } /* Insert segments for v1->v2. */ tstart = crossings[0].out; for (i = 1; i < BLI_array_len(crossings); i++) { cd = &crossings[i]; t = tnext = NULL; if (cd->lambda != 0.0) { if (is_constrained_edge(cd->in->edge)) { t = cd->vert->symedge; tnext = sym(t)->next; } } else if (cd->lambda == 0.0) { t = cd->in; tnext = cd->out; } if (t) { #ifdef DEBUG_CDT if (dbg_level > 1) { fprintf(stderr, "insert diagonal between\n"); dump_se(tstart, " "); dump_se(t, " "); dump_se_cycle(tstart, "tstart", 100); dump_se_cycle(t, "t", 100); } #endif if (tstart->next->vert == t->vert) { edge = tstart->edge; #ifdef DEBUG_CDT if (dbg_level > 1) { fprintf(stderr, "already there\n"); } #endif } else { edge = add_diagonal(cdt, tstart, t); } add_to_input_ids(&edge->input_ids, input_id, cdt); #ifdef DEBUG_CDT if (dbg_level > 1) { fprintf(stderr, "added\n"); } #endif if (r_edges != NULL) { BLI_linklist_append_pool(&edge_list, edge, cdt->listpool); } /* Now retriangulate upper and lower gaps. */ re_delaunay_triangulate(cdt, &edge->symedges[0]); re_delaunay_triangulate(cdt, &edge->symedges[1]); } if (i < BLI_array_len(crossings) - 1) { if (tnext != NULL) { tstart = tnext; #ifdef DEBUG_CDT if (dbg_level > 1) { fprintf(stderr, "now tstart = "); dump_se(tstart, ""); } #endif } } } } if (r_edges) { *r_edges = edge_list.list; } BLI_array_free(crossings); } /** * Add face_id to the input_ids lists of all #CDTFace's on the interior of the input face with that * id. face_symedge is on edge of the boundary of the input face, with assumption that interior is * on the left of that SymEdge. * * The algorithm is: starting from the #CDTFace for face_symedge, add the face_id and then * process all adjacent faces where the adjacency isn't across an edge that was a constraint added * for the boundary of the input face. * fedge_start..fedge_end is the inclusive range of edge input ids that are for the given face. * * Note: if the input face is not CCW oriented, we'll be labeling the outside, not the inside. * Note 2: if the boundary has self-crossings, this method will arbitrarily pick one of the * contiguous set of faces enclosed by parts of the boundary, leaving the other such untagged. This * may be a feature instead of a bug if the first contiguous section is most of the face and the * others are tiny self-crossing triangles at some parts of the boundary. On the other hand, if * decide we want to handle these in full generality, then will need a more complicated algorithm * (using "inside" tests and a parity rule) to decide on the interior. */ static void add_face_ids( CDT_state *cdt, SymEdge *face_symedge, int face_id, int fedge_start, int fedge_end) { Stack stack; SymEdge *se, *se_start, *se_sym; CDTFace *face, *face_other; int visit; /* Can't loop forever since eventually would visit every face. */ cdt->visit_count++; visit = cdt->visit_count; stack = NULL; push(&stack, face_symedge, cdt); while (!is_empty(&stack)) { se = pop(&stack, cdt); face = se->face; if (face->visit_index == visit) { continue; } face->visit_index = visit; add_to_input_ids(&face->input_ids, face_id, cdt); se_start = se; for (se = se->next; se != se_start; se = se->next) { if (!id_range_in_list(se->edge->input_ids, fedge_start, fedge_end)) { se_sym = sym(se); face_other = se_sym->face; if (face_other->visit_index != visit) { push(&stack, se_sym, cdt); } } } } } /* Delete_edge but try not to mess up outer face. * Also faces have symedges now, so make sure not * to mess those up either. */ static void dissolve_symedge(CDT_state *cdt, SymEdge *se) { SymEdge *symse = sym(se); if (symse->face == cdt->outer_face) { se = sym(se); symse = sym(se); } if (cdt->outer_face->symedge == se || cdt->outer_face->symedge == symse) { /* Advancing by 2 to get past possible 'sym(se)'. */ if (se->next->next == se) { cdt->outer_face->symedge = NULL; } else { cdt->outer_face->symedge = se->next->next; } } else { if (se->face->symedge == se) { se->face->symedge = se->next; } if (symse->face->symedge == symse) { symse->face->symedge = symse->next; } } delete_edge(cdt, se); } /* Return true if we can merge se's vert into se->next's vert * without making the area of any new triangle formed by doing * that into a zero or negative area triangle.*/ static bool can_collapse(const SymEdge *se) { SymEdge *loop_se; const double *co = se->next->vert->co; for (loop_se = se->rot; loop_se != se && loop_se->rot != se; loop_se = loop_se->rot) { if (orient2d(co, loop_se->next->vert->co, loop_se->rot->next->vert->co) <= 0.0) { return false; } } return true; } /* * Merge one end of e onto the other, fixing up surrounding faces. * * General situation looks something like: * * c-----e * / \ / \ * / \ / \ * a------b-----f * \ / \ / * \ / \ / * d-----g * * where ab is the tiny edge. We want to merge a and b and delete edge ab. * We don't want to change the coordinates of input vertices [We could revisit this * in the future, as API def doesn't prohibit this, but callers will appreciate if they * don't change.] * Sometimes the collapse shouldn't happen because the triangles formed by the changed * edges may end up with zero or negative area (see can_collapse, above). * So don't choose a collapse direction that is not allowed or one that has an original vertex * as origin and a non-original vertex as destination. * If both collapse directions are allowed by that rule, pick the one with the lower original * index. * * After merging, the faces abc and adb disappear (if they are not the outer face). * Suppose we merge b onto a. * Then edges cb and db are deleted. Face cbe becomes cae and face bdg becomes adg. * Any other faces attached to b now have a in their place. * We can do this by rotating edges round b, replacing their vert references with a. * Similar statements can be made about what happens when a merges into b; * in code below we'll swap a and b to make above lettering work for a b->a merge. * Return the vert at the collapsed edge, if a collapse happens. */ static CDTVert *collapse_tiny_edge(CDT_state *cdt, CDTEdge *e) { CDTVert *va, *vb; SymEdge *ab_se, *ba_se, *bd_se, *bc_se, *ad_se, *ac_se; SymEdge *bg_se, *be_se, *se, *gb_se, *ca_se; bool can_collapse_a_to_b, can_collapse_b_to_a; #ifdef DEBUG_CDT int dbg_level = 0; #endif ab_se = &e->symedges[0]; ba_se = &e->symedges[1]; va = ab_se->vert; vb = ba_se->vert; #ifdef DEBUG_CDT if (dbg_level > 0) { fprintf(stderr, "\ncollapse_tiny_edge\n"); dump_se(&e->symedges[0], "tiny edge"); fprintf(stderr, "a = [%d], b = [%d]\n", va->index, vb->index); validate_cdt(cdt, true, false, true); } #endif can_collapse_a_to_b = can_collapse(ab_se); can_collapse_b_to_a = can_collapse(ba_se); /* Now swap a and b if necessary and possible, so that from this point on we are collapsing b to * a. */ if (va->index > vb->index || !can_collapse_b_to_a) { if (can_collapse_a_to_b && !(is_original_vert(va, cdt) && !is_original_vert(vb, cdt))) { SWAP(CDTVert *, va, vb); ab_se = &e->symedges[1]; ba_se = &e->symedges[0]; #ifdef DEBUG_CDT if (dbg_level > 0) { fprintf(stderr, "swapped a and b\n"); } #endif } else { /* Neither collapse direction is OK. */ #ifdef DEBUG_CDT if (dbg_level > 0) { fprintf(stderr, "neither collapse direction ok\n"); } #endif return NULL; } } bc_se = ab_se->next; bd_se = ba_se->rot; if (bd_se == ba_se) { /* Shouldn't happen. Wire edge in outer face. */ fprintf(stderr, "unexpected wire edge\n"); return NULL; } vb->merge_to_index = va->merge_to_index == -1 ? va->index : va->merge_to_index; vb->symedge = NULL; #ifdef DEBUG_CDT if (dbg_level > 0) { fprintf(stderr, "vb = v[%d] merges to va = v[%d], vb->merge_to_index=%d\n", vb->index, va->index, vb->merge_to_index); } #endif /* First fix the vertex of intermediate triangles, like bgf. */ for (se = bd_se->rot; se != bc_se; se = se->rot) { #ifdef DEBUG_CDT if (dbg_level > 0) { dump_se(se, "intermediate tri edge, setting vert to va"); } #endif se->vert = va; } ad_se = sym(sym(bd_se)->rot); ca_se = bc_se->next; ac_se = sym(ca_se); if (bd_se->rot != bc_se) { bg_se = bd_se->rot; be_se = sym(bc_se)->next; gb_se = sym(bg_se); } else { bg_se = NULL; be_se = NULL; } #ifdef DEBUG_CDT if (dbg_level > 0) { fprintf(stderr, "delete bd, inputs to ad\n"); dump_se(bd_se, " bd"); dump_se(ad_se, " ad"); fprintf(stderr, "delete bc, inputs to ac\n"); dump_se(bc_se, " bc"); dump_se(ac_se, " ac"); fprintf(stderr, "delete ab\n"); dump_se(ab_se, " ab"); if (bg_se != NULL) { fprintf(stderr, "fix up bg, be\n"); dump_se(bg_se, " bg"); dump_se(be_se, " be"); } } #endif add_list_to_input_ids(&ad_se->edge->input_ids, bd_se->edge->input_ids, cdt); delete_edge(cdt, bd_se); add_list_to_input_ids(&ac_se->edge->input_ids, bc_se->edge->input_ids, cdt); delete_edge(cdt, sym(bc_se)); /* At this point we have this: * * c-----e * / / \ * / / \ * a------b-----f * \ \ / * \ \ / * d-----g * * Or, if there is not bg_se and be_se, like this: * * c * / * / * a------b * \ * \ * d * * (But we've already changed the vert field for bg, bf, ..., be to be va.) */ if (bg_se != NULL) { gb_se->next = ad_se; ad_se->rot = bg_se; ca_se->next = be_se; be_se->rot = ac_se; bg_se->vert = va; be_se->vert = va; } else { ca_se->next = ad_se; ad_se->rot = ac_se; } /* Don't use delete_edge as it changes too much. */ ab_se->next = ab_se->rot = NULL; ba_se->next = ba_se->rot = NULL; if (va->symedge == ab_se) { va->symedge = ac_se; } return va; } /* * Check to see if e is tiny (length <= epsilon) and queue it if so. */ static void maybe_enqueue_small_feature(CDT_state *cdt, CDTEdge *e, LinkNodePair *tiny_edge_queue) { SymEdge *se, *sesym; #ifdef DEBUG_CDT int dbg_level = 0; if (dbg_level > 0) { fprintf(stderr, "\nmaybe_enqueue_small_features\n"); dump_se(&e->symedges[0], " se0"); } #endif if (is_deleted_edge(e) || e->in_queue) { #ifdef DEBUG_CDT if (dbg_level > 0) { fprintf(stderr, "returning because of e conditions\n"); } #endif return; } se = &e->symedges[0]; sesym = &e->symedges[1]; if (len_squared_v2v2_db(se->vert->co, sesym->vert->co) <= cdt->epsilon_squared) { BLI_linklist_append_pool(tiny_edge_queue, e, cdt->listpool); e->in_queue = true; #ifdef DEBUG_CDT if (dbg_level > 0) { fprintf(stderr, "Queue tiny edge\n"); } #endif } } /* Consider all edges in rot ring around v for possible enqueing as small features .*/ static void maybe_enqueue_small_features(CDT_state *cdt, CDTVert *v, LinkNodePair *tiny_edge_queue) { SymEdge *se, *se_start; se = se_start = v->symedge; if (!se_start) { return; } do { maybe_enqueue_small_feature(cdt, se->edge, tiny_edge_queue); } while ((se = se->rot) != se_start); } /* Collapse small edges (length <= epsilon) until no more exist. */ static void remove_small_features(CDT_state *cdt) { double epsilon = cdt->epsilon; LinkNodePair tiny_edge_queue = {NULL, NULL}; BLI_mempool *pool = cdt->listpool; LinkNode *ln; CDTEdge *e; CDTVert *v_change; #ifdef DEBUG_CDT int dbg_level = 0; if (dbg_level > 0) { fprintf(stderr, "\nREMOVE_SMALL_FEATURES, epsilon=%g\n", epsilon); } #endif if (epsilon == 0.0) { return; } for (ln = cdt->edges; ln; ln = ln->next) { e = (CDTEdge *)ln->link; maybe_enqueue_small_feature(cdt, e, &tiny_edge_queue); } while (tiny_edge_queue.list != NULL) { e = (CDTEdge *)BLI_linklist_pop_pool(&tiny_edge_queue.list, pool); if (tiny_edge_queue.list == NULL) { tiny_edge_queue.last_node = NULL; } e->in_queue = false; if (is_deleted_edge(e)) { continue; } #ifdef DEBUG_CDT if (dbg_level > 0) { fprintf(stderr, "collapse tiny edge\n"); dump_se(&e->symedges[0], ""); } #endif v_change = collapse_tiny_edge(cdt, e); if (v_change) { maybe_enqueue_small_features(cdt, v_change, &tiny_edge_queue); } } } /* Remove all non-constraint edges. */ static void remove_non_constraint_edges(CDT_state *cdt) { LinkNode *ln; CDTEdge *e; SymEdge *se; for (ln = cdt->edges; ln; ln = ln->next) { e = (CDTEdge *)ln->link; se = &e->symedges[0]; if (!is_deleted_edge(e) && !is_constrained_edge(e)) { dissolve_symedge(cdt, se); } } } /* * Remove the non-constraint edges, but leave enough of them so that all of the * faces that would be bmesh faces (that is, the faces that have some input representative) * are valid: they can't have holes, they can't have repeated vertices, and they can't have * repeated edges. * * Not essential, but to make the result look more aesthetically nice, * remove the edges in order of decreasing length, so that it is more likely that the * final remaining support edges are short, and therefore likely to make a fairly * direct path from an outer face to an inner hole face. */ /* For sorting edges by decreasing length (squared). */ struct EdgeToSort { double len_squared; CDTEdge *e; }; static int edge_to_sort_cmp(const void *a, const void *b) { const struct EdgeToSort *e1 = a; const struct EdgeToSort *e2 = b; if (e1->len_squared > e2->len_squared) { return -1; } else if (e1->len_squared < e2->len_squared) { return 1; } return 0; } static void remove_non_constraint_edges_leave_valid_bmesh(CDT_state *cdt) { LinkNode *ln; CDTEdge *e; SymEdge *se, *se2; CDTFace *fleft, *fright; bool dissolve; size_t nedges; int i, ndissolvable; const double *co1, *co2; struct EdgeToSort *sorted_edges; nedges = 0; for (ln = cdt->edges; ln; ln = ln->next) { nedges++; } if (nedges == 0) { return; } sorted_edges = BLI_memarena_alloc(cdt->arena, nedges * sizeof(*sorted_edges)); i = 0; for (ln = cdt->edges; ln; ln = ln->next) { e = (CDTEdge *)ln->link; if (!is_deleted_edge(e) && !is_constrained_edge(e)) { sorted_edges[i].e = e; co1 = e->symedges[0].vert->co; co2 = e->symedges[1].vert->co; sorted_edges[i].len_squared = len_squared_v2v2_db(co1, co2); i++; } } ndissolvable = i; qsort(sorted_edges, ndissolvable, sizeof(*sorted_edges), edge_to_sort_cmp); for (i = 0; i < ndissolvable; i++) { e = sorted_edges[i].e; se = &e->symedges[0]; dissolve = true; if (true /*!edge_touches_frame(e)*/) { fleft = se->face; fright = sym(se)->face; if (fleft != cdt->outer_face && fright != cdt->outer_face && (fleft->input_ids != NULL || fright->input_ids != NULL)) { /* Is there another symedge with same left and right faces? * Or is there a vertex not part of e touching the same left and right faces? */ for (se2 = se->next; dissolve && se2 != se; se2 = se2->next) { if (sym(se2)->face == fright || (se2->vert != se->next->vert && vert_touches_face(se2->vert, fright))) { dissolve = false; } } } } if (dissolve) { dissolve_symedge(cdt, se); } } } static void remove_outer_edges_until_constraints(CDT_state *cdt) { LinkNode *fstack = NULL; SymEdge *se, *se_start; CDTFace *f, *fsym; int visit = ++cdt->visit_count; #ifdef DEBUG_CDT int dbg_level = 0; if (dbg_level > 0) { fprintf(stderr, "remove_outer_edges_until_constraints\n"); } #endif cdt->outer_face->visit_index = visit; /* Walk around outer face, adding faces on other side of dissolvable edges to stack. */ se_start = se = cdt->outer_face->symedge; do { if (!is_constrained_edge(se->edge)) { fsym = sym(se)->face; if (fsym->visit_index != visit) { #ifdef DEBUG_CDT if (dbg_level > 0) { fprintf(stderr, "pushing f=%p from symedge ", fsym); dump_se(se, "an outer edge"); } #endif BLI_linklist_prepend_pool(&fstack, fsym, cdt->listpool); } } } while ((se = se->next) != se_start); while (fstack != NULL) { LinkNode *to_dissolve = NULL; bool dissolvable; f = (CDTFace *)BLI_linklist_pop_pool(&fstack, cdt->listpool); if (f->visit_index == visit) { #ifdef DEBUG_CDT if (dbg_level > 0) { fprintf(stderr, "skipping f=%p, already visited\n", f); } #endif continue; } BLI_assert(f != cdt->outer_face); #ifdef DEBUG_CDT if (dbg_level > 0) { fprintf(stderr, "top of loop, f=%p\n", f); dump_se_cycle(f->symedge, "visit", 10000); if (dbg_level > 1) { dump_cdt(cdt, "cdt at top of loop"); cdt_draw(cdt, "top of dissolve loop"); } } #endif f->visit_index = visit; se_start = se = f->symedge; do { dissolvable = !is_constrained_edge(se->edge); #ifdef DEBUG_CDT if (dbg_level > 1) { dump_se(se, "edge in f"); fprintf(stderr, " dissolvable=%d\n", dissolvable); } #endif if (dissolvable) { fsym = sym(se)->face; #ifdef DEBUG_CDT if (dbg_level > 1) { dump_se_cycle(fsym->symedge, "fsym", 10000); fprintf(stderr, " visited=%d\n", fsym->visit_index == visit); } #endif if (fsym->visit_index != visit) { #ifdef DEBUG_CDT if (dbg_level > 0) { fprintf(stderr, "pushing face %p\n", fsym); dump_se_cycle(fsym->symedge, "pushed", 10000); } #endif BLI_linklist_prepend_pool(&fstack, fsym, cdt->listpool); } else { BLI_linklist_prepend_pool(&to_dissolve, se, cdt->listpool); } } se = se->next; } while (se != se_start); while (to_dissolve != NULL) { se = (SymEdge *)BLI_linklist_pop_pool(&to_dissolve, cdt->listpool); if (se->next != NULL) { dissolve_symedge(cdt, se); } } } } /** * Remove edges and merge faces to get desired output, as per options. * \note the cdt cannot be further changed after this. */ static void prepare_cdt_for_output(CDT_state *cdt, const CDT_output_type output_type) { CDTFace *f; CDTEdge *e; LinkNode *ln; cdt->output_prepared = true; if (cdt->edges == NULL) { return; } /* Make sure all non-deleted faces have a symedge. */ for (ln = cdt->edges; ln; ln = ln->next) { e = (CDTEdge *)ln->link; if (!is_deleted_edge(e)) { if (e->symedges[0].face->symedge == NULL) { e->symedges[0].face->symedge = &e->symedges[0]; } if (e->symedges[1].face->symedge == NULL) { e->symedges[1].face->symedge = &e->symedges[1]; } } } #ifdef DEBUG_CDT /* All non-deleted faces should have a symedge now. */ for (ln = cdt->faces; ln; ln = ln->next) { f = (CDTFace *)ln->link; if (!f->deleted) { BLI_assert(f->symedge != NULL); } } #else UNUSED_VARS(f); #endif if (output_type == CDT_CONSTRAINTS) { remove_non_constraint_edges(cdt); } else if (output_type == CDT_CONSTRAINTS_VALID_BMESH) { remove_non_constraint_edges_leave_valid_bmesh(cdt); } else if (output_type == CDT_INSIDE) { remove_outer_edges_until_constraints(cdt); } } static CDT_result *cdt_get_output(CDT_state *cdt, const CDT_input *input, const CDT_output_type output_type) { int i, j, nv, ne, nf, faces_len_total; int orig_map_size, orig_map_index; int *vert_to_output_map; CDT_result *result; CDTVert *v; LinkNode *lne, *lnf, *ln; SymEdge *se, *se_start; CDTEdge *e; CDTFace *f; prepare_cdt_for_output(cdt, output_type); result = (CDT_result *)MEM_callocN(sizeof(*result), __func__); if (cdt->vert_array_len == 0) { return result; } /* All verts without a merge_to_index will be output. * vert_to_output_map[i] will hold the output vertex index * corresponding to the vert in position i in cdt->vert_array. * Since merging picked the leftmost-lowermost representative, * that is not necessarily the same as the vertex with the lowest original * index (i.e., index in cdt->vert_array), so we need two passes: * one to get the non-merged-to vertices in vert_to_output_map, * and a second to put in the merge targets for merged-to vertices. */ vert_to_output_map = BLI_memarena_alloc(cdt->arena, (size_t)cdt->vert_array_len * sizeof(int *)); nv = 0; for (i = 0; i < cdt->vert_array_len; i++) { v = cdt->vert_array[i]; if (v->merge_to_index == -1) { vert_to_output_map[i] = nv; nv++; } } if (nv <= 0) { return result; } if (nv < cdt->vert_array_len) { for (i = 0; i < input->verts_len; i++) { v = cdt->vert_array[i]; if (v->merge_to_index != -1) { add_to_input_ids(&cdt->vert_array[v->merge_to_index]->input_ids, i, cdt); vert_to_output_map[i] = vert_to_output_map[v->merge_to_index]; } } } result->verts_len = nv; result->vert_coords = MEM_malloc_arrayN(nv, sizeof(result->vert_coords[0]), __func__); /* Make the vertex "orig" map arrays, mapping output verts to lists of input ones. */ orig_map_size = 0; for (i = 0; i < cdt->vert_array_len; i++) { if (cdt->vert_array[i]->merge_to_index == -1) { orig_map_size += 1 + BLI_linklist_count(cdt->vert_array[i]->input_ids); } } result->verts_orig_len_table = MEM_malloc_arrayN(nv, sizeof(int), __func__); result->verts_orig_start_table = MEM_malloc_arrayN(nv, sizeof(int), __func__); result->verts_orig = MEM_malloc_arrayN(orig_map_size, sizeof(int), __func__); orig_map_index = 0; i = 0; for (j = 0; j < cdt->vert_array_len; j++) { v = cdt->vert_array[j]; if (v->merge_to_index == -1) { result->vert_coords[i][0] = (float)v->co[0]; result->vert_coords[i][1] = (float)v->co[1]; result->verts_orig_start_table[i] = orig_map_index; result->verts_orig[orig_map_index++] = j; for (ln = v->input_ids; ln; ln = ln->next) { result->verts_orig[orig_map_index++] = POINTER_AS_INT(ln->link); } result->verts_orig_len_table[i] = orig_map_index - result->verts_orig_start_table[i]; i++; } } ne = 0; orig_map_size = 0; for (ln = cdt->edges; ln; ln = ln->next) { e = (CDTEdge *)ln->link; if (!is_deleted_edge(e)) { ne++; if (e->input_ids) { orig_map_size += BLI_linklist_count(e->input_ids); } } } if (ne != 0) { result->edges_len = ne; result->face_edge_offset = cdt->face_edge_offset; result->edges = MEM_malloc_arrayN(ne, sizeof(result->edges[0]), __func__); result->edges_orig_len_table = MEM_malloc_arrayN(ne, sizeof(int), __func__); result->edges_orig_start_table = MEM_malloc_arrayN(ne, sizeof(int), __func__); if (orig_map_size > 0) { result->edges_orig = MEM_malloc_arrayN(orig_map_size, sizeof(int), __func__); } orig_map_index = 0; i = 0; for (lne = cdt->edges; lne; lne = lne->next) { e = (CDTEdge *)lne->link; if (!is_deleted_edge(e)) { result->edges[i][0] = vert_to_output_map[e->symedges[0].vert->index]; result->edges[i][1] = vert_to_output_map[e->symedges[1].vert->index]; result->edges_orig_start_table[i] = orig_map_index; for (ln = e->input_ids; ln; ln = ln->next) { result->edges_orig[orig_map_index++] = POINTER_AS_INT(ln->link); } result->edges_orig_len_table[i] = orig_map_index - result->edges_orig_start_table[i]; i++; } } } nf = 0; faces_len_total = 0; orig_map_size = 0; for (ln = cdt->faces; ln; ln = ln->next) { f = (CDTFace *)ln->link; if (!f->deleted && f != cdt->outer_face) { nf++; se = se_start = f->symedge; BLI_assert(se != NULL); do { faces_len_total++; se = se->next; } while (se != se_start); if (f->input_ids) { orig_map_size += BLI_linklist_count(f->input_ids); } } } if (nf != 0) { result->faces_len = nf; result->faces_len_table = MEM_malloc_arrayN(nf, sizeof(int), __func__); result->faces_start_table = MEM_malloc_arrayN(nf, sizeof(int), __func__); result->faces = MEM_malloc_arrayN(faces_len_total, sizeof(int), __func__); result->faces_orig_len_table = MEM_malloc_arrayN(nf, sizeof(int), __func__); result->faces_orig_start_table = MEM_malloc_arrayN(nf, sizeof(int), __func__); if (orig_map_size > 0) { result->faces_orig = MEM_malloc_arrayN(orig_map_size, sizeof(int), __func__); } orig_map_index = 0; i = 0; j = 0; for (lnf = cdt->faces; lnf; lnf = lnf->next) { f = (CDTFace *)lnf->link; if (!f->deleted && f != cdt->outer_face) { result->faces_start_table[i] = j; se = se_start = f->symedge; do { result->faces[j++] = se->vert->index; se = se->next; } while (se != se_start); result->faces_len_table[i] = j - result->faces_start_table[i]; result->faces_orig_start_table[i] = orig_map_index; for (ln = f->input_ids; ln; ln = ln->next) { result->faces_orig[orig_map_index++] = POINTER_AS_INT(ln->link); } result->faces_orig_len_table[i] = orig_map_index - result->faces_orig_start_table[i]; i++; } } } return result; } /** * Calculate the Constrained Delaunay Triangulation of the 2d elements given in \a input. * * A Delaunay triangulation of a set of vertices is a triangulation where no triangle in the * triangulation has a circumcircle that strictly contains another vertex. Delaunay triangulations * are avoid long skinny triangles: they maximize the minimum angle of all triangles in the * triangulation. * * A Constrained Delaunay Triangulation adds the requirement that user-provided line segments must * appear as edges in the output (perhaps divided into several sub-segments). It is not required * that the input edges be non-intersecting: this routine will calculate the intersections. This * means that besides triangulating, this routine is also useful for general and robust 2d edge and * face intersection. * * This routine also takes an epsilon parameter in the \a input. Input vertices closer than epsilon * will be merged, and we collapse tiny edges (less than epsilon length) and skinny triangles * (having an altitude of less than epsilon). * * The current initial Deluanay triangulation algorithm is the Guibas-Stolfi Divide and Conquer * algorithm (see "Primitives for the Manipulation of General Subdivisions and the Computation of * Voronoi Diagrams"). and uses Shewchuk's exact predicates to issues where numeric errors cause * inconsistent geometric judgements. This is followed by inserting edge constraints (including the * edges implied by faces) using the algorithms discussed in "Fully Dynamic Constrained Delaunay * Triangulations" by Kallmann, Bieri, and Thalmann. * * \param input: points to a CDT_input struct which contains the vertices, edges, and faces to be * triangulated. \param output_type: specifies which edges to remove after doing the triangulation. * \return A pointer to an allocated CDT_result struct, which describes the triangulation in terms * of vertices, edges, and faces, and also has tables to map output elements back to input * elements. The caller must use BLI_delaunay_2d_cdt_free() on the result when done with it. * * See the header file BLI_delaunay_2d.h for details of the CDT_input and CDT_result structs and * the CDT_output_type enum. */ CDT_result *BLI_delaunay_2d_cdt_calc(const CDT_input *input, const CDT_output_type output_type) { int nv = input->verts_len; int ne = input->edges_len; int nf = input->faces_len; int i, iv1, iv2, f, fedge_start, fedge_end; CDT_state *cdt; CDTVert *v1, *v2; CDTEdge *face_edge; SymEdge *face_symedge; LinkNode *edge_list; CDT_result *result; static bool called_exactinit = false; #ifdef DEBUG_CDT int dbg_level = 0; #endif /* The exact orientation and incircle primitives need a one-time initialization of certain * constants. */ if (!called_exactinit) { exactinit(); called_exactinit = true; } #ifdef DEBUG_CDT if (dbg_level > 0) { fprintf(stderr, "\n\nCDT CALC, nv=%d, ne=%d, nf=%d, eps=%g\n", input->verts_len, input->edges_len, input->faces_len, input->epsilon); } if (dbg_level == -1) { write_cdt_input_to_file(input); } #endif if ((nv > 0 && input->vert_coords == NULL) || (ne > 0 && input->edges == NULL) || (nf > 0 && (input->faces == NULL || input->faces_start_table == NULL || input->faces_len_table == NULL))) { #ifdef DEBUG_CDT fprintf(stderr, "invalid input: unexpected NULL array(s)\n"); #endif return NULL; } cdt = new_cdt_init(input); initial_triangulation(cdt); #ifdef DEBUG_CDT if (dbg_level > 0) { validate_cdt(cdt, true, false, false); } #endif for (i = 0; i < ne; i++) { iv1 = input->edges[i][0]; iv2 = input->edges[i][1]; if (iv1 < 0 || iv1 >= nv || iv2 < 0 || iv2 >= nv) { #ifdef DEBUG_CDT fprintf(stderr, "edge indices for e%d not valid: v1=%d, v2=%d, nv=%d\n", i, iv1, iv2, nv); #endif continue; } v1 = cdt->vert_array[iv1]; v2 = cdt->vert_array[iv2]; if (v1->merge_to_index != -1) { v1 = cdt->vert_array[v1->merge_to_index]; } if (v2->merge_to_index != -1) { v2 = cdt->vert_array[v2->merge_to_index]; } add_edge_constraint(cdt, v1, v2, i, NULL); #ifdef DEBUG_CDT if (dbg_level > 3) { char namebuf[60]; sprintf(namebuf, "after edge constraint %d = (%d,%d)\n", i, iv1, iv2); cdt_draw(cdt, namebuf); // dump_cdt(cdt, namebuf); validate_cdt(cdt, true, true, false); } #endif } cdt->face_edge_offset = ne; for (f = 0; f < nf; f++) { int flen = input->faces_len_table[f]; int fstart = input->faces_start_table[f]; if (flen <= 2) { #ifdef DEBUG_CDT fprintf(stderr, "face %d has length %d; ignored\n", f, flen); #endif continue; } for (i = 0; i < flen; i++) { int face_edge_id = cdt->face_edge_offset + fstart + i; iv1 = input->faces[fstart + i]; iv2 = input->faces[fstart + ((i + 1) % flen)]; if (iv1 < 0 || iv1 >= nv || iv2 < 0 || iv2 >= nv) { #ifdef DEBUG_CDT fprintf(stderr, "face indices not valid: f=%d, iv1=%d, iv2=%d, nv=%d\n", f, iv1, iv2, nv); #endif continue; } v1 = cdt->vert_array[iv1]; v2 = cdt->vert_array[iv2]; if (v1->merge_to_index != -1) { v1 = cdt->vert_array[v1->merge_to_index]; } if (v2->merge_to_index != -1) { v2 = cdt->vert_array[v2->merge_to_index]; } add_edge_constraint(cdt, v1, v2, face_edge_id, &edge_list); #ifdef DEBUG_CDT if (dbg_level > 2) { fprintf(stderr, "edges for edge %d:\n", i); for (LinkNode *ln = edge_list; ln; ln = ln->next) { CDTEdge *cdt_e = (CDTEdge *)ln->link; fprintf(stderr, " (%.2f,%.2f)->(%.2f,%.2f)\n", F2(cdt_e->symedges[0].vert->co), F2(cdt_e->symedges[1].vert->co)); } } if (dbg_level > 2) { cdt_draw(cdt, "after a face edge"); if (dbg_level > 3) { dump_cdt(cdt, "after a face edge"); } validate_cdt(cdt, true, true, false); } #endif if (i == 0) { face_edge = (CDTEdge *)edge_list->link; face_symedge = &face_edge->symedges[0]; if (face_symedge->vert != v1) { face_symedge = &face_edge->symedges[1]; BLI_assert(face_symedge->vert == v1); } } BLI_linklist_free_pool(edge_list, NULL, cdt->listpool); } fedge_start = cdt->face_edge_offset + fstart; fedge_end = fedge_start + flen - 1; add_face_ids(cdt, face_symedge, f, fedge_start, fedge_end); } #ifdef DEBUG_CDT if (dbg_level > 0) { validate_cdt(cdt, true, true, false); } if (dbg_level > 1) { cdt_draw(cdt, "after adding edges and faces"); if (dbg_level > 2) { dump_cdt(cdt, "after adding edges and faces"); } } #endif if (cdt->epsilon > 0.0) { remove_small_features(cdt); #ifdef DEBUG_CDT if (dbg_level > 2) { cdt_draw(cdt, "after collapse skinny triangles\n"); if (dbg_level > 3) { dump_cdt(cdt, "after collapse skinny triangles\n"); } } #endif } result = cdt_get_output(cdt, input, output_type); #ifdef DEBUG_CDT if (dbg_level > 0) { cdt_draw(cdt, "final"); } #endif new_cdt_free(cdt); return result; } void BLI_delaunay_2d_cdt_free(CDT_result *result) { if (result == NULL) { return; } if (result->vert_coords) { MEM_freeN(result->vert_coords); } if (result->edges) { MEM_freeN(result->edges); } if (result->faces) { MEM_freeN(result->faces); } if (result->faces_start_table) { MEM_freeN(result->faces_start_table); } if (result->faces_len_table) { MEM_freeN(result->faces_len_table); } if (result->verts_orig) { MEM_freeN(result->verts_orig); } if (result->verts_orig_start_table) { MEM_freeN(result->verts_orig_start_table); } if (result->verts_orig_len_table) { MEM_freeN(result->verts_orig_len_table); } if (result->edges_orig) { MEM_freeN(result->edges_orig); } if (result->edges_orig_start_table) { MEM_freeN(result->edges_orig_start_table); } if (result->edges_orig_len_table) { MEM_freeN(result->edges_orig_len_table); } if (result->faces_orig) { MEM_freeN(result->faces_orig); } if (result->faces_orig_start_table) { MEM_freeN(result->faces_orig_start_table); } if (result->faces_orig_len_table) { MEM_freeN(result->faces_orig_len_table); } MEM_freeN(result); } #ifdef DEBUG_CDT ATTU static const char *vertname(const CDTVert *v) { static char vertnamebuf[20]; sprintf(vertnamebuf, "[%d]", v->index); return vertnamebuf; } ATTU static const char *sename(const SymEdge *se) { static char senamebuf[20]; sprintf(senamebuf, "{%x}", (POINTER_AS_UINT(se)) & 0xFFFF); return senamebuf; } static void dump_v(const CDTVert *v, const char *lab) { 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))", lab, vertname(se->vert), F2(se->vert->co), F2(se->next->vert->co)); fprintf(stderr, "%s\n", vertname(se->next->vert)); } 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) { int count = 0; const SymEdge *s = se; fprintf(stderr, "%s:\n", lab); do { dump_se(s, " "); s = s->next; count++; } while (s != se && count < limit); if (count == limit) { fprintf(stderr, " limit hit without cycle!\n"); } } static void dump_id_list(const LinkNode *id_list, const char *lab) { const LinkNode *ln; if (!id_list) { return; } fprintf(stderr, "%s", lab); for (ln = id_list; ln; ln = ln->next) { fprintf(stderr, "%d%c", POINTER_AS_INT(ln->link), ln->next ? ' ' : '\n'); } } /* 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_filtered(const CDT_state *cdt, bool (*filter_fn)(const CDT_state *, int, void *), void *filter_data, const char *lab) { LinkNode *ln; CDTVert *v, *vother; CDTEdge *e; CDTFace *f; SymEdge *se; int i, cnt; 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", vertname(v), PL(v), F2(v->co), PL(v->symedge)); if (v->merge_to_index == -1) { fprintf(stderr, "\n"); } else { fprintf(stderr, " merge to %s\n", vertname(cdt->vert_array[v->merge_to_index])); continue; } dump_id_list(v->input_ids, " "); se = v->symedge; cnt = 0; if (se) { fprintf(stderr, " edges out:\n"); do { if (se->next == NULL) { fprintf(stderr, " [NULL next/rot symedge, se=%x\n", PL(se)); break; } if (se->next->next == NULL) { fprintf(stderr, " [NULL next-next/rot symedge, se=%x\n", PL(se)); break; } vother = sym(se)->vert; fprintf(stderr, " %s (e=%x, se=%x)\n", vertname(vother), PL(se->edge), PL(se)); se = se->rot; cnt++; } while (se != v->symedge && cnt < 25); fprintf(stderr, "\n"); } } if (filter_fn) { return; } fprintf(stderr, "\nEDGES\n"); for (ln = cdt->edges; ln; ln = ln->next) { e = (CDTEdge *)ln->link; if (e->symedges[0].next == NULL) { continue; } fprintf(stderr, "%x:\n", PL(e)); for (i = 0; i < 2; i++) { se = &e->symedges[i]; fprintf(stderr, " se[%d] @%x: next=%x, rot=%x, vert=%x [%s] (%.2f,%.2f), edge=%x, face=%x\n", i, PL(se), PL(se->next), PL(se->rot), PL(se->vert), vertname(se->vert), F2(se->vert->co), PL(se->edge), PL(se->face)); } dump_id_list(e->input_ids, " "); } fprintf(stderr, "\nFACES\n"); for (ln = cdt->faces; ln; ln = ln->next) { f = (CDTFace *)ln->link; if (f->deleted) { continue; } if (f == cdt->outer_face) { fprintf(stderr, "%x: outer", PL(f)); } fprintf(stderr, " symedge=%x\n", PL(f->symedge)); dump_id_list(f->input_ids, " "); } fprintf(stderr, "\nOTHER\n"); fprintf(stderr, "outer_face=%x\n", PL(cdt->outer_face)); fprintf( stderr, "minx=%f, maxx=%f, miny=%f, maxy=%f\n", cdt->minx, cdt->maxx, cdt->miny, cdt->maxy); fprintf(stderr, "margin=%f\n", cdt->margin); } # 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 set_min_max(CDT_state *cdt) { int i; double minx, maxx, miny, maxy; double *co; minx = miny = DBL_MAX; maxx = maxy = -DBL_MAX; for (i = 0; i < cdt->vert_array_len; i++) { co = cdt->vert_array[i]->co; if (co[0] < minx) { minx = co[0]; } if (co[0] > maxx) { maxx = co[0]; } if (co[1] < miny) { miny = co[1]; } if (co[1] > maxy) { maxy = co[1]; } } if (minx != DBL_MAX) { cdt->minx = minx; cdt->miny = miny; cdt->maxx = maxx; cdt->maxy = maxy; } } 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. * Constraint edges are drawn thicker than non-constraint edges. * The first call creates DRAWFILE; subsequent calls append to it. */ # define DRAWFILE "/tmp/debug_draw.html" # 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 # define DRAW_EDGE_LABELS 0 static void cdt_draw_region( CDT_state *cdt, const char *lab, double minx, double miny, double maxx, double maxy) { static bool append = false; FILE *f = fopen(DRAWFILE, append ? "a" : "w"); int view_width, view_height; double width, height, aspect, scale; LinkNode *ln; CDTVert *v, *u; CDTEdge *e; int i, strokew; 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) { view_height = MAX_DRAW_HEIGHT; view_width = (int)(view_height / aspect); } scale = view_width / width; # define SX(x) ((x - minx) * scale) # define SY(y) ((maxy - y) * scale) if (!f) { printf("couldn't open file %s\n", DRAWFILE); return; } fprintf(f, "