Spline Arc-Length and Even Sampling
- To: mathgroup at smc.vnet.net
- Subject: [mg37649] Spline Arc-Length and Even Sampling
- From: "James Murray" <jmurray at litecycles.com>
- Date: Thu, 7 Nov 2002 06:42:34 -0500 (EST)
- Sender: owner-wri-mathgroup at wolfram.com
All,
In searching the archives for spline related topics, I found no direct
solution for evenly sampling along the arc-length of a SplineFit function.
Several related inquiries have been posted to the group. Below is a
straight-forward, non-elegant routine that solves the problem of pathlength
sampling of a SplineFit function for a set of points in a 2D plane. The
routine works by first sampling the spline function against the spline
parameter, u, which is nonlinear with pathlength. The differential
pathlength along the spline curve, ds/du, is determined by calculating the
distance between neighboring points. The pathlength function, s(u), is then
determined by interpolating the approximated cumulative integral from the
sampled ds/du points. An inverse interpolation routine (supplied by David
Park, MathGroup [mg29221] Inverse Interpolating Functions, or similar
routines, like Carl Woll's [mg29277] NInverse function) is used to determine
u(s), which is the spline parameter as a function of the pathlength, s. This
function is supplied to the spline function to generate a list of points
that are evenly spaced along the spline curve pathlength.
I am sure there are other more elegant and accurate methods of solving this
problem, but this will do in a pinch.
James
(********************Begin Functions Definition************************)
Clear[SplineEvenSamplePoints, InverseInterpolatingFunction, pairpts,
ptdist];
(*Define the InverseInterpolatingFunction written by David Park, MathGroup \
[mg29221] Inverse Interpolating Functions*)
InverseInterpolatingFunction[(f_InterpolatingFunction)? (And @@ Positive /@
\
ListConvolve[{1, -1}, #1[[4, 2]]] & )] :=
InterpolatingFunction[{{First[f[[4, \
2]]], Last[f[[4, 2]]]}}, f[[2]], {f[[4, 2]]}, {f[[4, 1]], f[[3, 1]]}] ;
(*Define the function pairpts[pts] to pair up neighboring points in the
point \
data list, pts.*)
pairpts[pts_List] := Drop[Transpose[{pts, RotateLeft[pts]}], -1];
(*Define the
function ptdist that calculates the distance between \
two points.*)
ptdist = Sqrt[(#[[1, 1]] - #[[2, 1]])^2 + (#[[1, 2]] - #[[2, 2]])^2] & ;
Options[SplineEvenSamplePoints] = {
NumberOfSmplePoints -> 100};
(*++++++++++++++++ SplineEvenSamplePoints ++++++++ + +++++++++++++++++*)
Clear[SplineEvenSamplePoints, InverseInterpolatingFunction, pairpts,
ptdist];
(*Define the InverseInterpolatingFunction written by David Park, MathGroup \
[mg29221] Inverse Interpolating Functions*)
InverseInterpolatingFunction[(f_InterpolatingFunction)? (And @@ Positive /@
\
ListConvolve[{1, -1}, #1[[4, 2]]] & )] :=
InterpolatingFunction[{{First[f[[4, \
2]]], Last[f[[4, 2]]]}}, f[[2]], {f[[4, 2]]}, {f[[4, 1]], f[[3, 1]]}] ;
(*Define the function pairpts[pts] to pair up neighboring points in the
point \
data list, pts.*)
pairpts[pts_List] := Drop[Transpose[{pts, RotateLeft[pts]}], -1];
(*Define the
function ptdist that calculates the distance between \
two points.*)
ptdist = Sqrt[(#[[1, 1]] - #[[2, 1]])^2 + (#[[1, 2]] - #[[2, 2]])^2] & ;
Options[SplineEvenSamplePoints] = {
NumberOfSamplePoints -> 100};
(*++++++++++++++++ SplineEvenSamplePoints ++++++++ + +++++++++++++++++*)
(*SplineEvenSamplePoints[spline,
opts] outputs a
list of sample points that are evenly spaced along the
pathlength \
of the spline function, spline. The option NumberOfSamplePoints -> 100 can
be \
used to set the number of desired output points.*)
SplineEvenSamplePoints[spline_SplineFunction, opts___Rule] := Module[{
u1, u2, du, num, cdata, dsdu, su, sofu, uofs},
(*Extract spline parameter domain limits*)
{u1, u2} = spline[[2]];
(*Define number of
points to sample along spline curve, making sure the number
is at least 2.*)
num = NumberOfSamplePoints /. {opts} /. Options[SplineEvenSamplePoints];
num = Max[{num - 1, 3}];
(*Spline parameter sampling size*)
du = (u2 - u1)/num;
(*Generate data sampled with spline parameter*)
cdata = Table[spline[u], {u, u1, u2, du}];
(*Generate a list containing the distance between sampled points. This is
an \
approximation of the sampled derivative of the pathlength, s, against the \
spline parameter, u.*)
dsdu = ptdist /@ pairpts[cdata];
(*Generate a list of accumulated
pathlength; i.e.,
the approximate sampled cumulative integral of dsdu.*)
su = FoldList[Plus, 0, dsdu];
(*Interpolate the
cumulative pathlength data indexed by the sampled u value.*)
sofu = Interpolation[
Transpose[{Table[u, {u, u1, u2, du}],
(su/Max[su])}]
];
(*Invert the interpolation to obtain u as a function of s, the
pathlength.*)
uofs = InverseInterpolatingFunction[sofu] ;
(*Evenly sample the spline curve using the
pathlength function as the spline parameter.*)
Table[spline[Chop[uofs[s]]], {s, 0, 1, num^-1}]
];
(********************End Functions Definition************************)
(********************Begin Simple Example************************)
(*Four points in a plane*)
pts1 = {{-0.351628, 0.235057}, {2.67602, 2.61574}, {2.13537,
2.69072}, {2.5061, 1.13484}};
(*Plot points*)
ListPlot[pts1,
PlotJoined -> True,
AspectRatio -> Automatic];
(*Load SplineFit*)
<< NumericalMath`SplineFit`;
(*Fit points to cubic spline*)
spline1 = SplineFit[pts1, Cubic];
(*Plot the spline1 curve*)
ParametricPlot[spline1[u], {u, 0, 3},
PlotRange -> All,
Compiled -> False,
AspectRatio -> Automatic];
(*Generate a list of 100 evenly spaced points alone the pathlength of the
spline1 curve.*)
ex1 = SplineEvenSamplePoints[spline1];
(*Plot the points*)
ListPlot[ex1]
(********************End Simple Example************************)
____________________________________________
James T. Murray, Ph.D.
Director of R&D
Lite Cycles, Inc.
2301 N. Forbes blvd. ste. 111
Tucson, AZ 85745
office: (520) 798-0652
mobile: (520) 991-5291
fax: (520) 798-0667
web site: www.litecycles.com
____________________________________________