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