self-hosted/postgis
self-hosted/postgis
1
0
Fork 0
mirror of https://git.osgeo.org/gitea/postgis/postgis synced 2024-10-24 17:12:35 +00:00
Code Issues Packages Projects Releases Wiki Activity

Test re-org and first cut at edge intersection.

Browse source 
git-svn-id: http://svn.osgeo.org/postgis/trunk@4537 b70326c6-7e19-0410-871a-916f4a2858ee
This commit is contained in:
Paul Ramsey 2009-09-28 22:45:37 +00:00
parent 986c9dccb9
commit 083f26fef4
6 changed files with 172 additions and 49 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

72
liblwgeom/cunit/cu_geodetic.c
View file

@ -28,9 +28,9 @@ CU_pSuite register_geodetic_suite(void)
if (
(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_gserialized_get_gbox_geocentric()", test_gserialized_get_gbox_geocentric)) ||
(NULL == CU_add_test(pSuite, "test_edge_intersection()", test_edge_intersection))
)
{
CU_cleanup_registry();
return NULL;
@ -172,7 +172,6 @@ void test_gserialized_get_gbox_geocentric(void)
/*
* Build LWGEOM on top of *aligned* structure so we can use the read-only
* point access methods on them.
*/
static LWGEOM* lwgeom_over_gserialized(char *wkt)
{
LWGEOM *lwg;
@ -183,50 +182,43 @@ static LWGEOM* lwgeom_over_gserialized(char *wkt)
lwgeom_free(lwg);
return lwgeom_from_gserialized(g);
}
*/
void test_gbox_calculation(void)
void test_edge_intersection(void)
{
LWGEOM *geom;
int i = 0;
GBOX *gbox = gbox_new(gflags(0,0,0));
BOX3D *box3d;
GEOGRAPHIC_EDGE e1, e2;
GEOGRAPHIC_POINT g;
int rv;
char ewkt[][512] = {
"POINT(0 0.2)",
"LINESTRING(-1 -1,-1 2.5,2 2,2 -1)",
"SRID=1;MULTILINESTRING((-1 -1,-1 2.5,2 2,2 -1),(-1 -1,-1 2.5,2 2,2 -1),(-1 -1,-1 2.5,2 2,2 -1),(-1 -1,-1 2.5,2 2,2 -1))",
"POLYGON((-1 -1,-1 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0))",
"SRID=4326;MULTIPOLYGON(((-1 -1,-1 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0),(-0.5 -0.5,-0.5 -0.4,-0.4 -0.4,-0.4 -0.5,-0.5 -0.5)),((-1 -1,-1 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0),(-0.5 -0.5,-0.5 -0.4,-0.4 -0.4,-0.4 -0.5,-0.5 -0.5)))",
"SRID=4326;GEOMETRYCOLLECTION(POINT(0 1),POLYGON((-1 -1,-1 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0)),MULTIPOLYGON(((-1 -1,-1 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0),(-0.5 -0.5,-0.5 -0.4,-0.4 -0.4,-0.4 -0.5,-0.5 -0.5))))",
"POINT(0 220.2)",
"LINESTRING(-1 -1,-1231 2.5,2 2,2 -1)",
"SRID=1;MULTILINESTRING((-1 -131,-1 2.5,2 2,2 -1),(-1 -1,-1 2.5,2 2,2 -1),(-1 -1,-1 2.5,2 2,2 -1),(-1 -1,-1 2.5,2 2,2 -1))",
"POLYGON((-1 -1,-1 2.5,2 2,2 -133,-1 -1),(0 0,0 1,1 1,1 0,0 0))",
"SRID=4326;MULTIPOLYGON(((-1 -1,-1 2.5,211 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0),(-0.5 -0.5,-0.5 -0.4,-0.4 -0.4,-0.4 -0.5,-0.5 -0.5)),((-1 -1,-1 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0),(-0.5 -0.5,-0.5 -0.4,-0.4 -0.4,-0.4 -0.5,-0.5 -0.5)))",
"SRID=4326;GEOMETRYCOLLECTION(POINT(0 1),POLYGON((-1 -1,-1111 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0)),MULTIPOLYGON(((-1 -1,-1 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0),(-0.5 -0.5,-0.5 -0.4,-0.4 -0.4,-0.4 -0.5,-0.5 -0.5))))",
};
for( i = 0; i < 6; i++ )
{
geom = lwgeom_over_gserialized(ewkt[i]);
lwgeom_calculate_gbox(geom, gbox);
box3d = lwgeom_compute_box3d(geom);
//printf("%g %g\n", gbox->xmin, box3d->xmin);
CU_ASSERT_EQUAL(gbox->xmin, box3d->xmin);
CU_ASSERT_EQUAL(gbox->xmax, box3d->xmax);
CU_ASSERT_EQUAL(gbox->ymin, box3d->ymin);
CU_ASSERT_EQUAL(gbox->ymax, box3d->ymax);
lwfree(box3d);
}
e1.start.lon = -1.0;
e1.start.lat = 45.0;
e1.end.lon = 1.0;
e1.end.lat = 45.0;
lwfree(gbox);
e2.start.lon = 0.0;
e2.start.lat = 44.0;
e2.end.lon = 0.0;
e2.end.lat = 46.0;
edge_deg2rad(&e1);
edge_deg2rad(&e2);
rv = edge_intersection(e1, e2, &g);
point_rad2deg(&g);
// printf("g = (%.9g %.9g)\n", g.lon, g.lat);
// printf("rv = %d\n", rv);
// CU_ASSERT_DOUBLE_EQUAL(g.lat, 45.0, 0.00001);
// CU_ASSERT_DOUBLE_EQUAL(g.lon, 0.0, 0.00001);
// CU_ASSERT_EQUAL(rv, LW_TRUE);
e2.end.lat = 42.0;
rv = edge_intersection(e1, e2, &g);
CU_ASSERT_EQUAL(rv, LW_FALSE);
}

1
liblwgeom/cunit/cu_geodetic.h
View file

@ -28,3 +28,4 @@ void test_signum(void);
void test_gbox_from_spherical_coordinates(void);
void test_gserialized_get_gbox_geocentric(void);
void test_gbox_calculation(void);
void test_edge_intersection(void);

42
liblwgeom/cunit/cu_libgeom.c
View file

@ -36,7 +36,9 @@ CU_pSuite register_libgeom_suite(void)
(NULL == CU_add_test(pSuite, "test_geometry_type_from_string()", test_geometry_type_from_string)) ||
(NULL == CU_add_test(pSuite, "test_lwgeom_check_geodetic()", test_lwgeom_check_geodetic)) ||
(NULL == CU_add_test(pSuite, "test_lwgeom_count_vertices()", test_lwgeom_count_vertices)) ||
(NULL == CU_add_test(pSuite, "test_on_gser_lwgeom_count_vertices()", test_on_gser_lwgeom_count_vertices))
(NULL == CU_add_test(pSuite, "test_on_gser_lwgeom_count_vertices()", test_on_gser_lwgeom_count_vertices)) ||
(NULL == CU_add_test(pSuite, "test_gbox_calculation()", test_gbox_calculation))
)
{
CU_cleanup_registry();
@ -417,3 +419,41 @@ void test_on_gser_lwgeom_count_vertices(void)
lwfree(g_ser1);
}
void test_gbox_calculation(void)
{
LWGEOM *geom;
int i = 0;
GBOX *gbox = gbox_new(gflags(0,0,0));
BOX3D *box3d;
char ewkt[][512] = {
"POINT(0 0.2)",
"LINESTRING(-1 -1,-1 2.5,2 2,2 -1)",
"SRID=1;MULTILINESTRING((-1 -1,-1 2.5,2 2,2 -1),(-1 -1,-1 2.5,2 2,2 -1),(-1 -1,-1 2.5,2 2,2 -1),(-1 -1,-1 2.5,2 2,2 -1))",
"POLYGON((-1 -1,-1 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0))",
"SRID=4326;MULTIPOLYGON(((-1 -1,-1 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0),(-0.5 -0.5,-0.5 -0.4,-0.4 -0.4,-0.4 -0.5,-0.5 -0.5)),((-1 -1,-1 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0),(-0.5 -0.5,-0.5 -0.4,-0.4 -0.4,-0.4 -0.5,-0.5 -0.5)))",
"SRID=4326;GEOMETRYCOLLECTION(POINT(0 1),POLYGON((-1 -1,-1 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0)),MULTIPOLYGON(((-1 -1,-1 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0),(-0.5 -0.5,-0.5 -0.4,-0.4 -0.4,-0.4 -0.5,-0.5 -0.5))))",
"POINT(0 220.2)",
"LINESTRING(-1 -1,-1231 2.5,2 2,2 -1)",
"SRID=1;MULTILINESTRING((-1 -131,-1 2.5,2 2,2 -1),(-1 -1,-1 2.5,2 2,2 -1),(-1 -1,-1 2.5,2 2,2 -1),(-1 -1,-1 2.5,2 2,2 -1))",
"POLYGON((-1 -1,-1 2.5,2 2,2 -133,-1 -1),(0 0,0 1,1 1,1 0,0 0))",
"SRID=4326;MULTIPOLYGON(((-1 -1,-1 2.5,211 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0),(-0.5 -0.5,-0.5 -0.4,-0.4 -0.4,-0.4 -0.5,-0.5 -0.5)),((-1 -1,-1 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0),(-0.5 -0.5,-0.5 -0.4,-0.4 -0.4,-0.4 -0.5,-0.5 -0.5)))",
"SRID=4326;GEOMETRYCOLLECTION(POINT(0 1),POLYGON((-1 -1,-1111 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0)),MULTIPOLYGON(((-1 -1,-1 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0),(-0.5 -0.5,-0.5 -0.4,-0.4 -0.4,-0.4 -0.5,-0.5 -0.5))))",
};
for( i = 0; i < 6; i++ )
{
geom = lwgeom_over_gserialized(ewkt[i]);
lwgeom_calculate_gbox(geom, gbox);
box3d = lwgeom_compute_box3d(geom);
//printf("%g %g\n", gbox->xmin, box3d->xmin);
CU_ASSERT_EQUAL(gbox->xmin, box3d->xmin);
CU_ASSERT_EQUAL(gbox->xmax, box3d->xmax);
CU_ASSERT_EQUAL(gbox->ymin, box3d->ymin);
CU_ASSERT_EQUAL(gbox->ymax, box3d->ymax);
lwfree(box3d);
}
lwfree(gbox);
}

1
liblwgeom/cunit/cu_libgeom.h
View file

@ -36,3 +36,4 @@ void test_lwgeom_check_geodetic(void);
void test_lwgeom_count_vertices(void);
void test_on_gser_lwgeom_count_vertices(void);
void test_gbox_serialized_size(void);
void test_gbox_calculation(void);

99
liblwgeom/lwgeodetic.c
View file

@ -11,6 +11,34 @@
#include "lwgeodetic.h"
void edge_deg2rad(GEOGRAPHIC_EDGE *e)
{
(e->start).lat = deg2rad((e->start).lat);
(e->end).lat = deg2rad((e->end).lat);
(e->start).lon = deg2rad((e->start).lon);
(e->end).lon = deg2rad((e->end).lon);
}
void edge_rad2deg(GEOGRAPHIC_EDGE *e)
{
(e->start).lat = rad2deg((e->start).lat);
(e->end).lat = rad2deg((e->end).lat);
(e->start).lon = rad2deg((e->start).lon);
(e->end).lon = rad2deg((e->end).lon);
}
void point_deg2rad(GEOGRAPHIC_POINT *p)
{
p->lat = deg2rad(p->lat);
p->lon = deg2rad(p->lon);
}
void point_rad2deg(GEOGRAPHIC_POINT *p)
{
p->lat = rad2deg(p->lat);
p->lon = rad2deg(p->lon);
}
/**
* 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
@ -73,6 +101,17 @@ static double inline dot_product(POINT3D p1, POINT3D p2)
return (p1.x*p2.x) + (p1.y*p2.y) + (p1.z*p2.z);
}
/**
* Calculate the cross product of two vectors
*/
static void inline cross_product(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;
n->z = a.x * b.y - a.y * b.x;
return;
}
/**
* Normalize to a unit vector.
*/
@ -92,16 +131,14 @@ static void inline normalize(POINT3D *p)
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;
n->z = a.x * b.y - a.y * b.x;
cross_product(a, b, n);
normalize(n);
return;
}
/**
* Computes the cross product of two unit vectors using their lat, lng representations.
* Good for small distances between p and q.
* Computes the cross product of two vectors using their lat, lng representations.
* Good even for small distances between p and q.
*/
void robust_cross_product(GEOGRAPHIC_POINT p, GEOGRAPHIC_POINT q, POINT3D *a)
{
@ -152,10 +189,20 @@ int edge_point_on_plane(GEOGRAPHIC_EDGE e, GEOGRAPHIC_POINT p)
*/
int edge_contains_longitude(GEOGRAPHIC_EDGE e, GEOGRAPHIC_POINT p)
{
LWDEBUGF(4, "e.start == GPOINT(%.6g %.6g) ", e.start.lat, e.start.lon);
LWDEBUGF(4, "e.end == GPOINT(%.6g %.6g) ", e.end.lat, e.end.lon);
LWDEBUGF(4, "p == GPOINT(%.6g %.6g) ", p.lat, p.lon);
if( e.start.lon < p.lon && p.lon < e.end.lon )
{
LWDEBUG(4, "returning TRUE");
return LW_TRUE;
}
if( e.end.lon < p.lon && p.lon < e.start.lon )
{
LWDEBUG(4, "returning TRUE");
return LW_TRUE;
}
LWDEBUG(4, "returning FALSE");
return LW_FALSE;
}
@ -195,8 +242,15 @@ double sphere_direction(GEOGRAPHIC_POINT s, GEOGRAPHIC_POINT e)
*/
int edge_contains_point(GEOGRAPHIC_EDGE e, GEOGRAPHIC_POINT p)
{
if( edge_point_on_plane(e, p) && edge_contains_longitude(e, p) )
return LW_TRUE;
if( edge_point_on_plane(e, p) )
{
LWDEBUG(4, "point is on plane");
if( edge_contains_longitude(e, p) )
{
LWDEBUG(4, "point is on edge");
return LW_TRUE;
}
}
return LW_FALSE;
}
@ -280,6 +334,37 @@ int clairaut_geographic(GEOGRAPHIC_POINT start, GEOGRAPHIC_POINT end, int top, G
return G_SUCCESS;
}
/**
* Returns true if an intersection can be calculated, and places it in *g.
* Returns false otherwise.
*/
int edge_intersection(GEOGRAPHIC_EDGE e1, GEOGRAPHIC_EDGE e2, GEOGRAPHIC_POINT *g)
{
POINT3D ea, eb, v;
robust_cross_product(e1.start, e1.end, &ea);
robust_cross_product(e2.start, e2.end, &eb);
cross_product(ea, eb, &v);
g->lat = atan2(v.z, sqrt(v.x * v.x + v.y * v.y));
g->lon = atan2(v.y, v.x);
if( edge_contains_point(e1, *g) && edge_contains_point(e2, *g) )
{
return LW_TRUE;
}
else
{
g->lat = -1.0 * g->lat;
g->lon = g->lon + M_PI;
if( g->lon > M_PI )
{
g->lon = -1.0 * (2.0 * M_PI - g->lon);
}
if( edge_contains_point(e1, *g) && edge_contains_point(e2, *g) )
{
return LW_TRUE;
}
}
return LW_FALSE;
}
/**
* Given a starting location r, a distance and an azimuth

6
liblwgeom/lwgeodetic.h
View file

@ -55,4 +55,8 @@ double sphere_distance(GEOGRAPHIC_POINT s, GEOGRAPHIC_POINT e);
double sphere_direction(GEOGRAPHIC_POINT s, GEOGRAPHIC_POINT e);
int sphere_project(GEOGRAPHIC_POINT r, double distance, double azimuth, GEOGRAPHIC_POINT *n);
int edge_calculate_gbox(GEOGRAPHIC_EDGE e, GBOX *gbox);
int edge_intersection(GEOGRAPHIC_EDGE e1, GEOGRAPHIC_EDGE e2, GEOGRAPHIC_POINT *g);
void edge_deg2rad(GEOGRAPHIC_EDGE *e);
void edge_rad2deg(GEOGRAPHIC_EDGE *e);
void point_deg2rad(GEOGRAPHIC_POINT *p);
void point_rad2deg(GEOGRAPHIC_POINT *p);