From a99ea5d04b8a8c7a842300d839496c1891e03392 Mon Sep 17 00:00:00 2001 From: Mark Cave-Ayland Date: Wed, 3 Dec 2008 16:27:44 +0000 Subject: [PATCH] Lightbulb moment: the fix for GBT#21: locate_along_measure: wrong values, invalid data required extra work as floating point errors could still be introduced by the removal of the memcpy(). In fact it was the clipping logic that was wrong, so this patch re-adds the memcpy() in the correct place(s) and corrects the clipping flags to remove this floating point error. With thanks to Stephen Davies. git-svn-id: http://svn.osgeo.org/postgis/branches/1.3@3358 b70326c6-7e19-0410-871a-916f4a2858ee --- lwgeom/lwgeom_functions_lrs.c | 116 ++++++++++++++++++---------------- lwgeom/ptarray.c | 2 +- 2 files changed, 64 insertions(+), 54 deletions(-) diff --git a/lwgeom/lwgeom_functions_lrs.c b/lwgeom/lwgeom_functions_lrs.c index cc31c435c..8a65883df 100644 --- a/lwgeom/lwgeom_functions_lrs.c +++ b/lwgeom/lwgeom_functions_lrs.c @@ -133,12 +133,29 @@ clip_seg_by_m_range(POINT4D *p1, POINT4D *p2, double m0, double m1) if ( p1->m < m0 ) { /* - * Interpolate coordinates + * To prevent rounding errors, then if m0==m1 and p2 lies within the range, copy + * p1 as a direct copy of p2 */ - p1->x += (dX*dM0); - p1->y += (dY*dM0); - p1->z += (dZ*dM0); - p1->m = m0; + if (m0 == m1 && p2->m <= m1) + { + memcpy(p1, p2, sizeof(POINT4D)); + +#if DEBUG_INTERPOLATION + lwnotice("Projected p1 on range (as copy of p2)"); +#endif + } + else + { + /* Otherwise interpolate coordinates */ + p1->x += (dX*dM0); + p1->y += (dY*dM0); + p1->z += (dZ*dM0); + p1->m = m0; + +#if DEBUG_INTERPOLATION + lwnotice("Projected p1 on range"); +#endif + } if ( swapped ) ret |= 0x0100; else ret |= 0x0010; @@ -151,12 +168,29 @@ clip_seg_by_m_range(POINT4D *p1, POINT4D *p2, double m0, double m1) if ( p2->m > m1 ) { /* - * Interpolate coordinates + * To prevent rounding errors, then if m0==m1 and p1 lies within the range, copy + * p2 as a direct copy of p1 */ - p2->x += (dX*dM1); - p2->y += (dY*dM1); - p2->z += (dZ*dM1); - p2->m = m1; + if (m0 == m1 && p1->m >= m0) + { + memcpy(p2, p1, sizeof(POINT4D)); + +#if DEBUG_INTERPOLATION + lwnotice("Projected p2 on range (as copy of p1)"); +#endif + } + else + { + /* Otherwise interpolate coordinates */ + p2->x += (dX*dM1); + p2->y += (dY*dM1); + p2->z += (dZ*dM1); + p2->m = m1; + +#if DEBUG_INTERPOLATION + lwnotice("Projected p2 on range"); +#endif + } if ( swapped ) ret |= 0x0010; else ret |= 0x0100; @@ -207,73 +241,49 @@ ptarray_locate_between_m(POINTARRAY *ipa, double m0, double m1) #endif - clipval=clip_seg_by_m_range(&p1, &p2, m0, m1); + clipval = clip_seg_by_m_range(&p1, &p2, m0, m1); /* segment completely outside, nothing to do */ if ( ! clipval ) continue; - /* - * second point is clipped, or this is last - * point in the array: close the pointarray - * (eventually opening it if none is defined) - */ - if ( clipval & 0x0100 || i == ipa->npoints-1 ) - { #if DEBUG_LRS - lwnotice(" second point clipped"); + lwnotice(" clipped to: [ %g %g %g %g - %g %g %g %g ] clipval: %x", p1.x, p1.y, p1.z, p1.m, + p2.x, p2.y, p2.z, p2.m, clipval); #endif - /* - * No output points defined, so - * we have to open a new one and add the first point - */ - if ( dpa == NULL ) - { -#if DEBUG_LRS - lwnotice(" 1 creating new POINARRAY with first point %g,%g,%g,%g", p1.x, p1.y, p1.z, p1.m); -#endif - dpa = dynptarray_create(2, ipa->dims); - dynptarray_addPoint4d(dpa, &p1, 1); - } - -#if DEBUG_LRS - lwnotice(" 1 adding new point %g,%g,%g,%g", p2.x, p2.y, p2.z, p2.m); -#endif - dynptarray_addPoint4d(dpa, &p2, 0); -#if DEBUG_LRS - lwnotice(" closing pointarray %x with %d points", dpa->pa, dpa->pa->npoints); -#endif - ret.ptarrays[ret.nptarrays++] = dpa->pa; - lwfree(dpa); dpa=NULL; - continue; - } - /* - * If dpa==NULL we create it and add the first segment - */ - if ( dpa==NULL ) + /* If no points have been accumulated so far, then if clipval != 0 the first point must be the + start of the intersection */ + if (dpa == NULL) { #if DEBUG_LRS - lwnotice(" 3 creating new POINARRAY with first point %g,%g,%g,%g", p1.x, p1.y, p1.z, p1.m); + lwnotice(" 1 creating new POINTARRAY with first point %g,%g,%g,%g", p1.x, p1.y, p1.z, p1.m); #endif dpa = dynptarray_create(ipa->npoints-i, ipa->dims); dynptarray_addPoint4d(dpa, &p1, 1); } + /* Otherwise always add the next point, avoiding duplicates */ + if (dpa) + dynptarray_addPoint4d(dpa, &p2, 0); + /* - * In any case we add the second point (w/out allowin duplicates) + * second point has been clipped */ + if ( clipval & 0x0100 || i == ipa->npoints-1 ) + { #if DEBUG_LRS - lwnotice(" 2 adding new point %g,%g,%g,%g", p2.x, p2.y, p2.z, p2.m); + lwnotice(" closing pointarray %x with %d points", dpa->pa, dpa->pa->npoints); #endif - dynptarray_addPoint4d(dpa, &p2, 0); - + ret.ptarrays[ret.nptarrays++] = dpa->pa; + lwfree(dpa); dpa=NULL; + } } /* * if dpa!=NULL it means we didn't close it yet. * this should never happen. */ - if ( dpa != NULL ) lwerror("Something wrong with algorightm"); + if ( dpa != NULL ) lwerror("Something wrong with algorithm"); return ret; } diff --git a/lwgeom/ptarray.c b/lwgeom/ptarray.c index 8e9b847e3..5c60650a4 100644 --- a/lwgeom/ptarray.c +++ b/lwgeom/ptarray.c @@ -852,7 +852,7 @@ dynptarray_addPoint4d(DYNPTARRAY *dpa, POINT4D *p4d, int allow_duplicates) * return 0 and do nothing else if previous point in list is * equal to this one (4D equality) */ - if ( ! memcmp(p4d, &tmp, sizeof(POINT4D)) ) return 0; + if (tmp.x == p4d->x && tmp.y == p4d->y && tmp.z == p4d->z && tmp.m == p4d->m) return 0; } ++pa->npoints; -- 2.49.0