From: Paul Ramsey Date: Thu, 1 Jan 2009 00:31:47 +0000 (+0000) Subject: Fix based on test case from MAC. X-Git-Tag: 1.4.0b1~351 X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=761363721f5f39e9330f8852cf072cc5bf4951d7;p=postgis 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 ef70971ca..e9648ecb4 100644 --- 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 28c37595c..78921683b 100644 --- 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 e9bb90470..b1c65c028 100644 --- 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); }