4
$\begingroup$

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?

2 Answers 2

1

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.

  • 0
    If my computations are right, these vectors are the same... – 2012-11-04
0

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.