triangle.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.
 */

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

static void triangle_destroy (GtsObject * object)
{
  GtsTriangle * triangle = GTS_TRIANGLE (object);
  GtsEdge * e1 = triangle->e1;
  GtsEdge * e2 = triangle->e2;
  GtsEdge * e3 = triangle->e3;

  e1->triangles = g_slist_remove (e1->triangles, triangle);
  if (!GTS_OBJECT_DESTROYED (e1) &&
      !gts_allow_floating_edges && e1->triangles == NULL)
    gts_object_destroy (GTS_OBJECT (e1));
  
  e2->triangles = g_slist_remove (e2->triangles, triangle);
  if (!GTS_OBJECT_DESTROYED (e2) &&
      !gts_allow_floating_edges && e2->triangles == NULL)
    gts_object_destroy (GTS_OBJECT (e2));
  
  e3->triangles = g_slist_remove (e3->triangles, triangle);
  if (!GTS_OBJECT_DESTROYED (e3) &&
      !gts_allow_floating_edges && e3->triangles == NULL)
    gts_object_destroy (GTS_OBJECT (e3));

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

static void triangle_class_init (GtsObjectClass * klass)
{
  klass->destroy = triangle_destroy;
}

static void triangle_init (GtsTriangle * triangle)
{
  triangle->e1 = triangle->e2 = triangle->e3 = NULL;
}

GtsTriangleClass * gts_triangle_class (void)
{
  static GtsTriangleClass * klass = NULL;

  if (klass == NULL) {
    GtsObjectClassInfo triangle_info = {
      "GtsTriangle",
      sizeof (GtsTriangle),
      sizeof (GtsTriangleClass),
      (GtsObjectClassInitFunc) triangle_class_init,
      (GtsObjectInitFunc) triangle_init,
      (GtsArgSetFunc) NULL,
      (GtsArgGetFunc) NULL
    };
    klass = gts_object_class_new (gts_object_class (), 
                                  &triangle_info);
  }

  return klass;
}

void gts_triangle_set (GtsTriangle * triangle, 
                       GtsEdge * e1, 
                       GtsEdge * e2,
                       GtsEdge * e3)
{
  g_return_if_fail (e1 != NULL);
  g_return_if_fail (e2 != NULL);
  g_return_if_fail (e3 != NULL);
  g_return_if_fail (e1 != e2 && e1 != e3 && e2 != e3);

  triangle->e1 = e1;
  triangle->e2 = e2;
  triangle->e3 = e3;

  if (GTS_SEGMENT (e1)->v1 == GTS_SEGMENT (e2)->v1)
    g_return_if_fail (gts_segment_connect (GTS_SEGMENT (e3), 
                                           GTS_SEGMENT (e1)->v2, 
                                           GTS_SEGMENT (e2)->v2));
  else if (GTS_SEGMENT (e1)->v2 == GTS_SEGMENT (e2)->v1)
    g_return_if_fail (gts_segment_connect (GTS_SEGMENT (e3), 
                                           GTS_SEGMENT (e1)->v1, 
                                           GTS_SEGMENT (e2)->v2));
  else if (GTS_SEGMENT (e1)->v2 == GTS_SEGMENT (e2)->v2)
    g_return_if_fail (gts_segment_connect (GTS_SEGMENT (e3), 
                                           GTS_SEGMENT (e1)->v1, 
                                           GTS_SEGMENT (e2)->v1));
  else if (GTS_SEGMENT (e1)->v1 == GTS_SEGMENT (e2)->v2)
    g_return_if_fail (gts_segment_connect (GTS_SEGMENT (e3), 
                                           GTS_SEGMENT (e1)->v2, 
                                           GTS_SEGMENT (e2)->v1));
  else
    g_assert_not_reached ();

  e1->triangles = g_slist_prepend (e1->triangles, triangle);
  e2->triangles = g_slist_prepend (e2->triangles, triangle);
  e3->triangles = g_slist_prepend (e3->triangles, triangle);
}

GtsTriangle * gts_triangle_new (GtsTriangleClass * klass,
                                GtsEdge * e1,
                                GtsEdge * e2,
                                GtsEdge * e3)
{
  GtsTriangle * t;

  t = GTS_TRIANGLE (gts_object_new (GTS_OBJECT_CLASS (klass)));
  gts_triangle_set (t, e1, e2, e3);
  
  return t;
}

GtsVertex * gts_triangle_vertex_opposite (GtsTriangle * t, GtsEdge * e)
{
  g_return_val_if_fail (t != NULL, NULL);
  g_return_val_if_fail (e != NULL, NULL);

  if (t->e1 == e) {
    GtsVertex * v = GTS_SEGMENT (t->e2)->v1;
    if (v != GTS_SEGMENT (e)->v1 && v != GTS_SEGMENT (e)->v2)
      return v;
    return GTS_SEGMENT (t->e2)->v2;
  }
  if (t->e2 == e) {
    GtsVertex * v = GTS_SEGMENT (t->e1)->v1;
    if (v != GTS_SEGMENT (e)->v1 && v != GTS_SEGMENT (e)->v2)
      return v;
    return GTS_SEGMENT (t->e1)->v2;
  }
  if (t->e3 == e) {
    GtsVertex * v = GTS_SEGMENT (t->e2)->v1;
    if (v != GTS_SEGMENT (e)->v1 && v != GTS_SEGMENT (e)->v2)
      return v;
    return GTS_SEGMENT (t->e2)->v2;
  }
  g_assert_not_reached ();
  return NULL;
}

GtsEdge * gts_triangle_edge_opposite (GtsTriangle * t, GtsVertex * v)
{
  GtsSegment * s1, * s2, * s3;

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

  s1 = GTS_SEGMENT (t->e1);
  s2 = GTS_SEGMENT (t->e2);

  if (s1->v1 != v && s1->v2 != v) {
    if (s2->v1 != v && s2->v2 != v)
      return NULL;
    return t->e1;
  }
  if (s2->v1 != v && s2->v2 != v)
    return t->e2;
  s3 = GTS_SEGMENT (t->e3);
  g_assert (s3->v1 != v && s3->v2 != v);
  return t->e3;
}

gdouble gts_triangles_angle (GtsTriangle * t1,
                             GtsTriangle * t2)
{
  gdouble nx1, ny1, nz1, nx2, ny2, nz2;
  gdouble pvx, pvy, pvz;
  gdouble theta;

  g_return_val_if_fail (t1 != NULL && t2 != NULL, 0.0);

  gts_triangle_normal (t1, &nx1, &ny1, &nz1);
  gts_triangle_normal (t2, &nx2, &ny2, &nz2);

  pvx = ny1*nz2 - nz1*ny2;
  pvy = nz1*nx2 - nx1*nz2;
  pvz = nx1*ny2 - ny1*nx2;

  theta = atan2 (sqrt (pvx*pvx + pvy*pvy + pvz*pvz), 
                 nx1*nx2 + ny1*ny2 + nz1*nz2) - M_PI;
  return theta < - M_PI ? theta + 2.*M_PI : theta;
}

gboolean gts_triangles_are_compatible (GtsTriangle * t1, 
                                       GtsTriangle * t2,
                                       GtsEdge * e)
{
  GtsEdge * e1 = NULL, * e2 = NULL;

  g_return_val_if_fail (t1 != NULL, FALSE);
  g_return_val_if_fail (t2 != NULL, FALSE);
  g_return_val_if_fail (e != NULL, FALSE);

  if (t1->e1 == e) e1 = t1->e2;
  else if (t1->e2 == e) e1 = t1->e3;
  else if (t1->e3 == e) e1 = t1->e1;
  else
    g_assert_not_reached ();
  if (t2->e1 == e) e2 = t2->e2;
  else if (t2->e2 == e) e2 = t2->e3;
  else if (t2->e3 == e) e2 = t2->e1;
  else
    g_assert_not_reached ();
  if (GTS_SEGMENT (e1)->v1 == GTS_SEGMENT (e2)->v1 || 
      GTS_SEGMENT (e1)->v1 == GTS_SEGMENT (e2)->v2 || 
      GTS_SEGMENT (e1)->v2 == GTS_SEGMENT (e2)->v1 || 
      GTS_SEGMENT (e1)->v2 == GTS_SEGMENT (e2)->v2)
    return FALSE;
  return TRUE;
}

gdouble gts_triangle_area (GtsTriangle * t)
{
  gdouble x, y, z;
  
  g_return_val_if_fail (t != NULL, 0.0);
  
  gts_triangle_normal (t, &x, &y, &z);
  
  return sqrt (x*x + y*y + z*z)/2.;
}

gdouble gts_triangle_perimeter (GtsTriangle * t)
{
  GtsVertex * v;

  g_return_val_if_fail (t != NULL, 0.0);

  v = gts_triangle_vertex (t);
  return 
    gts_point_distance (GTS_POINT (GTS_SEGMENT (t->e1)->v1), 
                        GTS_POINT (GTS_SEGMENT (t->e1)->v2)) +
    gts_point_distance (GTS_POINT (GTS_SEGMENT (t->e1)->v1), 
                        GTS_POINT (v)) +
    gts_point_distance (GTS_POINT (GTS_SEGMENT (t->e1)->v2), 
                        GTS_POINT (v));
}

/* perimeter of the equilateral triangle of area unity */
#define GOLDEN_PERIMETER 4.5590141139 

gdouble gts_triangle_quality (GtsTriangle * t)
{
  gdouble perimeter;

  g_return_val_if_fail (t != NULL, 0.0);

  perimeter = gts_triangle_perimeter (t);
  return perimeter > 0.0 ?
    GOLDEN_PERIMETER*sqrt (gts_triangle_area (t))/perimeter :
    0.0;
}

void gts_triangle_normal (GtsTriangle * t, 
                          gdouble * x, 
                          gdouble * y, 
                          gdouble * z)
{
  GtsVertex * v1, * v2 = NULL, * v3 = NULL;
  GtsPoint * p1, * p2, * p3;
  gdouble x1, y1, z1, x2, y2, z2;

  g_return_if_fail (t != NULL);

  v1 = GTS_SEGMENT (t->e1)->v1;
  if (GTS_SEGMENT (t->e1)->v1 == GTS_SEGMENT (t->e2)->v1) {
    v2 = GTS_SEGMENT (t->e2)->v2;
    v3 = GTS_SEGMENT (t->e1)->v2;
  }
  else if (GTS_SEGMENT (t->e1)->v2 == GTS_SEGMENT (t->e2)->v2) {
    v2 = GTS_SEGMENT (t->e1)->v2;
    v3 = GTS_SEGMENT (t->e2)->v1;
  }
  else if (GTS_SEGMENT (t->e1)->v1 == GTS_SEGMENT (t->e2)->v2) {
    v2 = GTS_SEGMENT (t->e2)->v1;
    v3 = GTS_SEGMENT (t->e1)->v2;
  }
  else if (GTS_SEGMENT (t->e1)->v2 == GTS_SEGMENT (t->e2)->v1) {
    v2 = GTS_SEGMENT (t->e1)->v2;
    v3 = GTS_SEGMENT (t->e2)->v2;
  }
  else {
    fprintf (stderr, "t: %p t->e1: %p t->e2: %p t->e3: %p t->e1->v1: %p t->e1->v2: %p t->e2->v1: %p t->e2->v2: %p t->e3->v1: %p t->e3->v2: %p\n",
         t, t->e1, t->e2, 
         t->e3, GTS_SEGMENT (t->e1)->v1, GTS_SEGMENT (t->e1)->v2, 
         GTS_SEGMENT (t->e2)->v1, GTS_SEGMENT (t->e2)->v2, 
         GTS_SEGMENT (t->e3)->v1, GTS_SEGMENT (t->e3)->v2);
    g_assert_not_reached ();
  }

  p1 = GTS_POINT (v1);
  p2 = GTS_POINT (v2);
  p3 = GTS_POINT (v3);

  x1 = p2->x - p1->x;
  y1 = p2->y - p1->y;
  z1 = p2->z - p1->z;

  x2 = p3->x - p1->x;
  y2 = p3->y - p1->y;
  z2 = p3->z - p1->z;

  *x = y1*z2 - z1*y2;
  *y = z1*x2 - x1*z2;
  *z = x1*y2 - y1*x2;
}

gdouble gts_triangle_orientation (GtsTriangle * t)
{
  GtsVertex * v1, * v2 = NULL, * v3 = NULL;

  g_return_val_if_fail (t != NULL, 0.0);

  v1 = GTS_SEGMENT (t->e1)->v1;
  if (GTS_SEGMENT (t->e1)->v1 == GTS_SEGMENT (t->e2)->v1) {
    v2 = GTS_SEGMENT (t->e2)->v2;
    v3 = GTS_SEGMENT (t->e1)->v2;
  }
  else if (GTS_SEGMENT (t->e1)->v2 == GTS_SEGMENT (t->e2)->v2) {
    v2 = GTS_SEGMENT (t->e1)->v2;
    v3 = GTS_SEGMENT (t->e2)->v1;
  }
  else if (GTS_SEGMENT (t->e1)->v1 == GTS_SEGMENT (t->e2)->v2) {
    v2 = GTS_SEGMENT (t->e2)->v1;
    v3 = GTS_SEGMENT (t->e1)->v2;
  }
  else if (GTS_SEGMENT (t->e1)->v2 == GTS_SEGMENT (t->e2)->v1) {
    v2 = GTS_SEGMENT (t->e1)->v2;
    v3 = GTS_SEGMENT (t->e2)->v2;
  }
  else
    g_assert_not_reached ();
  return gts_point_orientation (GTS_POINT (v1), 
                                GTS_POINT (v2), 
                                GTS_POINT (v3));
}

void gts_triangle_revert (GtsTriangle * t)
{
  GtsEdge * e;

  g_return_if_fail (t != NULL);

  e = t->e1;
  t->e1 = t->e2;
  t->e2 = e;
}

GSList * gts_triangles_from_edges (GSList * edges)
{
  GHashTable * hash;
  GSList * triangles = NULL, * i;

  hash = g_hash_table_new (NULL, NULL);
  i = edges;
  while (i) {
    GSList * j = GTS_EDGE (i->data)->triangles;
    while (j) {
      GtsTriangle * t = j->data;
      if (g_hash_table_lookup (hash, t) == NULL) {
        triangles = g_slist_prepend (triangles, t);
        g_hash_table_insert (hash, t, i);
      }
      j = j->next;
    }
    i = i->next;
  }
  g_hash_table_destroy (hash);

  return triangles;
}

void gts_triangle_vertices_edges (GtsTriangle * t, 
                                  GtsEdge * e,
                                  GtsVertex ** v1, 
                                  GtsVertex ** v2, 
                                  GtsVertex ** v3,
                                  GtsEdge ** e1,
                                  GtsEdge ** e2,
                                  GtsEdge ** e3)
{
  GtsEdge * ee1, * ee2;

  g_return_if_fail (t != NULL);
  
  if (e == t->e1 || e == NULL) {
    *e1 = ee1 = t->e1; *e2 = ee2 = t->e2; *e3 = t->e3;    
  }
  else if (e == t->e2) {
    *e1 = ee1 = e; *e2 = ee2 = t->e3; *e3 = t->e1;
  }
  else if (e == t->e3) {
    *e1 = ee1 = e; *e2 = ee2 = t->e1; *e3 = t->e2;
  }
  else {
    g_assert_not_reached ();
    ee1 = ee2 = NULL; /* to avoid complaints from the compiler */
  }
  if (GTS_SEGMENT (ee1)->v2 == GTS_SEGMENT (ee2)->v1) {
    *v1 = GTS_SEGMENT (ee1)->v1; 
    *v2 = GTS_SEGMENT (ee1)->v2; 
    *v3 = GTS_SEGMENT (ee2)->v2;
  }
  else if (GTS_SEGMENT (ee1)->v2 == GTS_SEGMENT (ee2)->v2) {
    *v1 = GTS_SEGMENT (ee1)->v1; 
    *v2 = GTS_SEGMENT (ee1)->v2; 
    *v3 = GTS_SEGMENT (ee2)->v1;
  }
  else if (GTS_SEGMENT (ee1)->v1 == GTS_SEGMENT (ee2)->v1) {
    *v1 = GTS_SEGMENT (ee1)->v2; 
    *v2 = GTS_SEGMENT (ee1)->v1; 
    *v3 = GTS_SEGMENT (ee2)->v2;
  }
  else if (GTS_SEGMENT (ee1)->v1 == GTS_SEGMENT (ee2)->v2) {
    *v1 = GTS_SEGMENT (ee1)->v2; 
    *v2 = GTS_SEGMENT (ee1)->v1; 
    *v3 = GTS_SEGMENT (ee2)->v1;
  }
  else
    g_assert_not_reached ();
}

/* sqrt(3) */
#define SQRT3 1.73205080757

GtsTriangle * gts_triangle_enclosing (GtsTriangleClass * klass,
                                      GSList * points, gdouble scale)
{
  gdouble xmax, xmin, ymax, ymin;
  gdouble xo, yo, r;
  GtsVertex * v1, * v2, * v3;
  GtsEdge * e1, * e2, * e3;

  if (points == NULL)
    return NULL;
  
  xmax = xmin = GTS_POINT (points->data)->x;
  ymax = ymin = GTS_POINT (points->data)->y;
  points = points->next;
  while (points) {
    GtsPoint * p = points->data;
    if (p->x > xmax) xmax = p->x;
    else if (p->x < xmin) xmin = p->x;
    if (p->y > ymax) ymax = p->y;
    else if (p->y < ymin) ymin = p->y;    
    points = points->next;
  }
  xo = (xmax + xmin)/2.;
  yo = (ymax + ymin)/2.;
  r = scale*sqrt((xmax - xo)*(xmax - xo) + (ymax - yo)*(ymax - yo));
  if (r == 0.0) r = scale;
  v1 = gts_vertex_new (gts_vertex_class (),
                       xo + r*SQRT3, yo - r, 0.0);
  v2 = gts_vertex_new (gts_vertex_class (),
                       xo, yo + 2.*r, 0.0);
  v3 = gts_vertex_new (gts_vertex_class (),
                       xo - r*SQRT3, yo - r, 0.0);
  e1 = gts_edge_new (gts_edge_class (), v1, v2);
  e2 = gts_edge_new (gts_edge_class (), v2, v3);
  e3 = gts_edge_new (gts_edge_class (), v3, v1);
  return gts_triangle_new (gts_triangle_class (), e1, e2, e3);
}

guint gts_triangle_neighbor_number (GtsTriangle * t)
{
  GSList * i;
  guint nn = 0;
  GtsEdge * ee[4], ** e = ee;
  
  g_return_val_if_fail (t != NULL, 0);

  ee[0] = t->e1; ee[1] = t->e2; ee[2] = t->e3; ee[3] = NULL;
  while (*e) {
    i = (*e++)->triangles;
    while (i) {
      GtsTriangle * t1 = i->data;
      if (t1 != t)
        nn++;
      i = i->next;
    }
  }
  return nn;
}

GSList * gts_triangle_neighbors (GtsTriangle * t)
{
  GSList * i, * list = NULL;
  GtsEdge * ee[4], ** e = ee;
  
  g_return_val_if_fail (t != NULL, NULL);

  ee[0] = t->e1; ee[1] = t->e2; ee[2] = t->e3; ee[3] = NULL;
  while (*e) {
    i = (*e++)->triangles;
    while (i) {
      GtsTriangle * t1 = i->data;
      if (t1 != t)
        list = g_slist_prepend (list, t1);
      i = i->next;
    }
  }
  return list;
}

GtsEdge * gts_triangles_common_edge (GtsTriangle * t1,
                                     GtsTriangle * t2)
{
  g_return_val_if_fail (t1 != NULL, NULL);
  g_return_val_if_fail (t2 != NULL, NULL);

  if (t1->e1 == t2->e1 || t1->e1 == t2->e2 || t1->e1 == t2->e3)
    return t1->e1;
  if (t1->e2 == t2->e1 || t1->e2 == t2->e2 || t1->e2 == t2->e3)
    return t1->e2;
  if (t1->e3 == t2->e1 || t1->e3 == t2->e2 || t1->e3 == t2->e3)
    return t1->e3;
  return NULL;
}

GtsTriangle * gts_triangle_is_duplicate (GtsTriangle * t)
{
  GSList * i;
  GtsEdge * e2, * e3;

  g_return_val_if_fail (t != NULL, NULL);

  e2 = t->e2;
  e3 = t->e3;
  i = t->e1->triangles;
  while (i) {
    GtsTriangle * t1 = i->data;
    if (t1 != t && 
        (t1->e1 == e2 || t1->e2 == e2 || t1->e3 == e2) &&
        (t1->e1 == e3 || t1->e2 == e3 || t1->e3 == e3))
      return t1;
    i = i->next;
  }
  
  return NULL;
}

GtsTriangle * gts_triangle_use_edges (GtsEdge * e1,
                                      GtsEdge * e2,
                                      GtsEdge * e3)
{
  GSList * i;
  
  g_return_val_if_fail (e1 != NULL, NULL);
  g_return_val_if_fail (e2 != NULL, NULL);
  g_return_val_if_fail (e3 != NULL, NULL);

  i = e1->triangles;
  while (i) {
    GtsTriangle * t = i->data;
    if ((t->e1 == e2 && (t->e2 == e3 || t->e3 == e3)) ||
        (t->e2 == e2 && (t->e1 == e3 || t->e3 == e3)) ||
        (t->e3 == e2 && (t->e1 == e3 || t->e2 == e3)))
      return t;
    i = i->next;
  }
  
  return NULL;
}

gboolean gts_triangle_is_ok (GtsTriangle * t)
{
  g_return_val_if_fail (t != NULL, FALSE);
  g_return_val_if_fail (t->e1 != NULL, FALSE);
  g_return_val_if_fail (t->e2 != NULL, FALSE);
  g_return_val_if_fail (t->e3 != NULL, FALSE);
  g_return_val_if_fail (t->e1 != t->e2 && t->e1 != t->e3 && t->e2 != t->e3, 
                        FALSE);
  g_return_val_if_fail (gts_segments_touch (GTS_SEGMENT (t->e1), 
                                            GTS_SEGMENT (t->e2)),
                        FALSE);
  g_return_val_if_fail (gts_segments_touch (GTS_SEGMENT (t->e1), 
                                            GTS_SEGMENT (t->e3)), 
                        FALSE);
  g_return_val_if_fail (gts_segments_touch (GTS_SEGMENT (t->e2), 
                                            GTS_SEGMENT (t->e3)), 
                        FALSE);
  g_return_val_if_fail (GTS_SEGMENT (t->e1)->v1 != GTS_SEGMENT (t->e1)->v2, 
                        FALSE);
  g_return_val_if_fail (GTS_SEGMENT (t->e2)->v1 != GTS_SEGMENT (t->e2)->v2, 
                        FALSE);
  g_return_val_if_fail (GTS_SEGMENT (t->e3)->v1 != GTS_SEGMENT (t->e3)->v2, 
                        FALSE);
  g_return_val_if_fail (GTS_OBJECT (t)->reserved == NULL, FALSE);
  g_return_val_if_fail (!gts_triangle_is_duplicate (t), FALSE);
  return TRUE;
}

void gts_triangle_vertices (GtsTriangle * t,
                            GtsVertex ** v1, GtsVertex ** v2, GtsVertex ** v3)
{
  GtsSegment * e1, * e2;

  g_return_if_fail (t != NULL);
  g_return_if_fail (v1 != NULL && v2 != NULL && v3 != NULL);

  e1 = GTS_SEGMENT (t->e1);
  e2 = GTS_SEGMENT (t->e2);

  if (GTS_SEGMENT (e1)->v2 == GTS_SEGMENT (e2)->v1) {
    *v1 = GTS_SEGMENT (e1)->v1; 
    *v2 = GTS_SEGMENT (e1)->v2; 
    *v3 = GTS_SEGMENT (e2)->v2;
  }
  else if (GTS_SEGMENT (e1)->v2 == GTS_SEGMENT (e2)->v2) {
    *v1 = GTS_SEGMENT (e1)->v1; 
    *v2 = GTS_SEGMENT (e1)->v2; 
    *v3 = GTS_SEGMENT (e2)->v1;
  }
  else if (GTS_SEGMENT (e1)->v1 == GTS_SEGMENT (e2)->v1) {
    *v1 = GTS_SEGMENT (e1)->v2; 
    *v2 = GTS_SEGMENT (e1)->v1; 
    *v3 = GTS_SEGMENT (e2)->v2;
  }
  else {
    *v1 = GTS_SEGMENT (e1)->v2; 
    *v2 = GTS_SEGMENT (e1)->v1; 
    *v3 = GTS_SEGMENT (e2)->v1;
  }
}

GtsPoint * gts_triangle_circumcircle_center (GtsTriangle * t,
                                             GtsPointClass * point_class)
{
  GtsVertex * v1, * v2, * v3;
  gdouble xa, ya, xb, yb, xc, yc;
  gdouble xd, yd, xe, ye;
  gdouble xad, yad, xae, yae;
  gdouble det;
  
  g_return_val_if_fail (t != NULL, NULL);
  g_return_val_if_fail (point_class != NULL, NULL);

  gts_triangle_vertices (t, &v1, &v2, &v3);

  xa = GTS_POINT (v1)->x; ya = GTS_POINT (v1)->y;
  xb = GTS_POINT (v2)->x; yb = GTS_POINT (v2)->y;
  xc = GTS_POINT (v3)->x; yc = GTS_POINT (v3)->y;
  xd = (xa + xb)/2.; yd = (ya + yb)/2.;
  xe = (xa + xc)/2.; ye = (ya + yc)/2.;
  xad = xd - xa; yad = yd - ya;
  xae = xe - xa; yae = ye - ya;
  det = xad*yae - xae*yad;
  if (det == 0.)
    return NULL;
  return gts_point_new (point_class,
                        (yae*yad*(yd - ye) + xad*yae*xd - xae*yad*xe)/det,
                        -(xae*xad*(xd - xe) + yad*xae*yd - yae*xad*ye)/det,
                        0.);
}

/* square of the maximum area ratio admissible */
#define AREA_RATIO_MAX2 1e8

static gboolean points_are_folded (GtsPoint * A,
                                   GtsPoint * B,
                                   GtsPoint * C,
                                   GtsPoint * D,
                                   gdouble max)
{
  GtsVector AB, AC, AD;
  GtsVector n1, n2;
  gdouble nn1, nn2, n1n2;

  gts_vector_init (AB, A, B);
  gts_vector_init (AC, A, C);
  gts_vector_init (AD, A, D);
  gts_vector_cross (n1, AB, AC);
  gts_vector_cross (n2, AD, AB);

  nn1 = gts_vector_scalar (n1, n1);
  nn2 = gts_vector_scalar (n2, n2);
  if (nn1 >= AREA_RATIO_MAX2*nn2 || nn2 >= AREA_RATIO_MAX2*nn1)
    return TRUE;
  n1n2 = gts_vector_scalar (n1, n2);
  if (n1n2 > 0.)
    return FALSE;
  if (n1n2*n1n2/(nn1*nn2) > max)
    return TRUE;
  return FALSE;
}

static GtsVertex * triangle_use_vertices (GtsTriangle * t,
                                          GtsVertex * A, 
                                          GtsVertex * B)
{
  GtsVertex 
    * v1 = GTS_SEGMENT (t->e1)->v1, 
    * v2 = GTS_SEGMENT (t->e1)->v2, 
    * v3 = gts_triangle_vertex (t);

  if (v1 == A) {
    if (v2 == B)
      return v3;
    g_assert (v3 == B);
    return v2;
  }
  if (v2 == A) {
    if (v1 == B)
      return v3;
    g_assert (v3 == B);
    return v1;
  }
  if (v3 == A) {
    if (v1 == B)
      return v2;
    g_assert (v2 == B);
    return v1;
  }
  g_assert_not_reached ();
  return NULL;
}

gboolean gts_triangles_are_folded (GSList * triangles,
                                   GtsVertex * A, GtsVertex * B,
                                   gdouble max)
{
  GSList * i;

  g_return_val_if_fail (A != NULL, TRUE);
  g_return_val_if_fail (B != NULL, TRUE);

  i = triangles;
  while (i) {
    GtsVertex * C = triangle_use_vertices (i->data, A, B);
    GSList * j = i->next;    
    while (j) {
      GtsVertex * D = triangle_use_vertices (j->data, A, B);
      if (points_are_folded (GTS_POINT (A), 
                             GTS_POINT (B), 
                             GTS_POINT (C), 
                             GTS_POINT (D), 
                             max))
        return TRUE;
      j = j->next;
    }
    i = i->next;
  }
  return FALSE;
}

GtsObject * gts_triangle_is_stabbed (GtsTriangle * t, 
                                     GtsPoint * p,
                                     gdouble * orientation)
{
  GtsVertex * v1, * v2, * v3, * inverted = NULL;
  GtsEdge * e1, * e2, * e3, * tmp;
  gdouble o, o1, o2, o3;

  g_return_val_if_fail (t != NULL, NULL);
  g_return_val_if_fail (p != NULL, NULL);

  gts_triangle_vertices_edges (t, NULL, &v1, &v2, &v3, &e1, &e2, &e3);
  o = gts_point_orientation (GTS_POINT (v1), GTS_POINT (v2), GTS_POINT (v3));
  if (o == 0.)
    return NULL;
  if (o < 0.) {
    inverted = v1;
    v1 = v2;
    v2 = inverted;
    tmp = e2;
    e2 = e3;
    e3 = tmp;
  }
  o = gts_point_orientation_3d (GTS_POINT (v1),
                                GTS_POINT (v2), 
                                GTS_POINT (v3), 
                                p);
  if (o < 0.)
    return NULL;
  o1 = gts_point_orientation (GTS_POINT (v1), GTS_POINT (v2), p);
  if (o1 < 0.)
    return NULL;
  o2 = gts_point_orientation (GTS_POINT (v2), GTS_POINT (v3), p);
  if (o2 < 0.)
    return NULL;
  o3 = gts_point_orientation (GTS_POINT (v3), GTS_POINT (v1), p);
  if (o3 < 0.)
    return NULL;  
  if (orientation) *orientation = inverted ? -o : o;
  if (o1 == 0.) {
    if (o2 == 0.)
      return GTS_OBJECT (v2);
    if (o3 == 0.)
      return GTS_OBJECT (v1);
    return GTS_OBJECT (e1);
  }
  if (o2 == 0.) {
    if (o3 == 0.)
      return GTS_OBJECT (v3);
    return GTS_OBJECT (e2);
  }
  if (o3 == 0.)
    return GTS_OBJECT (e3);
  return GTS_OBJECT (t);
}

void gts_triangle_interpolate_height (GtsTriangle * t, GtsPoint * p)
{
  GtsPoint * p1, * p2, * p3;
  gdouble x1, x2, y1, y2, det;

  g_return_if_fail (t != NULL);
  g_return_if_fail (p != NULL);

  p1 = GTS_POINT (GTS_SEGMENT (t->e1)->v1);
  p2 = GTS_POINT (GTS_SEGMENT (t->e1)->v2);
  p3 = GTS_POINT (gts_triangle_vertex (t));

  x1 = p2->x - p1->x;
  y1 = p2->y - p1->y;
  x2 = p3->x - p1->x;
  y2 = p3->y - p1->y;
  det = x1*y2 - x2*y1;
  if (det == 0.)
    p->z = (p1->z + p2->z + p3->z)/3.;
  else {
    gdouble x = p->x - p1->x;
    gdouble y = p->y - p1->y;
    gdouble a = (x*y2 - y*x2)/det;
    gdouble b = (y*x1 - x*y1)/det;

    p->z = (1. - a - b)*p1->z + a*p2->z + b*p3->z;
  }
}