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?