Partial work saved back for later.

git-svn-id: http://svn.osgeo.org/postgis/trunk@3440 b70326c6-7e19-0410-871a-916f4a2858ee

diff --git a/liblwgeom/Makefile.in b/liblwgeom/Makefile.in
index 288d4a9d39bc4bda1e853bf7d9caadcc4c49ec5b..7d81d3045e96e773d15f557867d21b734e459b15 100644 (file)
--- a/liblwgeom/Makefile.in
+++ b/liblwgeom/Makefile.in
@@ -17,7 +17,8 @@ YACC=@YACC@
 LEX=@LEX@
 
 # Standalone LWGEOM objects
-SA_OBJS=measures.o \
+SA_OBJS = \
+       measures.o \
        box2d.o \
        ptarray.o \
        lwgeom_api.o \
@@ -43,9 +44,13 @@ SA_OBJS=measures.o \
        lex.yy.o \
        vsprintf.o      
 
+SA_HEADERS = \
+       liblwgeom.h \
+       lwalgorithm.h
+
 all: liblwgeom.a
 
-liblwgeom.a: $(SA_OBJS)
+liblwgeom.a: $(SA_OBJS) $(SA_HEADERS) 
        ar rs liblwgeom.a $(SA_OBJS)    
 
 clean:
@@ -53,7 +58,7 @@ clean:
        rm -f liblwgeom.a 
 
 # Command to build each of the .o files
-$(SA_OBJS): %.o: %.c
+$(SA_OBJS): %.o: %.c 
        $(CC) $(CFLAGS) -c -o $@ $<
 
 # Commands to generate the lexer and parser from input files

diff --git a/liblwgeom/cunit/cu_algorithm.c b/liblwgeom/cunit/cu_algorithm.c
index e19c1101b06d759610af2ce0a7bc7e452b94b6ef..ac91437b472bdea61615edeb24cd298ffb4b2390 100644 (file)
--- a/liblwgeom/cunit/cu_algorithm.c
+++ b/liblwgeom/cunit/cu_algorithm.c
@@ -29,7 +29,10 @@ CU_pSuite register_cg_suite(void)
            (NULL == CU_add_test(pSuite, "test_lw_segment_side()", test_lw_segment_side)) ||
            (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_lwpoint_set_ordinate()", test_lwpoint_set_ordinate)) || 
+           (NULL == CU_add_test(pSuite, "test_lwpoint_get_ordinate()", test_lwpoint_get_ordinate)) ||
+           (NULL == CU_add_test(pSuite, "test_lwpoint_interpolate()", test_lwpoint_interpolate)) 
        )
        {
                CU_cleanup_registry();
@@ -47,6 +50,7 @@ POINT2D *p2 = NULL;
 POINT2D *q1 = NULL;
 POINT2D *q2 = NULL;
 POINT4D *p = NULL;
+POINT4D *q = NULL;
 /* Two-point objects */
 POINTARRAY *pa21 = NULL;
 POINTARRAY *pa22 = NULL;
@@ -66,6 +70,7 @@ LWLINE *l52 = NULL;
 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));
@@ -540,3 +545,60 @@ void test_lwline_crossing_long_lines(void)
 
 }
 
+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 );
+    
+}
+
+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 );
+        
+}
+
+void test_lwpoint_interpolate(void)
+{
+    POINT4D *r = NULL;
+    r = lwalloc(sizeof(POINT4D));
+    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);
+
+    lwfree(r);
+    
+}

diff --git a/liblwgeom/cunit/cu_algorithm.h b/liblwgeom/cunit/cu_algorithm.h
index a38f5d591502b3ede0c3bb6151c9d459f2fd9f8a..a2223c37f9cdeab4a7076d74361e429306b7f857 100644 (file)
--- a/liblwgeom/cunit/cu_algorithm.h
+++ b/liblwgeom/cunit/cu_algorithm.h
@@ -31,4 +31,6 @@ void test_lw_segment_side(void);
 void test_lw_segment_intersects(void);
 void test_lwline_crossing_short_lines(void);
 void test_lwline_crossing_long_lines(void);
-
+void test_lwpoint_set_ordinate(void);
+void test_lwpoint_get_ordinate(void);
+void test_lwpoint_interpolate(void);

diff --git a/liblwgeom/lwalgorithm.c b/liblwgeom/lwalgorithm.c
index 63e9b3c044aec899cdd418c7f6dd48af271b6645..33d39f960a991f2c855d517d8879c8994233e286 100644 (file)
--- a/liblwgeom/lwalgorithm.c
+++ b/liblwgeom/lwalgorithm.c
@@ -292,7 +292,6 @@ int lwline_crossing_direction(LWLINE *l1, LWLINE *l2)
        
 }
 
-#if 0 
 /*
 ** lwpoint_get_ordinate(point, ordinate) => double
 */
@@ -320,6 +319,69 @@ double lwpoint_get_ordinate(POINT4D *p, int ordinate)
        return p->x;
        
 }
+void lwpoint_set_ordinate(POINT4D *p, int ordinate, double value) 
+{
+       if( ! p ) 
+       {
+               lwerror("Null input geometry.");
+               return;
+       }
+       
+       if( ordinate > 3 || ordinate < 0 ) 
+       {
+               lwerror("Cannot extract ordinate %d.", ordinate);
+               return;
+       }
+
+       switch ( ordinate ) {
+        case 3:
+            p->m = value;
+            return;
+        case 2:
+            p->z = value;
+            return;
+        case 1:
+            p->y = value;
+            return;
+        case 0:
+            p->x = value;
+            return;
+    }          
+}
+
+
+int lwpoint_interpolate(POINT4D *p1, POINT4D *p2, POINT4D *p, int ndims, int ordinate, double interpolation_value)
+{
+    double p1_value = lwpoint_get_ordinate(p1, ordinate);
+    double p2_value = lwpoint_get_ordinate(p2, ordinate);
+    double proportion;
+    int i = 0;
+    
+    if( ordinate < 0 || ordinate >= ndims ) 
+    {
+        lwerror("Ordinate (%d) is not within ndims (%d).", ordinate, ndims);
+        return 0;
+    }
+    
+    if( FP_MIN(p1_value, p2_value) > interpolation_value || 
+        FP_MAX(p1_value, p2_value) < interpolation_value )
+    {
+        lwerror("Cannot interpolate to a value (%g) not between the input points.", interpolation_value);
+        return 0;
+    } 
+    
+    proportion = (interpolation_value - p1_value) / fabs(p2_value - p1_value);
+        
+    for( i = 0; i < ndims; i++ ) 
+    {
+        p1_value = lwpoint_get_ordinate(p1, i);
+        p2_value = lwpoint_get_ordinate(p2, i);
+        lwpoint_set_ordinate(p, i, p1_value + proportion * fabs(p2_value - p1_value));
+    }
+    
+    return 1;
+}
+
 
 /*
 ** lwline_clip_to_ordinate_range(line, ordinate, from, to) => lwmline
@@ -330,14 +392,14 @@ double lwpoint_get_ordinate(POINT4D *p, int ordinate)
 LWLINE *lwline_clip_to_ordinate_range(LWLINE *line, int ordinate, double from, double to) 
 {
        
-       POINTARRAY *pa_in;
-       LWMLINE *mline_out;
-       POINTARRAY *pa_out;
-       DYNPTARRAY *dp;
+       POINTARRAY *pa_in = NULL;
+       LWMLINE *mline_out = NULL;
+       POINTARRAY *pa_out = NULL;
+       DYNPTARRAY *dp = NULL;
        int i, rv;
-       int last_point = 0;
+       int added_last_point = 0;
        int nparts = 0;
-       POINT4D *p, *q;
+       POINT4D *p, *q, *r;
        double ordinate_value;
 
        
@@ -367,22 +429,61 @@ LWLINE *lwline_clip_to_ordinate_range(LWLINE *line, int ordinate, double from, d
        
        p = lwalloc(sizeof(POINT4D));
        q = lwalloc(sizeof(POINT4D));
+    r = lwalloc(sizeof(POINT4D));
 
        pa_in = (POINTARRAY*)line->points;
        
        dp = dynptarray_create(64, ndims);
 
-       
-       for ( i = 1; i < pa_in->npoints; i++ ) 
+       for ( i = 0; i < pa_in->npoints; i++ ) 
        {
                rv = getPoint4d_p(pa_in, i, p);
                ordinate_value = lwpoint_get_ordinate(p, ordinate);
+               /* Is this point inside the range? Yes. */
                if ( ordinate_value >= from && ordinate_value <= to )
                {
-                       rv =dynptarray_addPoint4d(dp, p, 1);
+                       if ( ! added_last_point ) 
+                       {
+                           /* TODO Make a new ptarray */
+                       if ( ordinate_value > from && ordinate_value < to &&
+                            i > 0 && i < pa_in->npoints - 1 )
+                   {
+                       /* We're transiting in so add an interpolated point */
+                    double interpolation_value;
+                    double last_value;
+                    rv = getPoint4d_p(pa_in, i-1, q);
+                    last_value = lwpoint_get_ordinate(q, ordinate);
+                    (last_value > to) ? (interpolation_value = to) : (interpolation_value = from);
+                    rv = lwpoint_interpolate(q, p, r, ndims, ordinate, interpolation_value);
+                    rv = dynptarray_addPoint4d(dp, r, 1);
+                    
+                   }
+                       }
+                   /* add the point */
+               rv = dynptarray_addPoint4d(dp, p, 1);
+            added_last_point = LW_TRUE;
                } 
+               /* Is this point inside the range? No. */
+               else 
+               {
+                   if( added_last_point ) 
+                   {
+                       /* We're transiting out, so add an interpolated point */
+                double interpolation_value;
+                rv = getPoint4d_p(pa_in, i-1, q);
+                (ordinate_value > to) ? (interpolation_value = to) : (interpolation_value = from);
+                rv = lwpoint_interpolate(q, p, r, ndims, ordinate, interpolation_value);
+                rv = dynptarray_addPoint4d(dp, r, 1);
+                
+                       /* TODO save back the current ptarray to a lwmline */
+               }
+            added_last_point = LW_FALSE;
+           }
        }
+       
+    lwfree(p);
+    lwfree(q);
+    lwfree(r);
 
        
 }
-#endif

diff --git a/liblwgeom/lwalgorithm.h b/liblwgeom/lwalgorithm.h
index b1ae14eb8f98606b5579fb9d5a3cb258db46354f..911c10cc6166ca3c34bd7b6a616c33ab5c5d98ea 100644 (file)
--- a/liblwgeom/lwalgorithm.h
+++ b/liblwgeom/lwalgorithm.h
@@ -10,6 +10,7 @@
  * 
  **********************************************************************/
 
+#include <math.h>
 #include "liblwgeom.h"
 
 enum CG_SEGMENT_INTERSECTION_TYPE { 
@@ -40,3 +41,6 @@ enum CG_LINE_CROSS_TYPE {
 int lwline_crossing_direction(LWLINE *l1, LWLINE *l2);
 
 double lwpoint_get_ordinate(POINT4D *p, int ordinate);
+void lwpoint_set_ordinate(POINT4D *p, int ordinate, double value);
+int lwpoint_interpolate(POINT4D *p1, POINT4D *p2, POINT4D *p, int ndims, int ordinate, double interpolation_value);
+LWLINE *lwline_clip_to_ordinate_range(LWLINE *line, int ordinate, double from, double to);