From: Paul Ramsey Date: Fri, 20 Nov 2009 18:36:38 +0000 (+0000) Subject: Fix graxing case and improve co-linear handling with help from Nicklas Aven (#314) X-Git-Tag: 1.5.0b1~199 X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=74dbaa24293a1f878960fbb5b7b729d9f45d9a7c;p=postgis Fix graxing case and improve co-linear handling with help from Nicklas Aven (#314) git-svn-id: http://svn.osgeo.org/postgis/trunk@4872 b70326c6-7e19-0410-871a-916f4a2858ee --- diff --git a/liblwgeom/cunit/cu_geodetic.c b/liblwgeom/cunit/cu_geodetic.c index 469aad704..09eca8597 100644 --- a/liblwgeom/cunit/cu_geodetic.c +++ b/liblwgeom/cunit/cu_geodetic.c @@ -39,9 +39,9 @@ CU_pSuite register_geodetic_suite(void) (NULL == CU_add_test(pSuite, "test_lwgeom_check_geodetic()", test_lwgeom_check_geodetic)) || (NULL == CU_add_test(pSuite, "test_gbox_calculation()", test_gbox_calculation)) || (NULL == CU_add_test(pSuite, "test_gserialized_from_lwgeom()", test_gserialized_from_lwgeom)) || - (NULL == CU_add_test(pSuite, "test_lwpoly_covers_point2d()", test_lwpoly_covers_point2d)) || (NULL == CU_add_test(pSuite, "test_spheroid_distance()", test_spheroid_distance)) || (NULL == CU_add_test(pSuite, "test_spheroid_area()", test_spheroid_area)) || + (NULL == CU_add_test(pSuite, "test_lwpoly_covers_point2d()", test_lwpoly_covers_point2d)) || (NULL == CU_add_test(pSuite, "test_ptarray_point_in_ring()", test_ptarray_point_in_ring)) ) { @@ -613,6 +613,83 @@ void test_ptarray_point_in_ring(void) POINT2D pt_outside; int result; + /* Simple containment case */ + lwg = lwgeom_from_ewkt("POLYGON((1.0 1.0, 1.0 1.1, 1.1 1.1, 1.1 1.0, 1.0 1.0))", PARSER_CHECK_NONE); + poly = (LWPOLY*)lwg; + pt_to_test.x = 1.05; + pt_to_test.y = 1.05; + pt_outside.x = 1.2; + pt_outside.y = 1.15; + result = ptarray_point_in_ring(poly->rings[0], pt_outside, pt_to_test); + CU_ASSERT_EQUAL(result, LW_TRUE); + lwgeom_free(lwg); + + /* Simple noncontainment case */ + lwg = lwgeom_from_ewkt("POLYGON((1.0 1.0, 1.0 1.1, 1.1 1.1, 1.1 1.0, 1.0 1.0))", PARSER_CHECK_NONE); + poly = (LWPOLY*)lwg; + pt_to_test.x = 1.05; + pt_to_test.y = 1.15; + pt_outside.x = 1.2; + pt_outside.y = 1.2; + result = ptarray_point_in_ring(poly->rings[0], pt_outside, pt_to_test); + CU_ASSERT_EQUAL(result, LW_FALSE); + lwgeom_free(lwg); + + /* Harder noncontainment case */ + lwg = lwgeom_from_ewkt("POLYGON((1.0 1.0, 1.0 1.1, 1.1 1.1, 1.1 1.0, 1.0 1.0))", PARSER_CHECK_NONE); + poly = (LWPOLY*)lwg; + pt_to_test.x = 1.05; + pt_to_test.y = 0.9; + pt_outside.x = 1.2; + pt_outside.y = 1.05; + result = ptarray_point_in_ring(poly->rings[0], pt_outside, pt_to_test); + CU_ASSERT_EQUAL(result, LW_FALSE); + lwgeom_free(lwg); + + /* Harder containment case */ + lwg = lwgeom_from_ewkt("POLYGON((0 0, 0 2, 1 2, 0 3, 2 3, 0 4, 3 5, 0 6, 6 10, 6 1, 0 0))", PARSER_CHECK_NONE); + poly = (LWPOLY*)lwg; + pt_to_test.x = 1.0; + pt_to_test.y = 1.0; + pt_outside.x = 1.0; + pt_outside.y = 10.0; + result = ptarray_point_in_ring(poly->rings[0], pt_outside, pt_to_test); + CU_ASSERT_EQUAL(result, LW_TRUE); + lwgeom_free(lwg); + + /* Point on ring at vertex case */ + lwg = lwgeom_from_ewkt("POLYGON((1.0 1.0, 1.0 1.1, 1.1 1.1, 1.1 1.0, 1.0 1.0))", PARSER_CHECK_NONE); + poly = (LWPOLY*)lwg; + pt_to_test.x = 1.1; + pt_to_test.y = 1.05; + pt_outside.x = 1.2; + pt_outside.y = 1.05; + result = ptarray_point_in_ring(poly->rings[0], pt_outside, pt_to_test); + CU_ASSERT_EQUAL(result, LW_TRUE); + lwgeom_free(lwg); + + /* Point on ring at first vertex case */ + lwg = lwgeom_from_ewkt("POLYGON((1.0 1.0, 1.0 1.1, 1.1 1.1, 1.1 1.0, 1.0 1.0))", PARSER_CHECK_NONE); + poly = (LWPOLY*)lwg; + pt_to_test.x = 1.0; + pt_to_test.y = 1.0; + pt_outside.x = 1.2; + pt_outside.y = 1.05; + result = ptarray_point_in_ring(poly->rings[0], pt_outside, pt_to_test); + CU_ASSERT_EQUAL(result, LW_TRUE); + lwgeom_free(lwg); + + /* Point on ring between vertexes case */ + lwg = lwgeom_from_ewkt("POLYGON((1.0 1.0, 1.0 1.1, 1.1 1.1, 1.1 1.0, 1.0 1.0))", PARSER_CHECK_NONE); + poly = (LWPOLY*)lwg; + pt_to_test.x = 1.0; + pt_to_test.y = 1.1; + pt_outside.x = 1.2; + pt_outside.y = 1.05; + result = ptarray_point_in_ring(poly->rings[0], pt_outside, pt_to_test); + CU_ASSERT_EQUAL(result, LW_TRUE); + lwgeom_free(lwg); + /* Co-linear crossing case for point-in-polygon test, should return LW_TRUE */ lwg = lwgeom_from_ewkt("POLYGON((1.0 1.0, 1.0 1.1, 1.1 1.1, 1.1 1.2, 1.2 1.2, 1.2 1.0, 1.0 1.0))", PARSER_CHECK_NONE); poly = (LWPOLY*)lwg; @@ -624,8 +701,8 @@ void test_ptarray_point_in_ring(void) CU_ASSERT_EQUAL(result, LW_TRUE); lwgeom_free(lwg); -#if 0 /* Grazing case for point-in-polygon test, should return LW_FALSE */ + lwg = lwgeom_from_ewkt("POLYGON((2.0 3.0, 2.0 0.0, 1.0 1.0, 2.0 3.0))", PARSER_CHECK_NONE); lwg = lwgeom_from_ewkt("POLYGON((1.0 1.0, 1.0 2.0, 1.5 1.5, 1.0 1.0))", PARSER_CHECK_NONE); poly = (LWPOLY*)lwg; pt_to_test.x = 1.5; @@ -635,8 +712,28 @@ void test_ptarray_point_in_ring(void) result = ptarray_point_in_ring(poly->rings[0], pt_outside, pt_to_test); CU_ASSERT_EQUAL(result, LW_FALSE); lwgeom_free(lwg); -#endif + /* Grazing case at first point for point-in-polygon test, should return LW_FALSE */ + lwg = lwgeom_from_ewkt("POLYGON((1.0 1.0, 2.0 3.0, 2.0 0.0, 1.0 1.0))", PARSER_CHECK_NONE); + poly = (LWPOLY*)lwg; + pt_to_test.x = 1.0; + pt_to_test.y = 0.0; + pt_outside.x = 1.0; + pt_outside.y = 2.0; + result = ptarray_point_in_ring(poly->rings[0], pt_outside, pt_to_test); + CU_ASSERT_EQUAL(result, LW_FALSE); + lwgeom_free(lwg); + + /* Point on vertex of ring */ + lwg = lwgeom_from_ewkt("POLYGON((-9 50,51 -11,-10 50,-9 50))", PARSER_CHECK_NONE); + poly = (LWPOLY*)lwg; + pt_to_test.x = -10.0; + pt_to_test.y = 50.0; + pt_outside.x = -10.2727799838316134; + pt_outside.y = -16.9370033133329976; + result = ptarray_point_in_ring(poly->rings[0], pt_outside, pt_to_test); + CU_ASSERT_EQUAL(result, LW_TRUE); + lwgeom_free(lwg); #if 0 /* Small polygon and huge distance between outside point and close-but-not-quite-inside point. Should return LW_FALSE. Pretty degenerate case. */ @@ -671,7 +768,6 @@ void test_lwpoly_covers_point2d(void) result = lwpoly_covers_point2d(poly, gbox, pt_to_test); CU_ASSERT_EQUAL(result, LW_TRUE); lwgeom_free(lwg); - } diff --git a/liblwgeom/lwgeodetic.c b/liblwgeom/lwgeodetic.c index d2fde7057..a698c4483 100644 --- a/liblwgeom/lwgeodetic.c +++ b/liblwgeom/lwgeodetic.c @@ -1453,6 +1453,7 @@ int ptarray_point_in_ring(POINTARRAY *pa, POINT2D pt_outside, POINT2D pt_to_test GEOGRAPHIC_EDGE crossing_edge, edge; POINT2D p; int count = 0; + int first_point = 0; int i; /* Null input, not enough points for a ring? You ain't closed! */ @@ -1463,37 +1464,69 @@ int ptarray_point_in_ring(POINTARRAY *pa, POINT2D pt_outside, POINT2D pt_to_test geographic_point_init(pt_to_test.x, pt_to_test.y, &(crossing_edge.start)); geographic_point_init(pt_outside.x, pt_outside.y, &(crossing_edge.end)); + /* Initialize first point */ + getPoint2d_p(pa, first_point, &p); + LWDEBUGF(4, "start point == POINT(%.12g %.12g)", p.x, p.y); + geographic_point_init(p.x, p.y, &(edge.start)); + + /* If the start point is on the stab line, back up until it isn't */ + while(edge_contains_point(crossing_edge, edge.start) && ! geographic_point_equals(crossing_edge.start, edge.start) ) + { + first_point--; + LWDEBUGF(4,"first point was on stab line, reversing %d points", first_point); + getPoint2d_p(pa, pa->npoints + first_point, &p); + LWDEBUGF(4, "start point == POINT(%.12g %.12g)", p.x, p.y); + geographic_point_init(p.x, p.y, &(edge.start)); + } + /* Walk every edge and see if the stab line hits it */ for( i = 1; i < pa->npoints; i++ ) { GEOGRAPHIC_POINT g; - int cross_type; - getPoint2d_p(pa, i-1, &p); - geographic_point_init(p.x, p.y, &(edge.start)); + + LWDEBUGF(4, "start point == POINT(%.12g %.12g)", p.x, p.y); + getPoint2d_p(pa, i, &p); geographic_point_init(p.x, p.y, &(edge.end)); - /* Does stab line cross, and if so, not on the first point. We except the - first point to avoid double counting crossings at vertices. */ + LWDEBUGF(4,"testing edge (%d)", i); + + /* Our test point is on an edge! Point is "in ring" by our definition */ + if( geographic_point_equals(crossing_edge.start, edge.start) || + geographic_point_equals(crossing_edge.start, edge.end) || + edge_contains_point(edge, crossing_edge.start) ) + { + LWDEBUGF(4,"edge (%d) contains the test point, returning true", i); + return LW_TRUE; + } + + /* If the end of our edge is on the stab line, extend the edge to the + next end point, by skipping the start->end assignment step at the + end of this loop */ + if(edge_contains_point(crossing_edge, edge.end)) + { + LWDEBUGF(4,"edge (%d) end point is on the stab line, continuing", i); + continue; + } + LWDEBUG(4,"testing edge crossing"); - cross_type = edge_intersection(edge, crossing_edge, &g); - LWDEBUGF(4,"edge(%d), edge_intersection == %d", i, cross_type); - if( cross_type != LW_FALSE ) + + if( edge_intersection(edge, crossing_edge, &g) ) { - /* Don't count crossings at start points or co-linear crossings. - Start point crossings will get counted by the pre- or succeeding end points crossings. - Co-linear crossings will also get counted by the crossings of the edges pre- or succeeding them. */ - if( ! ( geographic_point_equals(g, edge.start) || cross_type == 2 ) ) - { - count++; - LWDEBUGF(4,"edge crossing is not start point or co-linear, count == %d, end_point == %d", count, geographic_point_equals(g, edge.end) ); - } - else - { - LWDEBUG(4,"edge crossing is start point, disregarding"); - } + count++; + LWDEBUGF(4,"edge (%d) crossed, count == %d", i, count); + } + else + { + LWDEBUGF(4,"edge (%d) did not cross", i); } + + /* Increment to next edge */ + edge.start = edge.end; } + + LWDEBUGF(4,"final count == %d", count); + /* An odd number of crossings implies containment! */ if( count % 2 ) {