pcb 4.1.1
An interactive printed circuit board layout editor.

bbtree.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 #include <math.h>
00021 #include "gts.h"
00022 
00023 static void bbox_init (GtsBBox * bbox)
00024 {
00025   bbox->bounded = NULL;
00026 }
00027 
00033 GtsBBoxClass * gts_bbox_class (void)
00034 {
00035   static GtsBBoxClass * klass = NULL;
00036 
00037   if (klass == NULL) {
00038     GtsObjectClassInfo bbox_info = {
00039       "GtsBBox",
00040       sizeof (GtsBBox),
00041       sizeof (GtsBBoxClass),
00042       (GtsObjectClassInitFunc) NULL,
00043       (GtsObjectInitFunc) bbox_init,
00044       (GtsArgSetFunc) NULL,
00045       (GtsArgGetFunc) NULL
00046     };
00047     klass = gts_object_class_new (gts_object_class (), &bbox_info);
00048   }
00049 
00050   return klass;
00051 }
00052 
00066 void gts_bbox_set (GtsBBox * bbox,
00067                    gpointer bounded,
00068                    gdouble x1, gdouble y1, gdouble z1,
00069                    gdouble x2, gdouble y2, gdouble z2)
00070 {
00071   g_return_if_fail (bbox != NULL);
00072   g_return_if_fail (x2 >= x1 && y2 >= y1 && z2 >= z1);
00073 
00074   bbox->x1 = x1; bbox->y1 = y1; bbox->z1 = z1;
00075   bbox->x2 = x2; bbox->y2 = y2; bbox->z2 = z2;
00076   bbox->bounded = bounded;
00077 }
00078 
00092 GtsBBox * gts_bbox_new (GtsBBoxClass * klass,
00093                         gpointer bounded,
00094                         gdouble x1, gdouble y1, gdouble z1,
00095                         gdouble x2, gdouble y2, gdouble z2)
00096 {
00097   GtsBBox * bbox;
00098 
00099   g_return_val_if_fail (klass != NULL, NULL);
00100 
00101   bbox = GTS_BBOX (gts_object_new (GTS_OBJECT_CLASS (klass)));
00102   gts_bbox_set (bbox, bounded, x1, y1, z1, x2, y2, z2);
00103   return bbox;
00104 }
00105 
00113 GtsBBox * gts_bbox_triangle (GtsBBoxClass * klass,
00114                              GtsTriangle * t)
00115 {
00116   GtsBBox * bbox;
00117   GtsPoint * p;
00118 
00119   g_return_val_if_fail (t != NULL, NULL);
00120   g_return_val_if_fail (klass != NULL, NULL);
00121 
00122   p = GTS_POINT (GTS_SEGMENT (t->e1)->v1);
00123   bbox = gts_bbox_new (klass, t, p->x, p->y, p->z, p->x, p->y, p->z);
00124 
00125   p = GTS_POINT (GTS_SEGMENT (t->e1)->v2);
00126   if (p->x > bbox->x2) bbox->x2 = p->x;
00127   if (p->x < bbox->x1) bbox->x1 = p->x;
00128   if (p->y > bbox->y2) bbox->y2 = p->y;
00129   if (p->y < bbox->y1) bbox->y1 = p->y;
00130   if (p->z > bbox->z2) bbox->z2 = p->z;
00131   if (p->z < bbox->z1) bbox->z1 = p->z;
00132   p = GTS_POINT (gts_triangle_vertex (t));
00133   if (p->x > bbox->x2) bbox->x2 = p->x;
00134   if (p->x < bbox->x1) bbox->x1 = p->x;
00135   if (p->y > bbox->y2) bbox->y2 = p->y;
00136   if (p->y < bbox->y1) bbox->y1 = p->y;
00137   if (p->z > bbox->z2) bbox->z2 = p->z;
00138   if (p->z < bbox->z1) bbox->z1 = p->z;
00139   
00140   return bbox;
00141 }
00142 
00150 GtsBBox * gts_bbox_segment (GtsBBoxClass * klass, GtsSegment * s)
00151 {
00152   GtsBBox * bbox;
00153   GtsPoint * p1, * p2;
00154 
00155   g_return_val_if_fail (s != NULL, NULL);
00156   g_return_val_if_fail (klass != NULL, NULL);
00157 
00158   bbox = gts_bbox_new (klass, s, 0., 0., 0., 0., 0., 0.);
00159 
00160   p1 = GTS_POINT (s->v1); 
00161   p2 = GTS_POINT (s->v2);
00162   if (p1->x > p2->x) {
00163     bbox->x2 = p1->x; bbox->x1 = p2->x;
00164   }
00165   else {
00166     bbox->x1 = p1->x; bbox->x2 = p2->x;
00167   }
00168   if (p1->y > p2->y) {
00169     bbox->y2 = p1->y; bbox->y1 = p2->y;
00170   }
00171   else {
00172     bbox->y1 = p1->y; bbox->y2 = p2->y;
00173   }
00174   if (p1->z > p2->z) {
00175     bbox->z2 = p1->z; bbox->z1 = p2->z;
00176   }
00177   else {
00178     bbox->z1 = p1->z; bbox->z2 = p2->z;
00179   }
00180 
00181   return bbox;
00182 }
00183 
00184 static void bbox_foreach_vertex (GtsPoint * p, GtsBBox * bb)
00185 {
00186   if (p->x < bb->x1) bb->x1 = p->x;
00187   if (p->y < bb->y1) bb->y1 = p->y;
00188   if (p->z < bb->z1) bb->z1 = p->z;
00189   if (p->x > bb->x2) bb->x2 = p->x;
00190   if (p->y > bb->y2) bb->y2 = p->y;
00191   if (p->z > bb->z2) bb->z2 = p->z;
00192 }
00193 
00201 GtsBBox * gts_bbox_surface (GtsBBoxClass * klass, GtsSurface * surface)
00202 {
00203   GtsBBox * bbox;
00204 
00205   g_return_val_if_fail (klass != NULL, NULL);
00206   g_return_val_if_fail (surface != NULL, NULL);
00207 
00208   bbox = gts_bbox_new (klass, surface, 0., 0., 0., 0., 0., 0.);
00209   bbox->x1 = bbox->y1 = bbox->z1 = G_MAXDOUBLE;
00210   bbox->x2 = bbox->y2 = bbox->z2 = -G_MAXDOUBLE;
00211 
00212   gts_surface_foreach_vertex (surface, (GtsFunc) bbox_foreach_vertex, bbox);
00213 
00214   return bbox;
00215 }
00216 
00225 GtsBBox * gts_bbox_bboxes (GtsBBoxClass * klass, GSList * bboxes)
00226 {
00227   GtsBBox * bbox;
00228   GtsBBox * bb;
00229 
00230   g_return_val_if_fail (bboxes != NULL, NULL);
00231   g_return_val_if_fail (klass != NULL, NULL);
00232 
00233   bb = bboxes->data;
00234   bbox = gts_bbox_new (klass, bboxes, 
00235                        bb->x1, bb->y1, bb->z1, bb->x2, bb->y2, bb->z2);
00236   bboxes = bboxes->next;
00237   while (bboxes) {
00238     bb = bboxes->data;
00239     if (bb->x1 < bbox->x1) bbox->x1 = bb->x1;
00240     if (bb->y1 < bbox->y1) bbox->y1 = bb->y1;
00241     if (bb->z1 < bbox->z1) bbox->z1 = bb->z1;
00242     if (bb->x2 > bbox->x2) bbox->x2 = bb->x2;
00243     if (bb->y2 > bbox->y2) bbox->y2 = bb->y2;
00244     if (bb->z2 > bbox->z2) bbox->z2 = bb->z2;
00245     bboxes = bboxes->next;
00246   }
00247 
00248   return bbox;
00249 }
00250 
00258 GtsBBox * gts_bbox_points (GtsBBoxClass * klass, GSList * points)
00259 {
00260   GtsPoint * p;
00261   GtsBBox * bbox;
00262   GSList * i;
00263 
00264   if (points == NULL) 
00265     return NULL;
00266 
00267   p = points->data;  
00268   bbox = gts_bbox_new (klass, points, p->x, p->y, p->z, p->x, p->y, p->z);
00269 
00270   i = points->next;
00271   while (i) {
00272     p = i->data;
00273     if (p->x > bbox->x2) 
00274       bbox->x2 = p->x;
00275     else if (p->x < bbox->x1) 
00276       bbox->x1 = p->x;
00277     if (p->y > bbox->y2) 
00278       bbox->y2 = p->y;
00279     else if (p->y < bbox->y1) 
00280       bbox->y1 = p->y;
00281     if (p->z > bbox->z2) 
00282       bbox->z2 = p->z;
00283     else if (p->z < bbox->z1) 
00284       bbox->z1 = p->z;
00285     i = i->next;
00286   }
00287   
00288   return bbox;
00289 }
00290 
00299 gboolean gts_bboxes_are_overlapping (GtsBBox * bb1, GtsBBox * bb2)
00300 {
00301   if (bb1 == bb2)
00302     return TRUE;
00303   if (bb1->x1 > bb2->x2)
00304     return FALSE;
00305   if (bb2->x1 > bb1->x2)
00306     return FALSE;
00307   if (bb1->y1 > bb2->y2)
00308     return FALSE;
00309   if (bb2->y1 > bb1->y2)
00310     return FALSE;
00311   if (bb1->z1 > bb2->z2)
00312     return FALSE;
00313   if (bb2->z1 > bb1->z2)
00314     return FALSE;  
00315   return TRUE;
00316 }
00317 
00318 #define bbox_volume(bb) (((bb)->x2 -\
00319                           (bb)->x1)*\
00320                          ((bb)->y2 -\
00321                           (bb)->y1)*\
00322                          ((bb)->z2 -\
00323                           (bb)->z1))
00324 
00331 gdouble gts_bbox_diagonal2 (GtsBBox * bb)
00332 {
00333   gdouble x, y, z;
00334 
00335   g_return_val_if_fail (bb != NULL, 0.);
00336 
00337   x = bb->x2 - bb->x1;
00338   y = bb->y2 - bb->y1;
00339   z = bb->z2 - bb->z1;
00340 
00341   return x*x + y*y + z*z;
00342 }
00343 
00351 void gts_bbox_draw (GtsBBox * bb, FILE * fptr)
00352 {
00353   g_return_if_fail (bb != NULL);
00354 
00355   fprintf (fptr, "OFF 8 6 12\n");
00356   fprintf (fptr, "%g %g %g\n",
00357            bb->x1, bb->y1, bb->z1);
00358   fprintf (fptr, "%g %g %g\n",
00359            bb->x2, bb->y1, bb->z1);
00360   fprintf (fptr, "%g %g %g\n",
00361            bb->x2, bb->y2, bb->z1);
00362   fprintf (fptr, "%g %g %g\n",
00363            bb->x1, bb->y2, bb->z1);
00364   fprintf (fptr, "%g %g %g\n",
00365            bb->x1, bb->y1, bb->z2);
00366   fprintf (fptr, "%g %g %g\n",
00367            bb->x2, bb->y1, bb->z2);
00368   fprintf (fptr, "%g %g %g\n",
00369            bb->x2, bb->y2, bb->z2);
00370   fprintf (fptr, "%g %g %g\n",
00371            bb->x1, bb->y2, bb->z2);
00372   fputs ("4 3 2 1 0\n"
00373          "4 4 5 6 7\n"
00374          "4 2 3 7 6\n"
00375          "4 0 1 5 4\n"
00376          "4 0 4 7 3\n"
00377          "4 1 2 6 5\n",
00378          fptr);
00379 }
00380 
00381 #define MINMAX(x1, x2, xmin, xmax) { if (x1 < x2) { xmin = x1; xmax = x2; }\
00382                                      else { xmin = x2; xmax = x1; } }
00383 
00397 void gts_bbox_point_distance2 (GtsBBox * bb, GtsPoint * p,
00398                                gdouble * min, gdouble * max)
00399 {
00400   gdouble x1, y1, z1, x2, y2, z2, x, y, z;
00401   gdouble dmin, dmax, xd1, xd2, yd1, yd2, zd1, zd2;
00402   gdouble mx, Mx, my, My, mz, Mz;
00403     
00404   g_return_if_fail (bb != NULL);
00405   g_return_if_fail (p != NULL);
00406   g_return_if_fail (min != NULL);
00407   g_return_if_fail (max != NULL);
00408 
00409   x1 = bb->x1; y1 = bb->y1; z1 = bb->z1; 
00410   x2 = bb->x2; y2 = bb->y2; z2 = bb->z2;
00411   x = p->x; y = p->y; z = p->z;
00412 
00413   xd1 = (x1 - x)*(x1 - x);
00414   xd2 = (x - x2)*(x - x2);
00415   yd1 = (y1 - y)*(y1 - y);
00416   yd2 = (y - y2)*(y - y2);
00417   zd1 = (z1 - z)*(z1 - z);
00418   zd2 = (z - z2)*(z - z2);
00419   
00420   dmin = x < x1 ? xd1 : x > x2 ? xd2 : 0.0;
00421   dmin += y < y1 ? yd1 : y > y2 ? yd2 : 0.0;
00422   dmin += z < z1 ? zd1 : z > z2 ? zd2 : 0.0;
00423 
00424   MINMAX (xd1, xd2, mx, Mx);
00425   MINMAX (yd1, yd2, my, My);
00426   MINMAX (zd1, zd2, mz, Mz);
00427   
00428   dmax = mx + My + Mz;
00429   dmax = MIN (dmax, Mx + my + Mz);
00430   dmax = MIN (dmax, Mx + My + mz);
00431   
00432   *min = dmin;
00433   *max = dmax;
00434 }
00435 
00444 gboolean gts_bbox_is_stabbed (GtsBBox * bb, GtsPoint * p)
00445 {
00446   g_return_val_if_fail (bb != NULL, FALSE);
00447   g_return_val_if_fail (p != NULL, FALSE);
00448 
00449   if (p->x > bb->x2 ||
00450       p->y < bb->y1 || p->y > bb->y2 ||
00451       p->z < bb->z1 || p->z > bb->z2)
00452     return FALSE;
00453   return TRUE;
00454 }
00455 
00456 extern int triBoxOverlap (double boxcenter[3],
00457                           double boxhalfsize[3],
00458                           double triverts[3][3]);
00459 
00470 gboolean gts_bbox_overlaps_triangle (GtsBBox * bb, GtsTriangle * t)
00471 {
00472   double bc[3], bh[3], tv[3][3];
00473   GtsPoint * p1, * p2, * p3;
00474 
00475   g_return_val_if_fail (bb != NULL, FALSE);
00476   g_return_val_if_fail (t != NULL, FALSE);
00477 
00478   bc[0] = (bb->x2 + bb->x1)/2.;
00479   bh[0] = (bb->x2 - bb->x1)/2.;
00480   bc[1] = (bb->y2 + bb->y1)/2.;
00481   bh[1] = (bb->y2 - bb->y1)/2.;
00482   bc[2] = (bb->z2 + bb->z1)/2.;
00483   bh[2] = (bb->z2 - bb->z1)/2.;
00484   p1 = GTS_POINT (GTS_SEGMENT (t->e1)->v1);
00485   p2 = GTS_POINT (GTS_SEGMENT (t->e1)->v2);
00486   p3 = GTS_POINT (gts_triangle_vertex (t));
00487   tv[0][0] = p1->x; tv[0][1] = p1->y; tv[0][2] = p1->z;
00488   tv[1][0] = p2->x; tv[1][1] = p2->y; tv[1][2] = p2->z;
00489   tv[2][0] = p3->x; tv[2][1] = p3->y; tv[2][2] = p3->z;
00490 
00491   return triBoxOverlap (bc, bh, tv);
00492 }
00493 
00504 gboolean gts_bbox_overlaps_segment (GtsBBox * bb, GtsSegment * s)
00505 {
00506   double bc[3], bh[3], tv[3][3];
00507   GtsPoint * p1, * p2, * p3;
00508 
00509   g_return_val_if_fail (bb != NULL, FALSE);
00510   g_return_val_if_fail (s != NULL, FALSE);
00511 
00512   bc[0] = (bb->x2 + bb->x1)/2.;
00513   bh[0] = (bb->x2 - bb->x1)/2.;
00514   bc[1] = (bb->y2 + bb->y1)/2.;
00515   bh[1] = (bb->y2 - bb->y1)/2.;
00516   bc[2] = (bb->z2 + bb->z1)/2.;
00517   bh[2] = (bb->z2 - bb->z1)/2.;
00518   p1 = GTS_POINT (s->v1);
00519   p2 = GTS_POINT (s->v2);
00520   p3 = p1;
00521   tv[0][0] = p1->x; tv[0][1] = p1->y; tv[0][2] = p1->z;
00522   tv[1][0] = p2->x; tv[1][1] = p2->y; tv[1][2] = p2->z;
00523   tv[2][0] = p3->x; tv[2][1] = p3->y; tv[2][2] = p3->z;
00524 
00525   return triBoxOverlap (bc, bh, tv);
00526 }
00527 
00543 GNode * gts_bb_tree_new (GSList * bboxes)
00544 {
00545   GSList * i, * positive = NULL, * negative = NULL;
00546   GNode * node;
00547   GtsBBox * bbox;
00548   guint dir, np = 0, nn = 0;
00549   gdouble * p1, * p2;
00550   gdouble cut;
00551   
00552   g_return_val_if_fail (bboxes != NULL, NULL);
00553 
00554   if (bboxes->next == NULL) /* leaf node */
00555     return g_node_new (bboxes->data);
00556 
00557   bbox = gts_bbox_bboxes (gts_bbox_class (), bboxes);
00558   node = g_node_new (bbox);
00559 
00560   if (bbox->x2 - bbox->x1 > bbox->y2 - bbox->y1) {
00561     if (bbox->z2 - bbox->z1 > bbox->x2 - bbox->x1)
00562       dir = 2;
00563     else
00564       dir = 0;
00565   }
00566   else if (bbox->z2 - bbox->z1 > bbox->y2 - bbox->y1)
00567     dir = 2;
00568   else
00569     dir = 1;
00570 
00571   p1 = (gdouble *) &bbox->x1;
00572   p2 = (gdouble *) &bbox->x2;
00573   cut = (p1[dir] + p2[dir])/2.;
00574   i = bboxes;
00575   while (i) {
00576     bbox = i->data; 
00577     p1 = (gdouble *) &bbox->x1;
00578     p2 = (gdouble *) &bbox->x2;
00579     if ((p1[dir] + p2[dir])/2. > cut) {
00580       positive = g_slist_prepend (positive, bbox);
00581       np++;
00582     }
00583     else {
00584       negative = g_slist_prepend (negative, bbox);
00585       nn++;
00586     }
00587     i = i->next;
00588   }
00589   if (!positive) {
00590     GSList * last = g_slist_nth (negative, (nn - 1)/2);
00591     positive = last->next;
00592     last->next = NULL;
00593   }
00594   else if (!negative) {
00595     GSList * last = g_slist_nth (positive, (np - 1)/2);
00596     negative = last->next;
00597     last->next = NULL;
00598   }
00599   g_node_prepend (node, gts_bb_tree_new (positive));
00600   g_slist_free (positive);
00601   g_node_prepend (node, gts_bb_tree_new (negative));
00602   g_slist_free (negative);
00603   
00604   return node;
00605 }
00606 
00607 static void prepend_triangle_bbox (GtsTriangle * t, GSList ** bboxes)
00608 {
00609   *bboxes = g_slist_prepend (*bboxes, 
00610                              gts_bbox_triangle (gts_bbox_class (), t));
00611 }
00612 
00619 GNode * gts_bb_tree_surface (GtsSurface * s)
00620 {
00621   GSList * bboxes = NULL;
00622   GNode * tree;
00623 
00624   g_return_val_if_fail (s != NULL, NULL);
00625 
00626   gts_surface_foreach_face (s, (GtsFunc) prepend_triangle_bbox, &bboxes);
00627   tree = gts_bb_tree_new (bboxes);
00628   g_slist_free (bboxes);
00629 
00630   return tree;
00631 }
00632 
00641 GSList * gts_bb_tree_stabbed (GNode * tree, GtsPoint * p)
00642 {
00643   GSList * list = NULL;
00644   GtsBBox * bb;
00645   GNode * i;
00646 
00647   g_return_val_if_fail (tree != NULL, NULL);
00648   g_return_val_if_fail (p != NULL, NULL);
00649 
00650   bb = tree->data;
00651   if (!gts_bbox_is_stabbed (bb, p))
00652     return NULL;
00653   if (tree->children == NULL) /* leaf node */
00654     return g_slist_prepend (NULL, bb);
00655   i = tree->children;
00656   while (i) {
00657     list = g_slist_concat (list, gts_bb_tree_stabbed (i, p));
00658     i = i->next;
00659   }
00660   return list;
00661 }
00662 
00670 GSList * gts_bb_tree_overlap (GNode * tree, GtsBBox * bbox)
00671 {
00672   GSList * list = NULL;
00673   GtsBBox * bb;
00674   GNode * i;
00675 
00676   g_return_val_if_fail (tree != NULL, NULL);
00677   g_return_val_if_fail (bbox != NULL, NULL);
00678 
00679   bb = tree->data;
00680   if (!gts_bboxes_are_overlapping (bbox, bb))
00681     return NULL;
00682   if (tree->children == NULL) /* leaf node */
00683     return g_slist_prepend (NULL, bb);
00684   i = tree->children;
00685   while (i) {
00686     list = g_slist_concat (list, gts_bb_tree_overlap (i, bbox));
00687     i = i->next;
00688   }
00689   return list;
00690 }
00691 
00699 gboolean gts_bb_tree_is_overlapping (GNode * tree, GtsBBox * bbox)
00700 {
00701   GtsBBox * bb;
00702   GNode * i;
00703 
00704   g_return_val_if_fail (tree != NULL, FALSE);
00705   g_return_val_if_fail (bbox != NULL, FALSE);
00706 
00707   bb = tree->data;
00708   if (!gts_bboxes_are_overlapping (bbox, bb))
00709     return FALSE;
00710   if (tree->children == NULL) /* leaf node */
00711     return TRUE;
00712   i = tree->children;
00713   while (i) {
00714     if (gts_bb_tree_is_overlapping (i, bbox))
00715       return TRUE;
00716     i = i->next;
00717   }
00718   return FALSE;
00719 }
00720 
00730 void gts_bb_tree_traverse_overlapping (GNode * tree1, GNode * tree2,
00731                                        GtsBBTreeTraverseFunc func,
00732                                        gpointer data)
00733 {
00734   GtsBBox * bb1, * bb2;
00735 
00736   g_return_if_fail (tree1 != NULL && tree2 != NULL);
00737 
00738   bb1 = tree1->data; bb2 = tree2->data;
00739   if (!gts_bboxes_are_overlapping (bb1, bb2))
00740     return;
00741 
00742   if (tree1->children == NULL && tree2->children == NULL)
00743     (*func) (tree1->data, tree2->data, data);
00744   else if (tree2->children == NULL || 
00745            (tree1->children != NULL && 
00746             bbox_volume (bb1) > bbox_volume (bb2))) {
00747     GNode * i = tree1->children;
00748     while (i) {
00749       gts_bb_tree_traverse_overlapping (i, tree2, func, data);
00750       i = i->next;
00751     }
00752   }
00753   else {
00754     GNode * i = tree2->children;
00755     while (i) {
00756       gts_bb_tree_traverse_overlapping (tree1, i, func, data);
00757       i = i->next;
00758     }
00759   }
00760 }
00761 
00771 void gts_bb_tree_draw (GNode * tree, guint depth, FILE * fptr)
00772 {
00773   guint d;
00774 
00775   g_return_if_fail (tree != NULL);
00776   g_return_if_fail (fptr != NULL);
00777 
00778   d = g_node_depth (tree);
00779 
00780   if (d == 1)
00781     fprintf (fptr, "{ LIST");
00782 
00783   if (d == depth)
00784     gts_bbox_draw (tree->data, fptr);
00785   else if (d < depth) {
00786     GNode * i = tree->children;
00787     while (i) {
00788       gts_bb_tree_draw (i, depth, fptr);
00789       i = i->next;
00790     }
00791   }
00792 
00793   if (d == 1)
00794     fprintf (fptr, "}\n");
00795 }
00796 
00797 static void bb_tree_free (GNode * tree, gboolean free_leaves)
00798 {
00799   GNode * i;
00800 
00801   g_return_if_fail (tree != NULL);
00802 
00803   if (!free_leaves && tree->children == NULL) /* leaf node */
00804     return;
00805 
00806   gts_object_destroy (tree->data);
00807 
00808   i = tree->children;
00809   while (i) {
00810     bb_tree_free (i, free_leaves);
00811     i = i->next;
00812   }
00813 }
00814 
00824 void gts_bb_tree_destroy (GNode * tree, gboolean free_leaves)
00825 {
00826   g_return_if_fail (tree != NULL);
00827   
00828   bb_tree_free (tree, free_leaves);
00829   g_node_destroy (tree);
00830 }
00831 
00832 static gdouble bb_tree_min_max (GNode * tree,
00833                                 GtsPoint * p,
00834                                 gdouble min_max,
00835                                 GSList ** list)
00836 {
00837   GNode * tree1, * tree2;
00838   gdouble min1, max1, min2, max2;
00839 
00840   if (tree->children == NULL) {
00841     *list = g_slist_prepend (*list, tree->data);
00842     return min_max;
00843   }
00844   tree1 = tree->children;
00845   gts_bbox_point_distance2 (tree1->data, p, &min1, &max1);
00846   if (max1 < min_max)
00847     min_max = max1;
00848 
00849   tree2 = tree1->next;
00850   gts_bbox_point_distance2 (tree2->data, p, &min2, &max2);
00851   if (max2 < min_max)
00852     min_max = max2;
00853 
00854   if (min1 < min2) {
00855     if (min1 <= min_max) {
00856       min_max = bb_tree_min_max (tree1, p, min_max, list);
00857       if (min2 <= min_max)
00858         min_max = bb_tree_min_max (tree2, p, min_max, list);
00859     }
00860   }
00861   else {
00862     if (min2 <= min_max) {
00863       min_max = bb_tree_min_max (tree2, p, min_max, list);
00864       if (min1 <= min_max)
00865         min_max = bb_tree_min_max (tree1, p, min_max, list);
00866     }
00867   }
00868 
00869   return min_max;
00870 }
00871 
00880 GSList * gts_bb_tree_point_closest_bboxes (GNode * tree, 
00881                                            GtsPoint * p)
00882 {
00883   gdouble min, min_max;
00884   GSList * list = NULL, * i, * prev = NULL;
00885 
00886   g_return_val_if_fail (tree != NULL, NULL);
00887   g_return_val_if_fail (p != NULL, NULL);
00888 
00889   gts_bbox_point_distance2 (tree->data, p, &min, &min_max);
00890   min_max = bb_tree_min_max (tree, p, min_max, &list);
00891 
00892   i = list;
00893   while (i) {
00894     GSList * next = i->next;
00895     gdouble min, max;
00896 
00897     gts_bbox_point_distance2 (i->data, p, &min, &max);
00898 
00899     if (min > min_max) {
00900       if (prev == NULL)
00901         list = next;
00902       else
00903         prev->next = next;
00904       g_slist_free_1 (i);
00905     }
00906     else
00907       prev = i;
00908     i = next;
00909   }
00910 
00911   return list;
00912 }
00913 
00925 gdouble gts_bb_tree_point_distance (GNode * tree, 
00926                                     GtsPoint * p,
00927                                     GtsBBoxDistFunc distance,
00928                                     GtsBBox ** bbox)
00929 {
00930   GSList * list, * i;
00931   gdouble dmin = G_MAXDOUBLE;
00932 
00933   g_return_val_if_fail (tree != NULL, dmin);
00934   g_return_val_if_fail (p != NULL, dmin);
00935   g_return_val_if_fail (distance != NULL, dmin);
00936 
00937   i = list = gts_bb_tree_point_closest_bboxes (tree, p);
00938   while (i) {
00939     gdouble d = (*distance) (p, GTS_BBOX (i->data)->bounded);
00940 
00941     if (fabs (d) < fabs (dmin)) {
00942       dmin = d;
00943       if (bbox)
00944         *bbox = i->data;
00945     }
00946     i = i->next;
00947   }
00948   g_slist_free (list);
00949 
00950   return dmin;
00951 }
00952 
00964 GtsPoint * gts_bb_tree_point_closest (GNode * tree, 
00965                                       GtsPoint * p,
00966                                       GtsBBoxClosestFunc closest,
00967                                       gdouble * distance)
00968 {
00969   GSList * list, * i;
00970   gdouble dmin = G_MAXDOUBLE;
00971   GtsPoint * np = NULL;
00972 
00973   g_return_val_if_fail (tree != NULL, NULL);
00974   g_return_val_if_fail (p != NULL, NULL);
00975   g_return_val_if_fail (closest != NULL, NULL);
00976 
00977   i = list = gts_bb_tree_point_closest_bboxes (tree, p);
00978   while (i) {
00979     GtsPoint * tp = (*closest) (p, GTS_BBOX (i->data)->bounded);
00980     gdouble d = gts_point_distance2 (tp, p);
00981 
00982     if (d < dmin) {
00983       if (np)
00984         gts_object_destroy (GTS_OBJECT (np));
00985       np = tp;
00986       dmin = d;
00987     }
00988     else
00989       gts_object_destroy (GTS_OBJECT (tp));
00990     i = i->next;
00991   }
00992   g_slist_free (list);
00993 
00994   if (distance)
00995     *distance = dmin;
00996 
00997   return np;  
00998 }
00999 
01015 void gts_bb_tree_triangle_distance (GNode * tree,
01016                                     GtsTriangle * t,
01017                                     GtsBBoxDistFunc distance,
01018                                     gdouble delta,
01019                                     GtsRange * range)
01020 {
01021   GtsPoint * p1, * p2, * p3, * p;
01022   GtsVector p1p2, p1p3;
01023   gdouble l1, t1, dt1;
01024   guint i, n1;
01025 
01026   g_return_if_fail (tree != NULL);
01027   g_return_if_fail (t != NULL);
01028   g_return_if_fail (distance != NULL);
01029   g_return_if_fail (delta > 0.);
01030   g_return_if_fail (range != NULL);
01031 
01032   gts_triangle_vertices (t, 
01033                          (GtsVertex **) &p1, 
01034                          (GtsVertex **) &p2, 
01035                          (GtsVertex **) &p3);
01036 
01037   gts_vector_init (p1p2, p1, p2);
01038   gts_vector_init (p1p3, p1, p3);
01039   gts_range_init (range);
01040   p = GTS_POINT (gts_object_new (GTS_OBJECT_CLASS (gts_point_class ())));
01041 
01042   l1 = sqrt (gts_vector_scalar (p1p2, p1p2));
01043   n1 = l1/delta + 1;
01044   dt1 = 1.0/(gdouble) n1;
01045   t1 = 0.0;
01046   for (i = 0; i <= n1; i++, t1 += dt1) {
01047     gdouble t2 = 1. - t1;
01048     gdouble x = t2*p1p3[0];
01049     gdouble y = t2*p1p3[1];
01050     gdouble z = t2*p1p3[2];
01051     gdouble l2 = sqrt (x*x + y*y + z*z);
01052     guint j, n2 = (guint) (l2/delta + 1);
01053     gdouble dt2 = t2/(gdouble) n2;
01054 
01055     x = t2*p1->x + t1*p2->x;
01056     y = t2*p1->y + t1*p2->y;
01057     z = t2*p1->z + t1*p2->z;
01058     
01059     t2 = 0.0;
01060     for (j = 0; j <= n2; j++, t2 += dt2) {
01061       p->x = x + t2*p1p3[0];
01062       p->y = y + t2*p1p3[1];
01063       p->z = z + t2*p1p3[2];
01064 
01065       gts_range_add_value (range,
01066                     gts_bb_tree_point_distance (tree, p, distance, NULL));
01067     }
01068   }
01069 
01070   gts_object_destroy (GTS_OBJECT (p));
01071   gts_range_update (range);
01072 }
01073 
01089 void gts_bb_tree_segment_distance (GNode * tree,
01090                                    GtsSegment * s,
01091                                    gdouble (*distance) (GtsPoint *, 
01092                                                         gpointer),
01093                                    gdouble delta,
01094                                    GtsRange * range)
01095 {
01096   GtsPoint * p1, * p2, * p;
01097   GtsVector p1p2;
01098   gdouble l, t, dt;
01099   guint i, n;
01100 
01101   g_return_if_fail (tree != NULL);
01102   g_return_if_fail (s != NULL);
01103   g_return_if_fail (distance != NULL);
01104   g_return_if_fail (delta > 0.);
01105   g_return_if_fail (range != NULL);
01106 
01107   p1 = GTS_POINT (s->v1);
01108   p2 = GTS_POINT (s->v2);
01109 
01110   gts_vector_init (p1p2, p1, p2);
01111   gts_range_init (range);
01112   p = GTS_POINT (gts_object_new (GTS_OBJECT_CLASS (gts_point_class())));
01113 
01114   l = sqrt (gts_vector_scalar (p1p2, p1p2));
01115   n = (guint) (l/delta + 1);
01116   dt = 1.0/(gdouble) n;
01117   t = 0.0;
01118   for (i = 0; i <= n; i++, t += dt) {
01119     p->x = p1->x + t*p1p2[0];
01120     p->y = p1->y + t*p1p2[1];
01121     p->z = p1->z + t*p1p2[2];
01122     
01123     gts_range_add_value (range,
01124                          gts_bb_tree_point_distance (tree, p, distance, NULL));
01125   }
01126 
01127   gts_object_destroy (GTS_OBJECT (p));
01128   gts_range_update (range);
01129 }
01130 
01131 static void surface_distance_foreach_triangle (GtsTriangle * t, 
01132                                                gpointer * data)
01133 {
01134   gdouble * delta = data[1];
01135   GtsRange * range = data[2];
01136   gdouble * total_area = data[3], area;
01137   GtsRange range_triangle;
01138 
01139   gts_bb_tree_triangle_distance (data[0], t, data[4], *delta, &range_triangle);
01140 
01141   if (range_triangle.min < range->min)
01142     range->min = range_triangle.min;
01143   if (range_triangle.max > range->max)
01144     range->max = range_triangle.max;
01145   range->n += range_triangle.n;
01146 
01147   area = gts_triangle_area (t);
01148   *total_area += area;
01149   range->sum += area*range_triangle.mean;
01150   range->sum2 += area*range_triangle.mean*range_triangle.mean;
01151 }
01152 
01170 void gts_bb_tree_surface_distance (GNode * tree,
01171                                    GtsSurface * s,
01172                                    GtsBBoxDistFunc distance,
01173                                    gdouble delta,
01174                                    GtsRange * range)
01175 {
01176   gpointer data[5];
01177   gdouble total_area = 0.;
01178 
01179   g_return_if_fail (tree != NULL);
01180   g_return_if_fail (s != NULL);
01181   g_return_if_fail (delta > 0. && delta < 1.);
01182   g_return_if_fail (range != NULL);
01183 
01184   gts_range_init (range);
01185   delta *= sqrt (gts_bbox_diagonal2 (tree->data));
01186   data[0] = tree;
01187   data[1] = &delta;
01188   data[2] = range;
01189   data[3] = &total_area;
01190   data[4] = distance;
01191 
01192   gts_surface_foreach_face (s, 
01193                             (GtsFunc) surface_distance_foreach_triangle, 
01194                             data);
01195 
01196   if (total_area > 0.) {
01197     if (range->sum2 - range->sum*range->sum/total_area >= 0.)
01198       range->stddev = sqrt ((range->sum2 - range->sum*range->sum/total_area)
01199                             /total_area);
01200     else
01201       range->stddev = 0.;
01202     range->mean = range->sum/total_area;
01203   }
01204   else
01205     range->min = range->max = range->mean = range->stddev = 0.;
01206 }
01207 
01208 static void surface_distance_foreach_boundary (GtsEdge * e,
01209                                                gpointer * data)
01210 {
01211   gdouble * delta = data[1];
01212   GtsRange * range = data[2];
01213   gdouble * total_length = data[3], length;
01214   GtsRange range_edge;
01215 
01216   if (gts_edge_is_boundary (e, NULL)) {
01217     GtsSegment * s =  GTS_SEGMENT (e);
01218 
01219     gts_bb_tree_segment_distance (data[0], s, data[4], *delta, &range_edge);
01220 
01221     if (range_edge.min < range->min)
01222       range->min = range_edge.min;
01223     if (range_edge.max > range->max)
01224       range->max = range_edge.max;
01225     range->n += range_edge.n;
01226     
01227     length = gts_point_distance (GTS_POINT (s->v1), GTS_POINT (s->v2));
01228     *total_length += length;
01229     range->sum += length*range_edge.mean;
01230     range->sum2 += length*range_edge.mean*range_edge.mean;
01231   }
01232 }
01233 
01251 void gts_bb_tree_surface_boundary_distance (GNode * tree,
01252                                             GtsSurface * s,
01253                                             gdouble (*distance) (GtsPoint *,
01254                                                                  gpointer),
01255                                             gdouble delta,
01256                                             GtsRange * range)
01257 {
01258   gpointer data[5];
01259   gdouble total_length = 0.;
01260 
01261   g_return_if_fail (tree != NULL);
01262   g_return_if_fail (s != NULL);
01263   g_return_if_fail (delta > 0. && delta < 1.);
01264   g_return_if_fail (range != NULL);
01265 
01266   gts_range_init (range);
01267   delta *= sqrt (gts_bbox_diagonal2 (tree->data));
01268   data[0] = tree;
01269   data[1] = &delta;
01270   data[2] = range;
01271   data[3] = &total_length;
01272   data[4] = distance;
01273 
01274   gts_surface_foreach_edge (s, 
01275                             (GtsFunc) surface_distance_foreach_boundary, 
01276                             data);
01277 
01278   if (total_length > 0.) {
01279     if (range->sum2 - range->sum*range->sum/total_length >= 0.)
01280       range->stddev = sqrt ((range->sum2 - 
01281                              range->sum*range->sum/total_length)
01282                             /total_length);
01283     else
01284       range->stddev = 0.;
01285     range->mean = range->sum/total_length;
01286   }
01287   else
01288     range->min = range->max = range->mean = range->stddev = 0.;
01289 }