gts/boolean.c

   1 /* GTS - Library for the manipulation of triangulated surfaces
   2  * Copyright (C) 1999--2002 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 /*#define DEBUG*/
  24 /*#define DEBUG_BOOLEAN*/
  25 /*#define CHECK_ORIENTED*/
  26
  27 #ifdef DEBUG
  28 #  include "gts-private.h"
  29 #endif /* DEBUG */
  30
  31 static void surface_inter_destroy (GtsObject * object)
  32 {
  33   GtsSurfaceInter * si = GTS_SURFACE_INTER (object);
  34
  35   gts_object_destroy (GTS_OBJECT (si->s1));
  36   gts_object_destroy (GTS_OBJECT (si->s2));
  37   g_slist_free (si->edges);
  38
  39   (* GTS_OBJECT_CLASS (gts_surface_inter_class ())->parent_class->destroy)
  40     (object);
  41 }
  42
  43 static void surface_inter_class_init (GtsObjectClass * klass)
  44 {
  45   klass->destroy = surface_inter_destroy;
  46 }
  47
  48 static void surface_inter_init (GtsSurfaceInter * si)
  49 {
  50   si->s1 = si->s2 = NULL;
  51   si->edges = NULL;
  52 }
  53
  54 /**
  55  * gts_surface_inter_class:
  56  *
  57  * Returns: the #GtsSurfaceInterClass.
  58  */
  59 GtsSurfaceInterClass * gts_surface_inter_class (void)
  60 {
  61   static GtsSurfaceInterClass * klass = NULL;
  62
  63   if (klass == NULL) {
  64     GtsObjectClassInfo surface_inter_info = {
  65       "GtsSurfaceInter",
  66       sizeof (GtsSurfaceInter),
  67       sizeof (GtsSurfaceInterClass),
  68       (GtsObjectClassInitFunc) surface_inter_class_init,
  69       (GtsObjectInitFunc) surface_inter_init,
  70       (GtsArgSetFunc) NULL,
  71       (GtsArgGetFunc) NULL
  72     };
  73     klass = gts_object_class_new (gts_object_class (), &surface_inter_info);
  74   }
  75
  76   return klass;
  77 }
  78
  79 /* EdgeInter: Header */
  80
  81 typedef struct _EdgeInter         EdgeInter;
  82
  83 struct _EdgeInter {
  84   GtsEdge parent;
  85
  86   GtsTriangle * t1, * t2;
  87 };
  88
  89 #define EDGE_INTER(obj)            GTS_OBJECT_CAST (obj,\
  90                                                  EdgeInter,\
  91                                                  edge_inter_class ())
  92 #define IS_EDGE_INTER(obj)         (gts_object_is_from_class (obj,\
  93                                                  edge_inter_class ()))
  94
  95 static GtsEdgeClass * edge_inter_class  (void);
  96 static EdgeInter * edge_inter_new    (GtsVertex * v1, GtsVertex * v2,
  97                                       GtsTriangle * t1, GtsTriangle * t2);
  98
  99 /* EdgeInter: Object */
 100
 101 static GtsEdgeClass * edge_inter_class (void)
 102 {
 103   static GtsEdgeClass * klass = NULL;
 104
 105   if (klass == NULL) {
 106     GtsObjectClassInfo edge_inter_info = {
 107       "EdgeInter",
 108       sizeof (EdgeInter),
 109       sizeof (GtsEdgeClass),
 110       (GtsObjectClassInitFunc) NULL,
 111       (GtsObjectInitFunc) NULL,
 112       (GtsArgSetFunc) NULL,
 113       (GtsArgGetFunc) NULL
 114     };
 115     klass = gts_object_class_new (GTS_OBJECT_CLASS (gts_constraint_class ()),
 116                                   &edge_inter_info);
 117   }
 118
 119   return klass;
 120 }
 121
 122 static EdgeInter * edge_inter_new (GtsVertex * v1, GtsVertex * v2,
 123                                    GtsTriangle * t1, GtsTriangle * t2)
 124 {
 125   EdgeInter * object;
 126
 127   object = EDGE_INTER (gts_edge_new (GTS_EDGE_CLASS (edge_inter_class ()),
 128                                      v1, v2));
 129   object->t1 = t1;
 130   object->t2 = t2;
 131
 132   return object;
 133 }
 134
 135 #ifdef DEBUG
 136 static void write_surface_graph (GtsSurface * s, FILE * fp)
 137 {
 138   GSList * l = NULL;
 139   GtsGraph * g;
 140   static void add_to_list (gpointer data, GSList ** l) {
 141     *l = g_slist_prepend (*l, data);
 142   }
 143
 144   gts_surface_foreach_vertex (s, (GtsFunc) gts_object_reset_reserved, NULL);
 145   gts_surface_foreach_edge (s, (GtsFunc) gts_object_reset_reserved, NULL);
 146   gts_surface_foreach_edge (s, (GtsFunc) add_to_list, &l);
 147   g = gts_segments_graph_new (gts_graph_class (), l);
 148   gts_graph_write_dot (g, fp);
 149   gts_object_destroy (GTS_OBJECT (g));
 150   g_slist_free (l);
 151 }
 152 #endif /* DEBUG */
 153
 154 static GtsPoint * segment_triangle_intersection (GtsSegment * s,
 155                                                  GtsTriangle * t,
 156                                                  GtsPointClass * klass)
 157 {
 158   GtsPoint * A, * B, * C, * D, * E;
 159   gint ABCE, ABCD, ADCE, ABDE, BCDE;
 160   GtsEdge * AB, * BC, * CA;
 161   gdouble a, b, c;
 162
 163   g_return_val_if_fail (s != NULL, NULL);
 164   g_return_val_if_fail (t != NULL, NULL);
 165   g_return_val_if_fail (klass != NULL, NULL);
 166
 167   gts_triangle_vertices_edges (t, NULL,
 168                                (GtsVertex **) &A,
 169                                (GtsVertex **) &B,
 170                                (GtsVertex **) &C,
 171                                &AB, &BC, &CA);
 172   D = GTS_POINT (s->v1);
 173   E = GTS_POINT (s->v2);
 174
 175   ABCE = gts_point_orientation_3d_sos (A, B, C, E);
 176   ABCD = gts_point_orientation_3d_sos (A, B, C, D);
 177   if (ABCE < 0 || ABCD > 0) {
 178     GtsPoint * tmpp;
 179     gint tmp;
 180
 181     tmpp = E; E = D; D = tmpp;
 182     tmp = ABCE; ABCE = ABCD; ABCD = tmp;
 183   }
 184   if (ABCE < 0 || ABCD > 0)
 185     return NULL;
 186   ADCE = gts_point_orientation_3d_sos (A, D, C, E);
 187   if (ADCE < 0)
 188     return NULL;
 189   ABDE = gts_point_orientation_3d_sos (A, B, D, E);
 190   if (ABDE < 0)
 191     return NULL;
 192   BCDE = gts_point_orientation_3d_sos (B, C, D, E);
 193   if (BCDE < 0)
 194     return NULL;
 195   a = gts_point_orientation_3d (A, B, C, E);
 196   b = gts_point_orientation_3d (A, B, C, D);
 197   if (a != b) {
 198     c = a/(a - b);
 199     return gts_point_new (klass,
 200                           E->x + c*(D->x - E->x),
 201                           E->y + c*(D->y - E->y),
 202                           E->z + c*(D->z - E->z));
 203   }
 204   /* D and E are contained within ABC */
 205 #ifdef DEBUG
 206   fprintf (stderr,
 207            "segment: %p:%s triangle: %p:%s intersection\n"
 208            "D and E contained in ABC\n",
 209            s, GTS_NEDGE (s)->name, t, GTS_NFACE (t)->name);
 210 #endif /* DEBUG */
 211   g_assert (a == 0.);
 212   return gts_point_new (klass,
 213                         (E->x + D->x)/2.,
 214                         (E->y + D->y)/2.,
 215                         (E->z + D->z)/2.);
 216 }
 217
 218 static gint triangle_triangle_orientation (GtsPoint * p1,
 219                                            GtsPoint * p2, GtsPoint * p3,
 220                                            GtsPoint * p4, GtsPoint * p5,
 221                                            GtsPoint * p6)
 222 {
 223   gint o4 = 0, o5 = 0, o6 = 0;
 224
 225   if (p4 != p1 && p4 != p2 && p4 != p3)
 226     o4 = gts_point_orientation_3d_sos (p1, p2, p3, p4);
 227   if (p5 != p1 && p5 != p2 && p5 != p3)
 228     o5 = gts_point_orientation_3d_sos (p1, p2, p3, p5);
 229   if (o4*o5 < 0)
 230     return 0;
 231   if (p6 != p1 && p6 != p2 && p6 != p3)
 232     o6 = gts_point_orientation_3d_sos (p1, p2, p3, p6);
 233   if (o4*o6 < 0 || o5*o6 < 0)
 234     return 0;
 235   if (o4) return o4;
 236   if (o5) return o5;
 237   g_assert (o6);
 238   return o6;
 239 }
 240
 241 static gint triangle_point_orientation (GtsTriangle * t1,
 242                                         GtsTriangle * t2,
 243                                         gint o1,
 244                                         GtsPoint * p)
 245 {
 246   GtsPoint * p1 = GTS_POINT (GTS_SEGMENT (t1->e1)->v1);
 247   GtsPoint * p2 = GTS_POINT (GTS_SEGMENT (t1->e1)->v2);
 248   GtsPoint * p3 = GTS_POINT (gts_triangle_vertex (t1));
 249   GtsPoint * p4 = GTS_POINT (GTS_SEGMENT (t2->e1)->v1);
 250   GtsPoint * p5 = GTS_POINT (GTS_SEGMENT (t2->e1)->v2);
 251   GtsPoint * p6 = GTS_POINT (gts_triangle_vertex (t2));
 252   gint o = triangle_triangle_orientation (p1, p2, p3, p4, p5, p6);
 253
 254   if (o != 0)
 255     return o;
 256   o = triangle_triangle_orientation (p4, p5, p6, p1, p2, p3);
 257   if (o != 0) {
 258     gint o2 = gts_point_orientation_3d_sos (p4, p5, p6, p);
 259
 260     return - o*o1*o2;
 261   }
 262   return 0;
 263 }
 264
 265 static void add_edge_inter (GtsEdge * e,
 266                             GtsTriangle * t,
 267                             GtsVertex * v)
 268 {
 269   GtsVertex * ev1 = GTS_SEGMENT (e)->v1, * ev2 = GTS_SEGMENT (e)->v2;
 270   GList * i = GTS_OBJECT (e)->reserved;
 271
 272   GTS_OBJECT (v)->reserved = t;
 273   if (i == NULL) {
 274     GTS_OBJECT (e)->reserved = g_list_prepend (NULL, v);
 275 #ifdef DEBUG
 276     fprintf (stderr, "add_edge_inter: inserting %p (%p,%p)\n", v, e, t);
 277 #endif /* DEBUG */
 278   }
 279   else {
 280     GtsPoint * p1 = GTS_POINT (GTS_SEGMENT (t->e1)->v1);
 281     GtsPoint * p2 = GTS_POINT (GTS_SEGMENT (t->e1)->v2);
 282     GtsPoint * p3 = GTS_POINT (gts_triangle_vertex (t));
 283     gint o1, oref = gts_point_orientation_3d_sos (p1, p2, p3, GTS_POINT (ev1));
 284
 285     o1 = oref;
 286     while (i) {
 287       gint o2 = triangle_point_orientation (t, GTS_OBJECT (i->data)->reserved,
 288                                             oref, GTS_POINT (ev1));
 289
 290       if (o2 == 0) {
 291 #ifdef DEBUG
 292         g_warning ("add_edge_inter: safe sign evaluation failed\n");
 293 #endif /* DEBUG */
 294         o2 = gts_point_orientation_3d_sos (p1, p2, p3, i->data);
 295       }
 296
 297       if (o1*o2 < 0)
 298         break;
 299       ev1 = i->data;
 300       o1 = o2;
 301       i = i->next;
 302     }
 303     if (i != NULL) {
 304       GList * n = g_list_prepend (NULL, v);
 305
 306       ev2 = i->data;
 307       n->next = i;
 308       n->prev = i->prev;
 309       i->prev = n;
 310       if (n->prev == NULL)
 311         GTS_OBJECT (e)->reserved = n;
 312       else
 313         n->prev->next = n;
 314     }
 315     else {
 316       g_assert (o1*gts_point_orientation_3d_sos (p1, p2, p3, GTS_POINT (ev2))
 317                 < 0);
 318       GTS_OBJECT (e)->reserved = g_list_append (GTS_OBJECT (e)->reserved, v);
 319     }
 320 #ifdef DEBUG
 321     fprintf (stderr,
 322              "add_edge_inter: inserting %p (%p,%p) between %p and %p\n",
 323              v, e, t, ev1, ev2);
 324     i = GTS_OBJECT (e)->reserved;
 325     while (i) {
 326       fprintf (stderr, " %p", i->data);
 327       i = i->next;
 328     }
 329     fprintf (stderr, "\n");
 330 #endif /* DEBUG */
 331   }
 332 }
 333
 334 static GtsVertex * intersects (GtsEdge * e,
 335                                GtsTriangle * t,
 336                                GtsSurface * s)
 337 {
 338   GList * i = GTS_OBJECT (e)->reserved;
 339   GtsVertex * v;
 340
 341   while (i) {
 342     if (GTS_OBJECT (i->data)->reserved == t)
 343       return i->data;
 344     i = i->next;
 345   }
 346
 347   v = GTS_VERTEX (segment_triangle_intersection (GTS_SEGMENT (e), t,
 348                                     GTS_POINT_CLASS (s->vertex_class)));
 349   if (v != NULL) {
 350 #ifdef DEBUG
 351     if (GTS_IS_NVERTEX (v) && GTS_IS_NEDGE (e) && GTS_IS_NFACE (t) &&
 352         GTS_NVERTEX (v)->name[0] == '\0')
 353       g_snprintf (GTS_NVERTEX (v)->name, GTS_NAME_LENGTH, "%s|%s",
 354                   GTS_NEDGE (e)->name, GTS_NFACE (t)->name);
 355 #endif /* DEBUG */
 356     if (s->vertex_class->intersection_attributes)
 357       (*s->vertex_class->intersection_attributes)
 358         (v, GTS_OBJECT (e), GTS_OBJECT (t));
 359     add_edge_inter (e, t, v);
 360   }
 361   return v;
 362 }
 363
 364 /* see figure misc/orientation.fig */
 365 static gint intersection_orientation (GtsTriangle * t1,
 366                                       GtsEdge * e,
 367                                       GtsTriangle * t2)
 368 {
 369   GtsVertex * v1, * v2, * v3;
 370   GtsEdge * e2, * e3;
 371   GtsVertex * v4, * v5, * v6;
 372
 373   gts_triangle_vertices_edges (t1, e, &v1, &v2, &v3, &e, &e2, &e3);
 374   gts_triangle_vertices (t2, &v4, &v5, &v6);
 375
 376   return gts_point_orientation_3d_sos (GTS_POINT (v4),
 377                                        GTS_POINT (v5),
 378                                        GTS_POINT (v6),
 379                                        GTS_POINT (v2));
 380 }
 381
 382 #define UPDATE_ORIENTATION if (o > 0) { vi2 = v; /* e2 = e; */ } else { vi2 = vi1;\
 383                                                                         /* e2 = e1; */\
 384                                                                         vi1 = v;\
 385                                                                         /* e1 = e; */ }
 386
 387 static void intersect_edges (GtsBBox * bb1, GtsBBox * bb2,
 388                              GtsSurfaceInter * si)
 389 {
 390   GtsSurface * s1 = GTS_OBJECT (si->s1)->reserved;
 391   GtsTriangle * t1 = GTS_TRIANGLE (bb1->bounded);
 392   GtsTriangle * t2 = GTS_TRIANGLE (bb2->bounded);
 393   GtsVertex * v, * vi1 = NULL, * vi2 = NULL;
 394   //GtsEdge * e1 = NULL, * e2 = NULL, * e;
 395
 396   vi1 = intersects (t2->e1, t1, s1);
 397   //e1 = t2->e1;
 398   v = intersects (t2->e2, t1, s1);
 399   //e = t2->e2;
 400   if (!vi1) {
 401     vi1 = v;
 402     //e1 = e;
 403   }
 404   else if (v) {
 405     gint o = intersection_orientation (t2, t2->e2, t1);
 406     UPDATE_ORIENTATION;
 407   }
 408   if (!vi2) {
 409     v = intersects (t2->e3, t1, s1);
 410     //e = t2->e3;
 411     if (!vi1) {
 412       vi1 = v;
 413       //e1 = e;
 414     }
 415     else if (v) {
 416       gint o = intersection_orientation (t2, t2->e3, t1);
 417       UPDATE_ORIENTATION;
 418     }
 419   }
 420   if (!vi2) {
 421     v = intersects (t1->e1, t2, s1);
 422     //e = t1->e1;
 423     if (!vi1) {
 424       vi1 = v;
 425       //e1 = e;
 426     }
 427     else if (v) {
 428       gint o = - intersection_orientation (t1, t1->e1, t2);
 429       UPDATE_ORIENTATION;
 430     }
 431   }
 432   if (!vi2) {
 433     v = intersects (t1->e2, t2, s1);
 434     //e = t1->e2;
 435     if (!vi1) {
 436       vi1 = v;
 437       //e1 = e;
 438     }
 439     else if (v) {
 440       gint o = - intersection_orientation (t1, t1->e2, t2);
 441       UPDATE_ORIENTATION;
 442     }
 443   }
 444   if (!vi2) {
 445     v = intersects (t1->e3, t2, s1);
 446     //e = t1->e3;
 447     if (!vi1) {
 448       vi1 = v;
 449       //e1 = e;
 450     }
 451     else if (v) {
 452       gint o = - intersection_orientation (t1, t1->e3, t2);
 453       UPDATE_ORIENTATION;
 454     }
 455   }
 456
 457   g_assert ((!vi1 && !vi2) || (vi1 && vi2));
 458   if (vi1) {
 459     GtsEdge * e = GTS_EDGE (edge_inter_new (vi1, vi2, t1, t2));
 460
 461 #ifdef DEBUG
 462     fprintf (stderr, "creating constraint %p: %p->%p: %p/%p\n",
 463              e, vi1, vi2, t1, t2);
 464 #endif /* DEBUG */
 465     gts_surface_add_face (si->s1, GTS_FACE (t1));
 466     gts_surface_add_face (si->s2, GTS_FACE (t2));
 467     si->edges = g_slist_prepend (si->edges, e);
 468     GTS_OBJECT (t1)->reserved = g_slist_prepend (GTS_OBJECT (t1)->reserved, e);
 469     GTS_OBJECT (t2)->reserved = g_slist_prepend (GTS_OBJECT (t2)->reserved, e);
 470   }
 471 }
 472
 473 static GtsSurfaceInter * surface_inter_new (GtsSurfaceInterClass * klass,
 474                                             GtsSurface * s1,
 475                                             GtsSurface * s2,
 476                                             GNode * faces_tree1,
 477                                             GNode * faces_tree2)
 478 {
 479   GtsSurfaceInter * si;
 480
 481   si = GTS_SURFACE_INTER (gts_object_new (GTS_OBJECT_CLASS (klass)));
 482   si->s1 = gts_surface_new (gts_surface_class (),
 483                             s1->face_class,
 484                             s1->edge_class,
 485                             s1->vertex_class);
 486   GTS_OBJECT (si->s1)->reserved = s1;
 487   si->s2 = gts_surface_new (gts_surface_class (),
 488                             s2->face_class,
 489                             s2->edge_class,
 490                             s2->vertex_class);
 491   GTS_OBJECT (si->s2)->reserved = s2;
 492   gts_bb_tree_traverse_overlapping (faces_tree1, faces_tree2,
 493                                     (GtsBBTreeTraverseFunc) intersect_edges,
 494                                     si);
 495
 496   return si;
 497 }
 498
 499 static void free_slist (GtsObject * o)
 500 {
 501   g_slist_free (o->reserved);
 502   o->reserved = NULL;
 503 }
 504
 505 static void free_glist (GtsObject * o)
 506 {
 507   g_list_foreach (o->reserved, (GFunc) gts_object_reset_reserved, NULL);
 508   g_list_free (o->reserved);
 509   o->reserved = NULL;
 510 }
 511
 512 /**
 513  * gts_surface_intersection:
 514  * @s1: a #GtsSurface.
 515  * @s2: a #GtsSurface.
 516  * @faces_tree1: a bounding box tree (see gts_bb_tree_new()) for
 517  * the faces of @s1.
 518  * @faces_tree2: a bounding box tree for the faces of @s2.
 519  *
 520  * Returns: a list of #GtsEdge defining the curve intersection of the
 521  * two surfaces.
 522  */
 523 GSList * gts_surface_intersection (GtsSurface * s1,
 524                                    GtsSurface * s2,
 525                                    GNode * faces_tree1,
 526                                    GNode * faces_tree2)
 527 {
 528   GtsSurfaceInter * si;
 529   GSList * inter;
 530
 531   g_return_val_if_fail (s1 != NULL, NULL);
 532   g_return_val_if_fail (s2 != NULL, NULL);
 533   g_return_val_if_fail (faces_tree1 != NULL, NULL);
 534   g_return_val_if_fail (faces_tree2 != NULL, NULL);
 535
 536   si = surface_inter_new (gts_surface_inter_class (),
 537                           s1, s2, faces_tree1, faces_tree2);
 538
 539   gts_surface_foreach_face (si->s1, (GtsFunc) free_slist, NULL);
 540   gts_surface_foreach_face (si->s2, (GtsFunc) free_slist, NULL);
 541   gts_surface_foreach_edge (si->s1, (GtsFunc) free_glist, NULL);
 542   gts_surface_foreach_edge (si->s2, (GtsFunc) free_glist, NULL);
 543   inter = si->edges;
 544   si->edges = NULL;
 545   gts_object_destroy (GTS_OBJECT (si));
 546
 547   return inter;
 548 }
 549
 550 typedef enum {
 551   INTERIOR = 1 << (GTS_USER_FLAG),
 552   RELEVANT = 1 << (GTS_USER_FLAG + 1)
 553 } CurveFlag;
 554
 555 #define IS_SET(s, f) ((GTS_OBJECT_FLAGS (s) & (f)) != 0)
 556 #define SET(s, f)   (GTS_OBJECT_FLAGS (s) |= (f))
 557 #define UNSET(s, f) (GTS_OBJECT_FLAGS (s) &= ~(f))
 558 #define NEXT(s)  (GTS_OBJECT (s)->reserved)
 559
 560 #ifdef DEBUG
 561 static void print_segment (GtsSegment * s)
 562 {
 563   fprintf (stderr, "%p: %s->%s ", s,
 564            GTS_NVERTEX (s->v1)->name,
 565            GTS_NVERTEX (s->v2)->name);
 566   if (NEXT (s)) {
 567     GtsSegment * next = NEXT (s);
 568
 569     fprintf (stderr, "next %p: %s->%s\n", next,
 570              GTS_NVERTEX (next->v1)->name,
 571              GTS_NVERTEX (next->v2)->name);
 572   }
 573   else
 574     fprintf (stderr, "next NULL\n");
 575 }
 576
 577 static void write_nodes (GSList * i, GHashTable * hash, guint * nn,
 578                          FILE * fp)
 579 {
 580   while (i) {
 581     GtsSegment * s = i->data;
 582
 583     if (!g_hash_table_lookup (hash, s->v1)) {
 584       fprintf (fp, "  %u [ label = \"%p\" ];\n", *nn, s->v1);
 585       g_hash_table_insert (hash, s->v1, GUINT_TO_POINTER ((*nn)++));
 586     }
 587     if (!g_hash_table_lookup (hash, s->v2)) {
 588       fprintf (fp, "  %u [ label = \"%p\" ];\n", *nn, s->v2);
 589       g_hash_table_insert (hash, s->v2, GUINT_TO_POINTER ((*nn)++));
 590     }
 591     i = i->next;
 592   }
 593 }
 594
 595 static void write_edges (GSList * i, GHashTable * hash,
 596                          GtsSurface * surface,
 597                          FILE * fp)
 598 {
 599   while (i) {
 600     GtsSegment * s = i->data;
 601
 602     fprintf (fp, "  %u -> %u [ label = \"%p:%d\" ];\n",
 603              GPOINTER_TO_UINT (g_hash_table_lookup (hash, s->v1)),
 604              GPOINTER_TO_UINT (g_hash_table_lookup (hash, s->v2)),
 605              s,
 606              gts_edge_face_number (GTS_EDGE (s), surface));
 607     i = i->next;
 608   }
 609 }
 610
 611 static void write_graph (GSList * boundary, GSList * interior,
 612                          GtsSurface * surface,
 613                          FILE * fp)
 614 {
 615   GHashTable * hash = g_hash_table_new (NULL, NULL);
 616   guint nn = 1;
 617
 618   fprintf (fp, "digraph oriented_curve {\n");
 619   write_nodes (boundary, hash, &nn, fp);
 620   write_nodes (interior, hash, &nn, fp);
 621   write_edges (boundary, hash, surface, fp);
 622   fprintf (fp, "  edge [ color = red ];\n");
 623   write_edges (interior, hash, surface, fp);
 624   fprintf (fp, "}\n");
 625   g_hash_table_destroy (hash);
 626 }
 627
 628 static void write_graph1 (GtsSegment * start, GSList * i,
 629                           GtsSurface * surface,
 630                           FILE * fp)
 631 {
 632   GSList * boundary = NULL, * interior = NULL;
 633   GtsSegment * s = start;
 634
 635   do {
 636     boundary = g_slist_prepend (boundary, s);
 637     s = NEXT (s);
 638   } while (s != start);
 639   while (i) {
 640     if (IS_SET (i->data, INTERIOR))
 641       interior = g_slist_prepend (interior, i->data);
 642     i = i->next;
 643   }
 644   write_graph (boundary, interior, surface, fp);
 645   g_slist_free (boundary);
 646   g_slist_free (interior);
 647 }
 648
 649 static void print_loop (GtsSegment * start, FILE * fp)
 650 {
 651   GtsSegment * s = start;
 652
 653   do {
 654     fprintf (fp, "  %p: %p:%s -> %p:%s\n",
 655              s,
 656              s->v1, GTS_NVERTEX (s->v1)->name,
 657              s->v2, GTS_NVERTEX (s->v2)->name);
 658     s = NEXT (s);
 659   } while (s != start && s != NULL);
 660 }
 661
 662 static void draw_vector (GtsPoint * p1, GtsPoint * p2, FILE * fp)
 663 {
 664   gdouble x = p2->x - p1->x;
 665   gdouble y = p2->y - p1->y;
 666   gdouble z = p2->z - p1->z;
 667
 668   fprintf (fp, "VECT 1 3 0 3 0 %g %g %g %g %g %g %g %g %g\n",
 669            p1->x + x - (x - y/2.)/5.,
 670            p1->y + y - (x/2. + y)/5.,
 671            p1->z + z - (x/2. + z)/5.,
 672            p1->x + x,
 673            p1->y + y,
 674            p1->z + z,
 675            p1->x + x - (x + y/2.)/5.,
 676            p1->y + y + (x/2. - y)/5.,
 677            p1->z + z + (x/2. - z)/5.);
 678   fprintf (fp, "VECT 1 2 0 2 0 %g %g %g %g %g %g\n",
 679            p1->x, p1->y, p1->z,
 680            p1->x + x,
 681            p1->y + y,
 682            p1->z + z);
 683 }
 684
 685 static void draw_vector1 (GtsPoint * p1, GtsPoint * p2, GtsPoint * o,
 686                           FILE * fp)
 687 {
 688   gdouble x1 = o->x + 0.9*(p1->x - o->x);
 689   gdouble y1 = o->y + 0.9*(p1->y - o->y);
 690   gdouble z1 = o->z + 0.9*(p1->z - o->z);
 691   gdouble x2 = o->x + 0.9*(p2->x - o->x);
 692   gdouble y2 = o->y + 0.9*(p2->y - o->y);
 693   gdouble z2 = o->z + 0.9*(p2->z - o->z);
 694   gdouble x = x2 - x1;
 695   gdouble y = y2 - y1;
 696   gdouble z = z2 - z1;
 697
 698   fprintf (fp, "VECT 1 3 0 3 0 %g %g %g %g %g %g %g %g %g\n",
 699            x1 + x - (x - y/2.)/5.,
 700            y1 + y - (x/2. + y)/5.,
 701            z1 + z - (x/2. + z)/5.,
 702            x1 + x,
 703            y1 + y,
 704            z1 + z,
 705            x1 + x - (x + y/2.)/5.,
 706            y1 + y + (x/2. - y)/5.,
 707            z1 + z + (x/2. - z)/5.);
 708   fprintf (fp, "VECT 1 2 0 2 0 %g %g %g %g %g %g\n",
 709            x1, y1, z1,
 710            x1 + x,
 711            y1 + y,
 712            z1 + z);
 713 }
 714
 715 static void write_segments (GSList * boundary, GSList * interior,
 716                             FILE * fp)
 717 {
 718   GSList * i = boundary;
 719
 720   fprintf (fp, "LIST {\n");
 721   while (i) {
 722     GSList * inext = i->next;
 723     GtsSegment * s = i->data;
 724     GtsSegment * next = inext ? inext->data : boundary->data;
 725     GtsVertex * v1, * v2;
 726
 727     if (s->v1 != next->v1 && s->v1 != next->v2) {
 728       v1 = s->v1;
 729       v2 = s->v2;
 730     }
 731     else {
 732       v1 = s->v2;
 733       v2 = s->v1;
 734     }
 735     draw_vector (GTS_POINT (v1), GTS_POINT (v2), fp);
 736     i = inext;
 737   }
 738   i = interior;
 739   while (i) {
 740     GtsSegment * s = i->data;
 741
 742     draw_vector (GTS_POINT (s->v1), GTS_POINT (s->v2), fp);
 743     i = i->next;
 744   }
 745   fprintf (fp, "}\n");
 746 }
 747
 748 static void write_loops (GSList * i, FILE * fp)
 749 {
 750   guint nl = 0;
 751
 752   while (i) {
 753     GtsSegment * start = i->data, * s;
 754     GtsPoint os;
 755     guint n = 0;
 756
 757     fprintf (fp, "(geometry \"loop%d\" = LIST {\n", nl++);
 758
 759     os.x = os.y = os.z = 0.;
 760     s = start;
 761     do {
 762       GtsSegment * next = NEXT (s);
 763       GtsPoint * p;
 764
 765       if (s->v1 != next->v1 && s->v1 != next->v2)
 766         p = GTS_POINT (s->v1);
 767        else
 768          p = GTS_POINT (s->v2);
 769       os.x += p->x; os.y += p->y; os.z += p->z; n++;
 770       s = next;
 771      } while (s != start);
 772     os.x /= n; os.y /= n; os.z /= n;
 773
 774     s = start;
 775     do {
 776       GtsSegment * next = NEXT (s);
 777
 778       if (s->v1 != next->v1 && s->v1 != next->v2)
 779         draw_vector1 (GTS_POINT (s->v1), GTS_POINT (s->v2), &os, fp);
 780       else
 781          draw_vector1 (GTS_POINT (s->v2), GTS_POINT (s->v1), &os, fp);
 782       s = next;
 783     } while (s != start);
 784
 785     fprintf (fp, "})\n");
 786
 787     i = i->next;
 788   }
 789 }
 790
 791 #define NAME(v) (GTS_IS_NVERTEX (v) ? GTS_NVERTEX (v)->name : "")
 792 #endif /* DEBUG */
 793
 794 static GtsSegment * prev_flag (GtsSegment * s, CurveFlag flag)
 795 {
 796   GSList * i = s->v1->segments;
 797
 798   while (i) {
 799     if (i->data != s && IS_SET (i->data, flag))
 800       return i->data;
 801     i = i->next;
 802   }
 803   return NULL;
 804 }
 805
 806 static GtsSegment * next_flag (GtsSegment * s, CurveFlag flag)
 807 {
 808   GSList * i = s->v2->segments;
 809
 810   while (i) {
 811     if (i->data != s && IS_SET (i->data, flag))
 812       return i->data;
 813     i = i->next;
 814   }
 815   return NULL;
 816 }
 817
 818 static GtsSegment * next_interior (GtsVertex * v)
 819 {
 820   GSList * i = v->segments;
 821
 822   while (i) {
 823     GtsSegment * s = i->data;
 824
 825     if (s->v1 == v && IS_SET (s, INTERIOR))
 826       return s;
 827     i = i->next;
 828   }
 829   return NULL;
 830 }
 831
 832 static GtsSegment * prev_interior (GtsVertex * v)
 833 {
 834   GSList * i = v->segments;
 835
 836   while (i) {
 837     GtsSegment * s = i->data;
 838
 839     if (s->v2 == v && IS_SET (s, INTERIOR))
 840       return s;
 841     i = i->next;
 842   }
 843   return NULL;
 844 }
 845
 846 static GtsSegment * reverse (GtsSegment * start,
 847                              gboolean interior,
 848                              gboolean * isloop)
 849 {
 850   GtsSegment * s = start, * prev = NULL, * rprev = NULL;
 851   GtsSegment * rstart = NULL, * rstart1 = NULL;
 852
 853   do {
 854     GtsSegment * rs;
 855
 856     g_assert (IS_EDGE_INTER (s));
 857     rs = GTS_SEGMENT (edge_inter_new (s->v2, s->v1,
 858                                       EDGE_INTER (s)->t1, EDGE_INTER (s)->t2));
 859
 860     if (rstart == NULL)
 861       rstart = rs;
 862     else if (rstart1 == NULL)
 863       rstart1 = rs;
 864     if (interior)
 865       SET (rs, INTERIOR);
 866     NEXT (rs) = rprev;
 867     rprev = rs;
 868     prev = s;
 869     s = NEXT (s);
 870   } while (s != NULL && s != start);
 871   if (s == start) {
 872     NEXT (rstart) = rprev;
 873     *isloop = TRUE;
 874   }
 875   else {
 876     NEXT (rstart) = start;
 877     NEXT (prev) = rprev;
 878     *isloop = FALSE;
 879   }
 880   return rstart1;
 881 }
 882
 883 static GSList * interior_loops (GSList * interior)
 884 {
 885   GSList * i = interior;
 886   GSList * loops = NULL;
 887
 888   i = interior;
 889   while (i) {
 890     GtsSegment * s = i->data;
 891
 892     if (IS_SET (s, RELEVANT)) {
 893       GtsSegment * start = s, * end;
 894
 895       do {
 896         GtsSegment * next = next_flag (s, INTERIOR);
 897
 898         UNSET (s, RELEVANT);
 899         end = s;
 900         s = NEXT (s) = next;
 901       } while (s != NULL && s != start);
 902
 903       if (s == start)
 904         loops = g_slist_prepend (loops, start);
 905       else {
 906         GtsSegment * next, * prev;
 907         gboolean isloop;
 908
 909         s = prev_flag (start, INTERIOR);
 910         while (s) {
 911           UNSET (s, RELEVANT);
 912           NEXT (s) = start;
 913           start = s;
 914           s = prev_flag (s, INTERIOR);
 915         }
 916         next = next_flag (end, RELEVANT);
 917         prev = prev_flag (start, RELEVANT);
 918         if (prev != NULL)
 919           SET (start->v1, INTERIOR);
 920         if (next != NULL)
 921           SET (end->v2, INTERIOR);
 922         if (next == NULL && prev == NULL)
 923           loops = g_slist_prepend (loops, start);
 924         else
 925           reverse (start, TRUE, &isloop);
 926       }
 927     }
 928     i = i->next;
 929   }
 930   return loops;
 931 }
 932
 933 #define ORIENTATION(p1,p2,p3,o) (gts_point_orientation_3d (p1, p2, o, p3))
 934 #define ORIENTATION_SOS(p1,p2,p3,o) (gts_point_orientation_3d_sos (p1, p2, o, p3))
 935
 936 #define ORIENTED_VERTICES(s,next,w1,w2) {\
 937   if ((s)->v1 == (next)->v1 || (s)->v1 == (next)->v2) {\
 938     w1 = (s)->v2;\
 939     w2 = (s)->v1;\
 940   }\
 941   else {\
 942     w1 = (s)->v1;\
 943     w2 = (s)->v2;\
 944   }\
 945 }
 946
 947 #if 0
 948 static GtsSegment * segment_intersects (GtsPoint * p1, GtsPoint * p2,
 949                                         GSList * i,
 950                                         GtsPoint * o)
 951 {
 952   while (i) {
 953     GtsSegment * s = i->data;
 954     GtsPoint * p3 = GTS_POINT (s->v1);
 955     GtsPoint * p4 = GTS_POINT (s->v2);
 956
 957     if (p3 != p1 && p3 != p2 && p4 != p1 && p4 != p2) {
 958       gdouble o1 = ORIENTATION (p3, p4, p1, o);
 959       gdouble o2 = ORIENTATION (p3, p4, p2, o);
 960
 961       if ((o1 < 0. && o2 > 0.) || (o1 > 0. && o2 < 0.)) {
 962         o1 = ORIENTATION (p1, p2, p3, o);
 963         o2 = ORIENTATION (p1, p2, p4, o);
 964
 965         if ((o1 <= 0. && o2 >= 0.) || (o1 >= 0. && o2 <= 0.))
 966           return s;
 967       }
 968     }
 969     i = i->next;
 970   }
 971   return NULL;
 972 }
 973 #else
 974 static GtsSegment * segment_intersects (GtsPoint * p1, GtsPoint * p2,
 975                                         GSList * i,
 976                                         GtsPoint * o)
 977 {
 978   while (i) {
 979     GtsSegment * s = i->data;
 980     GtsPoint * p3 = GTS_POINT (s->v1);
 981     GtsPoint * p4 = GTS_POINT (s->v2);
 982
 983     if (p3 != p1 && p3 != p2 && p4 != p1 && p4 != p2) {
 984       gint o1 = ORIENTATION_SOS (p3, p4, p1, o);
 985       gint o2 = ORIENTATION_SOS (p3, p4, p2, o);
 986
 987       if (o1*o2 < 0) {
 988         o1 = ORIENTATION_SOS (p1, p2, p3, o);
 989         o2 = ORIENTATION_SOS (p1, p2, p4, o);
 990
 991         if (o1*o2 < 0)
 992           return s;
 993       }
 994     }
 995     i = i->next;
 996   }
 997   return NULL;
 998 }
 999 #endif
1000
1001 static gboolean is_inside_wedge (GtsSegment * s1, GtsSegment * s2,
1002                                  GtsPoint * p, GtsPoint * o)
1003 {
1004   GtsVertex * v1, * v2, * v3;
1005
1006   ORIENTED_VERTICES (s1, s2, v1, v2);
1007   v3 = s2->v1 != v2 ? s2->v1 : s2->v2;
1008
1009   if (ORIENTATION (GTS_POINT (v1), GTS_POINT (v2),
1010                    GTS_POINT (v3), o) >= 0.) {
1011     if (ORIENTATION (GTS_POINT (v1), GTS_POINT (v2), p, o) <= 0. ||
1012         ORIENTATION (GTS_POINT (v2), GTS_POINT (v3), p, o) <= 0.)
1013       return FALSE;
1014   }
1015   else if (ORIENTATION (GTS_POINT (v1), GTS_POINT (v2), p, o) <= 0. &&
1016            ORIENTATION (GTS_POINT (v2), GTS_POINT (v3), p, o) <= 0.)
1017     return FALSE;
1018   return TRUE;
1019 }
1020
1021 static GtsSegment * connection (GtsPoint * p,
1022                                 GSList * interior,
1023                                 GSList * bloops,
1024                                 GtsPoint * o)
1025 {
1026   while (bloops) {
1027     GtsSegment * start = bloops->data, * s = start;
1028
1029     do {
1030       GtsSegment * next = NEXT (s);
1031       GtsVertex * v2 = s->v1 == next->v1 || s->v1 == next->v2 ? s->v1 : s->v2;
1032
1033       if (is_inside_wedge (s, next, p, o) &&
1034           !segment_intersects (p, GTS_POINT (v2), interior, o))
1035         return s;
1036       s = next;
1037     } while (s != start);
1038     bloops = bloops->next;
1039   }
1040   return NULL;
1041 }
1042
1043 static gdouble loop_orientation (GtsSegment * start,
1044                                  GtsPoint * p, GtsPoint * o)
1045 {
1046   GtsSegment * s = start;
1047   gdouble or = 0.;
1048
1049   do {
1050     GtsSegment * next = NEXT (s);
1051     GtsVertex * v1, * v2;
1052
1053     ORIENTED_VERTICES (s, next, v1, v2);
1054     or += ORIENTATION (p, GTS_POINT (v1), GTS_POINT (v2), o);
1055     s = next;
1056   } while (s != start);
1057
1058 #ifdef DEBUG
1059   fprintf (stderr, "loop orientation: %g\n", or);
1060 #endif /* DEBUG */
1061
1062   return or;
1063 }
1064
1065 static void connect_interior_loop (GtsSegment * start,
1066                                    GSList ** interior,
1067                                    GSList ** bloops,
1068                                    GtsSurface * surface,
1069                                    GtsPoint * o)
1070 {
1071   GtsSegment * s = start, * c = NULL, * next, * s1, * rs1, * rs;
1072   GtsVertex * v, * cv;
1073   gboolean isloop;
1074
1075   do {
1076     if (!(c = connection (GTS_POINT (s->v2), *interior, *bloops, o)))
1077       s = NEXT (s);
1078   } while (s != start && !c);
1079   g_assert (c);
1080   next = NEXT (c);
1081   v = c->v1 == next->v1 || c->v1 == next->v2 ? c->v1 : c->v2;
1082   cv = s->v2;
1083 #ifdef DEBUG
1084   fprintf (stderr, "connecting %p:%s with %p:%s\n",
1085            cv, NAME (cv), v, NAME (v));
1086   fprintf (stderr, "  c: %p: %p:%s %p:%s\n", c,
1087            c->v1, NAME (c->v1),
1088            c->v2, NAME (c->v2));
1089   fprintf (stderr, "  next: %p: %p:%s %p:%s\n", next,
1090            next->v1, NAME (next->v1),
1091            next->v2, NAME (next->v2));
1092 #endif /* DEBUG */
1093   rs = reverse (s, FALSE, &isloop);
1094   if (isloop) {
1095     if (loop_orientation (rs, GTS_POINT (v), o) < 0.) {
1096       GtsSegment * tmp = s;
1097       s = rs;
1098       rs = tmp;
1099     }
1100     *bloops = g_slist_prepend (*bloops, rs);
1101   }
1102   s1 = GTS_SEGMENT (gts_edge_new (surface->edge_class, v, cv));
1103   rs1 = GTS_SEGMENT (gts_edge_new (surface->edge_class, cv, v));
1104   NEXT (c) = s1;
1105   NEXT (rs1) = next;
1106   *interior = g_slist_prepend (*interior, s1);
1107   NEXT (s1) = NEXT (s);
1108   NEXT (s) = rs1;
1109 }
1110
1111 static GSList * boundary_loops (GSList * boundary)
1112 {
1113   GSList * i = boundary;
1114   GtsSegment * start = i->data;
1115   GSList * loops = NULL;
1116
1117   while (i) {
1118     GtsSegment * s = i->data;
1119     GSList * inext = i->next;
1120     GtsSegment * next = inext ? inext->data : start;
1121     GtsVertex * v = s->v1 == next->v1 || s->v1 == next->v2 ? s->v1 : s->v2;
1122
1123     if (IS_SET (v, INTERIOR)) {
1124       GtsSegment * intprev = prev_interior (v);
1125
1126       NEXT (intprev) = next;
1127       NEXT (s) = next_interior (v);
1128       UNSET (v, INTERIOR);
1129     }
1130     else
1131       NEXT (s) = next;
1132     i = inext;
1133   }
1134
1135   i = boundary;
1136   while (i) {
1137     start = i->data;
1138
1139     if (IS_SET (start, RELEVANT)) {
1140       GtsSegment * s = start;
1141
1142       do {
1143         UNSET (s, RELEVANT);
1144         UNSET (s, INTERIOR);
1145         s = NEXT (s);
1146       } while (s != start);
1147       loops = g_slist_prepend (loops, start);
1148     }
1149     i = i->next;
1150   }
1151
1152   return loops;
1153 }
1154
1155 typedef struct _Ear    Ear;
1156
1157 struct _Ear {
1158   GtsVertex * v1, * v2, * v3;
1159   GtsSegment * s1, * s2, * s3;
1160 };
1161
1162 static gboolean point_in_wedge (GtsPoint * p1, GtsPoint * p2, GtsPoint * p3,
1163                                 GtsPoint * p, gboolean closed, GtsPoint * o)
1164 {
1165   gdouble o1;
1166
1167   if (p == p2 || p == p3)
1168     return FALSE;
1169   o1 = ORIENTATION (p1, p2, p, o);
1170   if ((closed && o1 < 0.) || (!closed && o1 <= 0.)) return FALSE;
1171   o1 = ORIENTATION (p3, p1, p, o);
1172   if ((closed && o1 < 0.) || (!closed && o1 <= 0.)) return FALSE;
1173   return TRUE;
1174 }
1175
1176 #if 0
1177 static gboolean segment_intersects1 (GtsPoint * p1, GtsPoint * p2,
1178                                      GtsPoint * p3, GtsPoint * p4,
1179                                      gboolean closed, GtsPoint * o)
1180 {
1181   gdouble o1 = ORIENTATION (p3, p4, p1, o);
1182   gdouble o2 = ORIENTATION (p3, p4, p2, o);
1183   gdouble o3, o4;
1184
1185   if ((closed && ((o1 > 0. && o2 > 0.) || (o1 < 0. && o2 < 0.))) ||
1186       (!closed && ((o1 >= 0. && o2 >= 0.) || (o1 <= 0. && o2 <= 0.))))
1187     return FALSE;
1188   o3 = ORIENTATION (p1, p2, p3, o);
1189   o4 = ORIENTATION (p1, p2, p4, o);
1190   if ((o3 > 0. && o4 > 0.) || (o3 < 0. && o4 < 0.))
1191     return FALSE;
1192   if (closed) return TRUE;
1193   if ((o3 == 0. && o4 > 0.) || (o4 == 0. && o3 > 0.))
1194     return TRUE;
1195   return FALSE;
1196 }
1197 #else
1198 static gboolean segment_intersects1 (GtsPoint * p1, GtsPoint * p2,
1199                                      GtsPoint * p3, GtsPoint * p4,
1200                                      gboolean closed, GtsPoint * o)
1201 {
1202   gint o1, o2;
1203
1204   o1 = ORIENTATION_SOS (p3, p4, p1, o);
1205   o2 = ORIENTATION_SOS (p3, p4, p2, o);
1206   if (o1*o2 > 0)
1207     return FALSE;
1208   o1 = ORIENTATION_SOS (p1, p2, p3, o);
1209   o2 = ORIENTATION_SOS (p1, p2, p4, o);
1210   if (o1*o2 > 0)
1211     return FALSE;
1212   return TRUE;
1213 }
1214 #endif
1215
1216 static GtsSegment * triangle_intersects_segments (GtsPoint * p1,
1217                                                   GtsPoint * p2,
1218                                                   GtsPoint * p3,
1219                                                   gboolean closed,
1220                                                   GtsSegment * start,
1221                                                   GtsPoint * o)
1222 {
1223   GtsSegment * s = start;
1224
1225   do {
1226     GtsPoint * p4 = GTS_POINT (s->v1);
1227     GtsPoint * p5 = GTS_POINT (s->v2);
1228
1229     if (p4 == p1) {
1230       if (point_in_wedge (p1, p2, p3, p5, closed, o))
1231         return s;
1232     }
1233     else if (p4 == p2) {
1234       if (point_in_wedge (p2, p3, p1, p5, closed, o))
1235         return s;
1236     }
1237     else if (p4 == p3) {
1238       if (point_in_wedge (p3, p1, p2, p5, closed, o))
1239         return s;
1240     }
1241     else if (p5 == p1) {
1242       if (point_in_wedge (p1, p2, p3, p4, closed, o))
1243         return s;
1244     }
1245     else if (p5 == p2) {
1246       if (point_in_wedge (p2, p3, p1, p4, closed, o))
1247         return s;
1248     }
1249     else if (p5 == p3) {
1250       if (point_in_wedge (p3, p1, p2, p4, closed, o))
1251         return s;
1252     }
1253     else if (segment_intersects1 (p1, p2, p4, p5, closed, o) ||
1254              segment_intersects1 (p2, p3, p4, p5, closed, o) ||
1255              segment_intersects1 (p3, p1, p4, p5, closed, o))
1256       return s;
1257     s = NEXT (s);
1258   } while (s != start);
1259   return NULL;
1260 }
1261
1262 static gboolean new_ear (GtsSegment * s,
1263                          Ear * e,
1264                          GtsSegment * start,
1265                          guint sloppy,
1266                          GtsPoint * o)
1267 {
1268   gdouble or;
1269
1270   e->s1 = s;
1271   e->s2 = NEXT (s);
1272
1273   g_return_val_if_fail (e->s2, FALSE);
1274   g_return_val_if_fail (e->s2 != e->s1, FALSE);
1275
1276   ORIENTED_VERTICES (e->s1, e->s2, e->v1, e->v2);
1277   e->v3 = e->s2->v1 != e->v2 ? e->s2->v1 : e->s2->v2;
1278   if (e->v3 == e->v1)
1279     return FALSE;
1280   e->s3 = NEXT (e->s2);
1281   if (gts_segment_connect (e->s3, e->v1, e->v3)) {
1282     if (NEXT (e->s3) != e->s1)
1283       return FALSE;
1284   }
1285   else if (gts_vertices_are_connected (e->v1, e->v3))
1286     return FALSE;
1287   else
1288     e->s3 = NULL;
1289   or = ORIENTATION (GTS_POINT (e->v1), GTS_POINT (e->v2), GTS_POINT (e->v3),o);
1290   switch (sloppy) {
1291   case 0:
1292     if (or <= 0. ||
1293         triangle_intersects_segments (GTS_POINT (e->v1), GTS_POINT (e->v2),
1294                                       GTS_POINT (e->v3), TRUE, start, o))
1295       return FALSE;
1296     break;
1297   case 1:
1298     if (or < 0. ||
1299         (or > 0. &&
1300          triangle_intersects_segments (GTS_POINT (e->v1), GTS_POINT (e->v2),
1301                                        GTS_POINT (e->v3), FALSE, start, o)))
1302       return FALSE;
1303     break;
1304   case 2:
1305     if ((or > 0. &&
1306          triangle_intersects_segments (GTS_POINT (e->v1), GTS_POINT (e->v2),
1307                                        GTS_POINT (e->v3), FALSE, start, o)) ||
1308         (or < 0. &&
1309          triangle_intersects_segments (GTS_POINT (e->v2), GTS_POINT (e->v1),
1310                                        GTS_POINT (e->v3), FALSE, start, o)))
1311       return FALSE;
1312     break;
1313   case 3:
1314     if (or < 0.)
1315       return FALSE;
1316     break;
1317   }
1318 #ifdef DEBUG
1319   if (or <= 0.)
1320     fprintf (stderr, "or: %g\n", or);
1321 #endif /* DEBUG */
1322   g_assert (or > -1e-6);
1323   return TRUE;
1324 }
1325
1326 static void triangulate_loop (GtsSegment * start,
1327                               GtsSurface * surface,
1328                               GtsPoint * o)
1329 {
1330   GtsSegment * prev = start, * s;
1331   guint sloppy = 0;
1332 #ifdef DEBUG
1333   guint nt = 0;
1334 #endif /* DEBUG */
1335
1336   s = NEXT (start);
1337   while (NEXT (s) != s) {
1338     GtsSegment * next = NEXT (s);
1339     Ear e;
1340
1341 #ifdef DEBUG
1342     fprintf (stderr, "prev: %p s: %p next: %p\n", prev, s, next);
1343 #endif /* DEBUG */
1344
1345     if (!new_ear (s, &e, start, sloppy, o)) {
1346       if (s == start) {
1347         sloppy++;
1348 #ifdef DEBUG
1349         fprintf (stderr, "sloppy: %u\n", sloppy);
1350 #endif /* DEBUG */
1351       }
1352       prev = s;
1353       s = next;
1354     }
1355     else {
1356       GtsFace * f;
1357
1358       if (!GTS_IS_EDGE (e.s3))
1359         e.s3 = GTS_SEGMENT (gts_edge_new (surface->edge_class, e.v1, e.v3));
1360       f = gts_face_new (surface->face_class,
1361                         GTS_EDGE (e.s1), GTS_EDGE (e.s2), GTS_EDGE (e.s3));
1362       gts_surface_add_face (surface, f);
1363       UNSET (e.s1, RELEVANT);
1364       UNSET (e.s1, INTERIOR);
1365       UNSET (e.s2, RELEVANT);
1366       UNSET (e.s2, INTERIOR);
1367       NEXT (prev) = e.s3;
1368       NEXT (e.s3) = NEXT (e.s2);
1369       NEXT (e.s1) = NEXT (e.s2) = NULL;
1370       start = prev;
1371       s = NEXT (prev);
1372       sloppy = 0;
1373 #ifdef DEBUG
1374       {
1375         gchar name[80];
1376         FILE * fp;
1377
1378         fprintf (stderr, " t.%u: (%p:%s,%p:%s,%p:%s)\n",
1379                  nt,
1380                  e.v1, NAME (e.v1),
1381                  e.v2, NAME (e.v2),
1382                  e.v3, NAME (e.v3));
1383         sprintf (name, "/tmp/t.%u", nt++);
1384         fp = fopen (name, "wt");
1385         //      gts_surface_write (surface, fp);
1386         gts_write_triangle (GTS_TRIANGLE (f), NULL, fp);
1387         //        write_graph1 (start, interior, surface, fp);
1388         fclose (fp);
1389         print_loop (start, stderr);
1390       }
1391 #endif /* DEBUG */
1392     }
1393   }
1394   UNSET (s, RELEVANT);
1395   UNSET (s, INTERIOR);
1396   NEXT (s) = NULL;
1397 }
1398
1399 #ifdef CHECK_ORIENTED
1400 static void check_object (GtsObject * o)
1401 {
1402   g_assert (o->reserved == NULL);
1403   g_assert (o->flags == 0);
1404 }
1405
1406 static void check_boundary (GtsEdge * e, GtsSurface * s)
1407 {
1408   check_object (GTS_OBJECT (e));
1409   check_object (GTS_OBJECT (GTS_SEGMENT (e)->v1));
1410   check_object (GTS_OBJECT (GTS_SEGMENT (e)->v2));
1411   g_assert (gts_edge_face_number (e, s) == 1);
1412 }
1413
1414 static void check_interior (GtsEdge * e, GtsSurface * s)
1415 {
1416   guint n;
1417   check_object (GTS_OBJECT (e));
1418   check_object (GTS_OBJECT (GTS_SEGMENT (e)->v1));
1419   check_object (GTS_OBJECT (GTS_SEGMENT (e)->v2));
1420
1421   n = gts_edge_face_number (e, s);
1422 #ifdef DEBUG
1423   if (n != 2)
1424     gts_surface_print_stats (s, stderr);
1425 #endif /* DEBUG */
1426   g_assert (n == 2);
1427 }
1428
1429 static void check_boundary_interior_triangulation (GSList * boundary,
1430                                                    GSList * interior,
1431                                                    GtsSurface * surface)
1432 {
1433   g_slist_foreach (boundary, (GFunc) check_boundary, surface);
1434   g_slist_foreach (interior, (GFunc) check_interior, surface);
1435 }
1436 #endif /*ifdef CHECK_ORIENTED */
1437
1438 static void merge_duplicate (GtsEdge * e)
1439 {
1440   GtsEdge * dup = gts_edge_is_duplicate (e);
1441
1442   g_assert (dup);
1443   gts_edge_replace (dup, e);
1444   gts_object_destroy (GTS_OBJECT (dup));
1445 }
1446
1447 static void triangulate_boundary_interior (GSList * boundary,
1448                                            GSList * interior,
1449                                            GtsSurface * s,
1450                                            GtsPoint * o)
1451 {
1452   GSList * iloops, * bloops, * i;
1453
1454   i = boundary;
1455   while (i) {
1456     SET (i->data, RELEVANT);
1457     i = i->next;
1458   }
1459   i = interior;
1460   while (i) {
1461     SET (i->data, RELEVANT);
1462     SET (i->data, INTERIOR);
1463     i = i->next;
1464   }
1465
1466   iloops = interior_loops (interior);
1467   bloops = boundary_loops (boundary);
1468
1469   i = iloops;
1470   while (i) {
1471 #ifdef DEBUG
1472     fprintf (stderr, "--- interior loop ---\n");
1473     print_loop (i->data, stderr);
1474 #endif /* DEBUG */
1475     connect_interior_loop (i->data, &interior, &bloops, s, o);
1476     i = i->next;
1477   }
1478
1479 #ifdef DEBUG
1480  {
1481    FILE * fp = fopen ("/tmp/bloops", "w");
1482    write_loops (bloops, fp);
1483    fclose (fp);
1484  }
1485 #endif /* DEBUG */
1486
1487   i = bloops;
1488   while (i) {
1489 #ifdef DEBUG
1490     fprintf (stderr, "--- boundary loop ---\n");
1491     print_loop (i->data, stderr);
1492 #endif /* DEBUG */
1493     triangulate_loop (i->data, s, o);
1494     i = i->next;
1495   }
1496
1497   g_slist_foreach (interior, (GFunc) merge_duplicate, NULL);
1498   g_slist_free (iloops);
1499   g_slist_free (bloops);
1500
1501 #ifdef CHECK_ORIENTED
1502   check_boundary_interior_triangulation (boundary, interior, s);
1503 #endif /* CHECK_ORIENTED */
1504 }
1505
1506 static void create_edges (GtsSegment * s, GtsSurface * surface)
1507 {
1508   if (GTS_OBJECT (s)->reserved) {
1509     GList * i = GTS_OBJECT (s)->reserved;
1510     GtsVertex * v1 = i->data;
1511
1512     GTS_OBJECT (s)->reserved = g_list_prepend (i,
1513                       gts_edge_new (surface->edge_class, s->v1, v1));
1514     while (i) {
1515       GList * next = i->next;
1516       GtsVertex * v2 = next ? next->data : s->v2;
1517
1518       GTS_OBJECT (i->data)->reserved = NULL;
1519       i->data = gts_edge_new (surface->edge_class, v1, v2);
1520       v1 = v2;
1521       i = next;
1522     }
1523   }
1524 }
1525
1526 static void add_boundary (GtsSegment * s, GtsSegment * next,
1527                           GSList ** boundary)
1528 {
1529   if (GTS_OBJECT (s)->reserved == NULL)
1530     *boundary = g_slist_prepend (*boundary, s);
1531   else {
1532     if (s->v2 == next->v2 || s->v2 == next->v1) {
1533       GList * i = g_list_last (GTS_OBJECT (s)->reserved);
1534
1535       while (i) {
1536         *boundary = g_slist_prepend (*boundary, i->data);
1537         i = i->prev;
1538       }
1539     }
1540     else {
1541       GList * i = GTS_OBJECT (s)->reserved;
1542
1543       while (i) {
1544         *boundary = g_slist_prepend (*boundary, i->data);
1545         i = i->next;
1546       }
1547     }
1548   }
1549 }
1550
1551 static void triangulate_face (GtsTriangle * t, GtsSurface * surface)
1552 {
1553   GSList * interior = GTS_OBJECT (t)->reserved;
1554   GSList * boundary = NULL;
1555   GtsSurface * s = gts_surface_new (gts_surface_class (),
1556                                     surface->face_class,
1557                                     surface->edge_class,
1558                                     surface->vertex_class);
1559   gdouble x, y, z;
1560   GtsPoint * p = GTS_POINT (GTS_SEGMENT (t->e1)->v1);
1561   GtsPoint * o;
1562
1563   GTS_OBJECT (t)->reserved = NULL;
1564   gts_triangle_normal (t, &x, &y, &z);
1565   g_assert (x != 0. || y != 0. || z != 0.);
1566   o = gts_point_new (gts_point_class (), p->x + x, p->y + y, p->z + z);
1567   add_boundary (GTS_SEGMENT (t->e3), GTS_SEGMENT (t->e1), &boundary);
1568   add_boundary (GTS_SEGMENT (t->e2), GTS_SEGMENT (t->e3), &boundary);
1569   add_boundary (GTS_SEGMENT (t->e1), GTS_SEGMENT (t->e2), &boundary);
1570 #ifdef DEBUG
1571   {
1572     static guint nt = 0;
1573     char name[80];
1574     FILE * fp;
1575
1576     fprintf (stderr, "%u: triangulating %p\n", nt, t);
1577 if (nt == 28)
1578   fprintf (stderr, "tintin!!!!\n");
1579     sprintf (name, "/tmp/oc.%u", nt++);
1580     fp = fopen (name, "wt");
1581     //    write_graph (boundary, interior, s, fp);
1582     write_segments (boundary, interior, fp);
1583     fclose (fp);
1584   }
1585 #endif /* DEBUG */
1586   triangulate_boundary_interior (boundary, interior, s, o);
1587   g_slist_free (interior);
1588   g_slist_free (boundary);
1589   if (GTS_OBJECT (t)->klass->attributes)
1590     gts_surface_foreach_face (s, (GtsFunc) gts_object_attributes, t);
1591   gts_surface_merge (surface, s);
1592   gts_object_destroy (GTS_OBJECT (s));
1593   gts_object_destroy (GTS_OBJECT (o));
1594 }
1595
1596 static void free_edge_list (GtsObject * o)
1597 {
1598   g_list_free (o->reserved);
1599   o->reserved = NULL;
1600 }
1601
1602 /**
1603  * gts_surface_inter_new:
1604  * @klass: a #GtsSurfaceInterClass.
1605  * @s1: a #GtsSurface.
1606  * @s2: a #GtsSurface.
1607  * @faces_tree1: a bounding box tree (see gts_bb_tree_new()) for
1608  * the faces of @s1.
1609  * @faces_tree2: a bounding box tree for the faces of @s2.
1610  * @is_open1: whether @s1 is an "open" surface.
1611  * @is_open2: whether @s2 is an "open" surface.
1612  *
1613  * When triangulating the cut faces, the new faces inherit the
1614  * attributes of these original faces through their attributes()
1615  * method.
1616  *
1617  * Returns: a new #GtsSurfaceInter describing the intersection of @s1
1618  * and @s2.
1619  */
1620 GtsSurfaceInter * gts_surface_inter_new (GtsSurfaceInterClass * klass,
1621                                          GtsSurface * s1,
1622                                          GtsSurface * s2,
1623                                          GNode * faces_tree1,
1624                                          GNode * faces_tree2,
1625                                          gboolean is_open1,
1626                                          gboolean is_open2)
1627 {
1628   GtsSurfaceInter * si;
1629   GtsSurface * s;
1630
1631   g_return_val_if_fail (klass != NULL, NULL);
1632   g_return_val_if_fail (s1 != NULL, NULL);
1633   g_return_val_if_fail (s2 != NULL, NULL);
1634   g_return_val_if_fail (faces_tree1 != NULL, NULL);
1635   g_return_val_if_fail (faces_tree2 != NULL, NULL);
1636
1637   si = surface_inter_new (klass, s1, s2, faces_tree1, faces_tree2);
1638
1639   gts_surface_foreach_edge (si->s1, (GtsFunc) create_edges, si->s1);
1640   gts_surface_foreach_edge (si->s2, (GtsFunc) create_edges, si->s2);
1641
1642 #ifdef DEBUG
1643   fprintf (stderr, "====== triangulating s1 ======\n");
1644 #endif /* DEBUG */
1645   s = gts_surface_new (gts_surface_class (),
1646                        s1->face_class,
1647                        s1->edge_class,
1648                        s1->vertex_class);
1649   gts_surface_foreach_face (si->s1, (GtsFunc) triangulate_face, s);
1650   gts_surface_foreach_edge (si->s1, (GtsFunc) free_edge_list, NULL);
1651   gts_object_destroy (GTS_OBJECT (si->s1));
1652   si->s1 = s;
1653   GTS_OBJECT (si->s1)->reserved = s1;
1654
1655 #ifdef DEBUG
1656   fprintf (stderr, "====== triangulating s2 ======\n");
1657 #endif /* DEBUG */
1658   s = gts_surface_new (gts_surface_class (),
1659                        s2->face_class,
1660                        s2->edge_class,
1661                        s2->vertex_class);
1662   gts_surface_foreach_face (si->s2, (GtsFunc) triangulate_face, s);
1663   gts_surface_foreach_edge (si->s2, (GtsFunc) free_edge_list, NULL);
1664   gts_object_destroy (GTS_OBJECT (si->s2));
1665   si->s2 = s;
1666   GTS_OBJECT (si->s2)->reserved = s2;
1667
1668   return si;
1669 }
1670
1671 static void check_surface_edge (GtsEdge * e, gpointer * data)
1672 {
1673   gboolean * ok = data[0];
1674   GtsSurface * s = data[1];
1675   GtsSurface * bs = GTS_OBJECT (s)->reserved;
1676   guint nf = gts_edge_face_number (e, s);
1677
1678   if (nf < 1 || nf > 2) {
1679     *ok = FALSE;
1680     g_return_if_fail (nf >= 1 && nf <= 2);
1681   }
1682   if (nf == 1 && gts_edge_face_number (e, bs) == 0) {
1683     *ok = FALSE;
1684     g_return_if_fail (gts_edge_face_number (e, bs) > 0);
1685   }
1686 }
1687
1688 static void mark_edge (GtsObject * o, gpointer data)
1689 {
1690   o->reserved = data;
1691 }
1692
1693 static gint triangle_orientation (GtsTriangle * t, GtsEdge * e)
1694 {
1695   GtsSegment * s = GTS_SEGMENT (t->e1 == e ? t->e2
1696                                 :
1697                                 t->e2 == e ? t->e3
1698                                 :
1699                                 t->e1);
1700   GtsVertex * v2 = GTS_SEGMENT (e)->v2;
1701
1702   if (s->v1 == v2 || s->v2 == v2)
1703     return 1;
1704   return -1;
1705 }
1706
1707 static gboolean check_orientation (GtsEdge * e, GtsSurface * s)
1708 {
1709   GtsTriangle * t1 = NULL, * t2 = NULL;
1710   GSList * i = e->triangles;
1711   gint o1 = 0, o2 = 0;
1712
1713   while (i) {
1714     if (GTS_IS_FACE (i->data) &&
1715         gts_face_has_parent_surface (i->data, s)) {
1716       if (t1 == NULL) {
1717         t1 = i->data;
1718         o1 = triangle_orientation (t1, e);
1719       }
1720       else if (t2 == NULL) {
1721         t2 = i->data;
1722         o2 = triangle_orientation (t2, e);
1723         g_return_val_if_fail (o1*o2 < 0, FALSE);
1724       }
1725       else
1726         g_assert_not_reached ();
1727     }
1728     i = i->next;
1729   }
1730   g_return_val_if_fail (t1 && t2, FALSE);
1731   return TRUE;
1732 }
1733
1734 static void check_edge (GtsSegment * s, gpointer * data)
1735 {
1736   gboolean * ok = data[0];
1737   GtsSurfaceInter * si = data[1];
1738   gboolean * closed = data[2];
1739   GSList * j;
1740   guint nn = 0;
1741
1742   j = s->v1->segments;
1743   while (j && *ok) {
1744     GtsSegment * s1 = j->data;
1745
1746     if (s1 != s && GTS_OBJECT (s1)->reserved == si) {
1747       if (s1->v2 != s->v1)
1748         *ok = FALSE;
1749       nn++;
1750     }
1751     j = j->next;
1752   }
1753   j = s->v2->segments;
1754   while (j && *ok) {
1755     GtsSegment * s1 = j->data;
1756
1757     if (s1 != s && GTS_OBJECT (s1)->reserved == si) {
1758       if (s1->v1 != s->v2)
1759         *ok = FALSE;
1760       nn++;
1761     }
1762     j = j->next;
1763   }
1764   if (nn != 2)
1765     *closed = FALSE;
1766
1767   if (!check_orientation (GTS_EDGE (s), si->s1))
1768     *ok = FALSE;
1769   if (!check_orientation (GTS_EDGE (s), si->s2))
1770     *ok = FALSE;
1771 }
1772
1773 /**
1774  * gts_surface_inter_check:
1775  * @si: a #GtsSurfaceInter.
1776  * @closed: is set to %TRUE if @si->edges is a closed curve, %FALSE
1777  * otherwise.
1778  *
1779  * Returns: %TRUE if the curve described by @si is an orientable
1780  * manifold, %FALSE otherwise.
1781  */
1782 gboolean gts_surface_inter_check (GtsSurfaceInter * si,
1783                                   gboolean * closed)
1784 {
1785   gboolean ok = TRUE;
1786   gpointer data[3];
1787
1788   g_return_val_if_fail (si != NULL, FALSE);
1789   g_return_val_if_fail (closed != NULL, FALSE);
1790
1791   *closed = si->edges ? TRUE : FALSE;
1792
1793   /* mark edges as used by si */
1794   g_slist_foreach (si->edges, (GFunc) mark_edge, si);
1795
1796   data[0] = &ok;
1797   data[1] = si;
1798   data[2] = closed;
1799   g_slist_foreach (si->edges, (GFunc) check_edge, data);
1800   g_slist_foreach (si->edges, (GFunc) gts_object_reset_reserved, NULL);
1801
1802   /* check connectivity of the faces of @si */
1803   if (*closed) {
1804     gpointer data[2];
1805
1806     data[0] = &ok;
1807     data[1] = si->s1;
1808     gts_surface_foreach_edge (si->s1, (GtsFunc) check_surface_edge, data);
1809     data[1] = si->s2;
1810     gts_surface_foreach_edge (si->s2, (GtsFunc) check_surface_edge, data);
1811   }
1812
1813   return ok;
1814 }
1815
1816 /* Given @e and @f returns a #GtsFace compatible with @f and belonging to
1817    @s1 or @s2 */
1818 static GtsFace * next_compatible_face (GtsEdge * e,
1819                                        GtsFace * f,
1820                                        GtsSurface * s1,
1821                                        GtsSurface * s2)
1822 {
1823   GSList * i = e->triangles;
1824   GtsFace * f2 = NULL, * f3 = NULL;
1825
1826   while (i) {
1827     GtsFace * f1 = i->data;
1828
1829     if (f1 != f && GTS_IS_FACE (f1)) {
1830       if (gts_face_has_parent_surface (f1, s1))
1831         return f1;
1832       if (gts_face_has_parent_surface (f1, s2)) {
1833         if (f2 == NULL) f2 = f1;
1834         else if (f3 == NULL) f3 = f1;
1835         else g_assert_not_reached (); /* s2 is a non-manifold surface */
1836       }
1837     }
1838     i = i->next;
1839   }
1840   if (f3 == NULL) {
1841     if (gts_edge_is_boundary (e, s2))
1842       return NULL;
1843     return f2;
1844   }
1845   g_assert (gts_face_has_parent_surface (f, s1));
1846   if (gts_triangles_are_compatible (GTS_TRIANGLE (f), GTS_TRIANGLE (f2), e))
1847     return f2;
1848   return f3;
1849 }
1850
1851 static void walk_faces (GtsEdge * e, GtsFace * f,
1852                         GtsSurface * s1,
1853                         GtsSurface * s2,
1854                         GtsSurface * s)
1855 {
1856   GtsFifo * faces = gts_fifo_new ();
1857   GtsFifo * edges = gts_fifo_new ();
1858
1859   gts_fifo_push (faces, f);
1860   gts_fifo_push (edges, e);
1861   while ((f = gts_fifo_pop (faces)) && (e = gts_fifo_pop (edges))) {
1862     if (!GTS_OBJECT (f)->reserved) {
1863       GtsTriangle * t = GTS_TRIANGLE (f);
1864       GtsFace * f1;
1865
1866       gts_surface_add_face (s, f);
1867       GTS_OBJECT (f)->reserved = s;
1868       if (t->e1 != e && !GTS_OBJECT (t->e1)->reserved &&
1869           (f1 = next_compatible_face (t->e1, f, s1, s2))) {
1870         gts_fifo_push (faces, f1);
1871         gts_fifo_push (edges, t->e1);
1872       }
1873       if (t->e2 != e && !GTS_OBJECT (t->e2)->reserved &&
1874           (f1 = next_compatible_face (t->e2, f, s1, s2))) {
1875         gts_fifo_push (faces, f1);
1876         gts_fifo_push (edges, t->e2);
1877       }
1878       if (t->e3 != e && !GTS_OBJECT (t->e3)->reserved &&
1879           (f1 = next_compatible_face (t->e3, f, s1, s2))) {
1880         gts_fifo_push (faces, f1);
1881         gts_fifo_push (edges, t->e3);
1882       }
1883     }
1884   }
1885   gts_fifo_destroy (faces);
1886   gts_fifo_destroy (edges);
1887 }
1888
1889 /**
1890  * gts_surface_inter_boolean:
1891  * @si: a #GtsSurfaceInter.
1892  * @surface: a #GtsSurface.
1893  * @op: a #GtsBooleanOperation.
1894  *
1895  * Adds to @surface the part of the surface described by @si and @op.
1896  */
1897 void gts_surface_inter_boolean (GtsSurfaceInter * si,
1898                                 GtsSurface * surface,
1899                                 GtsBooleanOperation op)
1900 {
1901   GtsSurface * s = NULL;
1902   gint orient = 1;
1903   GSList * i;
1904
1905   g_return_if_fail (si != NULL);
1906   g_return_if_fail (surface != NULL);
1907
1908   switch (op) {
1909   case GTS_1_OUT_2: s = si->s1; orient = 1; break;
1910   case GTS_1_IN_2: s = si->s1; orient = -1; break;
1911   case GTS_2_OUT_1: s = si->s2; orient = -1; break;
1912   case GTS_2_IN_1: s = si->s2; orient = 1; break;
1913   default: g_assert_not_reached ();
1914   }
1915
1916   /* mark edges as belonging to intersection */
1917   g_slist_foreach (si->edges, (GFunc) mark_edge, si);
1918
1919   i = si->edges;
1920   while (i) {
1921     GtsEdge * e = i->data;
1922     GSList * j = e->triangles;
1923
1924     while (j) {
1925       if (gts_face_has_parent_surface (j->data, s) &&
1926           orient*triangle_orientation (j->data, e) > 0) {
1927 #ifdef DEBUG_BOOLEAN
1928         GtsFace * boundary = gts_edge_is_boundary (e, surface);
1929
1930         g_assert (!boundary || boundary == j->data);
1931 #endif /* DEBUG_BOOLEAN */
1932         walk_faces (e, j->data, s, GTS_OBJECT (s)->reserved, surface);
1933         break;
1934       }
1935       j = j->next;
1936     }
1937     i = i->next;
1938   }
1939   g_slist_foreach (si->edges, (GFunc) gts_object_reset_reserved, NULL);
1940   gts_surface_foreach_face (surface,
1941                             (GtsFunc) gts_object_reset_reserved, NULL);
1942 }
1943
1944 static void self_intersecting (GtsBBox * bb1, GtsBBox * bb2,
1945                                gpointer * d)
1946 {
1947   GtsTriangle * t1 = bb1->bounded;
1948   GtsTriangle * t2 = bb2->bounded;
1949
1950   if (t1 != t2) {
1951     GtsSegment * s1 = GTS_SEGMENT (t1->e1);
1952     GtsSegment * s2 = GTS_SEGMENT (t1->e2);
1953     GtsSegment * s3 = GTS_SEGMENT (t1->e3);
1954     GtsSegment * s4 = GTS_SEGMENT (t2->e1);
1955     GtsSegment * s5 = GTS_SEGMENT (t2->e2);
1956     GtsSegment * s6 = GTS_SEGMENT (t2->e3);
1957     GtsPoint * pi;
1958
1959     if ((!gts_segments_touch (s4, s1) &&
1960          !gts_segments_touch (s4, s2) &&
1961          !gts_segments_touch (s4, s3) &&
1962          (pi = segment_triangle_intersection (s4, t1, gts_point_class ()))
1963          != NULL) ||
1964         (!gts_segments_touch (s5, s1) &&
1965          !gts_segments_touch (s5, s2) &&
1966          !gts_segments_touch (s5, s3) &&
1967          (pi = segment_triangle_intersection (s5, t1, gts_point_class ()))
1968          != NULL) ||
1969         (!gts_segments_touch (s6, s1) &&
1970          !gts_segments_touch (s6, s2) &&
1971          !gts_segments_touch (s6, s3) &&
1972          (pi = segment_triangle_intersection (s6, t1, gts_point_class ()))
1973          != NULL)) {
1974       GtsBBTreeTraverseFunc func = d[0];
1975       gpointer data = d[1];
1976       gboolean * self_inter = d[2];
1977
1978       gts_object_destroy (GTS_OBJECT (pi));
1979       *self_inter = TRUE;
1980       (* func) (bb1, bb2, data);
1981     }
1982   }
1983 }
1984
1985 /**
1986  * gts_surface_foreach_intersecting_face:
1987  * @s: a #GtsSurface.
1988  * @func: a #GtsBBTreeTraverseFunc.
1989  * @data: user data to pass to @func.
1990  *
1991  * Calls @func for each intersecting pair of faces of @s.
1992  *
1993  * Returns: %TRUE if @func was called at least once, %FALSE otherwise.
1994  */
1995 gboolean gts_surface_foreach_intersecting_face (GtsSurface * s,
1996                                                 GtsBBTreeTraverseFunc func,
1997                                                 gpointer data)
1998 {
1999   GNode * tree;
2000   gpointer d[3];
2001   gboolean self_inter = FALSE;
2002
2003   g_return_val_if_fail (s != NULL, FALSE);
2004   g_return_val_if_fail (func != NULL, FALSE);
2005
2006   tree = gts_bb_tree_surface (s);
2007   d[0] = func;
2008   d[1] = data;
2009   d[2] = &self_inter;
2010   gts_bb_tree_traverse_overlapping (tree, tree,
2011                                     (GtsBBTreeTraverseFunc) self_intersecting,
2012                                     d);
2013   gts_bb_tree_destroy (tree, TRUE);
2014
2015   return self_inter;
2016 }
2017
2018 static void add_intersecting (GtsBBox * bb1, GtsBBox * bb2,
2019                               GtsSurface * intersected)
2020 {
2021   gts_surface_add_face (intersected, bb1->bounded);
2022   gts_surface_add_face (intersected, bb2->bounded);
2023 }
2024
2025 /**
2026  * gts_surface_is_self_intersecting:
2027  * @s: a #GtsSurface.
2028  *
2029  * Returns: a new #GtsSurface containing the faces of @s which are
2030  * self-intersecting or %NULL if no faces of @s are self-intersecting.
2031  */
2032 GtsSurface * gts_surface_is_self_intersecting (GtsSurface * s)
2033 {
2034   GtsSurface * intersected;
2035
2036   g_return_val_if_fail (s != NULL, NULL);
2037
2038   intersected = gts_surface_new (GTS_SURFACE_CLASS (GTS_OBJECT (s)->klass),
2039                                  s->face_class,
2040                                  s->edge_class,
2041                                  s->vertex_class);
2042   if (!gts_surface_foreach_intersecting_face (s,
2043                       (GtsBBTreeTraverseFunc) add_intersecting, intersected)) {
2044     gts_object_destroy (GTS_OBJECT (intersected));
2045     intersected = NULL;
2046   }
2047   return intersected;
2048 }