reconstruct a 3D curve from discrete curvature and torsion

Question

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.

Answer 1

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 &,
    s0,
    Transpose[{κ, τ}]
    ]
   ];

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