I tried to reconstruct a 3D curve with given curvature and torsion. I saw some threads talking about using runge kutta. However, as far as I see that they required curvature and torsion were continuous. Unfortunately, my curvature and torsion are discrete. They are not functions of arc-length, but indexed by arc-length. Anyone can help to see how to solve the problem with discrete curvature and torsion. Thank you.

Hi here is an example. Assume we have a trajectory indexed with arc-length of $\delta s$ and total samples of 500. I have

$$y_{s} = \begin{bmatrix} r_{s}\\ T_{s}\\ N_{s}\\ B_{s} \end{bmatrix}$$

where $s = 1, 2, ..., 500$.

Now we can create a function used by runge-kutta as

$$y_{s}^{'} = \begin{bmatrix} 0 & 1 & 0 & 0\\ 0 & 0 & \kappa_{s} & 0\\ 0 & -\kappa_{s} & 0 & \tau_{s}\\ 0 & 0 & -\tau_{s} & 0 \end{bmatrix}y_{s}$$

Note that here $\kappa_{s}$ and $\tau_{s}$ are not function of $s$. Instead, they are indexed by $s$. Therefore, if $s = 1.5$, we are not able to compute $\kappa_{s}$ and $\tau_{s}$.

In terms of using runge-kutta, for

$$k_{1} = \delta s f(y_{s}, x_{s})$$

it is fine since $x_{s}/\delta s$ always is an integer. However for

$$k_{2} = \delta sf(y_{s}+\frac{k_{1}}{2}, x_{s}+\frac{\delta s}{2})$$

$\frac{x_{s}+\frac{\delta s}{2}}{\delta s}$ cannot be an integer, therefore we cannot find corresponding $\kappa$ and $\tau$ in this case.

Not sure if I explained my problem clearly.Thank you.

  • 1
    $\begingroup$ Welcome! If you were to add some (sample) data, that would make answering so much more comfortable! $\endgroup$
    – Yves Klett
    Sep 10 '14 at 7:06
  • $\begingroup$ It seems to me that Interpolate with NDSolve might solve the problem. But I can't say more without the data. Can you edit your post and add some example? $\endgroup$
    – ybeltukov
    Sep 10 '14 at 9:40
  • $\begingroup$ Related: Finding a 3d curve from torsion and curvature with NDSolve $\endgroup$
    – Michael E2
    Sep 10 '14 at 12:34

The curvature and torsion are rates of turning of the Frenet-Serret frame and can be used to integrate the frame using an Euler-type method. The unit tangent vector of the frame is the velocity and can be used to integrate the position.

Set up some initial data: the initial Frenet-Serret frame (randomly chosen below), initial point s0, the change in arclength per step ds, curvature data κ, and torsion data τ.

frame0 = If[Det[#] < 0, #[[{2, 1, 3}]], #] &@ Orthogonalize@RandomReal[{-1, 1}, {3, 3}];
s0 = {0., 0., 0.};
ds = 0.1;
κ = Table[3 Sin[ s/10.]^2, {s, 200}];
τ = Table[Sin[ s/10.], {s, 200}];

The unit tangent vector is estimated as the bisector of the tangent vector before and after rotating the frame. The order of rotation does make a difference, but it is of second-order. I ignored it, but one might think about how to incorporate it to improve the estimate.

path = Module[{frame = frame0},
   FoldList[#1 + Normalize[                                    (* #1 == current position *)
        frame[[All, 1]] +
        (frame = 
            RotationMatrix[Last[#2] ds, frame[[All, 1]]] .     (* Last[#2] == τ *)
             RotationMatrix[First[#2] ds, frame[[All, 3]]] .   (* First[#2] == κ *)
             frame)[[All, 1]]] ds &,
    Transpose[{κ, τ}]

  Line[path, VertexColors -> ColorData["Rainbow"] /@ Rescale@Range@Length@path]}]

Mathematica graphics


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.