Mercurial > repos > public > sbplib

comparison +grid/Curve.m @ 104:ffd9e68f63cc feature/arclen-param

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Splined in the other direction! Works great. Runs faster. Removed ugly fix outside [0,1] in spline(). Camelcased numerical_derivative().
author Martin Almquist <martin.almquist@it.uu.se>
date Mon, 07 Dec 2015 18:46:59 +0100
parents bc5db54f9efd
children 1df4f3704b76
comparison
equal deleted inserted replaced
103:bc5db54f9efd 104:ffd9e68f63cc
6 6
7 end 7 end
8 methods 8 methods
9 %TODO: 9 %TODO:
10 % Errors or FD if there is no derivative function added. 10 % Errors or FD if there is no derivative function added.
11 % -semi-done
12 11
13 % Concatenation of curves 12 % Concatenation of curves
14 % Subsections of curves 13 % Subsections of curves
15 % Stretching of curve parameter - semi-done 14 % Stretching of curve parameter - semi-done
16 % Curve to cell array of linesegments 15 % Curve to cell array of linesegments
96 default_arg('N',100); 95 default_arg('N',100);
97 96
98 % Construct arcLength function using splines 97 % Construct arcLength function using splines
99 tvec = linspace(0,1,N); 98 tvec = linspace(0,1,N);
100 arcVec = obj.arcLength(0,tvec); 99 arcVec = obj.arcLength(0,tvec);
101 arcLength = spline(tvec,arcVec); 100 tFunc = spline(arcVec,tvec); % t as a function of arcLength
102 101 L = obj.arcLength(0,1);
103 % Stretch the parameter, construct function with splines 102 arcPar = @(s) tFunc(s*L);
104 arcParVec = fzero_vec(@(s)arcLength(s) - tvec*arcLength(1),[0-10*eps,1+10*eps]);
105 arcPar = spline(tvec,arcParVec);
106 103
107 % New function and derivative 104 % New function and derivative
108 g_new = @(t)obj.g(arcPar(t)); 105 g_new = @(t)obj.g(arcPar(t));
109 gp_new = @(t) normalize(obj.gp(arcPar(t))); 106 gp_new = @(t) normalize(obj.gp(arcPar(t)));
110 curve = grid.Curve(g_new,gp_new); 107 curve = grid.Curve(g_new,gp_new);
247 end 244 end
248 end 245 end
249 246
250 methods (Static) 247 methods (Static)
251 248
252 function gp_out = numerical_derivative(g,D1) 249
250 function gp_out = numericalDerivative(g,D1)
253 m = length(D1); L = 1; % Assume curve parameter from 0 to 1. 251 m = length(D1); L = 1; % Assume curve parameter from 0 to 1.
254 t = linspace(0,L,m); 252 t = linspace(0,L,m);
255 gVec = g(t)'; 253 gVec = g(t)';
256 gpVec = (D1*gVec)'; 254 gpVec = (D1*gVec)';
257 255
374 I(i) = integral(f,a(i),b(i)); 372 I(i) = integral(f,a(i),b(i));
375 end 373 end
376 end 374 end
377 end 375 end
378 376
379 function I = fzero_vec(f,lim)
380 % Wrapper around the built-in function fzero that
381 % handles multiple functions (vector-valued f).
382
383 fval = f(lim(1));
384 I = zeros(size(fval));
385
386 for i = 1:length(fval)
387 e = zeros(size(fval));
388 e(i) = 1;
389 I(i) = fzero(@(t) sum(e.*f(t)),lim);
390 end
391
392 end
393
394 % Returns a function handle to the spline. 377 % Returns a function handle to the spline.
395 function f = spline(tval,fval,spline_order) 378 function f = spline(tval,fval,spline_order)
396 default_arg('spline_order',4); 379 default_arg('spline_order',4);
397 [m,~] = size(tval); 380 [m,~] = size(tval);
398 assert(m==1,'Need row vectors.'); 381 assert(m==1,'Need row vectors.');
399 382
400 % make vectors longer to be safe slightly beyond edges.
401 dt0 = tval(2)-tval(1); dt1 = tval(end)-tval(end-1);
402 df0 = fval(2)-fval(1); df1 = fval(end)-fval(end-1);
403 tval = [tval(1)-dt0, tval, tval(end)+dt1];
404 fval = [fval(1)-df0, fval, fval(end)+df1];
405
406 f_spline = spapi( optknt(tval,spline_order), tval, fval ); 383 f_spline = spapi( optknt(tval,spline_order), tval, fval );
407 f = @(t) fnval(f_spline,t); 384 f = @(t) fnval(f_spline,t);
408 end 385 end
409 386
410 387