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 469aad704097a356b457b9debc5c5ee23114a5f3..09eca8597a6097a4b5c7fded8a35631ef1a11bf6 100644 (file)
--- 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 d2fde7057e4ac131ea916fd0343b23ce89cdb05f..a698c4483049729a3ea5042fd5303e43eec7f6b4 100644 (file)
--- 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 )
        {