gts/refine.c

   1 /* GTS - Library for the manipulation of triangulated surfaces
   2  * Copyright (C) 1999 Stéphane Popinet
   3  *
   4  * This library is free software; you can redistribute it and/or
   5  * modify it under the terms of the GNU Library General Public
   6  * License as published by the Free Software Foundation; either
   7  * version 2 of the License, or (at your option) any later version.
   8  *
   9  * This library is distributed in the hope that it will be useful,
  10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  12  * Library General Public License for more details.
  13  *
  14  * You should have received a copy of the GNU Library General Public
  15  * License along with this library; if not, write to the
  16  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  17  * Boston, MA 02111-1307, USA.
  18  */
  19
  20 #include <math.h>
  21 #include "gts.h"
  22
  23 /**
  24  * gts_vertex_encroaches_edge:
  25  * @v: a #GtsVertex.
  26  * @e: a #GtsEdge.
  27  *
  28  * Returns: %TRUE if @v is strictly contained in the diametral circle of @e,
  29  * %FALSE otherwise.
  30  */
  31 gboolean gts_vertex_encroaches_edge (GtsVertex * v, GtsEdge * e)
  32 {
  33   GtsPoint * p, * p1, * p2;
  34
  35   g_return_val_if_fail (v != NULL, FALSE);
  36   g_return_val_if_fail (e != NULL, FALSE);
  37
  38   p = GTS_POINT (v);
  39   p1 = GTS_POINT (GTS_SEGMENT (e)->v1);
  40   p2 = GTS_POINT (GTS_SEGMENT (e)->v2);
  41
  42   if ((p1->x - p->x)*(p2->x - p->x) + (p1->y - p->y)*(p2->y - p->y) < 0.0)
  43     return TRUE;
  44   return FALSE;
  45 }
  46
  47 /**
  48  * gts_edge_is_encroached:
  49  * @e: a #GtsEdge.
  50  * @s: a #GtsSurface describing a (constrained) Delaunay triangulation.
  51  * @encroaches: a #GtsEncroachFunc.
  52  * @data: user data to be passed to @encroaches.
  53  *
  54  * Returns: a #GtsVertex belonging to @s and encroaching upon @e
  55  * (as defined by @encroaches) or %NULL if there is none.
  56  */
  57 GtsVertex * gts_edge_is_encroached (GtsEdge * e,
  58                                     GtsSurface * s,
  59                                     GtsEncroachFunc encroaches,
  60                                     gpointer data)
  61 {
  62   GSList * i;
  63
  64   g_return_val_if_fail (e != NULL, NULL);
  65   g_return_val_if_fail (s != NULL, NULL);
  66   g_return_val_if_fail (encroaches != NULL, NULL);
  67
  68   i = e->triangles;
  69   while (i) {
  70     GtsFace * f = i->data;
  71     if (GTS_IS_FACE (f) && gts_face_has_parent_surface (f, s)) {
  72       GtsVertex * v = gts_triangle_vertex_opposite (GTS_TRIANGLE (f), e);
  73       if ((* encroaches) (v, e, s, data))
  74         return v;
  75     }
  76     i = i->next;
  77   }
  78
  79   return NULL;
  80 }
  81
  82 #define ALREADY_ENCROACHED(c) (GTS_OBJECT (c)->reserved)
  83
  84 static void vertex_encroaches (GtsVertex * v,
  85                                GtsSurface * surface,
  86                                GtsFifo * encroached,
  87                                GtsEncroachFunc encroaches,
  88                                gpointer data)
  89 {
  90   GSList * triangles, * i;
  91
  92   g_return_if_fail (v != NULL);
  93   g_return_if_fail (surface != NULL);
  94   g_return_if_fail (encroached != NULL);
  95   g_return_if_fail (encroaches != NULL);
  96
  97   i = triangles = gts_vertex_triangles (v, NULL);
  98   while (i) {
  99     GtsFace * f = i->data;
 100     if (GTS_IS_FACE (f) && gts_face_has_parent_surface (f, surface)) {
 101       GtsEdge * e = gts_triangle_edge_opposite (i->data, v);
 102       if (!ALREADY_ENCROACHED (e) &&
 103           GTS_IS_CONSTRAINT (e) &&
 104           (* encroaches) (v, e, surface, data)) {
 105         gts_fifo_push (encroached, e);
 106         ALREADY_ENCROACHED (e) = encroached;
 107       }
 108     }
 109     i = i->next;
 110   }
 111   g_slist_free (triangles);
 112 }
 113
 114 static void make_encroached_fifo (GtsEdge * e, gpointer * datas)
 115 {
 116   GtsFifo * fifo = datas[0];
 117   GtsSurface * s = datas[1];
 118   GtsEncroachFunc encroaches = (GtsEncroachFunc) datas[2];
 119   gpointer data = datas[3];
 120
 121   if (GTS_IS_CONSTRAINT (e) &&
 122       gts_edge_is_encroached (e, s, encroaches, data)) {
 123     gts_fifo_push (fifo, e);
 124     ALREADY_ENCROACHED (e) = fifo;
 125   }
 126 }
 127
 128 #define SQUARE_ROOT_TWO 1.41421356237309504880168872420969807856967187
 129 #define DISTANCE_2D(v1, v2) (sqrt ((GTS_POINT (v2)->x - GTS_POINT (v1)->x)*\
 130                                    (GTS_POINT (v2)->x - GTS_POINT (v1)->x) +\
 131                                    (GTS_POINT (v2)->y - GTS_POINT (v1)->y)*\
 132                                    (GTS_POINT (v2)->y - GTS_POINT (v1)->y)))
 133
 134 /* finds where to split the given edge to avoid infinite cycles. (see
 135    Shewchuk's thesis for details */
 136 static GtsVertex * split_edge (GtsEdge * e,
 137                                GtsSurface * surface)
 138 {
 139   GSList * i = e->triangles;
 140   GtsEdge * c = NULL;
 141
 142   /* look for constraints touching e */
 143   while (i && !c) {
 144     GtsTriangle * t = i->data;
 145     if (GTS_IS_FACE (t) &&
 146         gts_face_has_parent_surface (GTS_FACE (t), surface)) {
 147       GtsEdge * e1, * e2;
 148       if (t->e1 == e) { e1 = t->e2; e2 = t->e3; }
 149       else if (t->e2 == e) { e1 = t->e1; e2 = t->e3; }
 150       else { e1 = t->e1; e2 = t->e2; }
 151       if (GTS_IS_CONSTRAINT (e1) && !GTS_IS_CONSTRAINT (e2))
 152         c = e1;
 153       else if (GTS_IS_CONSTRAINT (e2) && !GTS_IS_CONSTRAINT (e1))
 154         c = e2;
 155     }
 156     i = i->next;
 157   }
 158   if (c) {
 159     /* use power of two concentric shells */
 160     GtsVertex * v1 = GTS_SEGMENT (e)->v1;
 161     GtsVertex * v2 = GTS_SEGMENT (e)->v2;
 162     gdouble l = DISTANCE_2D (v1, v2);
 163     gdouble nearestpower = 1., split;
 164
 165     while (l > SQUARE_ROOT_TWO*nearestpower)
 166       nearestpower *= 2.;
 167     while (l < SQUARE_ROOT_TWO*nearestpower/2.)
 168       nearestpower /= 2.;
 169     split = nearestpower/l/2.;
 170
 171     if (GTS_SEGMENT (c)->v1 == v2 || GTS_SEGMENT (c)->v2 == v2)
 172       split = 1. - split;
 173     return gts_vertex_new (surface->vertex_class,
 174                            (1. - split)*GTS_POINT (v1)->x +
 175                            split*GTS_POINT (v2)->x,
 176                            (1. - split)*GTS_POINT (v1)->y +
 177                            split*GTS_POINT (v2)->y,
 178                            (1. - split)*GTS_POINT (v1)->z +
 179                            split*GTS_POINT (v2)->z);
 180   }
 181   else
 182     return gts_segment_midvertex (GTS_SEGMENT (e), surface->vertex_class);
 183 }
 184
 185 static gint split_encroached (GtsSurface * surface,
 186                               GtsFifo * encroached,
 187                               gint steiner_max,
 188                               GtsEncroachFunc encroaches,
 189                               gpointer data)
 190 {
 191   GtsSegment * s;
 192
 193   while (steiner_max-- != 0 && (s = gts_fifo_pop (encroached))) {
 194     GtsVertex *add_vertex_returned;
 195     GtsVertex * v = split_edge (GTS_EDGE (s), surface);
 196     GtsFace * boundary = gts_edge_is_boundary (GTS_EDGE (s), surface);
 197     GtsFace * f = boundary;
 198 #if 1
 199     GtsEdge * e1 = GTS_EDGE (gts_object_clone (GTS_OBJECT (s)));
 200     GtsEdge * e2 = GTS_EDGE (gts_object_clone (GTS_OBJECT (s)));
 201
 202     GTS_SEGMENT (e1)->v1 = s->v1;
 203     s->v1->segments = g_slist_prepend (s->v1->segments, e1);
 204     GTS_SEGMENT (e1)->v2 = v;
 205     v->segments = g_slist_prepend (v->segments, e1);
 206
 207     GTS_SEGMENT (e2)->v1 = v;
 208     v->segments = g_slist_prepend (v->segments, e2);
 209     GTS_SEGMENT (e2)->v2 = s->v2;
 210     s->v2->segments = g_slist_prepend (s->v2->segments, e2);
 211 #else
 212     GtsEdge * e1 = gts_edge_new (GTS_EDGE_CLASS (GTS_OBJECT (s)->klass),
 213                                  s->v1, v);
 214     GtsEdge * e2 = gts_edge_new (GTS_EDGE_CLASS (GTS_OBJECT (s)->klass),
 215                                  v, s->v2);
 216 #endif
 217
 218     GTS_OBJECT (s)->klass = GTS_OBJECT_CLASS (surface->edge_class);
 219
 220     if (f == NULL)
 221       f = gts_edge_has_parent_surface (GTS_EDGE (s), surface);
 222     g_assert (f != NULL);
 223     add_vertex_returned = gts_delaunay_add_vertex_to_face (surface, v, f);
 224     g_assert (add_vertex_returned == NULL);
 225
 226     if (boundary)
 227       gts_object_destroy (GTS_OBJECT (s));
 228
 229     vertex_encroaches (v, surface, encroached, encroaches, data);
 230
 231     if (gts_edge_is_encroached (e1, surface, encroaches, data)) {
 232       gts_fifo_push (encroached, e1);
 233       ALREADY_ENCROACHED (e1) = encroached;
 234     }
 235     if (gts_edge_is_encroached (e2, surface, encroaches, data)) {
 236       gts_fifo_push (encroached, e2);
 237       ALREADY_ENCROACHED (e2) = encroached;
 238     }
 239   }
 240
 241   return steiner_max;
 242 }
 243
 244 /**
 245  * gts_delaunay_conform:
 246  * @surface: a #GtsSurface describing a constrained Delaunay triangulation.
 247  * @steiner_max: maximum number of Steiner points.
 248  * @encroaches: a #GtsEncroachFunc.
 249  * @data: user-data to pass to @encroaches.
 250  *
 251  * Recursively split constraints of @surface which are encroached by
 252  * vertices of @surface (see Shewchuk 96 for details). The split
 253  * constraints are destroyed and replaced by a set of new constraints
 254  * of the same class. If gts_vertex_encroaches_edge() is used for
 255  * @encroaches, the resulting surface will be Delaunay conforming.
 256  *
 257  * If @steiner_max is positive or nul, the recursive splitting
 258  * procedure will stop when this maximum number of Steiner points is
 259  * reached. In that case the resulting surface will not necessarily be
 260  * Delaunay conforming.
 261  *
 262  * Returns: the number of remaining encroached edges. If @steiner_max
 263  * is set to a negative value and gts_vertex_encroaches_edge() is used
 264  * for @encroaches this should always be zero.
 265  */
 266 guint gts_delaunay_conform (GtsSurface * surface,
 267                             gint steiner_max,
 268                             GtsEncroachFunc encroaches,
 269                             gpointer data)
 270 {
 271   GtsFifo * encroached;
 272   gpointer datas[4];
 273   guint encroached_number;
 274
 275   g_return_val_if_fail (surface != NULL, 0);
 276   g_return_val_if_fail (surface != NULL, 0);
 277   g_return_val_if_fail (encroaches != NULL, 0);
 278
 279   datas[0] = encroached = gts_fifo_new ();
 280   datas[1] = surface;
 281   datas[2] = encroaches;
 282   datas[3] = data;
 283   gts_surface_foreach_edge (surface, (GtsFunc) make_encroached_fifo, datas);
 284
 285   split_encroached (surface,
 286                     encroached,
 287                     steiner_max,
 288                     encroaches, data);
 289   gts_fifo_foreach (encroached, (GtsFunc) gts_object_reset_reserved, NULL);
 290   encroached_number = gts_fifo_size (encroached);
 291   gts_fifo_destroy (encroached);
 292   return encroached_number;
 293 }
 294
 295 #define EHEAP_PAIR(f) (GTS_OBJECT (f)->reserved)
 296
 297 static void heap_surface_add_face (GtsSurface * s, GtsFace * f)
 298 {
 299   GtsEHeap * heap = GTS_OBJECT (s)->reserved;
 300   gdouble key = gts_eheap_key (heap, f);
 301
 302   if (key != 0.)
 303     EHEAP_PAIR (f) = gts_eheap_insert_with_key (heap, f, key);
 304
 305   if (GTS_SURFACE_CLASS (GTS_OBJECT (s)->klass->parent_class)->add_face)
 306     (* GTS_SURFACE_CLASS (GTS_OBJECT (s)->klass->parent_class)->add_face)
 307       (s, f);
 308 }
 309
 310 static void heap_surface_remove_face (GtsSurface * s, GtsFace * f)
 311 {
 312   GtsEHeap * heap = GTS_OBJECT (s)->reserved;
 313
 314   if (EHEAP_PAIR (f))
 315     gts_eheap_remove (heap, EHEAP_PAIR (f));
 316
 317   if (GTS_SURFACE_CLASS (GTS_OBJECT (s)->klass->parent_class)->remove_face)
 318     (* GTS_SURFACE_CLASS (GTS_OBJECT (s)->klass->parent_class)->remove_face)
 319       (s, f);
 320 }
 321
 322 static void heap_surface_class_init (GtsSurfaceClass * klass)
 323 {
 324   klass->add_face = heap_surface_add_face;
 325   klass->remove_face = heap_surface_remove_face;
 326 }
 327
 328 static GtsObjectClass * heap_surface_class_new (GtsObjectClass * parent_class)
 329 {
 330   GtsObjectClassInfo heap_surface_info;
 331
 332   heap_surface_info = parent_class->info;
 333   heap_surface_info.class_init_func = (GtsObjectClassInitFunc)
 334     heap_surface_class_init;
 335   return gts_object_class_new (parent_class,
 336                                &heap_surface_info);
 337 }
 338
 339 static void make_face_heap (GtsFace * f, GtsEHeap * heap)
 340 {
 341   gdouble key = gts_eheap_key (heap, f);
 342
 343   if (key != 0.)
 344     EHEAP_PAIR (f) = gts_eheap_insert_with_key (heap, f, key);
 345 }
 346
 347 /**
 348  * gts_delaunay_refine:
 349  * @surface: a #GtsSurface describing a conforming Delaunay triangulation.
 350  * @steiner_max: maximum number of Steiner points.
 351  * @encroaches: a #GtsEncroachFunc.
 352  * @encroach_data: user-data to pass to @encroaches.
 353  * @cost: a #GtsKeyFunc used to sort the faces during refinement.
 354  * @cost_data: user-data to pass to @cost.
 355  *
 356  * An implementation of the refinement algorithm described in Ruppert
 357  * (1995) and Shewchuk (1996).
 358  *
 359  * Returns: the number of unrefined faces of @surface left. Should be zero
 360  * if @steiner_max is set to a negative value.
 361  */
 362 guint gts_delaunay_refine (GtsSurface * surface,
 363                            gint steiner_max,
 364                            GtsEncroachFunc encroaches,
 365                            gpointer encroach_data,
 366                            GtsKeyFunc cost,
 367                            gpointer cost_data)
 368 {
 369   GtsObjectClass * heap_surface_class;
 370   GtsObjectClass * original_class;
 371   GtsEHeap * heap;
 372   GtsFifo * encroached;
 373   GtsFace * f;
 374   guint unrefined_number;
 375
 376   g_return_val_if_fail (surface != NULL, 0);
 377   g_return_val_if_fail (encroaches != NULL, 0);
 378   g_return_val_if_fail (cost != NULL, 0);
 379
 380   original_class = GTS_OBJECT (surface)->klass;
 381   heap_surface_class = heap_surface_class_new (original_class);
 382   GTS_OBJECT (surface)->klass = heap_surface_class;
 383
 384   heap = gts_eheap_new (cost, cost_data);
 385   gts_surface_foreach_face (surface, (GtsFunc) make_face_heap, heap);
 386   encroached = gts_fifo_new ();
 387
 388   GTS_OBJECT (surface)->reserved = heap;
 389
 390   while (steiner_max-- != 0 && (f = gts_eheap_remove_top (heap, NULL))) {
 391     GtsVertex *add_vertex_returned;
 392     GtsVertex * c =
 393       GTS_VERTEX (gts_triangle_circumcircle_center (GTS_TRIANGLE (f),
 394                   GTS_POINT_CLASS (surface->vertex_class)));
 395     EHEAP_PAIR (f) = NULL;
 396     g_assert (c != NULL);
 397     add_vertex_returned = gts_delaunay_add_vertex (surface, c, f);
 398     g_assert (add_vertex_returned == NULL);
 399
 400     vertex_encroaches (c, surface, encroached, encroaches, encroach_data);
 401     if (!gts_fifo_is_empty (encroached)) {
 402       gts_delaunay_remove_vertex (surface, c);
 403       steiner_max = split_encroached (surface,
 404                                       encroached,
 405                                       steiner_max,
 406                                       encroaches,
 407                                       encroach_data);
 408     }
 409   }
 410
 411   unrefined_number = gts_eheap_size (heap);
 412   gts_eheap_foreach (heap, (GFunc) gts_object_reset_reserved, NULL);
 413   gts_eheap_destroy (heap);
 414
 415   gts_fifo_foreach (encroached, (GtsFunc) gts_object_reset_reserved, NULL);
 416   gts_fifo_destroy (encroached);
 417
 418   GTS_OBJECT (surface)->klass = original_class;
 419   GTS_OBJECT (surface)->reserved = NULL;
 420   g_free (heap_surface_class);
 421
 422   return unrefined_number;
 423 }