Fix based on test case from MAC.

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

diff --git a/liblwgeom/cunit/cu_algorithm.c b/liblwgeom/cunit/cu_algorithm.c
index ef70971ca2a5ed8d5e6672bfd8cc7196eefaa517..e9648ecb47d3e6d33b432c7eefd2067f7cc10697 100644 (file)
--- a/liblwgeom/cunit/cu_algorithm.c
+++ b/liblwgeom/cunit/cu_algorithm.c
@@ -93,14 +93,14 @@ int init_cg_suite(void)
 */
 int clean_cg_suite(void)
 {
-       lwfree(p);
-       lwfree(p1);
-       lwfree(p2);
-       lwfree(q1);
-       lwfree(q2);
-       lwfree_line(l21);
-       lwfree_line(l22);
-       lwfree_line(l51);
+       if ( p ) lwfree(p);
+       if ( p1 )lwfree(p1);
+       if ( p2 ) lwfree(p2);
+       if ( q1 ) lwfree(q1);
+       if ( q2 ) lwfree(q2);
+       if ( l21 ) lwfree_line(l21);
+       if ( l22 ) lwfree_line(l22);
+       if ( l51 ) lwfree_line(l51);
        return 0;
 }
 
@@ -533,6 +533,7 @@ void test_lwmline_clip(void)
        LWCOLLECTION *c;
        char *ewkt;
        LWMLINE *mline = NULL;
+       LWLINE *line = NULL;
 
        /*
        ** Set up the input line. Trivial one-member case. 
@@ -579,6 +580,46 @@ void test_lwmline_clip(void)
 
        lwfree_mline(mline);
 
+       /* 
+       ** Set up input line from MAC
+       */
+       line = (LWLINE*)lwgeom_from_ewkt("LINESTRING(0 0 0 0,1 1 1 1,2 2 2 2,3 3 3 3,4 4 4 4,3 3 3 5,2 2 2 6,1 1 1 7,0 0 0 8)",0);
+
+       /* Clip from 3 to 3.5 */
+       c = lwline_clip_to_ordinate_range(line, 2, 3.0, 3.5);
+       ewkt = lwgeom_to_ewkt((LWGEOM*)c,0);
+       //printf("c = %s\n", ewkt);
+       CU_ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((3 3 3 3,3.5 3.5 3.5 3.5),(3.5 3.5 3.5 4.5,3 3 3 5))");
+       lwfree(ewkt);
+       lwfree_collection(c);
+
+       /* Clip from 2 to 3.5 */
+       c = lwline_clip_to_ordinate_range(line, 2, 2.0, 3.5);
+       ewkt = lwgeom_to_ewkt((LWGEOM*)c,0);
+       //printf("c = %s\n", ewkt);
+       CU_ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((2 2 2 2,3 3 3 3,3.5 3.5 3.5 3.5),(3.5 3.5 3.5 4.5,3 3 3 5,2 2 2 6))");
+       lwfree(ewkt);
+       lwfree_collection(c);
+
+       /* Clip from 3 to 4 */
+       c = lwline_clip_to_ordinate_range(line, 2, 3.0, 4.0);
+       ewkt = lwgeom_to_ewkt((LWGEOM*)c,0);
+       //printf("c = %s\n", ewkt);
+       CU_ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((3 3 3 3,4 4 4 4,3 3 3 5))");
+       lwfree(ewkt);
+       lwfree_collection(c);
+
+       /* Clip from 2 to 3 */
+       c = lwline_clip_to_ordinate_range(line, 2, 2.0, 3.0);
+       ewkt = lwgeom_to_ewkt((LWGEOM*)c,0);
+       //printf("c = %s\n", ewkt);
+       CU_ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((2 2 2 2,3 3 3 3),(3 3 3 5,2 2 2 6))");
+       lwfree(ewkt);
+       lwfree_collection(c);
+
+
+       lwfree_line(line);
+
 }
 
 

diff --git a/liblwgeom/lwalgorithm.c b/liblwgeom/lwalgorithm.c
index 28c37595cebfc76338041bb20b279ef7b17bc284..78921683ba885e3a78c22a299775ca8b323a1076 100644 (file)
--- a/liblwgeom/lwalgorithm.c
+++ b/liblwgeom/lwalgorithm.c
@@ -369,13 +369,14 @@ int lwpoint_interpolate(const POINT4D *p1, const POINT4D *p2, POINT4D *p, int nd
         return 0;
     } 
     
-    proportion = (interpolation_value - p1_value) / fabs(p2_value - p1_value);
-        
+    proportion = fabs((interpolation_value - p1_value) / (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));
+        lwpoint_set_ordinate(p, i, p1_value + proportion * (p2_value - p1_value));
+               LWDEBUGF(1, "   clip ordinate(%d) p1_value(%g) p2_value(%g) proportion(%g)", i, p1_value, p2_value, proportion );
     }
     
     return 1;
@@ -530,26 +531,13 @@ LWCOLLECTION *lwline_clip_to_ordinate_range(LWLINE *line, int ordinate, double f
                }
                rv = getPoint4d_p(pa_in, i, p);
                ordinate_value_p = lwpoint_get_ordinate(p, ordinate);
-               LWDEBUGF(1, "ordinate_value_p %g", ordinate_value_p);
-               LWDEBUGF(1, "ordinate_value_q %g", ordinate_value_q);
-
-               /* Is this point on the edges of our range? */
-               if( ordinate_value_p == from || ordinate_value_p == to ) 
-               {
-                       if( ! added_last_point ) 
-                       {       
-                               if( dp ) lwfree(dp);
-                               dp = dynptarray_create(64, line->type);
-                       }
-                       rv = dynptarray_addPoint4d(dp, p, 0);
-            added_last_point = 2; /* Added on a boundary. */
-                       continue;
-               }
+               LWDEBUGF(1, " ordinate_value_p %g (current)", ordinate_value_p);
+               LWDEBUGF(1, " ordinate_value_q %g (previous)", ordinate_value_q);
 
                /* Is this point inside the ordinate range? Yes. */
-               if ( ordinate_value_p > from && ordinate_value_p < to )
+               if ( ordinate_value_p >= from && ordinate_value_p <= to )
                {
-                       LWDEBUGF(1, "inside ordinate range (%g, %g)", from, to);
+                       LWDEBUGF(1, " inside ordinate range (%g, %g)", from, to);
 
                        if ( ! added_last_point ) 
                        {
@@ -559,22 +547,36 @@ LWCOLLECTION *lwline_clip_to_ordinate_range(LWLINE *line, int ordinate, double f
                                dp = dynptarray_create(64, line->type);
 
                        /* We're transiting into the range so add an interpolated 
-                               *  point at the range boundary. */
-                       if ( ordinate_value_p > from && ordinate_value_p < to && i > 0 )
+                               *  point at the range boundary. 
+                               *  If we're on a boundary and crossing from the far side, 
+                               *  we also need an interpolated point. */
+                       if ( i > 0 && ( /* Don't try to interpolate if this is the first point */
+                                    ( ordinate_value_p > from && ordinate_value_p < to ) || /* Inside */
+                                    ( ordinate_value_p == from && ordinate_value_q > to ) || /* Hopping from above */
+                                    ( ordinate_value_p == to && ordinate_value_q < from ) ) ) /* Hopping from below */
                    {
                     double interpolation_value;
                     (ordinate_value_q > to) ? (interpolation_value = to) : (interpolation_value = from);
                     rv = lwpoint_interpolate(q, p, r, dims, ordinate, interpolation_value);
                     rv = dynptarray_addPoint4d(dp, r, 0);
+                                       LWDEBUGF(1, " interpolating between (%g, %g) with interpolation point (%g)", ordinate_value_q, ordinate_value_p, interpolation_value);
                    }
                        }
                    /* Add the current vertex to the point array. */
                rv = dynptarray_addPoint4d(dp, p, 0);
-            added_last_point = 1; /* Added inside range. */
+                       if ( ordinate_value_p == from || ordinate_value_p == to ) 
+                       {
+               added_last_point = 2; /* Added on boundary. */
+                       }
+                       else
+                       {
+               added_last_point = 1; /* Added inside range. */
+                       }
                } 
                /* Is this point inside the ordinate range? No. */
                else 
                {
+                       LWDEBUGF(4, "  added_last_point (%d)", added_last_point);
                    if( added_last_point == 1 ) 
                    {
                        /* We're transiting out of the range, so add an interpolated point 
@@ -583,11 +585,22 @@ LWCOLLECTION *lwline_clip_to_ordinate_range(LWLINE *line, int ordinate, double f
                 (ordinate_value_p > to) ? (interpolation_value = to) : (interpolation_value = from);
                 rv = lwpoint_interpolate(q, p, r, dims, ordinate, interpolation_value);
                 rv = dynptarray_addPoint4d(dp, r, 0);
+                               LWDEBUGF(1, " interpolating between (%g, %g) with interpolation point (%g)", ordinate_value_q, ordinate_value_p, interpolation_value);
                }
                        else if ( added_last_point == 2 ) 
                        {
                                /* We're out and the last point was on the boundary.
-                               *  Nothing to do, the point was already added above. */
+                               *  If the last point was the near boundary, nothing to do.
+                               *  If it was the far boundary, we need an interpolated point. */
+                               if( (ordinate_value_q == from && ordinate_value_p > from) || 
+                                   (ordinate_value_q == to && ordinate_value_p < to) ) 
+                               {
+                       double interpolation_value;
+                       (ordinate_value_p > to) ? (interpolation_value = to) : (interpolation_value = from);
+                       rv = lwpoint_interpolate(q, p, r, dims, ordinate, interpolation_value);
+                       rv = dynptarray_addPoint4d(dp, r, 0);
+                                       LWDEBUGF(1, " interpolating between (%g, %g) with interpolation point (%g)", ordinate_value_q, ordinate_value_p, interpolation_value);
+                               }
                        }
                        else if ( ordinate_value_q < from && ordinate_value_p > to ) {
                                /* We just hopped over the whole range, from bottom to top, 

diff --git a/liblwgeom/lwgeom_api.c b/liblwgeom/lwgeom_api.c
index e9bb90470dae7ea158aff126bbf57be01e7a524a..b1c65c0280034cc010494c7b22ff4938e85528ec 100644 (file)
--- a/liblwgeom/lwgeom_api.c
+++ b/liblwgeom/lwgeom_api.c
@@ -795,8 +795,7 @@ pointArray_construct(uchar *points, char hasz, char hasm,
 int
 pointArray_ptsize(const POINTARRAY *pa)
 {
-       LWDEBUGF(2, "pointArray_ptsize: TYPE_NDIMS(pa->dims)=%x\n",
-               TYPE_NDIMS(pa->dims));
+       LWDEBUGF(2, "pointArray_ptsize: TYPE_NDIMS(pa->dims)=%x",TYPE_NDIMS(pa->dims));
 
        return sizeof(double)*TYPE_NDIMS(pa->dims);
 }