From: Paul Ramsey Date: Thu, 1 Oct 2009 05:53:07 +0000 (+0000) Subject: lwgeom sphere distance function and tests for point/linestring X-Git-Tag: 1.5.0b1~440 X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=68dab94cd07bccf4da933c5e4bdccf39d3070a9e;p=postgis lwgeom sphere distance function and tests for point/linestring git-svn-id: http://svn.osgeo.org/postgis/trunk@4570 b70326c6-7e19-0410-871a-916f4a2858ee --- diff --git a/liblwgeom/cunit/cu_geodetic.c b/liblwgeom/cunit/cu_geodetic.c index 56ed0f799..b76cbc24d 100644 --- a/liblwgeom/cunit/cu_geodetic.c +++ b/liblwgeom/cunit/cu_geodetic.c @@ -34,7 +34,10 @@ CU_pSuite register_geodetic_suite(void) (NULL == CU_add_test(pSuite, "test_clairaut()", test_clairaut)) || (NULL == CU_add_test(pSuite, "test_edge_intersection()", test_edge_intersection)) || (NULL == CU_add_test(pSuite, "test_edge_distance_to_point()", test_edge_distance_to_point)) || - (NULL == CU_add_test(pSuite, "test_edge_distance_to_edge()", test_edge_distance_to_edge)) + (NULL == CU_add_test(pSuite, "test_edge_distance_to_edge()", test_edge_distance_to_edge)) || + (NULL == CU_add_test(pSuite, "test_lwgeom_distance_sphere()", test_lwgeom_distance_sphere)) + + ) { CU_cleanup_registry(); @@ -409,7 +412,33 @@ void test_edge_distance_to_edge(void) CU_ASSERT_DOUBLE_EQUAL(c2.lon, 0.0, 0.00001); } +void test_lwgeom_distance_sphere(void) +{ + LWGEOM *lwg1, *lwg2; + double d; + + lwg1 = lwgeom_from_ewkt("LINESTRING(-30 10, -20 5, -10 3, 0 1)", PARSER_CHECK_NONE); + lwg2 = lwgeom_from_ewkt("LINESTRING(-10 -5, -5 0, 5 0, 10 -5)", PARSER_CHECK_NONE); + d = lwgeom_distance_sphere(lwg1, lwg2, 0.0); + CU_ASSERT_DOUBLE_EQUAL(d, M_PI / 180.0, 0.00001); + lwgeom_free(lwg1); + lwgeom_free(lwg2); + + lwg1 = lwgeom_from_ewkt("POINT(-4 1)", PARSER_CHECK_NONE); + lwg2 = lwgeom_from_ewkt("LINESTRING(-10 -5, -5 0, 5 0, 10 -5)", PARSER_CHECK_NONE); + d = lwgeom_distance_sphere(lwg1, lwg2, 0.0); + CU_ASSERT_DOUBLE_EQUAL(d, M_PI / 180.0, 0.00001); + lwgeom_free(lwg1); + lwgeom_free(lwg2); + lwg1 = lwgeom_from_ewkt("POINT(-4 1)", PARSER_CHECK_NONE); + lwg2 = lwgeom_from_ewkt("POINT(-4 -1)", PARSER_CHECK_NONE); + d = lwgeom_distance_sphere(lwg1, lwg2, 0.0); + CU_ASSERT_DOUBLE_EQUAL(d, M_PI / 90.0, 0.00001); + lwgeom_free(lwg1); + lwgeom_free(lwg2); + +} diff --git a/liblwgeom/cunit/cu_geodetic.h b/liblwgeom/cunit/cu_geodetic.h index b3c6738ca..ab6933685 100644 --- a/liblwgeom/cunit/cu_geodetic.h +++ b/liblwgeom/cunit/cu_geodetic.h @@ -32,3 +32,4 @@ void test_gbox_calculation(void); void test_edge_intersection(void); void test_edge_distance_to_point(void); void test_edge_distance_to_edge(void); +void test_lwgeom_distance_sphere(void); diff --git a/liblwgeom/libgeom.h b/liblwgeom/libgeom.h index 9fb7848ec..89e196d2b 100644 --- a/liblwgeom/libgeom.h +++ b/liblwgeom/libgeom.h @@ -374,7 +374,7 @@ extern int lwgeom_calculate_gbox(const LWGEOM *lwgeom, GBOX *gbox); * Calculate the geodetic distance from lwgeom1 to lwgeom2 on the unit sphere. * Pass in a tolerance in radians. */ -extern int lwgeom_distance_sphere(LWGEOM *lwgeom1, LWGEOM *lwgeom2, double tolerance, double *result); +extern double lwgeom_distance_sphere(LWGEOM *lwgeom1, LWGEOM *lwgeom2, double tolerance); /** * New function to read doubles directly from the double* coordinate array diff --git a/liblwgeom/lwgeodetic.c b/liblwgeom/lwgeodetic.c index 5ca9dbf86..4cfc86e06 100644 --- a/liblwgeom/lwgeodetic.c +++ b/liblwgeom/lwgeodetic.c @@ -1077,98 +1077,163 @@ int edge_calculate_gbox(GEOGRAPHIC_EDGE e, GBOX *gbox) return G_SUCCESS; } -static int lw_line_point_distance_sphere(LWLINE *lwline, LWPOINT *lwpt, double tolerance, double *result) +static double ptarray_distance_sphere(POINTARRAY *pa1, POINTARRAY *pa2, double tolerance) { - GEOGRAPHIC_POINT g; - GEOGRAPHIC_EDGE e; + GEOGRAPHIC_EDGE e1, e2; + GEOGRAPHIC_POINT g1, g2; POINT2D p; - int i; + double distance; + int i, j; + + /* Make result really big, so that everything will be smaller than it */ + distance = MAXFLOAT; - /* Initialize our point */ - getPoint2d_p(lwpt->point, 0, &p); - geographic_point_init(p.x, p.y, &g); + /* Empty point arrays? Return negative */ + if ( pa1->npoints == 0 || pa1->npoints == 0 ) + return -1.0; - /* Degenerate one-point line! We handle gracefully. */ - if ( lwline->points->npoints == 1 ) + /* Handle point/point case here */ + if ( pa1->npoints == 1 && pa2->npoints == 1 ) + { + getPoint2d_p(pa1, 0, &p); + geographic_point_init(p.x, p.y, &g1); + getPoint2d_p(pa2, 0, &p); + geographic_point_init(p.x, p.y, &g2); + return sphere_distance(g1, g2); + } + + /* Handle point/line case here */ + if ( pa1->npoints == 1 || pa2->npoints == 1 ) { - GEOGRAPHIC_POINT f; - getPoint2d_p(lwline->points, 0, &p); - geographic_point_init(p.x, p.y, &f); - *result = sphere_distance(g, f); - return LW_TRUE; + /* Handle one/many case here */ + int i; + POINTARRAY *pa_one, *pa_many; + + if( pa1->npoints == 1 ) + { + pa_one = pa1; + pa_many = pa2; + } + else + { + pa_one = pa2; + pa_many = pa1; + } + + /* Initialize our point */ + getPoint2d_p(pa_one, 0, &p); + geographic_point_init(p.x, p.y, &g1); + + /* Iterate through the edges in our line */ + for( i = 1; i < pa_many->npoints; i++ ) + { + double d; + getPoint2d_p(pa_many, i - 1, &p); + geographic_point_init(p.x, p.y, &(e1.start)); + getPoint2d_p(pa_many, i, &p); + geographic_point_init(p.x, p.y, &(e1.end)); + d = edge_distance_to_point(e1, g1, 0); + if( d < distance ) + distance = d; + if( d < tolerance ) + return distance; + } + return distance; } - /* Iterate through the edges in our line */ - for( i = 1; i < lwline->points->npoints; i++ ) + /* Handle line/line case */ + for( i = 1; i < pa1->npoints; i++ ) { - getPoint2d_p(lwline->points, i - 1, &p); - geographic_point_init(p.x, p.y, &(e.start)); - getPoint2d_p(lwline->points, i, &p); - geographic_point_init(p.x, p.y, &(e.end)); - *result = edge_distance_to_point(e, g, 0); - if( *result < tolerance ) - break; + getPoint2d_p(pa1, i - 1, &p); + geographic_point_init(p.x, p.y, &(e1.start)); + getPoint2d_p(pa1, i, &p); + geographic_point_init(p.x, p.y, &(e1.end)); + + for( j = 1; j < pa2->npoints; j++ ) + { + double d; + GEOGRAPHIC_POINT g; + + getPoint2d_p(pa2, j - 1, &p); + geographic_point_init(p.x, p.y, &(e2.start)); + getPoint2d_p(pa2, j, &p); + geographic_point_init(p.x, p.y, &(e2.end)); + + if ( edge_intersection(e1, e2, &g) ) + { + return 0.0; + } + d = edge_distance_to_edge(e1, e2, 0, 0); + if( d < distance ) + distance = d; + if( d < tolerance ) + return distance; + } } - return LW_TRUE; + return distance; } -static int lw_point_point_distance_sphere(LWPOINT *lwpt1, LWPOINT *lwpt2, double *result) -{ - GEOGRAPHIC_POINT g1, g2; - POINT2D p; - getPoint2d_p(lwpt1->point, 0, &p); - geographic_point_init(p.x, p.y, &g1); - getPoint2d_p(lwpt2->point, 0, &p); - geographic_point_init(p.x, p.y, &g2); - *result = sphere_distance(g1, g2); - return LW_TRUE; -} /** * Calculate the distance between two LWGEOMs, using the coordinates are * longitude and latitude. Return immediately when the calulated distance drops * below the tolerance (useful for dwithin calculations). */ -int lwgeom_distance_sphere(LWGEOM *lwgeom1, LWGEOM *lwgeom2, double tolerance, double *result) +double lwgeom_distance_sphere(LWGEOM *lwgeom1, LWGEOM *lwgeom2, double tolerance) { int type1, type2; assert(lwgeom1); assert(lwgeom2); + /* What's the distance to an empty geometry? We don't know. */ if( lwgeom_is_empty(lwgeom1) || lwgeom_is_empty(lwgeom2) ) { - *result = -1.0; - return LW_FALSE; + return -1.0; } - /* Ensure that for LINE/POINT, POLY/POINT, POLY/LINE cases - the lower dimensioned object is always type2 */ - if ( lwgeom_dimensionality(lwgeom1) < lwgeom_dimensionality(lwgeom2) ) - { - LWGEOM *tmp = lwgeom2; - lwgeom2 = lwgeom1; - lwgeom1 = tmp; - } - type1 = TYPE_GETTYPE(lwgeom1->type); type2 = TYPE_GETTYPE(lwgeom2->type); - if( type1 == POINTTYPE && type2 == POINTTYPE ) + /* Point/line combinations can all be handled with simple point array iterations */ + if( ( type1 == POINTTYPE || type1 == LINETYPE ) && + ( type2 == POINTTYPE || type2 == LINETYPE ) ) + { + POINTARRAY *pa1, *pa2; + + if( type1 == POINTTYPE ) + pa1 = ((LWPOINT*)lwgeom1)->point; + else + pa1 = ((LWLINE*)lwgeom1)->points; + + if( type2 == POINTTYPE ) + pa2 = ((LWPOINT*)lwgeom2)->point; + else + pa2 = ((LWLINE*)lwgeom2)->points; + + return ptarray_distance_sphere(pa1, pa2, tolerance); + } + + /* Point/Polygon cases, if point-in-poly, return zero, else return distance. */ + if( ( type1 == POLYGONTYPE && type2 == POINTTYPE ) || + ( type2 == POLYGONTYPE && type1 == POINTTYPE ) ) { - return lw_point_point_distance_sphere( - lwgeom_as_lwpoint(lwgeom1), - lwgeom_as_lwpoint(lwgeom2), result); } - if( type1 == LINETYPE && type2 == POINTTYPE ) + + /* Line/polygon case, if start point-in-poly, return zero, else return distance. */ + if( ( type1 == POLYGONTYPE && type2 == LINETYPE ) || + ( type2 == POLYGONTYPE && type1 == LINETYPE ) ) { - return lw_line_point_distance_sphere( - lwgeom_as_lwline(lwgeom1), - lwgeom_as_lwpoint(lwgeom2), - tolerance, result); + } + /* Polygon/polygon case, if start point-in-poly, return zero, else return distance. */ + if( ( type1 == POLYGONTYPE && type2 == LINETYPE ) || + ( type2 == POLYGONTYPE && type1 == LINETYPE ) ) + { } - return LW_FALSE; + + + return -1.0; } @@ -1212,8 +1277,7 @@ int ptarray_calculate_gbox_geodetic(POINTARRAY *pa, GBOX *gbox) POINT3D out_pt; GEOGRAPHIC_POINT gp; getPoint2d_p(pa, 0, &in_pt); - gp.lon = deg2rad(in_pt.x); - gp.lat = deg2rad(in_pt.y); + geographic_point_init(in_pt.x, in_pt.y, &gp); geog2cart(gp, &out_pt); gbox->xmin = gbox->xmax = out_pt.x; gbox->ymin = gbox->ymax = out_pt.y; @@ -1224,11 +1288,10 @@ int ptarray_calculate_gbox_geodetic(POINTARRAY *pa, GBOX *gbox) for ( i = 1; i < pa->npoints; i++ ) { getPoint2d_p(pa, i-1, &start_pt); + geographic_point_init(start_pt.x, start_pt.y, &(edge.start)); + getPoint2d_p(pa, i, &end_pt); - edge.start.lon = deg2rad(start_pt.x); - edge.start.lat = deg2rad(start_pt.y); - edge.end.lon = deg2rad(end_pt.x); - edge.end.lat = deg2rad(end_pt.y); + geographic_point_init(end_pt.x, end_pt.y, &(edge.end)); edge_calculate_gbox(edge, &edge_gbox);