Add untested new p-i-p approach for more testing later.

git-svn-id: http://svn.osgeo.org/postgis/trunk@4595 b70326c6-7e19-0410-871a-916f4a2858ee

diff --git a/liblwgeom/cunit/cu_geodetic.c b/liblwgeom/cunit/cu_geodetic.c
index bebddac55f983473c5fc999f30e54c0c12c0c61a..3e65445dbdb5a87ef7dc43bf904b6bff4fdc7887 100644 (file)
--- a/liblwgeom/cunit/cu_geodetic.c
+++ b/liblwgeom/cunit/cu_geodetic.c
@@ -34,10 +34,9 @@ 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_ptarray_point_in_ring_winding()", test_ptarray_point_in_ring_winding)) ||
            (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();
@@ -531,4 +530,6 @@ void test_lwgeom_distance_sphere(void)
 
 }
 
-
+void test_ptarray_point_in_ring_winding(void)
+{
+}

diff --git a/liblwgeom/cunit/cu_geodetic.h b/liblwgeom/cunit/cu_geodetic.h
index ab69336852679e3895a765c25c3d5c7709c5d723..364d61f9d623e14eaba1365f8ac6e34f0197b58c 100644 (file)
--- a/liblwgeom/cunit/cu_geodetic.h
+++ b/liblwgeom/cunit/cu_geodetic.h
@@ -32,4 +32,5 @@ 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_ptarray_point_in_ring_winding(void);
 void test_lwgeom_distance_sphere(void);

diff --git a/liblwgeom/lwgeodetic.c b/liblwgeom/lwgeodetic.c
index 56d16c49db5e9559bee2485aeffdf8c27d13af75..ad814ef3c6238c555f6b039180eafd40d4357e47 100644 (file)
--- a/liblwgeom/lwgeodetic.c
+++ b/liblwgeom/lwgeodetic.c
@@ -1260,6 +1260,58 @@ static int ptarray_point_in_ring(POINTARRAY *pa, POINT2D pt_outside, POINT2D pt_
        return LW_FALSE;
 }
 
+/**
+* This routine returns LW_TRUE if the point is inside the ring or on the boundary, LW_FALSE otherwise.
+* The pt_outside must be guaranteed to be outside the ring (use the geography_pt_outside() function
+* to derive one in postgis, or the gbox_pt_outside() function if you don't mind burning CPU cycles
+* building a gbox first).
+*/
+int ptarray_point_in_ring_winding(POINTARRAY *pa, POINT2D pt_to_test) 
+{
+       GEOGRAPHIC_POINT gpt;
+       GEOGRAPHIC_EDGE edge1, edge2;
+       POINT3D norm1, norm2;
+       POINT2D p;
+       double wind_number = 0;
+       int i;
+       
+       /* Null input, not enough points for a ring? You ain't closed! */
+       if( ! pa || pa->npoints < 4 ) 
+               return LW_FALSE;
+
+       /* Set up our test point */
+       geographic_point_init(pt_to_test.x, pt_to_test.y, &gpt);
+       edge1.start = gpt;
+       edge2.start = gpt;
+
+       /* Walk every edge and see if the stab line hits it */
+       for( i = 1; i < pa->npoints; i++ )
+       {
+               getPoint2d_p(pa, i-1, &p);
+               geographic_point_init(p.x, p.y, &(edge1.end));
+               getPoint2d_p(pa, i, &p);
+               geographic_point_init(p.x, p.y, &(edge2.end));
+
+               /* Calculate normals to the great circle planes */
+               robust_cross_product(edge1.start, edge1.end, &norm1);
+               robust_cross_product(edge2.start, edge2.end, &norm2);
+               normalize(&norm1);
+               normalize(&norm2);
+               
+               /* Add the wind */
+               wind_number += acos(dot_product(norm1, norm2));
+               
+       }
+
+       if( wind_number > M_PI )
+       {
+               return LW_TRUE;
+       }
+       
+       return LW_FALSE;
+}
+
+
 static double ptarray_distance_sphere(POINTARRAY *pa1, POINTARRAY *pa2, double tolerance)
 {
        GEOGRAPHIC_EDGE e1, e2;

diff --git a/liblwgeom/lwgeodetic.h b/liblwgeom/lwgeodetic.h
index 31d7e4704571197a3ccdbf86d14c878e6f6ddc30..42bc1a2de286c420e549ebf0663987cea49626a1 100644 (file)
--- a/liblwgeom/lwgeodetic.h
+++ b/liblwgeom/lwgeodetic.h
@@ -72,4 +72,5 @@ void edge_rad2deg(GEOGRAPHIC_EDGE *e);
 void point_deg2rad(GEOGRAPHIC_POINT *p);
 void point_rad2deg(GEOGRAPHIC_POINT *p);
 void geographic_point_init(double lon, double lat, GEOGRAPHIC_POINT *g);
+int ptarray_point_in_ring_winding(POINTARRAY *pa, POINT2D pt_to_test);
 int lwpoly_covers_point2d(const LWPOLY *poly, GBOX *gbox, POINT2D pt_to_test);