[postgis-commits] svn - r3358 - branches/1.3/lwgeom
postgis-commits at postgis.refractions.net
postgis-commits at postgis.refractions.net
Wed Dec 3 08:27:44 PST 2008
Author: mcayland
Date: 2008-12-03 08:27:44 -0800 (Wed, 03 Dec 2008)
New Revision: 3358
Modified:
branches/1.3/lwgeom/lwgeom_functions_lrs.c
branches/1.3/lwgeom/ptarray.c
Log:
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.
Modified: branches/1.3/lwgeom/lwgeom_functions_lrs.c
===================================================================
--- branches/1.3/lwgeom/lwgeom_functions_lrs.c 2008-12-03 16:02:32 UTC (rev 3357)
+++ branches/1.3/lwgeom/lwgeom_functions_lrs.c 2008-12-03 16:27:44 UTC (rev 3358)
@@ -133,13 +133,30 @@
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,13 +168,30 @@
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 @@
#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;
}
Modified: branches/1.3/lwgeom/ptarray.c
===================================================================
--- branches/1.3/lwgeom/ptarray.c 2008-12-03 16:02:32 UTC (rev 3357)
+++ branches/1.3/lwgeom/ptarray.c 2008-12-03 16:27:44 UTC (rev 3358)
@@ -852,7 +852,7 @@
* 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;
More information about the postgis-commits
mailing list