cdt.c

Go to the documentation of this file.

/* GTS - Library for the manipulation of triangulated surfaces
 * Copyright (C) 1999 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.
 */

#ifdef HAVE_CONFIG_H
# include <config.h>
#endif


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

#ifdef USE_SURFACE_BTREE

static gint find_closest (GtsTriangle * t, gpointer value, gpointer * data)
{
  guint * ns = data[2];
  guint * n = data[3];

  if (*n >= *ns)
    return TRUE;
  else {
    gdouble * dmin = data[0];
    gpointer * closest = data[1];
    GtsPoint * p = data[4];

    if (gts_triangle_orientation (t) > 0.) {
      GtsPoint * p1 = GTS_POINT (GTS_SEGMENT (t->e1)->v1);
      gdouble d = (p->x - p1->x)*(p->x - p1->x) + (p->y - p1->y)*(p->y - p1->y);
      
      if (d < *dmin) {
        *dmin = d;
        *closest = t;
      }
      (*n)++;
    }
  }
  return FALSE;
}

/* select the face closest to @p among n^1/3 randomly picked faces
 *  of @surface */
static GtsFace * closest_face (GtsSurface * s, GtsPoint * p)
{
  guint n = 0, nt, ns;
  gdouble dmin = G_MAXDOUBLE;
  GtsFace * closest = NULL;
  gpointer data[5];

  nt = gts_surface_face_number (s);
  if (!nt)
    return NULL;
  ns = exp (log ((gdouble) nt)/3.);

  data[0] = &dmin;
  data[1] = &closest;
  data[2] = &ns;
  data[3] = &n;
  data[4] = p;
  g_tree_traverse (s->faces, (GTraverseFunc) find_closest, G_IN_ORDER, data);

  return closest;
}

#else /* not USE_SURFACE_BTREE */

typedef struct _SFindClosest SFindClosest; 

struct _SFindClosest {
  gdouble dmin; 
  GtsFace *closest;
  GtsPoint * p;
  gint stop;
};

#  if GLIB_CHECK_VERSION(2,4,0)
/* finally, with g_hash_table_find we are able to stop iteration over the hash
   table in the middle */

static gboolean find_closest (gpointer key, gpointer value, gpointer user_data)
{
  SFindClosest * data = (SFindClosest *) user_data;
  GtsFace * f = GTS_FACE (value);
  
  if (gts_triangle_orientation (GTS_TRIANGLE (f)) > 0.) {
    GtsPoint * p1 = GTS_POINT (GTS_SEGMENT (GTS_TRIANGLE (f)->e1)->v1);
    gdouble d = ((data->p->x - p1->x)*(data->p->x - p1->x) + 
                 (data->p->y - p1->y)*(data->p->y - p1->y));

    if (d < data->dmin) {
      data->dmin = d;
      data->closest = f;
    }
  }
  data->stop--;
  return !(data->stop > 0);
}

static GtsFace * closest_face (GtsSurface * s, GtsPoint * p)
{
  SFindClosest fc;

  fc.dmin = G_MAXDOUBLE;
  fc.closest = NULL;
  fc.p = p;
  fc.stop = (gint) exp (log ((gdouble) g_hash_table_size (s->faces))/3.);
  g_hash_table_find (s->faces, find_closest, &fc);
  
  return fc.closest;
}

#  else /* VERSION < 2.4.0 */

static void
find_closest (gpointer key, gpointer value, gpointer user_data)
{
  SFindClosest * data = (SFindClosest *) user_data;
  GtsFace * f = GTS_FACE (value);

  if (gts_triangle_orientation (GTS_TRIANGLE (f)) > 0.) {
    GtsPoint * p1 = GTS_POINT (GTS_SEGMENT (GTS_TRIANGLE (f)->e1)->v1);
    gdouble d = ((data->p->x - p1->x)*(data->p->x - p1->x) +
                 (data->p->y - p1->y)*(data->p->y - p1->y));

    if (d < data->dmin) {
      data->dmin = d;
      data->closest = f;
    }
  }
  data->stop--;
}

/* select the face closest to @p among n^1/3 randomly picked faces
 * of @surface */
static GtsFace * closest_face (GtsSurface * s, GtsPoint * p)
{
  SFindClosest fc;

  if (!g_hash_table_size (s->faces))
    return NULL;

  fc.dmin = G_MAXDOUBLE;
  fc.closest = NULL;
  fc.p = p;
  fc.stop = (gint) exp (log ((gdouble) g_hash_table_size (s->faces))/3.);
  g_hash_table_foreach (s->faces, find_closest, &fc);
  return fc.closest;
}
#  endif /* VERSION < 2.4.0 */
#endif /* not USE_SURFACE_BTREE */

/* returns the face belonging to @surface and neighbor of @f via @e */
static GtsFace * neighbor (GtsFace * f,
                           GtsEdge * e,
                           GtsSurface * surface)
{
  GSList * i = e->triangles;
  GtsTriangle * t = GTS_TRIANGLE (f);

  while (i) {
    GtsTriangle * t1 = i->data;
    if (t1 != t &&
        GTS_IS_FACE (t1) &&
        gts_face_has_parent_surface (GTS_FACE (t1), surface))
      return GTS_FACE (t1);
    i = i->next;
  }
  return NULL;
}

/* given a triangle @t and a segment s (@o -> @p). 
   @o must be in @t. Returns the
   edge of @t which is intersected by s or %NULL if @p is also
   contained in @t (on_summit is set to %FALSE) or if s intersects @t 
   exactly on one of its summit (on_summit is set to %TRUE). */
static GtsEdge * triangle_next_edge (GtsTriangle * t,
                                     GtsPoint * o, GtsPoint * p,
                                     gboolean * on_summit)
{
  GtsVertex * v1, * v2, * v3;
  GtsEdge * e1, * e2, * e3;
  gdouble orient = 0.0;
  
  gts_triangle_vertices_edges (t, NULL,
                               &v1, &v2, &v3, 
                               &e1, &e2, &e3);

  *on_summit = FALSE;
  orient = gts_point_orientation (o, GTS_POINT (v1), p);
  if (orient > 0.0) {
    orient = gts_point_orientation (o, GTS_POINT (v2), p);
    if (orient > 0.0) {
      if (gts_point_orientation (GTS_POINT (v2), GTS_POINT (v3), p) >= 0.0)
        return NULL;
      return e2;
    }
    if (orient < 0.0) {
      if (gts_point_orientation (GTS_POINT (v1), GTS_POINT (v2), p) >= 0.0)
        return NULL;
      return e1;
    }
    if (gts_point_orientation (GTS_POINT (v1), GTS_POINT (v2), p) < 0.0)
      *on_summit = TRUE;
    return NULL;
  }

  if (orient < 0.0) {
    orient = gts_point_orientation (o, GTS_POINT (v3), p);
    if (orient > 0.0) {
      if (gts_point_orientation (GTS_POINT (v3), GTS_POINT (v1), p) >= 0.0)
        return NULL;
      return e3;
    }
    if (orient < 0.0) {
      if (gts_point_orientation (GTS_POINT (v2), GTS_POINT (v3), p) >= 0.0)
        return NULL;
      return e2;
    }
    if (gts_point_orientation (GTS_POINT (v3), GTS_POINT (v1), p) < 0.0)
      *on_summit = TRUE;
    return NULL;
  }
  
  if (gts_point_orientation (GTS_POINT (v2), GTS_POINT (v3), p) < 0.0)
    return e2;
  if (gts_point_orientation (GTS_POINT (v1), GTS_POINT (v2), p) < 0.0)
    *on_summit = TRUE;
  return NULL;
}

static void triangle_barycenter (GtsTriangle * t, GtsPoint * b)
{
  GtsPoint * p = GTS_POINT (gts_triangle_vertex (t));
  b->x = (p->x + 
          GTS_POINT (GTS_SEGMENT(t->e1)->v1)->x +
          GTS_POINT (GTS_SEGMENT(t->e1)->v2)->x)/3.;
  b->y = (p->y + 
          GTS_POINT (GTS_SEGMENT(t->e1)->v1)->y +
          GTS_POINT (GTS_SEGMENT(t->e1)->v2)->y)/3.;
}

static GtsFace * point_locate (GtsPoint * o,
                               GtsPoint * p,
                               GtsFace * f,
                               GtsSurface * surface)
{
  GtsEdge * prev;
  gboolean on_summit;
  GtsVertex * v1, * v2, * v3;
  GtsEdge * e2, * e3;    
  
  prev = triangle_next_edge (GTS_TRIANGLE (f), o, p, &on_summit);

  if (!prev) {
    GtsFace * f1;

    if (!on_summit)
      return f; /* p is inside f */

    /* s intersects f exactly on a summit: restarts from a neighbor of f */
    if ((f1 = neighbor (f, GTS_TRIANGLE (f)->e1, surface)) ||
        (f1 = neighbor (f, GTS_TRIANGLE (f)->e2, surface)) ||
        (f1 = neighbor (f, GTS_TRIANGLE (f)->e3, surface))) {
      triangle_barycenter (GTS_TRIANGLE (f1), o);
      return point_locate (o, p, f1, surface);
    }
    return NULL;
  }
  
  f = neighbor (f, prev, surface);
  if (f)
    gts_triangle_vertices_edges (GTS_TRIANGLE (f), prev, 
                                 &v1, &v2, &v3, &prev, &e2, &e3);
  while (f) {
    gdouble orient = gts_point_orientation (o, GTS_POINT (v3), p);

    if (orient < 0.0) {
      if (gts_point_orientation (GTS_POINT (v2), GTS_POINT (v3), p) >= 0.0)
        return f; /* p is inside f */
      f = neighbor (f, e2, surface);
      prev = e2;
      v1 = v3;      
    }
    else if (orient > 0.0) {
      if (gts_point_orientation (GTS_POINT (v3), GTS_POINT (v1), p) >= 0.0)
        return f; /* p is inside f */
      f = neighbor (f, e3, surface);
      prev = e3;
      v2 = v3;
    }
    else {
      GtsFace * f1;

      if (gts_point_orientation (GTS_POINT (v2), GTS_POINT (v3), p) >= 0.0)
        return f; /* p is inside f */

      /* s intersects f exactly on v3: restarts from a neighbor of f */
      if ((f1 = neighbor (f, e2, surface)) ||
          (f1 = neighbor (f, e3, surface))) {
        triangle_barycenter (GTS_TRIANGLE (f1), o);
        return point_locate (o, p, f1, surface);
      }
      return NULL;
    }
    /* update e2, e3, v3 for the new triangle */
    if (f) {
      if (prev == GTS_TRIANGLE (f)->e1) {
        e2 = GTS_TRIANGLE (f)->e2; e3 = GTS_TRIANGLE (f)->e3;
      }
      else if (prev == GTS_TRIANGLE (f)->e2) {
        e2 = GTS_TRIANGLE (f)->e3; e3 = GTS_TRIANGLE (f)->e1;
      }
      else {
        e2 = GTS_TRIANGLE (f)->e1; e3 = GTS_TRIANGLE (f)->e2;
      }
      if (GTS_SEGMENT (e2)->v1 == v1 || GTS_SEGMENT (e2)->v1 == v2)
        v3 = GTS_SEGMENT (e2)->v2;
      else
        v3 = GTS_SEGMENT (e2)->v1;
    }
  }
  return NULL;
}

GtsFace * gts_point_locate (GtsPoint * p, 
                            GtsSurface * surface,
                            GtsFace * guess)
{
  GtsFace * fr;
  GtsPoint * o;

  g_return_val_if_fail (p != NULL, NULL);
  g_return_val_if_fail (surface != NULL, NULL);
  g_return_val_if_fail (guess == NULL || 
                        gts_face_has_parent_surface (guess, surface), NULL);

  if (guess == NULL)
    guess = closest_face (surface, p);
  else
    g_return_val_if_fail (gts_triangle_orientation (GTS_TRIANGLE (guess)) > 0., NULL);

  if (guess == NULL)
    return NULL;

  o = GTS_POINT (gts_object_new (GTS_OBJECT_CLASS (gts_point_class ())));
  triangle_barycenter (GTS_TRIANGLE (guess), o);
  fr = point_locate (o, p, guess, surface);
  gts_object_destroy (GTS_OBJECT (o));

  return fr;
}


GtsConstraintClass * gts_constraint_class (void)
{
  static GtsConstraintClass * klass = NULL;

  if (klass == NULL) {
    GtsObjectClassInfo constraint_info = {
      "GtsConstraint",
      sizeof (GtsConstraint),
      sizeof (GtsConstraintClass),
      (GtsObjectClassInitFunc) NULL,
      (GtsObjectInitFunc) NULL,
      (GtsArgSetFunc) NULL,
      (GtsArgGetFunc) NULL
    };
    klass = gts_object_class_new (GTS_OBJECT_CLASS (gts_edge_class ()), 
                                  &constraint_info);
  }

  return klass;
}

static void split_list (GtsListFace * f, GtsListFace * f1, GtsListFace * f2, 
                        GtsPoint * p1, GtsPoint * p2,
                        GSList ** last1, GSList ** last2)
{
  GSList * i = f->points, * l1 = *last1, * l2 = *last2;

  while (i) {
    GtsPoint * p = i->data;
    
    if (gts_point_orientation (p1, p2, p) >= 0.) {
      if (l1) l1->next = i; else f1->points = i;
      l1 = i;
    }
    else {
      if (l2) l2->next = i; else f2->points = i;
      l2 = i;
    }
    i = i->next;
  }
  f->points = NULL;
  *last1 = l1;
  *last2 = l2;
}

/* cf. figure misc/swap.fig */
static void swap_if_in_circle (GtsFace * f1,
                               GtsVertex * v1, 
                               GtsVertex * v2, 
                               GtsVertex * v3,
                               GtsEdge * e1, 
                               GtsEdge * e2, 
                               GtsEdge * e3,
                               GtsSurface * surface)
{
  GtsFace * f2;
  GtsEdge * e4, *e5;
  GtsVertex * v4;

  if (GTS_IS_CONSTRAINT (e1)) /* @e1 is a constraint can not swap */
    return;

  f2 = neighbor (f1, e1, surface);
  if (f2 == NULL) /* @e1 is a boundary of @surface */
    return;

  if (GTS_TRIANGLE (f2)->e1 == e1) {
    e4 = GTS_TRIANGLE (f2)->e2; e5 = GTS_TRIANGLE (f2)->e3;
  }
  else if (GTS_TRIANGLE (f2)->e2 == e1) {
    e4 = GTS_TRIANGLE (f2)->e3; e5 = GTS_TRIANGLE (f2)->e1;
  }
  else {
    e4 = GTS_TRIANGLE (f2)->e1; e5 = GTS_TRIANGLE (f2)->e2;
  }
  if (GTS_SEGMENT (e4)->v1 == GTS_SEGMENT (e1)->v1 || 
      GTS_SEGMENT (e4)->v1 == GTS_SEGMENT (e1)->v2)
    v4 = GTS_SEGMENT (e4)->v2;
  else
    v4 = GTS_SEGMENT (e4)->v1;

  if (gts_point_in_circle (GTS_POINT (v4), GTS_POINT (v1), 
                           GTS_POINT (v2), GTS_POINT (v3)) > 0.0) {
    GtsEdge * en;
    GtsSegment * sn = gts_vertices_are_connected (v3, v4);
    GtsFace * f3, * f4;

    if (!GTS_IS_EDGE (sn))
      en = gts_edge_new (surface->edge_class, v3, v4);
    else
      en = GTS_EDGE (sn);

    f3 = gts_face_new (surface->face_class, en, e5, e2);
    gts_object_attributes (GTS_OBJECT (f3), GTS_OBJECT (f1));
    f4 = gts_face_new (surface->face_class, en, e3, e4);
    gts_object_attributes (GTS_OBJECT (f4), GTS_OBJECT (f2));
    
    if (GTS_IS_LIST_FACE (f3)) {
      GSList * last3 = NULL, * last4 = NULL;

      if (GTS_IS_LIST_FACE (f1))
        split_list (GTS_LIST_FACE (f1), GTS_LIST_FACE (f3), GTS_LIST_FACE (f4),
                    GTS_POINT (v3), GTS_POINT (v4), &last3, &last4);
      if (GTS_IS_LIST_FACE (f2))
        split_list (GTS_LIST_FACE (f2), GTS_LIST_FACE (f3), GTS_LIST_FACE (f4),
                    GTS_POINT (v3), GTS_POINT (v4), &last3, &last4);
      if (last3) last3->next = NULL;
      if (last4) last4->next = NULL;
    }

    gts_surface_remove_face (surface, f1);
    gts_surface_remove_face (surface, f2);
    gts_surface_add_face (surface, f3);
    gts_surface_add_face (surface, f4);

    swap_if_in_circle (f3, v4, v2, v3, e5, e2, en, surface);
    swap_if_in_circle (f4, v1, v4, v3, e4, en, e3, surface);
  }
}

GtsVertex * gts_delaunay_add_vertex_to_face (GtsSurface * surface, 
                                             GtsVertex * v,
                                             GtsFace * f)
{
  GtsEdge * e1, * e2, * e3;
  GtsSegment * s4, * s5, * s6;
  GtsEdge * e4, * e5, * e6;
  GtsVertex * v1, * v2, * v3;
  GtsFace * nf[3];

  g_return_val_if_fail (surface != NULL, v);
  g_return_val_if_fail (v != NULL, v);
  g_return_val_if_fail (f != NULL, v);

  gts_triangle_vertices_edges (GTS_TRIANGLE (f), NULL, 
                               &v1, &v2, &v3, &e1, &e2, &e3);
  if (v == v1 || v == v2 || v == v3) /* v already in @surface */
    return NULL;
  if (GTS_POINT (v)->x == GTS_POINT (v1)->x &&
      GTS_POINT (v)->y == GTS_POINT (v1)->y)
    return v1;
  if (GTS_POINT (v)->x == GTS_POINT (v2)->x &&
      GTS_POINT (v)->y == GTS_POINT (v2)->y)
    return v2;
  if (GTS_POINT (v)->x == GTS_POINT (v3)->x &&
      GTS_POINT (v)->y == GTS_POINT (v3)->y)
    return v3;

  s4 = gts_vertices_are_connected (v, v1);
  if (!GTS_IS_EDGE (s4))
    e4 = gts_edge_new (surface->edge_class, v, v1);
  else
    e4 = GTS_EDGE (s4);
  s5 = gts_vertices_are_connected (v, v2);
  if (!GTS_IS_EDGE (s5))
    e5 = gts_edge_new (surface->edge_class, v, v2);
  else
    e5 = GTS_EDGE (s5);
  s6 = gts_vertices_are_connected (v, v3);
  if (!GTS_IS_EDGE (s6))
    e6 = gts_edge_new (surface->edge_class, v, v3);
  else
    e6 = GTS_EDGE (s6);

  /* cf. figure misc/swap.fig */
  nf[0] = gts_face_new (surface->face_class, e4, e1, e5);
  gts_object_attributes (GTS_OBJECT (nf[0]), GTS_OBJECT (f));
  nf[1] = gts_face_new (surface->face_class, e5, e2, e6);
  gts_object_attributes (GTS_OBJECT (nf[1]), GTS_OBJECT (f));
  nf[2] = gts_face_new (surface->face_class, e6, e3, e4);
  gts_object_attributes (GTS_OBJECT (nf[2]), GTS_OBJECT (f));

  if (GTS_IS_LIST_FACE (f) && GTS_IS_LIST_FACE (nf[0])) {
    GSList * i = GTS_LIST_FACE (f)->points, * last[3] = { NULL, NULL, NULL };

    while (i) {
      GtsPoint * p = i->data;
      GSList * next = i->next;
      guint j;
      
      if (p != GTS_POINT (v)) {
        if (gts_point_orientation (GTS_POINT (v), GTS_POINT (v1), p) >= 0.) {
          gdouble o = gts_point_orientation (GTS_POINT (v), GTS_POINT (v2), p);

          if (o != 0.)
            j = o > 0. ? 1 : 0;
          else
            j = gts_point_orientation (GTS_POINT (v), GTS_POINT (v3), p) 
              > 0. ? 0 : 1;
        }
        else if (gts_point_orientation (GTS_POINT (v), GTS_POINT (v3), p) > 0.)
          j = 2;
        else
          j = 1;
        if (last[j])
          last[j]->next = i; 
        else 
          GTS_LIST_FACE (nf[j])->points = i;
        last[j] = i;
      }
      else
        g_slist_free_1 (i);
      i = next;
    }
    GTS_LIST_FACE (f)->points = NULL;
    if (last[0]) last[0]->next = NULL;
    if (last[1]) last[1]->next = NULL;
    if (last[2]) last[2]->next = NULL;
  }

  gts_surface_remove_face (surface, f);
  gts_surface_add_face (surface, nf[0]);
  gts_surface_add_face (surface, nf[1]);
  gts_surface_add_face (surface, nf[2]);

  swap_if_in_circle (nf[0], v1, v2, v, e1, e5, e4, surface);
  swap_if_in_circle (nf[1], v2, v3, v, e2, e6, e5, surface);
  swap_if_in_circle (nf[2], v3, v1, v, e3, e4, e6, surface);

  return NULL;
}

GtsVertex * gts_delaunay_add_vertex (GtsSurface * surface, 
                                     GtsVertex * v,
                                     GtsFace * guess)
{
  GtsFace * f;

  g_return_val_if_fail (surface != NULL, v);
  g_return_val_if_fail (v != NULL, v);

  if (!(f = gts_point_locate (GTS_POINT (v), surface, guess)))
    return v;
  return gts_delaunay_add_vertex_to_face (surface, v, f);
}

static gboolean polygon_in_circle (GSList * poly,
                                   GtsPoint * p1, 
                                   GtsPoint * p2,
                                   GtsPoint * p3)
{
  GtsVertex * v1 = NULL, * v2 = NULL;

  while (poly) {
    GtsSegment * s = poly->data;
    GtsVertex * v;
    v = s->v1;
    if (v != v1 && v != v2 &&
        v != GTS_VERTEX (p1) &&
        v != GTS_VERTEX (p2) &&
        v != GTS_VERTEX (p3) &&
        gts_point_in_circle (GTS_POINT (v), p1, p2, p3) > 0.)
      return TRUE;
    v = s->v2;
    if (v != v1 && v != v2 &&
        v != GTS_VERTEX (p1) &&
        v != GTS_VERTEX (p2) &&
        v != GTS_VERTEX (p3) &&
        gts_point_in_circle (GTS_POINT (v), p1, p2, p3) > 0.)
      return TRUE;
    v1 = s->v1;
    v2 = s->v2;
    poly = poly->next;
  }
  return FALSE;
}

static void triangulate_polygon (GSList * poly, 
                                 GtsSurface * surface,
                                 GtsFace * ref)
{
  GSList * i, * poly1, * poly2;
  GtsVertex * v1, * v2, * v3 = NULL;
  gboolean found = FALSE;
  GtsSegment * s, * s1, * s2;
  GtsEdge * e1, * e2;
  GtsFace * f;

  if (poly == NULL || poly->next == NULL) {
    g_slist_free (poly);
    return;
  }

  s = poly->data;
  s1 = poly->next->data;
  if (s->v1 == s1->v1 || s->v1 == s1->v2) {
    v1 = s->v2;
    v2 = s->v1;
  }
  else {
    g_assert (s->v2 == s1->v1 || s->v2 == s1->v2);
    v1 = s->v1;
    v2 = s->v2;
  }

  i = poly->next;
  v3 = v2;
  while (i && !found) {
    s1 = i->data;
    if (s1->v1 == v3)
      v3 = s1->v2;
    else {
      g_assert (s1->v2 == v3);
      v3 = s1->v1;
    }
    if (v3 != v1 &&
        gts_point_orientation (GTS_POINT (v1), 
                               GTS_POINT (v2), 
                               GTS_POINT (v3)) >= 0. &&
        !polygon_in_circle (poly, 
                            GTS_POINT (v1), 
                            GTS_POINT (v2), 
                            GTS_POINT (v3)))
      found = TRUE;
    else 
      i = i->next;
  }

  if (!found) {
    g_slist_free (poly);
    return;
  }

  s1 = gts_vertices_are_connected (v2, v3);
  if (!GTS_IS_EDGE (s1))
    e1 = gts_edge_new (surface->edge_class, v2, v3);
  else
    e1 = GTS_EDGE (s1);
  s2 = gts_vertices_are_connected (v3, v1);
  if (!GTS_IS_EDGE (s2))
    e2 = gts_edge_new (surface->edge_class, v3, v1);
  else
    e2 = GTS_EDGE (s2);
  f = gts_face_new (surface->face_class, GTS_EDGE (s), e1, e2);
  gts_object_attributes (GTS_OBJECT (f), GTS_OBJECT (ref));
  gts_surface_add_face (surface, f);

  poly1 = poly->next;
  g_slist_free_1 (poly);
  if (i->next && e2 != i->next->data)
    poly2 = g_slist_prepend (i->next, e2);
  else
    poly2 = i->next;
  if (e1 != i->data)
    i->next = g_slist_prepend (NULL, e1);
  else
    i->next = NULL;

 triangulate_polygon (poly1, surface, ref);
 triangulate_polygon (poly2, surface, ref);
}

void gts_delaunay_remove_vertex (GtsSurface * surface, GtsVertex * v)
{
  GSList * triangles, * i;
  GtsFace * ref = NULL;

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

  i = triangles = gts_vertex_triangles (v, NULL);
  while (i && !ref) {
    if (GTS_IS_FACE (i->data) &&
        gts_face_has_parent_surface (i->data, surface))
      ref = i->data;
    i = i->next;
  }
  if (!ref) {
    g_slist_free (triangles);
    g_return_if_fail (ref);
  }
  triangulate_polygon (gts_vertex_fan_oriented (v, surface), surface, ref);
  i = triangles;
  while (i) {
    if (GTS_IS_FACE (i->data) &&
        gts_face_has_parent_surface (i->data, surface))
      gts_surface_remove_face (surface, i->data);
    i = i->next;
  }
  g_slist_free (triangles);
}

#define NEXT_CUT(edge, edge1, list) { next = neighbor (f, edge, surface);\
                                      remove_triangles (e, surface);\
                                      if (!constraint && !e->triangles)\
                                        gts_object_destroy (GTS_OBJECT (e));\
                                      g_assert (next);\
                                      *list = g_slist_prepend (*list, edge1);\
                                      return g_slist_concat (constraint,\
                                        remove_intersected_edge (s, edge,\
                                               next, surface, left, right));\
                                    }

static void remove_triangles (GtsEdge * e, GtsSurface * s)
{
  GSList * i = e->triangles;

  while (i) {
    GSList * next = i->next;

    if (GTS_IS_FACE (i->data) && gts_face_has_parent_surface (i->data, s))
      gts_surface_remove_face (s, i->data);
    i = next;
  }
}

static GSList * 
remove_intersected_edge (GtsSegment * s,
                         GtsEdge * e,
                         GtsFace * f,
                         GtsSurface * surface,
                         GSList ** left, GSList ** right)
{
  GtsVertex * v1, * v2, * v3;
  GtsEdge * e1, * e2;
  gdouble o1, o2;
  GtsFace * next;
  GSList * constraint = NULL;

  if (GTS_IS_CONSTRAINT (e))
    constraint = g_slist_prepend (NULL, e);

  gts_triangle_vertices_edges (GTS_TRIANGLE (f), e, 
                               &v1, &v2, &v3, &e, &e1, &e2);
  
  o1 = gts_point_orientation (GTS_POINT (v2), GTS_POINT (v3), 
                              GTS_POINT (s->v2));
  o2 = gts_point_orientation (GTS_POINT (v3), GTS_POINT (v1), 
                              GTS_POINT (s->v2));

  if (o1 == 0. && o2 == 0.) {
/*    if(o2 != 0.) {
      fprintf(stderr, "o1 = %f o2 = %f\n", o1, o2);
      fprintf(stderr, "v1 = %f, %f\n", GTS_POINT(v1)->x, GTS_POINT(v1)->y);
      fprintf(stderr, "v2 = %f, %f\n", GTS_POINT(v2)->x, GTS_POINT(v2)->y);
      fprintf(stderr, "v3 = %f, %f\n", GTS_POINT(v3)->x, GTS_POINT(v3)->y);
      fprintf(stderr, "s->v2 = %f, %f\n", GTS_POINT(s->v2)->x, GTS_POINT(s->v2)->y);

      g_assert (o2 == 0.);
    }*/
 //   if(o2 == 0.) {
      remove_triangles (e, surface);
      if (!constraint && !e->triangles)
        gts_object_destroy (GTS_OBJECT (e));
      *left = g_slist_prepend (*left, e2);
      *right = g_slist_prepend (*right, e1);
//    }
  }
  else if (o1 > 0.) {
    g_assert (o2 <= 0.);
    NEXT_CUT (e2, e1, right)
  }
  else if (o2 >= 0.)
    NEXT_CUT (e1, e2, left)
  else {
    gdouble o3 = gts_point_orientation (GTS_POINT (s->v1), GTS_POINT (s->v2),
                                        GTS_POINT (v3));
    if (o3 > 0.)
      NEXT_CUT (e1, e2, left)
    else
      NEXT_CUT (e2, e1, right)
  }
  return constraint;
}

static GSList * 
remove_intersected_vertex (GtsSegment * s,
                           GtsVertex * v,
                           GtsSurface * surface,
                           GSList ** left,
                           GSList ** right,
                           GtsFace ** ref)
{
  GSList * triangles = gts_vertex_triangles (v, NULL);
  GSList * i;

  i = triangles;
  while (i) {
    GtsTriangle * t = i->data;
    if (GTS_IS_FACE (t) && 
        gts_face_has_parent_surface (GTS_FACE (t), surface)) {
      GtsVertex * v1, * v2, * v3;
      gdouble o1, o2;

      gts_triangle_vertices (t, &v1, &v2, &v3);
      if (v == v2) {
        v2 = v3;
        v3 = v1;
      }
      else if (v == v3) {
        v3 = v2;
        v2 = v1;        
      }
      else
        g_assert (v == v1);

      if ((o1 = gts_point_orientation (GTS_POINT (v), GTS_POINT (v2),
                                       GTS_POINT (s->v2))) >= 0. &&
          (o2 = gts_point_orientation (GTS_POINT (v3), GTS_POINT (v),
                                       GTS_POINT (s->v2))) >= 0.) {
        gdouble o3 = gts_point_orientation (GTS_POINT (v2), GTS_POINT (v3),
                                            GTS_POINT (s->v2));
        GtsEdge * e = gts_triangle_edge_opposite (t, v);
        GtsEdge * e1, * e2;
        GtsFace * next = neighbor (GTS_FACE (t), e, surface);

        *ref = GTS_FACE (t);
        gts_triangle_vertices_edges (t, e, &v2, &v3, &v, &e, &e2, &e1);

        g_slist_free (triangles);

        if (o3 >= 0.) /* @s->v2 is inside (or on the edge) of t */
          return NULL;

        gts_allow_floating_faces = TRUE;
        gts_surface_remove_face (surface, GTS_FACE (t));
        gts_allow_floating_faces = FALSE;

        *left = g_slist_prepend (*left, e2);
        *right = g_slist_prepend (*right, e1);

        g_assert (next);
        return remove_intersected_edge (s, e, next, surface, left, right);
      }
    }
    i = i->next;
  }

  g_assert_not_reached ();
  return NULL;
}

GSList * gts_delaunay_add_constraint (GtsSurface * surface,
                                      GtsConstraint * c)
{
  GSList * constraints;
  GtsVertex * v1; //, * v2;
  GSList * left = NULL, * right = NULL;
  GtsFace * ref = NULL;

  g_return_val_if_fail (surface != NULL, NULL);
  g_return_val_if_fail (c != NULL, NULL);
  g_return_val_if_fail (GTS_IS_CONSTRAINT (c), NULL);
  
  v1 = GTS_SEGMENT (c)->v1;
  //v2 = GTS_SEGMENT (c)->v2;
  
  gts_allow_floating_edges = TRUE;
  constraints = remove_intersected_vertex (GTS_SEGMENT (c), v1, surface,
                                           &left, &right, &ref);
  gts_allow_floating_edges = FALSE;
#if 1
  triangulate_polygon (g_slist_prepend (g_slist_reverse (right), c), 
                       surface, ref);
  triangulate_polygon (g_slist_prepend (left, c), 
                       surface, ref);
#else
  right = g_slist_prepend (g_slist_reverse (right), c);
  left = g_slist_prepend (left, c);
  {
    FILE * fp0 = fopen ("hole", "wt");
    FILE * fp1 = fopen ("right", "wt");
    FILE * fp2 = fopen ("left", "wt");
    GSList * i = left;

    gts_surface_write (surface, fp0);
    fclose (fp0);
 
    fprintf (fp2, "LIST {\n");
    while (i) {
      GtsSegment * s = i->data;
      fprintf (fp2, 
               "# %p: %p->%p\n"
               "VECT 1 2 0 2 0 %g %g 0 %g %g 0\n",
               s, s->v1, s->v2,
               GTS_POINT (s->v1)->x, GTS_POINT (s->v1)->y,
               GTS_POINT (s->v2)->x, GTS_POINT (s->v2)->y);
      i = i->next;
    }
    fprintf (fp2, "}\n");
    fprintf (fp1, "LIST {\n");
    i = right;
    while (i) {
      GtsSegment * s = i->data;
      fprintf (fp1, 
               "# %p: %p->%p\n"
               "VECT 1 2 0 2 0 %g %g 0 %g %g 0\n",
               s, s->v1, s->v2,
               GTS_POINT (s->v1)->x, GTS_POINT (s->v1)->y,
               GTS_POINT (s->v2)->x, GTS_POINT (s->v2)->y);
      i = i->next;
    }
    fprintf (fp1, "}\n");
    fclose (fp1);
    fclose (fp2);
  }
  triangulate_polygon (right, surface);
  triangulate_polygon (left, surface);
#endif
  if (ref && !ref->surfaces) {
    gts_allow_floating_edges = TRUE;
    gts_object_destroy (GTS_OBJECT (ref));
    gts_allow_floating_edges = FALSE;
  }
  return constraints;
}

static void delaunay_check (GtsTriangle * t, gpointer * data)
{
  GtsSurface * surface = data[0];
  GtsFace ** face = data[1];

  if (*face == NULL) {
    GSList * i, * list;
    GtsVertex * v1, * v2, * v3;

    gts_triangle_vertices (t, &v1, &v2, &v3);
    list = gts_vertex_neighbors (v1, NULL, surface);
    list = gts_vertex_neighbors (v2, list, surface);
    list = gts_vertex_neighbors (v3, list, surface);
    i = list;
    while (i && *face == NULL) {
      GtsVertex * v = i->data;
      if (v != v1 && v != v2 && v != v3 &&
          gts_point_in_circle (GTS_POINT (v), 
                               GTS_POINT (v1),
                               GTS_POINT (v2),  
                               GTS_POINT (v3)) > 0.)
        *face = GTS_FACE (t);
      i = i->next;
    }
    g_slist_free (list);
  }
}

GtsFace * gts_delaunay_check (GtsSurface * surface)
{
  GtsFace * face = NULL;
  gpointer data[2];

  g_return_val_if_fail (surface != NULL, FALSE);

  data[0] = surface;
  data[1] = &face;
  gts_surface_foreach_face (surface, (GtsFunc) delaunay_check, data);

  return face;
}

void gts_delaunay_remove_hull (GtsSurface * surface)
{
  GSList * boundary;

  g_return_if_fail (surface != NULL);

  boundary = gts_surface_boundary (surface);
  gts_allow_floating_edges = TRUE;
  while (boundary) {
    GSList * i = boundary;
    GtsEdge * e = i->data;

    boundary = i->next;
    g_slist_free_1 (i);
    if (!GTS_IS_CONSTRAINT (e)) {
      GtsTriangle * t = GTS_TRIANGLE (gts_edge_is_boundary (e, surface));

      if (t != NULL) {
        if (t->e1 != e && !GTS_IS_CONSTRAINT (t->e1) &&
            !gts_edge_is_boundary (t->e1, surface))
          boundary = g_slist_prepend (boundary, t->e1);
        if (t->e2 != e && !GTS_IS_CONSTRAINT (t->e2) &&
            !gts_edge_is_boundary (t->e2, surface))
          boundary = g_slist_prepend (boundary, t->e2);
        if (t->e3 != e && !GTS_IS_CONSTRAINT (t->e3) &&
            !gts_edge_is_boundary (t->e3, surface))
          boundary = g_slist_prepend (boundary, t->e3);
        gts_surface_remove_face (surface, GTS_FACE (t));
      }
      if (!e->triangles)
        gts_object_destroy (GTS_OBJECT (e));
    }
  }
  gts_allow_floating_edges = FALSE;
}

/* GtsListFace: Object */

static void gts_list_face_destroy (GtsObject * object)
{
  g_slist_free (GTS_LIST_FACE (object)->points);

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

static void gts_list_face_class_init (GtsFaceClass * klass)
{
  GTS_OBJECT_CLASS (klass)->destroy = gts_list_face_destroy;
}

GtsFaceClass * gts_list_face_class (void)
{
  static GtsFaceClass * klass = NULL;

  if (klass == NULL) {
    GtsObjectClassInfo gts_list_face_info = {
      "GtsListFace",
      sizeof (GtsListFace),
      sizeof (GtsFaceClass),
      (GtsObjectClassInitFunc) gts_list_face_class_init,
      (GtsObjectInitFunc) NULL,
      (GtsArgSetFunc) NULL,
      (GtsArgGetFunc) NULL
    };
    klass = gts_object_class_new (GTS_OBJECT_CLASS (gts_face_class ()),
                                  &gts_list_face_info);
  }

  return klass;
}