# Reconstructing space curves from its curvature and torsion

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?

-

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.

-
If my computations are right, these vectors are the same... –  Zsolt Nov 4 '12 at 3:46

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.

-