}
if (
- (NULL == CU_add_test(pSuite, "test_gbox_from_spherical_coordinates()", test_gbox_from_spherical_coordinates)) ||
+ (NULL == CU_add_test(pSuite, "test_signum()", test_signum)) ||
+ (NULL == CU_add_test(pSuite, "test_gbox_from_spherical_coordinates()", test_gbox_from_spherical_coordinates)) /* ||
(NULL == CU_add_test(pSuite, "test_gserialized_get_gbox_geocentric()", test_gserialized_get_gbox_geocentric)) ||
- (NULL == CU_add_test(pSuite, "test_gbox_calculation()", test_gbox_calculation))
+ (NULL == CU_add_test(pSuite, "test_gbox_calculation()", test_gbox_calculation)) */
)
{
CU_cleanup_registry();
}
+void test_signum(void)
+{
+ CU_ASSERT_EQUAL(signum(-5.0),-1);
+ CU_ASSERT_EQUAL(signum(5.0),1);
+}
void test_gbox_from_spherical_coordinates(void)
{
pa = pointArray_construct((uchar*)ll, 0, 0, 2);
rv = ptarray_calculate_gbox_geodetic(pa, box);
- printf("\n%s\n", gbox_to_string(box));
- printf("%s\n", "(0.9595879536661512, -0.052998127820754914, -0.09150161866340238) (0.9949586928172913, 0.2661172917357788, 0.17078275740560678)");
- CU_ASSERT_DOUBLE_EQUAL(box->xmin, 0.9949586928172913, 0.0001);
+// printf("\n%s\n", gbox_to_string(box));
+// printf("%s\n", "(0.95958795, -0.05299812, -0.09150161) (0.99495869, 0.26611729, 0.17078275)");
+ CU_ASSERT_DOUBLE_EQUAL(box->xmin, 0.95958795, 0.0001);
lwfree(pa);
pa = pointArray_construct((uchar*)ll, 0, 0, 2);
rv = ptarray_calculate_gbox_geodetic(pa, box);
- printf("\n%s\n", gbox_to_string(box));
- printf("%s\n", "(0.32139380484326974, -0.3491712110387923, 0.7933533402912352) (0.49866816905678146, 0.3830222215594891, 0.9014764513830344)");
- CU_ASSERT_DOUBLE_EQUAL(box->xmin, 0.9014764513830344, 0.0001);
+// printf("\n%s\n", gbox_to_string(box));
+// printf("%s\n", "(0.32139380, -0.34917121, 0.79335334) (0.49866816, 0.38302222, 0.90147645)");
+ CU_ASSERT_DOUBLE_EQUAL(box->xmin, 0.32139380, 0.0001);
lwfree(pa);
pa = pointArray_construct((uchar*)ll, 0, 0, 2);
rv = ptarray_calculate_gbox_geodetic(pa, box);
- printf("\n%s\n", gbox_to_string(box));
- printf("%s\n", "(-0.2985076686105858, -0.46856318207700426, 0.8314696123025452) (-0.19629681530223578, -0.2965721369531419, 0.9346189212108039)");
- CU_ASSERT_DOUBLE_EQUAL(box->zmin, 0.8314696123025452, 0.0001);
+// printf("\n%s\n", gbox_to_string(box));
+// printf("%s\n", "(-0.29850766, -0.46856318, 0.83146961) (-0.19629681, -0.29657213, 0.93461892)");
+ CU_ASSERT_DOUBLE_EQUAL(box->zmin, 0.83146961, 0.0001);
lwfree(pa);
lwfree(box);
/* Test functions */
+void test_signum(void);
void test_gbox_from_spherical_coordinates(void);
void test_gserialized_get_gbox_geocentric(void);
void test_gbox_calculation(void);
}
/* Add the libgeom suite to the registry */
- if (NULL == register_libgeom_suite())
- {
- CU_cleanup_registry();
- return CU_get_error();
- }
+// if (NULL == register_libgeom_suite())
+// {
+// CU_cleanup_registry();
+// return CU_get_error();
+// }
/* Add the geodetic suite to the registry */
if (NULL == register_geodetic_suite())
return g;
}
-int gbox_merge_point3d(GBOX *gbox, POINT3D *p)
+int gbox_merge_point3d(POINT3D p, GBOX *gbox)
{
- if( gbox->xmin > p->x ) gbox->xmin = p->x;
- if( gbox->ymin > p->y ) gbox->ymin = p->y;
- if( gbox->zmin > p->z ) gbox->zmin = p->z;
- if( gbox->xmax < p->x ) gbox->xmax = p->x;
- if( gbox->ymax < p->y ) gbox->ymax = p->y;
- if( gbox->zmax < p->z ) gbox->zmax = p->z;
+ if( gbox->xmin > p.x ) gbox->xmin = p.x;
+ if( gbox->ymin > p.y ) gbox->ymin = p.y;
+ if( gbox->zmin > p.z ) gbox->zmin = p.z;
+ if( gbox->xmax < p.x ) gbox->xmax = p.x;
+ if( gbox->ymax < p.y ) gbox->ymax = p.y;
+ if( gbox->zmax < p.z ) gbox->zmax = p.z;
return G_SUCCESS;
}
-int gbox_merge(GBOX *master_box, GBOX *new_box)
+
+int gbox_merge(GBOX new_box, GBOX *merge_box)
{
- assert(master_box);
- assert(new_box);
+ assert(merge_box);
- if( master_box->flags != new_box->flags )
+ if( merge_box->flags != new_box.flags )
return G_FAILURE;
- if( new_box->xmin < master_box->xmin) master_box->xmin = new_box->xmin;
- if( new_box->ymin < master_box->ymin) master_box->ymin = new_box->ymin;
- if( new_box->xmax > master_box->xmax) master_box->xmax = new_box->xmax;
- if( new_box->ymax > master_box->ymax) master_box->ymax = new_box->ymax;
+ if( new_box.xmin < merge_box->xmin) merge_box->xmin = new_box.xmin;
+ if( new_box.ymin < merge_box->ymin) merge_box->ymin = new_box.ymin;
+ if( new_box.xmax > merge_box->xmax) merge_box->xmax = new_box.xmax;
+ if( new_box.ymax > merge_box->ymax) merge_box->ymax = new_box.ymax;
- if( FLAGS_GET_Z(master_box->flags) || FLAGS_GET_GEODETIC(master_box->flags) )
+ if( FLAGS_GET_Z(merge_box->flags) || FLAGS_GET_GEODETIC(merge_box->flags) )
{
- if( new_box->zmin < master_box->zmin) master_box->zmin = new_box->zmin;
- if( new_box->zmax > master_box->zmax) master_box->zmax = new_box->zmax;
+ if( new_box.zmin < merge_box->zmin) merge_box->zmin = new_box.zmin;
+ if( new_box.zmax > merge_box->zmax) merge_box->zmax = new_box.zmax;
}
- if( FLAGS_GET_M(master_box->flags) )
+ if( FLAGS_GET_M(merge_box->flags) )
{
- if( new_box->mmin < master_box->mmin) master_box->mmin = new_box->mmin;
- if( new_box->mmax > master_box->mmax) master_box->mmax = new_box->mmax;
+ if( new_box.mmin < merge_box->mmin) merge_box->mmin = new_box.mmin;
+ if( new_box.mmax > merge_box->mmax) merge_box->mmax = new_box.mmax;
}
return G_SUCCESS;
return LW_TRUE;
}
+#if 0
+GBOX* gbox_from_string(char *str)
+{
+ int ndims = 0;
+ char *ptr = str;
+ char *nextptr;
+ double d;
+ char *gbox_start = strstr(str, "GBOX((");
+ GBOX *gbox = gbox_new(0);
+ if( ! gbox_start ) return NULL; /* No header found */
+ ptr += 5;
+ do {
+ ptr++;
+ d = strtod(ptr, &nextptr);
+ if( ptr == nextptr ) return NULL; /* No double found */
+ gbox->xmin = d;
+ ndims++;
+ ptr = nextptr;
+
+
+}
+#endif
+
char* gbox_to_string(GBOX *gbox)
{
static int sz = 128;
return copy;
}
-void gbox_duplicate(GBOX *original, GBOX *duplicate)
+void gbox_duplicate(GBOX original, GBOX *duplicate)
{
- assert(original);
assert(duplicate);
- if( original->flags != duplicate->flags )
+ if( original.flags != duplicate->flags )
lwerror("gbox_duplicate: geometries have inconsistent dimensionality");
- duplicate->xmin = original->xmin;
- duplicate->ymin = original->ymin;
- duplicate->xmax = original->xmax;
- duplicate->ymax = original->ymax;
+ duplicate->xmin = original.xmin;
+ duplicate->ymin = original.ymin;
+ duplicate->xmax = original.xmax;
+ duplicate->ymax = original.ymax;
- if( FLAGS_GET_GEODETIC(original->flags) || FLAGS_GET_Z(original->flags) )
+ if( FLAGS_GET_GEODETIC(original.flags) || FLAGS_GET_Z(original.flags) )
{
- duplicate->zmin = original->zmin;
- duplicate->zmax = original->zmax;
+ duplicate->zmin = original.zmin;
+ duplicate->zmax = original.zmax;
}
- if( FLAGS_GET_M(original->flags) )
+ if( FLAGS_GET_M(original.flags) )
{
- duplicate->mmin = original->mmin;
- duplicate->mmax = original->mmax;
+ duplicate->mmin = original.mmin;
+ duplicate->mmax = original.mmax;
}
return;
}
if (lwcircle_calculate_gbox(p1, p2, p3, &tmp) == G_FAILURE)
continue;
- gbox_merge(gbox, &tmp);
+ gbox_merge(tmp, gbox);
}
return G_SUCCESS;
{
if( first )
{
- gbox_duplicate(gbox, &subbox);
+ gbox_duplicate(subbox, gbox);
first = LW_FALSE;
}
else
{
- gbox_merge(gbox, &subbox);
+ gbox_merge(subbox, gbox);
}
result = G_SUCCESS;
}
{
if ( first )
{
- gbox_duplicate(&subbox, gbox);
+ gbox_duplicate(subbox, gbox);
first = LW_FALSE;
}
else
{
- gbox_merge(gbox, &subbox);
+ gbox_merge(subbox, gbox);
}
result = G_SUCCESS;
}
/**
* Update the merged #GBOX to be large enough to include itself and the new box.
*/
-extern int gbox_merge(GBOX *merged_box, GBOX *new_box);
+extern int gbox_merge(GBOX new_box, GBOX *merged_box);
/**
* Update the #GBOX to be large enough to include itself and the new point.
*/
-extern int gbox_merge_point3d(GBOX *gbox, POINT3D *p);
+extern int gbox_merge_point3d(POINT3D p, GBOX *gbox);
/**
* Allocate a string representation of the #GBOX, based on dimensionality of flags.
/**
* Copy the values of original #GBOX into duplicate.
*/
-extern void gbox_duplicate(GBOX *original, GBOX *duplicate);
+extern void gbox_duplicate(GBOX original, GBOX *duplicate);
/**
* Return the number of bytes necessary to hold a #GBOX of this dimension in
#include "lwgeodetic.h"
+/**
+* Check to see if this geocentric gbox is wrapped around the north pole.
+* Only makes sense if this gbox originated from a polygon, as it's assuming
+* the box is generated from external edges and there's an "interior" which
+* contains the pole.
+*/
+static int gbox_contains_north_pole(GBOX gbox)
+{
+ if( gbox.zmin > 0.0 &&
+ gbox.xmin < 0.0 && gbox.xmax > 0.0 &&
+ gbox.ymin < 0.0 && gbox.ymax > 0.0 )
+ {
+ return LW_TRUE;
+ }
+ return LW_FALSE;
+}
+
+/**
+* Check to see if this geocentric gbox is wrapped around the south pole.
+* Only makes sense if this gbox originated from a polygon, as it's assuming
+* the box is generated from external edges and there's an "interior" which
+* contains the pole.
+*/
+static int gbox_contains_south_pole(GBOX gbox)
+{
+ if( gbox.zmax < 0.0 &&
+ gbox.xmin < 0.0 && gbox.xmax > 0.0 &&
+ gbox.ymin < 0.0 && gbox.ymax > 0.0 )
+ {
+ return LW_TRUE;
+ }
+ return LW_FALSE;
+}
+
+
/**
* Convert spherical coordinates to cartesion coordinates on unit sphere
*/
-void inline geog2cart(GEOGRAPHIC_POINT g, POINT3D *p)
+static void inline geog2cart(GEOGRAPHIC_POINT g, POINT3D *p)
{
p->x = cos(g.lat) * cos(g.lon);
p->y = cos(g.lat) * sin(g.lon);
/**
* Convert cartesion coordinates to spherical coordinates on unit sphere
*/
-void inline cart2geog(POINT3D p, GEOGRAPHIC_POINT *g)
+static void inline cart2geog(POINT3D p, GEOGRAPHIC_POINT *g)
{
g->lon = atan2(p.y, p.x);
g->lat = asin(p.z);
/**
* Calculate the dot product of two unit vectors
*/
-double inline dot_product(POINT3D p1, POINT3D p2)
+static double inline dot_product(POINT3D p1, POINT3D p2)
{
return (p1.x*p2.x) + (p1.y*p2.y) + (p1.z*p2.z);
}
-
/**
* Normalize to a unit vector.
*/
-void inline normalize(POINT3D *p)
+static void inline normalize(POINT3D *p)
{
double d = sqrt(p->x*p->x + p->y*p->y + p->z*p->z);
if(FP_IS_ZERO(d))
return;
}
-void inline unit_normal(POINT3D a, POINT3D b, POINT3D *n)
+static void inline unit_normal(POINT3D a, POINT3D b, POINT3D *n)
{
n->x = a.y * b.z - a.z * b.y;
n->y = a.z * b.x - a.x * b.z;
return G_SUCCESS;
}
-
-
+/**
+* The magic function, given an edge in spherical coordinates, calculate a
+* 3D bounding box that fully contains it, taking into account the curvature
+* of the sphere on which it is inscribed. Note special case testing
+* for edges over poles and fully around the globe. An edge is assumed
+* to follow the shortest great circle route between the end points.
+*/
int edge_calculate_gbox(GEOGRAPHIC_EDGE e, GBOX *gbox)
{
double deltaLongitude;
{
geog2cart(vT1, &p);
LWDEBUGF(4, "p == POINT(%.8g %.8g %.8g)", p.x, p.y, p.z);
- gbox_merge_point3d(gbox, &p);
+ gbox_merge_point3d(p, gbox);
LWDEBUG(4, "edge contained vT1");
LWDEBUGF(4, "gbox: %s", gbox_to_string(gbox));
}
{
geog2cart(vT2, &p);
LWDEBUGF(4, "p == POINT(%.8g %.8g %.8g)", p.x, p.y, p.z);
- gbox_merge_point3d(gbox, &p);
+ gbox_merge_point3d(p, gbox);
LWDEBUG(4, "edge contained vT2");
LWDEBUGF(4, "gbox: %s", gbox_to_string(gbox));
}
{
geog2cart(vT1, &p);
LWDEBUGF(4, "p == POINT(%.8g %.8g %.8g)", p.x, p.y, p.z);
- gbox_merge_point3d(gbox, &p);
+ gbox_merge_point3d(p, gbox);
LWDEBUG(4, "edge contained vT1");
LWDEBUGF(4, "gbox: %s", gbox_to_string(gbox));
}
{
geog2cart(vT2, &p);
LWDEBUGF(4, "p == POINT(%.8g %.8g %.8g)", p.x, p.y, p.z);
- gbox_merge_point3d(gbox, &p);
+ gbox_merge_point3d(p, gbox);
LWDEBUG(4, "edge contained vT2");
LWDEBUGF(4, "gbox: %s", gbox_to_string(gbox));
}
{
geog2cart(vT1, &p);
LWDEBUGF(4, "p == POINT(%.8g %.8g %.8g)", p.x, p.y, p.z);
- gbox_merge_point3d(gbox, &p);
+ gbox_merge_point3d(p, gbox);
LWDEBUG(4, "edge contained vT1");
LWDEBUGF(4, "gbox: %s", gbox_to_string(gbox));
}
{
geog2cart(vT2, &p);
LWDEBUGF(4, "p == POINT(%.8g %.8g %.8g)", p.x, p.y, p.z);
- gbox_merge_point3d(gbox, &p);
+ gbox_merge_point3d(p, gbox);
LWDEBUG(4, "edge contained vT2");
LWDEBUGF(4, "gbox: %s", gbox_to_string(gbox));
}
LWDEBUGF(4, "gbox after: %s", gbox_to_string(gbox));
}
- LWDEBUG(4, "leaving function");
+ LWDEBUGF(4, "final gbox: %s", gbox_to_string(gbox));
return G_SUCCESS;
}
edge.end.lat = deg2rad(end_pt->y);
edge_calculate_gbox(edge, &edge_gbox);
+
+ LWDEBUGF(4, "edge_gbox: %s", gbox_to_string(&edge_gbox));
- /* Expand the box where necessary */
- if( ! first )
+ /* Initialize the box */
+ if( first )
{
- if( edge_gbox.xmin < gbox->xmin ) gbox->xmin = edge_gbox.xmin;
- if( edge_gbox.ymin < gbox->ymin ) gbox->ymin = edge_gbox.ymin;
- if( edge_gbox.zmin < gbox->zmin ) gbox->zmin = edge_gbox.zmin;
- if( edge_gbox.xmax > gbox->xmax ) gbox->xmax = edge_gbox.xmax;
- if( edge_gbox.ymax > gbox->ymax ) gbox->ymax = edge_gbox.ymax;
- if( edge_gbox.zmax > gbox->zmax ) gbox->zmax = edge_gbox.zmax;
+ gbox_duplicate(edge_gbox,gbox);
+ LWDEBUGF(4, "gbox_duplicate: %s", gbox_to_string(gbox));
+ first = LW_FALSE;
}
- /* Initialize the box */
+ /* Expand the box where necessary */
else
{
- gbox->xmin = edge_gbox.xmin;
- gbox->ymin = edge_gbox.ymin;
- gbox->zmin = edge_gbox.zmin;
- gbox->xmax = edge_gbox.xmax;
- gbox->ymax = edge_gbox.ymax;
- gbox->zmax = edge_gbox.zmax;
- first = LW_FALSE;
+ gbox_merge(edge_gbox, gbox);
+ LWDEBUGF(4, "gbox_merge: %s", gbox_to_string(gbox));
}
}
return G_FAILURE;
if( first )
{
- gbox_duplicate(gbox, &ringbox);
+ gbox_duplicate(ringbox, gbox);
first = LW_FALSE;
}
else
{
- gbox_merge(gbox, &ringbox);
+ gbox_merge(ringbox, gbox);
}
}
+
+ /* If the box wraps the south pole, set zmin to the absolute min. */
+ if( gbox_contains_south_pole(*gbox) )
+ gbox->zmin = -1.0;
+ /* If the box wraps the north pole, set zmax to the absolute max. */
+ if( gbox_contains_north_pole(*gbox) )
+ gbox->zmax = 1.0;
+
return G_SUCCESS;
}
{
if( first )
{
- gbox_duplicate(gbox, &subbox);
+ gbox_duplicate(subbox, gbox);
first = LW_FALSE;
}
else
{
- gbox_merge(gbox, &subbox);
+ gbox_merge(subbox, gbox);
}
result = G_SUCCESS;
}
/*
** Prototypes for internal functions.
*/
-void inline geog2cart(GEOGRAPHIC_POINT g, POINT3D *p);
-void inline cart2geog(POINT3D p, GEOGRAPHIC_POINT *g);
-double inline dot_product(POINT3D p1, POINT3D p2);
-void inline unit_normal(POINT3D a, POINT3D b, POINT3D *n);
-void normalize(POINT3D *p);
void robust_cross_product(GEOGRAPHIC_POINT p, GEOGRAPHIC_POINT q, POINT3D *a);
void x_to_z(POINT3D *p);
void y_to_z(POINT3D *p);
projects / postgis / commitdiff