From: Paul Ramsey <pramsey@cleverelephant.ca>
Date: Mon, 5 Oct 2009 19:38:50 +0000 (+0000)
Subject: Add untested new p-i-p approach for more testing later.
X-Git-Tag: 1.5.0b1~415
X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=91a2035bdfab64b3867054467d96ccac7102dc06;p=postgis

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 bebddac55..3e65445db 100644
--- 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 ab6933685..364d61f9d 100644
--- 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 56d16c49d..ad814ef3c 100644
--- 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 31d7e4704..42bc1a2de 100644
--- 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);