Changeset 4787

Timestamp:

11/11/09 11:02:49 (3 years ago)

Author:

pramsey

Message:

Simplify code and improve consistency of linecrossing results (#272)

Location:

branches/1.4/liblwgeom

Files:

• cunit/cu_algorithm.c (modified) (18 diffs)
• cunit/cu_algorithm.h (modified) (1 diff)
• liblwgeom.h (modified) (1 diff)
• lwalgorithm.c (modified) (10 diffs)
• lwalgorithm.h (modified) (1 diff)

branches/1.4/liblwgeom/cunit/cu_algorithm.c

@@ -30 +30 @@
             (NULL == CU_add_test(pSuite, "test_lw_segment_intersects()", test_lw_segment_intersects)) ||
             (NULL == CU_add_test(pSuite, "test_lwline_crossing_short_lines()", test_lwline_crossing_short_lines)) ||
-            (NULL == CU_add_test(pSuite, "test_lwline_crossing_long_lines()", test_lwline_crossing_long_lines)) ||
+            (NULL == CU_add_test(pSuite, "test_lwline_crossing_long_lines()", test_lwline_crossing_long_lines)) || 
+            (NULL == CU_add_test(pSuite, "test_lwline_crossing_bugs()", test_lwline_crossing_bugs)) ||
             (NULL == CU_add_test(pSuite, "test_lwpoint_set_ordinate()", test_lwpoint_set_ordinate)) ||
             (NULL == CU_add_test(pSuite, "test_lwpoint_get_ordinate()", test_lwpoint_get_ordinate)) ||
@@ -39 +40 @@
             (NULL == CU_add_test(pSuite, "test_geohash_point()", test_geohash_point)) ||
             (NULL == CU_add_test(pSuite, "test_geohash_precision()", test_geohash_precision)) ||
-            (NULL == CU_add_test(pSuite, "test_geohash()", test_geohash))
+            (NULL == CU_add_test(pSuite, "test_geohash()", test_geohash)) 
         )
         {
@@ -52 +53 @@
 ** Global variables used by tests below
 */
-POINT2D *p1 = NULL;
-POINT2D *p2 = NULL;
-POINT2D *q1 = NULL;
-POINT2D *q2 = NULL;
-POINT4D *p = NULL;
-POINT4D *q = NULL;
+
 /* Two-point objects */
 POINTARRAY *pa21 = NULL;
@@ -63 +59 @@
 LWLINE *l21 = NULL;
 LWLINE *l22 = NULL;
-/* Five-point objects */
-LWLINE *l51 = NULL;
-LWLINE *l52 = NULL;
 /* Parsing support */
 LWGEOM_PARSER_RESULT parse_result;
@@ -76 +69 @@
 int init_cg_suite(void)
 {
-        p = lwalloc(sizeof(POINT4D));
-        q = lwalloc(sizeof(POINT4D));
-        p1 = lwalloc(sizeof(POINT2D));
-        p2 = lwalloc(sizeof(POINT2D));
-        q1 = lwalloc(sizeof(POINT2D));
-        q2 = lwalloc(sizeof(POINT2D));
         pa21 = ptarray_construct(0, 0, 2);
         pa22 = ptarray_construct(0, 0, 2);
@@ -96 +83 @@
 int clean_cg_suite(void)
 {
-        if ( p ) lwfree(p);
-        if ( p1 )lwfree(p1);
-        if ( p2 ) lwfree(p2);
-        if ( q1 ) lwfree(q1);
-        if ( q2 ) lwfree(q2);
         if ( l21 ) lwline_free(l21);
         if ( l22 ) lwline_free(l22);
@@ -112 +94 @@
 {
         double rv = 0.0;
-        POINT2D *q = NULL;
-        q = lwalloc(sizeof(POINT2D));
+        POINT2D p1, p2, q;
 
         /* Vertical line at x=0 */
-        p1->x = 0.0;
-        p1->y = 0.0;
-        p2->x = 0.0;
-        p2->y = 1.0;
+        p1.x = 0.0;
+        p1.y = 0.0;
+        p2.x = 0.0;
+        p2.y = 1.0;
 
         /* On the left */
-        q->x = -2.0;
-        q->y = 1.5;
+        q.x = -2.0;
+        q.y = 1.5;
         rv = lw_segment_side(p1, p2, q);
         //printf("left %g\n",rv);
@@ -129 +110 @@
 
         /* On the right */
-        q->x = 2.0;
+        q.x = 2.0;
         rv = lw_segment_side(p1, p2, q);
         //printf("right %g\n",rv);
@@ -135 +116 @@
 
         /* On the line */
-        q->x = 0.0;
+        q.x = 0.0;
         rv = lw_segment_side(p1, p2, q);
         //printf("on line %g\n",rv);
         CU_ASSERT_EQUAL(rv, 0.0);
-
-        lwfree(q);
 
 }
@@ -150 +129 @@
 {
 
+#define setpoint(p, x1, y1) {(p).x = (x1); (p).y = (y1);}
+
+        POINT2D p1, p2, q1, q2;
+        
         /* P: Vertical line at x=0 */
-        p1->x = 0.0;
-        p1->y = 0.0;
-        p2->x = 0.0;
-        p2->y = 1.0;
+        setpoint(p1, 0.0, 0.0);
+        p1.x = 0.0;
+        p1.y = 0.0;
+        p2.x = 0.0;
+        p2.y = 1.0;
 
         /* Q: Horizontal line crossing left to right */
-        q1->x = -0.5;
-        q1->y = 0.5;
-        q2->x = 0.5;
-        q2->y = 0.5;
+        q1.x = -0.5;
+        q1.y = 0.5;
+        q2.x = 0.5;
+        q2.y = 0.5;
         CU_ASSERT( lw_segment_intersects(p1, p2, q1, q2) == SEG_CROSS_RIGHT );
 
         /* Q: Horizontal line crossing right to left */
-        q1->x = 0.5;
-        q1->y = 0.5;
-        q2->x = -0.5;
-        q2->y = 0.5;
+        q1.x = 0.5;
+        q1.y = 0.5;
+        q2.x = -0.5;
+        q2.y = 0.5;
         CU_ASSERT( lw_segment_intersects(p1, p2, q1, q2) == SEG_CROSS_LEFT );
 
         /* Q: Horizontal line not crossing right to left */
-        q1->x = 0.5;
-        q1->y = 1.5;
-        q2->x = -0.5;
-        q2->y = 1.5;
+        q1.x = 0.5;
+        q1.y = 1.5;
+        q2.x = -0.5;
+        q2.y = 1.5;
         CU_ASSERT( lw_segment_intersects(p1, p2, q1, q2) == SEG_NO_INTERSECTION );
 
-        /* Q: Horizontal line crossing at vertex right to left */
-        q1->x = 0.5;
-        q1->y = 1.0;
-        q2->x = -0.5;
-        q2->y = 1.0;
+        /* Q: Horizontal line crossing at second vertex right to left */
+        q1.x = 0.5;
+        q1.y = 1.0;
+        q2.x = -0.5;
+        q2.y = 1.0;
+        CU_ASSERT( lw_segment_intersects(p1, p2, q1, q2) == SEG_NO_INTERSECTION );
+
+        /* Q: Horizontal line crossing at first vertex right to left */
+        q1.x = 0.5;
+        q1.y = 0.0;
+        q2.x = -0.5;
+        q2.y = 0.0;
         CU_ASSERT( lw_segment_intersects(p1, p2, q1, q2) == SEG_CROSS_LEFT );
 
-        /* Q: Diagonal line with large range crossing at vertex right to left */
-        q1->x = 0.5;
-        q1->y = 10.0;
-        q2->x = -0.5;
-        q2->y = -10.0;
+        /* Q: Diagonal line with large range crossing at first vertex right to left */
+        q1.x = 0.5;
+        q1.y = 10.0;
+        q2.x = -0.5;
+        q2.y = -10.0;
         CU_ASSERT( lw_segment_intersects(p1, p2, q1, q2) == SEG_CROSS_LEFT );
 
-        /* Q: Diagonal line with large range crossing at vertex right to left */
-        q1->x = 0.5;
-        q1->y = 11.0;
-        q2->x = -0.5;
-        q2->y = -9.0;
+        /* Q: Diagonal line with large range crossing at second vertex right to left */
+        q1.x = 0.5;
+        q1.y = 11.0;
+        q2.x = -0.5;
+        q2.y = -9.0;
+        CU_ASSERT( lw_segment_intersects(p1, p2, q1, q2) == SEG_NO_INTERSECTION );
+
+        /* Q: Horizontal touching from left at second vertex*/
+        q1.x = -0.5;
+        q1.y = 0.5;
+        q2.x = 0.0;
+        q2.y = 0.5;
+        CU_ASSERT( lw_segment_intersects(p1, p2, q1, q2) == SEG_NO_INTERSECTION );
+
+        /* Q: Horizontal touching from right at first vertex */
+        q1.x = 0.0;
+        q1.y = 0.5;
+        q2.x = 0.5;
+        q2.y = 0.5;
+        CU_ASSERT( lw_segment_intersects(p1, p2, q1, q2) == SEG_CROSS_RIGHT );
+
+        /* Q: Horizontal touching from left and far below on second vertex */
+        q1.x = -0.5;
+        q1.y = -10.5;
+        q2.x = 0.0;
+        q2.y = 0.5;
+        CU_ASSERT( lw_segment_intersects(p1, p2, q1, q2) == SEG_NO_INTERSECTION );
+
+        /* Q: Horizontal touching from right and far above on second vertex */
+        q1.x = 0.5;
+        q1.y = 10.5;
+        q2.x = 0.0;
+        q2.y = 0.5;
+        CU_ASSERT( lw_segment_intersects(p1, p2, q1, q2) == SEG_NO_INTERSECTION );
+
+        /* Q: Co-linear from top */
+        q1.x = 0.0;
+        q1.y = 10.0;
+        q2.x = 0.0;
+        q2.y = 0.5;
+        CU_ASSERT( lw_segment_intersects(p1, p2, q1, q2) == SEG_COLINEAR );
+
+        /* Q: Co-linear from bottom */
+        q1.x = 0.0;
+        q1.y = -10.0;
+        q2.x = 0.0;
+        q2.y = 0.5;
+        CU_ASSERT( lw_segment_intersects(p1, p2, q1, q2) == SEG_COLINEAR );
+
+        /* Q: Co-linear contained */
+        q1.x = 0.0;
+        q1.y = 0.4;
+        q2.x = 0.0;
+        q2.y = 0.5;
+        CU_ASSERT( lw_segment_intersects(p1, p2, q1, q2) == SEG_COLINEAR );
+
+        /* Q: Horizontal touching at end point from left */
+        q1.x = -0.5;
+        q1.y = 1.0;
+        q2.x = 0.0;
+        q2.y = 1.0;
+        CU_ASSERT( lw_segment_intersects(p1, p2, q1, q2) == SEG_NO_INTERSECTION );
+
+        /* Q: Horizontal touching at end point from right */
+        q1.x = 0.0;
+        q1.y = 1.0;
+        q2.x = 0.0;
+        q2.y = 0.5;
+        CU_ASSERT( lw_segment_intersects(p1, p2, q1, q2) == SEG_COLINEAR );
+
+        /* Q: Horizontal touching at start point from left */
+        q1.x = 0.0;
+        q1.y = 0.0;
+        q2.x = -0.5;
+        q2.y = 0.0;
         CU_ASSERT( lw_segment_intersects(p1, p2, q1, q2) == SEG_CROSS_LEFT );
 
-        /* Q: Horizontal touching from left */
-        q1->x = -0.5;
-        q1->y = 0.5;
-        q2->x = 0.0;
-        q2->y = 0.5;
-        CU_ASSERT( lw_segment_intersects(p1, p2, q1, q2) == SEG_TOUCH_LEFT );
-
-        /* Q: Horizontal touching from right */
-        q1->x = 0.5;
-        q1->y = 0.5;
-        q2->x = 0.0;
-        q2->y = 0.5;
-        CU_ASSERT( lw_segment_intersects(p1, p2, q1, q2) == SEG_TOUCH_RIGHT );
-
-        /* Q: Horizontal touching from left and far below*/
-        q1->x = -0.5;
-        q1->y = -10.5;
-        q2->x = 0.0;
-        q2->y = 0.5;
-        CU_ASSERT( lw_segment_intersects(p1, p2, q1, q2) == SEG_TOUCH_LEFT );
-
-        /* Q: Horizontal touching from right and far above */
-        q1->x = 0.5;
-        q1->y = 10.5;
-        q2->x = 0.0;
-        q2->y = 0.5;
-        CU_ASSERT( lw_segment_intersects(p1, p2, q1, q2) == SEG_TOUCH_RIGHT );
-
-        /* Q: Co-linear from top */
-        q1->x = 0.0;
-        q1->y = 10.0;
-        q2->x = 0.0;
-        q2->y = 0.5;
-        CU_ASSERT( lw_segment_intersects(p1, p2, q1, q2) == SEG_COLINEAR );
-
-        /* Q: Co-linear from bottom */
-        q1->x = 0.0;
-        q1->y = -10.0;
-        q2->x = 0.0;
-        q2->y = 0.5;
-        CU_ASSERT( lw_segment_intersects(p1, p2, q1, q2) == SEG_COLINEAR );
-
-        /* Q: Co-linear contained */
-        q1->x = 0.0;
-        q1->y = 0.4;
-        q2->x = 0.0;
-        q2->y = 0.5;
-        CU_ASSERT( lw_segment_intersects(p1, p2, q1, q2) == SEG_COLINEAR );
-
-        /* Q: Horizontal touching at end point from left */
-        q1->x = -0.5;
-        q1->y = 1.0;
-        q2->x = 0.0;
-        q2->y = 1.0;
-        CU_ASSERT( lw_segment_intersects(p1, p2, q1, q2) == SEG_TOUCH_LEFT );
-
-        /* Q: Horizontal touching at end point from right */
-        q1->x = 0.5;
-        q1->y = 1.0;
-        q2->x = 0.0;
-        q2->y = 1.0;
-        CU_ASSERT( lw_segment_intersects(p1, p2, q1, q2) == SEG_TOUCH_RIGHT );
-
-        /* Q: Horizontal touching at start point from left */
-        q1->x = -0.5;
-        q1->y = 0.0;
-        q2->x = 0.0;
-        q2->y = 0.0;
-        CU_ASSERT( lw_segment_intersects(p1, p2, q1, q2) == SEG_TOUCH_LEFT );
-
         /* Q: Horizontal touching at start point from right */
-        q1->x = 0.5;
-        q1->y = 0.0;
-        q2->x = 0.0;
-        q2->y = 0.0;
-        CU_ASSERT( lw_segment_intersects(p1, p2, q1, q2) == SEG_TOUCH_RIGHT );
+        q1.x = 0.0;
+        q1.y = 0.0;
+        q2.x = 0.5;
+        q2.y = 0.0;
+        CU_ASSERT( lw_segment_intersects(p1, p2, q1, q2) == SEG_CROSS_RIGHT );
 
 }
@@ -280 +271 @@
 {
 
+        POINT4D p;
+        
         /*
         ** Simple test, two two-point lines 
@@ -285 +278 @@
 
         /* Vertical line from 0,0 to 1,1 */
-        p->x = 0.0;
-        p->y = 0.0;
-        setPoint4d(pa21, 0, p);
-        p->y = 1.0;
-        setPoint4d(pa21, 1, p);
+        p.x = 0.0;
+        p.y = 0.0;
+        setPoint4d(pa21, 0, &p);
+        p.y = 1.0;
+        setPoint4d(pa21, 1, &p);
 
         /* Horizontal, crossing mid-segment */
-        p->x = -0.5;
-        p->y = 0.5;
-        setPoint4d(pa22, 0, p);
-        p->x = 0.5;
-        setPoint4d(pa22, 1, p);
+        p.x = -0.5;
+        p.y = 0.5;
+        setPoint4d(pa22, 0, &p);
+        p.x = 0.5;
+        setPoint4d(pa22, 1, &p);
 
         CU_ASSERT( lwline_crossing_direction(l21, l22) == LINE_CROSS_RIGHT );
 
-        /* Horizontal, crossing at top end vertex */
-        p->x = -0.5;
-        p->y = 1.0;
-        setPoint4d(pa22, 0, p);
-        p->x = 0.5;
-        setPoint4d(pa22, 1, p);
+        /* Horizontal, crossing at top end vertex (end crossings don't count) */
+        p.x = -0.5;
+        p.y = 1.0;
+        setPoint4d(pa22, 0, &p);
+        p.x = 0.5;
+        setPoint4d(pa22, 1, &p);
+
+        CU_ASSERT( lwline_crossing_direction(l21, l22) == LINE_NO_CROSS );
+
+        /* Horizontal, crossing at bottom end vertex */
+        p.x = -0.5;
+        p.y = 0.0;
+        setPoint4d(pa22, 0, &p);
+        p.x = 0.5;
+        setPoint4d(pa22, 1, &p);
 
         CU_ASSERT( lwline_crossing_direction(l21, l22) == LINE_CROSS_RIGHT );
 
-        /* Horizontal, crossing at bottom end vertex */
-        p->x = -0.5;
-        p->y = 0.0;
-        setPoint4d(pa22, 0, p);
-        p->x = 0.5;
-        setPoint4d(pa22, 1, p);
-
-        CU_ASSERT( lwline_crossing_direction(l21, l22) == LINE_CROSS_RIGHT );
-
         /* Horizontal, no crossing */
-        p->x = -0.5;
-        p->y = 2.0;
-        setPoint4d(pa22, 0, p);
-        p->x = 0.5;
-        setPoint4d(pa22, 1, p);
+        p.x = -0.5;
+        p.y = 2.0;
+        setPoint4d(pa22, 0, &p);
+        p.x = 0.5;
+        setPoint4d(pa22, 1, &p);
 
         CU_ASSERT( lwline_crossing_direction(l21, l22) == LINE_NO_CROSS );
 
         /* Vertical, no crossing */
-        p->x = -0.5;
-        p->y = 0.0;
-        setPoint4d(pa22, 0, p);
-        p->y = 1.0;
-        setPoint4d(pa22, 1, p);
+        p.x = -0.5;
+        p.y = 0.0;
+        setPoint4d(pa22, 0, &p);
+        p.y = 1.0;
+        setPoint4d(pa22, 1, &p);
 
         CU_ASSERT( lwline_crossing_direction(l21, l22) == LINE_NO_CROSS );
@@ -345 +338 @@
         ** More complex test, longer lines and multiple crossings 
         */
-
         /* Vertical line with vertices at y integers */
         l51 = (LWLINE*)lwgeom_from_ewkt("LINESTRING(0 0, 0 1, 0 2, 0 3, 0 4)", PARSER_CHECK_NONE);
@@ -374 +366 @@
         lwline_free(l52);
 
-        /* Two crossings, one at the last vertex on the next interior vertex */
-        l52 = (LWLINE*)lwgeom_from_ewkt("LINESTRING(1 4, 0 4, -1 3, 0 3, 1 3)", PARSER_CHECK_NONE);
-        CU_ASSERT( lwline_crossing_direction(l51, l52) == LINE_MULTICROSS_END_SAME_FIRST_LEFT );
-        lwline_free(l52);
-
         /* Three crossings, two at midpoints, one at vertex */
         l52 = (LWLINE*)lwgeom_from_ewkt("LINESTRING(0.5 1, -1 0.5, 1 2, -1 2, -1 3)", PARSER_CHECK_NONE);
@@ -403 +390 @@
 }
 
+
+void test_lwline_crossing_bugs(void)
+{
+        LWLINE *l1;
+        LWLINE *l2;
+        
+        l1 = (LWLINE*)lwgeom_from_ewkt("LINESTRING(2.99 90.16,71 74,20 140,171 154)", PARSER_CHECK_NONE);
+        l2 = (LWLINE*)lwgeom_from_ewkt("LINESTRING(25 169,89 114,40 70,86 43)", PARSER_CHECK_NONE);
+
+        CU_ASSERT( lwline_crossing_direction(l1, l2) == LINE_MULTICROSS_END_RIGHT );
+        lwline_free(l1);
+        lwline_free(l2);
+        
+}
+
 void test_lwpoint_set_ordinate(void)
 {
-        p->x = 0.0;
-        p->y = 0.0;
-        p->z = 0.0;
-        p->m = 0.0;
-
-        lwpoint_set_ordinate(p, 0, 1.5);
-        CU_ASSERT_EQUAL( p->x, 1.5 );
-
-        lwpoint_set_ordinate(p, 3, 2.5);
-        CU_ASSERT_EQUAL( p->m, 2.5 );
-
-        lwpoint_set_ordinate(p, 2, 3.5);
-        CU_ASSERT_EQUAL( p->z, 3.5 );
+        POINT4D p;
+
+        p.x = 0.0;
+        p.y = 0.0;
+        p.z = 0.0;
+        p.m = 0.0;
+
+        lwpoint_set_ordinate(&p, 0, 1.5);
+        CU_ASSERT_EQUAL( p.x, 1.5 );
+
+        lwpoint_set_ordinate(&p, 3, 2.5);
+        CU_ASSERT_EQUAL( p.m, 2.5 );
+
+        lwpoint_set_ordinate(&p, 2, 3.5);
+        CU_ASSERT_EQUAL( p.z, 3.5 );
 
 }
@@ -423 +427 @@
 void test_lwpoint_get_ordinate(void)
 {
-
-        p->x = 10.0;
-        p->y = 20.0;
-        p->z = 30.0;
-        p->m = 40.0;
-
-        CU_ASSERT_EQUAL( lwpoint_get_ordinate(p, 0), 10.0 );
-        CU_ASSERT_EQUAL( lwpoint_get_ordinate(p, 1), 20.0 );
-        CU_ASSERT_EQUAL( lwpoint_get_ordinate(p, 2), 30.0 );
-        CU_ASSERT_EQUAL( lwpoint_get_ordinate(p, 3), 40.0 );
+        POINT4D p;
+
+        p.x = 10.0;
+        p.y = 20.0;
+        p.z = 30.0;
+        p.m = 40.0;
+
+        CU_ASSERT_EQUAL( lwpoint_get_ordinate(&p, 0), 10.0 );
+        CU_ASSERT_EQUAL( lwpoint_get_ordinate(&p, 1), 20.0 );
+        CU_ASSERT_EQUAL( lwpoint_get_ordinate(&p, 2), 30.0 );
+        CU_ASSERT_EQUAL( lwpoint_get_ordinate(&p, 3), 40.0 );
 
 }
@@ -438 +443 @@
 void test_lwpoint_interpolate(void)
 {
-        POINT4D *r = NULL;
-        r = lwalloc(sizeof(POINT4D));
+        POINT4D p, q, r;
         int rv = 0;
 
-        p->x = 10.0;
-        p->y = 20.0;
-        p->z = 30.0;
-        p->m = 40.0;
-        q->x = 20.0;
-        q->y = 30.0;
-        q->z = 40.0;
-        q->m = 50.0;
-
-        rv = lwpoint_interpolate(p, q, r, 4, 2, 35.0);
-        CU_ASSERT_EQUAL( r->x, 15.0);
-
-        rv = lwpoint_interpolate(p, q, r, 4, 3, 41.0);
-        CU_ASSERT_EQUAL( r->y, 21.0);
-
-        rv = lwpoint_interpolate(p, q, r, 4, 3, 50.0);
-        CU_ASSERT_EQUAL( r->y, 30.0);
-
-        rv = lwpoint_interpolate(p, q, r, 4, 3, 40.0);
-        CU_ASSERT_EQUAL( r->y, 20.0);
-
-        lwfree(r);
+        p.x = 10.0;
+        p.y = 20.0;
+        p.z = 30.0;
+        p.m = 40.0;
+        
+        q.x = 20.0;
+        q.y = 30.0;
+        q.z = 40.0;
+        q.m = 50.0;
+
+        rv = lwpoint_interpolate(&p, &q, &r, 4, 2, 35.0);
+        CU_ASSERT_EQUAL( r.x, 15.0);
+
+        rv = lwpoint_interpolate(&p, &q, &r, 4, 3, 41.0);
+        CU_ASSERT_EQUAL( r.y, 21.0);
+
+        rv = lwpoint_interpolate(&p, &q, &r, 4, 3, 50.0);
+        CU_ASSERT_EQUAL( r.y, 30.0);
+
+        rv = lwpoint_interpolate(&p, &q, &r, 4, 3, 40.0);
+        CU_ASSERT_EQUAL( r.y, 20.0);
 
 }
@@ -679 +682 @@
         LWCOLLECTION *c;
         char *ewkt;
-
-        p->x = 0.0;
-        p->y = 0.0;
-        p->z = 0.0;
-        setPoint4d(pa, 0, p);
-
-        p->x = 1.0;
-        p->y = 1.0;
-        p->z = 1.0;
-        setPoint4d(pa, 1, p);
-
-        p->x = 2.0;
-        p->y = 2.0;
-        p->z = 2.0;
-        setPoint4d(pa, 2, p);
+        POINT4D p;
+
+        p.x = 0.0;
+        p.y = 0.0;
+        p.z = 0.0;
+        setPoint4d(pa, 0, &p);
+
+        p.x = 1.0;
+        p.y = 1.0;
+        p.z = 1.0;
+        setPoint4d(pa, 1, &p);
+
+        p.x = 2.0;
+        p.y = 2.0;
+        p.z = 2.0;
+        setPoint4d(pa, 2, &p);
 
         c = lwline_clip_to_ordinate_range(line, 2, 0.5, 1.5);

branches/1.4/liblwgeom/cunit/cu_algorithm.h

@@ -41 +41 @@
 void test_geohash_point(void);
 void test_geohash(void);
+void test_lwline_crossing_bugs(void);

branches/1.4/liblwgeom/liblwgeom.h

@@ -54 +54 @@
 #define FP_CONTAINS_EXCL(A, X, B) (FP_LT(A, X) && FP_LT(X, B))
 #define FP_CONTAINS(A, X, B) FP_CONTAINS_EXCL(A, X, B)
+#define FP_IS_ZERO(A) (fabs(A) <= PGIS_EPSILON)
 #define LW_TRUE 1
 #define LW_FALSE 0

branches/1.4/liblwgeom/lwalgorithm.c

@@ -21 +21 @@
 ** Return = 0.0 if point Q in on segment P
 */
-double lw_segment_side(POINT2D *p1, POINT2D *p2, POINT2D *q)
-{
-        return ( (q->x - p1->x) * (p2->y - p1->y) - (p2->x - p1->x) * (q->y - p1->y) );
+double lw_segment_side(POINT2D p1, POINT2D p2, POINT2D q)
+{
+        return ( (q.x - p1.x) * (p2.y - p1.y) - (p2.x - p1.x) * (q.y - p1.y) );
 }
 
@@ -33 +33 @@
         double maxp=LW_MAX(p1.x,p2.x);
 
-        if (minp>maxq || maxp<minq)
+        if (FP_GT(minp,maxq) || FP_LT(maxp,minq))
                 return LW_FALSE;
 
@@ -41 +41 @@
         maxp=LW_MAX(p1.y,p2.y);
 
-        if (minp>maxq || maxp<minq)
+        if (FP_GT(minp,maxq) || FP_LT(maxp,minq))
                 return LW_FALSE;
 
         return LW_TRUE;
 }
+
 
 /**
@@ -60 +61 @@
 **              SEG_CROSS_LEFT = 2,
 **              SEG_CROSS_RIGHT = 3,
-**              SEG_TOUCH_LEFT = 4,
-**              SEG_TOUCH_RIGHT = 5
 */
-int lw_segment_intersects(POINT2D *p1, POINT2D *p2, POINT2D *q1, POINT2D *q2)
+int lw_segment_intersects(POINT2D p1, POINT2D p2, POINT2D q1, POINT2D q2)
 {
 
@@ -69 +68 @@
 
         /* No envelope interaction => we are done. */
-        if (!lw_segment_envelope_intersects(*p1, *p2, *q1, *p2))
+        if (!lw_segment_envelope_intersects(p1, p2, q1, p2))
         {
                 return SEG_NO_INTERSECTION;
@@ -91 +90 @@
 
         /* Nobody is on one side or another? Must be colinear. */
-        if (pq1 == 0.0 && pq2 == 0.0 && qp1 == 0.0 && qp2 == 0.0)
+        if (FP_IS_ZERO(pq1) && FP_IS_ZERO(pq2) && FP_IS_ZERO(qp1) && FP_IS_ZERO(qp2))
         {
                 return SEG_COLINEAR;
@@ -98 +97 @@
         /*
         ** When one end-point touches, the sidedness is determined by the
-        ** location of the other end-point.
+        ** location of the other end-point. Only touches by the first point
+        ** will be considered "real" to avoid double counting.
         */
-        if ( pq2 == 0.0 )
-        {
-                if ( pq1 < 0.0 )
-                        return SEG_TOUCH_LEFT;
+        LWDEBUGF(4, "pq1=%.15g pq2=%.15g", pq1, pq2);
+        LWDEBUGF(4, "qp1=%.15g qp2=%.15g", qp1, qp2);
+
+        /* Second point of p or q touches, it's not a crossing. */
+        if ( FP_IS_ZERO(pq2) || FP_IS_ZERO(qp2) )
+        {
+                return SEG_NO_INTERSECTION;
+        }
+
+        /* First point of p touches, it's a "crossing". */
+        if ( FP_IS_ZERO(pq1) )
+        {
+                if ( FP_GT(pq2,0.0) )
+                        return SEG_CROSS_RIGHT;
                 else
-                        return SEG_TOUCH_RIGHT;
-        }
-        if ( pq1 == 0.0 )
-        {
-                if ( pq2 < 0.0 )
-                        return SEG_TOUCH_LEFT;
+                        return SEG_CROSS_LEFT;
+        }
+
+        /* First point of q touches, it's a crossing. */
+        if ( FP_IS_ZERO(qp1) )
+        {
+                if ( FP_LT(pq1,pq2) )
+                        return SEG_CROSS_RIGHT;
                 else
-                        return SEG_TOUCH_RIGHT;
-        }
-
+                        return SEG_CROSS_LEFT;
+        }
+        
         /* The segments cross, what direction is the crossing? */
-        if ( pq1 < pq2 )
+        if ( FP_LT(pq1,pq2) )
                 return SEG_CROSS_RIGHT;
         else
@@ -142 +154 @@
 int lwline_crossing_direction(LWLINE *l1, LWLINE *l2)
 {
-
         int i = 0, j = 0, rv = 0;
-        POINT2D *p1;
-        POINT2D *p2;
-        POINT2D *q1;
-        POINT2D *q2;
-        POINTARRAY *pa1;
-        POINTARRAY *pa2;
+        POINT2D p1, p2, q1, q2;
+        POINTARRAY *pa1 = NULL, *pa2 = NULL;
         int cross_left = 0;
         int cross_right = 0;
         int first_cross = 0;
-        int final_cross = 0;
         int this_cross = 0;
-        int vertex_touch = -1;
-        int vertex_touch_type = 0;
 
         pa1 = (POINTARRAY*)l1->points;
         pa2 = (POINTARRAY*)l2->points;
-
-        p1 = lwalloc(sizeof(POINT2D));
-        p2 = lwalloc(sizeof(POINT2D));
-        q1 = lwalloc(sizeof(POINT2D));
-        q2 = lwalloc(sizeof(POINT2D));
 
         /* One-point lines can't intersect (and shouldn't exist). */
@@ -170 +169 @@
                 return LINE_NO_CROSS;
 
-        LWDEBUGF(4, "lineCrossingDirection: l1 = %s", lwgeom_to_ewkt((LWGEOM*)l1,0));
-        LWDEBUGF(4, "lineCrossingDirection: l2 = %s", lwgeom_to_ewkt((LWGEOM*)l2,0));
+        LWDEBUGF(4, "l1 = %s", lwgeom_to_ewkt((LWGEOM*)l1,0));
+        LWDEBUGF(4, "l2 = %s", lwgeom_to_ewkt((LWGEOM*)l2,0));
+
+        /* Initialize first point of q */
+        rv = getPoint2d_p(pa2, 0, &q1);
 
         for ( i = 1; i < pa2->npoints; i++ )
         {
 
-                rv = getPoint2d_p(pa2, i-1, q1);
-                rv = getPoint2d_p(pa2, i, q2);
-
+                /* Update second point of q to next value */
+                rv = getPoint2d_p(pa2, i, &q2);
+
+                /* Initialize first point of p */
+                rv = getPoint2d_p(pa1, 0, &p1);
+                
                 for ( j = 1; j < pa1->npoints; j++ )
                 {
 
-                        rv = getPoint2d_p(pa1, j-1, p1);
-                        rv = getPoint2d_p(pa1, j, p2);
-
-                        LWDEBUGF(4, "lineCrossingDirection: i=%d, j=%d", i, j);
-
+                        /* Update second point of p to next value */
+                        rv = getPoint2d_p(pa1, j, &p2);
+                        
                         this_cross = lw_segment_intersects(p1, p2, q1, q2);
 
-                        if ( ! first_cross && this_cross )
-                                first_cross = this_cross;
-                        if ( this_cross )
-                                final_cross = this_cross;
+                        LWDEBUGF(4, "i=%d, j=%d (%.8g %.8g, %.8g %.8g)", this_cross, i, j, p1.x, p1.y, p2.x, p2.y);
 
                         if ( this_cross == SEG_CROSS_LEFT )
                         {
+                                LWDEBUG(4,"this_cross == SEG_CROSS_LEFT");
                                 cross_left++;
-                                break;
+                                if ( ! first_cross ) 
+                                        first_cross = SEG_CROSS_LEFT;
                         }
 
                         if ( this_cross == SEG_CROSS_RIGHT )
                         {
+                                LWDEBUG(4,"this_cross == SEG_CROSS_RIGHT");
                                 cross_right++;
-                                break;
+                                if ( ! first_cross ) 
+                                        first_cross = SEG_CROSS_LEFT;
                         }
 
                         /*
-                        ** Crossing at a co-linearity can be turned into crossing at
-                        ** a vertex by pulling the original touch point forward along
-                        ** the co-linearity.
+                        ** Crossing at a co-linearity can be turned handled by extending
+                        ** segment to next vertext and seeing if the end points straddle
+                        ** the co-linear segment.
                         */
-                        if ( this_cross == SEG_COLINEAR && vertex_touch == (i-1) )
-                        {
-                                vertex_touch = i;
-                                break;
-                        }
-
-                        /*
-                        ** Crossing-at-vertex will cause four segment touch interactions, two at
-                        ** j-1 and two at j. We avoid incrementing our cross count by ignoring the
-                        ** second pair.
-                        */
-                        if ( this_cross == SEG_TOUCH_LEFT )
-                        {
-                                if ( vertex_touch == (i-1) && vertex_touch_type == SEG_TOUCH_RIGHT )
-                                {
-                                        cross_left++;
-                                        vertex_touch = -1;
-                                        vertex_touch_type = 0;
-                                }
-                                else
-                                {
-                                        vertex_touch = i;
-                                        vertex_touch_type = this_cross;
-                                }
-                                break;
-                        }
-                        if ( this_cross == SEG_TOUCH_RIGHT )
-                        {
-                                if ( vertex_touch == (i-1) && vertex_touch_type == SEG_TOUCH_LEFT )
-                                {
-                                        cross_right++;
-                                        vertex_touch = -1;
-                                        vertex_touch_type = 0;
-                                }
-                                else
-                                {
-                                        vertex_touch = i;
-                                        vertex_touch_type = this_cross;
-                                }
-                                break;
-                        }
-
-                        /*
-                        ** TODO Handle co-linear cases.
-                        */
-
-                        LWDEBUGF(4, "lineCrossingDirection: this_cross=%d, vertex_touch=%d, vertex_touch_type=%d", this_cross, vertex_touch, vertex_touch_type);
-
-                }
-
-        }
-
-        LWDEBUGF(4, "first_cross=%d, final_cross=%d, cross_left=%d, cross_right=%d", first_cross, final_cross, cross_left, cross_right);
-
-        lwfree(p1);
-        lwfree(p2);
-        lwfree(q1);
-        lwfree(q2);
+                        if ( this_cross == SEG_COLINEAR )
+                        {
+                                LWDEBUG(4,"this_cross == SEG_COLINEAR");
+                                /* TODO: Add logic here and in segment_intersects() 
+                                continue; 
+                                */
+                        }
+                        
+                        LWDEBUG(4,"this_cross == SEG_NO_INTERSECTION");
+                        
+                        /* Turn second point of p into first point */
+                        p1 = p2;
+
+                }
+                
+                /* Turn second point of q into first point */
+                q1 = q2;
+
+        }
+
+        LWDEBUGF(4, "first_cross=%d, cross_left=%d, cross_right=%d", first_cross, cross_left, cross_right);
 
         if ( !cross_left && !cross_right )
@@ -291 +258 @@
                 return LINE_MULTICROSS_END_SAME_FIRST_RIGHT;
 
-        if ( cross_left - cross_right == 0 && first_cross == SEG_TOUCH_LEFT )
-                return LINE_MULTICROSS_END_SAME_FIRST_RIGHT;
-
-        if ( cross_left - cross_right == 0 && first_cross == SEG_TOUCH_RIGHT )
-                return LINE_MULTICROSS_END_SAME_FIRST_LEFT;
-
         return LINE_NO_CROSS;
 
 }
+
+
 
 /*

branches/1.4/liblwgeom/lwalgorithm.h

@@ -24 +24 @@
 };
 
-double lw_segment_side(POINT2D *p1, POINT2D *p2, POINT2D *q);
-int lw_segment_intersects(POINT2D *p1, POINT2D *p2, POINT2D *q1, POINT2D *q2);
+double lw_segment_side(POINT2D p1, POINT2D p2, POINT2D q);
+int lw_segment_intersects(POINT2D p1, POINT2D p2, POINT2D q1, POINT2D q2);
 int lw_segment_envelope_intersects(POINT2D p1, POINT2D p2, POINT2D q1, POINT2D q2);