Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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).

share|improve this question

migrated from stackoverflow.com Jul 18 '12 at 15:33

This question came from our site for professional and enthusiast programmers.

    
It won't be easy; the arclength function isn't the simplest thing... –  J. M. Jul 18 '12 at 15:36
2  
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... –  Yves Klett Jul 18 '12 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 ;-) –  Szabolcs Jul 18 '12 at 16:40
    
thanks for the suggestions! and will do, i didn't know this existed, its very useful =] –  scallionpancake Jul 18 '12 at 17:08
    
This is (mathematically) related: mathematica.stackexchange.com/a/2894/12 –  Szabolcs Jul 18 '12 at 18:09

5 Answers 5

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]]
 ]

Mathematica graphics

share|improve this answer
    
Nicely explained, too! –  Yves Klett Jul 18 '12 at 15:59
    
Very well done! –  user21 Jul 18 '12 at 16:00
    
thanks for the answer! =] –  scallionpancake Jul 18 '12 at 17:07
5  
didn't know you could use Derivative[-1] for doing integrals of pure functions, nice! –  Thies Heidecke Jul 18 '12 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. –  Szabolcs Jul 18 '12 at 21:50

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}]

equispaced trefoil points


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[].

share|improve this answer
1  
I got old and slow ... beaten by 6 seconds ... –  Szabolcs Jul 18 '12 at 15:54
    
This answer thread more or less exploded ;-) –  Yves Klett Jul 18 '12 at 15:55
    
Nice answer!... –  user21 Jul 18 '12 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... –  J. M. Jul 18 '12 at 16:05
    
a great and clear answer. thanks for the help!! –  scallionpancake Jul 18 '12 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]

enter image description here

share|improve this answer

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]]]

enter image description here

share|improve this answer
1  
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. –  J. M. Jul 20 '12 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 '13 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}]
 ]

Mathematica graphics

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}]]
share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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