5
$\begingroup$

By fitting the data using spline, I have created a curve.

sp = SplineFit[data1, Cubic]

I am trying to divide this curve into small segments of equal length. To do so, I am trying to calculate the length of the curve using:

NIntegrate[Sqrt[1 + sp'[z]^2], {z, 0, 34}]

I get:

NIntegrate::inumr: "The integrand \!\(\*SqrtBox[
RowBox[{\"1\", \"+\", SuperscriptBox[
RowBox[{SuperscriptBox[TemplateBox[{\"\\\"SplineFunction[\\\"\",\"Cubic\",\"\\\", \\\"\",RowBox[{\"{\", RowBox[{\"0.`\", \",\", \"34.`\"}], \"}\"}],\"\\\", \\\"\",\"\\\"<>\\\"\",\"\\\"]\\\"\"},
\"RowDefault\"], \"\[Prime]\",
MultilineFunction->None], \"[\", \"z\", \"]\"}], \"2\"]}]]\) has evaluated to non-numerical values for all sampling points in the region with boundaries {{0,34}}."

I am new to Mathematica and I really do not understand how splines work. I am trying to use the technique mentioned in How do I split up a curve into chords of equal length? to divide the curve into equal segments.

| improve this question | |
$\endgroup$
  • $\begingroup$ What version do you have? SplineFit seems to be somewhat deprecated. We now have Interpolation, BezierFunction, BSplineFunction -- maybe others I don't recall. The SplineFunction does not seem to be symbolically differentiable. $\endgroup$ – Michael E2 May 11 '15 at 20:03
  • $\begingroup$ I have Mathematica 9 as well as 10. I tried to Interpolate but I couldn't do so because I had to interpolate both x and y coordinates. $\endgroup$ – psimeson May 11 '15 at 20:35
  • 2
    $\begingroup$ Like this? xydata = Table[{t Cos[t], t Sin[t]}, {t, 0., 10}]; pf = Interpolation[ MapThread[{{#1}, #2} &, {Rescale@Range@Length@xydata, xydata}], Method -> "Spline"]; ParametricPlot[pf[t], {t, 0, 1}]. Note the data for Interpolation should be in the form Table[{{t}, {x, y}}, {t,...}] -- note the braces carefully. $\endgroup$ – Michael E2 May 11 '15 at 20:56
5
$\begingroup$

Using the somewhat outdated Splines package:

Needs["Splines`"]

f = SplineFit[{{0, 0}, {2, 0}, {2, 2}, {0, 1}}, Cubic];

ParametricPlot[f[t], {t, 0, 3}]

Mathematica graphics

The SplineFunction cannot be differentiated symbolically:

f'[0]
(*  (SplineFunction[Cubic, {0.,3.}, <>]^\[Prime])[0]  *)

NDSolve can construct a numerical derivative (apparently):

{arclength} = 
  NDSolve[{y[t] == f[t], s'[t] == Sqrt[y'[t].y'[t]], x'[t] == 1, 
    y[0] == f[0], x[0] == 0, s[0] == 0}, s, {t, 0, 3}];

s[3] /. arclength
(*  6.59702  *)

Plot[s[t] /. arclength, {t, 0, 3}]

Mathematica graphics

There are better ways, I think, if you have version 7 or higher. But I love how robust NDSolve is. (Probably, if you have a version below 7, it won't be robust enough to do the above. I cannot check that myself.)

Addendum

If generating evenly spaced points is the goal, then Alexey Popkov's answer to "How to get a uniformly distributed list of points from a B-spline function" is the easiest way:

plot = ParametricPlot[f[t], {t, 0, 3}, MeshFunctions -> {"ArcLength"}, Mesh -> {15}, 
 MeshStyle -> {PointSize[0.02], Red}]

Mathematica graphics

One can get the points (not the parameter) this way:

Cases[Normal @ plot, Point[p_] :> p, Infinity]
(*
  {{0.395124, -0.117757}, {0.985196, 1.74711}, {0.648418, 1.50936},
   {0.321879, 1.25766}, {0.794663, -0.219251}, {1.34214, 1.9529},
   {1.96322, -0.0415362}, {1.73477, 2.07018}, {1.61147, -0.246155},
   {2.17169, 0.311035}, {2.09348, 1.90693}, {2.25271, 1.52986},
   {1.20225, -0.278889}, {2.27279, 0.709978}, {2.29956, 1.12093}}
*)
| improve this answer | |
$\endgroup$
  • $\begingroup$ I am sorry I couldn't understand your solution. How did you come up with that differential equation? $\endgroup$ – psimeson May 11 '15 at 22:06
  • 1
    $\begingroup$ @psimeson The arclength $s$ is the integral of $\sqrt{1+(dy/dx)^2}\,dx$ or of $\sqrt{(dx/dt)^2+(dy/dt)^2}\,dt$. The first case is for a curve of the form of the graph of $y=f(z)$; the second is for a parametrized curve $x=x(t),\,y=y(t)$, which is your case. In your question you seemed to use the first form, but the second form is the correct one. Hence the s'[t] == Sqrt[y'[t].y'[t]] equation. The rest defines y = f = your Spline parametrization, and the dummy ODE system x'[t] == 1, x[0] = 0 is a standard trick to force NDSolve to integrate the system (it's the only ODE of the DAE)... $\endgroup$ – Michael E2 May 12 '15 at 0:13
  • $\begingroup$ ...(DAE = "Differential Algebraic Equation"). One could use s'[t] = Norm[s'[t]] instead of s'[t] == Sqrt[y'[t].y'[t]] and it works (at first I thought it wouldn't). While NDSolve is integrating x'[t] = 1, it will also be setting y[t] == f[t] and numerically integrating the arclength s via the equation s'[t] == Sqrt[y'[t].y'[t]]. $\endgroup$ – Michael E2 May 12 '15 at 0:18
  • 1
    $\begingroup$ @psimeson sol = NDSolve[..] gives {{s -> InterpolatingFunction[..]}}. Then s[0.5] /. sol gives {value} in a list. This is because sol is treated as a list of solutions {sol1}. To get just value we need to get rid of an extra set of braces. One way is to use s[0.5] /. First[sol]; similarly, one can use sol = First[NDSolve[..]]. Another way is to use a list assignment as I did. Note {x, y} = {1, 2} is the same as x = 1; y = 2. Similarly {arclength} = NDSolve[..] is the same as {arclength} = {{s -> ..}}, which is the same as arclength = {s ->..}. $\endgroup$ – Michael E2 May 13 '15 at 16:30
  • 1
    $\begingroup$ @psimeson See the "Introduction" part of this answer for links to introductory materials. $\endgroup$ – Michael E2 May 13 '15 at 16:33
5
$\begingroup$

If you do have version 10, another possibility is to use the new geometric computation functionality:

ArcLength@DiscretizeGraphics@BSplineCurve@data1

The above give the total length. To find the length of the segment of the spline function sp between times t1 and t2, you could extend the same approach like so: 

length[{t1_, t2_}] := ArcLength@DiscretizeGraphics@ParametricPlot[sp[t], {t, t1, t2}]

Note that this is not a good solution if you want to use the length function to find roots, as ParametricPlot will throw out errors if a non-numeric limit is supplied. Perhaps in the future, though, Mathematica will offer more options for working with geometric regions.

I also found a few things out while investigating Michael E2's excellent answer:

  1. BSplineFunctions can be symbolically differentiated, but what is produced is a little different than normal:

    f = BSplineFunction[{{1, 1}, {2, 3}, {3, -1}, {4, 1}, {5, 0}}];

df = D[f[t], t] (* -> BSplineFunction[{{0., 1.}}, <>][t] *)

    that is, it is already acting on t. We can plot it

    ParametricPlot[df, {t, 0, 1}, AspectRatio -> 1]

    df

    and find arc length with NIntegrate:

    s[tf_] := NIntegrate[Norm@df, {t, 0, tf}];

s[1] (* -> 5.30301 *)

    (Note that the equation for the arc length is $ds(t) = \sqrt{(dx/dt)^2 + (dy/dt)^2} dt $, hence the use of Norm@df.)

  2. To find the points that split a BSplineFunction into equal parts, we can use NDSolve, which is robust enough to handle the curve directly:

    findcuts[sf_, n_] := Module[
  {ti, tf, length, sol, times},
  {{ti, tf}} = sf["Domain"];
  length = NIntegrate[Norm@D[sf[t], t], {t, ti, tf}];
  {{sol}, {times}} = Reap@NDSolve[
    {r'[t] == Sqrt[sf'[t].sf'[t]],
     r[0] == 0,
     WhenEvent[Evaluate[Mod[r[t], length/n]], Sow[t]]},
    r, {t, ti, tf}];
  sf@# & /@ Drop[times, -1]];

Show[ParametricPlot[f[t], {t, 0, 1}], ListPlot[findcuts[f, 7], PlotStyle -> Orange]]

    f

    I heartily agree with Michael E2 about the versatility of NDSolve!

| improve this answer | |
$\endgroup$
  • $\begingroup$ I do have access to Mathematica 10. Thanks Virgil, the first command worked. It's nice that Bspline uses cubic spline. Is there a similar powerful way to divide this curve into segments of equal length? $\endgroup$ – psimeson May 11 '15 at 21:02
  • $\begingroup$ @psimeson, not that I know of (but many people here know more than me). Probably your best bet is to go with the method outlined in your linked question. $\endgroup$ – Virgil May 11 '15 at 21:50
  • $\begingroup$ One more question: Using ArcLength@DiscretizeGraphics@BSplineCurve@data1 I got the length to be 55.3209 and using length[{0,34}] = 56.0843 (This the contour length of the curve that I got using Image J). I don't know why these two are different. $\endgroup$ – psimeson May 11 '15 at 21:55
  • $\begingroup$ @psimeson, I get the same kind of difference. I think is stems from the difference between the way DiscretizeGraphics treats BSplineCurve and BSplineFunction. It seems that BSplineCurve, which is a graphics primitive, can only be approximated, so the quality of the discretization cannot be improved by passing options like MaxCellMeasure or MeshQualityGoal to DiscretizeGraphics. The parametric plot of the BSplineFunction can, however, be improved, so that should give a better arc length measure. $\endgroup$ – Virgil May 11 '15 at 22:14
  • $\begingroup$ @psimeson, actually, all nonlinear primitives can only be approximated by DiscretizeGraphics, so the parametric plot of the BSplineFunction cannot be improved either. However, the discretization seems to go by plot points, so it will take advantage of ParametricPlot's adaptive sampling and should be a closer fit in areas of higher curvature. $\endgroup$ – Virgil May 11 '15 at 22:34
2
$\begingroup$

Something easy to do, not very efficient though.

dz = .0000001;
arc[z_]:=NIntegrate[Norm[(sp[z + dz] - sp[z])/dz], {z, 0, 34}]
| improve this answer | |
$\endgroup$
1
$\begingroup$

In this previous answer, I gave a function that reproduces the functionality of the old NumericalMath`SplineFit`​ package, and outputs a function in BSplineFunction[] format. Unlike the older function, it plays nice with differentiation, which allows us to use it with NIntegrate[] to find the overall arclength:

sf = splinefit[{{0, 0}, {2, 0}, {2, 2}, {0, 1}}];
NIntegrate[Norm[sf'[t]], {t, 0, 3}]
   6.597018140608958

Of course, if you want an arclength function, you could use NDSolve[] in the same manner as Michael did for his answer:

s = NDSolveValue[{sp'[t] == Norm[sf'[t]], sp[0] == 0}, sp, {t, 0, 3}]
| improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.