Mercurial > repos > public > sbplib

comparison +parametrization/dataSpline.m @ 1091:b4054942e277 feature/dataspline

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Rewrite dataSpline() avoiding the spline function in Curve and using fnder for the differentiation
author Jonatan Werpers <jonatan@werpers.com>
date Fri, 05 Apr 2019 10:17:04 +0200
parents d7f6c10eab13
children 47a72344db71
comparison
equal deleted inserted replaced
1090:3c187e639dfa 1091:b4054942e277
9 function C = dataSpline(t_data, f_data) 9 function C = dataSpline(t_data, f_data)
10 10
11 assert(length(t_data)==length(f_data),'Vectors must be same length'); 11 assert(length(t_data)==length(f_data),'Vectors must be same length');
12 m_data = length(t_data); 12 m_data = length(t_data);
13 13
14 % Create spline interpolant 14 pp_g = spapi(4, t_data, f_data); % equivalent to g = spapi(aptknt(t_data, 4), t_data, f_data)
15 f = parametrization.Curve.spline(t_data, f_data); 15 % or (not sure what the difference is?!)
16 % g = spapi(optknt(t_data, 4), t_data, f_data)
17 pp_gp = fnder(g);
16 18
17 % Reparametrize with a parameter s in [0, 1] 19 g = @(t)fnval(pp_g, t);
18 tmin = min(t_data); 20 pp_gp = @(t)fnval(pp_gp, t);
19 tmax = max(t_data);
20 t = @(s) tmin + s*(tmax-tmin);
21 21
22 % Create parameterized curve
23 g = @(s) [t(s); f(t(s))];
24
25 % Compute numerical derivative of curve using twice as many points as in data set
26 m = 2*m_data;
27 ops = sbp.D2Standard(m, {0, 1}, 6);
28 gp = parametrization.Curve.numericalDerivative(g, ops.D1);
29
30 % Create curve object
31 C = parametrization.Curve(g, gp); 22 C = parametrization.Curve(g, gp);
32 23
33 % Reparametrize with arclength parametrization 24 % Reparametrize with arclength parametrization
34 C = C.arcLengthParametrization(m_data); 25 C = C.arcLengthParametrization(m_data);
35
36 % To avoid nested function calls, evaluate curve and compute final spline.
37 tv = linspace(0, 1, m_data);
38 gv = C.g(tv);
39 g = parametrization.Curve.spline(tv, gv);
40 C = parametrization.Curve(g);
41
42 end 26 end