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;
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;
#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;
}