Lightbulb moment: the fix for GBT#21: locate_along_measure: wrong values, invalid...

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

diff --git a/lwgeom/lwgeom_functions_lrs.c b/lwgeom/lwgeom_functions_lrs.c
index cc31c435ce035a5faf67874f1db8235dc230400d..8a65883df9defd7d728155c7db9aac1a754f6d87 100644 (file)
--- 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 8e9b847e33ea7d99b8f8f2bea59b52de6ccb008b..5c60650a4c287834b95e5cf41ff72cfb8f9edde1 100644 (file)
--- 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;