Take the 2-minute tour ×
Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

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?

share|improve this question
add comment

2 Answers

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.

share|improve this answer
    
If my computations are right, these vectors are the same... –  Zsolt Nov 4 '12 at 3:46
add comment

The Frenet-Serret formulas assume that the curve is parametrised by arc length. In part, this is easy to fix up. If you have a parameter p, then, for example

dT/dp = ds/dp * dT/ds

and

ds/dp = sqrt( de/dp * de/dp + dn/dp * dn/dp + db/dp * db/dp)

However the tricky part will be reparametrising the curvature and torsion functions. That is, if you have these functions as a function of arc-length, you will need to changethem to be functions of your new parameter p.

share|improve this answer
add comment

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.