Point-in-polygon function for arc-based point-arrays. (#2018)

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

diff --git a/liblwgeom/cunit/cu_ptarray.c b/liblwgeom/cunit/cu_ptarray.c
index 2ede179a19dc62043ad553eb918137a1bea75556..b1428dd954493d0cc8ff488e996075abce43e23d 100644 (file)
--- a/liblwgeom/cunit/cu_ptarray.c
+++ b/liblwgeom/cunit/cu_ptarray.c
@@ -463,6 +463,114 @@ static void test_ptarray_contains_point()
        lwline_free(lwline);
 }
 
+
+static void test_ptarray_contains_point_arc() 
+{
+/* int ptarray_contains_point_arc(const POINTARRAY *pa, const POINT2D *pt) */
+
+       LWLINE *lwline;
+       POINTARRAY *pa;
+       POINT2D pt;
+       int rv;
+       
+       /* Collection of semi-circles surrounding unit square */
+       lwline = lwgeom_as_lwline(lwgeom_from_text("LINESTRING(-1 -1, -2 0, -1 1, 0 2, 1 1, 2 0, 1 -1, 0 -2, -1 -1)"));
+       pa = lwline->points;
+
+       /* Point in middle of square */
+       pt.x = 0;
+       pt.y = 0;
+       rv = ptarray_contains_point_arc(pa, &pt);
+       CU_ASSERT_EQUAL(rv, 1);
+       
+       /* Point in left lobe */
+       pt.x = -1.1;
+       pt.y = 0.1;
+       rv = ptarray_contains_point_arc(pa, &pt);
+       CU_ASSERT_EQUAL(rv, 1); 
+
+       /* Point on boundary of left lobe */
+       pt.x = -1;
+       pt.y = 0;
+       rv = ptarray_contains_point_arc(pa, &pt);
+       CU_ASSERT_EQUAL(rv, 1); 
+
+       /* Point on boundary vertex */
+       pt.x = -1;
+       pt.y = 1;
+       rv = ptarray_contains_point_arc(pa, &pt);
+       CU_ASSERT_EQUAL(rv, 0); 
+
+       /* Point outside */
+       pt.x = -1.5;
+       pt.y = 1.5;
+       rv = ptarray_contains_point_arc(pa, &pt);
+       CU_ASSERT_EQUAL(rv, -1);        
+
+       /* Two-edge ring made up of semi-circles (really, a circle) */
+       lwfree(lwline);
+       lwline = lwgeom_as_lwline(lwgeom_from_text("LINESTRING(-1 0, 0 1, 1 0, 0 -1, -1 0)"));
+       pa = lwline->points;
+       
+       /* Point outside */
+       pt.x = -1.5;
+       pt.y = 1.5;
+       rv = ptarray_contains_point_arc(pa, &pt);
+       CU_ASSERT_EQUAL(rv, -1);        
+
+       /* Point inside */
+       pt.x = -0.2;
+       pt.y = 0.2;
+       rv = ptarray_contains_point_arc(pa, &pt);
+       CU_ASSERT_EQUAL(rv, 1); 
+
+       /* Point on edge vertex */
+       pt.x = 0;
+       pt.y = 1;
+       rv = ptarray_contains_point_arc(pa, &pt);
+       CU_ASSERT_EQUAL(rv, 0); 
+
+       /* One-edge ring, closed circle */
+       lwfree(lwline);
+       lwline = lwgeom_as_lwline(lwgeom_from_text("LINESTRING(-1 0, 1 0, -1 0)"));
+       pa = lwline->points;
+
+       /* Point inside */
+       pt.x = 0;
+       pt.y = 0;
+       rv = ptarray_contains_point_arc(pa, &pt);
+       CU_ASSERT_EQUAL(rv, 1); 
+
+       /* Point outside */
+       pt.x = 0;
+       pt.y = 2;
+       rv = ptarray_contains_point_arc(pa, &pt);
+       CU_ASSERT_EQUAL(rv, -1);        
+
+       /* Point on boundary */
+       pt.x = 0;
+       pt.y = 1;
+       rv = ptarray_contains_point_arc(pa, &pt);
+       CU_ASSERT_EQUAL(rv, 0); 
+
+       /* Overshort ring */
+       lwfree(lwline);
+       lwline = lwgeom_as_lwline(lwgeom_from_text("LINESTRING(-1 0, 1 0)"));
+       pa = lwline->points;
+       cu_error_msg_reset();
+       rv = ptarray_contains_point_arc(pa, &pt);
+       CU_ASSERT_STRING_EQUAL("ptarray_contains_point_arc called on unclosed ring", cu_error_msg);
+
+       /* Unclosed ring */
+       lwfree(lwline);
+       lwline = lwgeom_as_lwline(lwgeom_from_text("LINESTRING(-1 0, 1 0, 1 0)"));
+       pa = lwline->points;
+       cu_error_msg_reset();
+       rv = ptarray_contains_point_arc(pa, &pt);
+       CU_ASSERT_STRING_EQUAL("ptarray_contains_point_arc called on unclosed ring", cu_error_msg);
+
+       lwline_free(lwline);
+}
 /*
 ** Used by the test harness to register the tests in this file.
 */
@@ -475,6 +583,7 @@ CU_TestInfo ptarray_tests[] =
        PG_TEST(test_ptarray_desegmentize),
        PG_TEST(test_ptarray_insert_point),
        PG_TEST(test_ptarray_contains_point),
+       PG_TEST(test_ptarray_contains_point_arc),
        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 fd74b652a77fc811c2c0b1b02a96bdc5d24bb0c9..43404a090b8f663111eef08aca05775d1f14a87a 100644 (file)
--- a/liblwgeom/liblwgeom_internal.h
+++ b/liblwgeom/liblwgeom_internal.h
@@ -249,6 +249,7 @@ char *geohash_point(double longitude, double latitude, int precision);
 * Point comparisons
 */
 int p4d_same(const POINT4D *p1, const POINT4D *p2);
+int p3d_same(const POINT3D *p1, const POINT3D *p2);
 int p2d_same(const POINT2D *p1, const POINT2D *p2);
 
 /*
@@ -362,7 +363,7 @@ double lw_arc_length(const POINT2D *A1, const POINT2D *A2, const POINT2D *A3);
 double lw_arc_side(const POINT2D *A1, const POINT2D *A2, const POINT2D *A3, const POINT2D *Q);
 int pt_in_ring_2d(const POINT2D *p, const POINTARRAY *ring);
 int ptarray_contains_point(const POINTARRAY *pa, const POINT2D *pt);
-int ptarray_arc_contains_point(const POINTARRAY *pa, const POINT2D *pt);
+int ptarray_contains_point_arc(const POINTARRAY *pa, const POINT2D *pt);
 
 /**
 * Split a line by a point and push components to the provided multiline.

diff --git a/liblwgeom/lwalgorithm.c b/liblwgeom/lwalgorithm.c
index dd7307a756a2b3d2be6494558017733ba3df3aa0..4dff6f86782b9c4edbaeb4ce9d62a9b14ea15cc8 100644 (file)
--- a/liblwgeom/lwalgorithm.c
+++ b/liblwgeom/lwalgorithm.c
@@ -24,6 +24,33 @@ int signum(double n)
        return 0;
 }
 
+int
+p4d_same(const POINT4D *p1, const POINT4D *p2)
+{
+       if( FP_EQUALS(p1->x,p2->x) && FP_EQUALS(p1->y,p2->y) && FP_EQUALS(p1->z,p2->z) && FP_EQUALS(p1->m,p2->m) )
+               return LW_TRUE;
+       else
+               return LW_FALSE;
+}
+
+int
+p3d_same(const POINT3D *p1, const POINT3D *p2)
+{
+       if( FP_EQUALS(p1->x,p2->x) && FP_EQUALS(p1->y,p2->y) && FP_EQUALS(p1->z,p2->z) )
+               return LW_TRUE;
+       else
+               return LW_FALSE;
+}
+
+int
+p2d_same(const POINT2D *p1, const POINT2D *p2)
+{
+       if( FP_EQUALS(p1->x,p2->x) && FP_EQUALS(p1->y,p2->y) )
+               return LW_TRUE;
+       else
+               return LW_FALSE;
+}
+
 /**
 * lw_segment_side()
 *
@@ -34,7 +61,7 @@ int signum(double n)
 int lw_segment_side(const POINT2D *p1, const POINT2D *p2, const POINT2D *q)
 {
        double side = ( (q->x - p1->x) * (p2->y - p1->y) - (p2->x - p1->x) * (q->y - p1->y) );
-       if ( FP_IS_ZERO(side) )
+       if ( side == 0.0 )
                return 0;
        else
                return signum(side);
@@ -144,7 +171,7 @@ lw_arc_length(const POINT2D *A1, const POINT2D *A2, const POINT2D *A3)
        }
 
        /* Length as proportion of circumference */
-       return circumference_A* (angle / (2*M_PI));
+       return circumference_A * (angle / (2*M_PI));
 }
 
 double lw_arc_side(const POINT2D *A1, const POINT2D *A2, const POINT2D *A3, const POINT2D *Q)
@@ -164,6 +191,16 @@ double lw_arc_side(const POINT2D *A1, const POINT2D *A2, const POINT2D *A3, cons
                
        d = distance2d_pt_pt(Q, &C);
        
+       /* Q is on the arc boundary */
+       if ( d == radius_A && side_Q == side_A2 )
+       {       
+               return 0;
+       }
+       
+       /* 
+       * Q is inside the arc boundary, so it's not on the side we 
+       * might think from examining only the end points
+       */
        if ( d < radius_A && side_Q == side_A2 )
        {
                side_Q *= -1;

diff --git a/liblwgeom/ptarray.c b/liblwgeom/ptarray.c
index 35008b7ca01933aac90fb6b560f72aed90c1ee83..0b76426680c552dec7818c094b4d45c32fdeb64f 100644 (file)
--- a/liblwgeom/ptarray.c
+++ b/liblwgeom/ptarray.c
@@ -690,7 +690,7 @@ ptarray_is_closed_z(const POINTARRAY *in)
 }
 
 /**
-* Return -1 if the point is inside the POINTARRAY, 1 if it is outside,
+* Return 1 if the point is inside the POINTARRAY, -1 if it is outside,
 * and 0 if it is on the boundary.
 */
 int 
@@ -737,7 +737,7 @@ ptarray_contains_point(const POINTARRAY *pa, const POINT2D *pt)
                */
                if ( (side == 0) && lw_pt_in_seg(pt, seg1, seg2) )
                {
-                               return 0;
+                       return 0;
                }
 
                /*
@@ -774,34 +774,62 @@ ptarray_contains_point(const POINTARRAY *pa, const POINT2D *pt)
 }
 
 /**
-* Return -1 if the point is inside the POINTARRAY, 1 if it is outside,
+* For POINTARRAYs representing CIRCULARSTRINGS. That is, linked triples
+* with each triple being control points of a circular arc. Such
+* POINTARRAYs have an odd number of vertices.
+*
+* Return 1 if the point is inside the POINTARRAY, -1 if it is outside,
 * and 0 if it is on the boundary.
 */
 int 
-ptarray_arc_contains_point(const POINTARRAY *pa, const POINT2D *pt)
+ptarray_contains_point_arc(const POINTARRAY *pa, const POINT2D *pt)
 {
        int wn = 0;
-       int i;
-       double side;
+       int i, side;
        const POINT2D *seg1;
        const POINT2D *seg2;
        const POINT2D *seg3;
        GBOX gbox;
 
        /* Check for not an arc ring (always have odd # of points) */
-       if ( pa->npoints % 2 )
-               lwerror("ptarray_arc_contains_point called with even number of points");
-       
-       /* Check for unclosed */
-       seg1 = getPoint2d_cp(pa, 0);
-       seg2 = getPoint2d_cp(pa, pa->npoints-1);
+       if ( (pa->npoints % 2) == 0 )
+               lwerror("ptarray_contains_point_arc called with even number of points");
 
-       if ( ! p2d_same(seg1, seg2) )
-               lwerror("ptarray_contains_point called on unclosed ring");
+       /* Check for not an arc ring (always have >= 3 points) */
+       if ( pa->npoints < 3 )
+               lwerror("ptarray_contains_point_arc called too-short pointarray");
 
-       if ( pa->npoints == 3 && p2d_same(seg1, seg2) )
+       /* Check for unclosed case */
+       seg1 = getPoint2d_cp(pa, 0);
+       seg3 = getPoint2d_cp(pa, pa->npoints-1);
+       if ( ! p2d_same(seg1, seg3) )
+       {
+               lwerror("ptarray_contains_point_arc called on unclosed ring");
+       } 
+       /* OK, it's closed. Is it just one circle? */
+       else if ( pa->npoints == 3 )
+       {
+               double radius, d;
+               POINT2D c;
+               seg2 = getPoint2d_cp(pa, 1);
+               
+               /* Wait, it's just a point, so it can't contain anything */
+               if ( lw_arc_is_pt(seg1, seg2, seg3) )
+                       return -1;
+                       
+               /* See if the point is within the circle radius */
+               radius = lw_arc_center(seg1, seg2, seg3, &c);
+               d = distance2d_pt_pt(pt, &c);
+               if ( FP_EQUALS(d, radius) )
+                       return 0; /* Boundary of circle */
+               else if ( d < radius ) 
+                       return 1; /* Inside circle */
+               else 
+                       return -1; /* Outside circle */
+       } 
+       else if ( p2d_same(seg1, pt) )
        {
-               /* TODO handle circlular case */
+               return 0; /* Boundary case */
        }
 
        /* Start on the ring */
@@ -810,6 +838,11 @@ ptarray_arc_contains_point(const POINTARRAY *pa, const POINT2D *pt)
        {
                seg3 = getPoint2d_cp(pa, i);
                
+               /* Catch an easy boundary case */
+               if( p2d_same(seg3, pt) )
+                       return 0;
+               
+               /* Skip arcs that have no size */
                if ( lw_arc_is_pt(seg1, seg2, seg3) )
                {
                        seg1 = seg2;
@@ -817,9 +850,9 @@ ptarray_arc_contains_point(const POINTARRAY *pa, const POINT2D *pt)
                        continue;
                }
                
-               lw_arc_calculate_gbox_cartesian_2d(seg1, seg2, seg3, &gbox);
                
                /* Only test segments in our vertical range */
+               lw_arc_calculate_gbox_cartesian_2d(seg1, seg2, seg3, &gbox);
                if ( pt->y > gbox.ymax || pt->y < gbox.ymin ) 
                {
                        seg1 = seg2;
@@ -827,48 +860,22 @@ ptarray_arc_contains_point(const POINTARRAY *pa, const POINT2D *pt)
                        continue;
                }
                
-               side = lw_segment_side(seg1, seg3, pt);
-               
+               side = lw_arc_side(seg1, seg2, seg3, pt);
                
-               /* Going "up"! */
-               if ( (seg1->y <= pt->y) && (pt->y < seg3->y) )
-               {
-               }
-
-               /* Going "down"! */
-               if ( (seg2->y <= pt->y) && (pt->y < seg1->y) )
-               {
-               }
-
-
-
-               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) )
+               /* On the boundary */
+               if ( (side == 0) && lw_pt_in_arc(pt, seg1, seg2, seg3) )
                {
-                               return 0;
+                       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) )
+               
+               /* Going "up"! Point to left of arc. */
+               if ( side < 0 && (seg1->y <= pt->y) && (pt->y < seg3->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) )
+
+               /* Going "down"! */
+               if ( side > 0 && (seg2->y <= pt->y) && (pt->y < seg1->y) )
                {
                        wn--;
                }
@@ -1476,24 +1483,6 @@ ptarray_length(const POINTARRAY *pts)
 }
 
 
-int
-p4d_same(const POINT4D *p1, const POINT4D *p2)
-{
-       if( FP_EQUALS(p1->x,p2->x) && FP_EQUALS(p1->y,p2->y) && FP_EQUALS(p1->z,p2->z) && FP_EQUALS(p1->m,p2->m) )
-               return LW_TRUE;
-       else
-               return LW_FALSE;
-}
-
-int
-p2d_same(const POINT2D *p1, const POINT2D *p2)
-{
-       if( FP_EQUALS(p1->x,p2->x) && FP_EQUALS(p1->y,p2->y) )
-               return LW_TRUE;
-       else
-               return LW_FALSE;
-}
-
 /*
  * Get a pointer to nth point of a POINTARRAY.
  * You cannot safely cast this to a real POINT, due to memory alignment