diff options
-rw-r--r-- | source/blender/blenlib/intern/delaunay_2d.c | 3495 | ||||
-rw-r--r-- | tests/gtests/blenlib/BLI_delaunay_2d_test.cc | 896 |
2 files changed, 2947 insertions, 1444 deletions
diff --git a/source/blender/blenlib/intern/delaunay_2d.c b/source/blender/blenlib/intern/delaunay_2d.c index 4faaf1605e0..118949d1c46 100644 --- a/source/blender/blenlib/intern/delaunay_2d.c +++ b/source/blender/blenlib/intern/delaunay_2d.c @@ -17,8 +17,7 @@ /** \file * \ingroup bli * - * Dynamic Constrained Delaunay Triangulation. - * See paper by Marcelo Kallmann, Hanspeter Bieri, and Daniel Thalmann + * Constrained 2d Delaunay Triangulation. */ #include "MEM_guardedalloc.h" @@ -29,12 +28,11 @@ #include "BLI_math.h" #include "BLI_memarena.h" #include "BLI_mempool.h" -#include "BLI_rand.h" #include "BLI_delaunay_2d.h" /* Uncomment this define to get helpful debugging functions etc. defined. */ -// #define DEBUG_CDT +#define DEBUG_CDT struct CDTEdge; struct CDTFace; @@ -53,20 +51,21 @@ 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. */ + 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 { - double centroid[2]; /* Average of vertex coords. */ 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 { @@ -76,6 +75,7 @@ typedef struct CDT_state { CDTVert **vert_array; int vert_array_len; int vert_array_len_alloc; + int input_vert_tot; double minx; double miny; double maxx; @@ -83,19 +83,13 @@ typedef struct CDT_state { double margin; int visit_count; int face_edge_offset; - RNG *rng; MemArena *arena; BLI_mempool *listpool; double epsilon; + double epsilon_squared; bool output_prepared; } CDT_state; -typedef struct LocateResult { - enum { OnVert, OnEdge, InFace } loc_kind; - SymEdge *se; - double edge_lambda; -} LocateResult; - #define DLNY_ARENASIZE 1 << 14 /** @@ -105,60 +99,57 @@ typedef struct LocateResult { #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] -static void dump_se(const SymEdge *se, const char *lab); -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); +# 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 -/** - * Return 1 if a,b,c forms CCW angle, -1 if a CW angle, 0 if straight. - * For straight test, allow b to be withing eps of line. - */ -static int CCW_test(const double a[2], const double b[2], const double c[2], const double eps) +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) { - double det; - double ab; + return se->next->rot; +} - /* This is twice the signed area of triangle abc. */ - det = (b[0] - a[0]) * (c[1] - a[1]) - (c[0] - a[0]) * (b[1] - a[1]); - if (eps == 0.0) { - if (det > 0) { - return 1; - } - else if (det < 0) { - return -1; - } - else { - return 0; - } - } - ab = len_v2v2_db(a, b); - if (ab <= eps) { - return 0; - } - det /= ab; - if (det > eps) { - return 1; - } - else if (det < -eps) { - return -1; - } - return 0; +/** 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 - * CCW_test. */ -static bool in_line(const double a[2], const double b[2], const double c[2], double eps) +/** 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]) { - return fabs(len_v2v2_db(a, c) - (len_v2v2_db(a, b) + len_v2v2_db(b, c))) <= eps; + 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 @@ -176,21 +167,6 @@ static bool reachable(SymEdge *s1, SymEdge *s2, int limit) } #endif -static void calc_face_centroid(SymEdge *se) -{ - SymEdge *senext; - double *centroidp = se->face->centroid; - int count; - copy_v2_v2_db(centroidp, se->vert->co); - count = 1; - for (senext = se->next; senext != se; senext = senext->next) { - add_v2_v2_db(centroidp, senext->vert->co); - count++; - } - centroidp[0] /= count; - centroidp[1] /= count; -} - /** 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) { @@ -208,7 +184,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; + v->merge_to_index = -1; cdt->vert_array[cdt->vert_array_len++] = v; return v; } @@ -220,6 +196,7 @@ static CDTEdge *add_cdtedge( 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; @@ -243,6 +220,7 @@ static CDTFace *add_cdtface(CDT_state *cdt) 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; } @@ -290,37 +268,24 @@ static void add_list_to_input_ids(LinkNode **dst, const LinkNode *src, CDT_state } } -/** Return other #SymEdge for same #CDTEdge as se. */ -static inline SymEdge *sym(const SymEdge *se) -{ - return se->next->rot; -} - -/** Return SymEdge whose next is se. */ -static inline SymEdge *prev(const SymEdge *se) -{ - return se->rot->next->rot; -} - -static inline bool is_border_edge(const CDTEdge *e, const CDT_state *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; } -/** Does one edge of this edge touch the frame? */ -static bool edge_touches_frame(const CDTEdge *e) +BLI_INLINE bool is_constrained_edge(const CDTEdge *e) { - return e->symedges[0].vert->index < 4 || e->symedges[1].vert->index < 4; + return e->input_ids != NULL; } -static inline bool is_constrained_edge(const CDTEdge *e) +BLI_INLINE bool is_deleted_edge(const CDTEdge *e) { - return e->input_ids != NULL; + return e->symedges[0].next == NULL; } -static inline bool is_deleted_edge(const CDTEdge *e) +BLI_INLINE bool is_original_vert(const CDTVert *v, CDT_state *cdt) { - return e->symedges[0].next == NULL; + return (v->index < cdt->input_vert_tot); } /** Is there already an edge between a and b? */ @@ -357,7 +322,6 @@ static bool vert_touches_face(const CDTVert *v, const CDTFace *f) * 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. - * The centroids of both faces are recalculated. * Return the new diagonal's CDTEdge *. */ static CDTEdge *add_diagonal(CDT_state *cdt, SymEdge *s1, SymEdge *s2) @@ -366,8 +330,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, 2000)); - BLI_assert(reachable(s2, s1, 2000)); + BLI_assert(reachable(s1, s2, 20000)); + BLI_assert(reachable(s2, s1, 20000)); fold = s1->face; fnew = add_cdtface(cdt); s1prev = prev(s1); @@ -392,12 +356,61 @@ static CDTEdge *add_diagonal(CDT_state *cdt, SymEdge *s1, SymEdge *s2) se->face = fnew; } add_list_to_input_ids(&fnew->input_ids, fold->input_ids, cdt); - calc_face_centroid(sdiag); - calc_face_centroid(sdiagsym); 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. @@ -435,14 +448,12 @@ static CDTEdge *split_edge(CDT_state *cdt, SymEdge *se, double lambda) newsesym->vert->symedge = newsesym; } add_list_to_input_ids(&e->input_ids, se->edge->input_ids, cdt); - calc_face_centroid(se); - calc_face_centroid(sesym); 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; the centroid is updated. + * 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. @@ -519,820 +530,486 @@ static void delete_edge(CDT_state *cdt, SymEdge *e) cdt->outer_face = aface; } } - if (aface != cdt->outer_face) { - calc_face_centroid(f); - } } -/** - * The initial structure will be the rectangle with opposite corners (minx,miny) - * and (maxx,maxy), and a diagonal going between those two corners. - * We keep track of the outer face (surrounding the entire structure; its boundary - * is the clockwise traversal of the bounding box rectangle initially) in cdt->outer_face. - * - * The vertices are kept as pointers in an array (which may need to be reallocated from - * time to time); the edges and faces are kept in lists. Sometimes edges and faces are deleted, - * marked by setting all pointers to NULL (for edges), or setting the deleted flag to true (for - * faces). - * - * A #MemArena is allocated to do all allocations from except for link list nodes; a listpool - * is created for link list node allocations. - * - * The epsilon argument is stored and used in "near enough" distance calculations. - * - * When done, caller must call BLI_constrained_delaunay_free to free - * the memory used by the returned #CDT_state. - */ -static CDT_state *cdt_init(double minx, double maxx, double miny, double maxy, double epsilon) +static CDT_state *new_cdt_init(const CDT_input *in) { - double x0, x1, y0, y1; - double margin; - CDTVert *v[4]; - CDTEdge *e[4]; - CDTFace *f0, *fouter; - int i, inext, iprev; + int i; MemArena *arena = BLI_memarena_new(DLNY_ARENASIZE, __func__); - CDT_state *cdt = BLI_memarena_alloc(arena, sizeof(CDT_state)); - cdt->edges = NULL; - cdt->faces = NULL; - cdt->vert_array_len = 0; - cdt->vert_array_len_alloc = 32; + 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->minx = minx; - cdt->miny = miny; - cdt->maxx = maxx; - cdt->maxy = maxy; - cdt->arena = arena; - cdt->listpool = BLI_mempool_create(sizeof(LinkNode), 128, 128, 0); - cdt->rng = BLI_rng_new(0); - cdt->epsilon = epsilon; - - /* Expand bounding box a bit and make initial CDT from it. */ - margin = DLNY_MARGIN_PCT * max_dd(maxx - minx, maxy - miny) / 100.0; - if (margin <= 0.0) { - margin = 1.0; - } - if (margin < epsilon) { - margin = 4 * epsilon; /* Make sure constraint verts don't merge with border verts. */ - } - cdt->margin = margin; - x0 = minx - margin; - y0 = miny - margin; - x1 = maxx + margin; - y1 = maxy + margin; - - /* Make a quad, then split it with a diagonal. */ - v[0] = add_cdtvert(cdt, x0, y0); - v[1] = add_cdtvert(cdt, x1, y0); - v[2] = add_cdtvert(cdt, x1, y1); - v[3] = add_cdtvert(cdt, x0, y1); - cdt->outer_face = fouter = add_cdtface(cdt); - f0 = add_cdtface(cdt); - for (i = 0; i < 4; i++) { - e[i] = add_cdtedge(cdt, v[i], v[(i + 1) % 4], f0, fouter); - } - for (i = 0; i < 4; i++) { - inext = (i + 1) % 4; - iprev = (i + 3) % 4; - e[i]->symedges[0].next = &e[inext]->symedges[0]; - e[inext]->symedges[1].next = &e[i]->symedges[1]; - e[i]->symedges[0].rot = &e[iprev]->symedges[1]; - e[iprev]->symedges[1].rot = &e[i]->symedges[0]; - } - calc_face_centroid(&e[0]->symedges[0]); - add_diagonal(cdt, &e[0]->symedges[0], &e[2]->symedges[0]); - fouter->centroid[0] = fouter->centroid[1] = 0.0; - - cdt->visit_count = 0; - cdt->output_prepared = false; - cdt->face_edge_offset = 0; + 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 cdt_free(CDT_state *cdt) +static void new_cdt_free(CDT_state *cdt) { - BLI_rng_free(cdt->rng); BLI_mempool_destroy(cdt->listpool); BLI_memarena_free(cdt->arena); } -static bool locate_point_final(const double p[2], - SymEdge *tri_se, - bool try_neighbors, - const double epsilon, - LocateResult *r_lr) -{ - /* 'p' should be in or on our just outside of 'cur_tri'. */ - double dist_inside[3]; - int i; - SymEdge *se; - const double *a, *b; - double lambda, close[2]; - bool done = false; -#ifdef DEBUG_CDT - int dbg_level = 0; - - if (dbg_level > 0) { - fprintf(stderr, "locate_point_final %d\n", try_neighbors); - dump_se(tri_se, "tri_se"); - fprintf(stderr, "\n"); - } -#endif - se = tri_se; - i = 0; - do { -#ifdef DEBUG_CDT - if (dbg_level > 1) { - fprintf(stderr, "%d: ", i); - dump_se(se, "search se"); - } -#endif - a = se->vert->co; - b = se->next->vert->co; - lambda = closest_to_line_v2_db(close, p, a, b); - double len_close_p = len_v2v2_db(close, p); - if (len_close_p < epsilon) { - if (len_v2v2_db(p, a) < epsilon) { -#ifdef DEBUG_CDT - if (dbg_level > 0) { - fprintf(stderr, "OnVert case a (%.2f,%.2f)\n", F2(a)); - } -#endif - r_lr->loc_kind = OnVert; - r_lr->se = se; - r_lr->edge_lambda = 0.0; - done = true; - } - else if (len_v2v2_db(p, b) < epsilon) { -#ifdef DEBUG_CDT - if (dbg_level > 0) { - fprintf(stderr, "OnVert case b (%.2f,%.2f)\n", F2(b)); - } -#endif - r_lr->loc_kind = OnVert; - r_lr->se = se->next; - r_lr->edge_lambda = 0.0; - done = true; - } - else if (lambda > 0.0 && lambda < 1.0) { -#ifdef DEBUG_CDT - if (dbg_level > 0) { - fprintf(stderr, "OnEdge case, lambda=%f\n", lambda); - dump_se(se, "se"); - } -#endif - r_lr->loc_kind = OnEdge; - r_lr->se = se; - r_lr->edge_lambda = lambda; - done = true; - } - } - else { - dist_inside[i] = len_close_p; - dist_inside[i] = CCW_test(a, b, p, epsilon) >= 0 ? len_close_p : -len_close_p; - } - i++; - se = se->next; - } while (se != tri_se && !done); - if (!done) { -#ifdef DEBUG_CDT - if (dbg_level > 1) { - fprintf(stderr, - "not done, dist_inside=%f %f %f\n", - dist_inside[0], - dist_inside[1], - dist_inside[2]); - } -#endif - if (dist_inside[0] >= 0.0 && dist_inside[1] >= 0.0 && dist_inside[2] >= 0.0) { -#ifdef DEBUG_CDT - if (dbg_level > 0) { - fprintf(stderr, "InFace case\n"); - dump_se_cycle(tri_se, "tri", 10); - } -#endif - r_lr->loc_kind = InFace; - r_lr->se = tri_se; - r_lr->edge_lambda = 0.0; - done = true; - } - else if (try_neighbors) { - for (se = tri_se->next; se != tri_se; se = se->next) { - if (locate_point_final(p, se, false, epsilon, r_lr)) { - done = true; - break; - } - } - if (!done) { - /* Shouldn't happen desperation mode: pick something. */ - se = NULL; - if (dist_inside[0] > 0) { - se = tri_se; - } - if (dist_inside[1] > 0 && (se == NULL || dist_inside[1] < dist_inside[i])) { - se = tri_se->next; - } - if (se == NULL) { - se = tri_se->next->next; - } - a = se->vert->co; - b = se->next->vert->co; - lambda = closest_to_line_v2_db(close, p, a, b); - if (lambda <= 0.0) { - r_lr->loc_kind = OnVert; - r_lr->se = se; - r_lr->edge_lambda = 0.0; - } - else if (lambda >= 1.0) { - r_lr->loc_kind = OnVert; - r_lr->se = se->next; - r_lr->edge_lambda = 0.0; - } - else { - r_lr->loc_kind = OnEdge; - r_lr->se = se->next; - r_lr->edge_lambda = lambda; - } -#ifdef DEBUG_CDT - if (dbg_level > 0) { - fprintf( - stderr, "desperation case kind=%u lambda=%f\n", r_lr->loc_kind, r_lr->edge_lambda); - dump_se(r_lr->se, "se"); - BLI_assert(0); /* While developing, catch these "should not happens" */ - } -#endif - fprintf(stderr, "desperation! point=(%g,%g)\n", p[0], p[1]); // TODO: remove - return true; - } - } - } - return done; -} +typedef struct SiteInfo { + CDTVert *v; + int orig_index; +} SiteInfo; -static LocateResult locate_point(CDT_state *cdt, const double p[2]) +static int site_lexicographic_cmp(const void *a, const void *b) { - LocateResult lr; - SymEdge *cur_se, *next_se, *next_se_sym; - CDTFace *cur_tri; - bool done; - int sample_n, i, k; - CDTVert *v, *best_start_vert; - double dist_squared, best_dist_squared; - double *a, *b, *c; - const double epsilon = cdt->epsilon; - int visit = ++cdt->visit_count; - int loop_count = 0; -#ifdef DEBUG_CDT - int dbg_level = 0; + const SiteInfo *s1 = a; + const SiteInfo *s2 = b; + const double *co1 = s1->v->co; + const double *co2 = s2->v->co; - if (dbg_level > 0) { - fprintf(stderr, "locate_point (%.2f,%.2f), visit_index=%d\n", F2(p), visit); - } -#endif - /* Starting point determined by closest to p in an n ** (1/3) sized sample of current points. */ - BLI_assert(cdt->vert_array_len > 0); - sample_n = (int)round(pow((double)cdt->vert_array_len, 0.33333)); - if (sample_n < 1) { - sample_n = 1; + if (co1[0] < co2[0]) { + return -1; } - best_start_vert = NULL; - best_dist_squared = DBL_MAX; - for (k = 0; k < sample_n; k++) { - /* Yes, this may try some i's more than once, - * but will still get about an n ** (1/3) size sample. */ - i = (int)(BLI_rng_get_uint(cdt->rng) % cdt->vert_array_len); - v = cdt->vert_array[i]; - dist_squared = len_squared_v2v2_db(p, v->co); -#ifdef DEBUG_CDT - if (dbg_level > 0) { - fprintf(stderr, "try start vert %d, dist_squared=%f\n", i, dist_squared); - dump_v(v, "v"); - } -#endif - if (dist_squared < best_dist_squared) { - best_dist_squared = dist_squared; - best_start_vert = v; - } + else if (co1[0] > co2[0]) { + return 1; } - cur_se = &best_start_vert->symedge[0]; - if (cur_se->face == cdt->outer_face) { - cur_se = cur_se->rot; - BLI_assert(cur_se->face != cdt->outer_face); + else if (co1[1] < co2[1]) { + return -1; } -#ifdef DEBUG_CDT - if (dbg_level > 0) { - dump_se(cur_se, "start vert edge"); + else if (co1[1] > co2[1]) { + return 1; } -#endif - done = false; - while (!done) { - /* Find edge of cur_tri that separates p and t's centroid, - * and where other tri over the edge is unvisited. */ -#ifdef DEBUG_CDT - if (dbg_level > 0) { - dump_se_cycle(cur_se, "cur search face", 5); - } -#endif - cur_tri = cur_se->face; - BLI_assert(cur_tri != cdt->outer_face); - cur_tri->visit_index = visit; - /* Is p in or on current triangle? */ - a = cur_se->vert->co; - b = cur_se->next->vert->co; - c = cur_se->next->next->vert->co; - if (CCW_test(a, b, p, 0.0) >= 0 && CCW_test(b, c, p, 0.0) >= 0 && - CCW_test(c, a, p, 0.0) >= 0) { -#ifdef DEBUG_CDT - if (dbg_level > 1) { - fprintf(stderr, "p in current triangle\n"); - } -#endif - done = locate_point_final(p, cur_se, false, epsilon, &lr); - BLI_assert(done == true); - break; - } - bool found_next = false; - next_se = cur_se; - do { - a = next_se->vert->co; - b = next_se->next->vert->co; - c = next_se->next->next->vert->co; -#ifdef DEBUG_CDT - if (dbg_level > 1) { - dump_se(next_se, "search edge"); - fprintf(stderr, "tri centroid=(%.3f,%.3f)\n", F2(cur_tri->centroid)); - validate_face_centroid(next_se); - } -#endif - next_se_sym = sym(next_se); - if (CCW_test(a, b, p, 0.0) <= 0 && next_se->face != cdt->outer_face) { -#ifdef DEBUG_CDT - if (dbg_level > 1) { - fprintf(stderr, "CCW_test(a, b, p) <= 0\n"); - } -#endif -#ifdef DEBUG_CDT - if (dbg_level > 0) { - dump_se(next_se_sym, "next_se_sym"); - fprintf(stderr, "next_se_sym face visit=%d\n", next_se_sym->face->visit_index); - } -#endif - if (next_se_sym->face->visit_index != visit) { -#ifdef DEBUG_CDT - if (dbg_level > 0) { - fprintf(stderr, "found edge to cross\n"); - } -#endif - found_next = true; - cur_se = next_se_sym; - break; - } - } - next_se = next_se->next; - } while (next_se != cur_se); - if (!found_next) { - done = locate_point_final(p, cur_se, true, epsilon, &lr); - BLI_assert(done = true); - done = true; - } - if (++loop_count > 1000000) { - fprintf(stderr, "infinite search loop?\n"); - done = locate_point_final(p, cur_se, true, epsilon, &lr); - } + else if (s1->orig_index < s2->orig_index) { + return -1; } - - return lr; -} - -/** - * Return true if circumcircle(v1, v2, v3) does not contain p. - * To avoid possible infinite flip loops, we will say true even if p is inside the circle - * but less than epsilon from the boundary; or if v1, v2, v3, form a straight line. - */ -static bool delaunay_check(CDTVert *v1, CDTVert *v2, CDTVert *v3, CDTVert *p, const double epsilon) -{ - double x1, y1, x2, y2, den, cenx, ceny, rad, pc, a, b, w, z, q, s; - - /* Find center and radius of circum-circle of v1,v2,v3. - * Transform coords so v3 is at origin to help reduce floating point error - * when coordinates are far from (0,0) but close together. - */ - x1 = v1->co[0] - v3->co[0]; - y1 = v1->co[1] - v3->co[1]; - x2 = v2->co[0] - v3->co[0]; - y2 = v2->co[1] - v3->co[1]; - den = 2.0 * (x1 * y2 - x2 * y1); - if (UNLIKELY(den == 0.0)) { - /* v1, v2, v3 are in a line. */ - return true; + else if (s1->orig_index > s2->orig_index) { + return 1; } - /* cen[0] = det(x1**2 + y1**2, y1, x2**2 + y2**2, y2) / den - * cen[1] = det(x1, x1**2 + y1**2, x2, x2**2 + y2**2) / den - * den = 2 * det(x1, y1, x2, y2) - */ - a = x1 * x1 + y1 * y1; - b = x2 * x2 + y2 * y2; - cenx = (a * y2 - b * y1) / den; - ceny = (x1 * b - x2 * a) / den; - w = x1 - cenx; - z = y1 - ceny; - rad = sqrt(w * w + z * z); - q = p->co[0] - v3->co[0] - cenx; - s = p->co[1] - v3->co[1] - ceny; - pc = sqrt(q * q + s * s); - return (pc >= rad - epsilon); -} - -/* Return true if we can flip edge v1-v3 to edge v2-v4 inside quad v1v2v3v4 (in CCW order). - * We can do this if angles v4-v1-v2 and v2-v3-v4 are both CCW or straight. - */ -static inline bool can_flip(CDTVert *v1, CDTVert *v2, CDTVert *v3, CDTVert *v4) -{ - return CCW_test(v4->co, v1->co, v2->co, 0.0) >= 0 && CCW_test(v2->co, v3->co, v4->co, 0.0) >= 0; + return 0; } -/** Use LinkNode linked list as stack of SymEdges, allocating from cdt->listpool. */ -typedef LinkNode *Stack; - -static inline void push(Stack *stack, SymEdge *se, CDT_state *cdt) +BLI_INLINE bool vert_left_of_symedge(CDTVert *v, SymEdge *se) { - BLI_linklist_prepend_pool(stack, se, cdt->listpool); + return orient2d(v->co, se->vert->co, se->next->vert->co) > 0.0; } -static inline SymEdge *pop(Stack *stack, CDT_state *cdt) +BLI_INLINE bool vert_right_of_symedge(CDTVert *v, SymEdge *se) { - return (SymEdge *)BLI_linklist_pop_pool(stack, cdt->listpool); + return orient2d(v->co, se->next->vert->co, se->vert->co) > 0.0; } -static inline bool is_empty(Stack *stack) +/* Is se above basel? */ +BLI_INLINE bool dc_tri_valid(SymEdge *se, SymEdge *basel, SymEdge *basel_sym) { - return *stack == NULL; + return orient2d(se->next->vert->co, basel_sym->vert->co, basel->vert->co) > 0.0; } -/** - * <pre> - * /\ /\ - * /a|\ / \ - * / | sesym / \ - * / | \ / \ - * . b | d . -> . se______ - * \ se| / \ / - * \ |c/ \ / - * \ |/ \ / - * </pre> +/* 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 flip(SymEdge *se, CDT_state *cdt) +static void dc_tri( + CDT_state *cdt, SiteInfo *sites, int start, int end, SymEdge **r_le, SymEdge **r_re) { - SymEdge *a, *b, *c, *d; - SymEdge *sesym, *asym, *bsym, *csym, *dsym; - CDTFace *t1, *t2; - CDTVert *v1, *v2; + 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 - const int dbg_level = 0; -#endif + char label_buf[100]; + int dbg_level = 0; - sesym = sym(se); -#ifdef DEBUG_CDT if (dbg_level > 0) { - fprintf(stderr, "flip\n"); - dump_se(se, "se"); - dump_se(sesym, "sesym"); - } -#endif - a = se->next; - b = a->next; - c = sesym->next; - d = c->next; - asym = sym(a); - bsym = sym(b); - csym = sym(c); - dsym = sym(d); -#ifdef DEBUG_CDT - if (dbg_level > 1) { - dump_se(a, "a"); - dump_se(b, "b"); - dump_se(c, "c"); - dump_se(d, "d"); + fprintf(stderr, "DC_TRI start=%d end=%d\n", start, end); } #endif - v1 = se->vert; - v2 = sesym->vert; - t1 = a->face; - t2 = c->face; - - se->vert = b->vert; - sesym->vert = d->vert; - - a->next = se; - se->next = d; - d->next = a; - - sesym->next = b; - b->next = c; - c->next = sesym; - a->rot = dsym; - b->rot = se; - se->rot = asym; - - c->rot = bsym; - d->rot = sesym; - sesym->rot = csym; - - a->face = se->face = d->face = t1; - sesym->face = b->face = c->face = t2; - - if (v1->symedge == se) { - v1->symedge = c; + BLI_assert(r_le != NULL && r_re != NULL); + if (n <= 1) { + *r_le = NULL; + *r_re = NULL; + return; } - if (v2->symedge == sesym) { - v2->symedge = a; + 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); - calc_face_centroid(a); - calc_face_centroid(sesym); - + /* 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, "after flip\n"); - dump_se_cycle(a, "a cycle", 5); - dump_se_cycle(sesym, "sesym cycle", 5); + 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 - if (cdt) { - /* Pass. */ - } -} -static void flip_edges(CDTVert *v, Stack *stack, CDT_state *cdt) -{ - SymEdge *se, *sesym; - CDTVert *a, *b, *c, *d; - SymEdge *tri_without_p; - bool is_delaunay; - const double epsilon = cdt->epsilon; - int count = 3; + /* 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 - const int dbg_level = 0; if (dbg_level > 0) { - fprintf(stderr, "flip_edges, v=(%.2f,%.2f)\n", F2(v->co)); + fprintf(stderr, "common lower tangent is between\n"); + dump_se(rdi, "rdi"); + dump_se(ldi, "ldi"); } #endif - while (!is_empty(stack)) { - if (++count > 10000) { - fprintf(stderr, "infinite flip loop?\n"); - return; - } - se = pop(stack, cdt); + ebasel = connect_separate_parts(cdt, sym(rdi)->next, ldi); + basel = &ebasel->symedges[0]; + basel_sym = &ebasel->symedges[1]; #ifdef DEBUG_CDT - if (dbg_level > 0) { - dump_se(se, "flip_edges popped"); - } + if (dbg_level > 1) { + dump_se(basel, "basel"); + cdt_draw(cdt, "after basel made"); + } #endif - if (!is_constrained_edge(se->edge)) { - /* Edge is not constrained; is it Delaunay? */ + 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) { - dump_se_cycle(se, "unconstrained edge", 5); - } - else if (dbg_level > 0) { - fprintf(stderr, "unconstrained edge\n"); - } + if (dbg_level > 1) { + fprintf(stderr, "\ntop of merge loop\n"); + dump_se(lcand, "lcand"); + dump_se(rcand, "rcand"); + dump_se(basel, "basel"); + } #endif - a = se->vert; - b = se->next->vert; - c = se->next->next->vert; - sesym = sym(se); - d = sesym->next->next->vert; + if (dc_tri_valid(lcand, basel, basel_sym)) { #ifdef DEBUG_CDT if (dbg_level > 1) { - fprintf(stderr, "a=(%.3f,%.3f) b=(%.3f,%.3f)\n", F2(a->co), F2(b->co)); - fprintf(stderr, "c=(%.3f,%.3f) d=(%.3f,%.3f)\n", F2(c->co), F2(d->co)); + fprintf(stderr, "found valid lcand\n"); + dump_se(lcand, " lcand"); } #endif - if (v == c) { - tri_without_p = sesym; - is_delaunay = delaunay_check(a, b, c, d, epsilon); + 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, "v==c, delaunay(a,b,c,d)=%d\n", is_delaunay); + fprintf(stderr, "incircle says to remove lcand\n"); + dump_se(lcand, " lcand"); } #endif + t = lcand->rot; + delete_edge(cdt, sym(lcand)); + lcand = t; } - else { - tri_without_p = se; - BLI_assert(d == v); - is_delaunay = delaunay_check(b, a, d, c, epsilon); + } + /* 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, "v!=c, delaunay(b,a,d,c)=%d\n", is_delaunay); - } -#endif + if (dbg_level > 1) { + fprintf(stderr, "found valid rcand\n"); + dump_se(rcand, " rcand"); } - if (!is_delaunay && can_flip(a, d, b, c)) { - /* Push two edges of tri without p that aren't se. */ +#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, "maybe pushing more edges\n"); + fprintf(stderr, "incircle says to remove rcand\n"); + dump_se(lcand, " rcand"); } #endif - if (!is_border_edge(tri_without_p->next->edge, cdt)) { + 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) { - dump_se(tri_without_p->next, "push1"); - } + 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 - push(stack, tri_without_p->next, cdt); - } - if (!is_border_edge(tri_without_p->next->next->edge, cdt)) { + 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) { - dump_se(tri_without_p->next->next, "\npush2"); - } + if (dbg_level > 0) { + fprintf(stderr, "connecting rcand\n"); + dump_se(basel_sym, " se1=basel_sym"); + dump_se(rcand->next, " se2=rcand->next"); + } #endif - push(stack, tri_without_p->next->next, cdt); - } - flip(se, cdt); + ebasel = add_diagonal(cdt, rcand->next, basel_sym); + } + else { #ifdef DEBUG_CDT - if (dbg_level > 2) { - dump_cdt(cdt, "after flip"); - cdt_draw(cdt, "afer flip"); - validate_cdt(cdt, true); - } -#endif + 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); } -/** - * Splits e at lambda and returns a #SymEdge with new vert as its vert. - * The two opposite triangle vertices to e are connect to new point. - * <pre> - * /\ /\ - * /f|\ / |\ - * / |j\ / | \ - * / | i\ / k| \ - * . | . -> . l_ p m_. - * \g | / \ | / - * \ |h/ \ | / - * \e|/ \ e|/ - * - * t1 = {e, f, g}; t2 = {h, i, j}; - * t1' = {e, l.sym, g}; t2' = {h, m.sym, e'.sym} - * t3 = {k, f, l}; t4 = {m, i, j} - * </pre> - */ -static CDTVert *insert_point_in_edge(CDT_state *cdt, SymEdge *e, double lambda) +/* Guibas-Stolfi Divide-and_Conquer algorithm. */ +static void dc_triangulate(CDT_state *cdt, SiteInfo *sites, int nsites) { - SymEdge *f, *g, *h, *i, *j, *k; - CDTEdge *ke; - CDTVert *p; - Stack stack; - /* Split e at lambda. */ - - f = e->next; - g = f->next; - BLI_assert(g->next == e); - j = sym(e); - h = j->next; - i = h->next; - BLI_assert(i->next == j); - - ke = split_edge(cdt, e, lambda); - k = &ke->symedges[0]; - p = k->vert; + int i, j, n; + SymEdge *le, *re; - add_diagonal(cdt, g, k); - add_diagonal(cdt, sym(e), i); - - stack = NULL; - if (!is_border_edge(f->edge, cdt)) { - push(&stack, f, cdt); - } - if (!is_border_edge(g->edge, cdt)) { - push(&stack, g, cdt); - } - if (!is_border_edge(h->edge, cdt)) { - push(&stack, h, cdt); + /* 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++; + } } - if (!is_border_edge(i->edge, cdt)) { - push(&stack, i, cdt); + n = i; + if (n == 0) { + return; } - flip_edges(k->vert, &stack, cdt); - return p; + dc_tri(cdt, sites, 0, n, &le, &re); } /** - * Inserts p inside e's triangle and connects the three cornders - * of the triangle to the new point. Returns a SymEdge that has - * new point as its point. - * <pre> - * * * - * *g * * .j* - * * * * . * - * * p * -> * 1. p . 3* - * * * * . . * - * * e f* * . h 2 i . * - * * * * * * * * * * * * * * * * * * * * * * * * * * * - * </pre> + * 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 CDTVert *insert_point_in_face(CDT_state *cdt, SymEdge *e, const double p[2]) +static void initial_triangulation(CDT_state *cdt) { - SymEdge *f, *g, *h, *i, *j; - SymEdge *esym, *fsym, *gsym, *hsym, *isym, *jsym; - CDTVert *v; - CDTEdge *he, *ie, *je; - CDTFace *t1, *t2, *t3; - Stack stack; + 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, "insert point in face, p=(%.3f,%.3f)\n", F2(p)); - dump_se_cycle(e, "insert face", 20); + fprintf(stderr, "\nINITIAL TRIANGULATION\n\n"); } #endif - f = e->next; - g = f->next; - esym = sym(e); - fsym = sym(f); - gsym = sym(g); - t1 = e->face; - t2 = add_cdtface(cdt); - t3 = add_cdtface(cdt); - - v = add_cdtvert(cdt, p[0], p[1]); - he = add_cdtedge(cdt, e->vert, v, t1, t2); - h = &he->symedges[0]; - hsym = &he->symedges[1]; - ie = add_cdtedge(cdt, f->vert, v, t2, t3); - i = &ie->symedges[0]; - isym = &ie->symedges[1]; - je = add_cdtedge(cdt, g->vert, v, t3, t1); - j = &je->symedges[0]; - jsym = &je->symedges[1]; - - e->next = i; - i->next = hsym; - hsym->next = e; - e->face = t2; - - f->next = j; - j->next = isym; - isym->next = f; - f->face = t3; - - g->next = h; - h->next = jsym; - jsym->next = g; - g->face = t1; - - e->rot = h; - i->rot = esym; - hsym->rot = isym; - - f->rot = i; - j->rot = fsym; - isym->rot = jsym; - - g->rot = j; - h->rot = gsym; - jsym->rot = hsym; - - calc_face_centroid(e); - calc_face_centroid(f); - calc_face_centroid(g); - + /* 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 > 1) { - fprintf(stderr, "after initial insert:\n"); - dump_se_cycle(e, "e", 20); - dump_se_cycle(f, "f", 20); - dump_se_cycle(g, "g", 20); - if (dbg_level > 2) { - dump_cdt(cdt, "after initial insert, before flip"); - cdt_draw(cdt, "after initial insert, before flip"); - validate_cdt(cdt, true); + 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 - stack = NULL; - if (!is_border_edge(e->edge, cdt)) { - push(&stack, e, cdt); - } - if (!is_border_edge(f->edge, cdt)) { - push(&stack, f, cdt); - } - if (!is_border_edge(g->edge, cdt)) { - push(&stack, g, cdt); + /* 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++; + } } - flip_edges(v, &stack, cdt); - return v; + /* 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; } /** @@ -1347,12 +1024,16 @@ static void re_delaunay_triangulate(CDT_state *cdt, SymEdge *se) SymEdge *ss, *first, *cse; CDTVert *a, *b, *c, *v; CDTEdge *ebc, *eca; - const double epsilon = cdt->epsilon; 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); @@ -1391,7 +1072,7 @@ static void re_delaunay_triangulate(CDT_state *cdt, SymEdge *se) #endif for (ss = first->next; ss != se; ss = ss->next) { v = ss->vert; - if (!delaunay_check(a, b, c, v, epsilon)) { + if (incircle(a->co, b->co, c->co, v->co) > 0.0) { c = v; cse = ss; #ifdef DEBUG_CDT @@ -1440,50 +1121,6 @@ static void re_delaunay_triangulate(CDT_state *cdt, SymEdge *se) } /** - * Add a constrained point to cdt structure, and return the corresponding CDTVert*. - * May not be at exact coords given, because it can be merged with an existing vertex - * or moved to an existing edge (which could be a triangulation edge, not just a constraint one) - * if the point is within cdt->epsilon of those other elements. - * - * input_id will be added to the list of input_ids for the returned CDTVert (don't use -1 for id). - * - * Assumes cdt has been initialized, with min/max bounds that contain coords. - * Assumes that #BLI_constrained_delaunay_get_output has not been called yet. - */ -static CDTVert *add_point_constraint(CDT_state *cdt, const double coords[2], int input_id) -{ - LocateResult lr; - CDTVert *v; -#ifdef DEBUG_CDT - const int dbg_level = 0; -#endif - - BLI_assert(!cdt->output_prepared); -#ifdef DEBUG_CDT - if (dbg_level > 0) { - fprintf(stderr, "add point constraint (%.3f,%.3f), id=%d\n", F2(coords), input_id); - } -#endif - lr = locate_point(cdt, coords); -#ifdef DEBUG_CDT - if (dbg_level > 0) { - fprintf(stderr, " locate result has loc_kind %u\n", lr.loc_kind); - } -#endif - if (lr.loc_kind == OnVert) { - v = lr.se->vert; - } - else if (lr.loc_kind == OnEdge) { - v = insert_point_in_edge(cdt, lr.se, lr.edge_lambda); - } - else { - v = insert_point_in_face(cdt, lr.se, coords); - } - add_to_input_ids(&v->input_ids, input_id, cdt); - return v; -} - -/** * 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. @@ -1505,11 +1142,11 @@ static void add_edge_constraint( SymEdge *t, *tstart, *tout, *tnext; SymEdge *se; CDTEdge *edge; - int ccw1, ccw2, isect; - int i, search_count, visit; + int isect; + double orient1, orient2; + int i, search_count; double curco[2]; double lambda; - const double epsilon = cdt->epsilon; bool done, state_through_vert; LinkNodePair edge_list = {NULL, NULL}; typedef struct CrossData { @@ -1523,7 +1160,7 @@ static void add_edge_constraint( CrossData *cd; BLI_array_staticdeclare(crossings, 128); #ifdef DEBUG_CDT - const int dbg_level = 0; + int dbg_level = 0; #endif /* Find path through structure from v1 to v2 and record how we got there in crossings. @@ -1564,6 +1201,10 @@ static void add_edge_constraint( #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"); @@ -1571,9 +1212,6 @@ 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) { @@ -1600,7 +1238,7 @@ static void add_edge_constraint( cdata.vert = v2; BLI_array_append(crossings, cdata); #ifdef DEBUG_CDT - if (dbg_level > 1) { + if (dbg_level > 0) { fprintf(stderr, "special one segment case\n"); dump_se(t, " "); } @@ -1611,10 +1249,6 @@ static void add_edge_constraint( 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; @@ -1636,9 +1270,6 @@ static void add_edge_constraint( dump_se_cycle(t, "current t ", 4); } #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; @@ -1648,7 +1279,7 @@ static void add_edge_constraint( BLI_array_append(crossings, cdata); if (t->vert == v2) { #ifdef DEBUG_CDT - if (dbg_level > 0) { + if (dbg_level > 1) { fprintf(stderr, "found v2, so done\n"); } #endif @@ -1661,20 +1292,18 @@ static void add_edge_constraint( 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); + 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, "ccw1=%d, ccw2=%d\n", ccw1, ccw2); + fprintf(stderr, "orient1=%g\n", orient1); } #endif - if (ccw1 == 0 && in_line(t->vert->co, va->co, v2->co, epsilon) && - va->visit_index != visit) { + if (orient1 == 0.0 && in_line(t->vert->co, va->co, v2->co)) { #ifdef DEBUG_CDT - if (dbg_level > 0) { + if (dbg_level > 1) { fprintf(stderr, "ray goes through va\n"); } #endif @@ -1683,30 +1312,35 @@ static void add_edge_constraint( t = t->next; break; } - else if (ccw2 == 0 && in_line(t->vert->co, vb->co, v2->co, epsilon) && - vb->visit_index != visit) { + else if (t->face != cdt->outer_face) { + orient2 = orient2d(t->vert->co, vb->co, v2->co); #ifdef DEBUG_CDT - if (dbg_level > 0) { - fprintf(stderr, "ray goes through vb\n"); + if (dbg_level > 1) { + fprintf(stderr, "orient2=%g\n", orient2); } #endif - 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)) { + if (orient2 == 0.0 && in_line(t->vert->co, vb->co, v2->co)) { #ifdef DEBUG_CDT - if (dbg_level > 0) { - fprintf(stderr, "segment intersects\n"); + 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; + state_through_vert = false; + tout = t; + t = t->next; + break; + } } t = t->rot; #ifdef DEBUG_CDT @@ -1715,57 +1349,13 @@ static void add_edge_constraint( } #endif } 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, "add_edge_constraint desperation search\n"); - } -#endif - 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 > 0) { - fprintf(stderr, "add_edge_constraint desperation search chose to go through\n"); - dump_v(bestt->vert, "desperation vert"); - } -#endif - tout = bestt; - t = t->next; - } - } + 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. - * Whatever we make t next should not have both vert and face visited. */ va = t->vert; vb = t->next->vert; @@ -1784,19 +1374,41 @@ static void add_edge_constraint( } #endif 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. */ + 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, "add_edge_constraint no intersect found, using lambda = 0.5\n"); + fprintf(stderr, "intersect is collinear, treating as through vert\n"); } + dump_v(va, "va"); #endif - lambda = 0.5; + 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 %f along va--vb\n", lambda); - if (dbg_level == 1) { + 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"); } @@ -1817,15 +1429,15 @@ static void add_edge_constraint( /* '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); + 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, "ccw(vcur, v2, vc) = %d\n", ccw1); + fprintf(stderr, "orient2d(vcur, v2, vc) = %g\n", orient1); } #endif - if (ccw1 == -1 && (vc->visit_index != visit || tout->next->face->visit_index != visit)) { + if (orient1 < 0.0) { /* vcur--v2 should intersect vb--vc. */ #ifdef DEBUG_CDT if (dbg_level > 1) { @@ -1835,7 +1447,7 @@ static void add_edge_constraint( t = tout->next; state_through_vert = false; } - else if (ccw1 == 1 && tout->face->visit_index != visit) { + else if (orient1 > 0.0) { /* vcur--v2 should intersect va--vc. */ #ifdef DEBUG_CDT if (dbg_level > 1) { @@ -1845,7 +1457,7 @@ static void add_edge_constraint( t = tout; state_through_vert = false; } - else if (ccw1 == 0 && vc->visit_index != visit) { + else if (orient1 == 0.0) { #ifdef DEBUG_CDT if (dbg_level > 1) { fprintf(stderr, "ccw==0 case, so going through or to vc\n"); @@ -1856,16 +1468,11 @@ static void add_edge_constraint( } else { #ifdef DEBUG_CDT - fprintf(stderr, "add_edge_constraint desperation search 2\n"); + fprintf(stderr, "add_edge_constraint desperation search needed\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) { + if (++search_count > 1000000) { fprintf(stderr, "infinite loop? bailing out\n"); BLI_assert(0); /* Catch these while developing. */ break; @@ -2091,6 +1698,304 @@ static void dissolve_symedge(CDT_state *cdt, SymEdge *se) 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, picke 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) { @@ -2176,7 +2081,7 @@ static void remove_non_constraint_edges_leave_valid_bmesh(CDT_state *cdt) e = sorted_edges[i].e; se = &e->symedges[0]; dissolve = true; - if (!edge_touches_frame(e)) { + if (true /*!edge_touches_frame(e)*/) { fleft = se->face; fright = sym(se)->face; if (fleft != cdt->outer_face && fright != cdt->outer_face && @@ -2197,7 +2102,7 @@ static void remove_non_constraint_edges_leave_valid_bmesh(CDT_state *cdt) } } -static void remove_outer_edges(CDT_state *cdt, const bool remove_until_constraints) +static void remove_outer_edges_until_constraints(CDT_state *cdt) { LinkNode *fstack = NULL; SymEdge *se, *se_start; @@ -2207,27 +2112,28 @@ static void remove_outer_edges(CDT_state *cdt, const bool remove_until_constrain int dbg_level = 0; if (dbg_level > 0) { - fprintf(stderr, "remove_outer_edges, until_constraints=%d\n", remove_until_constraints); + fprintf(stderr, "remove_outer_edges_until_constraints\n"); } #endif cdt->outer_face->visit_index = visit; - - /* Find an f, not outer face, but touching outer face. */ - f = NULL; - se_start = se = cdt->vert_array[0]->symedge; + /* Walk around outer face, adding faces on other side of dissolvable edges to stack. */ + se_start = se = cdt->outer_face->symedge; do { - if (se->face != cdt->outer_face) { - f = se->face; - break; + 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); + } } - se = se->rot; - } while (se != se_start); - BLI_assert(f != NULL && f->symedge != NULL); - if (f == NULL) { - return; - } - BLI_linklist_prepend_pool(&fstack, f, cdt->listpool); + } while ((se = se->next) != se_start); + while (fstack != NULL) { LinkNode *to_dissolve = NULL; bool dissolvable; @@ -2245,18 +2151,16 @@ static void remove_outer_edges(CDT_state *cdt, const bool remove_until_constrain if (dbg_level > 0) { fprintf(stderr, "top of loop, f=%p\n", f); dump_se_cycle(f->symedge, "visit", 10000); - dump_cdt(cdt, "cdt at top of loop"); + 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 { - if (remove_until_constraints) { - dissolvable = !is_constrained_edge(se->edge); - } - else { - dissolvable = edge_touches_frame(se->edge); - } + dissolvable = !is_constrained_edge(se->edge); #ifdef DEBUG_CDT if (dbg_level > 1) { dump_se(se, "edge in f"); @@ -2306,15 +2210,20 @@ static void prepare_cdt_for_output(CDT_state *cdt, const CDT_output_type output_ 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 (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]; + 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 @@ -2335,19 +2244,20 @@ static void prepare_cdt_for_output(CDT_state *cdt, const CDT_output_type output_ else if (output_type == CDT_CONSTRAINTS_VALID_BMESH) { remove_non_constraint_edges_leave_valid_bmesh(cdt); } - else if (output_type == CDT_FULL || output_type == CDT_INSIDE) { - remove_outer_edges(cdt, output_type == CDT_INSIDE); + else if (output_type == CDT_INSIDE) { + remove_outer_edges_until_constraints(cdt); } } -#define NUM_BOUND_VERTS 4 -#define VERT_OUT_INDEX(v) ((v)->index - NUM_BOUND_VERTS) - -static CDT_result *cdt_get_output(CDT_state *cdt, const CDT_output_type output_type) +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; @@ -2356,35 +2266,70 @@ static CDT_result *cdt_get_output(CDT_state *cdt, const CDT_output_type output_t 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 except first NUM_BOUND_VERTS will be output. */ - nv = cdt->vert_array_len - NUM_BOUND_VERTS; + /* 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 < nv; i++) { - orig_map_size += BLI_linklist_count(cdt->vert_array[i + 4]->input_ids); + 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; - for (i = 0; i < nv; i++) { - j = i + NUM_BOUND_VERTS; - result->vert_coords[i][0] = (float)cdt->vert_array[j]->co[0]; - result->vert_coords[i][1] = (float)cdt->vert_array[j]->co[1]; - result->verts_orig_start_table[i] = orig_map_index; - for (ln = cdt->vert_array[j]->input_ids; ln; ln = ln->next) { - result->verts_orig[orig_map_index++] = POINTER_AS_INT(ln->link); + 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++; } - result->verts_orig_len_table[i] = orig_map_index - result->verts_orig_start_table[i]; } ne = 0; @@ -2412,8 +2357,8 @@ static CDT_result *cdt_get_output(CDT_state *cdt, const CDT_output_type output_t for (lne = cdt->edges; lne; lne = lne->next) { e = (CDTEdge *)lne->link; if (!is_deleted_edge(e)) { - result->edges[i][0] = VERT_OUT_INDEX(e->symedges[0].vert); - result->edges[i][1] = VERT_OUT_INDEX(e->symedges[1].vert); + 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); @@ -2462,7 +2407,7 @@ static CDT_result *cdt_get_output(CDT_state *cdt, const CDT_output_type output_t result->faces_start_table[i] = j; se = se_start = f->symedge; do { - result->faces[j++] = VERT_OUT_INDEX(se->vert); + result->faces[j++] = se->vert->index; se = se->next; } while (se != se_start); result->faces_len_table[i] = j - result->faces_start_table[i]; @@ -2478,28 +2423,72 @@ static CDT_result *cdt_get_output(CDT_state *cdt, const CDT_output_type output_t 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; - double epsilon = (double)input->epsilon; - int i, f, v1, v2; - int fedge_start, fedge_end; - double minx, maxx, miny, maxy; - float *xy; - double vert_co[2]; + int i, iv1, iv2, f, fedge_start, fedge_end; CDT_state *cdt; - CDT_result *result; - CDTVert **verts; - LinkNode *edge_list; + 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); } @@ -2514,77 +2503,43 @@ CDT_result *BLI_delaunay_2d_cdt_calc(const CDT_input *input, const CDT_output_ty return NULL; } - if (nv > 0) { - minx = miny = DBL_MAX; - maxx = maxy = -DBL_MAX; - for (i = 0; i < nv; i++) { - xy = input->vert_coords[i]; - if (xy[0] < minx) { - minx = xy[0]; - } - if (xy[0] > maxx) { - maxx = xy[0]; - } - if (xy[1] < miny) { - miny = xy[1]; - } - if (xy[1] > maxy) { - maxy = xy[1]; - } - } - verts = (CDTVert **)MEM_mallocN(nv * sizeof(CDTVert *), "constrained delaunay"); - } - else { - minx = miny = maxx = maxy = 0; - verts = NULL; - } - - if (epsilon == 0.0) { - epsilon = 1e-8; - } - cdt = cdt_init(minx, maxx, miny, maxy, epsilon); - /* TODO: use a random permutation for order of adding the vertices. */ - for (i = 0; i < nv; i++) { - vert_co[0] = (double)input->vert_coords[i][0]; - vert_co[1] = (double)input->vert_coords[i][1]; - verts[i] = add_point_constraint(cdt, vert_co, i); + cdt = new_cdt_init(input); + initial_triangulation(cdt); #ifdef DEBUG_CDT - if (dbg_level > 3) { - char namebuf[60]; - sprintf(namebuf, "after point %d = (%f,%f)\n", i, vert_co[0], vert_co[1]); - cdt_draw(cdt, namebuf); - dump_cdt(cdt, namebuf); - validate_cdt(cdt, true); - } -#endif + if (dbg_level > 0) { + validate_cdt(cdt, true, false, false); } +#endif + for (i = 0; i < ne; i++) { - v1 = input->edges[i][0]; - v2 = input->edges[i][1]; - if (v1 < 0 || v1 >= nv || v2 < 0 || v2 >= nv) { + 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 not valid: v1=%d, v2=%d, nv=%d\n", v1, v2, nv); + fprintf(stderr, "edge indices for e%d not valid: v1=%d, v2=%d, nv=%d\n", i, iv1, iv2, nv); #endif continue; } - add_edge_constraint(cdt, verts[v1], verts[v2], i, NULL); + 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, v1, v2); + sprintf(namebuf, "after edge constraint %d = (%d,%d)\n", i, iv1, iv2); cdt_draw(cdt, namebuf); - dump_cdt(cdt, namebuf); - validate_cdt(cdt, true); + // dump_cdt(cdt, namebuf); + validate_cdt(cdt, true, true, false); } #endif } -#ifdef DEBUG_CDT - if (dbg_level > 2) { - cdt_draw(cdt, "after edge constraints"); - dump_cdt(cdt, "after edge constraints"); - validate_cdt(cdt, true); - } -#endif + cdt->face_edge_offset = ne; for (f = 0; f < nf; f++) { int flen = input->faces_len_table[f]; @@ -2597,17 +2552,25 @@ CDT_result *BLI_delaunay_2d_cdt_calc(const CDT_input *input, const CDT_output_ty } for (i = 0; i < flen; i++) { int face_edge_id = cdt->face_edge_offset + fstart + i; - v1 = input->faces[fstart + i]; - v2 = input->faces[fstart + ((i + 1) % flen)]; - if (v1 < 0 || v1 >= nv || v2 < 0 || v2 >= nv) { + 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, v1=%d, v2=%d, nv=%d\n", f, v1, v2, nv); + fprintf(stderr, "face indices not valid: f=%d, iv1=%d, iv2=%d, nv=%d\n", f, iv1, iv2, nv); #endif continue; } - add_edge_constraint(cdt, verts[v1], verts[v2], face_edge_id, &edge_list); + 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 > 1) { + 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; @@ -2619,16 +2582,18 @@ CDT_result *BLI_delaunay_2d_cdt_calc(const CDT_input *input, const CDT_output_ty } if (dbg_level > 2) { cdt_draw(cdt, "after a face edge"); - dump_cdt(cdt, "after a face edge"); - validate_cdt(cdt, true); + 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 != verts[v1]) { + if (face_symedge->vert != v1) { face_symedge = &face_edge->symedges[1]; - BLI_assert(face_symedge->vert == verts[v1]); + BLI_assert(face_symedge->vert == v1); } } BLI_linklist_free_pool(edge_list, NULL, cdt->listpool); @@ -2639,22 +2604,36 @@ CDT_result *BLI_delaunay_2d_cdt_calc(const CDT_input *input, const CDT_output_ty } #ifdef DEBUG_CDT if (dbg_level > 0) { - validate_cdt(cdt, true); + validate_cdt(cdt, true, true, false); } if (dbg_level > 1) { - cdt_draw(cdt, "before cdt_get_output"); + cdt_draw(cdt, "after adding edges and faces"); + if (dbg_level > 2) { + dump_cdt(cdt, "after adding edges and faces"); + } } #endif - result = cdt_get_output(cdt, output_type); + + 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 - if (verts) { - MEM_freeN(verts); - } - cdt_free(cdt); + + new_cdt_free(cdt); return result; } @@ -2710,19 +2689,22 @@ void BLI_delaunay_2d_cdt_free(CDT_result *result) #ifdef DEBUG_CDT -static const char *vertname(const CDTVert *v) +ATTU static const char *vertname(const CDTVert *v) { static char vertnamebuf[20]; - if (v->index < 4) { - sprintf(vertnamebuf, "[%c]", "ABCD"[v->index]); - } - else { - sprintf(vertnamebuf, "[%d]", v->index - 4); - } + 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)); @@ -2732,11 +2714,12 @@ static void dump_se(const SymEdge *se, const char *lab) { if (se->next) { fprintf(stderr, - "%s%s((%.3f,%.3f)->(%.3f,%.3f))\n", + "%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)); @@ -2791,13 +2774,28 @@ static void dump_cdt_filtered(const CDT_state *cdt, 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)); + 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; @@ -2839,15 +2837,13 @@ static void dump_cdt_filtered(const CDT_state *cdt, continue; } if (f == cdt->outer_face) { - fprintf(stderr, "outer"); - } - else { - fprintf(stderr, "%x: centroid (%f,%f)", PL(f), F2(f->centroid)); + 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); @@ -2895,6 +2891,37 @@ static bool reachable_filter(const CDT_state *cdt, int v_index, void *filter_dat 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; @@ -2903,7 +2930,7 @@ static void dump_cdt_vert_neighborhood(CDT_state *cdt, int v, int maxdist, const 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. @@ -2916,24 +2943,20 @@ static void dump_cdt_vert_neighborhood(CDT_state *cdt, int v, int maxdist, const # define THICK_LINE 3 # define VERT_RADIUS 3 # define DRAW_VERT_LABELS 1 -static void cdt_draw(CDT_state *cdt, const char *lab) +# 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"); - double draw_margin = (cdt->maxx - cdt->minx + cdt->maxy - cdt->miny + 1) * 0.05; - double minx = cdt->minx - draw_margin; - double maxx = cdt->maxx + draw_margin; - double miny = cdt->miny - draw_margin; - double maxy = cdt->maxy + draw_margin; - double width, height, aspect; int view_width, view_height; - double scale; + double width, height, aspect, scale; LinkNode *ln; CDTVert *v, *u; 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; @@ -2979,10 +3002,23 @@ static void cdt_draw(CDT_state *cdt, const char *lab) fprintf(f, " <title>%s", vertname(u)); fprintf(f, "%s</title>\n", vertname(v)); fprintf(f, "</line>\n"); +# if DRAW_EDGE_LABELS + fprintf(f, + "<text x=\"%f\" y=\"%f\" font-size=\"small\">", + SX(0.5 * (u->co[0] + v->co[0])), + SY(0.5 * (u->co[1] + v->co[1]))); + fprintf(f, "%s", vertname(u)); + fprintf(f, "%s", vertname(v)); + fprintf(f, "%s", sename(&e->symedges[0])); + fprintf(f, "%s</text>\n", sename(&e->symedges[1])); +# endif } - i = cdt->output_prepared ? NUM_BOUND_VERTS : 0; + i = 0; for (; i < cdt->vert_array_len; i++) { v = cdt->vert_array[i]; + if (v->merge_to_index != -1) { + continue; + } fprintf(f, "<circle fill=\"black\" cx=\"%f\" cy=\"%f\" r=\"%d\">\n", SX(v->co[0]), @@ -3006,6 +3042,25 @@ static void cdt_draw(CDT_state *cdt, const char *lab) # undef SY } +static void cdt_draw(CDT_state *cdt, const char *lab) +{ + double draw_margin, minx, maxx, miny, maxy; + + set_min_max(cdt); + draw_margin = (cdt->maxx - cdt->minx + cdt->maxy - cdt->miny + 1) * 0.05; + minx = cdt->minx - draw_margin; + maxx = cdt->maxx + draw_margin; + miny = cdt->miny - draw_margin; + maxy = cdt->maxy + draw_margin; + cdt_draw_region(cdt, lab, minx, miny, maxx, maxy); +} + +static void cdt_draw_vertex_region(CDT_state *cdt, int v, double dist, const char *lab) +{ + const double *co = cdt->vert_array[v]->co; + cdt_draw_region(cdt, lab, co[0] - dist, co[1] - dist, co[0] + dist, co[1] + dist); +} + # define CDTFILE "/tmp/cdtinput.txt" static void write_cdt_input_to_file(const CDT_input *inp) { @@ -3029,16 +3084,18 @@ static void write_cdt_input_to_file(const CDT_input *inp) } # ifndef NDEBUG /* Only used in assert. */ -/** +/* * Is a visible from b: i.e., ab crosses no edge of cdt? * If constrained is true, consider only constrained edges as possible crossers. * In any case, don't count an edge ab itself. + * Note: this is an expensive test if there are a lot of edges. */ static bool is_visible(const CDTVert *a, const CDTVert *b, bool constrained, const CDT_state *cdt) { const LinkNode *ln; const CDTEdge *e; const SymEdge *se, *senext; + double lambda, mu; int ikind; for (ln = cdt->edges; ln; ln = ln->next) { @@ -3055,13 +3112,17 @@ static bool is_visible(const CDTVert *a, const CDTVert *b, bool constrained, con continue; } ikind = isect_seg_seg_v2_lambda_mu_db( - a->co, b->co, se->vert->co, senext->vert->co, NULL, NULL); + a->co, b->co, se->vert->co, senext->vert->co, &lambda, &mu); if (ikind != ISECT_LINE_LINE_NONE) { if (ikind == ISECT_LINE_LINE_COLINEAR) { /* TODO: special test here for overlap. */ continue; } - return false; + /* Allow an intersection very near or at ends, to allow for numerical error. */ + if (lambda > FLT_EPSILON && (1.0 - lambda) > FLT_EPSILON && mu > FLT_EPSILON && + (1.0 - mu) > FLT_EPSILON) { + return false; + } } } return true; @@ -3069,7 +3130,7 @@ static bool is_visible(const CDTVert *a, const CDTVert *b, bool constrained, con # endif # ifndef NDEBUG /* Only used in assert. */ -/** +/* * Check that edge ab satisfies constrained delaunay condition: * That is, for all non-constraint, non-border edges ab, * (1) ab is visible in the constraint graph; and @@ -3078,7 +3139,7 @@ static bool is_visible(const CDTVert *a, const CDTVert *b, bool constrained, con * The argument 'se' specifies ab by: a is se's vert and b is se->next's vert. * Return true if check is OK. */ -static bool is_delaunay_edge(const SymEdge *se, const double epsilon) +static bool is_delaunay_edge(const SymEdge *se) { int i; CDTVert *a, *b, *c; @@ -3099,7 +3160,7 @@ static bool is_delaunay_edge(const SymEdge *se, const double epsilon) b = curse->next->vert; c = curse->next->next->vert; for (ss = curse->rot; ss != curse; ss = ss->rot) { - ok[i] |= delaunay_check(a, b, c, ss->next->vert, epsilon); + ok[i] |= incircle(a->co, b->co, c->co, ss->next->vert->co) <= 0.0; } } return ok[0] || ok[1]; @@ -3113,48 +3174,29 @@ static bool plausible_non_null_ptr(void *p) } # endif -static void validate_face_centroid(SymEdge *se) -{ - SymEdge *senext; -# ifndef NDEBUG - double *centroidp = se->face->centroid; -# endif - double c[2]; - int count; - copy_v2_v2_db(c, se->vert->co); - BLI_assert(reachable(se->next, se, 100)); - count = 1; - for (senext = se->next; senext != se; senext = senext->next) { - add_v2_v2_db(c, senext->vert->co); - count++; - } - c[0] /= count; - c[1] /= count; - BLI_assert(fabs(c[0] - centroidp[0]) < 1e-8 && fabs(c[1] - centroidp[1]) < 1e-8); -} - -static void validate_cdt(CDT_state *cdt, bool check_all_tris) +static void validate_cdt(CDT_state *cdt, + bool check_all_tris, + bool check_delaunay, + bool check_visibility) { - LinkNode *ln, *lne; - int totedges, totfaces, totverts, totborderedges; + LinkNode *ln; + int totedges, totfaces, totverts; CDTEdge *e; SymEdge *se, *sesym, *s; CDTVert *v, *v1, *v2, *v3; CDTFace *f; - double *p; - double margin; int i, limit; bool isborder; if (cdt->output_prepared) { return; } + if (cdt->edges == NULL || cdt->edges->next == NULL) { + return; + } BLI_assert(cdt != NULL); - BLI_assert(cdt->maxx >= cdt->minx); - BLI_assert(cdt->maxy >= cdt->miny); totedges = 0; - totborderedges = 0; for (ln = cdt->edges; ln; ln = ln->next) { e = (CDTEdge *)ln->link; se = &e->symedges[0]; @@ -3165,13 +3207,6 @@ static void validate_cdt(CDT_state *cdt, bool check_all_tris) } totedges++; isborder = is_border_edge(e, cdt); - if (isborder) { - totborderedges++; - BLI_assert((se->face == cdt->outer_face && sesym->face != cdt->outer_face) || - (se->face != cdt->outer_face && sesym->face == cdt->outer_face)); - } - /* BLI_assert(se->face != sesym->face); - * Not required because faces can have intruding wire edges. */ BLI_assert(se->vert != sesym->vert); BLI_assert(se->edge == sesym->edge && se->edge == e); BLI_assert(sym(se) == sesym && sym(sesym) == se); @@ -3179,8 +3214,6 @@ static void validate_cdt(CDT_state *cdt, bool check_all_tris) se = &e->symedges[i]; v = se->vert; f = se->face; - p = v->co; - UNUSED_VARS_NDEBUG(p); BLI_assert(plausible_non_null_ptr(v)); if (f != NULL) { BLI_assert(plausible_non_null_ptr(f)); @@ -3198,9 +3231,9 @@ static void validate_cdt(CDT_state *cdt, bool check_all_tris) v1 = se->vert; v2 = se->next->vert; v3 = se->next->next->vert; - BLI_assert(CCW_test(v1->co, v2->co, v3->co, 0.0)); - BLI_assert(CCW_test(v2->co, v3->co, v1->co, 0.0)); - BLI_assert(CCW_test(v3->co, v1->co, v2->co, 0.0)); + BLI_assert(orient2d(v1->co, v2->co, v3->co) >= 0.0); + BLI_assert(orient2d(v2->co, v3->co, v1->co) >= 0.0); + BLI_assert(orient2d(v3->co, v1->co, v2->co) >= 0.0); } UNUSED_VARS_NDEBUG(limit); BLI_assert(se->next->next != se); @@ -3211,19 +3244,23 @@ static void validate_cdt(CDT_state *cdt, bool check_all_tris) s = s->next; } while (s != se); } - BLI_assert(isborder || is_visible(se->vert, se->next->vert, false, cdt)); - BLI_assert(isborder || is_delaunay_edge(se, cdt->epsilon)); + if (check_visibility) { + BLI_assert(isborder || is_visible(se->vert, se->next->vert, false, cdt)); + } + if (!isborder && check_delaunay) { + BLI_assert(is_delaunay_edge(se)); + } } totverts = 0; - margin = cdt->margin; for (i = 0; i < cdt->vert_array_len; i++) { - totverts++; v = cdt->vert_array[i]; BLI_assert(plausible_non_null_ptr(v)); - p = v->co; - BLI_assert(p[0] >= cdt->minx - margin && p[0] <= cdt->maxx + margin); - UNUSED_VARS_NDEBUG(margin); - BLI_assert(v->symedge->vert == v); + if (v->merge_to_index != -1) { + BLI_assert(v->merge_to_index >= 0 && v->merge_to_index < cdt->vert_array_len); + continue; + } + totverts++; + BLI_assert(cdt->vert_array_len <= 1 || v->symedge->vert == v); } totfaces = 0; for (ln = cdt->faces; ln; ln = ln->next) { @@ -3236,23 +3273,1111 @@ static void validate_cdt(CDT_state *cdt, bool check_all_tris) if (f == cdt->outer_face) { continue; } - for (lne = cdt->edges; lne; lne = lne->next) { - e = (CDTEdge *)lne->link; - if (!is_deleted_edge(e)) { - for (i = 0; i < 2; i++) { - if (e->symedges[i].face == f) { - validate_face_centroid(&e->symedges[i]); - } - } - } - } - p = f->centroid; - BLI_assert(p[0] >= cdt->minx - margin && p[0] <= cdt->maxx + margin); - BLI_assert(p[1] >= cdt->miny - margin && p[1] <= cdt->maxy + margin); } /* Euler's formula for planar graphs. */ - if (check_all_tris) { + if (check_all_tris && totfaces > 1) { BLI_assert(totverts - totedges + totfaces == 2); } } #endif + +/* Jonathan Shewchuk's adaptive predicates, trimmed to those needed here. + * Permission obtained by private communication from Jonathan to include this code in Blender. + */ + +/* + * Routines for Arbitrary Precision Floating-point Arithmetic + * and Fast Robust Geometric Predicates + * (predicates.c) + * + * May 18, 1996 + * + * Placed in the public domain by + * Jonathan Richard Shewchuk + * School of Computer Science + * Carnegie Mellon University + * 5000 Forbes Avenue + * Pittsburgh, Pennsylvania 15213-3891 + * jrs@cs.cmu.edu + * + * This file contains C implementation of algorithms for exact addition + * and multiplication of floating-point numbers, and predicates for + * robustly performing the orientation and incircle tests used in + * computational geometry. The algorithms and underlying theory are + * described in Jonathan Richard Shewchuk. "Adaptive Precision Floating- + * Point Arithmetic and Fast Robust Geometric Predicates." Technical + * Report CMU-CS-96-140, School of Computer Science, Carnegie Mellon + * University, Pittsburgh, Pennsylvania, May 1996. (Submitted to + * Discrete & Computational Geometry.) + * + * This file, the paper listed above, and other information are available + * from the Web page http://www.cs.cmu.edu/~quake/robust.html . + * + * Using this code: + * + * First, read the short or long version of the paper (from the Web page + * above). + * + * Be sure to call exactinit() once, before calling any of the arithmetic + * functions or geometric predicates. Also be sure to turn on the + * optimizer when compiling this file. + * + * On some machines, the exact arithmetic routines might be defeated by the + * use of internal extended precision floating-point registers. Sometimes + * this problem can be fixed by defining certain values to be volatile, + * thus forcing them to be stored to memory and rounded off. This isn't + * a great solution, though, as it slows the arithmetic down. + * + * To try this out, write "#define INEXACT volatile" below. Normally, + * however, INEXACT should be defined to be nothing. ("#define INEXACT".) + */ + +#define INEXACT /* Nothing */ +/* #define INEXACT volatile */ + +/* Which of the following two methods of finding the absolute values is + * fastest is compiler-dependent. A few compilers can inline and optimize + * the fabs() call; but most will incur the overhead of a function call, + * which is disastrously slow. A faster way on IEEE machines might be to + * mask the appropriate bit, but that's difficult to do in C. + */ + +#define Absolute(a) ((a) >= 0.0 ? (a) : -(a)) +/* #define Absolute(a) fabs(a) */ + +/* Many of the operations are broken up into two pieces, a main part that + * performs an approximate operation, and a "tail" that computes the + * roundoff error of that operation. + * + * The operations Fast_Two_Sum(), Fast_Two_Diff(), Two_Sum(), Two_Diff(), + * Split(), and Two_Product() are all implemented as described in the + * reference. Each of these macros requires certain variables to be + * defined in the calling routine. The variables `bvirt', `c', `abig', + * `_i', `_j', `_k', `_l', `_m', and `_n' are declared `INEXACT' because + * they store the result of an operation that may incur roundoff error. + * The input parameter `x' (or the highest numbered `x_' parameter) must + * also be declared `INEXACT'. + */ + +#define Fast_Two_Sum_Tail(a, b, x, y) \ + bvirt = x - a; \ + y = b - bvirt + +#define Fast_Two_Sum(a, b, x, y) \ + x = (double)(a + b); \ + Fast_Two_Sum_Tail(a, b, x, y) + +#define Fast_Two_Diff_Tail(a, b, x, y) \ + bvirt = a - x; \ + y = bvirt - b + +#define Fast_Two_Diff(a, b, x, y) \ + x = (double)(a - b); \ + Fast_Two_Diff_Tail(a, b, x, y) + +#define Two_Sum_Tail(a, b, x, y) \ + bvirt = (double)(x - a); \ + avirt = x - bvirt; \ + bround = b - bvirt; \ + around = a - avirt; \ + y = around + bround + +#define Two_Sum(a, b, x, y) \ + x = (double)(a + b); \ + Two_Sum_Tail(a, b, x, y) + +#define Two_Diff_Tail(a, b, x, y) \ + bvirt = (double)(a - x); \ + avirt = x + bvirt; \ + bround = bvirt - b; \ + around = a - avirt; \ + y = around + bround + +#define Two_Diff(a, b, x, y) \ + x = (double)(a - b); \ + Two_Diff_Tail(a, b, x, y) + +#define Split(a, ahi, alo) \ + c = (double)(splitter * a); \ + abig = (double)(c - a); \ + ahi = c - abig; \ + alo = a - ahi + +#define Two_Product_Tail(a, b, x, y) \ + Split(a, ahi, alo); \ + Split(b, bhi, blo); \ + err1 = x - (ahi * bhi); \ + err2 = err1 - (alo * bhi); \ + err3 = err2 - (ahi * blo); \ + y = (alo * blo) - err3 + +#define Two_Product(a, b, x, y) \ + x = (double)(a * b); \ + Two_Product_Tail(a, b, x, y) + +#define Two_Product_Presplit(a, b, bhi, blo, x, y) \ + x = (double)(a * b); \ + Split(a, ahi, alo); \ + err1 = x - (ahi * bhi); \ + err2 = err1 - (alo * bhi); \ + err3 = err2 - (ahi * blo); \ + y = (alo * blo) - err3 + +#define Square_Tail(a, x, y) \ + Split(a, ahi, alo); \ + err1 = x - (ahi * ahi); \ + err3 = err1 - ((ahi + ahi) * alo); \ + y = (alo * alo) - err3 + +#define Square(a, x, y) \ + x = (double)(a * a); \ + Square_Tail(a, x, y) + +#define Two_One_Sum(a1, a0, b, x2, x1, x0) \ + Two_Sum(a0, b, _i, x0); \ + Two_Sum(a1, _i, x2, x1) + +#define Two_One_Diff(a1, a0, b, x2, x1, x0) \ + Two_Diff(a0, b, _i, x0); \ + Two_Sum(a1, _i, x2, x1) + +#define Two_Two_Sum(a1, a0, b1, b0, x3, x2, x1, x0) \ + Two_One_Sum(a1, a0, b0, _j, _0, x0); \ + Two_One_Sum(_j, _0, b1, x3, x2, x1) + +#define Two_Two_Diff(a1, a0, b1, b0, x3, x2, x1, x0) \ + Two_One_Diff(a1, a0, b0, _j, _0, x0); \ + Two_One_Diff(_j, _0, b1, x3, x2, x1) + +static double splitter; /* = 2^ceiling(p / 2) + 1. Used to split floats in half. */ +static double m_epsilon; /* = 2^(-p). Used to estimate roundoff errors. */ +/* A set of coefficients used to calculate maximum roundoff errors. */ +static double resulterrbound; +static double ccwerrboundA, ccwerrboundB, ccwerrboundC; +static double o3derrboundA, o3derrboundB, o3derrboundC; +static double iccerrboundA, iccerrboundB, iccerrboundC; +static double isperrboundA, isperrboundB, isperrboundC; + +/* exactinit() Initialize the variables used for exact arithmetic. + * + * `epsilon' is the largest power of two such that 1.0 + epsilon = 1.0 in + * floating-point arithmetic. `epsilon' bounds the relative roundoff + * error. It is used for floating-point error analysis. + * + * `splitter' is used to split floating-point numbers into two half- + * length significands for exact multiplication. + * + * I imagine that a highly optimizing compiler might be too smart for its + * own good, and somehow cause this routine to fail, if it pretends that + * floating-point arithmetic is too much like real arithmetic. + * + * Don't change this routine unless you fully understand it. + */ + +static void exactinit(void) +{ + double half; + double check, lastcheck; + int every_other; + + every_other = 1; + half = 0.5; + m_epsilon = 1.0; + splitter = 1.0; + check = 1.0; + /* Repeatedly divide `epsilon' by two until it is too small to add to + * one without causing roundoff. (Also check if the sum is equal to + * the previous sum, for machines that round up instead of using exact + * rounding. Not that this library will work on such machines anyway. + */ + do { + lastcheck = check; + m_epsilon *= half; + if (every_other) { + splitter *= 2.0; + } + every_other = !every_other; + check = 1.0 + m_epsilon; + } while ((check != 1.0) && (check != lastcheck)); + splitter += 1.0; + + /* Error bounds for orientation and incircle tests. */ + resulterrbound = (3.0 + 8.0 * m_epsilon) * m_epsilon; + ccwerrboundA = (3.0 + 16.0 * m_epsilon) * m_epsilon; + ccwerrboundB = (2.0 + 12.0 * m_epsilon) * m_epsilon; + ccwerrboundC = (9.0 + 64.0 * m_epsilon) * m_epsilon * m_epsilon; + o3derrboundA = (7.0 + 56.0 * m_epsilon) * m_epsilon; + o3derrboundB = (3.0 + 28.0 * m_epsilon) * m_epsilon; + o3derrboundC = (26.0 + 288.0 * m_epsilon) * m_epsilon * m_epsilon; + iccerrboundA = (10.0 + 96.0 * m_epsilon) * m_epsilon; + iccerrboundB = (4.0 + 48.0 * m_epsilon) * m_epsilon; + iccerrboundC = (44.0 + 576.0 * m_epsilon) * m_epsilon * m_epsilon; + isperrboundA = (16.0 + 224.0 * m_epsilon) * m_epsilon; + isperrboundB = (5.0 + 72.0 * m_epsilon) * m_epsilon; + isperrboundC = (71.0 + 1408.0 * m_epsilon) * m_epsilon * m_epsilon; +} + +/* fast_expansion_sum_zeroelim() Sum two expansions, eliminating zero + * components from the output expansion. + * + * Sets h = e + f. See the long version of my paper for details. + * + * If round-to-even is used (as with IEEE 754), maintains the strongly + * nonoverlapping property. (That is, if e is strongly nonoverlapping, h + * will be also.) Does NOT maintain the nonoverlapping or nonadjacent + * properties. + */ + +static int fast_expansion_sum_zeroelim( + int elen, double *e, int flen, double *f, double *h) /* h cannot be e or f. */ +{ + double Q; + INEXACT double Qnew; + INEXACT double hh; + INEXACT double bvirt; + double avirt, bround, around; + int eindex, findex, hindex; + double enow, fnow; + + enow = e[0]; + fnow = f[0]; + eindex = findex = 0; + if ((fnow > enow) == (fnow > -enow)) { + Q = enow; + enow = e[++eindex]; + } + else { + Q = fnow; + fnow = f[++findex]; + } + hindex = 0; + if ((eindex < elen) && (findex < flen)) { + if ((fnow > enow) == (fnow > -enow)) { + Fast_Two_Sum(enow, Q, Qnew, hh); + enow = e[++eindex]; + } + else { + Fast_Two_Sum(fnow, Q, Qnew, hh); + fnow = f[++findex]; + } + Q = Qnew; + if (hh != 0.0) { + h[hindex++] = hh; + } + while ((eindex < elen) && (findex < flen)) { + if ((fnow > enow) == (fnow > -enow)) { + Two_Sum(Q, enow, Qnew, hh); + enow = e[++eindex]; + } + else { + Two_Sum(Q, fnow, Qnew, hh); + fnow = f[++findex]; + } + Q = Qnew; + if (hh != 0.0) { + h[hindex++] = hh; + } + } + } + while (eindex < elen) { + Two_Sum(Q, enow, Qnew, hh); + enow = e[++eindex]; + Q = Qnew; + if (hh != 0.0) { + h[hindex++] = hh; + } + } + while (findex < flen) { + Two_Sum(Q, fnow, Qnew, hh); + fnow = f[++findex]; + Q = Qnew; + if (hh != 0.0) { + h[hindex++] = hh; + } + } + if ((Q != 0.0) || (hindex == 0)) { + h[hindex++] = Q; + } + return hindex; +} + +/* scale_expansion_zeroelim() Multiply an expansion by a scalar, + * eliminating zero components from the + * output expansion. + * + * Sets h = be. See either version of my paper for details. + * + * Maintains the nonoverlapping property. If round-to-even is used (as + * with IEEE 754), maintains the strongly nonoverlapping and nonadjacent + * properties as well. (That is, if e has one of these properties, so + * will h.) + */ + +static int scale_expansion_zeroelim(int elen, + double *e, + double b, + double *h) /* e and h cannot be the same. */ +{ + INEXACT double Q, sum; + double hh; + INEXACT double product1; + double product0; + int eindex, hindex; + double enow; + INEXACT double bvirt; + double avirt, bround, around; + INEXACT double c; + INEXACT double abig; + double ahi, alo, bhi, blo; + double err1, err2, err3; + + Split(b, bhi, blo); + Two_Product_Presplit(e[0], b, bhi, blo, Q, hh); + hindex = 0; + if (hh != 0) { + h[hindex++] = hh; + } + for (eindex = 1; eindex < elen; eindex++) { + enow = e[eindex]; + Two_Product_Presplit(enow, b, bhi, blo, product1, product0); + Two_Sum(Q, product0, sum, hh); + if (hh != 0) { + h[hindex++] = hh; + } + Fast_Two_Sum(product1, sum, Q, hh); + if (hh != 0) { + h[hindex++] = hh; + } + } + if ((Q != 0.0) || (hindex == 0)) { + h[hindex++] = Q; + } + return hindex; +} + +/* estimate() Produce a one-word estimate of an expansion's value. + * + * See either version of my paper for details. + */ + +static double estimate(int elen, double *e) +{ + double Q; + int eindex; + + Q = e[0]; + for (eindex = 1; eindex < elen; eindex++) { + Q += e[eindex]; + } + return Q; +} + +/* orient2d() Adaptive exact 2D orientation test. Robust. + * + * Return a positive value if the points pa, pb, and pc occur + * in counterclockwise order; a negative value if they occur + * in clockwise order; and zero if they are collinear. The + * result is also a rough approximation of twice the signed + * area of the triangle defined by the three points. + * + * This uses exact arithmetic to ensure a correct answer. The + * result returned is the determinant of a matrix. + * This determinant is computed adaptively, in the sense that exact + * arithmetic is used only to the degree it is needed to ensure that the + * returned value has the correct sign. Hence, orient2d() is usually quite + * fast, but will run more slowly when the input points are collinear or + * nearly so. + */ + +static double orient2dadapt(const double *pa, const double *pb, const double *pc, double detsum) +{ + INEXACT double acx, acy, bcx, bcy; + double acxtail, acytail, bcxtail, bcytail; + INEXACT double detleft, detright; + double detlefttail, detrighttail; + double det, errbound; + double B[4], C1[8], C2[12], D[16]; + INEXACT double B3; + int C1length, C2length, Dlength; + double u[4]; + INEXACT double u3; + INEXACT double s1, t1; + double s0, t0; + + INEXACT double bvirt; + double avirt, bround, around; + INEXACT double c; + INEXACT double abig; + double ahi, alo, bhi, blo; + double err1, err2, err3; + INEXACT double _i, _j; + double _0; + + acx = (double)(pa[0] - pc[0]); + bcx = (double)(pb[0] - pc[0]); + acy = (double)(pa[1] - pc[1]); + bcy = (double)(pb[1] - pc[1]); + + Two_Product(acx, bcy, detleft, detlefttail); + Two_Product(acy, bcx, detright, detrighttail); + + Two_Two_Diff(detleft, detlefttail, detright, detrighttail, B3, B[2], B[1], B[0]); + B[3] = B3; + + det = estimate(4, B); + errbound = ccwerrboundB * detsum; + if ((det >= errbound) || (-det >= errbound)) { + return det; + } + + Two_Diff_Tail(pa[0], pc[0], acx, acxtail); + Two_Diff_Tail(pb[0], pc[0], bcx, bcxtail); + Two_Diff_Tail(pa[1], pc[1], acy, acytail); + Two_Diff_Tail(pb[1], pc[1], bcy, bcytail); + + if ((acxtail == 0.0) && (acytail == 0.0) && (bcxtail == 0.0) && (bcytail == 0.0)) { + return det; + } + + errbound = ccwerrboundC * detsum + resulterrbound * Absolute(det); + det += (acx * bcytail + bcy * acxtail) - (acy * bcxtail + bcx * acytail); + if ((det >= errbound) || (-det >= errbound)) { + return det; + } + + Two_Product(acxtail, bcy, s1, s0); + Two_Product(acytail, bcx, t1, t0); + Two_Two_Diff(s1, s0, t1, t0, u3, u[2], u[1], u[0]); + u[3] = u3; + C1length = fast_expansion_sum_zeroelim(4, B, 4, u, C1); + + Two_Product(acx, bcytail, s1, s0); + Two_Product(acy, bcxtail, t1, t0); + Two_Two_Diff(s1, s0, t1, t0, u3, u[2], u[1], u[0]); + u[3] = u3; + C2length = fast_expansion_sum_zeroelim(C1length, C1, 4, u, C2); + + Two_Product(acxtail, bcytail, s1, s0); + Two_Product(acytail, bcxtail, t1, t0); + Two_Two_Diff(s1, s0, t1, t0, u3, u[2], u[1], u[0]); + u[3] = u3; + Dlength = fast_expansion_sum_zeroelim(C2length, C2, 4, u, D); + + return (D[Dlength - 1]); +} + +static double orient2d(const double *pa, const double *pb, const double *pc) +{ + double detleft, detright, det; + double detsum, errbound; + + detleft = (pa[0] - pc[0]) * (pb[1] - pc[1]); + detright = (pa[1] - pc[1]) * (pb[0] - pc[0]); + det = detleft - detright; + + if (detleft > 0.0) { + if (detright <= 0.0) { + return det; + } + else { + detsum = detleft + detright; + } + } + else if (detleft < 0.0) { + if (detright >= 0.0) { + return det; + } + else { + detsum = -detleft - detright; + } + } + else { + return det; + } + + errbound = ccwerrboundA * detsum; + if ((det >= errbound) || (-det >= errbound)) { + return det; + } + + return orient2dadapt(pa, pb, pc, detsum); +} + +/* incircle() Adaptive exact 2D incircle test. Robust. + * + * Return a positive value if the point pd lies inside the + * circle passing through pa, pb, and pc; a negative value if + * it lies outside; and zero if the four points are cocircular. + * The points pa, pb, and pc must be in counterclockwise + * order, or the sign of the result will be reversed. + * + * This uses exact arithmetic to ensure a correct answer. + * The result returned is the determinant of a matrix. + * This determinant is computed adaptively, in the sense that exact + * arithmetic is used only to the degree it is needed to ensure that the + * returned value has the correct sign. Hence, incircle() is usually quite + * fast, but will run more slowly when the input points are cocircular or + * nearly so. + */ + +static double incircleadapt( + const double *pa, const double *pb, const double *pc, const double *pd, double permanent) +{ + INEXACT double adx, bdx, cdx, ady, bdy, cdy; + double det, errbound; + + INEXACT double bdxcdy1, cdxbdy1, cdxady1, adxcdy1, adxbdy1, bdxady1; + double bdxcdy0, cdxbdy0, cdxady0, adxcdy0, adxbdy0, bdxady0; + double bc[4], ca[4], ab[4]; + INEXACT double bc3, ca3, ab3; + double axbc[8], axxbc[16], aybc[8], ayybc[16], adet[32]; + int axbclen, axxbclen, aybclen, ayybclen, alen; + double bxca[8], bxxca[16], byca[8], byyca[16], bdet[32]; + int bxcalen, bxxcalen, bycalen, byycalen, blen; + double cxab[8], cxxab[16], cyab[8], cyyab[16], cdet[32]; + int cxablen, cxxablen, cyablen, cyyablen, clen; + double abdet[64]; + int ablen; + double fin1[1152], fin2[1152]; + double *finnow, *finother, *finswap; + int finlength; + + double adxtail, bdxtail, cdxtail, adytail, bdytail, cdytail; + INEXACT double adxadx1, adyady1, bdxbdx1, bdybdy1, cdxcdx1, cdycdy1; + double adxadx0, adyady0, bdxbdx0, bdybdy0, cdxcdx0, cdycdy0; + double aa[4], bb[4], cc[4]; + INEXACT double aa3, bb3, cc3; + INEXACT double ti1, tj1; + double ti0, tj0; + double u[4], v[4]; + INEXACT double u3, v3; + double temp8[8], temp16a[16], temp16b[16], temp16c[16]; + double temp32a[32], temp32b[32], temp48[48], temp64[64]; + int temp8len, temp16alen, temp16blen, temp16clen; + int temp32alen, temp32blen, temp48len, temp64len; + double axtbb[8], axtcc[8], aytbb[8], aytcc[8]; + int axtbblen, axtcclen, aytbblen, aytcclen; + double bxtaa[8], bxtcc[8], bytaa[8], bytcc[8]; + int bxtaalen, bxtcclen, bytaalen, bytcclen; + double cxtaa[8], cxtbb[8], cytaa[8], cytbb[8]; + int cxtaalen, cxtbblen, cytaalen, cytbblen; + double axtbc[8], aytbc[8], bxtca[8], bytca[8], cxtab[8], cytab[8]; + int axtbclen, aytbclen, bxtcalen, bytcalen, cxtablen, cytablen; + double axtbct[16], aytbct[16], bxtcat[16], bytcat[16], cxtabt[16], cytabt[16]; + int axtbctlen, aytbctlen, bxtcatlen, bytcatlen, cxtabtlen, cytabtlen; + double axtbctt[8], aytbctt[8], bxtcatt[8]; + double bytcatt[8], cxtabtt[8], cytabtt[8]; + int axtbcttlen, aytbcttlen, bxtcattlen, bytcattlen, cxtabttlen, cytabttlen; + double abt[8], bct[8], cat[8]; + int abtlen, bctlen, catlen; + double abtt[4], bctt[4], catt[4]; + int abttlen, bcttlen, cattlen; + INEXACT double abtt3, bctt3, catt3; + double negate; + + INEXACT double bvirt; + double avirt, bround, around; + INEXACT double c; + INEXACT double abig; + double ahi, alo, bhi, blo; + double err1, err2, err3; + INEXACT double _i, _j; + double _0; + + adx = (double)(pa[0] - pd[0]); + bdx = (double)(pb[0] - pd[0]); + cdx = (double)(pc[0] - pd[0]); + ady = (double)(pa[1] - pd[1]); + bdy = (double)(pb[1] - pd[1]); + cdy = (double)(pc[1] - pd[1]); + + Two_Product(bdx, cdy, bdxcdy1, bdxcdy0); + Two_Product(cdx, bdy, cdxbdy1, cdxbdy0); + Two_Two_Diff(bdxcdy1, bdxcdy0, cdxbdy1, cdxbdy0, bc3, bc[2], bc[1], bc[0]); + bc[3] = bc3; + axbclen = scale_expansion_zeroelim(4, bc, adx, axbc); + axxbclen = scale_expansion_zeroelim(axbclen, axbc, adx, axxbc); + aybclen = scale_expansion_zeroelim(4, bc, ady, aybc); + ayybclen = scale_expansion_zeroelim(aybclen, aybc, ady, ayybc); + alen = fast_expansion_sum_zeroelim(axxbclen, axxbc, ayybclen, ayybc, adet); + + Two_Product(cdx, ady, cdxady1, cdxady0); + Two_Product(adx, cdy, adxcdy1, adxcdy0); + Two_Two_Diff(cdxady1, cdxady0, adxcdy1, adxcdy0, ca3, ca[2], ca[1], ca[0]); + ca[3] = ca3; + bxcalen = scale_expansion_zeroelim(4, ca, bdx, bxca); + bxxcalen = scale_expansion_zeroelim(bxcalen, bxca, bdx, bxxca); + bycalen = scale_expansion_zeroelim(4, ca, bdy, byca); + byycalen = scale_expansion_zeroelim(bycalen, byca, bdy, byyca); + blen = fast_expansion_sum_zeroelim(bxxcalen, bxxca, byycalen, byyca, bdet); + + Two_Product(adx, bdy, adxbdy1, adxbdy0); + Two_Product(bdx, ady, bdxady1, bdxady0); + Two_Two_Diff(adxbdy1, adxbdy0, bdxady1, bdxady0, ab3, ab[2], ab[1], ab[0]); + ab[3] = ab3; + cxablen = scale_expansion_zeroelim(4, ab, cdx, cxab); + cxxablen = scale_expansion_zeroelim(cxablen, cxab, cdx, cxxab); + cyablen = scale_expansion_zeroelim(4, ab, cdy, cyab); + cyyablen = scale_expansion_zeroelim(cyablen, cyab, cdy, cyyab); + clen = fast_expansion_sum_zeroelim(cxxablen, cxxab, cyyablen, cyyab, cdet); + + ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet); + finlength = fast_expansion_sum_zeroelim(ablen, abdet, clen, cdet, fin1); + + det = estimate(finlength, fin1); + errbound = iccerrboundB * permanent; + if ((det >= errbound) || (-det >= errbound)) { + return det; + } + + Two_Diff_Tail(pa[0], pd[0], adx, adxtail); + Two_Diff_Tail(pa[1], pd[1], ady, adytail); + Two_Diff_Tail(pb[0], pd[0], bdx, bdxtail); + Two_Diff_Tail(pb[1], pd[1], bdy, bdytail); + Two_Diff_Tail(pc[0], pd[0], cdx, cdxtail); + Two_Diff_Tail(pc[1], pd[1], cdy, cdytail); + if ((adxtail == 0.0) && (bdxtail == 0.0) && (cdxtail == 0.0) && (adytail == 0.0) && + (bdytail == 0.0) && (cdytail == 0.0)) { + return det; + } + + errbound = iccerrboundC * permanent + resulterrbound * Absolute(det); + det += ((adx * adx + ady * ady) * + ((bdx * cdytail + cdy * bdxtail) - (bdy * cdxtail + cdx * bdytail)) + + 2.0 * (adx * adxtail + ady * adytail) * (bdx * cdy - bdy * cdx)) + + ((bdx * bdx + bdy * bdy) * + ((cdx * adytail + ady * cdxtail) - (cdy * adxtail + adx * cdytail)) + + 2.0 * (bdx * bdxtail + bdy * bdytail) * (cdx * ady - cdy * adx)) + + ((cdx * cdx + cdy * cdy) * + ((adx * bdytail + bdy * adxtail) - (ady * bdxtail + bdx * adytail)) + + 2.0 * (cdx * cdxtail + cdy * cdytail) * (adx * bdy - ady * bdx)); + if ((det >= errbound) || (-det >= errbound)) { + return det; + } + + finnow = fin1; + finother = fin2; + + if ((bdxtail != 0.0) || (bdytail != 0.0) || (cdxtail != 0.0) || (cdytail != 0.0)) { + Square(adx, adxadx1, adxadx0); + Square(ady, adyady1, adyady0); + Two_Two_Sum(adxadx1, adxadx0, adyady1, adyady0, aa3, aa[2], aa[1], aa[0]); + aa[3] = aa3; + } + if ((cdxtail != 0.0) || (cdytail != 0.0) || (adxtail != 0.0) || (adytail != 0.0)) { + Square(bdx, bdxbdx1, bdxbdx0); + Square(bdy, bdybdy1, bdybdy0); + Two_Two_Sum(bdxbdx1, bdxbdx0, bdybdy1, bdybdy0, bb3, bb[2], bb[1], bb[0]); + bb[3] = bb3; + } + if ((adxtail != 0.0) || (adytail != 0.0) || (bdxtail != 0.0) || (bdytail != 0.0)) { + Square(cdx, cdxcdx1, cdxcdx0); + Square(cdy, cdycdy1, cdycdy0); + Two_Two_Sum(cdxcdx1, cdxcdx0, cdycdy1, cdycdy0, cc3, cc[2], cc[1], cc[0]); + cc[3] = cc3; + } + + if (adxtail != 0.0) { + axtbclen = scale_expansion_zeroelim(4, bc, adxtail, axtbc); + temp16alen = scale_expansion_zeroelim(axtbclen, axtbc, 2.0 * adx, temp16a); + + axtcclen = scale_expansion_zeroelim(4, cc, adxtail, axtcc); + temp16blen = scale_expansion_zeroelim(axtcclen, axtcc, bdy, temp16b); + + axtbblen = scale_expansion_zeroelim(4, bb, adxtail, axtbb); + temp16clen = scale_expansion_zeroelim(axtbblen, axtbb, -cdy, temp16c); + + temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a, temp16blen, temp16b, temp32a); + temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c, temp32alen, temp32a, temp48); + finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len, temp48, finother); + finswap = finnow; + finnow = finother; + finother = finswap; + } + if (adytail != 0.0) { + aytbclen = scale_expansion_zeroelim(4, bc, adytail, aytbc); + temp16alen = scale_expansion_zeroelim(aytbclen, aytbc, 2.0 * ady, temp16a); + + aytbblen = scale_expansion_zeroelim(4, bb, adytail, aytbb); + temp16blen = scale_expansion_zeroelim(aytbblen, aytbb, cdx, temp16b); + + aytcclen = scale_expansion_zeroelim(4, cc, adytail, aytcc); + temp16clen = scale_expansion_zeroelim(aytcclen, aytcc, -bdx, temp16c); + + temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a, temp16blen, temp16b, temp32a); + temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c, temp32alen, temp32a, temp48); + finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len, temp48, finother); + finswap = finnow; + finnow = finother; + finother = finswap; + } + if (bdxtail != 0.0) { + bxtcalen = scale_expansion_zeroelim(4, ca, bdxtail, bxtca); + temp16alen = scale_expansion_zeroelim(bxtcalen, bxtca, 2.0 * bdx, temp16a); + + bxtaalen = scale_expansion_zeroelim(4, aa, bdxtail, bxtaa); + temp16blen = scale_expansion_zeroelim(bxtaalen, bxtaa, cdy, temp16b); + + bxtcclen = scale_expansion_zeroelim(4, cc, bdxtail, bxtcc); + temp16clen = scale_expansion_zeroelim(bxtcclen, bxtcc, -ady, temp16c); + + temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a, temp16blen, temp16b, temp32a); + temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c, temp32alen, temp32a, temp48); + finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len, temp48, finother); + finswap = finnow; + finnow = finother; + finother = finswap; + } + if (bdytail != 0.0) { + bytcalen = scale_expansion_zeroelim(4, ca, bdytail, bytca); + temp16alen = scale_expansion_zeroelim(bytcalen, bytca, 2.0 * bdy, temp16a); + + bytcclen = scale_expansion_zeroelim(4, cc, bdytail, bytcc); + temp16blen = scale_expansion_zeroelim(bytcclen, bytcc, adx, temp16b); + + bytaalen = scale_expansion_zeroelim(4, aa, bdytail, bytaa); + temp16clen = scale_expansion_zeroelim(bytaalen, bytaa, -cdx, temp16c); + + temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a, temp16blen, temp16b, temp32a); + temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c, temp32alen, temp32a, temp48); + finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len, temp48, finother); + finswap = finnow; + finnow = finother; + finother = finswap; + } + if (cdxtail != 0.0) { + cxtablen = scale_expansion_zeroelim(4, ab, cdxtail, cxtab); + temp16alen = scale_expansion_zeroelim(cxtablen, cxtab, 2.0 * cdx, temp16a); + + cxtbblen = scale_expansion_zeroelim(4, bb, cdxtail, cxtbb); + temp16blen = scale_expansion_zeroelim(cxtbblen, cxtbb, ady, temp16b); + + cxtaalen = scale_expansion_zeroelim(4, aa, cdxtail, cxtaa); + temp16clen = scale_expansion_zeroelim(cxtaalen, cxtaa, -bdy, temp16c); + + temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a, temp16blen, temp16b, temp32a); + temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c, temp32alen, temp32a, temp48); + finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len, temp48, finother); + finswap = finnow; + finnow = finother; + finother = finswap; + } + if (cdytail != 0.0) { + cytablen = scale_expansion_zeroelim(4, ab, cdytail, cytab); + temp16alen = scale_expansion_zeroelim(cytablen, cytab, 2.0 * cdy, temp16a); + + cytaalen = scale_expansion_zeroelim(4, aa, cdytail, cytaa); + temp16blen = scale_expansion_zeroelim(cytaalen, cytaa, bdx, temp16b); + + cytbblen = scale_expansion_zeroelim(4, bb, cdytail, cytbb); + temp16clen = scale_expansion_zeroelim(cytbblen, cytbb, -adx, temp16c); + + temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a, temp16blen, temp16b, temp32a); + temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c, temp32alen, temp32a, temp48); + finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len, temp48, finother); + finswap = finnow; + finnow = finother; + finother = finswap; + } + + if ((adxtail != 0.0) || (adytail != 0.0)) { + if ((bdxtail != 0.0) || (bdytail != 0.0) || (cdxtail != 0.0) || (cdytail != 0.0)) { + Two_Product(bdxtail, cdy, ti1, ti0); + Two_Product(bdx, cdytail, tj1, tj0); + Two_Two_Sum(ti1, ti0, tj1, tj0, u3, u[2], u[1], u[0]); + u[3] = u3; + negate = -bdy; + Two_Product(cdxtail, negate, ti1, ti0); + negate = -bdytail; + Two_Product(cdx, negate, tj1, tj0); + Two_Two_Sum(ti1, ti0, tj1, tj0, v3, v[2], v[1], v[0]); + v[3] = v3; + bctlen = fast_expansion_sum_zeroelim(4, u, 4, v, bct); + + Two_Product(bdxtail, cdytail, ti1, ti0); + Two_Product(cdxtail, bdytail, tj1, tj0); + Two_Two_Diff(ti1, ti0, tj1, tj0, bctt3, bctt[2], bctt[1], bctt[0]); + bctt[3] = bctt3; + bcttlen = 4; + } + else { + bct[0] = 0.0; + bctlen = 1; + bctt[0] = 0.0; + bcttlen = 1; + } + + if (adxtail != 0.0) { + temp16alen = scale_expansion_zeroelim(axtbclen, axtbc, adxtail, temp16a); + axtbctlen = scale_expansion_zeroelim(bctlen, bct, adxtail, axtbct); + temp32alen = scale_expansion_zeroelim(axtbctlen, axtbct, 2.0 * adx, temp32a); + temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a, temp32alen, temp32a, temp48); + finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len, temp48, finother); + finswap = finnow; + finnow = finother; + finother = finswap; + if (bdytail != 0.0) { + temp8len = scale_expansion_zeroelim(4, cc, adxtail, temp8); + temp16alen = scale_expansion_zeroelim(temp8len, temp8, bdytail, temp16a); + finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen, temp16a, finother); + finswap = finnow; + finnow = finother; + finother = finswap; + } + if (cdytail != 0.0) { + temp8len = scale_expansion_zeroelim(4, bb, -adxtail, temp8); + temp16alen = scale_expansion_zeroelim(temp8len, temp8, cdytail, temp16a); + finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen, temp16a, finother); + finswap = finnow; + finnow = finother; + finother = finswap; + } + + temp32alen = scale_expansion_zeroelim(axtbctlen, axtbct, adxtail, temp32a); + axtbcttlen = scale_expansion_zeroelim(bcttlen, bctt, adxtail, axtbctt); + temp16alen = scale_expansion_zeroelim(axtbcttlen, axtbctt, 2.0 * adx, temp16a); + temp16blen = scale_expansion_zeroelim(axtbcttlen, axtbctt, adxtail, temp16b); + temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a, temp16blen, temp16b, temp32b); + temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a, temp32blen, temp32b, temp64); + finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len, temp64, finother); + finswap = finnow; + finnow = finother; + finother = finswap; + } + if (adytail != 0.0) { + temp16alen = scale_expansion_zeroelim(aytbclen, aytbc, adytail, temp16a); + aytbctlen = scale_expansion_zeroelim(bctlen, bct, adytail, aytbct); + temp32alen = scale_expansion_zeroelim(aytbctlen, aytbct, 2.0 * ady, temp32a); + temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a, temp32alen, temp32a, temp48); + finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len, temp48, finother); + finswap = finnow; + finnow = finother; + finother = finswap; + + temp32alen = scale_expansion_zeroelim(aytbctlen, aytbct, adytail, temp32a); + aytbcttlen = scale_expansion_zeroelim(bcttlen, bctt, adytail, aytbctt); + temp16alen = scale_expansion_zeroelim(aytbcttlen, aytbctt, 2.0 * ady, temp16a); + temp16blen = scale_expansion_zeroelim(aytbcttlen, aytbctt, adytail, temp16b); + temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a, temp16blen, temp16b, temp32b); + temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a, temp32blen, temp32b, temp64); + finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len, temp64, finother); + finswap = finnow; + finnow = finother; + finother = finswap; + } + } + if ((bdxtail != 0.0) || (bdytail != 0.0)) { + if ((cdxtail != 0.0) || (cdytail != 0.0) || (adxtail != 0.0) || (adytail != 0.0)) { + Two_Product(cdxtail, ady, ti1, ti0); + Two_Product(cdx, adytail, tj1, tj0); + Two_Two_Sum(ti1, ti0, tj1, tj0, u3, u[2], u[1], u[0]); + u[3] = u3; + negate = -cdy; + Two_Product(adxtail, negate, ti1, ti0); + negate = -cdytail; + Two_Product(adx, negate, tj1, tj0); + Two_Two_Sum(ti1, ti0, tj1, tj0, v3, v[2], v[1], v[0]); + v[3] = v3; + catlen = fast_expansion_sum_zeroelim(4, u, 4, v, cat); + + Two_Product(cdxtail, adytail, ti1, ti0); + Two_Product(adxtail, cdytail, tj1, tj0); + Two_Two_Diff(ti1, ti0, tj1, tj0, catt3, catt[2], catt[1], catt[0]); + catt[3] = catt3; + cattlen = 4; + } + else { + cat[0] = 0.0; + catlen = 1; + catt[0] = 0.0; + cattlen = 1; + } + + if (bdxtail != 0.0) { + temp16alen = scale_expansion_zeroelim(bxtcalen, bxtca, bdxtail, temp16a); + bxtcatlen = scale_expansion_zeroelim(catlen, cat, bdxtail, bxtcat); + temp32alen = scale_expansion_zeroelim(bxtcatlen, bxtcat, 2.0 * bdx, temp32a); + temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a, temp32alen, temp32a, temp48); + finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len, temp48, finother); + finswap = finnow; + finnow = finother; + finother = finswap; + if (cdytail != 0.0) { + temp8len = scale_expansion_zeroelim(4, aa, bdxtail, temp8); + temp16alen = scale_expansion_zeroelim(temp8len, temp8, cdytail, temp16a); + finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen, temp16a, finother); + finswap = finnow; + finnow = finother; + finother = finswap; + } + if (adytail != 0.0) { + temp8len = scale_expansion_zeroelim(4, cc, -bdxtail, temp8); + temp16alen = scale_expansion_zeroelim(temp8len, temp8, adytail, temp16a); + finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen, temp16a, finother); + finswap = finnow; + finnow = finother; + finother = finswap; + } + + temp32alen = scale_expansion_zeroelim(bxtcatlen, bxtcat, bdxtail, temp32a); + bxtcattlen = scale_expansion_zeroelim(cattlen, catt, bdxtail, bxtcatt); + temp16alen = scale_expansion_zeroelim(bxtcattlen, bxtcatt, 2.0 * bdx, temp16a); + temp16blen = scale_expansion_zeroelim(bxtcattlen, bxtcatt, bdxtail, temp16b); + temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a, temp16blen, temp16b, temp32b); + temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a, temp32blen, temp32b, temp64); + finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len, temp64, finother); + finswap = finnow; + finnow = finother; + finother = finswap; + } + if (bdytail != 0.0) { + temp16alen = scale_expansion_zeroelim(bytcalen, bytca, bdytail, temp16a); + bytcatlen = scale_expansion_zeroelim(catlen, cat, bdytail, bytcat); + temp32alen = scale_expansion_zeroelim(bytcatlen, bytcat, 2.0 * bdy, temp32a); + temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a, temp32alen, temp32a, temp48); + finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len, temp48, finother); + finswap = finnow; + finnow = finother; + finother = finswap; + + temp32alen = scale_expansion_zeroelim(bytcatlen, bytcat, bdytail, temp32a); + bytcattlen = scale_expansion_zeroelim(cattlen, catt, bdytail, bytcatt); + temp16alen = scale_expansion_zeroelim(bytcattlen, bytcatt, 2.0 * bdy, temp16a); + temp16blen = scale_expansion_zeroelim(bytcattlen, bytcatt, bdytail, temp16b); + temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a, temp16blen, temp16b, temp32b); + temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a, temp32blen, temp32b, temp64); + finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len, temp64, finother); + finswap = finnow; + finnow = finother; + finother = finswap; + } + } + if ((cdxtail != 0.0) || (cdytail != 0.0)) { + if ((adxtail != 0.0) || (adytail != 0.0) || (bdxtail != 0.0) || (bdytail != 0.0)) { + Two_Product(adxtail, bdy, ti1, ti0); + Two_Product(adx, bdytail, tj1, tj0); + Two_Two_Sum(ti1, ti0, tj1, tj0, u3, u[2], u[1], u[0]); + u[3] = u3; + negate = -ady; + Two_Product(bdxtail, negate, ti1, ti0); + negate = -adytail; + Two_Product(bdx, negate, tj1, tj0); + Two_Two_Sum(ti1, ti0, tj1, tj0, v3, v[2], v[1], v[0]); + v[3] = v3; + abtlen = fast_expansion_sum_zeroelim(4, u, 4, v, abt); + + Two_Product(adxtail, bdytail, ti1, ti0); + Two_Product(bdxtail, adytail, tj1, tj0); + Two_Two_Diff(ti1, ti0, tj1, tj0, abtt3, abtt[2], abtt[1], abtt[0]); + abtt[3] = abtt3; + abttlen = 4; + } + else { + abt[0] = 0.0; + abtlen = 1; + abtt[0] = 0.0; + abttlen = 1; + } + + if (cdxtail != 0.0) { + temp16alen = scale_expansion_zeroelim(cxtablen, cxtab, cdxtail, temp16a); + cxtabtlen = scale_expansion_zeroelim(abtlen, abt, cdxtail, cxtabt); + temp32alen = scale_expansion_zeroelim(cxtabtlen, cxtabt, 2.0 * cdx, temp32a); + temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a, temp32alen, temp32a, temp48); + finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len, temp48, finother); + finswap = finnow; + finnow = finother; + finother = finswap; + if (adytail != 0.0) { + temp8len = scale_expansion_zeroelim(4, bb, cdxtail, temp8); + temp16alen = scale_expansion_zeroelim(temp8len, temp8, adytail, temp16a); + finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen, temp16a, finother); + finswap = finnow; + finnow = finother; + finother = finswap; + } + if (bdytail != 0.0) { + temp8len = scale_expansion_zeroelim(4, aa, -cdxtail, temp8); + temp16alen = scale_expansion_zeroelim(temp8len, temp8, bdytail, temp16a); + finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen, temp16a, finother); + finswap = finnow; + finnow = finother; + finother = finswap; + } + + temp32alen = scale_expansion_zeroelim(cxtabtlen, cxtabt, cdxtail, temp32a); + cxtabttlen = scale_expansion_zeroelim(abttlen, abtt, cdxtail, cxtabtt); + temp16alen = scale_expansion_zeroelim(cxtabttlen, cxtabtt, 2.0 * cdx, temp16a); + temp16blen = scale_expansion_zeroelim(cxtabttlen, cxtabtt, cdxtail, temp16b); + temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a, temp16blen, temp16b, temp32b); + temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a, temp32blen, temp32b, temp64); + finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len, temp64, finother); + finswap = finnow; + finnow = finother; + finother = finswap; + } + if (cdytail != 0.0) { + temp16alen = scale_expansion_zeroelim(cytablen, cytab, cdytail, temp16a); + cytabtlen = scale_expansion_zeroelim(abtlen, abt, cdytail, cytabt); + temp32alen = scale_expansion_zeroelim(cytabtlen, cytabt, 2.0 * cdy, temp32a); + temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a, temp32alen, temp32a, temp48); + finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len, temp48, finother); + finswap = finnow; + finnow = finother; + finother = finswap; + + temp32alen = scale_expansion_zeroelim(cytabtlen, cytabt, cdytail, temp32a); + cytabttlen = scale_expansion_zeroelim(abttlen, abtt, cdytail, cytabtt); + temp16alen = scale_expansion_zeroelim(cytabttlen, cytabtt, 2.0 * cdy, temp16a); + temp16blen = scale_expansion_zeroelim(cytabttlen, cytabtt, cdytail, temp16b); + temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a, temp16blen, temp16b, temp32b); + temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a, temp32blen, temp32b, temp64); + finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len, temp64, finother); + finswap = finnow; + finnow = finother; + finother = finswap; + } + } + + return finnow[finlength - 1]; +} + +static double incircle(const double *pa, const double *pb, const double *pc, const double *pd) +{ + double adx, bdx, cdx, ady, bdy, cdy; + double bdxcdy, cdxbdy, cdxady, adxcdy, adxbdy, bdxady; + double alift, blift, clift; + double det; + double permanent, errbound; + + adx = pa[0] - pd[0]; + bdx = pb[0] - pd[0]; + cdx = pc[0] - pd[0]; + ady = pa[1] - pd[1]; + bdy = pb[1] - pd[1]; + cdy = pc[1] - pd[1]; + + bdxcdy = bdx * cdy; + cdxbdy = cdx * bdy; + alift = adx * adx + ady * ady; + + cdxady = cdx * ady; + adxcdy = adx * cdy; + blift = bdx * bdx + bdy * bdy; + + adxbdy = adx * bdy; + bdxady = bdx * ady; + clift = cdx * cdx + cdy * cdy; + + det = alift * (bdxcdy - cdxbdy) + blift * (cdxady - adxcdy) + clift * (adxbdy - bdxady); + + permanent = (Absolute(bdxcdy) + Absolute(cdxbdy)) * alift + + (Absolute(cdxady) + Absolute(adxcdy)) * blift + + (Absolute(adxbdy) + Absolute(bdxady)) * clift; + errbound = iccerrboundA * permanent; + if ((det > errbound) || (-det > errbound)) { + return det; + } + + return incircleadapt(pa, pb, pc, pd, permanent); +} diff --git a/tests/gtests/blenlib/BLI_delaunay_2d_test.cc b/tests/gtests/blenlib/BLI_delaunay_2d_test.cc index 946dc03eca3..15cb44f4431 100644 --- a/tests/gtests/blenlib/BLI_delaunay_2d_test.cc +++ b/tests/gtests/blenlib/BLI_delaunay_2d_test.cc @@ -17,8 +17,8 @@ extern "C" { #include <sstream> #define DO_REGULAR_TESTS 1 -#define DO_RANDOM_TESTS 0 -#define DO_FILE_TESTS 0 +#define DO_RANDOM_TESTS 1 +#define DO_FILE_TESTS 1 static void fill_input_verts(CDT_input *r_input, float (*vcos)[2], int nverts) { @@ -322,15 +322,20 @@ TEST(delaunay, OnePt) { CDT_input in; CDT_result *out; - float p[][2] = {{0.0f, 0.0f}}; + const char *spec = R"(1 0 0 + 0.0 0.0 + )"; - fill_input_verts(&in, p, 1); + fill_input_from_string(&in, spec); out = BLI_delaunay_2d_cdt_calc(&in, CDT_FULL); EXPECT_EQ(out->verts_len, 1); EXPECT_EQ(out->edges_len, 0); EXPECT_EQ(out->faces_len, 0); - EXPECT_EQ(out->vert_coords[0][0], 0.0f); - EXPECT_EQ(out->vert_coords[0][1], 0.0f); + if (out->verts_len >= 1) { + EXPECT_EQ(out->vert_coords[0][0], 0.0f); + EXPECT_EQ(out->vert_coords[0][1], 0.0f); + } + free_spec_arrays(&in); BLI_delaunay_2d_cdt_free(out); } @@ -339,9 +344,12 @@ TEST(delaunay, TwoPt) CDT_input in; CDT_result *out; int v0_out, v1_out, e0_out; - float p[][2] = {{0.0f, -0.75f}, {0.0f, 0.75f}}; + const char *spec = R"(2 0 0 + 0.0 -0.75 + 0.0 0.75 + )"; - fill_input_verts(&in, p, 2); + fill_input_from_string(&in, spec); out = BLI_delaunay_2d_cdt_calc(&in, CDT_FULL); EXPECT_EQ(out->verts_len, 2); EXPECT_EQ(out->edges_len, 1); @@ -351,12 +359,15 @@ TEST(delaunay, TwoPt) EXPECT_NE(v0_out, -1); EXPECT_NE(v1_out, -1); EXPECT_NE(v0_out, v1_out); - EXPECT_NEAR(out->vert_coords[v0_out][0], p[0][0], in.epsilon); - EXPECT_NEAR(out->vert_coords[v0_out][1], p[0][1], in.epsilon); - EXPECT_NEAR(out->vert_coords[v1_out][0], p[1][0], in.epsilon); - EXPECT_NEAR(out->vert_coords[v1_out][1], p[1][1], in.epsilon); + if (out->verts_len >= 2) { + EXPECT_NEAR(out->vert_coords[v0_out][0], 0.0, in.epsilon); + EXPECT_NEAR(out->vert_coords[v0_out][1], -0.75, in.epsilon); + EXPECT_NEAR(out->vert_coords[v1_out][0], 0.0, in.epsilon); + EXPECT_NEAR(out->vert_coords[v1_out][1], 0.75, in.epsilon); + } e0_out = get_edge(out, v0_out, v1_out); EXPECT_EQ(e0_out, 0); + free_spec_arrays(&in); BLI_delaunay_2d_cdt_free(out); } @@ -367,9 +378,13 @@ TEST(delaunay, ThreePt) int v0_out, v1_out, v2_out; int e0_out, e1_out, e2_out; int f0_out; - float p[][2] = {{-0.1f, -0.75f}, {0.1f, 0.75f}, {0.5f, 0.5f}}; + const char *spec = R"(3 0 0 + -0.1 -0.75 + 0.1 0.75 + 0.5 0.5 + )"; - fill_input_verts(&in, p, 3); + fill_input_from_string(&in, spec); out = BLI_delaunay_2d_cdt_calc(&in, CDT_FULL); EXPECT_EQ(out->verts_len, 3); EXPECT_EQ(out->edges_len, 3); @@ -386,6 +401,7 @@ TEST(delaunay, ThreePt) EXPECT_TRUE(e0_out != e1_out && e0_out != e2_out && e1_out != e2_out); f0_out = get_face_tri(out, v0_out, v2_out, v1_out); EXPECT_EQ(f0_out, 0); + free_spec_arrays(&in); BLI_delaunay_2d_cdt_free(out); } @@ -394,11 +410,14 @@ TEST(delaunay, ThreePtsMerge) CDT_input in; CDT_result *out; int v0_out, v1_out, v2_out; - /* equilateral triangle with side 0.1 */ - float p[][2] = {{-0.05f, -0.05f}, {0.05f, -0.05f}, {0.0f, 0.03660254f}}; + const char *spec = R"(3 0 0 + -0.05 -0.05 + 0.05 -0.05 + 0.0 0.03660254 + )"; /* First with epsilon such that points are within that distance of each other */ - fill_input_verts(&in, p, 3); + fill_input_from_string(&in, spec); in.epsilon = 0.21f; out = BLI_delaunay_2d_cdt_calc(&in, CDT_FULL); EXPECT_EQ(out->verts_len, 1); @@ -422,6 +441,7 @@ TEST(delaunay, ThreePtsMerge) EXPECT_EQ(out->verts_len, 3); EXPECT_EQ(out->edges_len, 3); EXPECT_EQ(out->faces_len, 1); + free_spec_arrays(&in); BLI_delaunay_2d_cdt_free(out); } @@ -429,13 +449,19 @@ TEST(delaunay, MixedPts) { CDT_input in; CDT_result *out; - float p[][2] = {{0.0f, 0.0f}, {-0.5f, -0.5f}, {-0.4f, -0.25f}, {-0.3f, 0.8}}; - int e[][2] = {{0, 1}, {1, 2}, {2, 3}}; int v0_out, v1_out, v2_out, v3_out; int e0_out, e1_out, e2_out; + const char *spec = R"(4 3 0 + 0.0 0.0 + -0.5 -0.5 + -0.4 -0.25 + -0.3 0.8 + 0 1 + 1 2 + 2 3 + )"; - fill_input_verts(&in, p, 4); - add_input_edges(&in, e, 3); + fill_input_from_string(&in, spec); out = BLI_delaunay_2d_cdt_calc(&in, CDT_FULL); EXPECT_EQ(out->verts_len, 4); EXPECT_EQ(out->edges_len, 6); @@ -451,6 +477,134 @@ TEST(delaunay, MixedPts) EXPECT_TRUE(out_edge_has_input_id(out, e0_out, 0)); EXPECT_TRUE(out_edge_has_input_id(out, e1_out, 1)); EXPECT_TRUE(out_edge_has_input_id(out, e2_out, 2)); + free_spec_arrays(&in); + BLI_delaunay_2d_cdt_free(out); +} + +TEST(delaunay, Quad0) +{ + CDT_input in; + CDT_result *out; + int e_diag_out; + const char *spec = R"(4 0 0 + 0.0 1.0 + 1,0. 0.0 + 2.0 0.1 + 2.25 0.5 + )"; + fill_input_from_string(&in, spec); + out = BLI_delaunay_2d_cdt_calc(&in, CDT_FULL); + EXPECT_EQ(out->verts_len, 4); + EXPECT_EQ(out->edges_len, 5); + e_diag_out = get_edge(out, 1, 3); + EXPECT_NE(e_diag_out, -1); + free_spec_arrays(&in); + BLI_delaunay_2d_cdt_free(out); +} + +TEST(delaunay, Quad1) +{ + CDT_input in; + CDT_result *out; + int e_diag_out; + const char *spec = R"(4 0 0 + 0.0 0.0 + 0.9 -1.0 + 2.0 0.0 + 0.9 3.0 + )"; + fill_input_from_string(&in, spec); + out = BLI_delaunay_2d_cdt_calc(&in, CDT_FULL); + EXPECT_EQ(out->verts_len, 4); + EXPECT_EQ(out->edges_len, 5); + e_diag_out = get_edge(out, 0, 2); + EXPECT_NE(e_diag_out, -1); + free_spec_arrays(&in); + BLI_delaunay_2d_cdt_free(out); +} + +TEST(delaunay, Quad2) +{ + CDT_input in; + CDT_result *out; + int e_diag_out; + const char *spec = R"(4 0 0 + 0.5 0.0 + 0.15 0.2 + 0.3 0.4 + .45 0.35 + )"; + fill_input_from_string(&in, spec); + out = BLI_delaunay_2d_cdt_calc(&in, CDT_FULL); + EXPECT_EQ(out->verts_len, 4); + EXPECT_EQ(out->edges_len, 5); + e_diag_out = get_edge(out, 1, 3); + EXPECT_NE(e_diag_out, -1); + free_spec_arrays(&in); + BLI_delaunay_2d_cdt_free(out); +} + +TEST(delaunay, Quad3) +{ + CDT_input in; + CDT_result *out; + int e_diag_out; + const char *spec = R"(4 0 0 + 0.5 0.0 + 0.0 0.0 + 0.3 0.4 + .45 0.35 + )"; + fill_input_from_string(&in, spec); + out = BLI_delaunay_2d_cdt_calc(&in, CDT_FULL); + EXPECT_EQ(out->verts_len, 4); + EXPECT_EQ(out->edges_len, 5); + e_diag_out = get_edge(out, 0, 2); + EXPECT_NE(e_diag_out, -1); + free_spec_arrays(&in); + BLI_delaunay_2d_cdt_free(out); +} + +TEST(delaunay, Quad4) +{ + CDT_input in; + CDT_result *out; + int e_diag_out; + const char *spec = R"(4 0 0 + 1.0 1.0 + 0.0 0.0 + 1.0 -3.0 + 0.0 1.0 + )"; + fill_input_from_string(&in, spec); + out = BLI_delaunay_2d_cdt_calc(&in, CDT_FULL); + EXPECT_EQ(out->verts_len, 4); + EXPECT_EQ(out->edges_len, 5); + e_diag_out = get_edge(out, 0, 1); + EXPECT_NE(e_diag_out, -1); + free_spec_arrays(&in); + BLI_delaunay_2d_cdt_free(out); +} + +TEST(delaunay, LineInSquare) +{ + CDT_input in; + CDT_result *out; + const char *spec = R"(6 1 1 + -0.5 -0.5 + 0.5 -0.5 + -0.5 0.5 + 0.5 0.5 + -0.25 0.0 + 0.25 0.0 + 4 5 + 0 1 3 2 + )"; + fill_input_from_string(&in, spec); + out = BLI_delaunay_2d_cdt_calc(&in, CDT_CONSTRAINTS); + EXPECT_EQ(out->verts_len, 6); + EXPECT_EQ(out->faces_len, 1); + free_spec_arrays(&in); BLI_delaunay_2d_cdt_free(out); } @@ -458,13 +612,18 @@ TEST(delaunay, CrossSegs) { CDT_input in; CDT_result *out; - float p[][2] = {{-0.5f, 0.0f}, {0.5f, 0.0f}, {-0.4f, -0.5f}, {0.4f, 0.5f}}; - int e[][2] = {{0, 1}, {2, 3}}; int v0_out, v1_out, v2_out, v3_out, v_intersect; int i; + const char *spec = R"(4 2 0 + -0.5 0.0 + 0.5 0.0 + -0.4 -0.5 + 0.4 0.5 + 0 1 + 2 3 + )"; - fill_input_verts(&in, p, 4); - add_input_edges(&in, e, 2); + fill_input_from_string(&in, spec); out = BLI_delaunay_2d_cdt_calc(&in, CDT_FULL); EXPECT_EQ(out->verts_len, 5); EXPECT_EQ(out->edges_len, 8); @@ -482,8 +641,11 @@ TEST(delaunay, CrossSegs) } } EXPECT_NE(v_intersect, -1); - EXPECT_NEAR(out->vert_coords[v_intersect][0], 0.0f, in.epsilon); - EXPECT_NEAR(out->vert_coords[v_intersect][1], 0.0f, in.epsilon); + if (v_intersect != -1) { + EXPECT_NEAR(out->vert_coords[v_intersect][0], 0.0f, in.epsilon); + EXPECT_NEAR(out->vert_coords[v_intersect][1], 0.0f, in.epsilon); + } + free_spec_arrays(&in); BLI_delaunay_2d_cdt_free(out); } @@ -491,23 +653,27 @@ TEST(delaunay, DiamondCross) { CDT_input in; CDT_result *out; - float p[][2] = { - {0.0f, 0.0f}, - {1.0f, 3.0f}, - {2.0f, 0.0f}, - {1.0f, -3.0f}, - {0.0f, 0.0f}, - {1.0f, -3.0f}, - {1.0f, 3.0f}, - }; - int e[][2] = {{0, 1}, {1, 2}, {2, 3}, {3, 4}, {5, 6}}; - - fill_input_verts(&in, p, 7); - add_input_edges(&in, e, 5); + const char *spec = R"(7 5 0 + 0.0 0.0 + 1.0 3.0 + 2.0 0.0 + 1.0 -3.0 + 0.0 0.0 + 1.0 -3.0 + 1.0 3.0 + 0 1 + 1 2 + 2 3 + 3 4 + 5 6 + )"; + + fill_input_from_string(&in, spec); out = BLI_delaunay_2d_cdt_calc(&in, CDT_FULL); EXPECT_EQ(out->verts_len, 4); EXPECT_EQ(out->edges_len, 5); EXPECT_EQ(out->faces_len, 2); + free_spec_arrays(&in); BLI_delaunay_2d_cdt_free(out); } @@ -516,27 +682,35 @@ TEST(delaunay, TwoDiamondsCrossed) CDT_input in; CDT_result *out; /* Input has some repetition of vertices, on purpose */ - float p[][2] = { - {0.0f, 0.0f}, - {1.0f, 2.0f}, - {2.0f, 0.0f}, - {1.0f, -2.0f}, - {0.0f, 0.0f}, - {3.0f, 0.0f}, - {4.0f, 2.0f}, - {5.0f, 0.0f}, - {4.0f, -2.0f}, - {3.0f, 0.0f}, - {0.0f, 0.0f}, - {5.0f, 0.0f}, - }; int e[][2] = {{0, 1}, {1, 2}, {2, 3}, {3, 4}, {5, 6}, {6, 7}, {7, 8}, {8, 9}, {10, 11}}; int v_out[12]; int e_out[9], e_cross_1, e_cross_2, e_cross_3; int i; + const char *spec = R"(12 9 0 + 0.0 0.0 + 1.0 2.0 + 2.0 0.0 + 1.0 -2.0 + 0.0 0.0 + 3.0 0.0 + 4.0 2.0 + 5.0 0.0 + 4.0 -2.0 + 3.0 0.0 + 0.0 0.0 + 5.0 0.0 + 0 1 + 1 2 + 2 3 + 3 4 + 5 6 + 6 7 + 7 8 + 8 9 + 10 11 + )"; - fill_input_verts(&in, p, 12); - add_input_edges(&in, e, 9); + fill_input_from_string(&in, spec); out = BLI_delaunay_2d_cdt_calc(&in, CDT_FULL); EXPECT_EQ(out->verts_len, 8); EXPECT_EQ(out->edges_len, 15); @@ -562,6 +736,7 @@ TEST(delaunay, TwoDiamondsCrossed) EXPECT_TRUE(out_edge_has_input_id(out, e_cross_1, 8)); EXPECT_TRUE(out_edge_has_input_id(out, e_cross_2, 8)); EXPECT_TRUE(out_edge_has_input_id(out, e_cross_3, 8)); + free_spec_arrays(&in); BLI_delaunay_2d_cdt_free(out); } @@ -570,53 +745,63 @@ TEST(delaunay, ManyCross) CDT_input in; CDT_result *out; /* Input has some repetition of vertices, on purpose */ - float p[][2] = { - /* upper: verts 0 to 10 */ - {0.0f, 0.0f}, - {6.0f, 9.0f}, - {15.0f, 18.0f}, - {35.0f, 13.0f}, - {43.0f, 18.0f}, - {57.0f, 12.0f}, - {69.0f, 10.0f}, - {78.0f, 0.0f}, - {91.0f, 0.0f}, - {107.0f, 22.0f}, - {123.0f, 0.0f}, - /* lower part 1: verts 11 to 16 */ - {0.0f, 0.0f}, - {10.0f, -14.0f}, - {35.0f, -8.0f}, - {43.0f, -12.0f}, - {64.0f, -13.0f}, - {78.0f, 0.0f}, - /* lower part 2: verts 17 to 20 */ - {91.0f, 0.0f}, - {102.0f, -9.0f}, - {116.0f, -9.0f}, - {123.0f, 0.0f}, - /* cross 1: verts 21, 22 */ - {43.0f, 18.0f}, - {43.0f, -12.0f}, - /* cross 2: verts 23, 24 */ - {107.0f, 22.0f}, - {102.0f, -9.0f}, - /* cross all: verts 25, 26 */ - {0.0f, 0.0f}, - {123.0f, 0.0f}, - }; - int e[][2] = { - {0, 1}, {1, 2}, {2, 3}, {3, 4}, {4, 5}, {5, 6}, {6, 7}, - {7, 8}, {8, 9}, {9, 10}, {11, 12}, {12, 13}, {13, 14}, {14, 15}, - {15, 16}, {17, 18}, {18, 19}, {19, 20}, {21, 22}, {23, 24}, {25, 26}, - }; - - fill_input_verts(&in, p, 27); - add_input_edges(&in, e, 21); + const char *spec = R"(27 21 0 + 0.0 0.0 + 6.0 9.0 + 15.0 18.0 + 35.0 13.0 + 43.0 18.0 + 57.0 12.0 + 69.0 10.0 + 78.0 0.0 + 91.0 0.0 + 107.0 22.0 + 123.0 0.0 + 0.0 0.0 + 10.0 -14.0 + 35.0 -8.0 + 43.0 -12.0 + 64.0 -13.0 + 78.0 0.0 + 91.0 0.0 + 102.0 -9.0 + 116.0 -9.0 + 123.0 0.0 + 43.0 18.0 + 43.0 -12.0 + 107.0 22.0 + 102.0 -9.0 + 0.0 0.0 + 123.0 0.0 + 0 1 + 1 2 + 2 3 + 3 4 + 4 5 + 5 6 + 6 7 + 7 8 + 8 9 + 9 10 + 11 12 + 12 13 + 13 14 + 14 15 + 15 16 + 17 18 + 18 19 + 19 20 + 21 22 + 23 24 + 25 26 + )"; + + fill_input_from_string(&in, spec); out = BLI_delaunay_2d_cdt_calc(&in, CDT_FULL); EXPECT_EQ(out->verts_len, 19); EXPECT_EQ(out->edges_len, 46); EXPECT_EQ(out->faces_len, 28); + free_spec_arrays(&in); BLI_delaunay_2d_cdt_free(out); } @@ -624,16 +809,20 @@ TEST(delaunay, TwoFace) { CDT_input in; CDT_result *out; - float p[][2] = { - {0.0f, 0.0f}, {1.0f, 0.0f}, {0.5f, 1.0f}, {1.1f, 1.0f}, {1.1f, 0.0f}, {1.6f, 1.0f}}; - int f[] = {/* 0 */ 0, 1, 2, /* 1 */ 3, 4, 5}; - int fstart[] = {0, 3}; - int flen[] = {3, 3}; int v_out[6], f0_out, f1_out, e0_out, e1_out, e2_out; int i; + const char *spec = R"(6 0 2 + 0.0 0.0 + 1.0 0.0 + 0.5 1.0 + 1.1 1.0 + 1.1 0.0 + 1.6 1.0 + 0 1 2 + 3 4 5 + )"; - fill_input_verts(&in, p, 6); - add_input_faces(&in, f, fstart, flen, 2); + fill_input_from_string(&in, spec); out = BLI_delaunay_2d_cdt_calc(&in, CDT_FULL); EXPECT_EQ(out->verts_len, 6); EXPECT_EQ(out->edges_len, 9); @@ -657,6 +846,7 @@ TEST(delaunay, TwoFace) EXPECT_TRUE(out_edge_has_input_id(out, e2_out, out->face_edge_offset + 2)); EXPECT_TRUE(out_face_has_input_id(out, f0_out, 0)); EXPECT_TRUE(out_face_has_input_id(out, f1_out, 1)); + free_spec_arrays(&in); BLI_delaunay_2d_cdt_free(out); } @@ -664,28 +854,27 @@ TEST(delaunay, OverlapFaces) { CDT_input in; CDT_result *out; - float p[][2] = { - {0.0f, 0.0f}, - {1.0f, 0.0f}, - {1.0f, 1.0f}, - {0.0f, 1.0f}, - {0.5f, 0.5f}, - {1.5f, 0.5f}, - {1.5f, 1.3f}, - {0.5f, 1.3f}, - {0.1f, 0.1f}, - {0.3f, 0.1f}, - {0.3f, 0.3f}, - {0.1f, 0.3f}, - }; - int f[] = {/* 0 */ 0, 1, 2, 3, /* 1 */ 4, 5, 6, 7, /* 2*/ 8, 9, 10, 11}; - int fstart[] = {0, 4, 8}; - int flen[] = {4, 4, 4}; int v_out[12], v_int1, v_int2, f0_out, f1_out, f2_out; int i; + const char *spec = R"(12 0 3 + 0.0 0.0 + 1.0 0.0 + 1.0 1.0 + 0.0 1.0 + 0.5 0.5 + 1.5 0.5 + 1.5 1.3 + 0.5 1.3 + 0.1 0.1 + 0.3 0.1 + 0.3 0.3 + 0.1 0.3 + 0 1 2 3 + 4 5 6 7 + 8 9 10 11 + )"; - fill_input_verts(&in, p, 12); - add_input_faces(&in, f, fstart, flen, 3); + fill_input_from_string(&in, spec); out = BLI_delaunay_2d_cdt_calc(&in, CDT_FULL); EXPECT_EQ(out->verts_len, 14); EXPECT_EQ(out->edges_len, 33); @@ -696,22 +885,26 @@ TEST(delaunay, OverlapFaces) } v_int1 = 12; v_int2 = 13; - if (fabsf(out->vert_coords[v_int1][0] - 1.0f) > in.epsilon) { - v_int1 = 13; - v_int2 = 12; + if (out->verts_len > 13) { + if (fabsf(out->vert_coords[v_int1][0] - 1.0f) > in.epsilon) { + v_int1 = 13; + v_int2 = 12; + } + EXPECT_NEAR(out->vert_coords[v_int1][0], 1.0, in.epsilon); + EXPECT_NEAR(out->vert_coords[v_int1][1], 0.5, in.epsilon); + EXPECT_NEAR(out->vert_coords[v_int2][0], 0.5, in.epsilon); + EXPECT_NEAR(out->vert_coords[v_int2][1], 1.0, in.epsilon); } - EXPECT_NEAR(out->vert_coords[v_int1][0], 1.0, in.epsilon); - EXPECT_NEAR(out->vert_coords[v_int1][1], 0.5, in.epsilon); - EXPECT_NEAR(out->vert_coords[v_int2][0], 0.5, in.epsilon); - EXPECT_NEAR(out->vert_coords[v_int2][1], 1.0, in.epsilon); f0_out = get_face_tri(out, v_out[1], v_int1, v_out[4]); EXPECT_NE(f0_out, -1); EXPECT_TRUE(out_face_has_input_id(out, f0_out, 0)); f1_out = get_face_tri(out, v_out[4], v_int1, v_out[2]); EXPECT_NE(f1_out, -1); EXPECT_TRUE(out_face_has_input_id(out, f1_out, 0)); - EXPECT_TRUE(out_face_has_input_id(out, f1_out, 0)); + EXPECT_TRUE(out_face_has_input_id(out, f1_out, 1)); f2_out = get_face_tri(out, v_out[8], v_out[9], v_out[10]); + if (f2_out == -1) + f2_out = get_face_tri(out, v_out[8], v_out[9], v_out[11]); EXPECT_NE(f2_out, -1); EXPECT_TRUE(out_face_has_input_id(out, f2_out, 0)); EXPECT_TRUE(out_face_has_input_id(out, f2_out, 2)); @@ -728,6 +921,7 @@ TEST(delaunay, OverlapFaces) out = BLI_delaunay_2d_cdt_calc(&in, CDT_CONSTRAINTS_VALID_BMESH); EXPECT_EQ(out->faces_len, 5); + free_spec_arrays(&in); BLI_delaunay_2d_cdt_free(out); } @@ -735,52 +929,25 @@ TEST(delaunay, TwoSquaresOverlap) { CDT_input in; CDT_result *out; - float p[][2] = { - {1.0f, -1.0f}, - {-1.0f, -1.0f}, - {-1.0f, 1.0f}, - {1.0f, 1.0f}, - {-1.5f, 1.5f}, - {0.5f, 1.5f}, - {0.5f, -0.5f}, - {-1.5f, -0.5f}, - }; - int f[] = {/* 0 */ 7, 6, 5, 4, /* 1 */ 3, 2, 1, 0}; - int fstart[] = {0, 4}; - int flen[] = {4, 4}; - - fill_input_verts(&in, p, 8); - add_input_faces(&in, f, fstart, flen, 2); + const char *spec = R"(8 0 2 + 1.0 -1.0 + -1.0 -1.0 + -1.0 1.0 + 1.0 1.0 + -1.5 1.5 + 0.5 1.5 + 0.5 -0.5 + -1.5 -0.5 + 7 6 5 4 + 3 2 1 0 + )"; + + fill_input_from_string(&in, spec); out = BLI_delaunay_2d_cdt_calc(&in, CDT_CONSTRAINTS_VALID_BMESH); EXPECT_EQ(out->verts_len, 10); EXPECT_EQ(out->edges_len, 12); EXPECT_EQ(out->faces_len, 3); - BLI_delaunay_2d_cdt_free(out); -} - -TEST(delaunay, TriCutoff) -{ - CDT_input in; - CDT_result *out; - float p[][2] = { - {-3.53009f, 1.29403f}, - {-4.11844f, -1.08375f}, - {1.56893f, 1.29403f}, - {0.621034f, 0.897734f}, - {0.549125f, 1.29403f}, - }; - int f[] = {0, 2, 1}; - int fstart[] = {0}; - int flen[] = {3}; - int e[][2] = {{3, 4}}; - - fill_input_verts(&in, p, 5); - add_input_faces(&in, f, fstart, flen, 1); - add_input_edges(&in, e, 1); - out = BLI_delaunay_2d_cdt_calc(&in, CDT_CONSTRAINTS_VALID_BMESH); - EXPECT_EQ(out->verts_len, 5); - EXPECT_EQ(out->edges_len, 6); - EXPECT_EQ(out->faces_len, 2); + free_spec_arrays(&in); BLI_delaunay_2d_cdt_free(out); } @@ -788,24 +955,23 @@ TEST(delaunay, TriInTri) { CDT_input in; CDT_result *out; - float p[][2] = { - {-5.65685f, 0.0f}, - {1.41421f, -5.83095f}, - {0.0f, 0.0f}, - {-2.47487f, -1.45774f}, - {-0.707107f, -2.91548f}, - {-1.06066f, -1.45774f}, - }; - int f[] = {0, 1, 2, 3, 4, 5}; - int fstart[] = {0, 3}; - int flen[] = {3, 3}; - - fill_input_verts(&in, p, 6); - add_input_faces(&in, f, fstart, flen, 2); + const char *spec = R"(6 0 2 + -5.65685 0.0 + 1.41421 -5.83095 + 0.0 0.0 + -2.47487 -1.45774 + -0.707107 -2.91548 + -1.06066 -1.45774 + 0 1 2 + 3 4 5 + )"; + + fill_input_from_string(&in, spec); out = BLI_delaunay_2d_cdt_calc(&in, CDT_CONSTRAINTS_VALID_BMESH); EXPECT_EQ(out->verts_len, 6); EXPECT_EQ(out->edges_len, 8); EXPECT_EQ(out->faces_len, 3); + free_spec_arrays(&in); BLI_delaunay_2d_cdt_free(out); } @@ -831,6 +997,7 @@ TEST(delaunay, DiamondInSquare) EXPECT_EQ(out->edges_len, 10); EXPECT_EQ(out->faces_len, 3); free_spec_arrays(&in); + BLI_delaunay_2d_cdt_free(out); } TEST(delaunay, DiamondInSquareWire) @@ -861,52 +1028,76 @@ TEST(delaunay, DiamondInSquareWire) EXPECT_EQ(out->edges_len, 8); EXPECT_EQ(out->faces_len, 2); free_spec_arrays(&in); + BLI_delaunay_2d_cdt_free(out); } -TEST(delaunay, ClosePts) +TEST(delaunay, TinyEdge) { CDT_input in; CDT_result *out; - const char *spec = R"(7 2 1 - 0.46876350045204163 0.06087132915854454 - 0.46865847706794739 0.03632887825369835 - 0.49176687002182007 0.03632888197898865 - 0.49166208505630493 0.06087132543325424 - 0.49171400070190430 0.04841339960694313 - 0.49171534180641174 0.04839951172471046 - 0.49045535922050476 0.06087132915854454 - 4 5 - 6 4 - 0 1 2 3 + /* An intersect with triangle would be at (0.8, 0.2). */ + const char *spec = R"(4 1 1 + 0.0 0.0 + 1.0 0.0 + 0.5 0.5 + 0.84 0.21 + 0 3 + 0 1 2 )"; fill_input_from_string(&in, spec); - out = BLI_delaunay_2d_cdt_calc(&in, CDT_FULL); - EXPECT_EQ(out->verts_len, 7); - EXPECT_EQ(out->edges_len, 12); - EXPECT_EQ(out->faces_len, 6); + in.epsilon = 0.1; + out = BLI_delaunay_2d_cdt_calc(&in, CDT_CONSTRAINTS); + EXPECT_EQ(out->verts_len, 4); + EXPECT_EQ(out->edges_len, 5); + EXPECT_EQ(out->faces_len, 2); free_spec_arrays(&in); + BLI_delaunay_2d_cdt_free(out); } -TEST(delaunay, ClosePts2) +TEST(delaunay, TinyEdge2) { CDT_input in; CDT_result *out; + /* An intersect with triangle would be at (0.8, 0.2). */ const char *spec = R"(6 1 1 - -0.17878936231136322 -0.44374340772628784 - -0.17871695756912231 -0.45601493120193481 - -0.17544384300708771 -0.45601493120193481 - -0.17537136375904083 -0.44374340772628784 - -0.17544738948345184 -0.45602506399154663 - -0.17872454226016998 -0.45472940802574158 - 4 5 - 0 1 2 3 + 0.0 0.0 + 0.2 -0.2 + 1.0 0.0 + 0.5 0.5 + 0.2 0.4 + 0.84 0.21 + 0 5 + 0 1 2 3 4 )"; fill_input_from_string(&in, spec); - out = BLI_delaunay_2d_cdt_calc(&in, CDT_FULL); + in.epsilon = 0.1; + out = BLI_delaunay_2d_cdt_calc(&in, CDT_CONSTRAINTS); EXPECT_EQ(out->verts_len, 6); - EXPECT_EQ(out->edges_len, 10); - EXPECT_EQ(out->faces_len, 5); + EXPECT_EQ(out->edges_len, 7); + EXPECT_EQ(out->faces_len, 2); free_spec_arrays(&in); + BLI_delaunay_2d_cdt_free(out); +} + +TEST(delaunay, repeatededge) +{ + CDT_input in; + CDT_result *out; + const char *spec = R"(5 3 0 + 0.0 0.0 + 0.0 1.0 + 1.0 1.1 + 0.5 -0.5 + 0.5 2.5 + 0 1 + 2 3 + 2 3 + )"; + fill_input_from_string(&in, spec); + out = BLI_delaunay_2d_cdt_calc(&in, CDT_CONSTRAINTS); + EXPECT_EQ(out->edges_len, 2); + free_spec_arrays(&in); + BLI_delaunay_2d_cdt_free(out); } #endif @@ -915,67 +1106,216 @@ enum { RANDOM_PTS, RANDOM_SEGS, RANDOM_POLY, + RANDOM_TILTED_GRID, + RANDOM_CIRCLE, + RANDOM_TRI_BETWEEN_CIRCLES, }; -// #define DO_TIMING -static void rand_delaunay_test(int test_kind, - int max_lg_size, - int reps_per_size, - CDT_output_type otype) +# define DO_TIMING +static void rand_delaunay_test( + int test_kind, int max_lg_size, int reps_per_size, double param, CDT_output_type otype) { CDT_input in; CDT_result *out; - int lg_size, size, rep, i, npts, nedges; + int lg_size, size, rep, i, j, size_max, npts_max, nedges_max, nfaces_max, npts, nedges, nfaces; + int ia, ib, ic; float(*p)[2]; int(*e)[2]; + int *faces, *faces_start_table, *faces_len_table; + double start_angle, angle_delta, angle1, angle2, angle3; + float orient; double tstart; double *times; RNG *rng; rng = BLI_rng_new(0); - npts = (1 << max_lg_size); - p = (float(*)[2])MEM_malloc_arrayN(npts, 2 * sizeof(float), "delaunay"); + e = NULL; + faces = NULL; + faces_start_table = NULL; + faces_len_table = NULL; + nedges_max = 0; + nfaces_max = 0; + + /* Set up npts, nedges, nfaces, and allocate needed arrays at max length needed. */ + size_max = 1 << max_lg_size; switch (test_kind) { case RANDOM_PTS: - nedges = 0; - e = NULL; - break; - case RANDOM_SEGS: case RANDOM_POLY: - /* TODO: use faces for poly case, but need to deal with winding parity issue */ - nedges = npts - 1 + (test_kind == RANDOM_POLY); - e = (int(*)[2])MEM_malloc_arrayN(nedges, 2 * sizeof(int), "delaunay"); + npts_max = size_max; + if (test_kind == RANDOM_SEGS) { + nedges_max = npts_max - 1; + } + else if (test_kind == RANDOM_POLY) { + nedges_max = npts_max; + } + break; + + case RANDOM_TILTED_GRID: + /* A 'size' x 'size' grid of points, tilted by angle 'param'. + * Edges will go from left ends to right ends and tops to bottoms, so 2 x size of them. + * Depending on epsilon, the vertical-ish edges may or may not go through the intermediate + * vertices, but the horizontal ones always should. + */ + npts_max = size_max * size_max; + nedges_max = 2 * size_max; + break; + + case RANDOM_CIRCLE: + /* A circle with 'size' points, a random start angle, and equal spacing thereafter. + * Will be input as one face. + */ + npts_max = size_max; + nfaces_max = 1; + break; + + case RANDOM_TRI_BETWEEN_CIRCLES: + /* A set of 'size' triangles, each has two random points on the unit circle, + * and the third point is a random point on the circle with radius 'param'. + * Each triangle will be input as a face. + */ + npts_max = 3 * size_max; + nfaces_max = size_max; break; default: fprintf(stderr, "unknown random delaunay test kind\n"); return; } - times = (double *)MEM_malloc_arrayN(max_lg_size + 1, sizeof(double), "delaunay"); + fprintf(stderr, + "size_max=%d, npts_max=%d, nedges_max=%d, nfaces_max=%d\n", + size_max, + npts_max, + nedges_max, + nfaces_max); /*DEBUG!!*/ + p = (float(*)[2])MEM_malloc_arrayN(npts_max, 2 * sizeof(float), __func__); + if (nedges_max > 0) { + e = (int(*)[2])MEM_malloc_arrayN(nedges_max, 2 * sizeof(int), __func__); + } + if (nfaces_max > 0) { + faces_start_table = (int *)MEM_malloc_arrayN(nfaces_max, sizeof(int), __func__); + faces_len_table = (int *)MEM_malloc_arrayN(nfaces_max, sizeof(int), __func__); + faces = (int *)MEM_malloc_arrayN(npts_max, sizeof(int), __func__); + } + + times = (double *)MEM_malloc_arrayN(max_lg_size + 1, sizeof(double), __func__); + + /* For powers of 2 sizes up to max_lg_size power of 2. */ for (lg_size = 0; lg_size <= max_lg_size; lg_size++) { size = 1 << lg_size; + nedges = 0; + nfaces = 0; times[lg_size] = 0.0; - if (size == 1 && test_kind != RANDOM_PTS) + if (size == 1 && test_kind != RANDOM_PTS) { continue; + } + /* Do 'rep' repetitions. */ for (rep = 0; rep < reps_per_size; rep++) { - for (i = 0; i < size; i++) { - p[i][0] = (float)BLI_rng_get_double(rng); /* will be in range in [0,1) */ - p[i][1] = (float)BLI_rng_get_double(rng); + /* Make vertices and edges or faces. */ + switch (test_kind) { + case RANDOM_PTS: + case RANDOM_SEGS: + case RANDOM_POLY: + npts = size; + if (test_kind == RANDOM_SEGS) { + nedges = npts - 1; + } + else if (test_kind == RANDOM_POLY) { + nedges = npts; + } + for (i = 0; i < size; i++) { + p[i][0] = (float)BLI_rng_get_double(rng); /* will be in range in [0,1) */ + p[i][1] = (float)BLI_rng_get_double(rng); + if (test_kind != RANDOM_PTS) { + if (i > 0) { + e[i - 1][0] = i - 1; + e[i - 1][1] = i; + } + } + } + if (test_kind == RANDOM_POLY) { + e[size - 1][0] = size - 1; + e[size - 1][1] = 0; + } + break; + + case RANDOM_TILTED_GRID: + /* 'param' is slope of tilt of vertical lines. */ + npts = size * size; + nedges = 2 * size; + for (i = 0; i < size; i++) { + for (j = 0; j < size; j++) { + p[i * size + j][0] = i * param + j; + p[i * size + j][1] = i; + } + } + for (i = 0; i < size; i++) { + /* Horizontal edges: connect p(i,0) to p(i,size-1). */ + e[i][0] = i * size; + e[i][1] = i * size + size - 1; + /* Vertical edges: conntect p(0,i) to p(size-1,i). */ + e[size + i][0] = i; + e[size + i][1] = (size - 1) * size + i; + } + break; + + case RANDOM_CIRCLE: + npts = size; + nfaces = 1; + faces_start_table[0] = 0; + faces_len_table[0] = npts; + start_angle = BLI_rng_get_double(rng) * 2.0 * M_PI; + angle_delta = 2.0 * M_PI / size; + for (i = 0; i < size; i++) { + p[i][0] = (float)cos(start_angle + i * angle_delta); + p[i][1] = (float)sin(start_angle + i * angle_delta); + faces[i] = i; + } + break; + + case RANDOM_TRI_BETWEEN_CIRCLES: + npts = 3 * size; + nfaces = size; + for (i = 0; i < size; i++) { + /* Get three random angles in [0, 2pi). */ + angle1 = BLI_rng_get_double(rng) * 2.0 * M_PI; + angle2 = BLI_rng_get_double(rng) * 2.0 * M_PI; + angle3 = BLI_rng_get_double(rng) * 2.0 * M_PI; + ia = 3 * i; + ib = 3 * i + 1; + ic = 3 * i + 2; + p[ia][0] = (float)cos(angle1); + p[ia][1] = (float)sin(angle1); + p[ib][0] = (float)cos(angle2); + p[ib][1] = (float)sin(angle2); + p[ic][0] = (float)(param * cos(angle3)); + p[ic][1] = (float)(param * sin(angle3)); + faces_start_table[i] = 3 * i; + faces_len_table[i] = 3; + /* Put the coordinates in ccw order. */ + faces[ia] = ia; + orient = (p[ia][0] - p[ic][0]) * (p[ib][1] - p[ic][1]) - + (p[ib][0] - p[ic][0]) * (p[ia][1] - p[ic][1]); + if (orient >= 0.0f) { + faces[ib] = ib; + faces[ic] = ic; + } + else { + faces[ib] = ic; + faces[ic] = ib; + } + } + break; } - fill_input_verts(&in, p, size); - - if (test_kind == RANDOM_SEGS || test_kind == RANDOM_POLY) { - for (i = 0; i < size - 1; i++) { - e[i][0] = i; - e[i][1] = i + 1; - } - if (test_kind == RANDOM_POLY) { - e[size - 1][0] = size - 1; - e[size - 1][1] = 0; - } - add_input_edges(&in, e, size - 1 + (test_kind == RANDOM_POLY)); + fill_input_verts(&in, p, npts); + if (nedges > 0) { + add_input_edges(&in, e, nedges); + } + if (nfaces > 0) { + add_input_faces(&in, faces, faces_start_table, faces_len_table, nfaces); } + + /* Run the test. */ tstart = PIL_check_seconds_timer(); out = BLI_delaunay_2d_cdt_calc(&in, otype); EXPECT_NE(out->verts_len, 0); @@ -990,46 +1330,82 @@ static void rand_delaunay_test(int test_kind, } # endif MEM_freeN(p); - if (e) + if (e) { MEM_freeN(e); + } + if (faces) { + MEM_freeN(faces); + MEM_freeN(faces_start_table); + MEM_freeN(faces_len_table); + } MEM_freeN(times); BLI_rng_free(rng); } TEST(delaunay, randompts) { - rand_delaunay_test(RANDOM_PTS, 7, 100, CDT_FULL); + rand_delaunay_test(RANDOM_PTS, 7, 1, 0.0, CDT_FULL); } TEST(delaunay, randomsegs) { - rand_delaunay_test(RANDOM_SEGS, 7, 100, CDT_FULL); + rand_delaunay_test(RANDOM_SEGS, 7, 1, 0.0, CDT_FULL); } TEST(delaunay, randompoly) { - rand_delaunay_test(RANDOM_POLY, 7, 100, CDT_FULL); + rand_delaunay_test(RANDOM_POLY, 7, 1, 0.0, CDT_FULL); } TEST(delaunay, randompoly_inside) { - rand_delaunay_test(RANDOM_POLY, 7, 1, CDT_INSIDE); + rand_delaunay_test(RANDOM_POLY, 7, 1, 0.0, CDT_INSIDE); } TEST(delaunay, randompoly_constraints) { - rand_delaunay_test(RANDOM_POLY, 7, 1, CDT_CONSTRAINTS); + rand_delaunay_test(RANDOM_POLY, 7, 1, 0.0, CDT_CONSTRAINTS); } TEST(delaunay, randompoly_validbmesh) { - rand_delaunay_test(RANDOM_POLY, 7, 1, CDT_CONSTRAINTS_VALID_BMESH); + rand_delaunay_test(RANDOM_POLY, 7, 1, 0.0, CDT_CONSTRAINTS_VALID_BMESH); +} + +TEST(delaunay, grid) +{ + rand_delaunay_test(RANDOM_TILTED_GRID, 6, 1, 0.0, CDT_FULL); +} + +TEST(delaunay, tilted_grid_a) +{ + rand_delaunay_test(RANDOM_TILTED_GRID, 6, 1, 1.0, CDT_FULL); +} + +TEST(delaunay, tilted_grid_b) +{ + rand_delaunay_test(RANDOM_TILTED_GRID, 6, 1, 0.01, CDT_FULL); +} + +TEST(delaunay, randomcircle) +{ + rand_delaunay_test(RANDOM_CIRCLE, 7, 1, 0.0, CDT_FULL); +} + +TEST(delaunay, random_tris_circle) +{ + rand_delaunay_test(RANDOM_TRI_BETWEEN_CIRCLES, 6, 1, 0.25, CDT_FULL); +} + +TEST(delaunay, random_tris_circle_b) +{ + rand_delaunay_test(RANDOM_TRI_BETWEEN_CIRCLES, 6, 1, 1e-4, CDT_FULL); } #endif #if DO_FILE_TESTS /* For timing large examples of points only. - * The given file should have one point per line, as a space-separated pair of floats. + * See fill_input_from_file for file format. */ static void points_from_file_test(const char *filename) { @@ -1045,17 +1421,19 @@ static void points_from_file_test(const char *filename) free_spec_arrays(&in); } +# if 0 TEST(delaunay, debug) { CDT_input in; CDT_result *out; fill_input_from_file(&in, "/tmp/cdtinput.txt"); - out = BLI_delaunay_2d_cdt_calc(&in, CDT_FULL); + out = BLI_delaunay_2d_cdt_calc(&in, CDT_CONSTRAINTS); BLI_delaunay_2d_cdt_free(out); free_spec_arrays(&in); } +# endif -# if 0 +# if 1 # define POINTFILEROOT "/tmp/" TEST(delaunay, terrain1) |