pcb 4.1.1
An interactive printed circuit board layout editor.

cdt.c

Go to the documentation of this file.
00001 /* GTS - Library for the manipulation of triangulated surfaces
00002  * Copyright (C) 1999 Stéphane Popinet
00003  *
00004  * This library is free software; you can redistribute it and/or
00005  * modify it under the terms of the GNU Library General Public
00006  * License as published by the Free Software Foundation; either
00007  * version 2 of the License, or (at your option) any later version.
00008  *
00009  * This library is distributed in the hope that it will be useful,
00010  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00011  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00012  * Library General Public License for more details.
00013  *
00014  * You should have received a copy of the GNU Library General Public
00015  * License along with this library; if not, write to the
00016  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
00017  * Boston, MA 02111-1307, USA.
00018  */
00019 
00020 #ifdef HAVE_CONFIG_H
00021 # include <config.h>
00022 #endif
00023 
00024 
00025 #include <math.h>
00026 #include "gts.h"
00027 
00028 #ifdef USE_SURFACE_BTREE
00029 
00030 static gint find_closest (GtsTriangle * t, gpointer value, gpointer * data)
00031 {
00032   guint * ns = data[2];
00033   guint * n = data[3];
00034 
00035   if (*n >= *ns)
00036     return TRUE;
00037   else {
00038     gdouble * dmin = data[0];
00039     gpointer * closest = data[1];
00040     GtsPoint * p = data[4];
00041 
00042     if (gts_triangle_orientation (t) > 0.) {
00043       GtsPoint * p1 = GTS_POINT (GTS_SEGMENT (t->e1)->v1);
00044       gdouble d = (p->x - p1->x)*(p->x - p1->x) + (p->y - p1->y)*(p->y - p1->y);
00045       
00046       if (d < *dmin) {
00047         *dmin = d;
00048         *closest = t;
00049       }
00050       (*n)++;
00051     }
00052   }
00053   return FALSE;
00054 }
00055 
00056 /* select the face closest to @p among n^1/3 randomly picked faces
00057  *  of @surface */
00058 static GtsFace * closest_face (GtsSurface * s, GtsPoint * p)
00059 {
00060   guint n = 0, nt, ns;
00061   gdouble dmin = G_MAXDOUBLE;
00062   GtsFace * closest = NULL;
00063   gpointer data[5];
00064 
00065   nt = gts_surface_face_number (s);
00066   if (!nt)
00067     return NULL;
00068   ns = exp (log ((gdouble) nt)/3.);
00069 
00070   data[0] = &dmin;
00071   data[1] = &closest;
00072   data[2] = &ns;
00073   data[3] = &n;
00074   data[4] = p;
00075   g_tree_traverse (s->faces, (GTraverseFunc) find_closest, G_IN_ORDER, data);
00076 
00077   return closest;
00078 }
00079 
00080 #else /* not USE_SURFACE_BTREE */
00081 
00082 typedef struct _SFindClosest SFindClosest; 
00083 
00084 struct _SFindClosest {
00085   gdouble dmin; 
00086   GtsFace *closest;
00087   GtsPoint * p;
00088   gint stop;
00089 };
00090 
00091 #  if GLIB_CHECK_VERSION(2,4,0)
00092 /* finally, with g_hash_table_find we are able to stop iteration over the hash
00093    table in the middle */
00094 
00095 static gboolean find_closest (gpointer key, gpointer value, gpointer user_data)
00096 {
00097   SFindClosest * data = (SFindClosest *) user_data;
00098   GtsFace * f = GTS_FACE (value);
00099   
00100   if (gts_triangle_orientation (GTS_TRIANGLE (f)) > 0.) {
00101     GtsPoint * p1 = GTS_POINT (GTS_SEGMENT (GTS_TRIANGLE (f)->e1)->v1);
00102     gdouble d = ((data->p->x - p1->x)*(data->p->x - p1->x) + 
00103                  (data->p->y - p1->y)*(data->p->y - p1->y));
00104 
00105     if (d < data->dmin) {
00106       data->dmin = d;
00107       data->closest = f;
00108     }
00109   }
00110   data->stop--;
00111   return !(data->stop > 0);
00112 }
00113 
00114 static GtsFace * closest_face (GtsSurface * s, GtsPoint * p)
00115 {
00116   SFindClosest fc;
00117 
00118   fc.dmin = G_MAXDOUBLE;
00119   fc.closest = NULL;
00120   fc.p = p;
00121   fc.stop = (gint) exp (log ((gdouble) g_hash_table_size (s->faces))/3.);
00122   g_hash_table_find (s->faces, find_closest, &fc);
00123   
00124   return fc.closest;
00125 }
00126 
00127 #  else /* VERSION < 2.4.0 */
00128 
00129 static void
00130 find_closest (gpointer key, gpointer value, gpointer user_data)
00131 {
00132   SFindClosest * data = (SFindClosest *) user_data;
00133   GtsFace * f = GTS_FACE (value);
00134 
00135   if (gts_triangle_orientation (GTS_TRIANGLE (f)) > 0.) {
00136     GtsPoint * p1 = GTS_POINT (GTS_SEGMENT (GTS_TRIANGLE (f)->e1)->v1);
00137     gdouble d = ((data->p->x - p1->x)*(data->p->x - p1->x) +
00138                  (data->p->y - p1->y)*(data->p->y - p1->y));
00139 
00140     if (d < data->dmin) {
00141       data->dmin = d;
00142       data->closest = f;
00143     }
00144   }
00145   data->stop--;
00146 }
00147 
00148 /* select the face closest to @p among n^1/3 randomly picked faces
00149  * of @surface */
00150 static GtsFace * closest_face (GtsSurface * s, GtsPoint * p)
00151 {
00152   SFindClosest fc;
00153 
00154   if (!g_hash_table_size (s->faces))
00155     return NULL;
00156 
00157   fc.dmin = G_MAXDOUBLE;
00158   fc.closest = NULL;
00159   fc.p = p;
00160   fc.stop = (gint) exp (log ((gdouble) g_hash_table_size (s->faces))/3.);
00161   g_hash_table_foreach (s->faces, find_closest, &fc);
00162   return fc.closest;
00163 }
00164 #  endif /* VERSION < 2.4.0 */
00165 #endif /* not USE_SURFACE_BTREE */
00166 
00167 /* returns the face belonging to @surface and neighbor of @f via @e */
00168 static GtsFace * neighbor (GtsFace * f,
00169                            GtsEdge * e,
00170                            GtsSurface * surface)
00171 {
00172   GSList * i = e->triangles;
00173   GtsTriangle * t = GTS_TRIANGLE (f);
00174 
00175   while (i) {
00176     GtsTriangle * t1 = i->data;
00177     if (t1 != t &&
00178         GTS_IS_FACE (t1) &&
00179         gts_face_has_parent_surface (GTS_FACE (t1), surface))
00180       return GTS_FACE (t1);
00181     i = i->next;
00182   }
00183   return NULL;
00184 }
00185 
00186 /* given a triangle @t and a segment s (@o -> @p). 
00187    @o must be in @t. Returns the
00188    edge of @t which is intersected by s or %NULL if @p is also
00189    contained in @t (on_summit is set to %FALSE) or if s intersects @t 
00190    exactly on one of its summit (on_summit is set to %TRUE). */
00191 static GtsEdge * triangle_next_edge (GtsTriangle * t,
00192                                      GtsPoint * o, GtsPoint * p,
00193                                      gboolean * on_summit)
00194 {
00195   GtsVertex * v1, * v2, * v3;
00196   GtsEdge * e1, * e2, * e3;
00197   gdouble orient = 0.0;
00198   
00199   gts_triangle_vertices_edges (t, NULL,
00200                                &v1, &v2, &v3, 
00201                                &e1, &e2, &e3);
00202 
00203   *on_summit = FALSE;
00204   orient = gts_point_orientation (o, GTS_POINT (v1), p);
00205   if (orient > 0.0) {
00206     orient = gts_point_orientation (o, GTS_POINT (v2), p);
00207     if (orient > 0.0) {
00208       if (gts_point_orientation (GTS_POINT (v2), GTS_POINT (v3), p) >= 0.0)
00209         return NULL;
00210       return e2;
00211     }
00212     if (orient < 0.0) {
00213       if (gts_point_orientation (GTS_POINT (v1), GTS_POINT (v2), p) >= 0.0)
00214         return NULL;
00215       return e1;
00216     }
00217     if (gts_point_orientation (GTS_POINT (v1), GTS_POINT (v2), p) < 0.0)
00218       *on_summit = TRUE;
00219     return NULL;
00220   }
00221 
00222   if (orient < 0.0) {
00223     orient = gts_point_orientation (o, GTS_POINT (v3), p);
00224     if (orient > 0.0) {
00225       if (gts_point_orientation (GTS_POINT (v3), GTS_POINT (v1), p) >= 0.0)
00226         return NULL;
00227       return e3;
00228     }
00229     if (orient < 0.0) {
00230       if (gts_point_orientation (GTS_POINT (v2), GTS_POINT (v3), p) >= 0.0)
00231         return NULL;
00232       return e2;
00233     }
00234     if (gts_point_orientation (GTS_POINT (v3), GTS_POINT (v1), p) < 0.0)
00235       *on_summit = TRUE;
00236     return NULL;
00237   }
00238   
00239   if (gts_point_orientation (GTS_POINT (v2), GTS_POINT (v3), p) < 0.0)
00240     return e2;
00241   if (gts_point_orientation (GTS_POINT (v1), GTS_POINT (v2), p) < 0.0)
00242     *on_summit = TRUE;
00243   return NULL;
00244 }
00245 
00246 static void triangle_barycenter (GtsTriangle * t, GtsPoint * b)
00247 {
00248   GtsPoint * p = GTS_POINT (gts_triangle_vertex (t));
00249   b->x = (p->x + 
00250           GTS_POINT (GTS_SEGMENT(t->e1)->v1)->x +
00251           GTS_POINT (GTS_SEGMENT(t->e1)->v2)->x)/3.;
00252   b->y = (p->y + 
00253           GTS_POINT (GTS_SEGMENT(t->e1)->v1)->y +
00254           GTS_POINT (GTS_SEGMENT(t->e1)->v2)->y)/3.;
00255 }
00256 
00257 static GtsFace * point_locate (GtsPoint * o,
00258                                GtsPoint * p,
00259                                GtsFace * f,
00260                                GtsSurface * surface)
00261 {
00262   GtsEdge * prev;
00263   gboolean on_summit;
00264   GtsVertex * v1, * v2, * v3;
00265   GtsEdge * e2, * e3;    
00266   
00267   prev = triangle_next_edge (GTS_TRIANGLE (f), o, p, &on_summit);
00268 
00269   if (!prev) {
00270     GtsFace * f1;
00271 
00272     if (!on_summit)
00273       return f; /* p is inside f */
00274 
00275     /* s intersects f exactly on a summit: restarts from a neighbor of f */
00276     if ((f1 = neighbor (f, GTS_TRIANGLE (f)->e1, surface)) ||
00277         (f1 = neighbor (f, GTS_TRIANGLE (f)->e2, surface)) ||
00278         (f1 = neighbor (f, GTS_TRIANGLE (f)->e3, surface))) {
00279       triangle_barycenter (GTS_TRIANGLE (f1), o);
00280       return point_locate (o, p, f1, surface);
00281     }
00282     return NULL;
00283   }
00284   
00285   f = neighbor (f, prev, surface);
00286   if (f)
00287     gts_triangle_vertices_edges (GTS_TRIANGLE (f), prev, 
00288                                  &v1, &v2, &v3, &prev, &e2, &e3);
00289   while (f) {
00290     gdouble orient = gts_point_orientation (o, GTS_POINT (v3), p);
00291 
00292     if (orient < 0.0) {
00293       if (gts_point_orientation (GTS_POINT (v2), GTS_POINT (v3), p) >= 0.0)
00294         return f; /* p is inside f */
00295       f = neighbor (f, e2, surface);
00296       prev = e2;
00297       v1 = v3;      
00298     }
00299     else if (orient > 0.0) {
00300       if (gts_point_orientation (GTS_POINT (v3), GTS_POINT (v1), p) >= 0.0)
00301         return f; /* p is inside f */
00302       f = neighbor (f, e3, surface);
00303       prev = e3;
00304       v2 = v3;
00305     }
00306     else {
00307       GtsFace * f1;
00308 
00309       if (gts_point_orientation (GTS_POINT (v2), GTS_POINT (v3), p) >= 0.0)
00310         return f; /* p is inside f */
00311 
00312       /* s intersects f exactly on v3: restarts from a neighbor of f */
00313       if ((f1 = neighbor (f, e2, surface)) ||
00314           (f1 = neighbor (f, e3, surface))) {
00315         triangle_barycenter (GTS_TRIANGLE (f1), o);
00316         return point_locate (o, p, f1, surface);
00317       }
00318       return NULL;
00319     }
00320     /* update e2, e3, v3 for the new triangle */
00321     if (f) {
00322       if (prev == GTS_TRIANGLE (f)->e1) {
00323         e2 = GTS_TRIANGLE (f)->e2; e3 = GTS_TRIANGLE (f)->e3;
00324       }
00325       else if (prev == GTS_TRIANGLE (f)->e2) {
00326         e2 = GTS_TRIANGLE (f)->e3; e3 = GTS_TRIANGLE (f)->e1;
00327       }
00328       else {
00329         e2 = GTS_TRIANGLE (f)->e1; e3 = GTS_TRIANGLE (f)->e2;
00330       }
00331       if (GTS_SEGMENT (e2)->v1 == v1 || GTS_SEGMENT (e2)->v1 == v2)
00332         v3 = GTS_SEGMENT (e2)->v2;
00333       else
00334         v3 = GTS_SEGMENT (e2)->v1;
00335     }
00336   }
00337   return NULL;
00338 }
00339 
00357 GtsFace * gts_point_locate (GtsPoint * p, 
00358                             GtsSurface * surface,
00359                             GtsFace * guess)
00360 {
00361   GtsFace * fr;
00362   GtsPoint * o;
00363 
00364   g_return_val_if_fail (p != NULL, NULL);
00365   g_return_val_if_fail (surface != NULL, NULL);
00366   g_return_val_if_fail (guess == NULL || 
00367                         gts_face_has_parent_surface (guess, surface), NULL);
00368 
00369   if (guess == NULL)
00370     guess = closest_face (surface, p);
00371   else
00372     g_return_val_if_fail (gts_triangle_orientation (GTS_TRIANGLE (guess)) > 0., NULL);
00373 
00374   if (guess == NULL)
00375     return NULL;
00376 
00377   o = GTS_POINT (gts_object_new (GTS_OBJECT_CLASS (gts_point_class ())));
00378   triangle_barycenter (GTS_TRIANGLE (guess), o);
00379   fr = point_locate (o, p, guess, surface);
00380   gts_object_destroy (GTS_OBJECT (o));
00381 
00382   return fr;
00383 }
00384 
00385 
00391 GtsConstraintClass * gts_constraint_class (void)
00392 {
00393   static GtsConstraintClass * klass = NULL;
00394 
00395   if (klass == NULL) {
00396     GtsObjectClassInfo constraint_info = {
00397       "GtsConstraint",
00398       sizeof (GtsConstraint),
00399       sizeof (GtsConstraintClass),
00400       (GtsObjectClassInitFunc) NULL,
00401       (GtsObjectInitFunc) NULL,
00402       (GtsArgSetFunc) NULL,
00403       (GtsArgGetFunc) NULL
00404     };
00405     klass = gts_object_class_new (GTS_OBJECT_CLASS (gts_edge_class ()), 
00406                                   &constraint_info);
00407   }
00408 
00409   return klass;
00410 }
00411 
00412 static void split_list (GtsListFace * f, GtsListFace * f1, GtsListFace * f2, 
00413                         GtsPoint * p1, GtsPoint * p2,
00414                         GSList ** last1, GSList ** last2)
00415 {
00416   GSList * i = f->points, * l1 = *last1, * l2 = *last2;
00417 
00418   while (i) {
00419     GtsPoint * p = i->data;
00420     
00421     if (gts_point_orientation (p1, p2, p) >= 0.) {
00422       if (l1) l1->next = i; else f1->points = i;
00423       l1 = i;
00424     }
00425     else {
00426       if (l2) l2->next = i; else f2->points = i;
00427       l2 = i;
00428     }
00429     i = i->next;
00430   }
00431   f->points = NULL;
00432   *last1 = l1;
00433   *last2 = l2;
00434 }
00435 
00436 /* cf. figure misc/swap.fig */
00437 static void swap_if_in_circle (GtsFace * f1,
00438                                GtsVertex * v1, 
00439                                GtsVertex * v2, 
00440                                GtsVertex * v3,
00441                                GtsEdge * e1, 
00442                                GtsEdge * e2, 
00443                                GtsEdge * e3,
00444                                GtsSurface * surface)
00445 {
00446   GtsFace * f2;
00447   GtsEdge * e4, *e5;
00448   GtsVertex * v4;
00449 
00450   if (GTS_IS_CONSTRAINT (e1)) /* @e1 is a constraint can not swap */
00451     return;
00452 
00453   f2 = neighbor (f1, e1, surface);
00454   if (f2 == NULL) /* @e1 is a boundary of @surface */
00455     return;
00456 
00457   if (GTS_TRIANGLE (f2)->e1 == e1) {
00458     e4 = GTS_TRIANGLE (f2)->e2; e5 = GTS_TRIANGLE (f2)->e3;
00459   }
00460   else if (GTS_TRIANGLE (f2)->e2 == e1) {
00461     e4 = GTS_TRIANGLE (f2)->e3; e5 = GTS_TRIANGLE (f2)->e1;
00462   }
00463   else {
00464     e4 = GTS_TRIANGLE (f2)->e1; e5 = GTS_TRIANGLE (f2)->e2;
00465   }
00466   if (GTS_SEGMENT (e4)->v1 == GTS_SEGMENT (e1)->v1 || 
00467       GTS_SEGMENT (e4)->v1 == GTS_SEGMENT (e1)->v2)
00468     v4 = GTS_SEGMENT (e4)->v2;
00469   else
00470     v4 = GTS_SEGMENT (e4)->v1;
00471 
00472   if (gts_point_in_circle (GTS_POINT (v4), GTS_POINT (v1), 
00473                            GTS_POINT (v2), GTS_POINT (v3)) > 0.0) {
00474     GtsEdge * en;
00475     GtsSegment * sn = gts_vertices_are_connected (v3, v4);
00476     GtsFace * f3, * f4;
00477 
00478     if (!GTS_IS_EDGE (sn))
00479       en = gts_edge_new (surface->edge_class, v3, v4);
00480     else
00481       en = GTS_EDGE (sn);
00482 
00483     f3 = gts_face_new (surface->face_class, en, e5, e2);
00484     gts_object_attributes (GTS_OBJECT (f3), GTS_OBJECT (f1));
00485     f4 = gts_face_new (surface->face_class, en, e3, e4);
00486     gts_object_attributes (GTS_OBJECT (f4), GTS_OBJECT (f2));
00487     
00488     if (GTS_IS_LIST_FACE (f3)) {
00489       GSList * last3 = NULL, * last4 = NULL;
00490 
00491       if (GTS_IS_LIST_FACE (f1))
00492         split_list (GTS_LIST_FACE (f1), GTS_LIST_FACE (f3), GTS_LIST_FACE (f4),
00493                     GTS_POINT (v3), GTS_POINT (v4), &last3, &last4);
00494       if (GTS_IS_LIST_FACE (f2))
00495         split_list (GTS_LIST_FACE (f2), GTS_LIST_FACE (f3), GTS_LIST_FACE (f4),
00496                     GTS_POINT (v3), GTS_POINT (v4), &last3, &last4);
00497       if (last3) last3->next = NULL;
00498       if (last4) last4->next = NULL;
00499     }
00500 
00501     gts_surface_remove_face (surface, f1);
00502     gts_surface_remove_face (surface, f2);
00503     gts_surface_add_face (surface, f3);
00504     gts_surface_add_face (surface, f4);
00505 
00506     swap_if_in_circle (f3, v4, v2, v3, e5, e2, en, surface);
00507     swap_if_in_circle (f4, v1, v4, v3, e4, en, e3, surface);
00508   }
00509 }
00510 
00524 GtsVertex * gts_delaunay_add_vertex_to_face (GtsSurface * surface, 
00525                                              GtsVertex * v,
00526                                              GtsFace * f)
00527 {
00528   GtsEdge * e1, * e2, * e3;
00529   GtsSegment * s4, * s5, * s6;
00530   GtsEdge * e4, * e5, * e6;
00531   GtsVertex * v1, * v2, * v3;
00532   GtsFace * nf[3];
00533 
00534   g_return_val_if_fail (surface != NULL, v);
00535   g_return_val_if_fail (v != NULL, v);
00536   g_return_val_if_fail (f != NULL, v);
00537 
00538   gts_triangle_vertices_edges (GTS_TRIANGLE (f), NULL, 
00539                                &v1, &v2, &v3, &e1, &e2, &e3);
00540   if (v == v1 || v == v2 || v == v3) /* v already in @surface */
00541     return NULL;
00542   if (GTS_POINT (v)->x == GTS_POINT (v1)->x &&
00543       GTS_POINT (v)->y == GTS_POINT (v1)->y)
00544     return v1;
00545   if (GTS_POINT (v)->x == GTS_POINT (v2)->x &&
00546       GTS_POINT (v)->y == GTS_POINT (v2)->y)
00547     return v2;
00548   if (GTS_POINT (v)->x == GTS_POINT (v3)->x &&
00549       GTS_POINT (v)->y == GTS_POINT (v3)->y)
00550     return v3;
00551 
00552   s4 = gts_vertices_are_connected (v, v1);
00553   if (!GTS_IS_EDGE (s4))
00554     e4 = gts_edge_new (surface->edge_class, v, v1);
00555   else
00556     e4 = GTS_EDGE (s4);
00557   s5 = gts_vertices_are_connected (v, v2);
00558   if (!GTS_IS_EDGE (s5))
00559     e5 = gts_edge_new (surface->edge_class, v, v2);
00560   else
00561     e5 = GTS_EDGE (s5);
00562   s6 = gts_vertices_are_connected (v, v3);
00563   if (!GTS_IS_EDGE (s6))
00564     e6 = gts_edge_new (surface->edge_class, v, v3);
00565   else
00566     e6 = GTS_EDGE (s6);
00567 
00568   /* cf. figure misc/swap.fig */
00569   nf[0] = gts_face_new (surface->face_class, e4, e1, e5);
00570   gts_object_attributes (GTS_OBJECT (nf[0]), GTS_OBJECT (f));
00571   nf[1] = gts_face_new (surface->face_class, e5, e2, e6);
00572   gts_object_attributes (GTS_OBJECT (nf[1]), GTS_OBJECT (f));
00573   nf[2] = gts_face_new (surface->face_class, e6, e3, e4);
00574   gts_object_attributes (GTS_OBJECT (nf[2]), GTS_OBJECT (f));
00575 
00576   if (GTS_IS_LIST_FACE (f) && GTS_IS_LIST_FACE (nf[0])) {
00577     GSList * i = GTS_LIST_FACE (f)->points, * last[3] = { NULL, NULL, NULL };
00578 
00579     while (i) {
00580       GtsPoint * p = i->data;
00581       GSList * next = i->next;
00582       guint j;
00583       
00584       if (p != GTS_POINT (v)) {
00585         if (gts_point_orientation (GTS_POINT (v), GTS_POINT (v1), p) >= 0.) {
00586           gdouble o = gts_point_orientation (GTS_POINT (v), GTS_POINT (v2), p);
00587 
00588           if (o != 0.)
00589             j = o > 0. ? 1 : 0;
00590           else
00591             j = gts_point_orientation (GTS_POINT (v), GTS_POINT (v3), p) 
00592               > 0. ? 0 : 1;
00593         }
00594         else if (gts_point_orientation (GTS_POINT (v), GTS_POINT (v3), p) > 0.)
00595           j = 2;
00596         else
00597           j = 1;
00598         if (last[j])
00599           last[j]->next = i; 
00600         else 
00601           GTS_LIST_FACE (nf[j])->points = i;
00602         last[j] = i;
00603       }
00604       else
00605         g_slist_free_1 (i);
00606       i = next;
00607     }
00608     GTS_LIST_FACE (f)->points = NULL;
00609     if (last[0]) last[0]->next = NULL;
00610     if (last[1]) last[1]->next = NULL;
00611     if (last[2]) last[2]->next = NULL;
00612   }
00613 
00614   gts_surface_remove_face (surface, f);
00615   gts_surface_add_face (surface, nf[0]);
00616   gts_surface_add_face (surface, nf[1]);
00617   gts_surface_add_face (surface, nf[2]);
00618 
00619   swap_if_in_circle (nf[0], v1, v2, v, e1, e5, e4, surface);
00620   swap_if_in_circle (nf[1], v2, v3, v, e2, e6, e5, surface);
00621   swap_if_in_circle (nf[2], v3, v1, v, e3, e4, e6, surface);
00622 
00623   return NULL;
00624 }
00625 
00642 GtsVertex * gts_delaunay_add_vertex (GtsSurface * surface, 
00643                                      GtsVertex * v,
00644                                      GtsFace * guess)
00645 {
00646   GtsFace * f;
00647 
00648   g_return_val_if_fail (surface != NULL, v);
00649   g_return_val_if_fail (v != NULL, v);
00650 
00651   if (!(f = gts_point_locate (GTS_POINT (v), surface, guess)))
00652     return v;
00653   return gts_delaunay_add_vertex_to_face (surface, v, f);
00654 }
00655 
00656 static gboolean polygon_in_circle (GSList * poly,
00657                                    GtsPoint * p1, 
00658                                    GtsPoint * p2,
00659                                    GtsPoint * p3)
00660 {
00661   GtsVertex * v1 = NULL, * v2 = NULL;
00662 
00663   while (poly) {
00664     GtsSegment * s = poly->data;
00665     GtsVertex * v;
00666     v = s->v1;
00667     if (v != v1 && v != v2 &&
00668         v != GTS_VERTEX (p1) &&
00669         v != GTS_VERTEX (p2) &&
00670         v != GTS_VERTEX (p3) &&
00671         gts_point_in_circle (GTS_POINT (v), p1, p2, p3) > 0.)
00672       return TRUE;
00673     v = s->v2;
00674     if (v != v1 && v != v2 &&
00675         v != GTS_VERTEX (p1) &&
00676         v != GTS_VERTEX (p2) &&
00677         v != GTS_VERTEX (p3) &&
00678         gts_point_in_circle (GTS_POINT (v), p1, p2, p3) > 0.)
00679       return TRUE;
00680     v1 = s->v1;
00681     v2 = s->v2;
00682     poly = poly->next;
00683   }
00684   return FALSE;
00685 }
00686 
00687 static void triangulate_polygon (GSList * poly, 
00688                                  GtsSurface * surface,
00689                                  GtsFace * ref)
00690 {
00691   GSList * i, * poly1, * poly2;
00692   GtsVertex * v1, * v2, * v3 = NULL;
00693   gboolean found = FALSE;
00694   GtsSegment * s, * s1, * s2;
00695   GtsEdge * e1, * e2;
00696   GtsFace * f;
00697 
00698   if (poly == NULL || poly->next == NULL) {
00699     g_slist_free (poly);
00700     return;
00701   }
00702 
00703   s = poly->data;
00704   s1 = poly->next->data;
00705   if (s->v1 == s1->v1 || s->v1 == s1->v2) {
00706     v1 = s->v2;
00707     v2 = s->v1;
00708   }
00709   else {
00710     g_assert (s->v2 == s1->v1 || s->v2 == s1->v2);
00711     v1 = s->v1;
00712     v2 = s->v2;
00713   }
00714 
00715   i = poly->next;
00716   v3 = v2;
00717   while (i && !found) {
00718     s1 = i->data;
00719     if (s1->v1 == v3)
00720       v3 = s1->v2;
00721     else {
00722       g_assert (s1->v2 == v3);
00723       v3 = s1->v1;
00724     }
00725     if (v3 != v1 &&
00726         gts_point_orientation (GTS_POINT (v1), 
00727                                GTS_POINT (v2), 
00728                                GTS_POINT (v3)) >= 0. &&
00729         !polygon_in_circle (poly, 
00730                             GTS_POINT (v1), 
00731                             GTS_POINT (v2), 
00732                             GTS_POINT (v3)))
00733       found = TRUE;
00734     else 
00735       i = i->next;
00736   }
00737 
00738   if (!found) {
00739     g_slist_free (poly);
00740     return;
00741   }
00742 
00743   s1 = gts_vertices_are_connected (v2, v3);
00744   if (!GTS_IS_EDGE (s1))
00745     e1 = gts_edge_new (surface->edge_class, v2, v3);
00746   else
00747     e1 = GTS_EDGE (s1);
00748   s2 = gts_vertices_are_connected (v3, v1);
00749   if (!GTS_IS_EDGE (s2))
00750     e2 = gts_edge_new (surface->edge_class, v3, v1);
00751   else
00752     e2 = GTS_EDGE (s2);
00753   f = gts_face_new (surface->face_class, GTS_EDGE (s), e1, e2);
00754   gts_object_attributes (GTS_OBJECT (f), GTS_OBJECT (ref));
00755   gts_surface_add_face (surface, f);
00756 
00757   poly1 = poly->next;
00758   g_slist_free_1 (poly);
00759   if (i->next && e2 != i->next->data)
00760     poly2 = g_slist_prepend (i->next, e2);
00761   else
00762     poly2 = i->next;
00763   if (e1 != i->data)
00764     i->next = g_slist_prepend (NULL, e1);
00765   else
00766     i->next = NULL;
00767 
00768  triangulate_polygon (poly1, surface, ref);
00769  triangulate_polygon (poly2, surface, ref);
00770 }
00771 
00782 void gts_delaunay_remove_vertex (GtsSurface * surface, GtsVertex * v)
00783 {
00784   GSList * triangles, * i;
00785   GtsFace * ref = NULL;
00786 
00787   g_return_if_fail (surface != NULL);
00788   g_return_if_fail (v != NULL);
00789 
00790   i = triangles = gts_vertex_triangles (v, NULL);
00791   while (i && !ref) {
00792     if (GTS_IS_FACE (i->data) &&
00793         gts_face_has_parent_surface (i->data, surface))
00794       ref = i->data;
00795     i = i->next;
00796   }
00797   if (!ref) {
00798     g_slist_free (triangles);
00799     g_return_if_fail (ref);
00800   }
00801   triangulate_polygon (gts_vertex_fan_oriented (v, surface), surface, ref);
00802   i = triangles;
00803   while (i) {
00804     if (GTS_IS_FACE (i->data) &&
00805         gts_face_has_parent_surface (i->data, surface))
00806       gts_surface_remove_face (surface, i->data);
00807     i = i->next;
00808   }
00809   g_slist_free (triangles);
00810 }
00811 
00812 #define NEXT_CUT(edge, edge1, list) { next = neighbor (f, edge, surface);\
00813                                       remove_triangles (e, surface);\
00814                                       if (!constraint && !e->triangles)\
00815                                         gts_object_destroy (GTS_OBJECT (e));\
00816                                       g_assert (next);\
00817                                       *list = g_slist_prepend (*list, edge1);\
00818                                       return g_slist_concat (constraint,\
00819                                         remove_intersected_edge (s, edge,\
00820                                                next, surface, left, right));\
00821                                     }
00822 
00823 static void remove_triangles (GtsEdge * e, GtsSurface * s)
00824 {
00825   GSList * i = e->triangles;
00826 
00827   while (i) {
00828     GSList * next = i->next;
00829 
00830     if (GTS_IS_FACE (i->data) && gts_face_has_parent_surface (i->data, s))
00831       gts_surface_remove_face (s, i->data);
00832     i = next;
00833   }
00834 }
00835 
00836 static GSList * 
00837 remove_intersected_edge (GtsSegment * s,
00838                          GtsEdge * e,
00839                          GtsFace * f,
00840                          GtsSurface * surface,
00841                          GSList ** left, GSList ** right)
00842 {
00843   GtsVertex * v1, * v2, * v3;
00844   GtsEdge * e1, * e2;
00845   gdouble o1, o2;
00846   GtsFace * next;
00847   GSList * constraint = NULL;
00848 
00849   if (GTS_IS_CONSTRAINT (e))
00850     constraint = g_slist_prepend (NULL, e);
00851 
00852   gts_triangle_vertices_edges (GTS_TRIANGLE (f), e, 
00853                                &v1, &v2, &v3, &e, &e1, &e2);
00854   
00855   o1 = gts_point_orientation (GTS_POINT (v2), GTS_POINT (v3), 
00856                               GTS_POINT (s->v2));
00857   o2 = gts_point_orientation (GTS_POINT (v3), GTS_POINT (v1), 
00858                               GTS_POINT (s->v2));
00859 
00860   if (o1 == 0. && o2 == 0.) {
00861 /*    if(o2 != 0.) {
00862       fprintf(stderr, "o1 = %f o2 = %f\n", o1, o2);
00863       fprintf(stderr, "v1 = %f, %f\n", GTS_POINT(v1)->x, GTS_POINT(v1)->y);
00864       fprintf(stderr, "v2 = %f, %f\n", GTS_POINT(v2)->x, GTS_POINT(v2)->y);
00865       fprintf(stderr, "v3 = %f, %f\n", GTS_POINT(v3)->x, GTS_POINT(v3)->y);
00866       fprintf(stderr, "s->v2 = %f, %f\n", GTS_POINT(s->v2)->x, GTS_POINT(s->v2)->y);
00867 
00868       g_assert (o2 == 0.);
00869     }*/
00870  //   if(o2 == 0.) {
00871       remove_triangles (e, surface);
00872       if (!constraint && !e->triangles)
00873         gts_object_destroy (GTS_OBJECT (e));
00874       *left = g_slist_prepend (*left, e2);
00875       *right = g_slist_prepend (*right, e1);
00876 //    }
00877   }
00878   else if (o1 > 0.) {
00879     g_assert (o2 <= 0.);
00880     NEXT_CUT (e2, e1, right)
00881   }
00882   else if (o2 >= 0.)
00883     NEXT_CUT (e1, e2, left)
00884   else {
00885     gdouble o3 = gts_point_orientation (GTS_POINT (s->v1), GTS_POINT (s->v2),
00886                                         GTS_POINT (v3));
00887     if (o3 > 0.)
00888       NEXT_CUT (e1, e2, left)
00889     else
00890       NEXT_CUT (e2, e1, right)
00891   }
00892   return constraint;
00893 }
00894 
00895 static GSList * 
00896 remove_intersected_vertex (GtsSegment * s,
00897                            GtsVertex * v,
00898                            GtsSurface * surface,
00899                            GSList ** left,
00900                            GSList ** right,
00901                            GtsFace ** ref)
00902 {
00903   GSList * triangles = gts_vertex_triangles (v, NULL);
00904   GSList * i;
00905 
00906   i = triangles;
00907   while (i) {
00908     GtsTriangle * t = i->data;
00909     if (GTS_IS_FACE (t) && 
00910         gts_face_has_parent_surface (GTS_FACE (t), surface)) {
00911       GtsVertex * v1, * v2, * v3;
00912       gdouble o1, o2;
00913 
00914       gts_triangle_vertices (t, &v1, &v2, &v3);
00915       if (v == v2) {
00916         v2 = v3;
00917         v3 = v1;
00918       }
00919       else if (v == v3) {
00920         v3 = v2;
00921         v2 = v1;        
00922       }
00923       else
00924         g_assert (v == v1);
00925 
00926       if ((o1 = gts_point_orientation (GTS_POINT (v), GTS_POINT (v2),
00927                                        GTS_POINT (s->v2))) >= 0. &&
00928           (o2 = gts_point_orientation (GTS_POINT (v3), GTS_POINT (v),
00929                                        GTS_POINT (s->v2))) >= 0.) {
00930         gdouble o3 = gts_point_orientation (GTS_POINT (v2), GTS_POINT (v3),
00931                                             GTS_POINT (s->v2));
00932         GtsEdge * e = gts_triangle_edge_opposite (t, v);
00933         GtsEdge * e1, * e2;
00934         GtsFace * next = neighbor (GTS_FACE (t), e, surface);
00935 
00936         *ref = GTS_FACE (t);
00937         gts_triangle_vertices_edges (t, e, &v2, &v3, &v, &e, &e2, &e1);
00938 
00939         g_slist_free (triangles);
00940 
00941         if (o3 >= 0.) /* @s->v2 is inside (or on the edge) of t */
00942           return NULL;
00943 
00944         gts_allow_floating_faces = TRUE;
00945         gts_surface_remove_face (surface, GTS_FACE (t));
00946         gts_allow_floating_faces = FALSE;
00947 
00948         *left = g_slist_prepend (*left, e2);
00949         *right = g_slist_prepend (*right, e1);
00950 
00951         g_assert (next);
00952         return remove_intersected_edge (s, e, next, surface, left, right);
00953       }
00954     }
00955     i = i->next;
00956   }
00957 
00958   g_assert_not_reached ();
00959   return NULL;
00960 }
00961 
00973 GSList * gts_delaunay_add_constraint (GtsSurface * surface,
00974                                       GtsConstraint * c)
00975 {
00976   GSList * constraints;
00977   GtsVertex * v1; //, * v2;
00978   GSList * left = NULL, * right = NULL;
00979   GtsFace * ref = NULL;
00980 
00981   g_return_val_if_fail (surface != NULL, NULL);
00982   g_return_val_if_fail (c != NULL, NULL);
00983   g_return_val_if_fail (GTS_IS_CONSTRAINT (c), NULL);
00984   
00985   v1 = GTS_SEGMENT (c)->v1;
00986   //v2 = GTS_SEGMENT (c)->v2;
00987   
00988   gts_allow_floating_edges = TRUE;
00989   constraints = remove_intersected_vertex (GTS_SEGMENT (c), v1, surface,
00990                                            &left, &right, &ref);
00991   gts_allow_floating_edges = FALSE;
00992 #if 1
00993   triangulate_polygon (g_slist_prepend (g_slist_reverse (right), c), 
00994                        surface, ref);
00995   triangulate_polygon (g_slist_prepend (left, c), 
00996                        surface, ref);
00997 #else
00998   right = g_slist_prepend (g_slist_reverse (right), c);
00999   left = g_slist_prepend (left, c);
01000   {
01001     FILE * fp0 = fopen ("hole", "wt");
01002     FILE * fp1 = fopen ("right", "wt");
01003     FILE * fp2 = fopen ("left", "wt");
01004     GSList * i = left;
01005 
01006     gts_surface_write (surface, fp0);
01007     fclose (fp0);
01008  
01009     fprintf (fp2, "LIST {\n");
01010     while (i) {
01011       GtsSegment * s = i->data;
01012       fprintf (fp2, 
01013                "# %p: %p->%p\n"
01014                "VECT 1 2 0 2 0 %g %g 0 %g %g 0\n",
01015                s, s->v1, s->v2,
01016                GTS_POINT (s->v1)->x, GTS_POINT (s->v1)->y,
01017                GTS_POINT (s->v2)->x, GTS_POINT (s->v2)->y);
01018       i = i->next;
01019     }
01020     fprintf (fp2, "}\n");
01021     fprintf (fp1, "LIST {\n");
01022     i = right;
01023     while (i) {
01024       GtsSegment * s = i->data;
01025       fprintf (fp1, 
01026                "# %p: %p->%p\n"
01027                "VECT 1 2 0 2 0 %g %g 0 %g %g 0\n",
01028                s, s->v1, s->v2,
01029                GTS_POINT (s->v1)->x, GTS_POINT (s->v1)->y,
01030                GTS_POINT (s->v2)->x, GTS_POINT (s->v2)->y);
01031       i = i->next;
01032     }
01033     fprintf (fp1, "}\n");
01034     fclose (fp1);
01035     fclose (fp2);
01036   }
01037   triangulate_polygon (right, surface);
01038   triangulate_polygon (left, surface);
01039 #endif
01040   if (ref && !ref->surfaces) {
01041     gts_allow_floating_edges = TRUE;
01042     gts_object_destroy (GTS_OBJECT (ref));
01043     gts_allow_floating_edges = FALSE;
01044   }
01045   return constraints;
01046 }
01047 
01048 static void delaunay_check (GtsTriangle * t, gpointer * data)
01049 {
01050   GtsSurface * surface = data[0];
01051   GtsFace ** face = data[1];
01052 
01053   if (*face == NULL) {
01054     GSList * i, * list;
01055     GtsVertex * v1, * v2, * v3;
01056 
01057     gts_triangle_vertices (t, &v1, &v2, &v3);
01058     list = gts_vertex_neighbors (v1, NULL, surface);
01059     list = gts_vertex_neighbors (v2, list, surface);
01060     list = gts_vertex_neighbors (v3, list, surface);
01061     i = list;
01062     while (i && *face == NULL) {
01063       GtsVertex * v = i->data;
01064       if (v != v1 && v != v2 && v != v3 &&
01065           gts_point_in_circle (GTS_POINT (v), 
01066                                GTS_POINT (v1),
01067                                GTS_POINT (v2),  
01068                                GTS_POINT (v3)) > 0.)
01069         *face = GTS_FACE (t);
01070       i = i->next;
01071     }
01072     g_slist_free (list);
01073   }
01074 }
01075 
01084 GtsFace * gts_delaunay_check (GtsSurface * surface)
01085 {
01086   GtsFace * face = NULL;
01087   gpointer data[2];
01088 
01089   g_return_val_if_fail (surface != NULL, FALSE);
01090 
01091   data[0] = surface;
01092   data[1] = &face;
01093   gts_surface_foreach_face (surface, (GtsFunc) delaunay_check, data);
01094 
01095   return face;
01096 }
01097 
01105 void gts_delaunay_remove_hull (GtsSurface * surface)
01106 {
01107   GSList * boundary;
01108 
01109   g_return_if_fail (surface != NULL);
01110 
01111   boundary = gts_surface_boundary (surface);
01112   gts_allow_floating_edges = TRUE;
01113   while (boundary) {
01114     GSList * i = boundary;
01115     GtsEdge * e = i->data;
01116 
01117     boundary = i->next;
01118     g_slist_free_1 (i);
01119     if (!GTS_IS_CONSTRAINT (e)) {
01120       GtsTriangle * t = GTS_TRIANGLE (gts_edge_is_boundary (e, surface));
01121 
01122       if (t != NULL) {
01123         if (t->e1 != e && !GTS_IS_CONSTRAINT (t->e1) &&
01124             !gts_edge_is_boundary (t->e1, surface))
01125           boundary = g_slist_prepend (boundary, t->e1);
01126         if (t->e2 != e && !GTS_IS_CONSTRAINT (t->e2) &&
01127             !gts_edge_is_boundary (t->e2, surface))
01128           boundary = g_slist_prepend (boundary, t->e2);
01129         if (t->e3 != e && !GTS_IS_CONSTRAINT (t->e3) &&
01130             !gts_edge_is_boundary (t->e3, surface))
01131           boundary = g_slist_prepend (boundary, t->e3);
01132         gts_surface_remove_face (surface, GTS_FACE (t));
01133       }
01134       if (!e->triangles)
01135         gts_object_destroy (GTS_OBJECT (e));
01136     }
01137   }
01138   gts_allow_floating_edges = FALSE;
01139 }
01140 
01141 /* GtsListFace: Object */
01142 
01143 static void gts_list_face_destroy (GtsObject * object)
01144 {
01145   g_slist_free (GTS_LIST_FACE (object)->points);
01146 
01147   (* GTS_OBJECT_CLASS (gts_list_face_class ())->parent_class->destroy) 
01148     (object);
01149 }
01150 
01151 static void gts_list_face_class_init (GtsFaceClass * klass)
01152 {
01153   GTS_OBJECT_CLASS (klass)->destroy = gts_list_face_destroy;
01154 }
01155 
01156 GtsFaceClass * gts_list_face_class (void)
01157 {
01158   static GtsFaceClass * klass = NULL;
01159 
01160   if (klass == NULL) {
01161     GtsObjectClassInfo gts_list_face_info = {
01162       "GtsListFace",
01163       sizeof (GtsListFace),
01164       sizeof (GtsFaceClass),
01165       (GtsObjectClassInitFunc) gts_list_face_class_init,
01166       (GtsObjectInitFunc) NULL,
01167       (GtsArgSetFunc) NULL,
01168       (GtsArgGetFunc) NULL
01169     };
01170     klass = gts_object_class_new (GTS_OBJECT_CLASS (gts_face_class ()),
01171                                   &gts_list_face_info);
01172   }
01173 
01174   return klass;
01175 }