From: Paul Ramsey Date: Tue, 9 Oct 2012 17:49:03 +0000 (+0000) Subject: Point-in-polygon function for arc-based point-arrays. (#2018) X-Git-Tag: 2.1.0beta2~579 X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=d89c7aff1fc8b062a70380f20c77f7dc1ebb9ce3;p=postgis 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 2ede179a1..b1428dd95 100644 --- 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 fd74b652a..43404a090 100644 --- 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 dd7307a75..4dff6f867 100644 --- 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 35008b7ca..0b7642668 100644 --- 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