boolean.c

Go to the documentation of this file.

/* GTS - Library for the manipulation of triangulated surfaces
 * Copyright (C) 1999--2002 Stéphane Popinet
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Library General Public
 * License as published by the Free Software Foundation; either
 * version 2 of the License, or (at your option) any later version.
 *
 * This library is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * Library General Public License for more details.
 *
 * You should have received a copy of the GNU Library General Public
 * License along with this library; if not, write to the
 * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
 * Boston, MA 02111-1307, USA.
 */

#include <math.h>
#include "gts.h"

/*#define DEBUG*/
/*#define DEBUG_BOOLEAN*/
/*#define CHECK_ORIENTED*/

#ifdef DEBUG
#ifdef DEBUG_FUNCTIONS
#include "gts-private.h"
#else
#define DEBUG_FUNCTIONS
#include "gts-private.h"
#undef DEBUG_FUNCTIONS
#endif
#endif

static void surface_inter_destroy (GtsObject * object)
{
  GtsSurfaceInter * si = GTS_SURFACE_INTER (object);

  gts_object_destroy (GTS_OBJECT (si->s1));
  gts_object_destroy (GTS_OBJECT (si->s2));
  g_slist_free (si->edges);

  (* GTS_OBJECT_CLASS (gts_surface_inter_class ())->parent_class->destroy)
    (object);
}

static void surface_inter_class_init (GtsObjectClass * klass)
{
  klass->destroy = surface_inter_destroy;
}

static void surface_inter_init (GtsSurfaceInter * si)
{
  si->s1 = si->s2 = NULL;
  si->edges = NULL;
}

GtsSurfaceInterClass * gts_surface_inter_class (void)
{
  static GtsSurfaceInterClass * klass = NULL;

  if (klass == NULL) {
    GtsObjectClassInfo surface_inter_info = {
      "GtsSurfaceInter",
      sizeof (GtsSurfaceInter),
      sizeof (GtsSurfaceInterClass),
      (GtsObjectClassInitFunc) surface_inter_class_init,
      (GtsObjectInitFunc) surface_inter_init,
      (GtsArgSetFunc) NULL,
      (GtsArgGetFunc) NULL
    };
    klass = gts_object_class_new (gts_object_class (), &surface_inter_info);
  }

  return klass;
}

/* EdgeInter: Header */

typedef struct _EdgeInter         EdgeInter;

struct _EdgeInter {
  GtsEdge parent;

  GtsTriangle * t1, * t2;
};

#define EDGE_INTER(obj)            GTS_OBJECT_CAST (obj,\
                                                 EdgeInter,\
                                                 edge_inter_class ())
#define IS_EDGE_INTER(obj)         (gts_object_is_from_class (obj,\
                                                 edge_inter_class ()))

static GtsEdgeClass * edge_inter_class  (void);
static EdgeInter * edge_inter_new    (GtsVertex * v1, GtsVertex * v2,
                                      GtsTriangle * t1, GtsTriangle * t2);

/* EdgeInter: Object */

static GtsEdgeClass * edge_inter_class (void)
{
  static GtsEdgeClass * klass = NULL;

  if (klass == NULL) {
    GtsObjectClassInfo edge_inter_info = {
      "EdgeInter",
      sizeof (EdgeInter),
      sizeof (GtsEdgeClass),
      (GtsObjectClassInitFunc) NULL,
      (GtsObjectInitFunc) NULL,
      (GtsArgSetFunc) NULL,
      (GtsArgGetFunc) NULL
    };
    klass = gts_object_class_new (GTS_OBJECT_CLASS (gts_constraint_class ()),
                                  &edge_inter_info);
  }

  return klass;
}

static EdgeInter * edge_inter_new (GtsVertex * v1, GtsVertex * v2,
                                   GtsTriangle * t1, GtsTriangle * t2)
{
  EdgeInter * object;

  object = EDGE_INTER (gts_edge_new (GTS_EDGE_CLASS (edge_inter_class ()), 
                                     v1, v2));
  object->t1 = t1;
  object->t2 = t2;

  return object;
}

static GtsPoint * segment_triangle_intersection (GtsSegment * s,
                                                 GtsTriangle * t,
                                                 GtsPointClass * klass)
{
  GtsPoint * A, * B, * C, * D, * E;
  gint ABCE, ABCD, ADCE, ABDE, BCDE;
  GtsEdge * AB, * BC, * CA;
  gdouble a, b, c;

  g_return_val_if_fail (s != NULL, NULL);
  g_return_val_if_fail (t != NULL, NULL);
  g_return_val_if_fail (klass != NULL, NULL);

  gts_triangle_vertices_edges (t, NULL, 
                               (GtsVertex **) &A, 
                               (GtsVertex **) &B, 
                               (GtsVertex **) &C, 
                               &AB, &BC, &CA);
  D = GTS_POINT (s->v1);
  E = GTS_POINT (s->v2);

  ABCE = gts_point_orientation_3d_sos (A, B, C, E);
  ABCD = gts_point_orientation_3d_sos (A, B, C, D);
  if (ABCE < 0 || ABCD > 0) {
    GtsPoint * tmpp;
    gint tmp;

    tmpp = E; E = D; D = tmpp;
    tmp = ABCE; ABCE = ABCD; ABCD = tmp;
  }
  if (ABCE < 0 || ABCD > 0)
    return NULL;
  ADCE = gts_point_orientation_3d_sos (A, D, C, E);
  if (ADCE < 0)
    return NULL;
  ABDE = gts_point_orientation_3d_sos (A, B, D, E);
  if (ABDE < 0)
    return NULL;
  BCDE = gts_point_orientation_3d_sos (B, C, D, E);
  if (BCDE < 0)
    return NULL;
  a = gts_point_orientation_3d (A, B, C, E);
  b = gts_point_orientation_3d (A, B, C, D);
  if (a != b) {
    c = a/(a - b);
    return gts_point_new (klass,
                          E->x + c*(D->x - E->x),
                          E->y + c*(D->y - E->y),
                          E->z + c*(D->z - E->z));
  }
  /* D and E are contained within ABC */
#ifdef DEBUG
  fprintf (stderr, 
           "segment: %p:%s triangle: %p:%s intersection\n"
           "D and E contained in ABC\n",
           s, GTS_NEDGE (s)->name, t, GTS_NFACE (t)->name);
#endif /* DEBUG */  
  g_assert (a == 0.); 
  return gts_point_new (klass,
                        (E->x + D->x)/2.,
                        (E->y + D->y)/2.,
                        (E->z + D->z)/2.);
}

static gint triangle_triangle_orientation (GtsPoint * p1, 
                                           GtsPoint * p2, GtsPoint * p3,
                                           GtsPoint * p4, GtsPoint * p5,
                                           GtsPoint * p6)
{
  gint o4 = 0, o5 = 0, o6 = 0;

  if (p4 != p1 && p4 != p2 && p4 != p3)
    o4 = gts_point_orientation_3d_sos (p1, p2, p3, p4);
  if (p5 != p1 && p5 != p2 && p5 != p3)
    o5 = gts_point_orientation_3d_sos (p1, p2, p3, p5);
  if (o4*o5 < 0)
    return 0;
  if (p6 != p1 && p6 != p2 && p6 != p3)
    o6 = gts_point_orientation_3d_sos (p1, p2, p3, p6);
  if (o4*o6 < 0 || o5*o6 < 0)
    return 0;
  if (o4) return o4;
  if (o5) return o5;
  g_assert (o6);
  return o6;
}

static gint triangle_point_orientation (GtsTriangle * t1, 
                                        GtsTriangle * t2,
                                        gint o1,
                                        GtsPoint * p)
{
  GtsPoint * p1 = GTS_POINT (GTS_SEGMENT (t1->e1)->v1);
  GtsPoint * p2 = GTS_POINT (GTS_SEGMENT (t1->e1)->v2);
  GtsPoint * p3 = GTS_POINT (gts_triangle_vertex (t1));
  GtsPoint * p4 = GTS_POINT (GTS_SEGMENT (t2->e1)->v1);
  GtsPoint * p5 = GTS_POINT (GTS_SEGMENT (t2->e1)->v2);
  GtsPoint * p6 = GTS_POINT (gts_triangle_vertex (t2));
  gint o = triangle_triangle_orientation (p1, p2, p3, p4, p5, p6);

  if (o != 0)
    return o;
  o = triangle_triangle_orientation (p4, p5, p6, p1, p2, p3);
  if (o != 0) {
    gint o2 = gts_point_orientation_3d_sos (p4, p5, p6, p);

    return - o*o1*o2;
  }
  return 0;
}

static void add_edge_inter (GtsEdge * e,
                            GtsTriangle * t,
                            GtsVertex * v)
{
  GtsVertex * ev1 = GTS_SEGMENT (e)->v1, * ev2 = GTS_SEGMENT (e)->v2;
  GList * i = GTS_OBJECT (e)->reserved;

  GTS_OBJECT (v)->reserved = t;
  if (i == NULL) {
    GTS_OBJECT (e)->reserved = g_list_prepend (NULL, v);
#ifdef DEBUG
    fprintf (stderr, "add_edge_inter: inserting %p (%p,%p)\n", v, e, t);
#endif /* DEBUG */
  }
  else {
    GtsPoint * p1 = GTS_POINT (GTS_SEGMENT (t->e1)->v1);
    GtsPoint * p2 = GTS_POINT (GTS_SEGMENT (t->e1)->v2);
    GtsPoint * p3 = GTS_POINT (gts_triangle_vertex (t));
    gint o1, oref = gts_point_orientation_3d_sos (p1, p2, p3, GTS_POINT (ev1));
    
    o1 = oref;
    while (i) {
      gint o2 = triangle_point_orientation (t, GTS_OBJECT (i->data)->reserved,
                                            oref, GTS_POINT (ev1));

      if (o2 == 0) {
#ifdef DEBUG
        g_warning ("add_edge_inter: safe sign evaluation failed\n");
#endif /* DEBUG */
        o2 = gts_point_orientation_3d_sos (p1, p2, p3, i->data);
      }

      if (o1*o2 < 0)
        break;
      ev1 = i->data;
      o1 = o2;
      i = i->next;
    }
    if (i != NULL) {
      GList * n = g_list_prepend (NULL, v);

      ev2 = i->data;
      n->next = i;
      n->prev = i->prev;
      i->prev = n;
      if (n->prev == NULL)
        GTS_OBJECT (e)->reserved = n;
      else
        n->prev->next = n;
    }
    else {
      g_assert (o1*gts_point_orientation_3d_sos (p1, p2, p3, GTS_POINT (ev2))
                < 0);
      GTS_OBJECT (e)->reserved = g_list_append (GTS_OBJECT (e)->reserved, v);
    }
#ifdef DEBUG
    fprintf (stderr, 
             "add_edge_inter: inserting %p (%p,%p) between %p and %p\n", 
             v, e, t, ev1, ev2);
    i = GTS_OBJECT (e)->reserved;
    while (i) {
      fprintf (stderr, " %p", i->data);
      i = i->next;
    }
    fprintf (stderr, "\n");
#endif /* DEBUG */
  }
}

static GtsVertex * intersects (GtsEdge * e,
                               GtsTriangle * t,
                               GtsSurface * s)
{
  GList * i = GTS_OBJECT (e)->reserved;
  GtsVertex * v;

  while (i) {
    if (GTS_OBJECT (i->data)->reserved == t)
      return i->data;
    i = i->next;
  }

  v = GTS_VERTEX (segment_triangle_intersection (GTS_SEGMENT (e), t, 
                                    GTS_POINT_CLASS (s->vertex_class)));
  if (v != NULL) {
#ifdef DEBUG
    if (GTS_IS_NVERTEX (v) && GTS_IS_NEDGE (e) && GTS_IS_NFACE (t) &&
        GTS_NVERTEX (v)->name[0] == '\0')
      g_snprintf (GTS_NVERTEX (v)->name, GTS_NAME_LENGTH, "%s|%s",
                  GTS_NEDGE (e)->name, GTS_NFACE (t)->name);
#endif /* DEBUG */
    if (s->vertex_class->intersection_attributes)
      (*s->vertex_class->intersection_attributes)
        (v, GTS_OBJECT (e), GTS_OBJECT (t));
    add_edge_inter (e, t, v);
  }
  return v;
}

/* see figure misc/orientation.fig */
static gint intersection_orientation (GtsTriangle * t1, 
                                      GtsEdge * e,
                                      GtsTriangle * t2)
{
  GtsVertex * v1, * v2, * v3;
  GtsEdge * e2, * e3;
  GtsVertex * v4, * v5, * v6;

  gts_triangle_vertices_edges (t1, e, &v1, &v2, &v3, &e, &e2, &e3);
  gts_triangle_vertices (t2, &v4, &v5, &v6);

  return gts_point_orientation_3d_sos (GTS_POINT (v4), 
                                       GTS_POINT (v5), 
                                       GTS_POINT (v6),
                                       GTS_POINT (v2));
}

#define UPDATE_ORIENTATION if (o > 0) { vi2 = v; /* e2 = e; */ } else { vi2 = vi1;\
                                                                        /* e2 = e1; */\
                                                                        vi1 = v;\
                                                                        /* e1 = e; */ }

static void intersect_edges (GtsBBox * bb1, GtsBBox * bb2,
                             GtsSurfaceInter * si)
{
  GtsSurface * s1 = GTS_OBJECT (si->s1)->reserved;
  GtsTriangle * t1 = GTS_TRIANGLE (bb1->bounded);
  GtsTriangle * t2 = GTS_TRIANGLE (bb2->bounded);
  GtsVertex * v, * vi1 = NULL, * vi2 = NULL;
  //GtsEdge * e1 = NULL, * e2 = NULL, * e;

  vi1 = intersects (t2->e1, t1, s1);
  //e1 = t2->e1;
  v = intersects (t2->e2, t1, s1);
  //e = t2->e2;
  if (!vi1) {
    vi1 = v;
    //e1 = e;
  }
  else if (v) {
    gint o = intersection_orientation (t2, t2->e2, t1);
    UPDATE_ORIENTATION;
  }
  if (!vi2) {
    v = intersects (t2->e3, t1, s1);
    //e = t2->e3;
    if (!vi1) {
      vi1 = v;
      //e1 = e;
    }
    else if (v) {
      gint o = intersection_orientation (t2, t2->e3, t1);
      UPDATE_ORIENTATION;
    }
  }
  if (!vi2) {
    v = intersects (t1->e1, t2, s1);
    //e = t1->e1;
    if (!vi1) {
      vi1 = v;
      //e1 = e;
    }
    else if (v) {
      gint o = - intersection_orientation (t1, t1->e1, t2);
      UPDATE_ORIENTATION;
    }
  }
  if (!vi2) {
    v = intersects (t1->e2, t2, s1);
    //e = t1->e2;
    if (!vi1) {
      vi1 = v;
      //e1 = e;
    }
    else if (v) {
      gint o = - intersection_orientation (t1, t1->e2, t2);
      UPDATE_ORIENTATION;
    }
  }
  if (!vi2) {
    v = intersects (t1->e3, t2, s1);
    //e = t1->e3;
    if (!vi1) {
      vi1 = v;
      //e1 = e;
    }
    else if (v) {
      gint o = - intersection_orientation (t1, t1->e3, t2);
      UPDATE_ORIENTATION;
    }
  }

  g_assert ((!vi1 && !vi2) || (vi1 && vi2));
  if (vi1) {
    GtsEdge * e = GTS_EDGE (edge_inter_new (vi1, vi2, t1, t2));

#ifdef DEBUG
    fprintf (stderr, "creating constraint %p: %p->%p: %p/%p\n", 
             e, vi1, vi2, t1, t2);
#endif /* DEBUG */
    gts_surface_add_face (si->s1, GTS_FACE (t1));
    gts_surface_add_face (si->s2, GTS_FACE (t2));
    si->edges = g_slist_prepend (si->edges, e);
    GTS_OBJECT (t1)->reserved = g_slist_prepend (GTS_OBJECT (t1)->reserved, e);
    GTS_OBJECT (t2)->reserved = g_slist_prepend (GTS_OBJECT (t2)->reserved, e);
  }
}

static GtsSurfaceInter * surface_inter_new (GtsSurfaceInterClass * klass,
                                            GtsSurface * s1,
                                            GtsSurface * s2,
                                            GNode * faces_tree1,
                                            GNode * faces_tree2)
{
  GtsSurfaceInter * si;

  si = GTS_SURFACE_INTER (gts_object_new (GTS_OBJECT_CLASS (klass)));
  si->s1 = gts_surface_new (gts_surface_class (),
                            s1->face_class,
                            s1->edge_class,
                            s1->vertex_class);
  GTS_OBJECT (si->s1)->reserved = s1;
  si->s2 = gts_surface_new (gts_surface_class (),
                            s2->face_class,
                            s2->edge_class,
                            s2->vertex_class);
  GTS_OBJECT (si->s2)->reserved = s2;
  gts_bb_tree_traverse_overlapping (faces_tree1, faces_tree2,
                                    (GtsBBTreeTraverseFunc) intersect_edges, 
                                    si);

  return si;
}

static void free_slist (GtsObject * o)
{
  g_slist_free (o->reserved);
  o->reserved = NULL;
}

static void free_glist (GtsObject * o)
{
  g_list_foreach (o->reserved, (GFunc) gts_object_reset_reserved, NULL);
  g_list_free (o->reserved);
  o->reserved = NULL;
}

GSList * gts_surface_intersection (GtsSurface * s1,
                                   GtsSurface * s2,
                                   GNode * faces_tree1,
                                   GNode * faces_tree2)
{
  GtsSurfaceInter * si;
  GSList * inter;

  g_return_val_if_fail (s1 != NULL, NULL);
  g_return_val_if_fail (s2 != NULL, NULL);
  g_return_val_if_fail (faces_tree1 != NULL, NULL);
  g_return_val_if_fail (faces_tree2 != NULL, NULL);

  si = surface_inter_new (gts_surface_inter_class (),
                          s1, s2, faces_tree1, faces_tree2);

  gts_surface_foreach_face (si->s1, (GtsFunc) free_slist, NULL);
  gts_surface_foreach_face (si->s2, (GtsFunc) free_slist, NULL);
  gts_surface_foreach_edge (si->s1, (GtsFunc) free_glist, NULL);
  gts_surface_foreach_edge (si->s2, (GtsFunc) free_glist, NULL);
  inter = si->edges;
  si->edges = NULL;
  gts_object_destroy (GTS_OBJECT (si));

  return inter;  
}

typedef enum {
  INTERIOR = 1 << (GTS_USER_FLAG),
  RELEVANT = 1 << (GTS_USER_FLAG + 1)
} CurveFlag;

#define IS_SET(s, f) ((GTS_OBJECT_FLAGS (s) & (f)) != 0)
#define SET(s, f)   (GTS_OBJECT_FLAGS (s) |= (f))
#define UNSET(s, f) (GTS_OBJECT_FLAGS (s) &= ~(f))
#define NEXT(s)  (GTS_OBJECT (s)->reserved)

#ifdef DEBUG
static void write_nodes (GSList * i, GHashTable * hash, guint * nn,
                         FILE * fp)
{
  while (i) {
    GtsSegment * s = i->data;

    if (!g_hash_table_lookup (hash, s->v1)) {
      fprintf (fp, "  %u [ label = \"%p\" ];\n", *nn, s->v1);
      g_hash_table_insert (hash, s->v1, GUINT_TO_POINTER ((*nn)++));
    }
    if (!g_hash_table_lookup (hash, s->v2)) {
      fprintf (fp, "  %u [ label = \"%p\" ];\n", *nn, s->v2);
      g_hash_table_insert (hash, s->v2, GUINT_TO_POINTER ((*nn)++));
    }
    i = i->next;
  }
}

static void write_edges (GSList * i, GHashTable * hash, 
                         GtsSurface * surface,
                         FILE * fp)
{
  while (i) {
    GtsSegment * s = i->data;

    fprintf (fp, "  %u -> %u [ label = \"%p:%d\" ];\n",
             GPOINTER_TO_UINT (g_hash_table_lookup (hash, s->v1)),
             GPOINTER_TO_UINT (g_hash_table_lookup (hash, s->v2)),
             s,
             gts_edge_face_number (GTS_EDGE (s), surface));
    i = i->next;
  }
}

static void write_graph (GSList * boundary, GSList * interior,
                         GtsSurface * surface,
                         FILE * fp)
{
  GHashTable * hash = g_hash_table_new (NULL, NULL);
  guint nn = 1;
  
  fprintf (fp, "digraph oriented_curve {\n");
  write_nodes (boundary, hash, &nn, fp);
  write_nodes (interior, hash, &nn, fp);
  write_edges (boundary, hash, surface, fp);
  fprintf (fp, "  edge [ color = red ];\n");
  write_edges (interior, hash, surface, fp);
  fprintf (fp, "}\n");
  g_hash_table_destroy (hash);
}

static void print_loop (GtsSegment * start, FILE * fp)
{
  GtsSegment * s = start;

  do {
    fprintf (fp, "  %p: %p:%s -> %p:%s\n",
             s, 
             s->v1, GTS_NVERTEX (s->v1)->name, 
             s->v2, GTS_NVERTEX (s->v2)->name);
    s = NEXT (s);
  } while (s != start && s != NULL);
}

static void draw_vector (GtsPoint * p1, GtsPoint * p2, FILE * fp)
{
  gdouble x = p2->x - p1->x;
  gdouble y = p2->y - p1->y;
  gdouble z = p2->z - p1->z;

  fprintf (fp, "VECT 1 3 0 3 0 %g %g %g %g %g %g %g %g %g\n",
           p1->x + x - (x - y/2.)/5.,
           p1->y + y - (x/2. + y)/5.,
           p1->z + z - (x/2. + z)/5.,
           p1->x + x,
           p1->y + y,
           p1->z + z,
           p1->x + x - (x + y/2.)/5.,
           p1->y + y + (x/2. - y)/5.,
           p1->z + z + (x/2. - z)/5.);
  fprintf (fp, "VECT 1 2 0 2 0 %g %g %g %g %g %g\n",
           p1->x, p1->y, p1->z,
           p1->x + x,
           p1->y + y,
           p1->z + z);
}

static void draw_vector1 (GtsPoint * p1, GtsPoint * p2, GtsPoint * o,
                          FILE * fp)
{
  gdouble x1 = o->x + 0.9*(p1->x - o->x);
  gdouble y1 = o->y + 0.9*(p1->y - o->y);
  gdouble z1 = o->z + 0.9*(p1->z - o->z);
  gdouble x2 = o->x + 0.9*(p2->x - o->x);
  gdouble y2 = o->y + 0.9*(p2->y - o->y);
  gdouble z2 = o->z + 0.9*(p2->z - o->z);
  gdouble x = x2 - x1;
  gdouble y = y2 - y1;
  gdouble z = z2 - z1;

  fprintf (fp, "VECT 1 3 0 3 0 %g %g %g %g %g %g %g %g %g\n",
           x1 + x - (x - y/2.)/5.,
           y1 + y - (x/2. + y)/5.,
           z1 + z - (x/2. + z)/5.,
           x1 + x,
           y1 + y,
           z1 + z,
           x1 + x - (x + y/2.)/5.,
           y1 + y + (x/2. - y)/5.,
           z1 + z + (x/2. - z)/5.);
  fprintf (fp, "VECT 1 2 0 2 0 %g %g %g %g %g %g\n",
           x1, y1, z1,
           x1 + x,
           y1 + y,
           z1 + z);
}

static void write_segments (GSList * boundary, GSList * interior,
                            FILE * fp)
{
  GSList * i = boundary;

  fprintf (fp, "LIST {\n");
  while (i) {
    GSList * inext = i->next;
    GtsSegment * s = i->data;
    GtsSegment * next = inext ? inext->data : boundary->data;
    GtsVertex * v1, * v2;

    if (s->v1 != next->v1 && s->v1 != next->v2) {
      v1 = s->v1;
      v2 = s->v2;
    }
    else {
      v1 = s->v2;
      v2 = s->v1;
    }
    draw_vector (GTS_POINT (v1), GTS_POINT (v2), fp);
    i = inext;
  }
  i = interior;
  while (i) {
    GtsSegment * s = i->data;

    draw_vector (GTS_POINT (s->v1), GTS_POINT (s->v2), fp);
    i = i->next;
  }
  fprintf (fp, "}\n");
}

static void write_loops (GSList * i, FILE * fp)
{
  guint nl = 0;

  while (i) {
    GtsSegment * start = i->data, * s;
    GtsPoint os;
    guint n = 0;

    fprintf (fp, "(geometry \"loop%d\" = LIST {\n", nl++);    

    os.x = os.y = os.z = 0.;
    s = start;
    do {
      GtsSegment * next = NEXT (s);
      GtsPoint * p;
      
      if (s->v1 != next->v1 && s->v1 != next->v2)
        p = GTS_POINT (s->v1);
       else
         p = GTS_POINT (s->v2);
      os.x += p->x; os.y += p->y; os.z += p->z; n++;
      s = next;
     } while (s != start);
    os.x /= n; os.y /= n; os.z /= n;
    
    s = start;
    do {
      GtsSegment * next = NEXT (s);
      
      if (s->v1 != next->v1 && s->v1 != next->v2)
        draw_vector1 (GTS_POINT (s->v1), GTS_POINT (s->v2), &os, fp);
      else
         draw_vector1 (GTS_POINT (s->v2), GTS_POINT (s->v1), &os, fp);
      s = next;
    } while (s != start);
    
    fprintf (fp, "})\n");

    i = i->next;
  }
}

#define NAME(v) (GTS_IS_NVERTEX (v) ? GTS_NVERTEX (v)->name : "")
#endif /* DEBUG */

static GtsSegment * prev_flag (GtsSegment * s, CurveFlag flag)
{
  GSList * i = s->v1->segments;

  while (i) {
    if (i->data != s && IS_SET (i->data, flag))
      return i->data;
    i = i->next;
  }
  return NULL;
}

static GtsSegment * next_flag (GtsSegment * s, CurveFlag flag)
{
  GSList * i = s->v2->segments;

  while (i) {
    if (i->data != s && IS_SET (i->data, flag))
      return i->data;
    i = i->next;
  }
  return NULL;
}

static GtsSegment * next_interior (GtsVertex * v)
{
  GSList * i = v->segments;

  while (i) {
    GtsSegment * s = i->data;

    if (s->v1 == v && IS_SET (s, INTERIOR))
      return s;
    i = i->next;
  }
  return NULL;
}

static GtsSegment * prev_interior (GtsVertex * v)
{
  GSList * i = v->segments;

  while (i) {
    GtsSegment * s = i->data;

    if (s->v2 == v && IS_SET (s, INTERIOR))
      return s;
    i = i->next;
  }
  return NULL;
}

static GtsSegment * reverse (GtsSegment * start,
                             gboolean interior,
                             gboolean * isloop)
{
  GtsSegment * s = start, * prev = NULL, * rprev = NULL;
  GtsSegment * rstart = NULL, * rstart1 = NULL;

  do {
    GtsSegment * rs;

    g_assert (IS_EDGE_INTER (s));
    rs = GTS_SEGMENT (edge_inter_new (s->v2, s->v1,
                                      EDGE_INTER (s)->t1, EDGE_INTER (s)->t2));

    if (rstart == NULL)
      rstart = rs;
    else if (rstart1 == NULL)
      rstart1 = rs;
    if (interior)
      SET (rs, INTERIOR);
    NEXT (rs) = rprev;
    rprev = rs;
    prev = s;
    s = NEXT (s);
  } while (s != NULL && s != start);
  if (s == start) {
    NEXT (rstart) = rprev;
    *isloop = TRUE;
  }
  else {
    NEXT (rstart) = start;
    NEXT (prev) = rprev;
    *isloop = FALSE;
  }    
  return rstart1;
}

static GSList * interior_loops (GSList * interior)
{
  GSList * i = interior;
  GSList * loops = NULL;

  i = interior;
  while (i) {
    GtsSegment * s = i->data;

    if (IS_SET (s, RELEVANT)) {
      GtsSegment * start = s, * end;

      do {
        GtsSegment * next = next_flag (s, INTERIOR);

        UNSET (s, RELEVANT);
        end = s; 
        s = NEXT (s) = next;
      } while (s != NULL && s != start);

      if (s == start)
        loops = g_slist_prepend (loops, start);
      else {
        GtsSegment * next, * prev;
        gboolean isloop;

        s = prev_flag (start, INTERIOR);
        while (s) {
          UNSET (s, RELEVANT);
          NEXT (s) = start;
          start = s;
          s = prev_flag (s, INTERIOR);
        }
        next = next_flag (end, RELEVANT);
        prev = prev_flag (start, RELEVANT);
        if (prev != NULL)
          SET (start->v1, INTERIOR);
        if (next != NULL)
          SET (end->v2, INTERIOR);
        if (next == NULL && prev == NULL)
          loops = g_slist_prepend (loops, start);
        else
          reverse (start, TRUE, &isloop);
      }
    }
    i = i->next;
  }
  return loops;
}

#define ORIENTATION(p1,p2,p3,o) (gts_point_orientation_3d (p1, p2, o, p3))
#define ORIENTATION_SOS(p1,p2,p3,o) (gts_point_orientation_3d_sos (p1, p2, o, p3))

#define ORIENTED_VERTICES(s,next,w1,w2) {\
  if ((s)->v1 == (next)->v1 || (s)->v1 == (next)->v2) {\
    w1 = (s)->v2;\
    w2 = (s)->v1;\
  }\
  else {\
    w1 = (s)->v1;\
    w2 = (s)->v2;\
  }\
}

#if 0
static GtsSegment * segment_intersects (GtsPoint * p1, GtsPoint * p2,
                                        GSList * i,
                                        GtsPoint * o)
{
  while (i) {
    GtsSegment * s = i->data;
    GtsPoint * p3 = GTS_POINT (s->v1);
    GtsPoint * p4 = GTS_POINT (s->v2);

    if (p3 != p1 && p3 != p2 && p4 != p1 && p4 != p2) {
      gdouble o1 = ORIENTATION (p3, p4, p1, o);
      gdouble o2 = ORIENTATION (p3, p4, p2, o);

      if ((o1 < 0. && o2 > 0.) || (o1 > 0. && o2 < 0.)) {
        o1 = ORIENTATION (p1, p2, p3, o);
        o2 = ORIENTATION (p1, p2, p4, o);

        if ((o1 <= 0. && o2 >= 0.) || (o1 >= 0. && o2 <= 0.))
          return s;
      }
    }
    i = i->next;
  }
  return NULL;
}
#else
static GtsSegment * segment_intersects (GtsPoint * p1, GtsPoint * p2,
                                        GSList * i,
                                        GtsPoint * o)
{
  while (i) {
    GtsSegment * s = i->data;
    GtsPoint * p3 = GTS_POINT (s->v1);
    GtsPoint * p4 = GTS_POINT (s->v2);

    if (p3 != p1 && p3 != p2 && p4 != p1 && p4 != p2) {
      gint o1 = ORIENTATION_SOS (p3, p4, p1, o);
      gint o2 = ORIENTATION_SOS (p3, p4, p2, o);

      if (o1*o2 < 0) {
        o1 = ORIENTATION_SOS (p1, p2, p3, o);
        o2 = ORIENTATION_SOS (p1, p2, p4, o);

        if (o1*o2 < 0)
          return s;
      }
    }
    i = i->next;
  }
  return NULL;
}
#endif

static gboolean is_inside_wedge (GtsSegment * s1, GtsSegment * s2,
                                 GtsPoint * p, GtsPoint * o)
{
  GtsVertex * v1, * v2, * v3;

  ORIENTED_VERTICES (s1, s2, v1, v2);
  v3 = s2->v1 != v2 ? s2->v1 : s2->v2;

  if (ORIENTATION (GTS_POINT (v1), GTS_POINT (v2), 
                   GTS_POINT (v3), o) >= 0.) {
    if (ORIENTATION (GTS_POINT (v1), GTS_POINT (v2), p, o) <= 0. ||
        ORIENTATION (GTS_POINT (v2), GTS_POINT (v3), p, o) <= 0.)
      return FALSE;
  }
  else if (ORIENTATION (GTS_POINT (v1), GTS_POINT (v2), p, o) <= 0. &&
           ORIENTATION (GTS_POINT (v2), GTS_POINT (v3), p, o) <= 0.)
    return FALSE;
  return TRUE;
}

static GtsSegment * connection (GtsPoint * p, 
                                GSList * interior,
                                GSList * bloops,
                                GtsPoint * o)
{
  while (bloops) {
    GtsSegment * start = bloops->data, * s = start;

    do {
      GtsSegment * next = NEXT (s);
      GtsVertex * v2 = s->v1 == next->v1 || s->v1 == next->v2 ? s->v1 : s->v2;

      if (is_inside_wedge (s, next, p, o) &&
          !segment_intersects (p, GTS_POINT (v2), interior, o))
        return s;
      s = next;
    } while (s != start);
    bloops = bloops->next;
  }
  return NULL;
}

static gdouble loop_orientation (GtsSegment * start,
                                 GtsPoint * p, GtsPoint * o)
{
  GtsSegment * s = start;
  gdouble or = 0.;

  do {
    GtsSegment * next = NEXT (s);
    GtsVertex * v1, * v2;

    ORIENTED_VERTICES (s, next, v1, v2);
    or += ORIENTATION (p, GTS_POINT (v1), GTS_POINT (v2), o);
    s = next;
  } while (s != start);

#ifdef DEBUG
  fprintf (stderr, "loop orientation: %g\n", or);
#endif /* DEBUG */

  return or;
}

static void connect_interior_loop (GtsSegment * start,
                                   GSList ** interior,
                                   GSList ** bloops,
                                   GtsSurface * surface,
                                   GtsPoint * o)
{
  GtsSegment * s = start, * c = NULL, * next, * s1, * rs1, * rs;
  GtsVertex * v, * cv;
  gboolean isloop;

  do {
    if (!(c = connection (GTS_POINT (s->v2), *interior, *bloops, o)))
      s = NEXT (s);
  } while (s != start && !c);
  g_assert (c);
  next = NEXT (c);
  v = c->v1 == next->v1 || c->v1 == next->v2 ? c->v1 : c->v2;
  cv = s->v2;
#ifdef DEBUG
  fprintf (stderr, "connecting %p:%s with %p:%s\n", 
           cv, NAME (cv), v, NAME (v));
  fprintf (stderr, "  c: %p: %p:%s %p:%s\n", c, 
           c->v1, NAME (c->v1),
           c->v2, NAME (c->v2));
  fprintf (stderr, "  next: %p: %p:%s %p:%s\n", next,
           next->v1, NAME (next->v1),
           next->v2, NAME (next->v2));
#endif /* DEBUG */
  rs = reverse (s, FALSE, &isloop);
  if (isloop) {
    if (loop_orientation (rs, GTS_POINT (v), o) < 0.) {
      GtsSegment * tmp = s;
      s = rs;
      rs = tmp;
    }
    *bloops = g_slist_prepend (*bloops, rs);
  }
  s1 = GTS_SEGMENT (gts_edge_new (surface->edge_class, v, cv));
  rs1 = GTS_SEGMENT (gts_edge_new (surface->edge_class, cv, v));
  NEXT (c) = s1;
  NEXT (rs1) = next;
  *interior = g_slist_prepend (*interior, s1);
  NEXT (s1) = NEXT (s);
  NEXT (s) = rs1;
}

static GSList * boundary_loops (GSList * boundary)
{
  GSList * i = boundary;  
  GtsSegment * start = i->data;
  GSList * loops = NULL;

  while (i) {
    GtsSegment * s = i->data;
    GSList * inext = i->next;
    GtsSegment * next = inext ? inext->data : start;
    GtsVertex * v = s->v1 == next->v1 || s->v1 == next->v2 ? s->v1 : s->v2;

    if (IS_SET (v, INTERIOR)) {
      GtsSegment * intprev = prev_interior (v);

      NEXT (intprev) = next;
      NEXT (s) = next_interior (v);
      UNSET (v, INTERIOR);
    }
    else
      NEXT (s) = next;
    i = inext;
  }

  i = boundary;
  while (i) {
    start = i->data;
    
    if (IS_SET (start, RELEVANT)) {
      GtsSegment * s = start;

      do {
        UNSET (s, RELEVANT);
        UNSET (s, INTERIOR);
        s = NEXT (s);
      } while (s != start);
      loops = g_slist_prepend (loops, start);
    }
    i = i->next;
  }

  return loops;
}

typedef struct _Ear    Ear;

struct _Ear {
  GtsVertex * v1, * v2, * v3;
  GtsSegment * s1, * s2, * s3;
};

static gboolean point_in_wedge (GtsPoint * p1, GtsPoint * p2, GtsPoint * p3,
                                GtsPoint * p, gboolean closed, GtsPoint * o)
{
  gdouble o1;

  if (p == p2 || p == p3)
    return FALSE;
  o1 = ORIENTATION (p1, p2, p, o);
  if ((closed && o1 < 0.) || (!closed && o1 <= 0.)) return FALSE;
  o1 = ORIENTATION (p3, p1, p, o);
  if ((closed && o1 < 0.) || (!closed && o1 <= 0.)) return FALSE;
  return TRUE;
}

#if 0
static gboolean segment_intersects1 (GtsPoint * p1, GtsPoint * p2, 
                                     GtsPoint * p3, GtsPoint * p4,
                                     gboolean closed, GtsPoint * o)
{
  gdouble o1 = ORIENTATION (p3, p4, p1, o);
  gdouble o2 = ORIENTATION (p3, p4, p2, o);
  gdouble o3, o4;

  if ((closed && ((o1 > 0. && o2 > 0.) || (o1 < 0. && o2 < 0.))) ||
      (!closed && ((o1 >= 0. && o2 >= 0.) || (o1 <= 0. && o2 <= 0.))))
    return FALSE;
  o3 = ORIENTATION (p1, p2, p3, o);
  o4 = ORIENTATION (p1, p2, p4, o);
  if ((o3 > 0. && o4 > 0.) || (o3 < 0. && o4 < 0.))
    return FALSE;
  if (closed) return TRUE;
  if ((o3 == 0. && o4 > 0.) || (o4 == 0. && o3 > 0.))
    return TRUE;
  return FALSE;
}
#else
static gboolean segment_intersects1 (GtsPoint * p1, GtsPoint * p2, 
                                     GtsPoint * p3, GtsPoint * p4,
                                     gboolean closed, GtsPoint * o)
{
  gint o1, o2;

  o1 = ORIENTATION_SOS (p3, p4, p1, o);
  o2 = ORIENTATION_SOS (p3, p4, p2, o);
  if (o1*o2 > 0)
    return FALSE;
  o1 = ORIENTATION_SOS (p1, p2, p3, o);
  o2 = ORIENTATION_SOS (p1, p2, p4, o);
  if (o1*o2 > 0)
    return FALSE;
  return TRUE;
}
#endif

static GtsSegment * triangle_intersects_segments (GtsPoint * p1,
                                                  GtsPoint * p2,
                                                  GtsPoint * p3,
                                                  gboolean closed,
                                                  GtsSegment * start,
                                                  GtsPoint * o)
{
  GtsSegment * s = start;

  do {
    GtsPoint * p4 = GTS_POINT (s->v1);
    GtsPoint * p5 = GTS_POINT (s->v2);

    if (p4 == p1) {
      if (point_in_wedge (p1, p2, p3, p5, closed, o))
        return s;
    }
    else if (p4 == p2) {
      if (point_in_wedge (p2, p3, p1, p5, closed, o))
        return s;
    }
    else if (p4 == p3) {
      if (point_in_wedge (p3, p1, p2, p5, closed, o))
        return s;
    }
    else if (p5 == p1) {
      if (point_in_wedge (p1, p2, p3, p4, closed, o))
        return s;
    }
    else if (p5 == p2) {
      if (point_in_wedge (p2, p3, p1, p4, closed, o))
        return s;
    }
    else if (p5 == p3) {
      if (point_in_wedge (p3, p1, p2, p4, closed, o))
        return s;
    }
    else if (segment_intersects1 (p1, p2, p4, p5, closed, o) ||
             segment_intersects1 (p2, p3, p4, p5, closed, o) ||
             segment_intersects1 (p3, p1, p4, p5, closed, o))
      return s;
    s = NEXT (s);
  } while (s != start);
  return NULL;
}

static gboolean new_ear (GtsSegment * s, 
                         Ear * e, 
                         GtsSegment * start,
                         guint sloppy,
                         GtsPoint * o)
{
  gdouble or;

  e->s1 = s;
  e->s2 = NEXT (s);

  g_return_val_if_fail (e->s2, FALSE);
  g_return_val_if_fail (e->s2 != e->s1, FALSE);

  ORIENTED_VERTICES (e->s1, e->s2, e->v1, e->v2);
  e->v3 = e->s2->v1 != e->v2 ? e->s2->v1 : e->s2->v2;
  if (e->v3 == e->v1)
    return FALSE;
  e->s3 = NEXT (e->s2);
  if (gts_segment_connect (e->s3, e->v1, e->v3)) {
    if (NEXT (e->s3) != e->s1)
      return FALSE;
  }
  else if (gts_vertices_are_connected (e->v1, e->v3))
    return FALSE;
  else
    e->s3 = NULL;
  or = ORIENTATION (GTS_POINT (e->v1), GTS_POINT (e->v2), GTS_POINT (e->v3),o);
  switch (sloppy) {
  case 0: 
    if (or <= 0. ||
        triangle_intersects_segments (GTS_POINT (e->v1), GTS_POINT (e->v2),
                                      GTS_POINT (e->v3), TRUE, start, o))
      return FALSE;
    break;
  case 1:
    if (or < 0. || 
        (or > 0. && 
         triangle_intersects_segments (GTS_POINT (e->v1), GTS_POINT (e->v2),
                                       GTS_POINT (e->v3), FALSE, start, o)))
      return FALSE;
    break;
  case 2:
    if ((or > 0. && 
         triangle_intersects_segments (GTS_POINT (e->v1), GTS_POINT (e->v2),
                                       GTS_POINT (e->v3), FALSE, start, o)) ||
        (or < 0. && 
         triangle_intersects_segments (GTS_POINT (e->v2), GTS_POINT (e->v1),
                                       GTS_POINT (e->v3), FALSE, start, o)))
      return FALSE;
    break;
  case 3:
    if (or < 0.)
      return FALSE;
    break;
  }
#ifdef DEBUG
  if (or <= 0.)
    fprintf (stderr, "or: %g\n", or);
#endif /* DEBUG */
  g_assert (or > -1e-6);
  return TRUE;
}

static void triangulate_loop (GtsSegment * start,
                              GtsSurface * surface,
                              GtsPoint * o)
{
  GtsSegment * prev = start, * s;
  guint sloppy = 0;
#ifdef DEBUG
  guint nt = 0;
#endif /* DEBUG */

  s = NEXT (start);
  while (NEXT (s) != s) {
    GtsSegment * next = NEXT (s);
    Ear e;

#ifdef DEBUG
    fprintf (stderr, "prev: %p s: %p next: %p\n", prev, s, next);
#endif /* DEBUG */
  
    if (!new_ear (s, &e, start, sloppy, o)) {
      if (s == start) {
        sloppy++;
#ifdef DEBUG
        fprintf (stderr, "sloppy: %u\n", sloppy);
#endif /* DEBUG */
      }
      prev = s;
      s = next;
    }
    else {
      GtsFace * f;

      if (!GTS_IS_EDGE (e.s3))
        e.s3 = GTS_SEGMENT (gts_edge_new (surface->edge_class, e.v1, e.v3));
      f = gts_face_new (surface->face_class, 
                        GTS_EDGE (e.s1), GTS_EDGE (e.s2), GTS_EDGE (e.s3));
      gts_surface_add_face (surface, f);
      UNSET (e.s1, RELEVANT);
      UNSET (e.s1, INTERIOR);
      UNSET (e.s2, RELEVANT);
      UNSET (e.s2, INTERIOR);
      NEXT (prev) = e.s3;
      NEXT (e.s3) = NEXT (e.s2);
      NEXT (e.s1) = NEXT (e.s2) = NULL;
      start = prev;
      s = NEXT (prev);
      sloppy = 0;
#ifdef DEBUG
      {
        gchar name[80];
        FILE * fp;
        
        fprintf (stderr, " t.%u: (%p:%s,%p:%s,%p:%s)\n",
                 nt, 
                 e.v1, NAME (e.v1),
                 e.v2, NAME (e.v2),
                 e.v3, NAME (e.v3));
        sprintf (name, "/tmp/t.%u", nt++);
        fp = fopen (name, "wt");
        //      gts_surface_write (surface, fp);
        gts_write_triangle (GTS_TRIANGLE (f), NULL, fp);
        //        write_graph1 (start, interior, surface, fp);
        fclose (fp);
        print_loop (start, stderr);
      }
#endif /* DEBUG */
    }
  }
  UNSET (s, RELEVANT);
  UNSET (s, INTERIOR);
  NEXT (s) = NULL;
}

#ifdef CHECK_ORIENTED
static void check_object (GtsObject * o)
{
  g_assert (o->reserved == NULL);
  g_assert (o->flags == 0);  
}

static void check_boundary (GtsEdge * e, GtsSurface * s)
{
  check_object (GTS_OBJECT (e));
  check_object (GTS_OBJECT (GTS_SEGMENT (e)->v1));
  check_object (GTS_OBJECT (GTS_SEGMENT (e)->v2));
  g_assert (gts_edge_face_number (e, s) == 1);
}

static void check_interior (GtsEdge * e, GtsSurface * s)
{
  guint n;
  check_object (GTS_OBJECT (e));
  check_object (GTS_OBJECT (GTS_SEGMENT (e)->v1));
  check_object (GTS_OBJECT (GTS_SEGMENT (e)->v2));

  n = gts_edge_face_number (e, s);
#ifdef DEBUG
  if (n != 2)
    gts_surface_print_stats (s, stderr);
#endif /* DEBUG */
  g_assert (n == 2);
}

static void check_boundary_interior_triangulation (GSList * boundary,
                                                   GSList * interior,
                                                   GtsSurface * surface)
{
  g_slist_foreach (boundary, (GFunc) check_boundary, surface);
  g_slist_foreach (interior, (GFunc) check_interior, surface);
}
#endif /*ifdef CHECK_ORIENTED */

static void merge_duplicate (GtsEdge * e)
{
  GtsEdge * dup = gts_edge_is_duplicate (e);

  g_assert (dup);
  gts_edge_replace (dup, e);
  gts_object_destroy (GTS_OBJECT (dup));
}

static void triangulate_boundary_interior (GSList * boundary, 
                                           GSList * interior,
                                           GtsSurface * s,
                                           GtsPoint * o)
{
  GSList * iloops, * bloops, * i;

  i = boundary;
  while (i) {
    SET (i->data, RELEVANT);
    i = i->next;
  }
  i = interior;
  while (i) {
    SET (i->data, RELEVANT);
    SET (i->data, INTERIOR);
    i = i->next;
  }

  iloops = interior_loops (interior);
  bloops = boundary_loops (boundary);

  i = iloops;
  while (i) {
#ifdef DEBUG
    fprintf (stderr, "--- interior loop ---\n");
    print_loop (i->data, stderr);
#endif /* DEBUG */
    connect_interior_loop (i->data, &interior, &bloops, s, o);
    i = i->next;
  }
  
#ifdef DEBUG
 {
   FILE * fp = fopen ("/tmp/bloops", "w");
   write_loops (bloops, fp);
   fclose (fp);
 }
#endif /* DEBUG */

  i = bloops;
  while (i) {
#ifdef DEBUG
    fprintf (stderr, "--- boundary loop ---\n");
    print_loop (i->data, stderr);
#endif /* DEBUG */
    triangulate_loop (i->data, s, o);
    i = i->next;
  }
  
  g_slist_foreach (interior, (GFunc) merge_duplicate, NULL);
  g_slist_free (iloops);
  g_slist_free (bloops);

#ifdef CHECK_ORIENTED
  check_boundary_interior_triangulation (boundary, interior, s);
#endif /* CHECK_ORIENTED */
}

static void create_edges (GtsSegment * s, GtsSurface * surface)
{
  if (GTS_OBJECT (s)->reserved) {
    GList * i = GTS_OBJECT (s)->reserved;
    GtsVertex * v1 = i->data;

    GTS_OBJECT (s)->reserved = g_list_prepend (i, 
                      gts_edge_new (surface->edge_class, s->v1, v1));
    while (i) {
      GList * next = i->next;
      GtsVertex * v2 = next ? next->data : s->v2;

      GTS_OBJECT (i->data)->reserved = NULL;
      i->data = gts_edge_new (surface->edge_class, v1, v2);
      v1 = v2;
      i = next;
    }
  }
}

static void add_boundary (GtsSegment * s, GtsSegment * next, 
                          GSList ** boundary)
{
  if (GTS_OBJECT (s)->reserved == NULL)
    *boundary = g_slist_prepend (*boundary, s);
  else {
    if (s->v2 == next->v2 || s->v2 == next->v1) {
      GList * i = g_list_last (GTS_OBJECT (s)->reserved);

      while (i) {
        *boundary = g_slist_prepend (*boundary, i->data);
        i = i->prev;
      }
    }
    else {
      GList * i = GTS_OBJECT (s)->reserved;

      while (i) {
        *boundary = g_slist_prepend (*boundary, i->data);
        i = i->next;
      }
    }
  }
}

static void triangulate_face (GtsTriangle * t, GtsSurface * surface)
{
  GSList * interior = GTS_OBJECT (t)->reserved;
  GSList * boundary = NULL;
  GtsSurface * s = gts_surface_new (gts_surface_class (),
                                    surface->face_class,
                                    surface->edge_class,
                                    surface->vertex_class);
  gdouble x, y, z;
  GtsPoint * p = GTS_POINT (GTS_SEGMENT (t->e1)->v1);
  GtsPoint * o;

  GTS_OBJECT (t)->reserved = NULL;  
  gts_triangle_normal (t, &x, &y, &z);
  g_assert (x != 0. || y != 0. || z != 0.);
  o = gts_point_new (gts_point_class (), p->x + x, p->y + y, p->z + z);
  add_boundary (GTS_SEGMENT (t->e3), GTS_SEGMENT (t->e1), &boundary);
  add_boundary (GTS_SEGMENT (t->e2), GTS_SEGMENT (t->e3), &boundary);
  add_boundary (GTS_SEGMENT (t->e1), GTS_SEGMENT (t->e2), &boundary);
#ifdef DEBUG
  {
    static guint nt = 0;
    char name[80];
    FILE * fp;

    fprintf (stderr, "%u: triangulating %p\n", nt, t);
if (nt == 28)
  fprintf (stderr, "tintin!!!!\n");
    sprintf (name, "/tmp/oc.%u", nt++);
    fp = fopen (name, "wt");
    //    write_graph (boundary, interior, s, fp);
    write_segments (boundary, interior, fp);
    fclose (fp);
  }
#endif /* DEBUG */
  triangulate_boundary_interior (boundary, interior, s, o);
  g_slist_free (interior);
  g_slist_free (boundary);
  if (GTS_OBJECT (t)->klass->attributes)
    gts_surface_foreach_face (s, (GtsFunc) gts_object_attributes, t);
  gts_surface_merge (surface, s);
  gts_object_destroy (GTS_OBJECT (s));
  gts_object_destroy (GTS_OBJECT (o));
}

static void free_edge_list (GtsObject * o)
{
  g_list_free (o->reserved);
  o->reserved = NULL;
}

GtsSurfaceInter * gts_surface_inter_new (GtsSurfaceInterClass * klass,
                                         GtsSurface * s1,
                                         GtsSurface * s2,
                                         GNode * faces_tree1,
                                         GNode * faces_tree2,
                                         gboolean is_open1,
                                         gboolean is_open2)
{
  GtsSurfaceInter * si;
  GtsSurface * s;

  g_return_val_if_fail (klass != NULL, NULL);
  g_return_val_if_fail (s1 != NULL, NULL);
  g_return_val_if_fail (s2 != NULL, NULL);
  g_return_val_if_fail (faces_tree1 != NULL, NULL);
  g_return_val_if_fail (faces_tree2 != NULL, NULL);

  si = surface_inter_new (klass, s1, s2, faces_tree1, faces_tree2);

  gts_surface_foreach_edge (si->s1, (GtsFunc) create_edges, si->s1);
  gts_surface_foreach_edge (si->s2, (GtsFunc) create_edges, si->s2);

#ifdef DEBUG
  fprintf (stderr, "====== triangulating s1 ======\n");
#endif /* DEBUG */
  s = gts_surface_new (gts_surface_class (),
                       s1->face_class,
                       s1->edge_class,
                       s1->vertex_class);
  gts_surface_foreach_face (si->s1, (GtsFunc) triangulate_face, s);
  gts_surface_foreach_edge (si->s1, (GtsFunc) free_edge_list, NULL);
  gts_object_destroy (GTS_OBJECT (si->s1));
  si->s1 = s;
  GTS_OBJECT (si->s1)->reserved = s1;
  
#ifdef DEBUG
  fprintf (stderr, "====== triangulating s2 ======\n");
#endif /* DEBUG */
  s = gts_surface_new (gts_surface_class (),
                       s2->face_class,
                       s2->edge_class,
                       s2->vertex_class);
  gts_surface_foreach_face (si->s2, (GtsFunc) triangulate_face, s);
  gts_surface_foreach_edge (si->s2, (GtsFunc) free_edge_list, NULL);
  gts_object_destroy (GTS_OBJECT (si->s2));
  si->s2 = s;
  GTS_OBJECT (si->s2)->reserved = s2;

  return si;
}

static void check_surface_edge (GtsEdge * e, gpointer * data)
{
  gboolean * ok = data[0];
  GtsSurface * s = data[1];
  GtsSurface * bs = GTS_OBJECT (s)->reserved;
  guint nf = gts_edge_face_number (e, s);

  if (nf < 1 || nf > 2) {
    *ok = FALSE;
    g_return_if_fail (nf >= 1 && nf <= 2);
  }
  if (nf == 1 && gts_edge_face_number (e, bs) == 0) {
    *ok = FALSE;
    g_return_if_fail (gts_edge_face_number (e, bs) > 0);
  }
}

static void mark_edge (GtsObject * o, gpointer data)
{
  o->reserved = data;
}

static gint triangle_orientation (GtsTriangle * t, GtsEdge * e)
{
  GtsSegment * s = GTS_SEGMENT (t->e1 == e ? t->e2 
                                : 
                                t->e2 == e ? t->e3 
                                : 
                                t->e1);
  GtsVertex * v2 = GTS_SEGMENT (e)->v2;

  if (s->v1 == v2 || s->v2 == v2)
    return 1;
  return -1;
}

static gboolean check_orientation (GtsEdge * e, GtsSurface * s)
{
  GtsTriangle * t1 = NULL, * t2 = NULL;
  GSList * i = e->triangles;
  gint o1 = 0, o2 = 0;

  while (i) {
    if (GTS_IS_FACE (i->data) && 
        gts_face_has_parent_surface (i->data, s)) {
      if (t1 == NULL) {
        t1 = i->data;
        o1 = triangle_orientation (t1, e);
      }
      else if (t2 == NULL) {
        t2 = i->data;
        o2 = triangle_orientation (t2, e);
        g_return_val_if_fail (o1*o2 < 0, FALSE);
      }
      else
        g_assert_not_reached ();
    }
    i = i->next;
  }
  g_return_val_if_fail (t1 && t2, FALSE);
  return TRUE;
}

static void check_edge (GtsSegment * s, gpointer * data)
{
  gboolean * ok = data[0];
  GtsSurfaceInter * si = data[1];
  gboolean * closed = data[2];
  GSList * j;
  guint nn = 0;
  
  j = s->v1->segments;
  while (j && *ok) {
    GtsSegment * s1 = j->data;
    
    if (s1 != s && GTS_OBJECT (s1)->reserved == si) {
      if (s1->v2 != s->v1)
        *ok = FALSE;
      nn++;
    }
    j = j->next;
  }
  j = s->v2->segments;
  while (j && *ok) {
    GtsSegment * s1 = j->data;
    
    if (s1 != s && GTS_OBJECT (s1)->reserved == si) {
      if (s1->v1 != s->v2)
        *ok = FALSE;
      nn++;
    }
    j = j->next;
  }
  if (nn != 2)
    *closed = FALSE;

  if (!check_orientation (GTS_EDGE (s), si->s1))
    *ok = FALSE;
  if (!check_orientation (GTS_EDGE (s), si->s2))
    *ok = FALSE;
}

gboolean gts_surface_inter_check (GtsSurfaceInter * si,
                                  gboolean * closed)
{
  gboolean ok = TRUE;
  gpointer data[3];

  g_return_val_if_fail (si != NULL, FALSE);
  g_return_val_if_fail (closed != NULL, FALSE);

  *closed = si->edges ? TRUE : FALSE;

  /* mark edges as used by si */
  g_slist_foreach (si->edges, (GFunc) mark_edge, si);

  data[0] = &ok;
  data[1] = si;
  data[2] = closed;
  g_slist_foreach (si->edges, (GFunc) check_edge, data);
  g_slist_foreach (si->edges, (GFunc) gts_object_reset_reserved, NULL);

  /* check connectivity of the faces of @si */
  if (*closed) {
    gpointer data[2];

    data[0] = &ok;
    data[1] = si->s1;
    gts_surface_foreach_edge (si->s1, (GtsFunc) check_surface_edge, data);
    data[1] = si->s2;
    gts_surface_foreach_edge (si->s2, (GtsFunc) check_surface_edge, data);
  }

  return ok;
}

/* Given @e and @f returns a #GtsFace compatible with @f and belonging to
   @s1 or @s2 */
static GtsFace * next_compatible_face (GtsEdge * e, 
                                       GtsFace * f, 
                                       GtsSurface * s1,
                                       GtsSurface * s2)
{
  GSList * i = e->triangles;
  GtsFace * f2 = NULL, * f3 = NULL;

  while (i) {
    GtsFace * f1 = i->data;

    if (f1 != f && GTS_IS_FACE (f1)) {
      if (gts_face_has_parent_surface (f1, s1))
        return f1;
      if (gts_face_has_parent_surface (f1, s2)) {
        if (f2 == NULL) f2 = f1;
        else if (f3 == NULL) f3 = f1;
        else g_assert_not_reached (); /* s2 is a non-manifold surface */
      }
    }
    i = i->next;
  }
  if (f3 == NULL) {
    if (gts_edge_is_boundary (e, s2))
      return NULL;
    return f2; 
  }
  g_assert (gts_face_has_parent_surface (f, s1));
  if (gts_triangles_are_compatible (GTS_TRIANGLE (f), GTS_TRIANGLE (f2), e))
    return f2;
  return f3;
}

static void walk_faces (GtsEdge * e, GtsFace * f, 
                        GtsSurface * s1,
                        GtsSurface * s2,
                        GtsSurface * s)
{
  GtsFifo * faces = gts_fifo_new ();
  GtsFifo * edges = gts_fifo_new ();

  gts_fifo_push (faces, f);
  gts_fifo_push (edges, e);
  while ((f = gts_fifo_pop (faces)) && (e = gts_fifo_pop (edges))) {
    if (!GTS_OBJECT (f)->reserved) {
      GtsTriangle * t = GTS_TRIANGLE (f);
      GtsFace * f1;

      gts_surface_add_face (s, f);
      GTS_OBJECT (f)->reserved = s;
      if (t->e1 != e && !GTS_OBJECT (t->e1)->reserved &&
          (f1 = next_compatible_face (t->e1, f, s1, s2))) {
        gts_fifo_push (faces, f1);
        gts_fifo_push (edges, t->e1);
      } 
      if (t->e2 != e && !GTS_OBJECT (t->e2)->reserved &&
          (f1 = next_compatible_face (t->e2, f, s1, s2))) {
        gts_fifo_push (faces, f1);
        gts_fifo_push (edges, t->e2);
      } 
      if (t->e3 != e && !GTS_OBJECT (t->e3)->reserved &&
          (f1 = next_compatible_face (t->e3, f, s1, s2))) {
        gts_fifo_push (faces, f1);
        gts_fifo_push (edges, t->e3);
      } 
    }
  }
  gts_fifo_destroy (faces);
  gts_fifo_destroy (edges);
}

void gts_surface_inter_boolean (GtsSurfaceInter * si,
                                GtsSurface * surface,
                                GtsBooleanOperation op)
{
  GtsSurface * s = NULL;
  gint orient = 1;
  GSList * i;

  g_return_if_fail (si != NULL);
  g_return_if_fail (surface != NULL);

  switch (op) {
  case GTS_1_OUT_2: s = si->s1; orient = 1; break;
  case GTS_1_IN_2: s = si->s1; orient = -1; break;
  case GTS_2_OUT_1: s = si->s2; orient = -1; break;
  case GTS_2_IN_1: s = si->s2; orient = 1; break;
  default: g_assert_not_reached ();
  }

  /* mark edges as belonging to intersection */
  g_slist_foreach (si->edges, (GFunc) mark_edge, si);

  i = si->edges;
  while (i) {
    GtsEdge * e = i->data;
    GSList * j = e->triangles;
    
    while (j) {
      if (gts_face_has_parent_surface (j->data, s) &&
          orient*triangle_orientation (j->data, e) > 0) {
#ifdef DEBUG_BOOLEAN
        GtsFace * boundary = gts_edge_is_boundary (e, surface);

        g_assert (!boundary || boundary == j->data);
#endif /* DEBUG_BOOLEAN */
        walk_faces (e, j->data, s, GTS_OBJECT (s)->reserved, surface);
        break;
      }
      j = j->next;
    }
    i = i->next;
  }
  g_slist_foreach (si->edges, (GFunc) gts_object_reset_reserved, NULL);
  gts_surface_foreach_face (surface, 
                            (GtsFunc) gts_object_reset_reserved, NULL);
}

static void self_intersecting (GtsBBox * bb1, GtsBBox * bb2, 
                               gpointer * d)
{
  GtsTriangle * t1 = bb1->bounded;
  GtsTriangle * t2 = bb2->bounded;

  if (t1 != t2) {
    GtsSegment * s1 = GTS_SEGMENT (t1->e1);
    GtsSegment * s2 = GTS_SEGMENT (t1->e2);
    GtsSegment * s3 = GTS_SEGMENT (t1->e3);
    GtsSegment * s4 = GTS_SEGMENT (t2->e1);
    GtsSegment * s5 = GTS_SEGMENT (t2->e2);
    GtsSegment * s6 = GTS_SEGMENT (t2->e3);
    GtsPoint * pi;

    if ((!gts_segments_touch (s4, s1) && 
         !gts_segments_touch (s4, s2) &&
         !gts_segments_touch (s4, s3) &&
         (pi = segment_triangle_intersection (s4, t1, gts_point_class ()))
         != NULL) ||
        (!gts_segments_touch (s5, s1) && 
         !gts_segments_touch (s5, s2) &&
         !gts_segments_touch (s5, s3) &&
         (pi = segment_triangle_intersection (s5, t1, gts_point_class ())) 
         != NULL) ||
        (!gts_segments_touch (s6, s1) && 
         !gts_segments_touch (s6, s2) &&
         !gts_segments_touch (s6, s3) &&
         (pi = segment_triangle_intersection (s6, t1, gts_point_class ())) 
         != NULL)) {
      GtsBBTreeTraverseFunc func = d[0];
      gpointer data = d[1];
      gboolean * self_inter = d[2];

      gts_object_destroy (GTS_OBJECT (pi));
      *self_inter = TRUE;
      (* func) (bb1, bb2, data);
    }
  }
}

gboolean gts_surface_foreach_intersecting_face (GtsSurface * s,
                                                GtsBBTreeTraverseFunc func,
                                                gpointer data)
{
  GNode * tree;
  gpointer d[3];
  gboolean self_inter = FALSE;

  g_return_val_if_fail (s != NULL, FALSE);
  g_return_val_if_fail (func != NULL, FALSE);

  tree = gts_bb_tree_surface (s);
  d[0] = func;
  d[1] = data;
  d[2] = &self_inter;
  gts_bb_tree_traverse_overlapping (tree, tree, 
                                    (GtsBBTreeTraverseFunc) self_intersecting,
                                    d);
  gts_bb_tree_destroy (tree, TRUE);

  return self_inter;
}

static void add_intersecting (GtsBBox * bb1, GtsBBox * bb2, 
                              GtsSurface * intersected)
{
  gts_surface_add_face (intersected, bb1->bounded);
  gts_surface_add_face (intersected, bb2->bounded);
}

GtsSurface * gts_surface_is_self_intersecting (GtsSurface * s)
{
  GtsSurface * intersected;

  g_return_val_if_fail (s != NULL, NULL);

  intersected = gts_surface_new (GTS_SURFACE_CLASS (GTS_OBJECT (s)->klass),
                                 s->face_class,
                                 s->edge_class,
                                 s->vertex_class);
  if (!gts_surface_foreach_intersecting_face (s,
                      (GtsBBTreeTraverseFunc) add_intersecting, intersected)) {
    gts_object_destroy (GTS_OBJECT (intersected));
    intersected = NULL;
  }
  return intersected;
}