gts/point.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 <stdlib.h>
  22 #include "gts.h"
  23 #include "gts-private.h"
  24 #include "predicates.h"
  25
  26 static void point_read (GtsObject ** o, GtsFile * f)
  27 {
  28   GtsPoint * p = GTS_POINT (*o);
  29
  30   if (GTS_POINT_CLASS ((*o)->klass)->binary) {
  31     if (gts_file_read (f, &(p->x), sizeof (gdouble), 1) != 1) {
  32       gts_file_error (f, "expecting a binary number (x coordinate)");
  33       return;
  34     }
  35     if (gts_file_read (f, &(p->y), sizeof (gdouble), 1) != 1) {
  36       gts_file_error (f, "expecting a binary number (y coordinate)");
  37       return;
  38     }
  39     if (gts_file_read (f, &(p->z), sizeof (gdouble), 1) != 1) {
  40       gts_file_error (f, "expecting a binary number (z coordinate)");
  41       return;
  42     }
  43   }
  44   else {
  45     if (f->type != GTS_INT && f->type != GTS_FLOAT) {
  46       gts_file_error (f, "expecting a number (x coordinate)");
  47       return;
  48     }
  49     p->x = atof (f->token->str);
  50
  51     gts_file_next_token (f);
  52     if (f->type != GTS_INT && f->type != GTS_FLOAT) {
  53       gts_file_error (f, "expecting a number (y coordinate)");
  54       return;
  55     }
  56     p->y = atof (f->token->str);
  57
  58     gts_file_next_token (f);
  59     if (f->type != GTS_INT && f->type != GTS_FLOAT) {
  60       gts_file_error (f, "expecting a number (z coordinate)");
  61       return;
  62     }
  63     p->z = atof (f->token->str);
  64
  65     gts_file_next_token (f);
  66   }
  67 }
  68
  69 static void point_write (GtsObject * o, FILE * fptr)
  70 {
  71   GtsPoint * p = GTS_POINT (o);
  72
  73   if (GTS_POINT_CLASS ((o)->klass)->binary) {
  74     fwrite (&(p->x), sizeof (gdouble), 1, fptr);
  75     fwrite (&(p->y), sizeof (gdouble), 1, fptr);
  76     fwrite (&(p->z), sizeof (gdouble), 1, fptr);
  77   }
  78   else
  79     fprintf (fptr, "%.10g %.10g %.10g", p->x, p->y, p->z);
  80 }
  81
  82 static void point_class_init (GtsObjectClass * klass)
  83 {
  84   klass->read = point_read;
  85   klass->write = point_write;
  86 }
  87
  88 /**
  89  * gts_point_class:
  90  *
  91  * Returns: the #GtsPointClass.
  92  */
  93 GtsPointClass * gts_point_class (void)
  94 {
  95   static GtsPointClass * klass = NULL;
  96
  97   if (klass == NULL) {
  98     GtsObjectClassInfo point_info = {
  99       "GtsPoint",
 100       sizeof (GtsPoint),
 101       sizeof (GtsPointClass),
 102       (GtsObjectClassInitFunc) point_class_init,
 103       (GtsObjectInitFunc) NULL,
 104       (GtsArgSetFunc) NULL,
 105       (GtsArgGetFunc) NULL
 106     };
 107     klass = gts_object_class_new (gts_object_class (),
 108                                   &point_info);
 109   }
 110
 111   return klass;
 112 }
 113
 114 /**
 115  * gts_point_new:
 116  * @klass: a #GtsPointClass.
 117  * @x: the x-coordinate.
 118  * @y: the y-coordinate.
 119  * @z: the z-coordinate.
 120  *
 121  * Returns: a new #GtsPoint.
 122  */
 123 GtsPoint * gts_point_new (GtsPointClass * klass,
 124                           gdouble x, gdouble y, gdouble z)
 125 {
 126   GtsPoint * p;
 127
 128   p = GTS_POINT (gts_object_new (GTS_OBJECT_CLASS (klass)));
 129   p->x = x;
 130   p->y = y;
 131   p->z = z;
 132
 133   return p;
 134 }
 135
 136 /**
 137  * gts_point_set:
 138  * @p: a #GtsPoint.
 139  * @x: the x-coordinate.
 140  * @y: the y-coordinate.
 141  * @z: the z-coordinate.
 142  *
 143  * Sets the coordinates of @p.
 144  */
 145 void gts_point_set (GtsPoint * p, gdouble x, gdouble y, gdouble z)
 146 {
 147   g_return_if_fail (p != NULL);
 148
 149   p->x = x;
 150   p->y = y;
 151   p->z = z;
 152 }
 153
 154 /**
 155  * gts_point_distance:
 156  * @p1: a #GtsPoint.
 157  * @p2: another #GtsPoint.
 158  *
 159  * Returns: the Euclidean distance between @p1 and @p2.
 160  */
 161 gdouble gts_point_distance (GtsPoint * p1, GtsPoint * p2)
 162 {
 163   g_return_val_if_fail (p1 != NULL && p2 != NULL, 0.0);
 164
 165   return sqrt ((p1->x - p2->x)*(p1->x - p2->x) +
 166                (p1->y - p2->y)*(p1->y - p2->y) +
 167                (p1->z - p2->z)*(p1->z - p2->z));
 168 }
 169
 170 /**
 171  * gts_point_distance2:
 172  * @p1: a #GtsPoint.
 173  * @p2: another #GtsPoint.
 174  *
 175  * Returns: the square of the Euclidean distance between @p1 and @p2.
 176  */
 177 gdouble gts_point_distance2 (GtsPoint * p1, GtsPoint * p2)
 178 {
 179   g_return_val_if_fail (p1 != NULL && p2 != NULL, 0.0);
 180
 181   return
 182     (p1->x - p2->x)*(p1->x - p2->x) +
 183     (p1->y - p2->y)*(p1->y - p2->y) +
 184     (p1->z - p2->z)*(p1->z - p2->z);
 185 }
 186
 187 /**
 188  * gts_point_orientation_3d:
 189  * @p1: a #GtsPoint.
 190  * @p2: a #GtsPoint.
 191  * @p3: a #GtsPoint.
 192  * @p4: a #GtsPoint.
 193  *
 194  * Checks if @p4 lies above, below or on the plane passing through the
 195  * points @p1, @p2 and @p3. Below is defined so that @p1, @p2 and @p3
 196  * appear in counterclockwise order when viewed from above the
 197  * plane. The returned value is an approximation of six times the
 198  * signed volume of the tetrahedron defined by the four points. This
 199  * function uses adaptive floating point arithmetic and is
 200  * consequently geometrically robust.
 201  *
 202  * Returns: a positive value if @p4 lies below, a negative value if
 203  * @p4 lies above the plane, zero if the four points are coplanar.
 204  */
 205 gdouble gts_point_orientation_3d (GtsPoint * p1,
 206                                   GtsPoint * p2,
 207                                   GtsPoint * p3,
 208                                   GtsPoint * p4)
 209 {
 210   g_return_val_if_fail (p1 != NULL && p2 != NULL &&
 211                         p3 != NULL && p4 != NULL, 0.0);
 212   return orient3d ((gdouble *) &p1->x,
 213                    (gdouble *) &p2->x,
 214                    (gdouble *) &p3->x,
 215                    (gdouble *) &p4->x);
 216 }
 217
 218 /**
 219  * gts_point_is_in_triangle:
 220  * @p: a #GtsPoint.
 221  * @t: a #GtsTriangle.
 222  *
 223  * Tests if the planar projection (x, y) of @p is inside, outside or
 224  * on the boundary of the planar projection of @t.  This function is
 225  * geometrically robust.
 226  *
 227  * Returns: %GTS_IN if @p is inside @t, %GTS_ON if @p is on the boundary of
 228  * @t, %GTS_OUT otherwise.
 229  */
 230 GtsIntersect gts_point_is_in_triangle (GtsPoint * p, GtsTriangle * t)
 231 {
 232   GtsVertex * v1, * v2, * v3;
 233   gdouble d1, d2, d3;
 234
 235   g_return_val_if_fail (p != NULL && t != NULL, FALSE);
 236
 237   gts_triangle_vertices (t, &v1, &v2, &v3);
 238
 239   d1 = gts_point_orientation (GTS_POINT (v1), GTS_POINT (v2), p);
 240   if (d1 < 0.0)
 241     return GTS_OUT;
 242   d2 = gts_point_orientation (GTS_POINT (v2), GTS_POINT (v3), p);
 243   if (d2 < 0.0)
 244     return GTS_OUT;
 245   d3 = gts_point_orientation (GTS_POINT (v3), GTS_POINT (v1), p);
 246   if (d3 < 0.0)
 247     return GTS_OUT;
 248   if (d1 == 0.0 || d2 == 0.0 || d3 == 0.0)
 249     return GTS_ON;
 250   return GTS_IN;
 251 }
 252
 253 /**
 254  * gts_point_in_triangle_circle:
 255  * @p: a #GtsPoint.
 256  * @t: a #GtsTriangle.
 257  *
 258  * Tests if the planar projection (x, y) of @p is inside or outside
 259  * the circumcircle of the planar projection of @t. This function is
 260  * geometrically robust.
 261  *
 262  * Returns: a positive number if @p lies inside,
 263  * a negative number if @p lies outside and zero if @p lies on
 264  * the circumcircle of @t.
 265  */
 266 gdouble gts_point_in_triangle_circle (GtsPoint * p, GtsTriangle * t)
 267 {
 268   GtsPoint * p1, * p2, * p3;
 269
 270   g_return_val_if_fail (p != NULL && t != NULL, 0.0);
 271
 272   gts_triangle_vertices (t,
 273                          (GtsVertex **) &p1,
 274                          (GtsVertex **) &p2,
 275                          (GtsVertex **) &p3);
 276
 277   return incircle ((gdouble *) &p1->x,
 278                    (gdouble *) &p2->x,
 279                    (gdouble *) &p3->x,
 280                    (gdouble *) &p->x);
 281 }
 282
 283 /**
 284  * gts_point_in_circle:
 285  * @p: a #GtsPoint.
 286  * @p1: a #GtsPoint.
 287  * @p2: a #GtsPoint.
 288  * @p3: a #GtsPoint.
 289  *
 290  * Tests if the planar projection (x, y) of @p is inside or outside the
 291  * circle defined by the planar projection of @p1, @p2 and @p3.
 292  *
 293  * Returns: a positive number if @p lies inside,
 294  * a negative number if @p lies outside and zero if @p lies on
 295  * the circle.
 296  */
 297 gdouble gts_point_in_circle (GtsPoint * p,
 298                              GtsPoint * p1, GtsPoint * p2, GtsPoint * p3)
 299 {
 300   g_return_val_if_fail (p != NULL && p1 != NULL && p2 != NULL && p3 != NULL,
 301                         0.0);
 302
 303   return incircle ((gdouble *) &p1->x,
 304                    (gdouble *) &p2->x,
 305                    (gdouble *) &p3->x,
 306                    (gdouble *) &p->x);
 307 }
 308
 309 /**
 310  * gts_point_in_sphere:
 311  * @p: a #GtsPoint.
 312  * @p1: a #GtsPoint.
 313  * @p2: a #GtsPoint.
 314  * @p3: a #GtsPoint.
 315  * @p4: a #GtsPoint.
 316  *
 317  * Tests if @p is inside or outside the sphere defined by @p1, @p2,
 318  * @p3 and @p4.
 319  *
 320  * Returns: a positive number if @p lies inside,
 321  * a negative number if @p lies outside and zero if @p lies on
 322  * the sphere.
 323  */
 324 gdouble gts_point_in_sphere (GtsPoint * p,
 325                              GtsPoint * p1, GtsPoint * p2, GtsPoint * p3, GtsPoint * p4)
 326 {
 327   g_return_val_if_fail (p != NULL && p1 != NULL && p2 != NULL && p3 != NULL && p4 != NULL,
 328                         0.0);
 329
 330   return insphere ((gdouble *) &p1->x,
 331                    (gdouble *) &p2->x,
 332                    (gdouble *) &p3->x,
 333                    (gdouble *) &p4->x,
 334                    (gdouble *) &p->x);
 335 }
 336
 337 /**
 338  * gts_point_segment_distance2:
 339  * @p: a #GtsPoint.
 340  * @s: a #GtsSegment.
 341  *
 342  * Returns: the square of the minimun Euclidean distance between @p and @s.
 343  */
 344 gdouble gts_point_segment_distance2 (GtsPoint * p, GtsSegment * s)
 345 {
 346   gdouble t, ns2, x, y, z;
 347   GtsPoint * p1, * p2;
 348
 349   g_return_val_if_fail (p != NULL, 0.0);
 350   g_return_val_if_fail (s != NULL, 0.0);
 351
 352   p1 = GTS_POINT (s->v1);
 353   p2 = GTS_POINT (s->v2);
 354   ns2 = gts_point_distance2 (p1, p2);
 355   if (ns2 == 0.0)
 356     return gts_point_distance2 (p, p1);
 357   t = ((p2->x - p1->x)*(p->x - p1->x) +
 358        (p2->y - p1->y)*(p->y - p1->y) +
 359        (p2->z - p1->z)*(p->z - p1->z))/ns2;
 360   if (t > 1.0)
 361     return gts_point_distance2 (p, p2);
 362   if (t < 0.0)
 363     return gts_point_distance2 (p, p1);
 364   x = (1. - t)*p1->x + t*p2->x - p->x;
 365   y = (1. - t)*p1->y + t*p2->y - p->y;
 366   z = (1. - t)*p1->z + t*p2->z - p->z;
 367   return x*x + y*y + z*z;
 368 }
 369
 370 /**
 371  * gts_point_segment_distance:
 372  * @p: a #GtsPoint.
 373  * @s: a #GtsSegment.
 374  *
 375  * Returns: the minimun Euclidean distance between @p and @s.
 376  */
 377 gdouble gts_point_segment_distance (GtsPoint * p, GtsSegment * s)
 378 {
 379   g_return_val_if_fail (p != NULL, 0.0);
 380   g_return_val_if_fail (s != NULL, 0.0);
 381
 382   return sqrt (gts_point_segment_distance2 (p, s));
 383 }
 384
 385 /**
 386  * gts_point_segment_closest:
 387  * @p: a #GtsPoint.
 388  * @s: a #GtsSegment.
 389  * @closest: a #GtsPoint.
 390  *
 391  * Set the coordinates of @closest to the coordinates of the point belonging
 392  * to @s closest to @p.
 393  */
 394 void gts_point_segment_closest (GtsPoint * p,
 395                                 GtsSegment * s,
 396                                 GtsPoint * closest)
 397 {
 398   gdouble t, ns2;
 399   GtsPoint * p1, * p2;
 400
 401   g_return_if_fail (p != NULL);
 402   g_return_if_fail (s != NULL);
 403   g_return_if_fail (closest != NULL);
 404
 405   p1 = GTS_POINT (s->v1);
 406   p2 = GTS_POINT (s->v2);
 407   ns2 = gts_point_distance2 (p1, p2);
 408
 409   if (ns2 == 0.0) {
 410     gts_point_set (closest, p1->x, p1->y, p1->z);
 411     return;
 412   }
 413
 414   t = ((p2->x - p1->x)*(p->x - p1->x) +
 415        (p2->y - p1->y)*(p->y - p1->y) +
 416        (p2->z - p1->z)*(p->z - p1->z))/ns2;
 417
 418   if (t > 1.0)
 419     gts_point_set (closest, p2->x, p2->y, p2->z);
 420   else if (t < 0.0)
 421     gts_point_set (closest, p1->x, p1->y, p1->z);
 422   else
 423     gts_point_set (closest,
 424                    (1. - t)*p1->x + t*p2->x,
 425                    (1. - t)*p1->y + t*p2->y,
 426                    (1. - t)*p1->z + t*p2->z);
 427 }
 428
 429 /**
 430  * gts_point_triangle_distance2:
 431  * @p: a #GtsPoint.
 432  * @t: a #GtsTriangle.
 433  *
 434  * Returns: the square of the minimun Euclidean distance between @p and @t.
 435  */
 436 gdouble gts_point_triangle_distance2 (GtsPoint * p, GtsTriangle * t)
 437 {
 438   GtsPoint * p1, * p2, * p3;
 439   GtsEdge * e1, * e2, * e3;
 440   GtsVector p1p2, p1p3, pp1;
 441   gdouble A, B, C, D, E, det;
 442   gdouble t1, t2;
 443   gdouble x, y, z;
 444
 445   g_return_val_if_fail (p != NULL, 0.0);
 446   g_return_val_if_fail (t != NULL, 0.0);
 447
 448   gts_triangle_vertices_edges (t, NULL,
 449                                (GtsVertex **) &p1,
 450                                (GtsVertex **) &p2,
 451                                (GtsVertex **) &p3,
 452                                &e1, &e2, &e3);
 453
 454   gts_vector_init (p1p2, p1, p2);
 455   gts_vector_init (p1p3, p1, p3);
 456   gts_vector_init (pp1, p, p1);
 457
 458   B = gts_vector_scalar (p1p3, p1p2);
 459   E = gts_vector_scalar (p1p2, p1p2);
 460   C = gts_vector_scalar (p1p3, p1p3);
 461
 462   det = B*B - E*C;
 463   if (det == 0.) { /* p1p2 and p1p3 are colinear */
 464     gdouble d1 = gts_point_segment_distance2 (p, GTS_SEGMENT (e1));
 465     gdouble d2 = gts_point_segment_distance2 (p, GTS_SEGMENT (e3));
 466     if (d1 < d2)
 467       return d1;
 468     return d2;
 469   }
 470
 471   A = gts_vector_scalar (p1p3, pp1);
 472   D = gts_vector_scalar (p1p2, pp1);
 473
 474   t1 = (D*C - A*B)/det;
 475   t2 = (A*E - D*B)/det;
 476
 477   if (t1 < 0.)
 478     return gts_point_segment_distance2 (p, GTS_SEGMENT (e3));
 479   if (t2 < 0.)
 480     return gts_point_segment_distance2 (p, GTS_SEGMENT (e1));
 481   if (t1 + t2 > 1.)
 482     return gts_point_segment_distance2 (p, GTS_SEGMENT (e2));
 483
 484   x = pp1[0] + t1*p1p2[0] + t2*p1p3[0];
 485   y = pp1[1] + t1*p1p2[1] + t2*p1p3[1];
 486   z = pp1[2] + t1*p1p2[2] + t2*p1p3[2];
 487
 488   return x*x + y*y + z*z;
 489 }
 490
 491 /**
 492  * gts_point_triangle_distance:
 493  * @p: a #GtsPoint.
 494  * @t: a #GtsTriangle.
 495  *
 496  * Returns: the minimun Euclidean distance between @p and @t.
 497  */
 498 gdouble gts_point_triangle_distance (GtsPoint * p, GtsTriangle * t)
 499 {
 500   g_return_val_if_fail (p != NULL, 0.0);
 501   g_return_val_if_fail (t != NULL, 0.0);
 502
 503   return sqrt (gts_point_triangle_distance2 (p, t));
 504 }
 505
 506 /**
 507  * gts_point_triangle_closest:
 508  * @p: a #GtsPoint.
 509  * @t: a #GtsTriangle.
 510  * @closest: a #GtsPoint.
 511  *
 512  * Set the coordinates of @closest to those of the point belonging to @t and
 513  * closest to @p.
 514  */
 515 void gts_point_triangle_closest (GtsPoint * p,
 516                                  GtsTriangle * t,
 517                                  GtsPoint * closest)
 518 {
 519   GtsPoint * p1, * p2, * p3;
 520   GtsEdge * e1, * e2, * e3;
 521   GtsVector p1p2, p1p3, pp1;
 522   gdouble A, B, C, D, E, det;
 523   gdouble t1, t2;
 524
 525   g_return_if_fail (p != NULL);
 526   g_return_if_fail (t != NULL);
 527   g_return_if_fail (closest != NULL);
 528
 529   gts_triangle_vertices_edges (t, NULL,
 530                                (GtsVertex **) &p1,
 531                                (GtsVertex **) &p2,
 532                                (GtsVertex **) &p3,
 533                                &e1, &e2, &e3);
 534
 535   gts_vector_init (p1p2, p1, p2);
 536   gts_vector_init (p1p3, p1, p3);
 537   gts_vector_init (pp1, p, p1);
 538
 539   B = gts_vector_scalar (p1p3, p1p2);
 540   E = gts_vector_scalar (p1p2, p1p2);
 541   C = gts_vector_scalar (p1p3, p1p3);
 542
 543   det = B*B - E*C;
 544   if (det == 0.) { /* p1p2 and p1p3 are colinear */
 545     GtsPoint * cp =
 546       GTS_POINT (gts_object_new (GTS_OBJECT_CLASS (gts_point_class ())));
 547     gts_point_segment_closest (p, GTS_SEGMENT (e1), cp);
 548     gts_point_segment_closest (p, GTS_SEGMENT (e3), closest);
 549
 550     if (gts_point_distance2 (cp, p) < gts_point_distance2 (closest, p))
 551       gts_point_set (closest, cp->x, cp->y, cp->z);
 552     gts_object_destroy (GTS_OBJECT (cp));
 553     return;
 554   }
 555
 556   A = gts_vector_scalar (p1p3, pp1);
 557   D = gts_vector_scalar (p1p2, pp1);
 558
 559   t1 = (D*C - A*B)/det;
 560   t2 = (A*E - D*B)/det;
 561
 562   if (t1 < 0.)
 563     gts_point_segment_closest (p, GTS_SEGMENT (e3), closest);
 564   else if (t2 < 0.)
 565     gts_point_segment_closest (p, GTS_SEGMENT (e1), closest);
 566   else if (t1 + t2 > 1.)
 567     gts_point_segment_closest (p, GTS_SEGMENT (e2), closest);
 568   else
 569     gts_point_set (closest,
 570                    p1->x + t1*p1p2[0] + t2*p1p3[0],
 571                    p1->y + t1*p1p2[1] + t2*p1p3[1],
 572                    p1->z + t1*p1p2[2] + t2*p1p3[2]);
 573 }
 574
 575 /**
 576  * gts_segment_triangle_intersection:
 577  * @s: a #GtsSegment.
 578  * @t: a #GtsTriangle.
 579  * @boundary: if %TRUE, the boundary of @t is taken into account.
 580  * @klass: a #GtsPointClass to be used for the new point.
 581  *
 582  * Checks if @s intersects @t. If this is the case, creates a new
 583  * point pi intersection of @s with @t.
 584  *
 585  * This function is geometrically robust in the sense that it will not
 586  * return a point if @s and @t do not intersect and will return a
 587  * point if @s and @t do intersect. However, the point coordinates are
 588  * subject to round-off errors.
 589  *
 590  * Note that this function will not return any point if @s is contained in
 591  * the plane defined by @t.
 592  *
 593  * Returns: a summit of @t (if @boundary is set to %TRUE), one of the endpoints
 594  * of @s or a new #GtsPoint, intersection of @s with @t or %NULL if @s
 595  * and @t don't intersect.
 596  */
 597 GtsPoint * gts_segment_triangle_intersection (GtsSegment * s,
 598                                               GtsTriangle * t,
 599                                               gboolean boundary,
 600                                               GtsPointClass * klass)
 601 {
 602   GtsPoint * A, * B, * C, * D, * E, * I;
 603   gdouble ABCE, ABCD, ADCE, ABDE, BCDE;
 604   gdouble c;
 605
 606   g_return_val_if_fail (s != NULL, NULL);
 607   g_return_val_if_fail (t != NULL, NULL);
 608   g_return_val_if_fail (klass != NULL, NULL);
 609
 610   A = GTS_POINT (GTS_SEGMENT (t->e1)->v1);
 611   B = GTS_POINT (GTS_SEGMENT (t->e1)->v2);
 612   C = GTS_POINT (gts_triangle_vertex (t));
 613   D = GTS_POINT (s->v1);
 614   E = GTS_POINT (s->v2);
 615
 616   ABCE = gts_point_orientation_3d (A, B, C, E);
 617   ABCD = gts_point_orientation_3d (A, B, C, D);
 618   if (ABCE < 0.0 || ABCD > 0.0) {
 619     GtsPoint * tmpp;
 620     gdouble tmp;
 621     tmpp = E; E = D; D = tmpp;
 622     tmp = ABCE; ABCE = ABCD; ABCD = tmp;
 623   }
 624   if (ABCE < 0.0 || ABCD > 0.0)
 625     return NULL;
 626   ADCE = gts_point_orientation_3d (A, D, C, E);
 627   if ((boundary && ADCE < 0.) || (!boundary && ADCE <= 0.))
 628     return NULL;
 629   ABDE = gts_point_orientation_3d (A, B, D, E);
 630   if ((boundary && ABDE < 0.) || (!boundary && ABDE <= 0.))
 631     return NULL;
 632   BCDE = gts_point_orientation_3d (B, C, D, E);
 633   if ((boundary && BCDE < 0.) || (!boundary && BCDE <= 0.))
 634     return NULL;
 635   if (ABCE == 0.0) {
 636     if (ABCD == 0.0)
 637       /* s is contained in the plane defined by t*/
 638       return NULL;
 639     return E;
 640   }
 641   if (ABCD == 0.0)
 642     return D;
 643   if (boundary) { /* corners of @t */
 644     if (ABDE == 0.) {
 645       if (ADCE == 0.)
 646         return A;
 647       if (BCDE == 0.)
 648         return B;
 649     }
 650     else if (BCDE == 0. && ADCE == 0.)
 651       return C;
 652   }
 653   c = ABCE/(ABCE - ABCD);
 654   I = GTS_POINT (gts_object_new (GTS_OBJECT_CLASS (klass)));
 655   gts_point_set (I,
 656                  E->x + c*(D->x - E->x),
 657                  E->y + c*(D->y - E->y),
 658                  E->z + c*(D->z - E->z));
 659   return I;
 660 }
 661
 662 /**
 663  * gts_point_transform:
 664  * @p: a #GtsPoint.
 665  * @m: the #GtsMatrix representing the transformation to
 666  * apply to the coordinates of @p.
 667  *
 668  * Transform the coordinates of @p according to @m. (p[] becomes m[][].p[]).
 669  */
 670 void gts_point_transform (GtsPoint * p, GtsMatrix * m)
 671 {
 672   gdouble x, y, z;
 673   g_return_if_fail (p != NULL && m != NULL);
 674   x = m[0][0]*p->x + m[0][1]*p->y + m[0][2]*p->z + m[0][3];
 675   y = m[1][0]*p->x + m[1][1]*p->y + m[1][2]*p->z + m[1][3];
 676   z = m[2][0]*p->x + m[2][1]*p->y + m[2][2]*p->z + m[2][3];
 677   p->x = x; p->y = y; p->z = z;
 678 }
 679
 680 /**
 681  * gts_point_orientation:
 682  * @p1: a #GtsPoint.
 683  * @p2: a #GtsPoint.
 684  * @p3: a #GtsPoint.
 685  *
 686  * Checks for orientation of the projection of three points on the
 687  * (x,y) plane. The result is also an approximation of twice the
 688  * signed area of the triangle defined by the three points. This
 689  * function uses adaptive floating point arithmetic and is
 690  * consequently geometrically robust.
 691  *
 692  * Returns: a positive value if @p1, @p2 and @p3 appear in
 693  * counterclockwise order, a negative value if they appear in
 694  * clockwise order and zero if they are colinear.
 695  */
 696 gdouble gts_point_orientation (GtsPoint * p1, GtsPoint * p2, GtsPoint * p3)
 697 {
 698   g_return_val_if_fail (p1 != NULL && p2 != NULL && p3 != NULL, 0.0);
 699
 700   return orient2d ((gdouble *) &p1->x,
 701                    (gdouble *) &p2->x,
 702                    (gdouble *) &p3->x);
 703 }
 704
 705 static gboolean ray_intersects_triangle (GtsPoint * D, GtsPoint * E,
 706                                          GtsTriangle * t)
 707 {
 708   GtsPoint * A, * B, * C;
 709   gint ABCE, ABCD, ADCE, ABDE, BCDE;
 710
 711   gts_triangle_vertices (t, (GtsVertex **) &A,
 712                          (GtsVertex **) &B,
 713                          (GtsVertex **) &C);
 714
 715   ABCE = gts_point_orientation_3d_sos (A, B, C, E);
 716   ABCD = gts_point_orientation_3d_sos (A, B, C, D);
 717   if (ABCE < 0 || ABCD > 0) {
 718     GtsPoint * tmpp;
 719     gint tmp;
 720
 721     tmpp = E; E = D; D = tmpp;
 722     tmp = ABCE; ABCE = ABCD; ABCD = tmp;
 723   }
 724   if (ABCE < 0 || ABCD > 0)
 725     return FALSE;
 726   ADCE = gts_point_orientation_3d_sos (A, D, C, E);
 727   if (ADCE < 0)
 728     return FALSE;
 729   ABDE = gts_point_orientation_3d_sos (A, B, D, E);
 730   if (ABDE < 0)
 731     return FALSE;
 732   BCDE = gts_point_orientation_3d_sos (B, C, D, E);
 733   if (BCDE < 0)
 734     return FALSE;
 735   return TRUE;
 736 }
 737
 738 /**
 739  * gts_point_is_inside_surface:
 740  * @p: a #GtsPoint.
 741  * @tree: a bounding box tree of the faces of a closed, orientable
 742  * surface (see gts_bb_tree_surface()).
 743  * @is_open: %TRUE if the surface defined by @tree is "open" i.e. its volume
 744  * is negative, %FALSE otherwise.
 745  *
 746  * Returns: %TRUE if @p is inside the surface defined by @tree, %FALSE
 747  * otherwise.
 748  */
 749 gboolean gts_point_is_inside_surface (GtsPoint * p,
 750                                       GNode * tree,
 751                                       gboolean is_open)
 752 {
 753   GSList * list, * i;
 754   guint nc = 0;
 755   GtsPoint * p1;
 756   GtsBBox * bb;
 757
 758   g_return_val_if_fail (p != NULL, FALSE);
 759   g_return_val_if_fail (tree != NULL, FALSE);
 760
 761   bb = tree->data;
 762   p1 = gts_point_new (gts_point_class (), bb->x2 + fabs (bb->x2)/10., p->y, p->z);
 763   i = list = gts_bb_tree_stabbed (tree, p);
 764   while (i) {
 765     GtsTriangle * t = GTS_TRIANGLE (GTS_BBOX (i->data)->bounded);
 766
 767     if (ray_intersects_triangle (p, p1, t))
 768       nc++;
 769     i = i->next;
 770   }
 771   g_slist_free (list);
 772   gts_object_destroy (GTS_OBJECT (p1));
 773
 774   return is_open ? (nc % 2 == 0) : (nc % 2 != 0);
 775 }
 776
 777 #define SIGN(x) ((x) > 0. ? 1 : -1)
 778 #define ORIENT1D(a,b) ((a) > (b) ? 1 : (a) < (b) ? -1 : 0)
 779
 780 static gint sortp (gpointer * p, guint n)
 781 {
 782   gint sign = 1;
 783   guint i, j;
 784
 785   for (i = 0; i < n - 1; i++)
 786     for (j = 0; j < n - 1 - i; j++)
 787       if (GPOINTER_TO_UINT (p[j+1]) < GPOINTER_TO_UINT (p[j])) {
 788         gpointer tmp = p[j];
 789
 790         p[j] = p[j+1];
 791         p[j+1] = tmp;
 792         sign = - sign;
 793       }
 794   return sign;
 795 }
 796
 797 /**
 798  * gts_point_orientation_3d_sos:
 799  * @p1: a #GtsPoint.
 800  * @p2: a #GtsPoint.
 801  * @p3: a #GtsPoint.
 802  * @p4: a #GtsPoint.
 803  *
 804  * Checks if @p4 lies above or below the plane passing through the
 805  * points @p1, @p2 and @p3. Below is defined so that @p1, @p2 and @p3
 806  * appear in counterclockwise order when viewed from above the
 807  * plane. This function uses adaptive floating point arithmetic and is
 808  * consequently geometrically robust.
 809  *
 810  * Simulation of Simplicity (SoS) is used to break ties when the
 811  * orientation is degenerate (i.e. @p4 lies on the plane defined by
 812  * @p1, @p2 and @p3).
 813  *
 814  * Returns: +1 if @p4 lies below, -1 if @p4 lies above the plane.
 815  */
 816 gint gts_point_orientation_3d_sos (GtsPoint * p1,
 817                                    GtsPoint * p2,
 818                                    GtsPoint * p3,
 819                                    GtsPoint * p4)
 820 {
 821   gdouble o;
 822
 823   g_return_val_if_fail (p1 != NULL && p2 != NULL &&
 824                         p3 != NULL && p4 != NULL, 0);
 825
 826   o = orient3d ((gdouble *) &p1->x,
 827                 (gdouble *) &p2->x,
 828                 (gdouble *) &p3->x,
 829                 (gdouble *) &p4->x);
 830   if (o != 0.)
 831     return SIGN (o);
 832   else {
 833     GtsPoint * p[4];
 834     gdouble a[2], b[2], c[2];
 835     gint sign;
 836
 837     p[0] = p1; p[1] = p2; p[2] = p3; p[3] = p4;
 838     sign = sortp ((gpointer *) p, 4);
 839
 840     /* epsilon^1/8 */
 841     a[0] = p[1]->x; a[1] = p[1]->y;
 842     b[0] = p[2]->x; b[1] = p[2]->y;
 843     c[0] = p[3]->x; c[1] = p[3]->y;
 844     o = orient2d (a, b, c);
 845     if (o != 0.)
 846       return SIGN (o)*sign;
 847
 848     /* epsilon^1/4 */
 849     a[0] = p[1]->x; a[1] = p[1]->z;
 850     b[0] = p[2]->x; b[1] = p[2]->z;
 851     c[0] = p[3]->x; c[1] = p[3]->z;
 852     o = orient2d (a, b, c);
 853     if (o != 0.)
 854       return - SIGN (o)*sign;
 855
 856     /* epsilon^1/2 */
 857     a[0] = p[1]->y; a[1] = p[1]->z;
 858     b[0] = p[2]->y; b[1] = p[2]->z;
 859     c[0] = p[3]->y; c[1] = p[3]->z;
 860     o = orient2d (a, b, c);
 861     if (o != 0.)
 862       return SIGN (o)*sign;
 863
 864     /* epsilon */
 865     a[0] = p[0]->x; a[1] = p[0]->y;
 866     b[0] = p[2]->x; b[1] = p[2]->y;
 867     c[0] = p[3]->x; c[1] = p[3]->y;
 868     o = orient2d (a, b, c);
 869     if (o != 0.)
 870       return - SIGN (o)*sign;
 871
 872     /* epsilon^5/4 */
 873     o = ORIENT1D (p[2]->x, p[3]->x);
 874     if (o != 0.)
 875       return SIGN (o)*sign;
 876
 877     /* epsilon^3/2 */
 878     o = ORIENT1D (p[2]->y, p[3]->y);
 879     if (o != 0.)
 880       return - SIGN (o)*sign;
 881
 882     /* epsilon^2 */
 883     a[0] = p[0]->x; a[1] = p[0]->z;
 884     b[0] = p[2]->x; b[1] = p[2]->z;
 885     c[0] = p[3]->x; c[1] = p[3]->z;
 886     o = orient2d (a, b, c);
 887     if (o != 0.)
 888       return SIGN (o)*sign;
 889
 890     /* epsilon^5/2 */
 891     o = ORIENT1D (p[2]->z, p[3]->z);
 892     if (o != 0.)
 893       return SIGN (o)*sign;
 894
 895     /* epsilon^4 */
 896     a[0] = p[0]->y; a[1] = p[0]->z;
 897     b[0] = p[2]->y; b[1] = p[2]->z;
 898     c[0] = p[3]->y; c[1] = p[3]->z;
 899     o = orient2d (a, b, c);
 900     if (o != 0.)
 901       return - SIGN (o)*sign;
 902
 903     /* epsilon^8 */
 904     a[0] = p[0]->x; a[1] = p[0]->y;
 905     b[0] = p[1]->x; b[1] = p[1]->y;
 906     c[0] = p[3]->x; c[1] = p[3]->y;
 907     o = orient2d (a, b, c);
 908     if (o != 0.)
 909       return SIGN (o)*sign;
 910
 911     /* epsilon^33/4 */
 912     o = ORIENT1D (p[1]->x, p[3]->x);
 913     if (o != 0.)
 914       return - SIGN (o)*sign;
 915
 916     /* epsilon^17/2 */
 917     o = ORIENT1D (p[1]->y, p[3]->y);
 918     if (o != 0.)
 919       return SIGN (o)*sign;
 920
 921     /* epsilon^10 */
 922     o = ORIENT1D (p[0]->x, p[3]->x);
 923     if (o != 0.)
 924       return SIGN (o)*sign;
 925
 926     /* epsilon^21/2 */
 927     return sign;
 928   }
 929 }
 930
 931 /**
 932  * gts_point_orientation_sos:
 933  * @p1: a #GtsPoint.
 934  * @p2: a #GtsPoint.
 935  * @p3: a #GtsPoint.
 936  *
 937  * Checks for orientation of the projection of three points on the
 938  * (x,y) plane.
 939  *
 940  * Simulation of Simplicity (SoS) is used to break ties when the
 941  * orientation is degenerate (i.e. @p3 lies on the line defined by
 942  * @p1 and @p2).
 943  *
 944  * Returns: a positive value if @p1, @p2 and @p3 appear in
 945  * counterclockwise order or a negative value if they appear in
 946  * clockwise order.
 947  */
 948 gint gts_point_orientation_sos (GtsPoint * p1,
 949                                 GtsPoint * p2,
 950                                 GtsPoint * p3)
 951 {
 952   gdouble o;
 953
 954   g_return_val_if_fail (p1 != NULL && p2 != NULL && p3 != NULL, 0);
 955
 956   o = orient2d ((gdouble *) &p1->x,
 957                 (gdouble *) &p2->x,
 958                 (gdouble *) &p3->x);
 959   if (o != 0.)
 960     return SIGN (o);
 961   else {
 962     GtsPoint * p[3];
 963     gint sign;
 964
 965     p[0] = p1; p[1] = p2; p[2] = p3;
 966     sign = sortp ((gpointer *) p, 3);
 967
 968     /* epsilon^1/4 */
 969     o = ORIENT1D (p[1]->x, p[2]->x);
 970     if (o != 0.)
 971       return - SIGN (o)*sign;
 972
 973     /* epsilon^1/2 */
 974     o = ORIENT1D (p[1]->y, p[2]->y);
 975     if (o != 0.)
 976       return SIGN (o)*sign;
 977
 978     /* epsilon */
 979     o = ORIENT1D (p[0]->x, p[2]->x);
 980     if (o != 0.)
 981       return SIGN (o)*sign;
 982
 983     /* epsilon^3/2 */
 984     return sign;
 985   }
 986 }