Add ptarray_contains_point to ptarray file, so that all other liblwgeom functions...

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

diff --git a/liblwgeom/cunit/cu_ptarray.c b/liblwgeom/cunit/cu_ptarray.c
index 7ecec2dfd82c730f7f7cf1f51158ccdce448b771..2ede179a19dc62043ad553eb918137a1bea75556 100644 (file)
--- a/liblwgeom/cunit/cu_ptarray.c
+++ b/liblwgeom/cunit/cu_ptarray.c
@@ -374,6 +374,95 @@ static void test_ptarray_desegmentize()
        
 }
 
+static void test_ptarray_contains_point() 
+{
+/* int ptarray_contains_point(const POINTARRAY *pa, const POINT2D *pt) */
+
+       LWLINE *lwline;
+       POINTARRAY *pa;
+       POINT2D pt;
+       int rv;
+       
+       lwline = lwgeom_as_lwline(lwgeom_from_text("LINESTRING(0 0, 0 1, 1 1, 1 0, 0 0)"));
+       pa = lwline->points;
+
+       /* Point in middle of square */
+       pt.x = 0.5;
+       pt.y = 0.5;
+       rv = ptarray_contains_point(pa, &pt);
+       CU_ASSERT_EQUAL(rv, 1);
+       
+       /* Point on left edge of square */
+       pt.x = 0;
+       pt.y = 0.5;
+       rv = ptarray_contains_point(pa, &pt);
+       CU_ASSERT_EQUAL(rv, 0);
+
+       /* Point on top edge of square */
+       pt.x = 0.5;
+       pt.y = 1;
+       rv = ptarray_contains_point(pa, &pt);
+       CU_ASSERT_EQUAL(rv, 0);
+
+       /* Point on bottom left corner of square */
+       pt.x = 0;
+       pt.y = 0;
+       rv = ptarray_contains_point(pa, &pt);
+       CU_ASSERT_EQUAL(rv, 0);
+
+       /* Point on top left corner of square */
+       pt.x = 0;
+       pt.y = 1;
+       rv = ptarray_contains_point(pa, &pt);
+       CU_ASSERT_EQUAL(rv, 0);
+
+       /* Point outside top left corner of square */
+       pt.x = -0.1;
+       pt.y = 1;
+       rv = ptarray_contains_point(pa, &pt);
+       CU_ASSERT_EQUAL(rv, -1);
+
+       /* Point outside top left corner of square */
+       pt.x = 0;
+       pt.y = 1.1;
+       rv = ptarray_contains_point(pa, &pt);
+       CU_ASSERT_EQUAL(rv, -1);
+
+       /* Point outside left side of square */
+       pt.x = -0.2;
+       pt.y = 0.5;
+       rv = ptarray_contains_point(pa, &pt);
+       CU_ASSERT_EQUAL(rv, -1);
+       
+       lwline_free(lwline);
+       lwline = lwgeom_as_lwline(lwgeom_from_text("LINESTRING(0 0, 1 1, 2 0, 0 0)"));
+       pa = lwline->points;
+       
+       /* Point outside grazing top of triangle */
+       pt.x = 0;
+       pt.y = 1;
+       rv = ptarray_contains_point(pa, &pt);
+       CU_ASSERT_EQUAL(rv, -1);
+
+       lwline_free(lwline);
+       lwline = lwgeom_as_lwline(lwgeom_from_text("LINESTRING(0 0, 0 4, 1 4, 2 2, 3 4, 4 4, 4 0, 0 0)"));
+       pa = lwline->points;
+
+       /* Point outside grazing top of triangle */
+       pt.x = 1;
+       pt.y = 2;
+       rv = ptarray_contains_point(pa, &pt);
+       CU_ASSERT_EQUAL(rv, 1);
+
+       /* Point outside grazing top of triangle */
+       pt.x = 3;
+       pt.y = 2;
+       rv = ptarray_contains_point(pa, &pt);
+       CU_ASSERT_EQUAL(rv, 1);
+
+       lwline_free(lwline);
+}
+
 /*
 ** Used by the test harness to register the tests in this file.
 */
@@ -385,6 +474,7 @@ CU_TestInfo ptarray_tests[] =
        PG_TEST(test_ptarray_isccw),
        PG_TEST(test_ptarray_desegmentize),
        PG_TEST(test_ptarray_insert_point),
+       PG_TEST(test_ptarray_contains_point),
        CU_TEST_INFO_NULL
 };
 CU_SuiteInfo ptarray_suite = {"ptarray", NULL, NULL, ptarray_tests };

diff --git a/liblwgeom/liblwgeom_internal.h b/liblwgeom/liblwgeom_internal.h
index e56cdb9e2edc1e152150ba1e6abd602fcad14685..df90ead34071d7a629442ef61c9caee144546f12 100644 (file)
--- a/liblwgeom/liblwgeom_internal.h
+++ b/liblwgeom/liblwgeom_internal.h
@@ -353,6 +353,10 @@ int lwpsurface_is_closed(const LWPSURFACE *psurface);
 int lwtin_is_closed(const LWTIN *tin);
 
 
+
+int lw_pt_in_seg(const POINT2D *P, const POINT2D *A1, const POINT2D *A2);
+int ptarray_contains_point(const POINTARRAY *pa, const POINT2D *pt);
+
 /**
 * Split a line by a point and push components to the provided multiline.
 * If the point doesn't split the line, push nothing to the container.

diff --git a/liblwgeom/lwalgorithm.c b/liblwgeom/lwalgorithm.c
index 01f043caefecd2dafa5ba2110c9d316c4201daf1..5b0dc465c96792e09ab9d48139d5dcd31fb17ed3 100644 (file)
--- a/liblwgeom/lwalgorithm.c
+++ b/liblwgeom/lwalgorithm.c
@@ -25,11 +25,11 @@ int signum(double n)
 }
 
 /**
-** lw_segment_side()
-**
-** Return < 0.0 if point Q is left of segment P
-** Return > 0.0 if point Q is right of segment P
-** Return = 0.0 if point Q in on segment P
+* lw_segment_side()
+*
+* Return < 0.0 if point Q is left of segment P
+* Return > 0.0 if point Q is right of segment P
+* Return = 0.0 if point Q in on segment P
 */
 double lw_segment_side(const POINT2D *p1, const POINT2D *p2, const POINT2D *q)
 {
@@ -41,6 +41,7 @@ double lw_segment_side(const POINT2D *p1, const POINT2D *p2, const POINT2D *q)
 }
 
 
+
 int lw_segment_envelope_intersects(const POINT2D *p1, const POINT2D *p2, const POINT2D *q1, const POINT2D *q2)
 {
        double minq=FP_MIN(q1->x,q2->x);

diff --git a/liblwgeom/measures.c b/liblwgeom/measures.c
index e7410df6f4b4d4faf5d09bd78e43ab97fb978e50..f42cf2eb0c74936005330e325e80fe4559d1e432 100644 (file)
--- a/liblwgeom/measures.c
+++ b/liblwgeom/measures.c
@@ -840,11 +840,11 @@ lw_pt_in_arc(const POINT2D *P, const POINT2D *A1, const POINT2D *A2, const POINT
 * Returns true if P is between A1/A2. Only makes sense if P has already been
 * deterined to be on the line defined by A1/A2.
 */
-static int
+int
 lw_pt_in_seg(const POINT2D *P, const POINT2D *A1, const POINT2D *A2)
 {
-       return ((A1->x < P->x && P->x < A2->x) || (A1->x > P->x && P->x > A2->x)) ||
-              ((A1->y < P->y && P->y < A2->y) || (A1->y > P->y && P->y > A2->y));
+       return ((A1->x <= P->x && P->x < A2->x) || (A1->x >= P->x && P->x > A2->x)) ||
+              ((A1->y <= P->y && P->y < A2->y) || (A1->y >= P->y && P->y > A2->y));
 }
 
 /**

diff --git a/liblwgeom/ptarray.c b/liblwgeom/ptarray.c
index b7a4504acf570a62435d775c184a97742e730019..e50270252196d23b4b37a01e0a7037b0578f393d 100644 (file)
--- a/liblwgeom/ptarray.c
+++ b/liblwgeom/ptarray.c
@@ -689,6 +689,86 @@ ptarray_is_closed_z(const POINTARRAY *in)
                return ptarray_is_closed_2d(in);
 }
 
+/**
+* Return -1 if the point is inside the POINTARRAY, 1 if it is outside,
+* and 0 if it is on the boundary.
+*/
+int 
+ptarray_contains_point(const POINTARRAY *pa, const POINT2D *pt)
+{
+       int wn = 0;
+       int i;
+       double side;
+       const POINT2D *seg1;
+       const POINT2D *seg2;
+       double ymin, ymax;
+
+       seg1 = getPoint2d_cp(pa, 0);
+       seg2 = getPoint2d_cp(pa, pa->npoints-1);
+       if ( ! p2d_same(seg1, seg2) )
+               lwerror("ptarray_contains_point called on unclosed ring");
+       
+       for ( i=1; i < pa->npoints; i++ )
+       {
+               seg2 = getPoint2d_cp(pa, i);
+               
+               /* Zero length segments are ignored. */
+               if ( seg1->x == seg2->x && seg1->y == seg2->y )
+                       continue;
+                       
+               ymin = FP_MIN(seg1->y, seg2->y);
+               ymax = FP_MAX(seg1->y, seg2->y);
+               
+               /* Only test segments in our vertical range */
+               if ( pt->y > ymax || pt->y < ymin ) 
+               {
+                       seg1 = seg2;
+                       continue;
+               }
+
+               side = lw_segment_side(seg1, seg2, pt);
+
+               /* 
+               * A point on the boundary of a ring is not contained. 
+               * WAS: if (fabs(side) < 1e-12), see #852 
+               */
+               if ( (side == 0.0) && lw_pt_in_seg(pt, seg1, seg2) )
+               {
+                               return 0;
+               }
+
+               /*
+               * If the point is to the left of the line, and it's rising,
+               * then the line is to the right of the point and
+               * circling counter-clockwise, so incremement.
+               */
+               if ( (side < 0) && (seg1->y <= pt->y) && (pt->y < seg2->y) )
+               {
+                       wn++;
+               }
+               
+               /*
+               * If the point is to the right of the line, and it's falling,
+               * then the line is to the right of the point and circling
+               * clockwise, so decrement.
+               */
+               else if ( (side > 0) && (seg2->y <= pt->y) && (pt->y < seg1->y) )
+               {
+                       wn--;
+               }
+               
+               seg1 = seg2;
+       }
+
+       /* Outside */
+       if (wn == 0)
+       {
+               return -1;
+       }
+       
+       /* Inside */
+       return 1;
+}
 
 POINTARRAY*
 ptarray_force_dims(const POINTARRAY *pa, int hasz, int hasm)