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 56ed0f799b801dbdc64732a2479aef0c28cf3115..b76cbc24dee1656307732277e2993ce8a09dd266 100644 (file)
--- 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 b3c6738cafd1b8704514a5c6bc67978cde034c07..ab69336852679e3895a765c25c3d5c7709c5d723 100644 (file)
--- 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 9fb7848ec87e5910bbbc2482296d875209057480..89e196d2b6943be408bbc8e21aa1f33b80f800b2 100644 (file)
--- 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 5ca9dbf8614a6c307998bc3a77852b8135fd8a4a..4cfc86e0664ff8eb0fb6a7c66cf701c40a8c65fa 100644 (file)
--- 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);