[Math] Reconstructing space curves from its curvature and torsion

curvaturedifferential-geometryMATLABoctaveplane-curves

I have to write a program which gets two functions (curvature and torsion) and 3 vectors of the Frenet-Serret "frame" at the starting point – and I have to reconstruct the space curve from this given input data.

So I've written some Octave code, which works fine, for example on the curve

 Phi(t) = (cos(t/sqrt(2)), sin(t/sqrt(2)), t/sqrt(2))  % natural parametrization

(it's some kind of spiral) and I was happy because when I plotted results, the original curve was almost the same as the reconstructed one. But when I tried out some other parametrization, there was huge difference between the reconstructed and the original one:

 Phi(t) = (cos(t), sin(t), t)  % other parametrization of the same curve

I can't really find what I did wrong in my code… Maybe my idea was completely wrong, and I can only reconstruct curves which is naturally parametrized?

% note: this returns a vector
% K : curvature
% T : torsion
function equation = f(s, gvect, K, T, velocity)
    gs = toStruct(gvect);
    e_v = velocity(s)*K(s).* gs.n;
    n_v = -velocity(s)*K(s).* gs.e + velocity(s)*T(s).* gs.b;
    b_v = -velocity(s)*T(s).* gs.n; 
    equation = [e_v; n_v; b_v];
    return;
endfunction

% returns a struct
function str = toStruct(vect)
    str.e = vect(1:3);
    str.n = vect(4:6);
    str.b = vect(7:9);
endfunction

% returns a vector
function vect = toVect(str)
    vect = [str.e; str.n; str.b];
endfunction

% numerical integration
function p = nintegrate(x, fx)
    acc = 0;
    for i = [1:length(x)-1]
            p(i) = acc;
            delta_x = x(i+1) - x(i);
            acc += delta_x * fx(i);
    endfor
endfunction

% K : curvature
% T : torsion
function gvals = RungeKutta4(S, g, K, T, velocity, delta)
    for i = [1:length(S)-1]
            g_i = toVect(g(i));
            k1 = f(S(i), g_i, K, T, velocity);
            k2 = f(S(i)+delta/2, g_i+delta/2*k1, K, T, velocity);
            k3 = f(S(i)+delta/2, g_i+delta/2*k2, K, T, velocity);
            k4 = f(S(i)+delta, g_i+delta*k3, K, T, velocity);
            g(i+1) = toStruct( g_i + delta/6*(k1+2*k2+2*k3+k4) );
    endfor
    gvals = g;
endfunction

function drawPlots(xs, ys, zs, orig_xs, orig_ys, orig_zs)
    plot3(xs, ys, zs, ";reconstructed;", orig_xs, orig_ys, orig_zs, ";original;");
endfunction

function Phi = reconstruct(K, T, e0, n0, b0, alpha, beta, start, delta = 0.1, velocity = @()1)
    S = [alpha:delta:beta];
    g(1).e = e0;    % unit tangent vector 
    g(1).n = n0;    % principal normal vector 
    g(1).b = b0;    % binormal vector
    g = RungeKutta4(S, g, K, T, velocity, delta);

    e_x = cellfun(@(v)v(1), {g.e});
    e_y = cellfun(@(v)v(2), {g.e});
    e_z = cellfun(@(v)v(3), {g.e});

    Phi.x = nintegrate(S, e_x) + start(1);
    Phi.y = nintegrate(S, e_y) + start(2);
    Phi.z = nintegrate(S, e_z) + start(3);
    return;
endfunction

% works OK for this
function spiralNaturalParametrized()
    L = sqrt(2)*4*pi;
    e0 = [ 0,  1/sqrt(2), 1/sqrt(2)]';  % unit tangent vector 
    n0 = [-1,          0,         0]';  % principal normal vector 
    b0 = [ 0, -1/sqrt(2), 1/sqrt(2)]';  % binormal vector
    Phi = reconstruct(@()1/2, @()1/2, e0, n0, b0, 0, L, [1 0 0], 0.05);

    t = [0:0.05:L];

    drawPlots(Phi.x, Phi.y, Phi.z, 
          cos(t./sqrt(2)), sin(t./sqrt(2)), t./sqrt(2));              
endfunction

% another parametrization: huge error
function spiralAnother()
    L = 4*pi;
    e0 = [ 0,  1/sqrt(2), 1/sqrt(2)]'; 
    n0 = [-1,          0,         0]'; 
    b0 = [ 0, -1/sqrt(2), 1/sqrt(2)]';
    Phi = reconstruct(@()1/2, @()1/2, e0, n0, b0, 0, L, [1 0 0], 0.05, @(v)sqrt(2));

    t = [0:0.05:L];
    drawPlots(Phi.x, Phi.y, Phi.z, cos(t), sin(t), t);
endfunction

I included my code.
Any help is appreciated!

Update:

I have created another test which reveals the fact that the algorithm has some bug inside even when it deals with the natural parametric case:

function YetAnotherTest()
    L = 10*pi;
    e0 = [       0,    2/sqrt(13),   3/sqrt(13)]'; 
    n0 = [      -1,             0,            0]'; 
    b0 = [       0,   -3/sqrt(13),   2/sqrt(13)]';
    Phi = reconstruct(@()(2/13),        % curvature
                      @()sqrt(3/13),    % torsion
                      e0, n0, b0,       % initial frame
                      0, L,             % interval
                      [2 0 0], 0.05);   % start point and step size

    t = [0:0.05:L];

    drawPlots(Phi.x, Phi.y, Phi.z, 
              2*cos(t./sqrt(13)), 2*sin(t./sqrt(13)), t*(3/sqrt(13)));
endfunction

Note: this is a test for $$Phi(t) = ( 2*\cos(t/sqrt(13)), 2*\sin(t/sqrt(13)), 3t/sqrt(13) )$$

which is a natural parametric curve.
If you take a look at this plot, you can see that the reconstructed curve is totally different than the original one.

Update2:

I have a suspicion that maybe I solve the differential equations wrong, and that causes the error. I have the differential equation system (my solver can be found in the RungeKutta4 and the f functions above):

Frenet–Serret formulas:

e'(t) =  K(t)*velocity(t)*n(t)
n'(t) = -K(t)*velocity(t)*e(t) + T(t)*velocity(t)*b(t)
b'(t) = -T(t)*velocity(t)*n(t)

given the initial values:
  e(0) = e0
  n(0) = n0
  b(0) = b0

where
K : curvature
T : torsion

Maybe I misuse Runge-Kutta…

Update3:

If I change my own Runge-Kutta algorithm to the built-in Octave function lsode and I rewrite the numerical integration code to use trapz (trapezoid rule), nothing changes. So I think the bug is not in these functions…But where?

Best Answer

Your TNB frame for the second parametrization is not going to be the same as it is for the first, but it looks like you are assuming that it will be.

Compute the Tangent, Normal and Binormal vectors for the parametrization of the helix given by $$X(t) = ( \cos(t), \sin(t), t )$$ and then set the values of L, e0, n0, b0 accordingly in your spiralAnother() function.

Related Question