gts/cdt.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 #ifdef HAVE_CONFIG_H
  21 # include <config.h>
  22 #endif
  23
  24
  25 #include <math.h>
  26 #include "gts.h"
  27
  28 #ifdef USE_SURFACE_BTREE
  29
  30 static gint find_closest (GtsTriangle * t, gpointer value, gpointer * data)
  31 {
  32   guint * ns = data[2];
  33   guint * n = data[3];
  34
  35   if (*n >= *ns)
  36     return TRUE;
  37   else {
  38     gdouble * dmin = data[0];
  39     gpointer * closest = data[1];
  40     GtsPoint * p = data[4];
  41
  42     if (gts_triangle_orientation (t) > 0.) {
  43       GtsPoint * p1 = GTS_POINT (GTS_SEGMENT (t->e1)->v1);
  44       gdouble d = (p->x - p1->x)*(p->x - p1->x) + (p->y - p1->y)*(p->y - p1->y);
  45
  46       if (d < *dmin) {
  47         *dmin = d;
  48         *closest = t;
  49       }
  50       (*n)++;
  51     }
  52   }
  53   return FALSE;
  54 }
  55
  56 /* select the face closest to @p among n^1/3 randomly picked faces
  57  *  of @surface */
  58 static GtsFace * closest_face (GtsSurface * s, GtsPoint * p)
  59 {
  60   guint n = 0, nt, ns;
  61   gdouble dmin = G_MAXDOUBLE;
  62   GtsFace * closest = NULL;
  63   gpointer data[5];
  64
  65   nt = gts_surface_face_number (s);
  66   if (!nt)
  67     return NULL;
  68   ns = exp (log ((gdouble) nt)/3.);
  69
  70   data[0] = &dmin;
  71   data[1] = &closest;
  72   data[2] = &ns;
  73   data[3] = &n;
  74   data[4] = p;
  75   g_tree_traverse (s->faces, (GTraverseFunc) find_closest, G_IN_ORDER, data);
  76
  77   return closest;
  78 }
  79
  80 #else /* not USE_SURFACE_BTREE */
  81
  82 typedef struct _SFindClosest SFindClosest;
  83
  84 struct _SFindClosest {
  85   gdouble dmin;
  86   GtsFace *closest;
  87   GtsPoint * p;
  88   gint stop;
  89 };
  90
  91 #  if GLIB_CHECK_VERSION(2,4,0)
  92 /* finally, with g_hash_table_find we are able to stop iteration over the hash
  93    table in the middle */
  94
  95 static gboolean find_closest (gpointer key, gpointer value, gpointer user_data)
  96 {
  97   SFindClosest * data = (SFindClosest *) user_data;
  98   GtsFace * f = GTS_FACE (value);
  99
 100   if (gts_triangle_orientation (GTS_TRIANGLE (f)) > 0.) {
 101     GtsPoint * p1 = GTS_POINT (GTS_SEGMENT (GTS_TRIANGLE (f)->e1)->v1);
 102     gdouble d = ((data->p->x - p1->x)*(data->p->x - p1->x) +
 103                  (data->p->y - p1->y)*(data->p->y - p1->y));
 104
 105     if (d < data->dmin) {
 106       data->dmin = d;
 107       data->closest = f;
 108     }
 109   }
 110   data->stop--;
 111   return !(data->stop > 0);
 112 }
 113
 114 static GtsFace * closest_face (GtsSurface * s, GtsPoint * p)
 115 {
 116   SFindClosest fc;
 117
 118   fc.dmin = G_MAXDOUBLE;
 119   fc.closest = NULL;
 120   fc.p = p;
 121   fc.stop = (gint) exp (log ((gdouble) g_hash_table_size (s->faces))/3.);
 122   g_hash_table_find (s->faces, find_closest, &fc);
 123
 124   return fc.closest;
 125 }
 126
 127 #  else /* VERSION < 2.4.0 */
 128
 129 static void
 130 find_closest (gpointer key, gpointer value, gpointer user_data)
 131 {
 132   SFindClosest * data = (SFindClosest *) user_data;
 133   GtsFace * f = GTS_FACE (value);
 134
 135   if (gts_triangle_orientation (GTS_TRIANGLE (f)) > 0.) {
 136     GtsPoint * p1 = GTS_POINT (GTS_SEGMENT (GTS_TRIANGLE (f)->e1)->v1);
 137     gdouble d = ((data->p->x - p1->x)*(data->p->x - p1->x) +
 138                  (data->p->y - p1->y)*(data->p->y - p1->y));
 139
 140     if (d < data->dmin) {
 141       data->dmin = d;
 142       data->closest = f;
 143     }
 144   }
 145   data->stop--;
 146 }
 147
 148 /* select the face closest to @p among n^1/3 randomly picked faces
 149  * of @surface */
 150 static GtsFace * closest_face (GtsSurface * s, GtsPoint * p)
 151 {
 152   SFindClosest fc;
 153
 154   if (!g_hash_table_size (s->faces))
 155     return NULL;
 156
 157   fc.dmin = G_MAXDOUBLE;
 158   fc.closest = NULL;
 159   fc.p = p;
 160   fc.stop = (gint) exp (log ((gdouble) g_hash_table_size (s->faces))/3.);
 161   g_hash_table_foreach (s->faces, find_closest, &fc);
 162   return fc.closest;
 163 }
 164 #  endif /* VERSION < 2.4.0 */
 165 #endif /* not USE_SURFACE_BTREE */
 166
 167 /* returns the face belonging to @surface and neighbor of @f via @e */
 168 static GtsFace * neighbor (GtsFace * f,
 169                            GtsEdge * e,
 170                            GtsSurface * surface)
 171 {
 172   GSList * i = e->triangles;
 173   GtsTriangle * t = GTS_TRIANGLE (f);
 174
 175   while (i) {
 176     GtsTriangle * t1 = i->data;
 177     if (t1 != t &&
 178         GTS_IS_FACE (t1) &&
 179         gts_face_has_parent_surface (GTS_FACE (t1), surface))
 180       return GTS_FACE (t1);
 181     i = i->next;
 182   }
 183   return NULL;
 184 }
 185
 186 /* given a triangle @t and a segment s (@o -> @p).
 187    @o must be in @t. Returns the
 188    edge of @t which is intersected by s or %NULL if @p is also
 189    contained in @t (on_summit is set to %FALSE) or if s intersects @t
 190    exactly on one of its summit (on_summit is set to %TRUE). */
 191 static GtsEdge * triangle_next_edge (GtsTriangle * t,
 192                                      GtsPoint * o, GtsPoint * p,
 193                                      gboolean * on_summit)
 194 {
 195   GtsVertex * v1, * v2, * v3;
 196   GtsEdge * e1, * e2, * e3;
 197   gdouble orient = 0.0;
 198
 199   gts_triangle_vertices_edges (t, NULL,
 200                                &v1, &v2, &v3,
 201                                &e1, &e2, &e3);
 202
 203   *on_summit = FALSE;
 204   orient = gts_point_orientation (o, GTS_POINT (v1), p);
 205   if (orient > 0.0) {
 206     orient = gts_point_orientation (o, GTS_POINT (v2), p);
 207     if (orient > 0.0) {
 208       if (gts_point_orientation (GTS_POINT (v2), GTS_POINT (v3), p) >= 0.0)
 209         return NULL;
 210       return e2;
 211     }
 212     if (orient < 0.0) {
 213       if (gts_point_orientation (GTS_POINT (v1), GTS_POINT (v2), p) >= 0.0)
 214         return NULL;
 215       return e1;
 216     }
 217     if (gts_point_orientation (GTS_POINT (v1), GTS_POINT (v2), p) < 0.0)
 218       *on_summit = TRUE;
 219     return NULL;
 220   }
 221
 222   if (orient < 0.0) {
 223     orient = gts_point_orientation (o, GTS_POINT (v3), p);
 224     if (orient > 0.0) {
 225       if (gts_point_orientation (GTS_POINT (v3), GTS_POINT (v1), p) >= 0.0)
 226         return NULL;
 227       return e3;
 228     }
 229     if (orient < 0.0) {
 230       if (gts_point_orientation (GTS_POINT (v2), GTS_POINT (v3), p) >= 0.0)
 231         return NULL;
 232       return e2;
 233     }
 234     if (gts_point_orientation (GTS_POINT (v3), GTS_POINT (v1), p) < 0.0)
 235       *on_summit = TRUE;
 236     return NULL;
 237   }
 238
 239   if (gts_point_orientation (GTS_POINT (v2), GTS_POINT (v3), p) < 0.0)
 240     return e2;
 241   if (gts_point_orientation (GTS_POINT (v1), GTS_POINT (v2), p) < 0.0)
 242     *on_summit = TRUE;
 243   return NULL;
 244 }
 245
 246 static void triangle_barycenter (GtsTriangle * t, GtsPoint * b)
 247 {
 248   GtsPoint * p = GTS_POINT (gts_triangle_vertex (t));
 249   b->x = (p->x +
 250           GTS_POINT (GTS_SEGMENT(t->e1)->v1)->x +
 251           GTS_POINT (GTS_SEGMENT(t->e1)->v2)->x)/3.;
 252   b->y = (p->y +
 253           GTS_POINT (GTS_SEGMENT(t->e1)->v1)->y +
 254           GTS_POINT (GTS_SEGMENT(t->e1)->v2)->y)/3.;
 255 }
 256
 257 static GtsFace * point_locate (GtsPoint * o,
 258                                GtsPoint * p,
 259                                GtsFace * f,
 260                                GtsSurface * surface)
 261 {
 262   GtsEdge * prev;
 263   gboolean on_summit;
 264   GtsVertex * v1, * v2, * v3;
 265   GtsEdge * e2, * e3;
 266
 267   prev = triangle_next_edge (GTS_TRIANGLE (f), o, p, &on_summit);
 268
 269   if (!prev) {
 270     GtsFace * f1;
 271
 272     if (!on_summit)
 273       return f; /* p is inside f */
 274
 275     /* s intersects f exactly on a summit: restarts from a neighbor of f */
 276     if ((f1 = neighbor (f, GTS_TRIANGLE (f)->e1, surface)) ||
 277         (f1 = neighbor (f, GTS_TRIANGLE (f)->e2, surface)) ||
 278         (f1 = neighbor (f, GTS_TRIANGLE (f)->e3, surface))) {
 279       triangle_barycenter (GTS_TRIANGLE (f1), o);
 280       return point_locate (o, p, f1, surface);
 281     }
 282     return NULL;
 283   }
 284
 285   f = neighbor (f, prev, surface);
 286   if (f)
 287     gts_triangle_vertices_edges (GTS_TRIANGLE (f), prev,
 288                                  &v1, &v2, &v3, &prev, &e2, &e3);
 289   while (f) {
 290     gdouble orient = gts_point_orientation (o, GTS_POINT (v3), p);
 291
 292     if (orient < 0.0) {
 293       if (gts_point_orientation (GTS_POINT (v2), GTS_POINT (v3), p) >= 0.0)
 294         return f; /* p is inside f */
 295       f = neighbor (f, e2, surface);
 296       prev = e2;
 297       v1 = v3;
 298     }
 299     else if (orient > 0.0) {
 300       if (gts_point_orientation (GTS_POINT (v3), GTS_POINT (v1), p) >= 0.0)
 301         return f; /* p is inside f */
 302       f = neighbor (f, e3, surface);
 303       prev = e3;
 304       v2 = v3;
 305     }
 306     else {
 307       GtsFace * f1;
 308
 309       if (gts_point_orientation (GTS_POINT (v2), GTS_POINT (v3), p) >= 0.0)
 310         return f; /* p is inside f */
 311
 312       /* s intersects f exactly on v3: restarts from a neighbor of f */
 313       if ((f1 = neighbor (f, e2, surface)) ||
 314           (f1 = neighbor (f, e3, surface))) {
 315         triangle_barycenter (GTS_TRIANGLE (f1), o);
 316         return point_locate (o, p, f1, surface);
 317       }
 318       return NULL;
 319     }
 320     /* update e2, e3, v3 for the new triangle */
 321     if (f) {
 322       if (prev == GTS_TRIANGLE (f)->e1) {
 323         e2 = GTS_TRIANGLE (f)->e2; e3 = GTS_TRIANGLE (f)->e3;
 324       }
 325       else if (prev == GTS_TRIANGLE (f)->e2) {
 326         e2 = GTS_TRIANGLE (f)->e3; e3 = GTS_TRIANGLE (f)->e1;
 327       }
 328       else {
 329         e2 = GTS_TRIANGLE (f)->e1; e3 = GTS_TRIANGLE (f)->e2;
 330       }
 331       if (GTS_SEGMENT (e2)->v1 == v1 || GTS_SEGMENT (e2)->v1 == v2)
 332         v3 = GTS_SEGMENT (e2)->v2;
 333       else
 334         v3 = GTS_SEGMENT (e2)->v1;
 335     }
 336   }
 337   return NULL;
 338 }
 339
 340 /**
 341  * gts_point_locate:
 342  * @p: a #GtsPoint.
 343  * @surface: a #GtsSurface.
 344  * @guess: %NULL or a face of @surface close to @p.
 345  *
 346  * Locates the face of the planar projection of @surface containing
 347  * @p. The planar projection of @surface must define a connected set
 348  * of triangles without holes and bounded by a convex boundary. The
 349  * algorithm is randomized and performs in O(n^1/3) expected time
 350  * where n is the number of triangles of @surface.
 351  *
 352  * If a good @guess is given the point location can be significantly faster.
 353  *
 354  * Returns: a #GtsFace of @surface containing @p or %NULL if @p is not
 355  * contained within the boundary of @surface.
 356  */
 357 GtsFace * gts_point_locate (GtsPoint * p,
 358                             GtsSurface * surface,
 359                             GtsFace * guess)
 360 {
 361   GtsFace * fr;
 362   GtsPoint * o;
 363
 364   g_return_val_if_fail (p != NULL, NULL);
 365   g_return_val_if_fail (surface != NULL, NULL);
 366   g_return_val_if_fail (guess == NULL ||
 367                         gts_face_has_parent_surface (guess, surface), NULL);
 368
 369   if (guess == NULL)
 370     guess = closest_face (surface, p);
 371   else
 372     g_return_val_if_fail (gts_triangle_orientation (GTS_TRIANGLE (guess)) > 0., NULL);
 373
 374   if (guess == NULL)
 375     return NULL;
 376
 377   o = GTS_POINT (gts_object_new (GTS_OBJECT_CLASS (gts_point_class ())));
 378   triangle_barycenter (GTS_TRIANGLE (guess), o);
 379   fr = point_locate (o, p, guess, surface);
 380   gts_object_destroy (GTS_OBJECT (o));
 381
 382   return fr;
 383 }
 384
 385
 386 /**
 387  * gts_constraint_class:
 388  *
 389  * Returns: the #GtsConstraintClass.
 390  */
 391 GtsConstraintClass * gts_constraint_class (void)
 392 {
 393   static GtsConstraintClass * klass = NULL;
 394
 395   if (klass == NULL) {
 396     GtsObjectClassInfo constraint_info = {
 397       "GtsConstraint",
 398       sizeof (GtsConstraint),
 399       sizeof (GtsConstraintClass),
 400       (GtsObjectClassInitFunc) NULL,
 401       (GtsObjectInitFunc) NULL,
 402       (GtsArgSetFunc) NULL,
 403       (GtsArgGetFunc) NULL
 404     };
 405     klass = gts_object_class_new (GTS_OBJECT_CLASS (gts_edge_class ()),
 406                                   &constraint_info);
 407   }
 408
 409   return klass;
 410 }
 411
 412 static void split_list (GtsListFace * f, GtsListFace * f1, GtsListFace * f2,
 413                         GtsPoint * p1, GtsPoint * p2,
 414                         GSList ** last1, GSList ** last2)
 415 {
 416   GSList * i = f->points, * l1 = *last1, * l2 = *last2;
 417
 418   while (i) {
 419     GtsPoint * p = i->data;
 420
 421     if (gts_point_orientation (p1, p2, p) >= 0.) {
 422       if (l1) l1->next = i; else f1->points = i;
 423       l1 = i;
 424     }
 425     else {
 426       if (l2) l2->next = i; else f2->points = i;
 427       l2 = i;
 428     }
 429     i = i->next;
 430   }
 431   f->points = NULL;
 432   *last1 = l1;
 433   *last2 = l2;
 434 }
 435
 436 /* cf. figure misc/swap.fig */
 437 static void swap_if_in_circle (GtsFace * f1,
 438                                GtsVertex * v1,
 439                                GtsVertex * v2,
 440                                GtsVertex * v3,
 441                                GtsEdge * e1,
 442                                GtsEdge * e2,
 443                                GtsEdge * e3,
 444                                GtsSurface * surface)
 445 {
 446   GtsFace * f2;
 447   GtsEdge * e4, *e5;
 448   GtsVertex * v4;
 449
 450   if (GTS_IS_CONSTRAINT (e1)) /* @e1 is a constraint can not swap */
 451     return;
 452
 453   f2 = neighbor (f1, e1, surface);
 454   if (f2 == NULL) /* @e1 is a boundary of @surface */
 455     return;
 456
 457   if (GTS_TRIANGLE (f2)->e1 == e1) {
 458     e4 = GTS_TRIANGLE (f2)->e2; e5 = GTS_TRIANGLE (f2)->e3;
 459   }
 460   else if (GTS_TRIANGLE (f2)->e2 == e1) {
 461     e4 = GTS_TRIANGLE (f2)->e3; e5 = GTS_TRIANGLE (f2)->e1;
 462   }
 463   else {
 464     e4 = GTS_TRIANGLE (f2)->e1; e5 = GTS_TRIANGLE (f2)->e2;
 465   }
 466   if (GTS_SEGMENT (e4)->v1 == GTS_SEGMENT (e1)->v1 ||
 467       GTS_SEGMENT (e4)->v1 == GTS_SEGMENT (e1)->v2)
 468     v4 = GTS_SEGMENT (e4)->v2;
 469   else
 470     v4 = GTS_SEGMENT (e4)->v1;
 471
 472   if (gts_point_in_circle (GTS_POINT (v4), GTS_POINT (v1),
 473                            GTS_POINT (v2), GTS_POINT (v3)) > 0.0) {
 474     GtsEdge * en;
 475     GtsSegment * sn = gts_vertices_are_connected (v3, v4);
 476     GtsFace * f3, * f4;
 477
 478     if (!GTS_IS_EDGE (sn))
 479       en = gts_edge_new (surface->edge_class, v3, v4);
 480     else
 481       en = GTS_EDGE (sn);
 482
 483     f3 = gts_face_new (surface->face_class, en, e5, e2);
 484     gts_object_attributes (GTS_OBJECT (f3), GTS_OBJECT (f1));
 485     f4 = gts_face_new (surface->face_class, en, e3, e4);
 486     gts_object_attributes (GTS_OBJECT (f4), GTS_OBJECT (f2));
 487
 488     if (GTS_IS_LIST_FACE (f3)) {
 489       GSList * last3 = NULL, * last4 = NULL;
 490
 491       if (GTS_IS_LIST_FACE (f1))
 492         split_list (GTS_LIST_FACE (f1), GTS_LIST_FACE (f3), GTS_LIST_FACE (f4),
 493                     GTS_POINT (v3), GTS_POINT (v4), &last3, &last4);
 494       if (GTS_IS_LIST_FACE (f2))
 495         split_list (GTS_LIST_FACE (f2), GTS_LIST_FACE (f3), GTS_LIST_FACE (f4),
 496                     GTS_POINT (v3), GTS_POINT (v4), &last3, &last4);
 497       if (last3) last3->next = NULL;
 498       if (last4) last4->next = NULL;
 499     }
 500
 501     gts_surface_remove_face (surface, f1);
 502     gts_surface_remove_face (surface, f2);
 503     gts_surface_add_face (surface, f3);
 504     gts_surface_add_face (surface, f4);
 505
 506     swap_if_in_circle (f3, v4, v2, v3, e5, e2, en, surface);
 507     swap_if_in_circle (f4, v1, v4, v3, e4, en, e3, surface);
 508   }
 509 }
 510
 511 /**
 512  * gts_delaunay_add_vertex_to_face:
 513  * @surface: a #GtsSurface.
 514  * @v: a #GtsVertex.
 515  * @f: a #GtsFace belonging to @surface.
 516  *
 517  * Adds vertex @v to the face @f of the Delaunay triangulation defined
 518  * by @surface.
 519  *
 520  * Returns: %NULL is @v has been successfully added to @surface or was
 521  * already contained in @surface or a #GtsVertex having the same x and
 522  * y coordinates as @v.
 523  */
 524 GtsVertex * gts_delaunay_add_vertex_to_face (GtsSurface * surface,
 525                                              GtsVertex * v,
 526                                              GtsFace * f)
 527 {
 528   GtsEdge * e1, * e2, * e3;
 529   GtsSegment * s4, * s5, * s6;
 530   GtsEdge * e4, * e5, * e6;
 531   GtsVertex * v1, * v2, * v3;
 532   GtsFace * nf[3];
 533
 534   g_return_val_if_fail (surface != NULL, v);
 535   g_return_val_if_fail (v != NULL, v);
 536   g_return_val_if_fail (f != NULL, v);
 537
 538   gts_triangle_vertices_edges (GTS_TRIANGLE (f), NULL,
 539                                &v1, &v2, &v3, &e1, &e2, &e3);
 540   if (v == v1 || v == v2 || v == v3) /* v already in @surface */
 541     return NULL;
 542   if (GTS_POINT (v)->x == GTS_POINT (v1)->x &&
 543       GTS_POINT (v)->y == GTS_POINT (v1)->y)
 544     return v1;
 545   if (GTS_POINT (v)->x == GTS_POINT (v2)->x &&
 546       GTS_POINT (v)->y == GTS_POINT (v2)->y)
 547     return v2;
 548   if (GTS_POINT (v)->x == GTS_POINT (v3)->x &&
 549       GTS_POINT (v)->y == GTS_POINT (v3)->y)
 550     return v3;
 551
 552   s4 = gts_vertices_are_connected (v, v1);
 553   if (!GTS_IS_EDGE (s4))
 554     e4 = gts_edge_new (surface->edge_class, v, v1);
 555   else
 556     e4 = GTS_EDGE (s4);
 557   s5 = gts_vertices_are_connected (v, v2);
 558   if (!GTS_IS_EDGE (s5))
 559     e5 = gts_edge_new (surface->edge_class, v, v2);
 560   else
 561     e5 = GTS_EDGE (s5);
 562   s6 = gts_vertices_are_connected (v, v3);
 563   if (!GTS_IS_EDGE (s6))
 564     e6 = gts_edge_new (surface->edge_class, v, v3);
 565   else
 566     e6 = GTS_EDGE (s6);
 567
 568   /* cf. figure misc/swap.fig */
 569   nf[0] = gts_face_new (surface->face_class, e4, e1, e5);
 570   gts_object_attributes (GTS_OBJECT (nf[0]), GTS_OBJECT (f));
 571   nf[1] = gts_face_new (surface->face_class, e5, e2, e6);
 572   gts_object_attributes (GTS_OBJECT (nf[1]), GTS_OBJECT (f));
 573   nf[2] = gts_face_new (surface->face_class, e6, e3, e4);
 574   gts_object_attributes (GTS_OBJECT (nf[2]), GTS_OBJECT (f));
 575
 576   if (GTS_IS_LIST_FACE (f) && GTS_IS_LIST_FACE (nf[0])) {
 577     GSList * i = GTS_LIST_FACE (f)->points, * last[3] = { NULL, NULL, NULL };
 578
 579     while (i) {
 580       GtsPoint * p = i->data;
 581       GSList * next = i->next;
 582       guint j;
 583
 584       if (p != GTS_POINT (v)) {
 585         if (gts_point_orientation (GTS_POINT (v), GTS_POINT (v1), p) >= 0.) {
 586           gdouble o = gts_point_orientation (GTS_POINT (v), GTS_POINT (v2), p);
 587
 588           if (o != 0.)
 589             j = o > 0. ? 1 : 0;
 590           else
 591             j = gts_point_orientation (GTS_POINT (v), GTS_POINT (v3), p)
 592               > 0. ? 0 : 1;
 593         }
 594         else if (gts_point_orientation (GTS_POINT (v), GTS_POINT (v3), p) > 0.)
 595           j = 2;
 596         else
 597           j = 1;
 598         if (last[j])
 599           last[j]->next = i;
 600         else
 601           GTS_LIST_FACE (nf[j])->points = i;
 602         last[j] = i;
 603       }
 604       else
 605         g_slist_free_1 (i);
 606       i = next;
 607     }
 608     GTS_LIST_FACE (f)->points = NULL;
 609     if (last[0]) last[0]->next = NULL;
 610     if (last[1]) last[1]->next = NULL;
 611     if (last[2]) last[2]->next = NULL;
 612   }
 613
 614   gts_surface_remove_face (surface, f);
 615   gts_surface_add_face (surface, nf[0]);
 616   gts_surface_add_face (surface, nf[1]);
 617   gts_surface_add_face (surface, nf[2]);
 618
 619   swap_if_in_circle (nf[0], v1, v2, v, e1, e5, e4, surface);
 620   swap_if_in_circle (nf[1], v2, v3, v, e2, e6, e5, surface);
 621   swap_if_in_circle (nf[2], v3, v1, v, e3, e4, e6, surface);
 622
 623   return NULL;
 624 }
 625
 626 /**
 627  * gts_delaunay_add_vertex:
 628  * @surface: a #GtsSurface.
 629  * @v: a #GtsVertex.
 630  * @guess: %NULL or a #GtsFace belonging to @surface to be used as an initial
 631  * guess for point location.
 632  *
 633  * Adds vertex @v to the Delaunay triangulation defined by
 634  * @surface. If @v is not contained in the convex hull bounding
 635  * @surface, @v is not added to the triangulation.
 636  *
 637  * Returns: %NULL is @v has been successfully added to @surface or was
 638  * already contained in @surface, @v if @v is not contained in the
 639  * convex hull bounding surface or a #GtsVertex having the same x and
 640  * y coordinates as @v.
 641  */
 642 GtsVertex * gts_delaunay_add_vertex (GtsSurface * surface,
 643                                      GtsVertex * v,
 644                                      GtsFace * guess)
 645 {
 646   GtsFace * f;
 647
 648   g_return_val_if_fail (surface != NULL, v);
 649   g_return_val_if_fail (v != NULL, v);
 650
 651   if (!(f = gts_point_locate (GTS_POINT (v), surface, guess)))
 652     return v;
 653   return gts_delaunay_add_vertex_to_face (surface, v, f);
 654 }
 655
 656 static gboolean polygon_in_circle (GSList * poly,
 657                                    GtsPoint * p1,
 658                                    GtsPoint * p2,
 659                                    GtsPoint * p3)
 660 {
 661   GtsVertex * v1 = NULL, * v2 = NULL;
 662
 663   while (poly) {
 664     GtsSegment * s = poly->data;
 665     GtsVertex * v;
 666     v = s->v1;
 667     if (v != v1 && v != v2 &&
 668         v != GTS_VERTEX (p1) &&
 669         v != GTS_VERTEX (p2) &&
 670         v != GTS_VERTEX (p3) &&
 671         gts_point_in_circle (GTS_POINT (v), p1, p2, p3) > 0.)
 672       return TRUE;
 673     v = s->v2;
 674     if (v != v1 && v != v2 &&
 675         v != GTS_VERTEX (p1) &&
 676         v != GTS_VERTEX (p2) &&
 677         v != GTS_VERTEX (p3) &&
 678         gts_point_in_circle (GTS_POINT (v), p1, p2, p3) > 0.)
 679       return TRUE;
 680     v1 = s->v1;
 681     v2 = s->v2;
 682     poly = poly->next;
 683   }
 684   return FALSE;
 685 }
 686
 687 static void triangulate_polygon (GSList * poly,
 688                                  GtsSurface * surface,
 689                                  GtsFace * ref)
 690 {
 691   GSList * i, * poly1, * poly2;
 692   GtsVertex * v1, * v2, * v3 = NULL;
 693   gboolean found = FALSE;
 694   GtsSegment * s, * s1, * s2;
 695   GtsEdge * e1, * e2;
 696   GtsFace * f;
 697
 698   if (poly == NULL || poly->next == NULL) {
 699     g_slist_free (poly);
 700     return;
 701   }
 702
 703   s = poly->data;
 704   s1 = poly->next->data;
 705   if (s->v1 == s1->v1 || s->v1 == s1->v2) {
 706     v1 = s->v2;
 707     v2 = s->v1;
 708   }
 709   else {
 710     g_assert (s->v2 == s1->v1 || s->v2 == s1->v2);
 711     v1 = s->v1;
 712     v2 = s->v2;
 713   }
 714
 715   i = poly->next;
 716   v3 = v2;
 717   while (i && !found) {
 718     s1 = i->data;
 719     if (s1->v1 == v3)
 720       v3 = s1->v2;
 721     else {
 722       g_assert (s1->v2 == v3);
 723       v3 = s1->v1;
 724     }
 725     if (v3 != v1 &&
 726         gts_point_orientation (GTS_POINT (v1),
 727                                GTS_POINT (v2),
 728                                GTS_POINT (v3)) >= 0. &&
 729         !polygon_in_circle (poly,
 730                             GTS_POINT (v1),
 731                             GTS_POINT (v2),
 732                             GTS_POINT (v3)))
 733       found = TRUE;
 734     else
 735       i = i->next;
 736   }
 737
 738   if (!found) {
 739     g_slist_free (poly);
 740     return;
 741   }
 742
 743   s1 = gts_vertices_are_connected (v2, v3);
 744   if (!GTS_IS_EDGE (s1))
 745     e1 = gts_edge_new (surface->edge_class, v2, v3);
 746   else
 747     e1 = GTS_EDGE (s1);
 748   s2 = gts_vertices_are_connected (v3, v1);
 749   if (!GTS_IS_EDGE (s2))
 750     e2 = gts_edge_new (surface->edge_class, v3, v1);
 751   else
 752     e2 = GTS_EDGE (s2);
 753   f = gts_face_new (surface->face_class, GTS_EDGE (s), e1, e2);
 754   gts_object_attributes (GTS_OBJECT (f), GTS_OBJECT (ref));
 755   gts_surface_add_face (surface, f);
 756
 757   poly1 = poly->next;
 758   g_slist_free_1 (poly);
 759   if (i->next && e2 != i->next->data)
 760     poly2 = g_slist_prepend (i->next, e2);
 761   else
 762     poly2 = i->next;
 763   if (e1 != i->data)
 764     i->next = g_slist_prepend (NULL, e1);
 765   else
 766     i->next = NULL;
 767
 768  triangulate_polygon (poly1, surface, ref);
 769  triangulate_polygon (poly2, surface, ref);
 770 }
 771
 772 /**
 773  * gts_delaunay_remove_vertex:
 774  * @surface: a #GtsSurface.
 775  * @v: a #GtsVertex.
 776  *
 777  * Removes @v from the Delaunay triangulation defined by @surface and
 778  * restores the Delaunay property. Vertex @v must not be used by any
 779  * constrained edge otherwise the triangulation is not guaranteed to
 780  * be Delaunay.
 781  */
 782 void gts_delaunay_remove_vertex (GtsSurface * surface, GtsVertex * v)
 783 {
 784   GSList * triangles, * i;
 785   GtsFace * ref = NULL;
 786
 787   g_return_if_fail (surface != NULL);
 788   g_return_if_fail (v != NULL);
 789
 790   i = triangles = gts_vertex_triangles (v, NULL);
 791   while (i && !ref) {
 792     if (GTS_IS_FACE (i->data) &&
 793         gts_face_has_parent_surface (i->data, surface))
 794       ref = i->data;
 795     i = i->next;
 796   }
 797   if (!ref) {
 798     g_slist_free (triangles);
 799     g_return_if_fail (ref);
 800   }
 801   triangulate_polygon (gts_vertex_fan_oriented (v, surface), surface, ref);
 802   i = triangles;
 803   while (i) {
 804     if (GTS_IS_FACE (i->data) &&
 805         gts_face_has_parent_surface (i->data, surface))
 806       gts_surface_remove_face (surface, i->data);
 807     i = i->next;
 808   }
 809   g_slist_free (triangles);
 810 }
 811
 812 #define NEXT_CUT(edge, edge1, list) { next = neighbor (f, edge, surface);\
 813                                       remove_triangles (e, surface);\
 814                                       if (!constraint && !e->triangles)\
 815                                         gts_object_destroy (GTS_OBJECT (e));\
 816                                       g_assert (next);\
 817                                       *list = g_slist_prepend (*list, edge1);\
 818                                       return g_slist_concat (constraint,\
 819                                         remove_intersected_edge (s, edge,\
 820                                                next, surface, left, right));\
 821                                     }
 822
 823 static void remove_triangles (GtsEdge * e, GtsSurface * s)
 824 {
 825   GSList * i = e->triangles;
 826
 827   while (i) {
 828     GSList * next = i->next;
 829
 830     if (GTS_IS_FACE (i->data) && gts_face_has_parent_surface (i->data, s))
 831       gts_surface_remove_face (s, i->data);
 832     i = next;
 833   }
 834 }
 835
 836 static GSList *
 837 remove_intersected_edge (GtsSegment * s,
 838                          GtsEdge * e,
 839                          GtsFace * f,
 840                          GtsSurface * surface,
 841                          GSList ** left, GSList ** right)
 842 {
 843   GtsVertex * v1, * v2, * v3;
 844   GtsEdge * e1, * e2;
 845   gdouble o1, o2;
 846   GtsFace * next;
 847   GSList * constraint = NULL;
 848
 849   if (GTS_IS_CONSTRAINT (e))
 850     constraint = g_slist_prepend (NULL, e);
 851
 852   gts_triangle_vertices_edges (GTS_TRIANGLE (f), e,
 853                                &v1, &v2, &v3, &e, &e1, &e2);
 854
 855   o1 = gts_point_orientation (GTS_POINT (v2), GTS_POINT (v3),
 856                               GTS_POINT (s->v2));
 857   o2 = gts_point_orientation (GTS_POINT (v3), GTS_POINT (v1),
 858                               GTS_POINT (s->v2));
 859
 860   if (o1 == 0. && o2 == 0.) {
 861 /*    if(o2 != 0.) {
 862       fprintf(stderr, "o1 = %f o2 = %f\n", o1, o2);
 863       fprintf(stderr, "v1 = %f, %f\n", GTS_POINT(v1)->x, GTS_POINT(v1)->y);
 864       fprintf(stderr, "v2 = %f, %f\n", GTS_POINT(v2)->x, GTS_POINT(v2)->y);
 865       fprintf(stderr, "v3 = %f, %f\n", GTS_POINT(v3)->x, GTS_POINT(v3)->y);
 866       fprintf(stderr, "s->v2 = %f, %f\n", GTS_POINT(s->v2)->x, GTS_POINT(s->v2)->y);
 867
 868       g_assert (o2 == 0.);
 869     }*/
 870  //   if(o2 == 0.) {
 871       remove_triangles (e, surface);
 872       if (!constraint && !e->triangles)
 873         gts_object_destroy (GTS_OBJECT (e));
 874       *left = g_slist_prepend (*left, e2);
 875       *right = g_slist_prepend (*right, e1);
 876 //    }
 877   }
 878   else if (o1 > 0.) {
 879     g_assert (o2 <= 0.);
 880     NEXT_CUT (e2, e1, right)
 881   }
 882   else if (o2 >= 0.)
 883     NEXT_CUT (e1, e2, left)
 884   else {
 885     gdouble o3 = gts_point_orientation (GTS_POINT (s->v1), GTS_POINT (s->v2),
 886                                         GTS_POINT (v3));
 887     if (o3 > 0.)
 888       NEXT_CUT (e1, e2, left)
 889     else
 890       NEXT_CUT (e2, e1, right)
 891   }
 892   return constraint;
 893 }
 894
 895 static GSList *
 896 remove_intersected_vertex (GtsSegment * s,
 897                            GtsVertex * v,
 898                            GtsSurface * surface,
 899                            GSList ** left,
 900                            GSList ** right,
 901                            GtsFace ** ref)
 902 {
 903   GSList * triangles = gts_vertex_triangles (v, NULL);
 904   GSList * i;
 905
 906   i = triangles;
 907   while (i) {
 908     GtsTriangle * t = i->data;
 909     if (GTS_IS_FACE (t) &&
 910         gts_face_has_parent_surface (GTS_FACE (t), surface)) {
 911       GtsVertex * v1, * v2, * v3;
 912       gdouble o1, o2;
 913
 914       gts_triangle_vertices (t, &v1, &v2, &v3);
 915       if (v == v2) {
 916         v2 = v3;
 917         v3 = v1;
 918       }
 919       else if (v == v3) {
 920         v3 = v2;
 921         v2 = v1;
 922       }
 923       else
 924         g_assert (v == v1);
 925
 926       if ((o1 = gts_point_orientation (GTS_POINT (v), GTS_POINT (v2),
 927                                        GTS_POINT (s->v2))) >= 0. &&
 928           (o2 = gts_point_orientation (GTS_POINT (v3), GTS_POINT (v),
 929                                        GTS_POINT (s->v2))) >= 0.) {
 930         gdouble o3 = gts_point_orientation (GTS_POINT (v2), GTS_POINT (v3),
 931                                             GTS_POINT (s->v2));
 932         GtsEdge * e = gts_triangle_edge_opposite (t, v);
 933         GtsEdge * e1, * e2;
 934         GtsFace * next = neighbor (GTS_FACE (t), e, surface);
 935
 936         *ref = GTS_FACE (t);
 937         gts_triangle_vertices_edges (t, e, &v2, &v3, &v, &e, &e2, &e1);
 938
 939         g_slist_free (triangles);
 940
 941         if (o3 >= 0.) /* @s->v2 is inside (or on the edge) of t */
 942           return NULL;
 943
 944         gts_allow_floating_faces = TRUE;
 945         gts_surface_remove_face (surface, GTS_FACE (t));
 946         gts_allow_floating_faces = FALSE;
 947
 948         *left = g_slist_prepend (*left, e2);
 949         *right = g_slist_prepend (*right, e1);
 950
 951         g_assert (next);
 952         return remove_intersected_edge (s, e, next, surface, left, right);
 953       }
 954     }
 955     i = i->next;
 956   }
 957
 958   g_assert_not_reached ();
 959   return NULL;
 960 }
 961
 962 /**
 963  * gts_delaunay_add_constraint:
 964  * @surface: a #GtsSurface.
 965  * @c: a #GtsConstraint.
 966  *
 967  * Add constraint @c to the constrained Delaunay triangulation defined by
 968  * @surface.
 969  *
 970  * Returns: a list of #GtsConstraint conflicting (i.e. intersecting) with @c
 971  * which were removed from @surface (%NULL if there was none).
 972  */
 973 GSList * gts_delaunay_add_constraint (GtsSurface * surface,
 974                                       GtsConstraint * c)
 975 {
 976   GSList * constraints;
 977   GtsVertex * v1; //, * v2;
 978   GSList * left = NULL, * right = NULL;
 979   GtsFace * ref = NULL;
 980
 981   g_return_val_if_fail (surface != NULL, NULL);
 982   g_return_val_if_fail (c != NULL, NULL);
 983   g_return_val_if_fail (GTS_IS_CONSTRAINT (c), NULL);
 984
 985   v1 = GTS_SEGMENT (c)->v1;
 986   //v2 = GTS_SEGMENT (c)->v2;
 987
 988   gts_allow_floating_edges = TRUE;
 989   constraints = remove_intersected_vertex (GTS_SEGMENT (c), v1, surface,
 990                                            &left, &right, &ref);
 991   gts_allow_floating_edges = FALSE;
 992 #if 1
 993   triangulate_polygon (g_slist_prepend (g_slist_reverse (right), c),
 994                        surface, ref);
 995   triangulate_polygon (g_slist_prepend (left, c),
 996                        surface, ref);
 997 #else
 998   right = g_slist_prepend (g_slist_reverse (right), c);
 999   left = g_slist_prepend (left, c);
1000   {
1001     FILE * fp0 = fopen ("hole", "wt");
1002     FILE * fp1 = fopen ("right", "wt");
1003     FILE * fp2 = fopen ("left", "wt");
1004     GSList * i = left;
1005
1006     gts_surface_write (surface, fp0);
1007     fclose (fp0);
1008
1009     fprintf (fp2, "LIST {\n");
1010     while (i) {
1011       GtsSegment * s = i->data;
1012       fprintf (fp2,
1013                "# %p: %p->%p\n"
1014                "VECT 1 2 0 2 0 %g %g 0 %g %g 0\n",
1015                s, s->v1, s->v2,
1016                GTS_POINT (s->v1)->x, GTS_POINT (s->v1)->y,
1017                GTS_POINT (s->v2)->x, GTS_POINT (s->v2)->y);
1018       i = i->next;
1019     }
1020     fprintf (fp2, "}\n");
1021     fprintf (fp1, "LIST {\n");
1022     i = right;
1023     while (i) {
1024       GtsSegment * s = i->data;
1025       fprintf (fp1,
1026                "# %p: %p->%p\n"
1027                "VECT 1 2 0 2 0 %g %g 0 %g %g 0\n",
1028                s, s->v1, s->v2,
1029                GTS_POINT (s->v1)->x, GTS_POINT (s->v1)->y,
1030                GTS_POINT (s->v2)->x, GTS_POINT (s->v2)->y);
1031       i = i->next;
1032     }
1033     fprintf (fp1, "}\n");
1034     fclose (fp1);
1035     fclose (fp2);
1036   }
1037   triangulate_polygon (right, surface);
1038   triangulate_polygon (left, surface);
1039 #endif
1040   if (ref && !ref->surfaces) {
1041     gts_allow_floating_edges = TRUE;
1042     gts_object_destroy (GTS_OBJECT (ref));
1043     gts_allow_floating_edges = FALSE;
1044   }
1045   return constraints;
1046 }
1047
1048 static void delaunay_check (GtsTriangle * t, gpointer * data)
1049 {
1050   GtsSurface * surface = data[0];
1051   GtsFace ** face = data[1];
1052
1053   if (*face == NULL) {
1054     GSList * i, * list;
1055     GtsVertex * v1, * v2, * v3;
1056
1057     gts_triangle_vertices (t, &v1, &v2, &v3);
1058     list = gts_vertex_neighbors (v1, NULL, surface);
1059     list = gts_vertex_neighbors (v2, list, surface);
1060     list = gts_vertex_neighbors (v3, list, surface);
1061     i = list;
1062     while (i && *face == NULL) {
1063       GtsVertex * v = i->data;
1064       if (v != v1 && v != v2 && v != v3 &&
1065           gts_point_in_circle (GTS_POINT (v),
1066                                GTS_POINT (v1),
1067                                GTS_POINT (v2),
1068                                GTS_POINT (v3)) > 0.)
1069         *face = GTS_FACE (t);
1070       i = i->next;
1071     }
1072     g_slist_free (list);
1073   }
1074 }
1075
1076 /**
1077  * gts_delaunay_check:
1078  * @surface: a #GtsSurface.
1079  *
1080  * Returns: %NULL if the planar projection of @surface is a Delaunay
1081  * triangulation (unconstrained), a #GtsFace violating the Delaunay
1082  * property otherwise.
1083  */
1084 GtsFace * gts_delaunay_check (GtsSurface * surface)
1085 {
1086   GtsFace * face = NULL;
1087   gpointer data[2];
1088
1089   g_return_val_if_fail (surface != NULL, FALSE);
1090
1091   data[0] = surface;
1092   data[1] = &face;
1093   gts_surface_foreach_face (surface, (GtsFunc) delaunay_check, data);
1094
1095   return face;
1096 }
1097
1098 /**
1099  * gts_delaunay_remove_hull:
1100  * @surface: a #GtsSurface.
1101  *
1102  * Removes all the edges of the boundary of @surface which are not
1103  * constraints.
1104  */
1105 void gts_delaunay_remove_hull (GtsSurface * surface)
1106 {
1107   GSList * boundary;
1108
1109   g_return_if_fail (surface != NULL);
1110
1111   boundary = gts_surface_boundary (surface);
1112   gts_allow_floating_edges = TRUE;
1113   while (boundary) {
1114     GSList * i = boundary;
1115     GtsEdge * e = i->data;
1116
1117     boundary = i->next;
1118     g_slist_free_1 (i);
1119     if (!GTS_IS_CONSTRAINT (e)) {
1120       GtsTriangle * t = GTS_TRIANGLE (gts_edge_is_boundary (e, surface));
1121
1122       if (t != NULL) {
1123         if (t->e1 != e && !GTS_IS_CONSTRAINT (t->e1) &&
1124             !gts_edge_is_boundary (t->e1, surface))
1125           boundary = g_slist_prepend (boundary, t->e1);
1126         if (t->e2 != e && !GTS_IS_CONSTRAINT (t->e2) &&
1127             !gts_edge_is_boundary (t->e2, surface))
1128           boundary = g_slist_prepend (boundary, t->e2);
1129         if (t->e3 != e && !GTS_IS_CONSTRAINT (t->e3) &&
1130             !gts_edge_is_boundary (t->e3, surface))
1131           boundary = g_slist_prepend (boundary, t->e3);
1132         gts_surface_remove_face (surface, GTS_FACE (t));
1133       }
1134       if (!e->triangles)
1135         gts_object_destroy (GTS_OBJECT (e));
1136     }
1137   }
1138   gts_allow_floating_edges = FALSE;
1139 }
1140
1141 /* GtsListFace: Object */
1142
1143 static void gts_list_face_destroy (GtsObject * object)
1144 {
1145   g_slist_free (GTS_LIST_FACE (object)->points);
1146
1147   (* GTS_OBJECT_CLASS (gts_list_face_class ())->parent_class->destroy)
1148     (object);
1149 }
1150
1151 static void gts_list_face_class_init (GtsFaceClass * klass)
1152 {
1153   GTS_OBJECT_CLASS (klass)->destroy = gts_list_face_destroy;
1154 }
1155
1156 GtsFaceClass * gts_list_face_class (void)
1157 {
1158   static GtsFaceClass * klass = NULL;
1159
1160   if (klass == NULL) {
1161     GtsObjectClassInfo gts_list_face_info = {
1162       "GtsListFace",
1163       sizeof (GtsListFace),
1164       sizeof (GtsFaceClass),
1165       (GtsObjectClassInitFunc) gts_list_face_class_init,
1166       (GtsObjectInitFunc) NULL,
1167       (GtsArgSetFunc) NULL,
1168       (GtsArgGetFunc) NULL
1169     };
1170     klass = gts_object_class_new (GTS_OBJECT_CLASS (gts_face_class ()),
1171                                   &gts_list_face_info);
1172   }
1173
1174   return klass;
1175 }