/* GTS - Library for the manipulation of triangulated surfaces
* Copyright (C) 1999 Stéphane Popinet
*
* This library is free software; you can redistribute it and/or
* modify it under the terms of the GNU Library General Public
* License as published by the Free Software Foundation; either
* version 2 of the License, or (at your option) any later version.
*
* This library is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* Library General Public License for more details.
*
* You should have received a copy of the GNU Library General Public
* License along with this library; if not, write to the
* Free Software Foundation, Inc., 59 Temple Place - Suite 330,
* Boston, MA 02111-1307, USA.
*/
#include <math.h>
#include "gts.h"
static void bbox_init (GtsBBox * bbox)
{
bbox->bounded = NULL;
}
/**
* gts_bbox_class:
*
* Returns: the #GtsBBoxClass.
*/
GtsBBoxClass * gts_bbox_class (void)
{
static GtsBBoxClass * klass = NULL;
if (klass == NULL) {
GtsObjectClassInfo bbox_info = {
"GtsBBox",
sizeof (GtsBBox),
sizeof (GtsBBoxClass),
(GtsObjectClassInitFunc) NULL,
(GtsObjectInitFunc) bbox_init,
(GtsArgSetFunc) NULL,
(GtsArgGetFunc) NULL
};
klass = gts_object_class_new (gts_object_class (), &bbox_info);
}
return klass;
}
/**
* gts_bbox_set:
* @bbox: a #GtsBBox.
* @bounded: the object to be bounded.
* @x1: x-coordinate of the lower left corner.
* @y1: y-coordinate of the lower left corner.
* @z1: z-coordinate of the lower left corner.
* @x2: x-coordinate of the upper right corner.
* @y2: y-coordinate of the upper right corner.
* @z2: z-coordinate of the upper right corner.
*
* Sets fields of @bbox.
*/
void gts_bbox_set (GtsBBox * bbox,
gpointer bounded,
gdouble x1, gdouble y1, gdouble z1,
gdouble x2, gdouble y2, gdouble z2)
{
g_return_if_fail (bbox != NULL);
g_return_if_fail (x2 >= x1 && y2 >= y1 && z2 >= z1);
bbox->x1 = x1; bbox->y1 = y1; bbox->z1 = z1;
bbox->x2 = x2; bbox->y2 = y2; bbox->z2 = z2;
bbox->bounded = bounded;
}
/**
* gts_bbox_new:
* @klass: a #GtsBBoxClass.
* @bounded: the object to be bounded.
* @x1: x-coordinate of the lower left corner.
* @y1: y-coordinate of the lower left corner.
* @z1: z-coordinate of the lower left corner.
* @x2: x-coordinate of the upper right corner.
* @y2: y-coordinate of the upper right corner.
* @z2: z-coordinate of the upper right corner.
*
* Returns: a new #GtsBBox.
*/
GtsBBox * gts_bbox_new (GtsBBoxClass * klass,
gpointer bounded,
gdouble x1, gdouble y1, gdouble z1,
gdouble x2, gdouble y2, gdouble z2)
{
GtsBBox * bbox;
g_return_val_if_fail (klass != NULL, NULL);
bbox = GTS_BBOX (gts_object_new (GTS_OBJECT_CLASS (klass)));
gts_bbox_set (bbox, bounded, x1, y1, z1, x2, y2, z2);
return bbox;
}
/**
* gts_bbox_triangle:
* @klass: a #GtsBBoxClass.
* @t: a #GtsTriangle.
*
* Returns: a new #GtsBBox bounding box of @t.
*/
GtsBBox * gts_bbox_triangle (GtsBBoxClass * klass,
GtsTriangle * t)
{
GtsBBox * bbox;
GtsPoint * p;
g_return_val_if_fail (t != NULL, NULL);
g_return_val_if_fail (klass != NULL, NULL);
p = GTS_POINT (GTS_SEGMENT (t->e1)->v1);
bbox = gts_bbox_new (klass, t, p->x, p->y, p->z, p->x, p->y, p->z);
p = GTS_POINT (GTS_SEGMENT (t->e1)->v2);
if (p->x > bbox->x2) bbox->x2 = p->x;
if (p->x < bbox->x1) bbox->x1 = p->x;
if (p->y > bbox->y2) bbox->y2 = p->y;
if (p->y < bbox->y1) bbox->y1 = p->y;
if (p->z > bbox->z2) bbox->z2 = p->z;
if (p->z < bbox->z1) bbox->z1 = p->z;
p = GTS_POINT (gts_triangle_vertex (t));
if (p->x > bbox->x2) bbox->x2 = p->x;
if (p->x < bbox->x1) bbox->x1 = p->x;
if (p->y > bbox->y2) bbox->y2 = p->y;
if (p->y < bbox->y1) bbox->y1 = p->y;
if (p->z > bbox->z2) bbox->z2 = p->z;
if (p->z < bbox->z1) bbox->z1 = p->z;
return bbox;
}
/**
* gts_bbox_segment:
* @klass: a #GtsBBoxClass.
* @s: a #GtsSegment.
*
* Returns: a new #GtsBBox bounding box of @s.
*/
GtsBBox * gts_bbox_segment (GtsBBoxClass * klass, GtsSegment * s)
{
GtsBBox * bbox;
GtsPoint * p1, * p2;
g_return_val_if_fail (s != NULL, NULL);
g_return_val_if_fail (klass != NULL, NULL);
bbox = gts_bbox_new (klass, s, 0., 0., 0., 0., 0., 0.);
p1 = GTS_POINT (s->v1);
p2 = GTS_POINT (s->v2);
if (p1->x > p2->x) {
bbox->x2 = p1->x; bbox->x1 = p2->x;
}
else {
bbox->x1 = p1->x; bbox->x2 = p2->x;
}
if (p1->y > p2->y) {
bbox->y2 = p1->y; bbox->y1 = p2->y;
}
else {
bbox->y1 = p1->y; bbox->y2 = p2->y;
}
if (p1->z > p2->z) {
bbox->z2 = p1->z; bbox->z1 = p2->z;
}
else {
bbox->z1 = p1->z; bbox->z2 = p2->z;
}
return bbox;
}
static void bbox_foreach_vertex (GtsPoint * p, GtsBBox * bb)
{
if (p->x < bb->x1) bb->x1 = p->x;
if (p->y < bb->y1) bb->y1 = p->y;
if (p->z < bb->z1) bb->z1 = p->z;
if (p->x > bb->x2) bb->x2 = p->x;
if (p->y > bb->y2) bb->y2 = p->y;
if (p->z > bb->z2) bb->z2 = p->z;
}
/**
* gts_bbox_surface:
* @klass: a #GtsBBoxClass.
* @surface: a #GtsSurface.
*
* Returns: a new #GtsBBox bounding box of @surface.
*/
GtsBBox * gts_bbox_surface (GtsBBoxClass * klass, GtsSurface * surface)
{
GtsBBox * bbox;
g_return_val_if_fail (klass != NULL, NULL);
g_return_val_if_fail (surface != NULL, NULL);
bbox = gts_bbox_new (klass, surface, 0., 0., 0., 0., 0., 0.);
bbox->x1 = bbox->y1 = bbox->z1 = G_MAXDOUBLE;
bbox->x2 = bbox->y2 = bbox->z2 = -G_MAXDOUBLE;
gts_surface_foreach_vertex (surface, (GtsFunc) bbox_foreach_vertex, bbox);
return bbox;
}
/**
* gts_bbox_bboxes:
* @klass: a #GtsBBoxClass.
* @bboxes: a list of #GtsBBox.
*
* Returns: a new #GtsBBox bounding box of all the bounding boxes in
* @bboxes.
*/
GtsBBox * gts_bbox_bboxes (GtsBBoxClass * klass, GSList * bboxes)
{
GtsBBox * bbox;
GtsBBox * bb;
g_return_val_if_fail (bboxes != NULL, NULL);
g_return_val_if_fail (klass != NULL, NULL);
bb = bboxes->data;
bbox = gts_bbox_new (klass, bboxes,
bb->x1, bb->y1, bb->z1, bb->x2, bb->y2, bb->z2);
bboxes = bboxes->next;
while (bboxes) {
bb = bboxes->data;
if (bb->x1 < bbox->x1) bbox->x1 = bb->x1;
if (bb->y1 < bbox->y1) bbox->y1 = bb->y1;
if (bb->z1 < bbox->z1) bbox->z1 = bb->z1;
if (bb->x2 > bbox->x2) bbox->x2 = bb->x2;
if (bb->y2 > bbox->y2) bbox->y2 = bb->y2;
if (bb->z2 > bbox->z2) bbox->z2 = bb->z2;
bboxes = bboxes->next;
}
return bbox;
}
/**
* gts_bbox_points:
* @klass: a #GtsBBoxClass.
* @points: a list of #GtsPoint.
*
* Returns: a new #GtsBBox bounding box of @points.
*/
GtsBBox * gts_bbox_points (GtsBBoxClass * klass, GSList * points)
{
GtsPoint * p;
GtsBBox * bbox;
GSList * i;
if (points == NULL)
return NULL;
p = points->data;
bbox = gts_bbox_new (klass, points, p->x, p->y, p->z, p->x, p->y, p->z);
i = points->next;
while (i) {
p = i->data;
if (p->x > bbox->x2)
bbox->x2 = p->x;
else if (p->x < bbox->x1)
bbox->x1 = p->x;
if (p->y > bbox->y2)
bbox->y2 = p->y;
else if (p->y < bbox->y1)
bbox->y1 = p->y;
if (p->z > bbox->z2)
bbox->z2 = p->z;
else if (p->z < bbox->z1)
bbox->z1 = p->z;
i = i->next;
}
return bbox;
}
/**
* gts_bboxes_are_overlapping:
* @bb1: a #GtsBBox.
* @bb2: a #GtsBBox.
*
* Returns: %TRUE if the bounding boxes @bb1 and @bb2 are overlapping
* (including just touching), %FALSE otherwise.
*/
gboolean gts_bboxes_are_overlapping (GtsBBox * bb1, GtsBBox * bb2)
{
if (bb1 == bb2)
return TRUE;
if (bb1->x1 > bb2->x2)
return FALSE;
if (bb2->x1 > bb1->x2)
return FALSE;
if (bb1->y1 > bb2->y2)
return FALSE;
if (bb2->y1 > bb1->y2)
return FALSE;
if (bb1->z1 > bb2->z2)
return FALSE;
if (bb2->z1 > bb1->z2)
return FALSE;
return TRUE;
}
#define bbox_volume(bb) (((bb)->x2 -\
(bb)->x1)*\
((bb)->y2 -\
(bb)->y1)*\
((bb)->z2 -\
(bb)->z1))
/**
* gts_bbox_diagonal2:
* @bb: a #GtsBBox.
*
* Returns: the squared length of the diagonal of @bb.
*/
gdouble gts_bbox_diagonal2 (GtsBBox * bb)
{
gdouble x, y, z;
g_return_val_if_fail (bb != NULL, 0.);
x = bb->x2 - bb->x1;
y = bb->y2 - bb->y1;
z = bb->z2 - bb->z1;
return x*x + y*y + z*z;
}
/**
* gts_bbox_draw:
* @bb: a #GtsBBox.
* @fptr: a file pointer.
*
* Writes in file @fptr an OOGL (Geomview) description of @bb.
*/
void gts_bbox_draw (GtsBBox * bb, FILE * fptr)
{
g_return_if_fail (bb != NULL);
fprintf (fptr, "OFF 8 6 12\n");
fprintf (fptr, "%g %g %g\n",
bb->x1, bb->y1, bb->z1);
fprintf (fptr, "%g %g %g\n",
bb->x2, bb->y1, bb->z1);
fprintf (fptr, "%g %g %g\n",
bb->x2, bb->y2, bb->z1);
fprintf (fptr, "%g %g %g\n",
bb->x1, bb->y2, bb->z1);
fprintf (fptr, "%g %g %g\n",
bb->x1, bb->y1, bb->z2);
fprintf (fptr, "%g %g %g\n",
bb->x2, bb->y1, bb->z2);
fprintf (fptr, "%g %g %g\n",
bb->x2, bb->y2, bb->z2);
fprintf (fptr, "%g %g %g\n",
bb->x1, bb->y2, bb->z2);
fputs ("4 3 2 1 0\n"
"4 4 5 6 7\n"
"4 2 3 7 6\n"
"4 0 1 5 4\n"
"4 0 4 7 3\n"
"4 1 2 6 5\n",
fptr);
}
#define MINMAX(x1, x2, xmin, xmax) { if (x1 < x2) { xmin = x1; xmax = x2; }\
else { xmin = x2; xmax = x1; } }
/**
* gts_bbox_point_distance2:
* @bb: a #GtsBBox.
* @p: a #GtsPoint.
* @min: a pointer on a gdouble.
* @max: a pointer on a gdouble.
*
* Sets @min and @max to lower and upper bounds for the square of the
* Euclidean distance between the object contained in @bb and @p. For these
* bounds to make any sense the bounding box must be "tight" i.e. each of the
* 6 faces of the box must at least be touched by one point of the bounded
* object.
*/
void gts_bbox_point_distance2 (GtsBBox * bb, GtsPoint * p,
gdouble * min, gdouble * max)
{
gdouble x1, y1, z1, x2, y2, z2, x, y, z;
gdouble dmin, dmax, xd1, xd2, yd1, yd2, zd1, zd2;
gdouble mx, Mx, my, My, mz, Mz;
g_return_if_fail (bb != NULL);
g_return_if_fail (p != NULL);
g_return_if_fail (min != NULL);
g_return_if_fail (max != NULL);
x1 = bb->x1; y1 = bb->y1; z1 = bb->z1;
x2 = bb->x2; y2 = bb->y2; z2 = bb->z2;
x = p->x; y = p->y; z = p->z;
xd1 = (x1 - x)*(x1 - x);
xd2 = (x - x2)*(x - x2);
yd1 = (y1 - y)*(y1 - y);
yd2 = (y - y2)*(y - y2);
zd1 = (z1 - z)*(z1 - z);
zd2 = (z - z2)*(z - z2);
dmin = x < x1 ? xd1 : x > x2 ? xd2 : 0.0;
dmin += y < y1 ? yd1 : y > y2 ? yd2 : 0.0;
dmin += z < z1 ? zd1 : z > z2 ? zd2 : 0.0;
MINMAX (xd1, xd2, mx, Mx);
MINMAX (yd1, yd2, my, My);
MINMAX (zd1, zd2, mz, Mz);
dmax = mx + My + Mz;
dmax = MIN (dmax, Mx + my + Mz);
dmax = MIN (dmax, Mx + My + mz);
*min = dmin;
*max = dmax;
}
/**
* gts_bbox_is_stabbed:
* @bb: a #GtsBBox.
* @p: a #GtsPoint.
*
* Returns: %TRUE if the ray starting at @p and ending at (+infty,
* @p->y, @p->z) intersects with @bb, %FALSE otherwise.
*/
gboolean gts_bbox_is_stabbed (GtsBBox * bb, GtsPoint * p)
{
g_return_val_if_fail (bb != NULL, FALSE);
g_return_val_if_fail (p != NULL, FALSE);
if (p->x > bb->x2 ||
p->y < bb->y1 || p->y > bb->y2 ||
p->z < bb->z1 || p->z > bb->z2)
return FALSE;
return TRUE;
}
extern int triBoxOverlap (double boxcenter[3],
double boxhalfsize[3],
double triverts[3][3]);
/**
* gts_bbox_overlaps_triangle:
* @bb: a #GtsBBox.
* @t: a #GtsTriangle.
*
* This is a wrapper around the fast overlap test of Tomas
* Akenine-Moller (http://www.cs.lth.se/home/Tomas_Akenine_Moller/).
*
* Returns: %TRUE if @bb overlaps with @t, %FALSE otherwise.
*/
gboolean gts_bbox_overlaps_triangle (GtsBBox * bb, GtsTriangle * t)
{
double bc[3], bh[3], tv[3][3];
GtsPoint * p1, * p2, * p3;
g_return_val_if_fail (bb != NULL, FALSE);
g_return_val_if_fail (t != NULL, FALSE);
bc[0] = (bb->x2 + bb->x1)/2.;
bh[0] = (bb->x2 - bb->x1)/2.;
bc[1] = (bb->y2 + bb->y1)/2.;
bh[1] = (bb->y2 - bb->y1)/2.;
bc[2] = (bb->z2 + bb->z1)/2.;
bh[2] = (bb->z2 - bb->z1)/2.;
p1 = GTS_POINT (GTS_SEGMENT (t->e1)->v1);
p2 = GTS_POINT (GTS_SEGMENT (t->e1)->v2);
p3 = GTS_POINT (gts_triangle_vertex (t));
tv[0][0] = p1->x; tv[0][1] = p1->y; tv[0][2] = p1->z;
tv[1][0] = p2->x; tv[1][1] = p2->y; tv[1][2] = p2->z;
tv[2][0] = p3->x; tv[2][1] = p3->y; tv[2][2] = p3->z;
return triBoxOverlap (bc, bh, tv);
}
/**
* gts_bbox_overlaps_segment:
* @bb: a #GtsBBox.
* @s: a #GtsSegment.
*
* This functions uses gts_bbox_overlaps_triangle() with a degenerate
* triangle.
*
* Returns: %TRUE if @bb overlaps with @s, %FALSE otherwise.
*/
gboolean gts_bbox_overlaps_segment (GtsBBox * bb, GtsSegment * s)
{
double bc[3], bh[3], tv[3][3];
GtsPoint * p1, * p2, * p3;
g_return_val_if_fail (bb != NULL, FALSE);
g_return_val_if_fail (s != NULL, FALSE);
bc[0] = (bb->x2 + bb->x1)/2.;
bh[0] = (bb->x2 - bb->x1)/2.;
bc[1] = (bb->y2 + bb->y1)/2.;
bh[1] = (bb->y2 - bb->y1)/2.;
bc[2] = (bb->z2 + bb->z1)/2.;
bh[2] = (bb->z2 - bb->z1)/2.;
p1 = GTS_POINT (s->v1);
p2 = GTS_POINT (s->v2);
p3 = p1;
tv[0][0] = p1->x; tv[0][1] = p1->y; tv[0][2] = p1->z;
tv[1][0] = p2->x; tv[1][1] = p2->y; tv[1][2] = p2->z;
tv[2][0] = p3->x; tv[2][1] = p3->y; tv[2][2] = p3->z;
return triBoxOverlap (bc, bh, tv);
}
/**
* gts_bb_tree_new:
* @bboxes: a list of #GtsBBox.
*
* Builds a new hierarchy of bounding boxes for @bboxes. At each
* level, the GNode->data field contains a #GtsBBox bounding box of
* all the children. The tree is binary and is built by repeatedly
* cutting in two approximately equal halves the bounding boxes at
* each level until a leaf node (i.e. a bounding box given in @bboxes)
* is reached. In order to minimize the depth of the tree, the cutting
* direction is always chosen as perpendicular to the longest
* dimension of the bounding box.
*
* Returns: a new hierarchy of bounding boxes.
*/
GNode * gts_bb_tree_new (GSList * bboxes)
{
GSList * i, * positive = NULL, * negative = NULL;
GNode * node;
GtsBBox * bbox;
guint dir, np = 0, nn = 0;
gdouble * p1, * p2;
gdouble cut;
g_return_val_if_fail (bboxes != NULL, NULL);
if (bboxes->next == NULL) /* leaf node */
return g_node_new (bboxes->data);
bbox = gts_bbox_bboxes (gts_bbox_class (), bboxes);
node = g_node_new (bbox);
if (bbox->x2 - bbox->x1 > bbox->y2 - bbox->y1) {
if (bbox->z2 - bbox->z1 > bbox->x2 - bbox->x1)
dir = 2;
else
dir = 0;
}
else if (bbox->z2 - bbox->z1 > bbox->y2 - bbox->y1)
dir = 2;
else
dir = 1;
p1 = (gdouble *) &bbox->x1;
p2 = (gdouble *) &bbox->x2;
cut = (p1[dir] + p2[dir])/2.;
i = bboxes;
while (i) {
bbox = i->data;
p1 = (gdouble *) &bbox->x1;
p2 = (gdouble *) &bbox->x2;
if ((p1[dir] + p2[dir])/2. > cut) {
positive = g_slist_prepend (positive, bbox);
np++;
}
else {
negative = g_slist_prepend (negative, bbox);
nn++;
}
i = i->next;
}
if (!positive) {
GSList * last = g_slist_nth (negative, (nn - 1)/2);
positive = last->next;
last->next = NULL;
}
else if (!negative) {
GSList * last = g_slist_nth (positive, (np - 1)/2);
negative = last->next;
last->next = NULL;
}
g_node_prepend (node, gts_bb_tree_new (positive));
g_slist_free (positive);
g_node_prepend (node, gts_bb_tree_new (negative));
g_slist_free (negative);
return node;
}
static void prepend_triangle_bbox (GtsTriangle * t, GSList ** bboxes)
{
*bboxes = g_slist_prepend (*bboxes,
gts_bbox_triangle (gts_bbox_class (), t));
}
/**
* gts_bb_tree_surface:
* @s: a #GtsSurface.
*
* Returns: a new hierarchy of bounding boxes bounding the faces of @s.
*/
GNode * gts_bb_tree_surface (GtsSurface * s)
{
GSList * bboxes = NULL;
GNode * tree;
g_return_val_if_fail (s != NULL, NULL);
gts_surface_foreach_face (s, (GtsFunc) prepend_triangle_bbox, &bboxes);
tree = gts_bb_tree_new (bboxes);
g_slist_free (bboxes);
return tree;
}
/**
* gts_bb_tree_stabbed:
* @tree: a bounding box tree.
* @p: a #GtsPoint.
*
* Returns: a list of bounding boxes, leaves of @tree which are
* stabbed by the ray defined by @p (see gts_bbox_is_stabbed()).
*/
GSList * gts_bb_tree_stabbed (GNode * tree, GtsPoint * p)
{
GSList * list = NULL;
GtsBBox * bb;
GNode * i;
g_return_val_if_fail (tree != NULL, NULL);
g_return_val_if_fail (p != NULL, NULL);
bb = tree->data;
if (!gts_bbox_is_stabbed (bb, p))
return NULL;
if (tree->children == NULL) /* leaf node */
return g_slist_prepend (NULL, bb);
i = tree->children;
while (i) {
list = g_slist_concat (list, gts_bb_tree_stabbed (i, p));
i = i->next;
}
return list;
}
/**
* gts_bb_tree_overlap:
* @tree: a bounding box tree.
* @bbox: a #GtsBBox.
*
* Returns: a list of bounding boxes, leaves of @tree which overlap @bbox.
*/
GSList * gts_bb_tree_overlap (GNode * tree, GtsBBox * bbox)
{
GSList * list = NULL;
GtsBBox * bb;
GNode * i;
g_return_val_if_fail (tree != NULL, NULL);
g_return_val_if_fail (bbox != NULL, NULL);
bb = tree->data;
if (!gts_bboxes_are_overlapping (bbox, bb))
return NULL;
if (tree->children == NULL) /* leaf node */
return g_slist_prepend (NULL, bb);
i = tree->children;
while (i) {
list = g_slist_concat (list, gts_bb_tree_overlap (i, bbox));
i = i->next;
}
return list;
}
/**
* gts_bb_tree_is_overlapping:
* @tree: a bounding box tree.
* @bbox: a #GtsBBox.
*
* Returns: %TRUE if any leaf of @tree overlaps @bbox, %FALSE otherwise.
*/
gboolean gts_bb_tree_is_overlapping (GNode * tree, GtsBBox * bbox)
{
GtsBBox * bb;
GNode * i;
g_return_val_if_fail (tree != NULL, FALSE);
g_return_val_if_fail (bbox != NULL, FALSE);
bb = tree->data;
if (!gts_bboxes_are_overlapping (bbox, bb))
return FALSE;
if (tree->children == NULL) /* leaf node */
return TRUE;
i = tree->children;
while (i) {
if (gts_bb_tree_is_overlapping (i, bbox))
return TRUE;
i = i->next;
}
return FALSE;
}
/**
* gts_bb_tree_traverse_overlapping:
* @tree1: a bounding box tree.
* @tree2: a bounding box tree.
* @func: a #GtsBBTreeTraverseFunc.
* @data: user data to be passed to @func.
*
* Calls @func for each overlapping pair of leaves of @tree1 and @tree2.
*/
void gts_bb_tree_traverse_overlapping (GNode * tree1, GNode * tree2,
GtsBBTreeTraverseFunc func,
gpointer data)
{
GtsBBox * bb1, * bb2;
g_return_if_fail (tree1 != NULL && tree2 != NULL);
bb1 = tree1->data; bb2 = tree2->data;
if (!gts_bboxes_are_overlapping (bb1, bb2))
return;
if (tree1->children == NULL && tree2->children == NULL)
(*func) (tree1->data, tree2->data, data);
else if (tree2->children == NULL ||
(tree1->children != NULL &&
bbox_volume (bb1) > bbox_volume (bb2))) {
GNode * i = tree1->children;
while (i) {
gts_bb_tree_traverse_overlapping (i, tree2, func, data);
i = i->next;
}
}
else {
GNode * i = tree2->children;
while (i) {
gts_bb_tree_traverse_overlapping (tree1, i, func, data);
i = i->next;
}
}
}
/**
* gts_bb_tree_draw:
* @tree: a bounding box tree.
* @depth: a specified depth.
* @fptr: a file pointer.
*
* Write in @fptr an OOGL (Geomview) description of @tree for the
* depth specified by @depth.
*/
void gts_bb_tree_draw (GNode * tree, guint depth, FILE * fptr)
{
guint d;
g_return_if_fail (tree != NULL);
g_return_if_fail (fptr != NULL);
d = g_node_depth (tree);
if (d == 1)
fprintf (fptr, "{ LIST");
if (d == depth)
gts_bbox_draw (tree->data, fptr);
else if (d < depth) {
GNode * i = tree->children;
while (i) {
gts_bb_tree_draw (i, depth, fptr);
i = i->next;
}
}
if (d == 1)
fprintf (fptr, "}\n");
}
static void bb_tree_free (GNode * tree, gboolean free_leaves)
{
GNode * i;
g_return_if_fail (tree != NULL);
if (!free_leaves && tree->children == NULL) /* leaf node */
return;
gts_object_destroy (tree->data);
i = tree->children;
while (i) {
bb_tree_free (i, free_leaves);
i = i->next;
}
}
/**
* gts_bb_tree_destroy:
* @tree: a bounding box tree.
* @free_leaves: if %TRUE the bounding boxes given by the user are freed.
*
* Destroys all the bounding boxes created by @tree and destroys the
* tree itself. If @free_leaves is set to %TRUE, destroys boxes given
* by the user when creating the tree (i.e. leaves of the tree).
*/
void gts_bb_tree_destroy (GNode * tree, gboolean free_leaves)
{
g_return_if_fail (tree != NULL);
bb_tree_free (tree, free_leaves);
g_node_destroy (tree);
}
static gdouble bb_tree_min_max (GNode * tree,
GtsPoint * p,
gdouble min_max,
GSList ** list)
{
GNode * tree1, * tree2;
gdouble min1, max1, min2, max2;
if (tree->children == NULL) {
*list = g_slist_prepend (*list, tree->data);
return min_max;
}
tree1 = tree->children;
gts_bbox_point_distance2 (tree1->data, p, &min1, &max1);
if (max1 < min_max)
min_max = max1;
tree2 = tree1->next;
gts_bbox_point_distance2 (tree2->data, p, &min2, &max2);
if (max2 < min_max)
min_max = max2;
if (min1 < min2) {
if (min1 <= min_max) {
min_max = bb_tree_min_max (tree1, p, min_max, list);
if (min2 <= min_max)
min_max = bb_tree_min_max (tree2, p, min_max, list);
}
}
else {
if (min2 <= min_max) {
min_max = bb_tree_min_max (tree2, p, min_max, list);
if (min1 <= min_max)
min_max = bb_tree_min_max (tree1, p, min_max, list);
}
}
return min_max;
}
/**
* gts_bb_tree_point_closest_bboxes:
* @tree: a bounding box tree.
* @p: a #GtsPoint.
*
* Returns: a list of #GtsBBox. One of the bounding boxes is assured to contain
* the object of @tree closest to @p.
*/
GSList * gts_bb_tree_point_closest_bboxes (GNode * tree,
GtsPoint * p)
{
gdouble min, min_max;
GSList * list = NULL, * i, * prev = NULL;
g_return_val_if_fail (tree != NULL, NULL);
g_return_val_if_fail (p != NULL, NULL);
gts_bbox_point_distance2 (tree->data, p, &min, &min_max);
min_max = bb_tree_min_max (tree, p, min_max, &list);
i = list;
while (i) {
GSList * next = i->next;
gdouble min, max;
gts_bbox_point_distance2 (i->data, p, &min, &max);
if (min > min_max) {
if (prev == NULL)
list = next;
else
prev->next = next;
g_slist_free_1 (i);
}
else
prev = i;
i = next;
}
return list;
}
/**
* gts_bb_tree_point_distance:
* @tree: a bounding box tree.
* @p: a #GtsPoint.
* @distance: a #GtsBBoxDistFunc.
* @bbox: if not %NULL is set to the bounding box containing the closest
* object.
*
* Returns: the distance as evaluated by @distance between @p and the closest
* object in @tree.
*/
gdouble gts_bb_tree_point_distance (GNode * tree,
GtsPoint * p,
GtsBBoxDistFunc distance,
GtsBBox ** bbox)
{
GSList * list, * i;
gdouble dmin = G_MAXDOUBLE;
g_return_val_if_fail (tree != NULL, dmin);
g_return_val_if_fail (p != NULL, dmin);
g_return_val_if_fail (distance != NULL, dmin);
i = list = gts_bb_tree_point_closest_bboxes (tree, p);
while (i) {
gdouble d = (*distance) (p, GTS_BBOX (i->data)->bounded);
if (fabs (d) < fabs (dmin)) {
dmin = d;
if (bbox)
*bbox = i->data;
}
i = i->next;
}
g_slist_free (list);
return dmin;
}
/**
* gts_bb_tree_point_closest:
* @tree: a bounding box tree.
* @p: a #GtsPoint.
* @closest: a #GtsBBoxClosestFunc.
* @distance: if not %NULL is set to the distance between @p and the
* new #GtsPoint.
*
* Returns: a new #GtsPoint, closest point to @p and belonging to an object of
* @tree.
*/
GtsPoint * gts_bb_tree_point_closest (GNode * tree,
GtsPoint * p,
GtsBBoxClosestFunc closest,
gdouble * distance)
{
GSList * list, * i;
gdouble dmin = G_MAXDOUBLE;
GtsPoint * np = NULL;
g_return_val_if_fail (tree != NULL, NULL);
g_return_val_if_fail (p != NULL, NULL);
g_return_val_if_fail (closest != NULL, NULL);
i = list = gts_bb_tree_point_closest_bboxes (tree, p);
while (i) {
GtsPoint * tp = (*closest) (p, GTS_BBOX (i->data)->bounded);
gdouble d = gts_point_distance2 (tp, p);
if (d < dmin) {
if (np)
gts_object_destroy (GTS_OBJECT (np));
np = tp;
dmin = d;
}
else
gts_object_destroy (GTS_OBJECT (tp));
i = i->next;
}
g_slist_free (list);
if (distance)
*distance = dmin;
return np;
}
/**
* gts_bb_tree_triangle_distance:
* @tree: a bounding box tree.
* @t: a #GtsTriangle.
* @distance: a #GtsBBoxDistFunc.
* @delta: spatial scale of the sampling to be used.
* @range: a #GtsRange to be filled with the results.
*
* Given a triangle @t, points are sampled regularly on its surface
* using @delta as increment. The distance from each of these points
* to the closest object of @tree is computed using @distance and the
* gts_bb_tree_point_distance() function. The fields of @range are
* filled with the number of points sampled, the minimum, average and
* maximum value and the standard deviation.
*/
void gts_bb_tree_triangle_distance (GNode * tree,
GtsTriangle * t,
GtsBBoxDistFunc distance,
gdouble delta,
GtsRange * range)
{
GtsPoint * p1, * p2, * p3, * p;
GtsVector p1p2, p1p3;
gdouble l1, t1, dt1;
guint i, n1;
g_return_if_fail (tree != NULL);
g_return_if_fail (t != NULL);
g_return_if_fail (distance != NULL);
g_return_if_fail (delta > 0.);
g_return_if_fail (range != NULL);
gts_triangle_vertices (t,
(GtsVertex **) &p1,
(GtsVertex **) &p2,
(GtsVertex **) &p3);
gts_vector_init (p1p2, p1, p2);
gts_vector_init (p1p3, p1, p3);
gts_range_init (range);
p = GTS_POINT (gts_object_new (GTS_OBJECT_CLASS (gts_point_class ())));
l1 = sqrt (gts_vector_scalar (p1p2, p1p2));
n1 = l1/delta + 1;
dt1 = 1.0/(gdouble) n1;
t1 = 0.0;
for (i = 0; i <= n1; i++, t1 += dt1) {
gdouble t2 = 1. - t1;
gdouble x = t2*p1p3[0];
gdouble y = t2*p1p3[1];
gdouble z = t2*p1p3[2];
gdouble l2 = sqrt (x*x + y*y + z*z);
guint j, n2 = (guint) (l2/delta + 1);
gdouble dt2 = t2/(gdouble) n2;
x = t2*p1->x + t1*p2->x;
y = t2*p1->y + t1*p2->y;
z = t2*p1->z + t1*p2->z;
t2 = 0.0;
for (j = 0; j <= n2; j++, t2 += dt2) {
p->x = x + t2*p1p3[0];
p->y = y + t2*p1p3[1];
p->z = z + t2*p1p3[2];
gts_range_add_value (range,
gts_bb_tree_point_distance (tree, p, distance, NULL));
}
}
gts_object_destroy (GTS_OBJECT (p));
gts_range_update (range);
}
/**
* gts_bb_tree_segment_distance:
* @tree: a bounding box tree.
* @s: a #GtsSegment.
* @distance: a #GtsBBoxDistFunc.
* @delta: spatial scale of the sampling to be used.
* @range: a #GtsRange to be filled with the results.
*
* Given a segment @s, points are sampled regularly on its length
* using @delta as increment. The distance from each of these points
* to the closest object of @tree is computed using @distance and the
* gts_bb_tree_point_distance() function. The fields of @range are
* filled with the number of points sampled, the minimum, average and
* maximum value and the standard deviation.
*/
void gts_bb_tree_segment_distance (GNode * tree,
GtsSegment * s,
gdouble (*distance) (GtsPoint *,
gpointer),
gdouble delta,
GtsRange * range)
{
GtsPoint * p1, * p2, * p;
GtsVector p1p2;
gdouble l, t, dt;
guint i, n;
g_return_if_fail (tree != NULL);
g_return_if_fail (s != NULL);
g_return_if_fail (distance != NULL);
g_return_if_fail (delta > 0.);
g_return_if_fail (range != NULL);
p1 = GTS_POINT (s->v1);
p2 = GTS_POINT (s->v2);
gts_vector_init (p1p2, p1, p2);
gts_range_init (range);
p = GTS_POINT (gts_object_new (GTS_OBJECT_CLASS (gts_point_class())));
l = sqrt (gts_vector_scalar (p1p2, p1p2));
n = (guint) (l/delta + 1);
dt = 1.0/(gdouble) n;
t = 0.0;
for (i = 0; i <= n; i++, t += dt) {
p->x = p1->x + t*p1p2[0];
p->y = p1->y + t*p1p2[1];
p->z = p1->z + t*p1p2[2];
gts_range_add_value (range,
gts_bb_tree_point_distance (tree, p, distance, NULL));
}
gts_object_destroy (GTS_OBJECT (p));
gts_range_update (range);
}
static void surface_distance_foreach_triangle (GtsTriangle * t,
gpointer * data)
{
gdouble * delta = data[1];
GtsRange * range = data[2];
gdouble * total_area = data[3], area;
GtsRange range_triangle;
gts_bb_tree_triangle_distance (data[0], t, data[4], *delta, &range_triangle);
if (range_triangle.min < range->min)
range->min = range_triangle.min;
if (range_triangle.max > range->max)
range->max = range_triangle.max;
range->n += range_triangle.n;
area = gts_triangle_area (t);
*total_area += area;
range->sum += area*range_triangle.mean;
range->sum2 += area*range_triangle.mean*range_triangle.mean;
}
/**
* gts_bb_tree_surface_distance:
* @tree: a bounding box tree.
* @s: a #GtsSurface.
* @distance: a #GtsBBoxDistFunc.
* @delta: a sampling increment defined as the percentage of the diagonal
* of the root bounding box of @tree.
* @range: a #GtsRange to be filled with the results.
*
* Calls gts_bb_tree_triangle_distance() for each face of @s. The
* fields of @range are filled with the minimum, maximum and average
* distance. The average distance is defined as the sum of the average
* distances for each triangle weighthed by their area and divided by
* the total area of the surface. The standard deviation is defined
* accordingly. The @n field of @range is filled with the number of
* sampled points used.
*/
void gts_bb_tree_surface_distance (GNode * tree,
GtsSurface * s,
GtsBBoxDistFunc distance,
gdouble delta,
GtsRange * range)
{
gpointer data[5];
gdouble total_area = 0.;
g_return_if_fail (tree != NULL);
g_return_if_fail (s != NULL);
g_return_if_fail (delta > 0. && delta < 1.);
g_return_if_fail (range != NULL);
gts_range_init (range);
delta *= sqrt (gts_bbox_diagonal2 (tree->data));
data[0] = tree;
data[1] = δ
data[2] = range;
data[3] = &total_area;
data[4] = distance;
gts_surface_foreach_face (s,
(GtsFunc) surface_distance_foreach_triangle,
data);
if (total_area > 0.) {
if (range->sum2 - range->sum*range->sum/total_area >= 0.)
range->stddev = sqrt ((range->sum2 - range->sum*range->sum/total_area)
/total_area);
else
range->stddev = 0.;
range->mean = range->sum/total_area;
}
else
range->min = range->max = range->mean = range->stddev = 0.;
}
static void surface_distance_foreach_boundary (GtsEdge * e,
gpointer * data)
{
gdouble * delta = data[1];
GtsRange * range = data[2];
gdouble * total_length = data[3], length;
GtsRange range_edge;
if (gts_edge_is_boundary (e, NULL)) {
GtsSegment * s = GTS_SEGMENT (e);
gts_bb_tree_segment_distance (data[0], s, data[4], *delta, &range_edge);
if (range_edge.min < range->min)
range->min = range_edge.min;
if (range_edge.max > range->max)
range->max = range_edge.max;
range->n += range_edge.n;
length = gts_point_distance (GTS_POINT (s->v1), GTS_POINT (s->v2));
*total_length += length;
range->sum += length*range_edge.mean;
range->sum2 += length*range_edge.mean*range_edge.mean;
}
}
/**
* gts_bb_tree_surface_boundary_distance:
* @tree: a bounding box tree.
* @s: a #GtsSurface.
* @distance: a #GtsBBoxDistFunc.
* @delta: a sampling increment defined as the percentage of the diagonal
* of the root bounding box of @tree.
* @range: a #GtsRange to be filled with the results.
*
* Calls gts_bb_tree_segment_distance() for each edge boundary of @s.
* The fields of @range are filled with the minimum, maximum and
* average distance. The average distance is defined as the sum of the
* average distances for each boundary edge weighthed by their length
* and divided by the total length of the boundaries. The standard
* deviation is defined accordingly. The @n field of @range is filled
* with the number of sampled points used.
*/
void gts_bb_tree_surface_boundary_distance (GNode * tree,
GtsSurface * s,
gdouble (*distance) (GtsPoint *,
gpointer),
gdouble delta,
GtsRange * range)
{
gpointer data[5];
gdouble total_length = 0.;
g_return_if_fail (tree != NULL);
g_return_if_fail (s != NULL);
g_return_if_fail (delta > 0. && delta < 1.);
g_return_if_fail (range != NULL);
gts_range_init (range);
delta *= sqrt (gts_bbox_diagonal2 (tree->data));
data[0] = tree;
data[1] = δ
data[2] = range;
data[3] = &total_length;
data[4] = distance;
gts_surface_foreach_edge (s,
(GtsFunc) surface_distance_foreach_boundary,
data);
if (total_length > 0.) {
if (range->sum2 - range->sum*range->sum/total_length >= 0.)
range->stddev = sqrt ((range->sum2 -
range->sum*range->sum/total_length)
/total_length);
else
range->stddev = 0.;
range->mean = range->sum/total_length;
}
else
range->min = range->max = range->mean = range->stddev = 0.;
}
syntax highlighted by Code2HTML, v. 0.9.1