# Generating evenly spaced points on a curve

In the KnotData package a simple command such as

points = Table[KnotData[{3, 1}, "SpaceCurve"][t], {t, 0, 2 Pi, 0.1}];


will generate a series of points that comprise a knot (here, a trefoil, if successive points are connected).

However, the distance between points is not constant; is there a way to generate evenly spaced points? (I.e. evenly spaced arc lengths between points).

• It won't be easy; the arclength function isn't the simplest thing... Jul 18, 2012 at 15:36
• KnotData[{3, 1}, "SpaceCurve"]will give you the space curve as a pure function, e.g. {Sin[#1] + 2 Sin[2 #1], Cos[#1] - 2 Cos[2 #1], -Sin[3 #1]} &. You should be able to at least numerically determine the arc length to get your points... Jul 18, 2012 at 15:45
• Welcome to Mathematica.SE! You can consider asking future Mathematica-related questions here instead of StackOverflow. As you can see, answers come quickly and in abundance here ;-) Jul 18, 2012 at 16:40
• thanks for the suggestions! and will do, i didn't know this existed, its very useful =] Jul 18, 2012 at 17:08
• This is (mathematically) related: mathematica.stackexchange.com/a/2894/12 Jul 18, 2012 at 18:09

I'm going to brute force it numerically.

First, let's define the function we're interested in:

fun = KnotData[{3, 1}, "SpaceCurve"]


Imagine that this function fun[t] describes the position of a moving point in time. The the magnitude of its velocity as a function of the time t is

Sqrt[#.#] & [fun'[t]]


I'm going to make an interpolating function out of this to speed up and simplify the numerical calculations:

v = FunctionInterpolation[Sqrt[#.#] & [fun'[t]], {t, 0, 2 Pi}]


The integral (antiderivative) of this function will give us the distance covered as a function of time.

dist = Derivative[-1][v]


This is a monotonically increasing function, so it can be inverted:

invdist = InverseFunction[dist]


Using the inverse we can generate the times at which the point passes through the equally spaced points. Let's divide the curve into 20 equal-length segments:

times = Table[invdist[x], {x, 0, dist[2 Pi], dist[2 Pi]/20}]


Now we can easily plot the equally spaced points:

Show[
ParametricPlot3D[fun[t], {t, 0, 2 Pi}, PlotStyle -> Black],
ListPointPlot3D[fun /@ times, PlotStyle -> Directive[PointSize[0.02], Red]]
]


• Nicely explained, too! Jul 18, 2012 at 15:59
• Very well done!
– user21
Jul 18, 2012 at 16:00
• didn't know you could use Derivative[-1] for doing integrals of pure functions, nice! Jul 18, 2012 at 21:23
• @ThiesHeidecke It was essential to this solution that Integrate can handle interpolating functions. Derivative[-1] uses Integrate, but as you noted it's much more convenient to use on pure functions. Jul 18, 2012 at 21:50
• In terms of points being equidistant (making evenly spaced arc lengths), what does ParametricPlot3D[f[t], {t, 0, 2 Pi}, Mesh -> 20, MeshFunctions -> (#4 &)] produce differently? Because I can Cases them out there.
– BoLe
Mar 25, 2013 at 13:30

Just wanted to update everyone that things are much simpler, - there is built in support for this:

MeshFunctions -> {"ArcLength"}


So for our case:

Show[{
ParametricPlot3D[KnotData[{3, 1}, "SpaceCurve"][t], {t, 0, 2 Pi},

(* the trick *)
Mesh -> 15,
MeshFunctions -> {"ArcLength"},

(* styles *)
MeshStyle -> Directive[Red, PointSize[.03]],

(* the boundary split point *)
Graphics3D[{Red, PointSize[.03],
Point[KnotData[{3, 1}, "SpaceCurve"][0]]}]}]


• this certainly does make things simpler...thank you for update :) May 26, 2015 at 1:25
• @Vitaliy Probably you'll be interested to learn about important feature of MeshFunctions -> {"ArcLength"} that I recently discovered. Feb 11, 2019 at 23:40

Here's one way to implement Yves's suggestion:

(* arclength function *)
trefarc = \[FormalS] /. First[NDSolve[
{\[FormalS]'[t] == Norm[KnotData[{3, 1}, "SpaceCurve"]'[t]], \[FormalS][0] == 0},
\[FormalS], {t, 0, 2 Pi}, Method -> "Extrapolation"]]

(* length of trefoil *)
end = trefarc[2 Pi];

With[{n = 25}, (* n - number of points to generate *)
pts = (KnotData[{3, 1}, "SpaceCurve"][\[FormalT]] /.
FindRoot[trefarc[\[FormalT]] == #,
{\[FormalT], Rescale[#, {0, end}, {0, 2 Pi}], 0, 2 Pi}]) & /@
(end*Range[0, 1, 1/n])]

Graphics3D[{Line[pts], AbsolutePointSize[6], Point /@ pts}]


In case the code given above was not sufficiently transparent to readers:

1. The arclength function trefarc is produced through the use of Mathematica's numerical differential equation solver NDSolve[] (thus yielding an InterpolatingFunction[] suitable for evaluation); recall that the function $s(t)=\int_0^t g(u)\;\mathrm du$ is the solution to the initial value problem $s^\prime(t)=g(t)$, where $g(t)$ is the norm of the derivative of your (plane or space) curve with respect to the parameter.

2. We can use FindRoot[] for inverting monotonically increasing functions like the arclength function $s(t)$. I used the bracketed form of FindRoot[], FindRoot[eqn, {var, start, min, max}] to ensure that the algorithms within FindRoot[] do not go outside the domain of the InterpolatingFunction[]. For the starting point, since the arclength function looks almost like a line, I used Rescale[] to produce the inverse function of the line joining $(0,0)$ and $(2\pi,s(2\pi))$.

Taking a page from Szabolcs's code, the snippet for generating pts can be shortened a fair bit:

With[{n = 25},
pts = Composition[KnotData[{3, 1}, "SpaceCurve"],
InverseFunction[trefarc]] /@ (end*Range[0, 1, 1/n])]


This hinges on the fact that the InterpolatingFunction[] output by NDSolve[] can be inverted by InverseFunction[].

• I got old and slow ... beaten by 6 seconds ... Jul 18, 2012 at 15:54
• This answer thread more or less exploded ;-) Jul 18, 2012 at 15:55
– user21
Jul 18, 2012 at 16:02
• FWIW: one could have removed the Method -> "Extrapolation" part in the NDSolve[] for producing the arclength function, but I quite like using Bulirsch-Stoer myself... Jul 18, 2012 at 16:05
• a great and clear answer. thanks for the help!! Jul 18, 2012 at 17:07

Minimalist:

k[t_]   := KnotData[{3, 1}, "SpaceCurve"][t];
len[r_] := N@Integrate[Total[D[k[t], t]^2], {t, 0, r}];
pts[n_] := Solve[len[t1] == len[2 Pi] #/n, t1, Reals] & /@ Range[n + 1];

Graphics3D@Tube[k[t1] /. # & /@ (Flatten@pts@30), .1]


This numerical approach is based on length of chord (not arc), so it is a good approximation as long as the curve is smooth and you have close points. You have a parametric curve f of variable t. Define a numerical function that given ti finds such tf that Norm[f[tf]-f[ti]] stays constant. You need 2 functions to keep both negative and positive roots which will show different directions of going arond the knot and will draw knot's different halves.

falMax[ti_] := Max[t /. NSolve[Norm[KnotData[{3, 1}, "SpaceCurve"][t] -
KnotData[{3, 1}, "SpaceCurve"][ti]] == .2, t]] // Quiet

falMin[ti_] := Min[t /. NSolve[Norm[KnotData[{3, 1}, "SpaceCurve"][t] -
KnotData[{3, 1}, "SpaceCurve"][ti]] == .2, t]] // Quiet


Now the only thing you need is to nest that function - recursively feed it its own output, so it goes in even steps around the knot. Amount of nesting steps you get approximately after a few tries. You do not want to over-draw the extra steps.

Show[Graphics3D[{Orange, Specularity[White, 20],

Sphere[#, .1] & /@ (KnotData[{3, 1}, "SpaceCurve"] /@
NestList[falMax[.5], 0, 30]),

{Sphere[#, .1] & /@ (KnotData[{3, 1}, "SpaceCurve"] /@
NestList[falMin[.5], 0, 30])}}],

ParametricPlot3D[KnotData[{3, 1}, "SpaceCurve"][t], {t, 0, 2 Pi},
PlotStyle -> Directive[Green, Thick]]]


• Chord-length parametrization, as used here by Vitaliy, is a very popular method for generating almost-equispaced points when computing the arclength function is impractical. Sometimes, one also uses the centripetal parametrization for the purpose. See E.T.Y. Lee's paper for more details. Jul 20, 2012 at 8:32
• Very equi points, indeed! I don't think one needs to worry about nesting steps. NestWhileList[falMax, 0, #2 > #1 &, 2]
– BoLe
Apr 12, 2013 at 23:20

Here is another way to get points equally spaced by chord-length. It will give a good result if the chord length is relatively small compared to the maximum radius of curvature along the curve. (If, say, the chord length is no greater than the maximum radius, then between successive points, the turning will be less 60 deg., and the difference between the chord and arc lengths will be less than 5%).

The method is fast because the objective function, which is the sum of the squares of the differences in the lengths of adjacent chords, is quick to be minimized (provided the initial values are ordered). Another reason is that the speed of this particular parametrization does not vary much, which leads to excellent initial points.

The first and last point are set equal to the end points of the curve. The code below assumes that the parametrization is periodic. The end points are assumed to be equal and one of them is discarded.

xFN = KnotData[{3, 1}, "SpaceCurve"];
nPts = 80;                                                (* number of points *)

Block[{t},
t[1] = 0; t[nPts + 1] = 2 Pi;                            (* initialize end points *)
sols = Array[t, nPts] /.
Last@FindMinimum[
Total[Differences[Norm /@ Differences[xFN /@ Array[t, nPts + 1]]]^2],
Table[{t[i], (i - 1)/(nPts) 2 Pi}, {i, 2, nPts}]];   (* omit end points *)
]

Show[
Graphics3D[{Red, Translate[Sphere[{0, 0, 0}, 0.08], xFN /@ sols]}],
ParametricPlot3D[xFN[t], {t, 0, 2 Pi}]
]


P.S. If you desire relatively few points so that the chord length is greater than the radius of curvature, then generate many points and down-sample. To get the 80 points above took slightly over 1/20 sec. One can get 8 equally spaced points from the above calculation by using sols[[;; ;; 10]]:

Show[Graphics3D[{Red, Translate[Sphere[{0, 0, 0}, 0.08], xFN /@ sols[[;; ;; 10]]]}],
ParametricPlot3D[xFN[t], {t, 0, 2 Pi}]]