From: Barak Itkin Date: Fri, 27 Apr 2012 17:04:49 +0000 (+0300) Subject: Commit the code for the Delaunay Terminator Algorithm X-Git-Tag: p2tc-0.1.0~90 X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=598b0d57a906f26bfeb223cecc3a82acc829fb82;p=poly2tri-c Commit the code for the Delaunay Terminator Algorithm --- diff --git a/refine/delaunay-terminator.c b/refine/delaunay-terminator.c index a5c5853..0e93029 100644 --- a/refine/delaunay-terminator.c +++ b/refine/delaunay-terminator.c @@ -1,3 +1,5 @@ +#include +#include #include #include "utils.h" @@ -9,9 +11,21 @@ #include "mesh.h" #include "cdt.h" +#include "cluster.h" #include "delaunay-terminator.h" +/* The code in this file is based on the "Delaunay Terminator" algorithm + * - an algorithm for refining constrained delaunay triangulations. The + * algorithm itself appears in a paper by Jonathan Richard Shewchuk as + * described below: + * + * Delaunay Refinement Algorithms for Triangular Mesh Generation + * Computational Geometry: Theory and Applications 22(1–3):21–74, May 2002 + * Jonathan Richard Shewchuk + * http://www.cs.berkeley.edu/~jrs/papers/2dj.pdf + */ + gboolean p2tr_cdt_test_encroachment_ignore_visibility (const P2trVector2 *w, P2trEdge *e) @@ -66,4 +80,346 @@ p2tr_cdt_is_encroached (P2trEdge *E) return (T1 != NULL && p2tr_cdt_test_encroachment_ignore_visibility (&p2tr_triangle_get_opposite_point (T1, E)->c, E)) || (T2 != NULL && p2tr_cdt_test_encroachment_ignore_visibility (&p2tr_triangle_get_opposite_point (T2, E)->c, E)); -} \ No newline at end of file +} + +/* ****************************************************************** */ +/* Now for the algorithm itself */ +/* ****************************************************************** */ + +static gboolean +SplitPermitted (P2trDelaunayTerminator *self, P2trEdge *s, gdouble d); + +static void +SplitEncroachedSubsegments (P2trDelaunayTerminator *self, gdouble theta, P2trTriangleTooBig delta); + +static void +NewVertex (P2trDelaunayTerminator *self, P2trPoint *v, gdouble theta, P2trTriangleTooBig delta); + +static gdouble +ShortestEdgeLength (P2trTriangle *tri); + +static gboolean +TolerantIsShorter (P2trEdge *toTest, P2trEdge *reference); + +static inline gdouble +LOG2 (gdouble v); + +static gboolean +TolerantIsPowerOfTwoLength (gdouble length); + +static void +ChooseSplitVertex(P2trEdge *e, P2trVector2 *dst); + + + +static inline gint +triangle_quality_compare (P2trTriangle *t1, P2trTriangle *t2) +{ + gdouble a1 = p2tr_triangle_smallest_non_constrained_angle (t1); + gdouble a2 = p2tr_triangle_smallest_non_constrained_angle (t2); + + return (a1 < a2) ? -1 : ((a1 == a2) ? 0 : 1); +} + +P2trDelaunayTerminator* +p2tr_dt_new (gdouble theta, P2trTriangleTooBig delta, P2trCDT *cdt) +{ + P2trDelaunayTerminator *self = g_slice_new (P2trDelaunayTerminator); + self->Qt = g_sequence_new (NULL); + g_queue_init (&self->Qs); + self->delta = delta; + self->theta = theta; + self->mesh = cdt; + return self; +} + +static void +p2tr_dt_enqueue_tri (P2trDelaunayTerminator *self, + P2trTriangle *tri) +{ + p2tr_triangle_ref (tri); + g_sequence_insert_sorted (self->Qt, tri, (GCompareDataFunc)triangle_quality_compare, NULL); +} + +static inline gboolean +p2tr_dt_tri_queue_is_empty (P2trDelaunayTerminator *self) +{ + return g_sequence_iter_is_end (g_sequence_get_begin_iter (self->Qt)); +} + +static P2trTriangle* +p2tr_dt_dequeue_tri (P2trDelaunayTerminator *self) +{ + GSequenceIter *first = g_sequence_get_begin_iter (self->Qt); + + /* If we have an empty sequence, return NULL */ + if (p2tr_dt_tri_queue_is_empty (self)) + return NULL; + else + { + P2trTriangle *ret = g_sequence_get (first); + g_sequence_remove (first); + return ret; + } +} + +static void +p2tr_dt_enqueue_segment (P2trDelaunayTerminator *self, + P2trEdge *E) +{ + if (! E->constrained) + p2tr_exception_programmatic ("Tried to append a non-segment!"); + + p2tr_edge_ref (E); + g_queue_push_tail (&self->Qs, E); +} + +static P2trEdge* +p2tr_dt_dequeue_segment (P2trDelaunayTerminator *self) +{ + if (g_queue_is_empty (&self->Qs)) + return NULL; + else + return g_queue_pop_head (&self->Qs); +} + +static gboolean +p2tr_dt_segment_queue_is_empty (P2trDelaunayTerminator *self) +{ + return g_queue_is_empty (&self->Qs); +} + +static gboolean +false_too_big (P2trTriangle *tri) +{ + return FALSE; +} + +void p2tr_dt_refine (P2trDelaunayTerminator *self, + gint max_steps) +{ + P2trHashSetIter hs_iter; + P2trEdge *s; + P2trTriangle *t; + gint steps; + + p2tr_hash_set_iter_init (&hs_iter, self->mesh->mesh->edges); + while (p2tr_hash_set_iter_next (&hs_iter, (gpointer*)&s)) + if (s->constrained && p2tr_cdt_is_encroached (s)) + p2tr_dt_enqueue_segment (self, s); + + SplitEncroachedSubsegments (self, 0, false_too_big); + + p2tr_hash_set_iter_init (&hs_iter, self->mesh->mesh->triangles); + while (p2tr_hash_set_iter_next (&hs_iter, (gpointer*)&t)) + if (p2tr_triangle_smallest_non_constrained_angle (t) < self->theta) + p2tr_dt_enqueue_tri (self, t); + + steps = 0; + while (! p2tr_dt_tri_queue_is_empty (self) && steps < max_steps) + { + t = p2tr_dt_dequeue_tri (self); + + if (p2tr_hash_set_contains (self->mesh->mesh->triangles, t)) + { + P2trCircle tCircum; + P2trVector2 *c; + P2trTriangle *triContaining_c_NOREF; + P2trHashSet *E; + + p2tr_triangle_get_circum_circle (t, &tCircum); + c = &tCircum.center; + + triContaining_c_NOREF = p2tr_mesh_find_point2 (self->mesh->mesh, c, t); + + /* If no edge is encroached, then this must be + * inside the triangulation domain!!! */ + if (triContaining_c_NOREF == NULL) + p2tr_exception_geometric ("Should not happen!"); + + /* Now, check if this point would encroach any edge + * of the triangulation */ + E = p2tr_cdt_get_segments_encroached_by (self->mesh, c); + + if (p2tr_hash_set_size (E) == 0) + { + P2trPoint *cPoint; + cPoint = p2tr_cdt_insert_point (self->mesh, c, triContaining_c_NOREF); + NewVertex (self, cPoint, self->theta, self->delta); + + p2tr_point_unref (cPoint); + } + else + { + gdouble d = ShortestEdgeLength(t); + + p2tr_hash_set_iter_init (&hs_iter, E); + while (p2tr_hash_set_iter_next (&hs_iter, (gpointer*)&s)) + if (self->delta (t) || SplitPermitted(self, s, d)) + p2tr_dt_enqueue_segment (self, s); + + if (! p2tr_dt_segment_queue_is_empty (self)) + { + p2tr_dt_enqueue_tri (self, t); + SplitEncroachedSubsegments(self, self->theta, self->delta); + } + } + p2tr_hash_set_free (E); + } + + p2tr_triangle_unref (t); + } +} + +static gboolean +SplitPermitted (P2trDelaunayTerminator *self, P2trEdge *s, gdouble d) +{ + P2trCluster *startCluster = p2tr_cluster_get_for (P2TR_EDGE_START(s), s); + P2trCluster *endCluster = p2tr_cluster_get_for (s->end, s); + P2trCluster *S_NOREF; + GList *iter; + + gboolean permitted = FALSE; + + if (! TolerantIsPowerOfTwoLength (p2tr_edge_get_length (s)) + /* True when different, meaning both null or both exist */ + || ((startCluster != NULL) ^ (endCluster == NULL))) + { + permitted = TRUE; + } + + if (! permitted) + { + S_NOREF = (startCluster != NULL) ? startCluster : endCluster; + + for (iter = g_queue_peek_head_link (&S_NOREF->edges); iter != NULL; iter = iter->next) + if (TolerantIsShorter((P2trEdge*) iter->data, s)) /* e shorter than s */ + { + permitted = TRUE; + break; + } + } + + if (! permitted) + { + gdouble rmin = p2tr_edge_get_length(s) * sin(S_NOREF->min_angle / 2); + if (rmin >= d) + permitted = TRUE; + } + + if (startCluster) p2tr_cluster_free (startCluster); + if (endCluster) p2tr_cluster_free (endCluster); + + return permitted; +} + +static void +SplitEncroachedSubsegments (P2trDelaunayTerminator *self, gdouble theta, P2trTriangleTooBig delta) +{ + while (! p2tr_dt_segment_queue_is_empty (self)) + { + P2trEdge *s = p2tr_dt_dequeue_segment (self); + if (p2tr_hash_set_contains (self->mesh->mesh->edges, s)) + { + P2trVector2 v; + P2trPoint *Pv; + GList *parts, *iter; + + ChooseSplitVertex (s, &v); + Pv = p2tr_mesh_new_point (self->mesh->mesh, &v); + + /* Update here if using diametral lenses */ + + parts = p2tr_cdt_split_edge (self->mesh, s, Pv); + + NewVertex (self, Pv, theta, delta); + + for (iter = parts; iter != NULL; iter = iter->next) + { + P2trEdge *e = (P2trEdge*)iter->data; + if (p2tr_cdt_is_encroached (e)) + p2tr_dt_enqueue_segment (self, e); + p2tr_edge_unref (e); + } + } + p2tr_edge_unref (s); + } +} + +static void +NewVertex (P2trDelaunayTerminator *self, P2trPoint *v, gdouble theta, P2trTriangleTooBig delta) +{ + GList *iter; + for (iter = v->outgoing_edges; iter != NULL; iter = iter->next) + { + P2trEdge *outEdge = (P2trEdge*) iter->data; + P2trTriangle *t = outEdge->tri; + P2trEdge *e; + + if (t == NULL) + continue; + + e = p2tr_triangle_get_opposite_edge (t, v); + + /* we want the fast check and for new points we don't + * use that check... So let's go on the full check + * since it's still faster */ + if (e->constrained && p2tr_cdt_is_encroached (e)) + p2tr_dt_enqueue_segment (self, e); + else if (delta (t) || p2tr_triangle_smallest_non_constrained_angle (t) < theta) + p2tr_dt_enqueue_tri (self, t); + } +} + +static gdouble +ShortestEdgeLength (P2trTriangle *tri) +{ + gdouble a1 = p2tr_edge_get_length_squared (tri->edges[0]); + gdouble a2 = p2tr_edge_get_length_squared (tri->edges[1]); + gdouble a3 = p2tr_edge_get_length_squared (tri->edges[2]); + + return sqrt (MIN (a1, MIN (a2, a3))); +} + +static gboolean +TolerantIsShorter (P2trEdge *toTest, P2trEdge *reference) +{ + return p2tr_edge_get_length(toTest) < p2tr_edge_get_length(reference) * 1.01; +} + +static inline gdouble +LOG2 (gdouble v) +{ + return log10 (v) / G_LOG_2_BASE_10; +} + +static gboolean +TolerantIsPowerOfTwoLength (gdouble length) +{ + gdouble exp = LOG2 (length); + gdouble intpart, frac = modf (exp, &intpart); + + return ABS(frac) < 0.05; +} + +static void +ChooseSplitVertex(P2trEdge *e, P2trVector2 *dst) +{ + gdouble sourceLength = p2tr_edge_get_length(e); + gdouble newLength = pow(2, round(LOG2(sourceLength))); + gdouble ratio, resultLength; + + /* IMPORTANT! DIVIDE BY 2! */ + newLength /= 2; + + ratio = newLength / sourceLength; + + dst->x = (1 - ratio) * P2TR_EDGE_START(e)->c.x + (ratio) * e->end->c.x; + dst->y = (1 - ratio) * P2TR_EDGE_START(e)->c.y + (ratio) * e->end->c.y; + + /* now let's avoid consistency problems */ + resultLength = sqrt(P2TR_VECTOR2_DISTANCE_SQ(&P2TR_EDGE_START(e)->c, dst)); + + if (! TolerantIsPowerOfTwoLength(resultLength)) + p2tr_exception_numeric ("Bad rounding!"); +} diff --git a/refine/delaunay-terminator.h b/refine/delaunay-terminator.h index 3a9b96d..e536bd6 100644 --- a/refine/delaunay-terminator.h +++ b/refine/delaunay-terminator.h @@ -4,10 +4,15 @@ #include #include "cdt.h" +typedef gboolean (*P2trTriangleTooBig) (P2trTriangle *tri); + typedef struct { - P2trCDT *mesh; - GQueue Qs, Qt; + P2trCDT *mesh; + GQueue Qs; + GSequence *Qt; + gdouble theta; + P2trTriangleTooBig delta; } P2trDelaunayTerminator; gboolean p2tr_cdt_test_encroachment_ignore_visibility (const P2trVector2 *w, @@ -22,5 +27,7 @@ P2trHashSet* p2tr_cdt_get_segments_encroached_by (P2trCDT *self, gboolean p2tr_cdt_is_encroached (P2trEdge *E); +P2trDelaunayTerminator* +p2tr_dt_new (gdouble theta, P2trTriangleTooBig delta, P2trCDT *cdt); #endif \ No newline at end of file diff --git a/refine/utils.h b/refine/utils.h index 93be8b1..efbf967 100755 --- a/refine/utils.h +++ b/refine/utils.h @@ -57,6 +57,7 @@ extern "C" #define p2tr_hash_set_remove(set,element) g_hash_table_remove ((set), (element)) #define p2tr_hash_set_remove_all(set) g_hash_table_remove_all ((set)) #define p2tr_hash_set_free(set) g_hash_table_destroy(set) +#define p2tr_hash_set_size(set) g_hash_table_size(set) #define p2tr_hash_set_iter_init(iter,hash_set) g_hash_table_iter_init ((iter),(hash_set)) #define p2tr_hash_set_iter_next(iter,val) g_hash_table_iter_next ((iter),(val),NULL)