Skip to main content
added 154 characters in body
Source Link
parametrizeCurve[pts_List, a : (_?NumericQ) : 1/2] /; MatrixQ[pts, NumericQ] :=   FoldList[Plus, 0, Normalize[(Norm /@ Differences[pts])^a, Total]] /; MatrixQ[pts, NumericQ] 
periodicSplineSlopes[pts_?MatrixQ] := Module[{n = Length[pts], dy, ha, xa, ya}, {xa, ya} = Transpose[pts];   ha = {##, #1} & @@ Differences[xa];   dy = ({##, #1} & @@ Differences[ya])/ha; dy = LinearSolve[   dy = SparseArray[LinearSolve[SparseArray[{Band[{2, 1}] -> Drop[ha, 2],   Band[{1, 1}] -> ListConvolve[{2, 2}, ha],   Band[{1, 2}] -> Drop[ha, -2],   {1, n - 1} -> ha[[2]], {n - 1, 1} -> ha[[-2]]}],   3 MapThread[Dot[#1, Reverse[#2]] &,  Partition[#, 2, 1] & /@ {ha, dy}]];   Prepend[dy, Last[dy]]] 
sc = Table[{Sin[t], Cos[t], Cos[t] Sin[t]}, {t, 0, 2 π, π/5}] // N; tvals = parametrizeCurve[sc]   {0, 0.102805, 0.196242, 0.303758, 0.397195, 0.5, 0.602805, 0.696242,  0.803758, 0.897195, 1.} cmps = Transpose[sc]; slopes = periodicSplineSlopes[Transpose[{tvals, #}]] & /@ cmps; {f1, f2, f3} = MapThread[Interpolation[Transpose[{List /@ tvals, #1, #2}],   InterpolationOrder -> 3, Method -> "Hermite",   PeriodicInterpolation -> True] &, {cmps, slopes}]; 
parametrizeCurve[pts_List, a : (_?NumericQ) : 1/2] := FoldList[Plus, 0, Normalize[(Norm /@ Differences[pts])^a, Total]] /; MatrixQ[pts, NumericQ] 
periodicSplineSlopes[pts_?MatrixQ] := Module[{n = Length[pts], dy, ha, xa, ya}, {xa, ya} = Transpose[pts]; ha = {##, #1} & @@ Differences[xa]; dy = ({##, #1} & @@ Differences[ya])/ha; dy = LinearSolve[   SparseArray[{Band[{2, 1}] -> Drop[ha, 2], Band[{1, 1}] -> ListConvolve[{2, 2}, ha], Band[{1, 2}] -> Drop[ha, -2], {1, n - 1} -> ha[[2]], {n - 1, 1} -> ha[[-2]]}], 3 MapThread[Dot[#1, Reverse[#2]] &, Partition[#, 2, 1] & /@ {ha, dy}]]; Prepend[dy, Last[dy]]] 
sc = Table[{Sin[t], Cos[t], Cos[t] Sin[t]}, {t, 0, 2 π, π/5}] // N; tvals = parametrizeCurve[sc] {0, 0.102805, 0.196242, 0.303758, 0.397195, 0.5, 0.602805, 0.696242, 0.803758, 0.897195, 1.} cmps = Transpose[sc]; slopes = periodicSplineSlopes[Transpose[{tvals, #}]] & /@ cmps; {f1, f2, f3} = MapThread[Interpolation[Transpose[{List /@ tvals, #1, #2}], InterpolationOrder -> 3, Method -> "Hermite", PeriodicInterpolation -> True] &, {cmps, slopes}]; 
parametrizeCurve[pts_List, a : (_?NumericQ) : 1/2] /; MatrixQ[pts, NumericQ] :=   FoldList[Plus, 0, Normalize[(Norm /@ Differences[pts])^a, Total]] 
periodicSplineSlopes[pts_?MatrixQ] := Module[{n = Length[pts], dy, ha, xa, ya}, {xa, ya} = Transpose[pts];   ha = {##, #1} & @@ Differences[xa];   dy = ({##, #1} & @@ Differences[ya])/ha; dy = LinearSolve[SparseArray[{Band[{2, 1}] -> Drop[ha, 2],   Band[{1, 1}] -> ListConvolve[{2, 2}, ha],   Band[{1, 2}] -> Drop[ha, -2],   {1, n - 1} -> ha[[2]], {n - 1, 1} -> ha[[-2]]}],   3 MapThread[Dot[#1, Reverse[#2]] &,  Partition[#, 2, 1] & /@ {ha, dy}]];   Prepend[dy, Last[dy]]] 
sc = Table[{Sin[t], Cos[t], Cos[t] Sin[t]}, {t, 0, 2 π, π/5}] // N; tvals = parametrizeCurve[sc]   {0, 0.102805, 0.196242, 0.303758, 0.397195, 0.5, 0.602805, 0.696242,  0.803758, 0.897195, 1.} cmps = Transpose[sc]; slopes = periodicSplineSlopes[Transpose[{tvals, #}]] & /@ cmps; {f1, f2, f3} = MapThread[Interpolation[Transpose[{List /@ tvals, #1, #2}],   InterpolationOrder -> 3, Method -> "Hermite",   PeriodicInterpolation -> True] &, {cmps, slopes}]; 
deleted 33 characters in body
Source Link
sc = Table[{Sin[t], Cos[t], Cos[t] Sin[t]}, {t, 0, 2 π, π/5}] // N; tvals = parametrize[sc]parametrizeCurve[sc] {0, 0.102805, 0.196242, 0.303758, 0.397195, 0.5, 0.602805, 0.696242, 0.803758, 0.897195, 1.} {p1, p2, p3}cmps = Transpose[sc]; {d1, d2, d3}slopes = periodicSplineSlopes[Transpose[{tvals, #}]] & /@ {p1, p2, p3};cmps; {f1, f2, f3} = MapThread[Interpolation[Transpose[{List /@ tvals, #1, #2}], InterpolationOrder -> 3, Method -> "Hermite", PeriodicInterpolation -> True] &,   {{p1, p2, p3}, {d1, d2cmps, d3}slopes}]; 
sc = Table[{Sin[t], Cos[t], Cos[t] Sin[t]}, {t, 0, 2 π, π/5}] // N; tvals = parametrize[sc] {0, 0.102805, 0.196242, 0.303758, 0.397195, 0.5, 0.602805, 0.696242, 0.803758, 0.897195, 1.} {p1, p2, p3} = Transpose[sc]; {d1, d2, d3} = periodicSplineSlopes[Transpose[{tvals, #}]] & /@ {p1, p2, p3}; {f1, f2, f3} = MapThread[Interpolation[Transpose[{List /@ tvals, #1, #2}], InterpolationOrder -> 3, Method -> "Hermite", PeriodicInterpolation -> True] &,   {{p1, p2, p3}, {d1, d2, d3}}]; 
sc = Table[{Sin[t], Cos[t], Cos[t] Sin[t]}, {t, 0, 2 π, π/5}] // N; tvals = parametrizeCurve[sc] {0, 0.102805, 0.196242, 0.303758, 0.397195, 0.5, 0.602805, 0.696242, 0.803758, 0.897195, 1.} cmps = Transpose[sc]; slopes = periodicSplineSlopes[Transpose[{tvals, #}]] & /@ cmps; {f1, f2, f3} = MapThread[Interpolation[Transpose[{List /@ tvals, #1, #2}], InterpolationOrder -> 3, Method -> "Hermite", PeriodicInterpolation -> True] &, {cmps, slopes}]; 
Source Link

It does seem that the options PeriodicInterpolation -> True and Method -> "Spline" are incompatible, so I'll give a method for implementing a genuine cubic periodic spline for curves. First, let's talk about parametrizing the curve.

Eugene Lee, in this paper, introduced what is known as centripetal parametrization that can be used when one wants to interpolate across an arbitrary curve in $\mathbb R^n$. Here's a Mathematica implementation of his method:

parametrizeCurve[pts_List, a : (_?NumericQ) : 1/2] := FoldList[Plus, 0, Normalize[(Norm /@ Differences[pts])^a, Total]] /; MatrixQ[pts, NumericQ] 

The default setting of the second parameter for parametrizeCurve[] gives the centripetal parametrization. Other popular settings include a == 0 (uniform parametrization) and a == 1 (chord length parametrization).

Now we turn to generating the derivatives needed for periodic cubic spline interpolation. This is done through the solution of an appropriate cyclic tridiagonal system, for which the functions LinearSolve[], SparseArray[], and Band[] come in handy:

periodicSplineSlopes[pts_?MatrixQ] := Module[{n = Length[pts], dy, ha, xa, ya}, {xa, ya} = Transpose[pts]; ha = {##, #1} & @@ Differences[xa]; dy = ({##, #1} & @@ Differences[ya])/ha; dy = LinearSolve[ SparseArray[{Band[{2, 1}] -> Drop[ha, 2], Band[{1, 1}] -> ListConvolve[{2, 2}, ha], Band[{1, 2}] -> Drop[ha, -2], {1, n - 1} -> ha[[2]], {n - 1, 1} -> ha[[-2]]}], 3 MapThread[Dot[#1, Reverse[#2]] &, Partition[#, 2, 1] & /@ {ha, dy}]]; Prepend[dy, Last[dy]]] 

Using Sjoerd's example:

sc = Table[{Sin[t], Cos[t], Cos[t] Sin[t]}, {t, 0, 2 π, π/5}] // N; tvals = parametrize[sc] {0, 0.102805, 0.196242, 0.303758, 0.397195, 0.5, 0.602805, 0.696242, 0.803758, 0.897195, 1.} {p1, p2, p3} = Transpose[sc]; {d1, d2, d3} = periodicSplineSlopes[Transpose[{tvals, #}]] & /@ {p1, p2, p3}; {f1, f2, f3} = MapThread[Interpolation[Transpose[{List /@ tvals, #1, #2}], InterpolationOrder -> 3, Method -> "Hermite", PeriodicInterpolation -> True] &, {{p1, p2, p3}, {d1, d2, d3}}]; 

Plot the space curve:

Show[ParametricPlot3D[{f1[u], f2[u], f3[u]}, {u, 0, 1}], Graphics3D[{AbsolutePointSize[6], Point[sc]}]] 

interpolated space curve

Individually plotting the components and their respective derivatives verify that the interpolating functions are $C^2$, as expected of a cubic spline:

plots of components and derivatives