]> granicus.if.org Git - postgis/commitdiff
First cut of line crossing function, and associated cunit tests.
authorPaul Ramsey <pramsey@cleverelephant.ca>
Tue, 16 Dec 2008 22:29:42 +0000 (22:29 +0000)
committerPaul Ramsey <pramsey@cleverelephant.ca>
Tue, 16 Dec 2008 22:29:42 +0000 (22:29 +0000)
git-svn-id: http://svn.osgeo.org/postgis/trunk@3424 b70326c6-7e19-0410-871a-916f4a2858ee

cunit/Makefile [new file with mode: 0644]
cunit/cu_algorithm.c [new file with mode: 0644]
cunit/cu_algorithm.h [new file with mode: 0644]
cunit/cu_tester.c [new file with mode: 0644]
liblwgeom/Makefile.in
liblwgeom/liblwgeom.h
liblwgeom/lwalgorithm.c [new file with mode: 0644]
liblwgeom/lwalgorithm.h [new file with mode: 0644]

diff --git a/cunit/Makefile b/cunit/Makefile
new file mode 100644 (file)
index 0000000..4e08915
--- /dev/null
@@ -0,0 +1,25 @@
+
+LWGEOM=../liblwgeom
+LIBCUNIT=/usr/local/lib/libcunit.a
+LIBLWGEOM=$(LWGEOM)/liblwgeom.a
+
+OBJS = cu_tester.o \
+       cu_algorithm.o
+
+#---------------------------------------------------------------
+
+INCLUDES = -I/usr/local/include -I$(LWGEOM)
+CFLAGS = -g $(MSDEFINES) $(INCLUDES)
+LIBS = $(LIBCUNIT) $(LIBLWGEOM)
+LDFLAGS = 
+
+cu_tester: $(OBJS) $(LIBS)
+
+test: cu_tester
+       ./cu_tester
+
+check: test
+
+clean:
+       @rm $(OBJS)
+       @rm cu_tester
diff --git a/cunit/cu_algorithm.c b/cunit/cu_algorithm.c
new file mode 100644 (file)
index 0000000..10a2929
--- /dev/null
@@ -0,0 +1,476 @@
+
+#include "cu_algorithm.h"
+
+/*
+** Called from test harness to register the tests in this file.
+*/
+CU_pSuite register_cg_suite(void)
+{
+       CU_pSuite pSuite;
+       pSuite = CU_add_suite("PostGIS Computational Geometry Suite", init_cg_suite, clean_cg_suite);
+       if (NULL == pSuite)
+       {
+               CU_cleanup_registry();
+               return NULL;
+       }
+
+       if (
+           (NULL == CU_add_test(pSuite, "testSegmentSide()", testSegmentSide)) ||
+           (NULL == CU_add_test(pSuite, "testSegmentIntersects()", testSegmentIntersects)) ||
+           (NULL == CU_add_test(pSuite, "testLineCrossingShortLines()", testLineCrossingShortLines)) ||
+           (NULL == CU_add_test(pSuite, "testLineCrossingLongLines()", testLineCrossingLongLines)) 
+       )
+       {
+               CU_cleanup_registry();
+               return NULL;
+       }
+       return pSuite;
+}
+
+
+/*
+** Global variables used by tests below
+*/
+POINT2D *p1 = NULL;
+POINT2D *p2 = NULL;
+POINT2D *q1 = NULL;
+POINT2D *q2 = NULL;
+POINT4D *p = NULL;
+/* Two-point objects */
+POINTARRAY *pa21 = NULL;
+POINTARRAY *pa22 = NULL;
+LWLINE *l21 = NULL;
+LWLINE *l22 = NULL;
+/* Five-point objects */
+POINTARRAY *pa51 = NULL;
+POINTARRAY *pa52 = NULL;
+LWLINE *l51 = NULL;
+LWLINE *l52 = NULL;
+
+
+/*
+** The suite initialization function.
+** Create any re-used objects.
+*/
+int init_cg_suite(void)
+{
+       p = 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);
+       l21 = lwline_construct(-1, NULL, pa21);
+       l22 = lwline_construct(-1, NULL, pa22);
+       pa51 = ptarray_construct(0, 0, 5);
+       pa52 = ptarray_construct(0, 0, 5);
+       l51 = lwline_construct(-1, NULL, pa51);
+       l52 = lwline_construct(-1, NULL, pa52);
+       return 0;
+}
+
+/*
+** The suite cleanup function.
+** Frees any global objects.
+*/
+int clean_cg_suite(void)
+{
+       lwfree(p);
+       lwfree(p1);
+       lwfree(p2);
+       lwfree(q1);
+       lwfree(q2);
+       lwfree(pa21);
+       lwfree(pa22);
+       lwfree(pa51);
+       lwfree(pa52);
+       lwfree(l21);
+       lwfree(l22);
+       lwfree(l51);
+       lwfree(l52);
+       return 0;
+}
+
+/*
+** Test left/right side.
+*/
+void testSegmentSide(void)
+{
+       double rv = 0.0;
+       POINT2D *q = NULL;
+       q = lwalloc(sizeof(POINT2D));
+       
+       /* Vertical line at x=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;
+       rv = segmentSide(p1, p2, q);
+       //printf("left %g\n",rv);
+       CU_ASSERT(rv < 0.0);
+       
+       /* On the right */
+       q->x = 2.0;
+       rv = segmentSide(p1, p2, q);
+       //printf("right %g\n",rv);
+       CU_ASSERT(rv > 0.0);
+       
+       /* On the line */
+       q->x = 0.0;
+       rv = segmentSide(p1, p2, q);
+       //printf("on line %g\n",rv);
+       CU_ASSERT_EQUAL(rv, 0.0);
+       
+       lwfree(q);
+       
+}
+
+/*
+** Test crossings side.
+*/
+void testSegmentIntersects(void)
+{
+       
+       /* P: Vertical line at x=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;
+       CU_ASSERT( segmentIntersects(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;
+       CU_ASSERT( segmentIntersects(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;
+       CU_ASSERT( segmentIntersects(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;
+       CU_ASSERT( segmentIntersects(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;
+       CU_ASSERT( segmentIntersects(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;
+       CU_ASSERT( segmentIntersects(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( segmentIntersects(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( segmentIntersects(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( segmentIntersects(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( segmentIntersects(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( segmentIntersects(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( segmentIntersects(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( segmentIntersects(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( segmentIntersects(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( segmentIntersects(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( segmentIntersects(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( segmentIntersects(p1, p2, q1, q2) == SEG_TOUCH_RIGHT );
+
+}
+
+void testLineCrossingShortLines(void) 
+{
+
+       int rv = 0;
+       
+       /*
+       ** Simple test, two two-point lines 
+       */
+
+       /* 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);
+       
+       /* Horizontal, crossing mid-segment */
+       p->x = -0.5;
+       p->y = 0.5;
+       setPoint4d(pa22, 0, p);
+       p->x = 0.5;
+       setPoint4d(pa22, 1, p);
+
+       CU_ASSERT( lineCrossingDirection(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);
+
+       CU_ASSERT( lineCrossingDirection(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( lineCrossingDirection(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);
+
+       CU_ASSERT( lineCrossingDirection(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);
+
+       CU_ASSERT( lineCrossingDirection(l21, l22) == LINE_NO_CROSS );
+
+}
+       
+void testLineCrossingLongLines(void) 
+{
+
+       /* 
+       ** More complex test, longer lines and multiple crossings 
+       */
+
+       /* Vertical line with vertices at y integers */
+       p->x = 0.0;
+       p->y = 0.0;
+       setPoint4d(pa51, 0, p);
+       p->y = 1.0;
+       setPoint4d(pa51, 1, p);
+       p->y = 2.0;
+       setPoint4d(pa51, 2, p);
+       p->y = 3.0;
+       setPoint4d(pa51, 3, p);
+       p->y = 4.0;
+       setPoint4d(pa51, 4, p);
+
+       /* Two crossings at segment midpoints */
+       p->x = 1.0;
+       p->y = 1.0;
+       setPoint4d(pa52, 0, p);
+       p->x = -1.0;
+       p->y = 1.5;
+       setPoint4d(pa52, 1, p);
+       p->x = 1.0;
+       p->y = 3.0;
+       setPoint4d(pa52, 2, p);
+       p->x = 1.0;
+       p->y = 4.0;
+       setPoint4d(pa52, 3, p);
+       p->x = 1.0;
+       p->y = 5.0;
+       setPoint4d(pa52, 4, p);
+       CU_ASSERT( lineCrossingDirection(l51, l52) == LINE_MULTICROSS_END_SAME_FIRST_LEFT );
+
+       /* One crossing at interior vertex */
+       p->x = 1.0;
+       p->y = 1.0;
+       setPoint4d(pa52, 0, p);
+       p->x = 0.0;
+       p->y = 1.0;
+       setPoint4d(pa52, 1, p);
+       p->x = -1.0;
+       p->y = 1.0;
+       setPoint4d(pa52, 2, p);
+       p->x = -1.0;
+       p->y = 2.0;
+       setPoint4d(pa52, 3, p);
+       p->x = -1.0;
+       p->y = 3.0;
+       setPoint4d(pa52, 4, p);
+       CU_ASSERT( lineCrossingDirection(l51, l52) == LINE_CROSS_LEFT );
+
+       /* Two crossings at interior vertices */
+       p->x = 1.0;
+       p->y = 1.0;
+       setPoint4d(pa52, 0, p);
+       p->x = 0.0;
+       p->y = 1.0;
+       setPoint4d(pa52, 1, p);
+       p->x = -1.0;
+       p->y = 1.0;
+       setPoint4d(pa52, 2, p);
+       p->x = 0.0;
+       p->y = 3.0;
+       setPoint4d(pa52, 3, p);
+       p->x = 1.0;
+       p->y = 3.0;
+       setPoint4d(pa52, 4, p);
+       CU_ASSERT( lineCrossingDirection(l51, l52) == LINE_MULTICROSS_END_SAME_FIRST_LEFT );
+
+       /* Two crossings, one at the first vertex on at interior vertex */
+       p->x = 1.0;
+       p->y = 0.0;
+       setPoint4d(pa52, 0, p);
+       p->x = 0.0;
+       p->y = 0.0;
+       setPoint4d(pa52, 1, p);
+       p->x = -1.0;
+       p->y = 1.0;
+       setPoint4d(pa52, 2, p);
+       p->x = 0.0;
+       p->y = 3.0;
+       setPoint4d(pa52, 3, p);
+       p->x = 1.0;
+       p->y = 3.0;
+       setPoint4d(pa52, 4, p);
+       CU_ASSERT( lineCrossingDirection(l51, l52) == LINE_MULTICROSS_END_SAME_FIRST_LEFT );
+
+       /* Two crossings, one at the first vertex on the next interior vertex */
+       p->x = 1.0;
+       p->y = 0.0;
+       setPoint4d(pa52, 0, p);
+       p->x = 0.0;
+       p->y = 0.0;
+       setPoint4d(pa52, 1, p);
+       p->x = -1.0;
+       p->y = 1.0;
+       setPoint4d(pa52, 2, p);
+       p->x = 0.0;
+       p->y = 1.0;
+       setPoint4d(pa52, 3, p);
+       p->x = 1.0;
+       p->y = 2.0;
+       setPoint4d(pa52, 4, p);
+       CU_ASSERT( lineCrossingDirection(l51, l52) == LINE_MULTICROSS_END_SAME_FIRST_LEFT );
+
+       /* Two crossings, one at the last vertex on the next interior vertex */
+       p->x = 1.0;
+       p->y = 4.0;
+       setPoint4d(pa52, 0, p);
+       p->x = 0.0;
+       p->y = 4.0;
+       setPoint4d(pa52, 1, p);
+       p->x = -1.0;
+       p->y = 3.0;
+       setPoint4d(pa52, 2, p);
+       p->x = 0.0;
+       p->y = 3.0;
+       setPoint4d(pa52, 3, p);
+       p->x = 1.0;
+       p->y = 3.0;
+       setPoint4d(pa52, 4, p);
+       CU_ASSERT( lineCrossingDirection(l51, l52) == LINE_MULTICROSS_END_SAME_FIRST_LEFT );
+
+       /* Three crossings, two at midpoints, one at vertex */
+       p->x = 0.5;
+       p->y = 1.0;
+       setPoint4d(pa52, 0, p);
+       p->x = -1.0;
+       p->y = 0.5;
+       setPoint4d(pa52, 1, p);
+       p->x = 1.0;
+       p->y = 2.0;
+       setPoint4d(pa52, 2, p);
+       p->x = -1.0;
+       p->y = 2.0;
+       setPoint4d(pa52, 3, p);
+       p->x = -1.0;
+       p->y = 3.0;
+       setPoint4d(pa52, 4, p);
+       CU_ASSERT( lineCrossingDirection(l51, l52) == LINE_MULTICROSS_END_LEFT );
+
+}
+
diff --git a/cunit/cu_algorithm.h b/cunit/cu_algorithm.h
new file mode 100644 (file)
index 0000000..fbac994
--- /dev/null
@@ -0,0 +1,23 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include "CUnit/Basic.h"
+
+#include "lwalgorithm.h"
+
+/***********************************************************************
+** for Computational Geometry Suite
+*/
+
+/* Set-up / clean-up functions */
+CU_pSuite register_cg_suite(void);
+int init_cg_suite(void);
+int clean_cg_suite(void);
+
+/* Test functions */
+void testSegmentSide(void);
+void testSegmentIntersects(void);
+void testLineCrossingShortLines(void);
+void testLineCrossingLongLines(void);
+
+
diff --git a/cunit/cu_tester.c b/cunit/cu_tester.c
new file mode 100644 (file)
index 0000000..21e7854
--- /dev/null
@@ -0,0 +1,56 @@
+
+#include <stdio.h>
+#include <string.h>
+#include "CUnit/Basic.h"
+
+#include "cu_algorithm.h"
+
+
+/*
+** Set up liblwgeom to run in stand-alone mode using the 
+** usual system memory handling functions.
+*/
+void cunit_lwerr(const char *fmt, va_list ap) {
+       printf("Got an LWERR.\n");
+}
+
+void cunit_lwnotice(const char *fmt, va_list ap) {
+       printf("Got an LWNOTICE.\n");
+}
+
+void lwgeom_init_allocators(void) {
+        /* liblwgeom callback - install PostgreSQL handlers */
+        lwalloc_var = malloc;
+        lwrealloc_var = realloc;
+        lwfree_var = free;
+        lwerror_var = default_errorreporter;
+        lwnotice_var = default_noticereporter;
+}
+
+/*
+** The main() function for setting up and running the tests.
+** Returns a CUE_SUCCESS on successful running, another
+** CUnit error code on failure.
+*/
+int main()
+{
+
+       /* initialize the CUnit test registry */
+       if (CUE_SUCCESS != CU_initialize_registry())
+               return CU_get_error();
+
+       /* Add the PostGIS suite to the registry */
+       if (NULL == register_cg_suite())
+       {
+               CU_cleanup_registry();
+               return CU_get_error();
+       }
+
+       /* Run all tests using the CUnit Basic interface */
+       CU_basic_set_mode(CU_BRM_VERBOSE);
+       CU_basic_run_tests();
+       CU_cleanup_registry();
+
+       return CU_get_error();
+
+}
index 854f4a30ca5e547790150e1ae6128594d4044de6..288d4a9d39bc4bda1e853bf7d9caadcc4c49ec5b 100644 (file)
@@ -35,6 +35,7 @@ SA_OBJS=measures.o \
        lwmcurve.o \
        lwmsurface.o \
        lwutil.o \
+       lwalgorithm.o \
        lwgunparse.o \
        lwgparse.o \
        lwsegmentize.o \
index f7b84e8fba794a4fb7f9b791f0637f32d0e2ae87..b6b4eaaf70953c7c0fced9c2559e8aced4ada25e 100644 (file)
@@ -22,6 +22,9 @@
 #define FP_CONTAINS_INCL(A, X, B) (FP_LTEQ(A, X) && FP_LTEQ(X, B))
 #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 LW_TRUE 1
+#define LW_FALSE 0
+
 
 /*
  * Memory management function types
diff --git a/liblwgeom/lwalgorithm.c b/liblwgeom/lwalgorithm.c
new file mode 100644 (file)
index 0000000..d23eaaf
--- /dev/null
@@ -0,0 +1,262 @@
+#include "lwalgorithm.h"
+
+/*
+** segmentSide()
+**
+** Return < 0.0 if point Q is left of segment P
+** Return > 0.0 if point Q is right of segment P
+** Return = 0.0 if point Q in on segment P
+*/
+double segmentSide(POINT2D *p1, POINT2D *p2, POINT2D *q)
+{
+       return ( (q->x - p1->x) * (p2->y - p1->y) - (p2->x - p1->x) * (q->y - p1->y) );
+}
+
+int segmentEnvelopeIntersects(POINT2D *p1, POINT2D *p2, POINT2D *q1, POINT2D *q2) {
+       double minq=LW_MIN(q1->x,q2->x);
+       double maxq=LW_MAX(q1->x,q2->x);
+       double minp=LW_MIN(p1->x,p2->x);
+       double maxp=LW_MAX(p1->x,p2->x);
+
+       if(minp>maxq || maxp<minq) 
+               return LW_FALSE;
+
+       minq=LW_MIN(q1->y,q2->y);
+       maxq=LW_MAX(q1->y,q2->y);
+       minp=LW_MIN(p1->y,p2->y);
+       maxp=LW_MAX(p1->y,p2->y);
+
+       if(minp>maxq || maxp<minq) 
+               return LW_FALSE;
+
+       return LW_TRUE;
+}
+
+/*
+** segmentIntersects()
+** 
+** Returns one of
+**     SEG_ERROR = -1,
+**     SEG_NO_INTERSECTION = 0, 
+**     SEG_COLINEAR = 1, 
+**     SEG_CROSS_LEFT = 2, 
+**     SEG_CROSS_RIGHT = 3, 
+**     SEG_TOUCH_LEFT = 4, 
+**     SEG_TOUCH_RIGHT = 5
+*/
+int segmentIntersects(POINT2D *p1, POINT2D *p2, POINT2D *q1, POINT2D *q2) {
+       
+       double pq1, pq2, qp1, qp2;
+       
+       /* No envelope interaction => we are done. */
+       if (!segmentEnvelopeIntersects(p1, p2, q1, p2)) {
+               return SEG_NO_INTERSECTION;
+       }
+
+       /* Are the start and end points of q on the same side of p? */
+       pq1=segmentSide(p1,p2,q1);
+       pq2=segmentSide(p1,p2,q2);
+       if ((pq1>0 && pq2>0) || (pq1<0 && pq2<0)) {
+               return SEG_NO_INTERSECTION;
+       }
+
+       /* Are the start and end points of p on the same side of q? */
+       qp1=segmentSide(q1,q2,p1);
+       qp2=segmentSide(q1,q2,p2);
+       if ((qp1>0 && qp2>0) || (qp1<0 && qp2<0)) {
+               return SEG_NO_INTERSECTION;
+       }
+
+       /* Nobody is on one side or another? Must be colinear. */
+       if (pq1 == 0.0 && pq2 == 0.0 && qp1 == 0.0 && qp2 == 0.0) {
+               return SEG_COLINEAR;
+       }
+
+       /* 
+       ** When one end-point touches, the sidedness is determined by the 
+       ** location of the other end-point.
+       */
+       if ( pq2 == 0.0 ) {
+               if ( pq1 < 0.0 ) 
+                       return SEG_TOUCH_LEFT;
+               else 
+                       return SEG_TOUCH_RIGHT;
+       }
+       if ( pq1 == 0.0 ) {
+               if ( pq2 < 0.0 ) 
+                       return SEG_TOUCH_LEFT;
+               else 
+                       return SEG_TOUCH_RIGHT;
+       }
+
+       /* The segments cross, what direction is the crossing? */
+       if ( pq1 < pq2 ) 
+               return SEG_CROSS_RIGHT;
+       else 
+               return SEG_CROSS_LEFT;
+       
+       /* This should never happen! */
+       return SEG_ERROR;
+               
+}
+
+/*
+** lineCrossingDirection()
+** 
+** Returns one of
+**   LINE_NO_CROSS = 0
+**   LINE_CROSS_LEFT = -1
+**   LINE_CROSS_RIGHT = 1
+**   LINE_MULTICROSS_END_LEFT = -2
+**   LINE_MULTICROSS_END_RIGHT = 2
+**   LINE_MULTICROSS_END_SAME_FIRST_LEFT = -3
+**   LINE_MULTICROSS_END_SAME_FIRST_RIGHT = 3
+**   LINE_TOUCH_LEFT = -4
+**   LINE_TOUCH_RIGHT = 4
+**
+*/
+int lineCrossingDirection(LWLINE *l1, LWLINE *l2) {
+       
+       int i = 0, j = 0, rv = 0;
+       POINT2D *p1;
+       POINT2D *p2;
+       POINT2D *q1;
+       POINT2D *q2;
+       POINTARRAY *pa1;
+       POINTARRAY *pa2;
+       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). */
+       if( pa1->npoints < 2 || pa2->npoints < 2 ) 
+               return LINE_NO_CROSS;
+
+       LWDEBUGF(5, "lineCrossingDirection: l1 = %s", lwgeom_to_ewkt((LWGEOM*)l1,0));
+       LWDEBUGF(5, "lineCrossingDirection: l2 = %s", lwgeom_to_ewkt((LWGEOM*)l2,0));
+
+       for ( i = 1; i < pa2->npoints; i++ ) {
+               
+               rv = getPoint2d_p(pa2, i-1, q1);
+               rv = getPoint2d_p(pa2, i, q2);
+               
+               for ( j = 1; j < pa1->npoints; j++ ) {   
+                       
+                       rv = getPoint2d_p(pa1, j-1, p1);
+                       rv = getPoint2d_p(pa1, j, p2);
+                       
+                       LWDEBUGF(5, "lineCrossingDirection: i=%d, j=%d", i, j);
+                       
+                       this_cross = segmentIntersects(p1, p2, q1, q2);
+               
+                       if( ! first_cross && this_cross )
+                               first_cross = this_cross;
+                       if( this_cross ) 
+                               final_cross = this_cross;
+                               
+                       if( this_cross == SEG_CROSS_LEFT ) {
+                               cross_left++;
+                               break;
+                       }
+                               
+                       if( this_cross == SEG_CROSS_RIGHT ) {
+                               cross_right++;
+                               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(5, "lineCrossingDirection: this_cross=%d, prev_cross=%d, vertex_touch=%d, vertex_touch_type=%d", this_cross, prev_cross, vertex_touch, vertex_touch_type);
+                               
+               }
+               
+       }
+
+       LWDEBUGF(5, "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( !cross_left && !cross_right && !first_cross && !final_cross ) 
+               return LINE_NO_CROSS;
+
+       if( !cross_left && !cross_right && ( first_cross == SEG_TOUCH_RIGHT || final_cross == SEG_TOUCH_RIGHT ) )
+               return LINE_TOUCH_RIGHT;
+
+       if( !cross_left && !cross_right && ( first_cross == SEG_TOUCH_LEFT || final_cross == SEG_TOUCH_LEFT ) )
+               return LINE_TOUCH_LEFT;
+               
+       if( !cross_left && cross_right == 1 )
+               return LINE_CROSS_RIGHT;
+               
+       if( !cross_right && cross_left == 1 )
+               return LINE_CROSS_LEFT;
+
+       if( cross_left - cross_right == 1 )
+               return LINE_MULTICROSS_END_LEFT;
+
+       if( cross_left - cross_right == -1 )
+               return LINE_MULTICROSS_END_RIGHT;
+
+       if( cross_left - cross_right == 0 && first_cross == SEG_CROSS_LEFT )
+               return LINE_MULTICROSS_END_SAME_FIRST_LEFT;
+
+       if( cross_left - cross_right == 0 && first_cross == SEG_CROSS_RIGHT )
+               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;
+       
+}
+
diff --git a/liblwgeom/lwalgorithm.h b/liblwgeom/lwalgorithm.h
new file mode 100644 (file)
index 0000000..7e2329c
--- /dev/null
@@ -0,0 +1,30 @@
+#include "liblwgeom.h"
+
+enum CG_SEGMENT_INTERSECTION_TYPE { 
+       SEG_ERROR = -1,
+       SEG_NO_INTERSECTION = 0, 
+       SEG_COLINEAR = 1, 
+       SEG_CROSS_LEFT = 2, 
+       SEG_CROSS_RIGHT = 3, 
+       SEG_TOUCH_LEFT = 4, 
+       SEG_TOUCH_RIGHT = 5
+};
+
+double segmentSide(POINT2D *p1, POINT2D *p2, POINT2D *q);
+int segmentIntersects(POINT2D *p1, POINT2D *p2, POINT2D *q1, POINT2D *q2);
+int segmentEnvelopeIntersects(POINT2D *p1, POINT2D *p2, POINT2D *q1, POINT2D *q2);
+
+
+enum CG_LINE_CROSS_TYPE {
+       LINE_NO_CROSS = 0,
+       LINE_CROSS_LEFT = -1,
+       LINE_CROSS_RIGHT = 1,
+       LINE_MULTICROSS_END_LEFT = -2,
+       LINE_MULTICROSS_END_RIGHT = 2,
+       LINE_MULTICROSS_END_SAME_FIRST_LEFT = -3,
+       LINE_MULTICROSS_END_SAME_FIRST_RIGHT = 3,
+       LINE_TOUCH_LEFT = -4,
+       LINE_TOUCH_RIGHT = 4
+};
+
+int lineCrossingDirection(LWLINE *l1, LWLINE *l2);