gts/triangle.c

   1 /* GTS - Library for the manipulation of triangulated surfaces
   2  * Copyright (C) 1999 Stéphane Popinet
   3  *
   4  * This library is free software; you can redistribute it and/or
   5  * modify it under the terms of the GNU Library General Public
   6  * License as published by the Free Software Foundation; either
   7  * version 2 of the License, or (at your option) any later version.
   8  *
   9  * This library is distributed in the hope that it will be useful,
  10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  12  * Library General Public License for more details.
  13  *
  14  * You should have received a copy of the GNU Library General Public
  15  * License along with this library; if not, write to the
  16  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  17  * Boston, MA 02111-1307, USA.
  18  */
  19
  20 #include <math.h>
  21 #include "gts.h"
  22
  23 static void triangle_destroy (GtsObject * object)
  24 {
  25   GtsTriangle * triangle = GTS_TRIANGLE (object);
  26   GtsEdge * e1 = triangle->e1;
  27   GtsEdge * e2 = triangle->e2;
  28   GtsEdge * e3 = triangle->e3;
  29
  30   e1->triangles = g_slist_remove (e1->triangles, triangle);
  31   if (!GTS_OBJECT_DESTROYED (e1) &&
  32       !gts_allow_floating_edges && e1->triangles == NULL)
  33     gts_object_destroy (GTS_OBJECT (e1));
  34
  35   e2->triangles = g_slist_remove (e2->triangles, triangle);
  36   if (!GTS_OBJECT_DESTROYED (e2) &&
  37       !gts_allow_floating_edges && e2->triangles == NULL)
  38     gts_object_destroy (GTS_OBJECT (e2));
  39
  40   e3->triangles = g_slist_remove (e3->triangles, triangle);
  41   if (!GTS_OBJECT_DESTROYED (e3) &&
  42       !gts_allow_floating_edges && e3->triangles == NULL)
  43     gts_object_destroy (GTS_OBJECT (e3));
  44
  45   (* GTS_OBJECT_CLASS (gts_triangle_class ())->parent_class->destroy) (object);
  46 }
  47
  48 static void triangle_class_init (GtsObjectClass * klass)
  49 {
  50   klass->destroy = triangle_destroy;
  51 }
  52
  53 static void triangle_init (GtsTriangle * triangle)
  54 {
  55   triangle->e1 = triangle->e2 = triangle->e3 = NULL;
  56 }
  57
  58 /**
  59  * gts_triangle_class:
  60  *
  61  * Returns: the #GtsTriangleClass.
  62  */
  63 GtsTriangleClass * gts_triangle_class (void)
  64 {
  65   static GtsTriangleClass * klass = NULL;
  66
  67   if (klass == NULL) {
  68     GtsObjectClassInfo triangle_info = {
  69       "GtsTriangle",
  70       sizeof (GtsTriangle),
  71       sizeof (GtsTriangleClass),
  72       (GtsObjectClassInitFunc) triangle_class_init,
  73       (GtsObjectInitFunc) triangle_init,
  74       (GtsArgSetFunc) NULL,
  75       (GtsArgGetFunc) NULL
  76     };
  77     klass = gts_object_class_new (gts_object_class (),
  78                                   &triangle_info);
  79   }
  80
  81   return klass;
  82 }
  83
  84 /**
  85  * gts_triangle_set:
  86  * @triangle: a #GtsTriangle.
  87  * @e1: a #GtsEdge.
  88  * @e2: another #GtsEdge touching @e1.
  89  * @e3: another #GtsEdge touching both @e1 and @e2.
  90  *
  91  * Sets the edge of @triangle to @e1, @e2 and @e3 while checking that they
  92  * define a valid triangle.
  93  */
  94 void gts_triangle_set (GtsTriangle * triangle,
  95                        GtsEdge * e1,
  96                        GtsEdge * e2,
  97                        GtsEdge * e3)
  98 {
  99   g_return_if_fail (e1 != NULL);
 100   g_return_if_fail (e2 != NULL);
 101   g_return_if_fail (e3 != NULL);
 102   g_return_if_fail (e1 != e2 && e1 != e3 && e2 != e3);
 103
 104   triangle->e1 = e1;
 105   triangle->e2 = e2;
 106   triangle->e3 = e3;
 107
 108   if (GTS_SEGMENT (e1)->v1 == GTS_SEGMENT (e2)->v1)
 109     g_return_if_fail (gts_segment_connect (GTS_SEGMENT (e3),
 110                                            GTS_SEGMENT (e1)->v2,
 111                                            GTS_SEGMENT (e2)->v2));
 112   else if (GTS_SEGMENT (e1)->v2 == GTS_SEGMENT (e2)->v1)
 113     g_return_if_fail (gts_segment_connect (GTS_SEGMENT (e3),
 114                                            GTS_SEGMENT (e1)->v1,
 115                                            GTS_SEGMENT (e2)->v2));
 116   else if (GTS_SEGMENT (e1)->v2 == GTS_SEGMENT (e2)->v2)
 117     g_return_if_fail (gts_segment_connect (GTS_SEGMENT (e3),
 118                                            GTS_SEGMENT (e1)->v1,
 119                                            GTS_SEGMENT (e2)->v1));
 120   else if (GTS_SEGMENT (e1)->v1 == GTS_SEGMENT (e2)->v2)
 121     g_return_if_fail (gts_segment_connect (GTS_SEGMENT (e3),
 122                                            GTS_SEGMENT (e1)->v2,
 123                                            GTS_SEGMENT (e2)->v1));
 124   else
 125     g_assert_not_reached ();
 126
 127   e1->triangles = g_slist_prepend (e1->triangles, triangle);
 128   e2->triangles = g_slist_prepend (e2->triangles, triangle);
 129   e3->triangles = g_slist_prepend (e3->triangles, triangle);
 130 }
 131
 132 /**
 133  * gts_triangle_new:
 134  * @klass: a #GtsTriangleClass.
 135  * @e1: a #GtsEdge.
 136  * @e2: another #GtsEdge touching @e1.
 137  * @e3: another #GtsEdge touching both @e1 and @e2.
 138  *
 139  * Returns: a new #GtsTriangle having @e1, @e2 and @e3 as edges.
 140  */
 141 GtsTriangle * gts_triangle_new (GtsTriangleClass * klass,
 142                                 GtsEdge * e1,
 143                                 GtsEdge * e2,
 144                                 GtsEdge * e3)
 145 {
 146   GtsTriangle * t;
 147
 148   t = GTS_TRIANGLE (gts_object_new (GTS_OBJECT_CLASS (klass)));
 149   gts_triangle_set (t, e1, e2, e3);
 150
 151   return t;
 152 }
 153
 154 /**
 155  * gts_triangle_vertex_opposite:
 156  * @t: a #GtsTriangle.
 157  * @e: a #GtsEdge used by @t.
 158  *
 159  * This function fails if @e is not an edge of @t.
 160  *
 161  * Returns: a #GtsVertex, vertex of @t which does not belong to @e.
 162  */
 163 GtsVertex * gts_triangle_vertex_opposite (GtsTriangle * t, GtsEdge * e)
 164 {
 165   g_return_val_if_fail (t != NULL, NULL);
 166   g_return_val_if_fail (e != NULL, NULL);
 167
 168   if (t->e1 == e) {
 169     GtsVertex * v = GTS_SEGMENT (t->e2)->v1;
 170     if (v != GTS_SEGMENT (e)->v1 && v != GTS_SEGMENT (e)->v2)
 171       return v;
 172     return GTS_SEGMENT (t->e2)->v2;
 173   }
 174   if (t->e2 == e) {
 175     GtsVertex * v = GTS_SEGMENT (t->e1)->v1;
 176     if (v != GTS_SEGMENT (e)->v1 && v != GTS_SEGMENT (e)->v2)
 177       return v;
 178     return GTS_SEGMENT (t->e1)->v2;
 179   }
 180   if (t->e3 == e) {
 181     GtsVertex * v = GTS_SEGMENT (t->e2)->v1;
 182     if (v != GTS_SEGMENT (e)->v1 && v != GTS_SEGMENT (e)->v2)
 183       return v;
 184     return GTS_SEGMENT (t->e2)->v2;
 185   }
 186   g_assert_not_reached ();
 187   return NULL;
 188 }
 189
 190 /**
 191  * gts_triangle_edge_opposite:
 192  * @t: a #GtsTriangle.
 193  * @v: a #GtsVertex of @t.
 194  *
 195  * Returns: the edge of @t opposite @v or %NULL if @v is not a vertice of @t.
 196  */
 197 GtsEdge * gts_triangle_edge_opposite (GtsTriangle * t, GtsVertex * v)
 198 {
 199   GtsSegment * s1, * s2, * s3;
 200
 201   g_return_val_if_fail (t != NULL, NULL);
 202   g_return_val_if_fail (v != NULL, NULL);
 203
 204   s1 = GTS_SEGMENT (t->e1);
 205   s2 = GTS_SEGMENT (t->e2);
 206
 207   if (s1->v1 != v && s1->v2 != v) {
 208     if (s2->v1 != v && s2->v2 != v)
 209       return NULL;
 210     return t->e1;
 211   }
 212   if (s2->v1 != v && s2->v2 != v)
 213     return t->e2;
 214   s3 = GTS_SEGMENT (t->e3);
 215   g_assert (s3->v1 != v && s3->v2 != v);
 216   return t->e3;
 217 }
 218
 219 /**
 220  * gts_triangles_angle:
 221  * @t1: a #GtsTriangle.
 222  * @t2: a #GtsTriangle.
 223  *
 224  * Returns: the value (in radians) of the angle between @t1 and @t2.
 225  */
 226 gdouble gts_triangles_angle (GtsTriangle * t1,
 227                              GtsTriangle * t2)
 228 {
 229   gdouble nx1, ny1, nz1, nx2, ny2, nz2;
 230   gdouble pvx, pvy, pvz;
 231   gdouble theta;
 232
 233   g_return_val_if_fail (t1 != NULL && t2 != NULL, 0.0);
 234
 235   gts_triangle_normal (t1, &nx1, &ny1, &nz1);
 236   gts_triangle_normal (t2, &nx2, &ny2, &nz2);
 237
 238   pvx = ny1*nz2 - nz1*ny2;
 239   pvy = nz1*nx2 - nx1*nz2;
 240   pvz = nx1*ny2 - ny1*nx2;
 241
 242   theta = atan2 (sqrt (pvx*pvx + pvy*pvy + pvz*pvz),
 243                  nx1*nx2 + ny1*ny2 + nz1*nz2) - M_PI;
 244   return theta < - M_PI ? theta + 2.*M_PI : theta;
 245 }
 246
 247 /**
 248  * gts_triangles_are_compatible:
 249  * @t1: a #GtsTriangle.
 250  * @t2: a #GtsTriangle.
 251  * @e: a #GtsEdge used by both @t1 and @t2.
 252  *
 253  * Checks if @t1 and @t2 have compatible orientations i.e. if @t1 and
 254  * @t2 can be part of the same surface without conflict in the surface
 255  * normal orientation.
 256  *
 257  * Returns: %TRUE if @t1 and @t2 are compatible, %FALSE otherwise.
 258  */
 259 gboolean gts_triangles_are_compatible (GtsTriangle * t1,
 260                                        GtsTriangle * t2,
 261                                        GtsEdge * e)
 262 {
 263   GtsEdge * e1 = NULL, * e2 = NULL;
 264
 265   g_return_val_if_fail (t1 != NULL, FALSE);
 266   g_return_val_if_fail (t2 != NULL, FALSE);
 267   g_return_val_if_fail (e != NULL, FALSE);
 268
 269   if (t1->e1 == e) e1 = t1->e2;
 270   else if (t1->e2 == e) e1 = t1->e3;
 271   else if (t1->e3 == e) e1 = t1->e1;
 272   else
 273     g_assert_not_reached ();
 274   if (t2->e1 == e) e2 = t2->e2;
 275   else if (t2->e2 == e) e2 = t2->e3;
 276   else if (t2->e3 == e) e2 = t2->e1;
 277   else
 278     g_assert_not_reached ();
 279   if (GTS_SEGMENT (e1)->v1 == GTS_SEGMENT (e2)->v1 ||
 280       GTS_SEGMENT (e1)->v1 == GTS_SEGMENT (e2)->v2 ||
 281       GTS_SEGMENT (e1)->v2 == GTS_SEGMENT (e2)->v1 ||
 282       GTS_SEGMENT (e1)->v2 == GTS_SEGMENT (e2)->v2)
 283     return FALSE;
 284   return TRUE;
 285 }
 286
 287 /**
 288  * gts_triangle_area:
 289  * @t: a #GtsTriangle.
 290  *
 291  * Returns: the area of the triangle @t.
 292  */
 293 gdouble gts_triangle_area (GtsTriangle * t)
 294 {
 295   gdouble x, y, z;
 296
 297   g_return_val_if_fail (t != NULL, 0.0);
 298
 299   gts_triangle_normal (t, &x, &y, &z);
 300
 301   return sqrt (x*x + y*y + z*z)/2.;
 302 }
 303
 304 /**
 305  * gts_triangle_perimeter:
 306  * @t: a #GtsTriangle.
 307  *
 308  * Returns: the perimeter of the triangle @t.
 309  */
 310 gdouble gts_triangle_perimeter (GtsTriangle * t)
 311 {
 312   GtsVertex * v;
 313
 314   g_return_val_if_fail (t != NULL, 0.0);
 315
 316   v = gts_triangle_vertex (t);
 317   return
 318     gts_point_distance (GTS_POINT (GTS_SEGMENT (t->e1)->v1),
 319                         GTS_POINT (GTS_SEGMENT (t->e1)->v2)) +
 320     gts_point_distance (GTS_POINT (GTS_SEGMENT (t->e1)->v1),
 321                         GTS_POINT (v)) +
 322     gts_point_distance (GTS_POINT (GTS_SEGMENT (t->e1)->v2),
 323                         GTS_POINT (v));
 324 }
 325
 326 /* perimeter of the equilateral triangle of area unity */
 327 #define GOLDEN_PERIMETER 4.5590141139
 328
 329 /**
 330  * gts_triangle_quality:
 331  * @t: a #GtsTriangle.
 332  *
 333  * The quality of a triangle is defined as the ratio of the square
 334  * root of its surface area to its perimeter relative to this same
 335  * ratio for an equilateral triangle with the same area. The quality
 336  * is then one for an equilateral triangle and tends to zero for a
 337  * very stretched triangle.
 338  *
 339  * Returns: the quality of the triangle @t.
 340  */
 341 gdouble gts_triangle_quality (GtsTriangle * t)
 342 {
 343   gdouble perimeter;
 344
 345   g_return_val_if_fail (t != NULL, 0.0);
 346
 347   perimeter = gts_triangle_perimeter (t);
 348   return perimeter > 0.0 ?
 349     GOLDEN_PERIMETER*sqrt (gts_triangle_area (t))/perimeter :
 350     0.0;
 351 }
 352
 353 /**
 354  * gts_triangle_normal:
 355  * @t: a #GtsTriangle.
 356  * @x: the x coordinate of the normal.
 357  * @y: the y coordinate of the normal.
 358  * @z: the z coordinate of the normal.
 359  *
 360  * Computes the coordinates of the oriented normal of @t as the
 361  * cross-product of two edges, using the left-hand rule. The normal is
 362  * not normalized.  If this triangle is part of a closed and oriented
 363  * surface, the normal points to the outside of the surface.
 364  */
 365 void gts_triangle_normal (GtsTriangle * t,
 366                           gdouble * x,
 367                           gdouble * y,
 368                           gdouble * z)
 369 {
 370   GtsVertex * v1, * v2 = NULL, * v3 = NULL;
 371   GtsPoint * p1, * p2, * p3;
 372   gdouble x1, y1, z1, x2, y2, z2;
 373
 374   g_return_if_fail (t != NULL);
 375
 376   v1 = GTS_SEGMENT (t->e1)->v1;
 377   if (GTS_SEGMENT (t->e1)->v1 == GTS_SEGMENT (t->e2)->v1) {
 378     v2 = GTS_SEGMENT (t->e2)->v2;
 379     v3 = GTS_SEGMENT (t->e1)->v2;
 380   }
 381   else if (GTS_SEGMENT (t->e1)->v2 == GTS_SEGMENT (t->e2)->v2) {
 382     v2 = GTS_SEGMENT (t->e1)->v2;
 383     v3 = GTS_SEGMENT (t->e2)->v1;
 384   }
 385   else if (GTS_SEGMENT (t->e1)->v1 == GTS_SEGMENT (t->e2)->v2) {
 386     v2 = GTS_SEGMENT (t->e2)->v1;
 387     v3 = GTS_SEGMENT (t->e1)->v2;
 388   }
 389   else if (GTS_SEGMENT (t->e1)->v2 == GTS_SEGMENT (t->e2)->v1) {
 390     v2 = GTS_SEGMENT (t->e1)->v2;
 391     v3 = GTS_SEGMENT (t->e2)->v2;
 392   }
 393   else {
 394     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",
 395          t, t->e1, t->e2,
 396          t->e3, GTS_SEGMENT (t->e1)->v1, GTS_SEGMENT (t->e1)->v2,
 397          GTS_SEGMENT (t->e2)->v1, GTS_SEGMENT (t->e2)->v2,
 398          GTS_SEGMENT (t->e3)->v1, GTS_SEGMENT (t->e3)->v2);
 399     g_assert_not_reached ();
 400   }
 401
 402   p1 = GTS_POINT (v1);
 403   p2 = GTS_POINT (v2);
 404   p3 = GTS_POINT (v3);
 405
 406   x1 = p2->x - p1->x;
 407   y1 = p2->y - p1->y;
 408   z1 = p2->z - p1->z;
 409
 410   x2 = p3->x - p1->x;
 411   y2 = p3->y - p1->y;
 412   z2 = p3->z - p1->z;
 413
 414   *x = y1*z2 - z1*y2;
 415   *y = z1*x2 - x1*z2;
 416   *z = x1*y2 - y1*x2;
 417 }
 418
 419 /**
 420  * gts_triangle_orientation:
 421  * @t: a #GtsTriangle.
 422  *
 423  * Checks for the orientation of the plane (x,y) projection of a
 424  * triangle. See gts_point_orientation() for details. This function
 425  * is geometrically robust.
 426  *
 427  * Returns: a number depending on the orientation of the vertices of @t.
 428  */
 429 gdouble gts_triangle_orientation (GtsTriangle * t)
 430 {
 431   GtsVertex * v1, * v2 = NULL, * v3 = NULL;
 432
 433   g_return_val_if_fail (t != NULL, 0.0);
 434
 435   v1 = GTS_SEGMENT (t->e1)->v1;
 436   if (GTS_SEGMENT (t->e1)->v1 == GTS_SEGMENT (t->e2)->v1) {
 437     v2 = GTS_SEGMENT (t->e2)->v2;
 438     v3 = GTS_SEGMENT (t->e1)->v2;
 439   }
 440   else if (GTS_SEGMENT (t->e1)->v2 == GTS_SEGMENT (t->e2)->v2) {
 441     v2 = GTS_SEGMENT (t->e1)->v2;
 442     v3 = GTS_SEGMENT (t->e2)->v1;
 443   }
 444   else if (GTS_SEGMENT (t->e1)->v1 == GTS_SEGMENT (t->e2)->v2) {
 445     v2 = GTS_SEGMENT (t->e2)->v1;
 446     v3 = GTS_SEGMENT (t->e1)->v2;
 447   }
 448   else if (GTS_SEGMENT (t->e1)->v2 == GTS_SEGMENT (t->e2)->v1) {
 449     v2 = GTS_SEGMENT (t->e1)->v2;
 450     v3 = GTS_SEGMENT (t->e2)->v2;
 451   }
 452   else
 453     g_assert_not_reached ();
 454   return gts_point_orientation (GTS_POINT (v1),
 455                                 GTS_POINT (v2),
 456                                 GTS_POINT (v3));
 457 }
 458
 459 /**
 460  * gts_triangle_revert:
 461  * @t: a #GtsTriangle.
 462  *
 463  * Changes the orientation of triangle @t, turning it inside out.
 464  */
 465 void gts_triangle_revert (GtsTriangle * t)
 466 {
 467   GtsEdge * e;
 468
 469   g_return_if_fail (t != NULL);
 470
 471   e = t->e1;
 472   t->e1 = t->e2;
 473   t->e2 = e;
 474 }
 475
 476 /**
 477  * gts_triangles_from_edges:
 478  * @edges: a list of #GtsEdge.
 479  *
 480  * Builds a list of unique triangles which have one of their edges in @edges.
 481  *
 482  * Returns: the list of triangles.
 483  */
 484 GSList * gts_triangles_from_edges (GSList * edges)
 485 {
 486   GHashTable * hash;
 487   GSList * triangles = NULL, * i;
 488
 489   hash = g_hash_table_new (NULL, NULL);
 490   i = edges;
 491   while (i) {
 492     GSList * j = GTS_EDGE (i->data)->triangles;
 493     while (j) {
 494       GtsTriangle * t = j->data;
 495       if (g_hash_table_lookup (hash, t) == NULL) {
 496         triangles = g_slist_prepend (triangles, t);
 497         g_hash_table_insert (hash, t, i);
 498       }
 499       j = j->next;
 500     }
 501     i = i->next;
 502   }
 503   g_hash_table_destroy (hash);
 504
 505   return triangles;
 506 }
 507
 508 /**
 509  * gts_triangle_vertices_edges:
 510  * @t: a #GtsTriangle.
 511  * @e: a #GtsEdge belonging to the edges of @t or %NULL.
 512  * @v1: a #GtsVertex used by @t.
 513  * @v2: a #GtsVertex used by @t.
 514  * @v3: a #GtsVertex used by @t.
 515  * @e1: a #GtsEdge used by @t.
 516  * @e2: a #GtsEdge used by @t.
 517  * @e3: a #GtsEdge used by @t.
 518  *
 519  * Given @t and @e, returns @v1, @v2, @v3, @e1, @e2 and @e3. @e1
 520  * has @v1 and @v2 as vertices, @e2 has @v2 and @v3 as vertices
 521  * and @e3 has @v3 and @v1 as vertices. @v1, @v2 and @v3 respects
 522  * the orientation of @t. If @e is not NULL, @e1 and @e are
 523  * identical.
 524  */
 525 void gts_triangle_vertices_edges (GtsTriangle * t,
 526                                   GtsEdge * e,
 527                                   GtsVertex ** v1,
 528                                   GtsVertex ** v2,
 529                                   GtsVertex ** v3,
 530                                   GtsEdge ** e1,
 531                                   GtsEdge ** e2,
 532                                   GtsEdge ** e3)
 533 {
 534   GtsEdge * ee1, * ee2;
 535
 536   g_return_if_fail (t != NULL);
 537
 538   if (e == t->e1 || e == NULL) {
 539     *e1 = ee1 = t->e1; *e2 = ee2 = t->e2; *e3 = t->e3;
 540   }
 541   else if (e == t->e2) {
 542     *e1 = ee1 = e; *e2 = ee2 = t->e3; *e3 = t->e1;
 543   }
 544   else if (e == t->e3) {
 545     *e1 = ee1 = e; *e2 = ee2 = t->e1; *e3 = t->e2;
 546   }
 547   else {
 548     g_assert_not_reached ();
 549     ee1 = ee2 = NULL; /* to avoid complaints from the compiler */
 550   }
 551   if (GTS_SEGMENT (ee1)->v2 == GTS_SEGMENT (ee2)->v1) {
 552     *v1 = GTS_SEGMENT (ee1)->v1;
 553     *v2 = GTS_SEGMENT (ee1)->v2;
 554     *v3 = GTS_SEGMENT (ee2)->v2;
 555   }
 556   else if (GTS_SEGMENT (ee1)->v2 == GTS_SEGMENT (ee2)->v2) {
 557     *v1 = GTS_SEGMENT (ee1)->v1;
 558     *v2 = GTS_SEGMENT (ee1)->v2;
 559     *v3 = GTS_SEGMENT (ee2)->v1;
 560   }
 561   else if (GTS_SEGMENT (ee1)->v1 == GTS_SEGMENT (ee2)->v1) {
 562     *v1 = GTS_SEGMENT (ee1)->v2;
 563     *v2 = GTS_SEGMENT (ee1)->v1;
 564     *v3 = GTS_SEGMENT (ee2)->v2;
 565   }
 566   else if (GTS_SEGMENT (ee1)->v1 == GTS_SEGMENT (ee2)->v2) {
 567     *v1 = GTS_SEGMENT (ee1)->v2;
 568     *v2 = GTS_SEGMENT (ee1)->v1;
 569     *v3 = GTS_SEGMENT (ee2)->v1;
 570   }
 571   else
 572     g_assert_not_reached ();
 573 }
 574
 575 /* sqrt(3) */
 576 #define SQRT3 1.73205080757
 577
 578 /**
 579  * gts_triangle_enclosing:
 580  * @klass: the class of the new triangle.
 581  * @points: a list of #GtsPoint.
 582  * @scale: a scaling factor (must be larger than one).
 583  *
 584  * Builds a new triangle (including new vertices and edges) enclosing
 585  * the plane projection of all the points in @points. This triangle is
 586  * equilateral and encloses a rectangle defined by the maximum and
 587  * minimum x and y coordinates of the points. @scale is an homothetic
 588  * scaling factor. If equal to one, the triangle encloses exactly the
 589  * enclosing rectangle.
 590  *
 591  * Returns: a new #GtsTriangle.
 592  */
 593 GtsTriangle * gts_triangle_enclosing (GtsTriangleClass * klass,
 594                                       GSList * points, gdouble scale)
 595 {
 596   gdouble xmax, xmin, ymax, ymin;
 597   gdouble xo, yo, r;
 598   GtsVertex * v1, * v2, * v3;
 599   GtsEdge * e1, * e2, * e3;
 600
 601   if (points == NULL)
 602     return NULL;
 603
 604   xmax = xmin = GTS_POINT (points->data)->x;
 605   ymax = ymin = GTS_POINT (points->data)->y;
 606   points = points->next;
 607   while (points) {
 608     GtsPoint * p = points->data;
 609     if (p->x > xmax) xmax = p->x;
 610     else if (p->x < xmin) xmin = p->x;
 611     if (p->y > ymax) ymax = p->y;
 612     else if (p->y < ymin) ymin = p->y;
 613     points = points->next;
 614   }
 615   xo = (xmax + xmin)/2.;
 616   yo = (ymax + ymin)/2.;
 617   r = scale*sqrt((xmax - xo)*(xmax - xo) + (ymax - yo)*(ymax - yo));
 618   if (r == 0.0) r = scale;
 619   v1 = gts_vertex_new (gts_vertex_class (),
 620                        xo + r*SQRT3, yo - r, 0.0);
 621   v2 = gts_vertex_new (gts_vertex_class (),
 622                        xo, yo + 2.*r, 0.0);
 623   v3 = gts_vertex_new (gts_vertex_class (),
 624                        xo - r*SQRT3, yo - r, 0.0);
 625   e1 = gts_edge_new (gts_edge_class (), v1, v2);
 626   e2 = gts_edge_new (gts_edge_class (), v2, v3);
 627   e3 = gts_edge_new (gts_edge_class (), v3, v1);
 628   return gts_triangle_new (gts_triangle_class (), e1, e2, e3);
 629 }
 630
 631 /**
 632  * gts_triangle_neighbor_number:
 633  * @t: a #GtsTriangle.
 634  *
 635  * Returns: the number of triangles neighbors of @t.
 636  */
 637 guint gts_triangle_neighbor_number (GtsTriangle * t)
 638 {
 639   GSList * i;
 640   guint nn = 0;
 641   GtsEdge * ee[4], ** e = ee;
 642
 643   g_return_val_if_fail (t != NULL, 0);
 644
 645   ee[0] = t->e1; ee[1] = t->e2; ee[2] = t->e3; ee[3] = NULL;
 646   while (*e) {
 647     i = (*e++)->triangles;
 648     while (i) {
 649       GtsTriangle * t1 = i->data;
 650       if (t1 != t)
 651         nn++;
 652       i = i->next;
 653     }
 654   }
 655   return nn;
 656 }
 657
 658 /**
 659  * gts_triangle_neighbors:
 660  * @t: a #GtsTriangle.
 661  *
 662  * Returns: a list of #GtsTriangle neighbors of @t.
 663  */
 664 GSList * gts_triangle_neighbors (GtsTriangle * t)
 665 {
 666   GSList * i, * list = NULL;
 667   GtsEdge * ee[4], ** e = ee;
 668
 669   g_return_val_if_fail (t != NULL, NULL);
 670
 671   ee[0] = t->e1; ee[1] = t->e2; ee[2] = t->e3; ee[3] = NULL;
 672   while (*e) {
 673     i = (*e++)->triangles;
 674     while (i) {
 675       GtsTriangle * t1 = i->data;
 676       if (t1 != t)
 677         list = g_slist_prepend (list, t1);
 678       i = i->next;
 679     }
 680   }
 681   return list;
 682 }
 683
 684 /**
 685  * gts_triangles_common_edge:
 686  * @t1: a #GtsTriangle.
 687  * @t2: a #GtsTriangle.
 688  *
 689  * Returns: a #GtsEdge common to both @t1 and @t2 or %NULL if @t1 and @t2
 690  * do not share any edge.
 691  */
 692 GtsEdge * gts_triangles_common_edge (GtsTriangle * t1,
 693                                      GtsTriangle * t2)
 694 {
 695   g_return_val_if_fail (t1 != NULL, NULL);
 696   g_return_val_if_fail (t2 != NULL, NULL);
 697
 698   if (t1->e1 == t2->e1 || t1->e1 == t2->e2 || t1->e1 == t2->e3)
 699     return t1->e1;
 700   if (t1->e2 == t2->e1 || t1->e2 == t2->e2 || t1->e2 == t2->e3)
 701     return t1->e2;
 702   if (t1->e3 == t2->e1 || t1->e3 == t2->e2 || t1->e3 == t2->e3)
 703     return t1->e3;
 704   return NULL;
 705 }
 706
 707 /**
 708  * gts_triangle_is_duplicate:
 709  * @t: a #GtsTriangle.
 710  *
 711  * Returns: a #GtsTriangle different from @t but sharing all its edges
 712  * with @t or %NULL if there is none.
 713  */
 714 GtsTriangle * gts_triangle_is_duplicate (GtsTriangle * t)
 715 {
 716   GSList * i;
 717   GtsEdge * e2, * e3;
 718
 719   g_return_val_if_fail (t != NULL, NULL);
 720
 721   e2 = t->e2;
 722   e3 = t->e3;
 723   i = t->e1->triangles;
 724   while (i) {
 725     GtsTriangle * t1 = i->data;
 726     if (t1 != t &&
 727         (t1->e1 == e2 || t1->e2 == e2 || t1->e3 == e2) &&
 728         (t1->e1 == e3 || t1->e2 == e3 || t1->e3 == e3))
 729       return t1;
 730     i = i->next;
 731   }
 732
 733   return NULL;
 734 }
 735
 736 /**
 737  * gts_triangle_use_edges:
 738  * @e1: a #GtsEdge.
 739  * @e2: a #GtsEdge.
 740  * @e3: a #GtsEdge.
 741  *
 742  * Returns: a #GtsTriangle having @e1, @e2 and @e3 as edges or %NULL if @e1,
 743  * @e2 and @e3 are not part of any triangle.
 744  */
 745 GtsTriangle * gts_triangle_use_edges (GtsEdge * e1,
 746                                       GtsEdge * e2,
 747                                       GtsEdge * e3)
 748 {
 749   GSList * i;
 750
 751   g_return_val_if_fail (e1 != NULL, NULL);
 752   g_return_val_if_fail (e2 != NULL, NULL);
 753   g_return_val_if_fail (e3 != NULL, NULL);
 754
 755   i = e1->triangles;
 756   while (i) {
 757     GtsTriangle * t = i->data;
 758     if ((t->e1 == e2 && (t->e2 == e3 || t->e3 == e3)) ||
 759         (t->e2 == e2 && (t->e1 == e3 || t->e3 == e3)) ||
 760         (t->e3 == e2 && (t->e1 == e3 || t->e2 == e3)))
 761       return t;
 762     i = i->next;
 763   }
 764
 765   return NULL;
 766 }
 767
 768 /**
 769  * gts_triangle_is_ok:
 770  * @t: a #GtsTriangle.
 771  *
 772  * Returns: %TRUE if @t is a non-degenerate, non-duplicate triangle,
 773  * %FALSE otherwise.
 774  */
 775 gboolean gts_triangle_is_ok (GtsTriangle * t)
 776 {
 777   g_return_val_if_fail (t != NULL, FALSE);
 778   g_return_val_if_fail (t->e1 != NULL, FALSE);
 779   g_return_val_if_fail (t->e2 != NULL, FALSE);
 780   g_return_val_if_fail (t->e3 != NULL, FALSE);
 781   g_return_val_if_fail (t->e1 != t->e2 && t->e1 != t->e3 && t->e2 != t->e3,
 782                         FALSE);
 783   g_return_val_if_fail (gts_segments_touch (GTS_SEGMENT (t->e1),
 784                                             GTS_SEGMENT (t->e2)),
 785                         FALSE);
 786   g_return_val_if_fail (gts_segments_touch (GTS_SEGMENT (t->e1),
 787                                             GTS_SEGMENT (t->e3)),
 788                         FALSE);
 789   g_return_val_if_fail (gts_segments_touch (GTS_SEGMENT (t->e2),
 790                                             GTS_SEGMENT (t->e3)),
 791                         FALSE);
 792   g_return_val_if_fail (GTS_SEGMENT (t->e1)->v1 != GTS_SEGMENT (t->e1)->v2,
 793                         FALSE);
 794   g_return_val_if_fail (GTS_SEGMENT (t->e2)->v1 != GTS_SEGMENT (t->e2)->v2,
 795                         FALSE);
 796   g_return_val_if_fail (GTS_SEGMENT (t->e3)->v1 != GTS_SEGMENT (t->e3)->v2,
 797                         FALSE);
 798   g_return_val_if_fail (GTS_OBJECT (t)->reserved == NULL, FALSE);
 799   g_return_val_if_fail (!gts_triangle_is_duplicate (t), FALSE);
 800   return TRUE;
 801 }
 802
 803 /**
 804  * gts_triangle_vertices:
 805  * @t: a #GtsTriangle.
 806  * @v1: a pointer on a #GtsVertex.
 807  * @v2: a pointer on a #GtsVertex.
 808  * @v3: a pointer on a #GtsVertex.
 809  *
 810  * Fills @v1, @v2 and @v3 with the oriented set of vertices, summits of @t.
 811  */
 812 void gts_triangle_vertices (GtsTriangle * t,
 813                             GtsVertex ** v1, GtsVertex ** v2, GtsVertex ** v3)
 814 {
 815   GtsSegment * e1, * e2;
 816
 817   g_return_if_fail (t != NULL);
 818   g_return_if_fail (v1 != NULL && v2 != NULL && v3 != NULL);
 819
 820   e1 = GTS_SEGMENT (t->e1);
 821   e2 = GTS_SEGMENT (t->e2);
 822
 823   if (GTS_SEGMENT (e1)->v2 == GTS_SEGMENT (e2)->v1) {
 824     *v1 = GTS_SEGMENT (e1)->v1;
 825     *v2 = GTS_SEGMENT (e1)->v2;
 826     *v3 = GTS_SEGMENT (e2)->v2;
 827   }
 828   else if (GTS_SEGMENT (e1)->v2 == GTS_SEGMENT (e2)->v2) {
 829     *v1 = GTS_SEGMENT (e1)->v1;
 830     *v2 = GTS_SEGMENT (e1)->v2;
 831     *v3 = GTS_SEGMENT (e2)->v1;
 832   }
 833   else if (GTS_SEGMENT (e1)->v1 == GTS_SEGMENT (e2)->v1) {
 834     *v1 = GTS_SEGMENT (e1)->v2;
 835     *v2 = GTS_SEGMENT (e1)->v1;
 836     *v3 = GTS_SEGMENT (e2)->v2;
 837   }
 838   else {
 839     *v1 = GTS_SEGMENT (e1)->v2;
 840     *v2 = GTS_SEGMENT (e1)->v1;
 841     *v3 = GTS_SEGMENT (e2)->v1;
 842   }
 843 }
 844
 845 /**
 846  * gts_triangle_circumcircle_center:
 847  * @t: a #GtsTriangle.
 848  * @point_class: a #GtsPointClass.
 849  *
 850  * Returns: a new #GtsPoint, center of the circumscribing circle of @t or
 851  * %NULL if the circumscribing circle is not defined.
 852  */
 853 GtsPoint * gts_triangle_circumcircle_center (GtsTriangle * t,
 854                                              GtsPointClass * point_class)
 855 {
 856   GtsVertex * v1, * v2, * v3;
 857   gdouble xa, ya, xb, yb, xc, yc;
 858   gdouble xd, yd, xe, ye;
 859   gdouble xad, yad, xae, yae;
 860   gdouble det;
 861
 862   g_return_val_if_fail (t != NULL, NULL);
 863   g_return_val_if_fail (point_class != NULL, NULL);
 864
 865   gts_triangle_vertices (t, &v1, &v2, &v3);
 866
 867   xa = GTS_POINT (v1)->x; ya = GTS_POINT (v1)->y;
 868   xb = GTS_POINT (v2)->x; yb = GTS_POINT (v2)->y;
 869   xc = GTS_POINT (v3)->x; yc = GTS_POINT (v3)->y;
 870   xd = (xa + xb)/2.; yd = (ya + yb)/2.;
 871   xe = (xa + xc)/2.; ye = (ya + yc)/2.;
 872   xad = xd - xa; yad = yd - ya;
 873   xae = xe - xa; yae = ye - ya;
 874   det = xad*yae - xae*yad;
 875   if (det == 0.)
 876     return NULL;
 877   return gts_point_new (point_class,
 878                         (yae*yad*(yd - ye) + xad*yae*xd - xae*yad*xe)/det,
 879                         -(xae*xad*(xd - xe) + yad*xae*yd - yae*xad*ye)/det,
 880                         0.);
 881 }
 882
 883 /* square of the maximum area ratio admissible */
 884 #define AREA_RATIO_MAX2 1e8
 885
 886 static gboolean points_are_folded (GtsPoint * A,
 887                                    GtsPoint * B,
 888                                    GtsPoint * C,
 889                                    GtsPoint * D,
 890                                    gdouble max)
 891 {
 892   GtsVector AB, AC, AD;
 893   GtsVector n1, n2;
 894   gdouble nn1, nn2, n1n2;
 895
 896   gts_vector_init (AB, A, B);
 897   gts_vector_init (AC, A, C);
 898   gts_vector_init (AD, A, D);
 899   gts_vector_cross (n1, AB, AC);
 900   gts_vector_cross (n2, AD, AB);
 901
 902   nn1 = gts_vector_scalar (n1, n1);
 903   nn2 = gts_vector_scalar (n2, n2);
 904   if (nn1 >= AREA_RATIO_MAX2*nn2 || nn2 >= AREA_RATIO_MAX2*nn1)
 905     return TRUE;
 906   n1n2 = gts_vector_scalar (n1, n2);
 907   if (n1n2 > 0.)
 908     return FALSE;
 909   if (n1n2*n1n2/(nn1*nn2) > max)
 910     return TRUE;
 911   return FALSE;
 912 }
 913
 914 static GtsVertex * triangle_use_vertices (GtsTriangle * t,
 915                                           GtsVertex * A,
 916                                           GtsVertex * B)
 917 {
 918   GtsVertex
 919     * v1 = GTS_SEGMENT (t->e1)->v1,
 920     * v2 = GTS_SEGMENT (t->e1)->v2,
 921     * v3 = gts_triangle_vertex (t);
 922
 923   if (v1 == A) {
 924     if (v2 == B)
 925       return v3;
 926     g_assert (v3 == B);
 927     return v2;
 928   }
 929   if (v2 == A) {
 930     if (v1 == B)
 931       return v3;
 932     g_assert (v3 == B);
 933     return v1;
 934   }
 935   if (v3 == A) {
 936     if (v1 == B)
 937       return v2;
 938     g_assert (v2 == B);
 939     return v1;
 940   }
 941   g_assert_not_reached ();
 942   return NULL;
 943 }
 944
 945 /**
 946  * gts_triangles_are_folded:
 947  * @triangles: a list of #GtsTriangle.
 948  * @A: a #GtsVertex.
 949  * @B: another #GtsVertex.
 950  * @max: the maximum value of the square of the cosine of the angle between
 951  * two triangles.
 952  *
 953  * Given a list of triangles sharing @A and @B as vertices, checks if any
 954  * two triangles in the list make an angle larger than a given value defined
 955  * by @max.
 956  *
 957  * Returns: %TRUE if any pair of triangles in @triangles makes an angle larger
 958  * than the maximum value, %FALSE otherwise.
 959  */
 960 gboolean gts_triangles_are_folded (GSList * triangles,
 961                                    GtsVertex * A, GtsVertex * B,
 962                                    gdouble max)
 963 {
 964   GSList * i;
 965
 966   g_return_val_if_fail (A != NULL, TRUE);
 967   g_return_val_if_fail (B != NULL, TRUE);
 968
 969   i = triangles;
 970   while (i) {
 971     GtsVertex * C = triangle_use_vertices (i->data, A, B);
 972     GSList * j = i->next;
 973     while (j) {
 974       GtsVertex * D = triangle_use_vertices (j->data, A, B);
 975       if (points_are_folded (GTS_POINT (A),
 976                              GTS_POINT (B),
 977                              GTS_POINT (C),
 978                              GTS_POINT (D),
 979                              max))
 980         return TRUE;
 981       j = j->next;
 982     }
 983     i = i->next;
 984   }
 985   return FALSE;
 986 }
 987
 988 /**
 989  * gts_triangle_is_stabbed:
 990  * @t: a #GtsTriangle.
 991  * @p: a #GtsPoint.
 992  * @orientation: a pointer or %NULL.
 993  *
 994  * Returns: one of the vertices of @t, one of the edges of @t or @t if
 995  * any of these are stabbed by the ray starting at @p (included) and
 996  * ending at (@p->x, @p->y, +infty), %NULL otherwise. If the ray is
 997  * contained in the plane of the triangle %NULL is also returned. If
 998  * @orientation is not %NULL, it is set to the value of the
 999  * orientation of @p relative to @t (as given by
1000  * gts_point_orientation_3d()).
1001  */
1002 GtsObject * gts_triangle_is_stabbed (GtsTriangle * t,
1003                                      GtsPoint * p,
1004                                      gdouble * orientation)
1005 {
1006   GtsVertex * v1, * v2, * v3, * inverted = NULL;
1007   GtsEdge * e1, * e2, * e3, * tmp;
1008   gdouble o, o1, o2, o3;
1009
1010   g_return_val_if_fail (t != NULL, NULL);
1011   g_return_val_if_fail (p != NULL, NULL);
1012
1013   gts_triangle_vertices_edges (t, NULL, &v1, &v2, &v3, &e1, &e2, &e3);
1014   o = gts_point_orientation (GTS_POINT (v1), GTS_POINT (v2), GTS_POINT (v3));
1015   if (o == 0.)
1016     return NULL;
1017   if (o < 0.) {
1018     inverted = v1;
1019     v1 = v2;
1020     v2 = inverted;
1021     tmp = e2;
1022     e2 = e3;
1023     e3 = tmp;
1024   }
1025   o = gts_point_orientation_3d (GTS_POINT (v1),
1026                                 GTS_POINT (v2),
1027                                 GTS_POINT (v3),
1028                                 p);
1029   if (o < 0.)
1030     return NULL;
1031   o1 = gts_point_orientation (GTS_POINT (v1), GTS_POINT (v2), p);
1032   if (o1 < 0.)
1033     return NULL;
1034   o2 = gts_point_orientation (GTS_POINT (v2), GTS_POINT (v3), p);
1035   if (o2 < 0.)
1036     return NULL;
1037   o3 = gts_point_orientation (GTS_POINT (v3), GTS_POINT (v1), p);
1038   if (o3 < 0.)
1039     return NULL;
1040   if (orientation) *orientation = inverted ? -o : o;
1041   if (o1 == 0.) {
1042     if (o2 == 0.)
1043       return GTS_OBJECT (v2);
1044     if (o3 == 0.)
1045       return GTS_OBJECT (v1);
1046     return GTS_OBJECT (e1);
1047   }
1048   if (o2 == 0.) {
1049     if (o3 == 0.)
1050       return GTS_OBJECT (v3);
1051     return GTS_OBJECT (e2);
1052   }
1053   if (o3 == 0.)
1054     return GTS_OBJECT (e3);
1055   return GTS_OBJECT (t);
1056 }
1057
1058 /**
1059  * gts_triangle_interpolate_height:
1060  * @t: a #GtsTriangle.
1061  * @p: a #GtsPoint.
1062  *
1063  * Fills the z-coordinate of point @p belonging to the plane
1064  * projection of triangle @t with the linearly interpolated value of
1065  * the z-coordinates of the vertices of @t.
1066  */
1067 void gts_triangle_interpolate_height (GtsTriangle * t, GtsPoint * p)
1068 {
1069   GtsPoint * p1, * p2, * p3;
1070   gdouble x1, x2, y1, y2, det;
1071
1072   g_return_if_fail (t != NULL);
1073   g_return_if_fail (p != NULL);
1074
1075   p1 = GTS_POINT (GTS_SEGMENT (t->e1)->v1);
1076   p2 = GTS_POINT (GTS_SEGMENT (t->e1)->v2);
1077   p3 = GTS_POINT (gts_triangle_vertex (t));
1078
1079   x1 = p2->x - p1->x;
1080   y1 = p2->y - p1->y;
1081   x2 = p3->x - p1->x;
1082   y2 = p3->y - p1->y;
1083   det = x1*y2 - x2*y1;
1084   if (det == 0.)
1085     p->z = (p1->z + p2->z + p3->z)/3.;
1086   else {
1087     gdouble x = p->x - p1->x;
1088     gdouble y = p->y - p1->y;
1089     gdouble a = (x*y2 - y*x2)/det;
1090     gdouble b = (y*x1 - x*y1)/det;
1091
1092     p->z = (1. - a - b)*p1->z + a*p2->z + b*p3->z;
1093   }
1094 }