Mercurial > repos > public > sbplib

comparison +grid/Curve.m @ 103:bc5db54f9efd feature/arclen-param

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
fzero_vec, integral_vec and spline are now local functions in Curve. Renamed arcLengthStretch to arcLengthParametrization. Removed plot_derivative. Added some comments and extra lines + removed unneccesary lines. arcLength is now a method and not static. Constructor does not accept difference operator anymore.
author Martin Almquist <martin.almquist@it.uu.se>
date Mon, 07 Dec 2015 17:24:28 +0100
parents ecf77a50d4fe
children ffd9e68f63cc
comparison
equal deleted inserted replaced
99:ecf77a50d4fe 103:bc5db54f9efd
13 % Concatenation of curves 13 % Concatenation of curves
14 % Subsections of curves 14 % Subsections of curves
15 % Stretching of curve parameter - semi-done 15 % Stretching of curve parameter - semi-done
16 % Curve to cell array of linesegments 16 % Curve to cell array of linesegments
17 17
18 % Can supply either derivative or a difference operator D1 18 % The curve parameter t must be in [0,1].
19 function obj = Curve(g,gp,transformation,D1) 19 % Can supply derivative.
20 function obj = Curve(g,gp,transformation)
20 default_arg('gp',[]); 21 default_arg('gp',[]);
21 default_arg('transformation',[]); 22 default_arg('transformation',[]);
22 default_arg('D1',[]);
23 p_test = g(0); 23 p_test = g(0);
24 assert(all(size(p_test) == [2,1]), 'A curve parametrization must return a 2x1 vector.'); 24 assert(all(size(p_test) == [2,1]), 'A curve parametrization must return a 2x1 vector.');
25
26 if(isempty(gp) && ~isempty(D1))
27 gp = grid.Curve.numerical_derivative(g,D1);
28 end
29 25
30 if ~isempty(transformation) 26 if ~isempty(transformation)
31 transformation.base_g = g; 27 transformation.base_g = g;
32 transformation.base_gp = gp; 28 transformation.base_gp = gp;
33 [g,gp] = grid.Curve.transform_g(g,gp,transformation); 29 [g,gp] = grid.Curve.transform_g(g,gp,transformation);
51 function h = plot(obj,n) 47 function h = plot(obj,n)
52 default_arg('n',100); 48 default_arg('n',100);
53 49
54 t = linspace(0,1,n); 50 t = linspace(0,1,n);
55 p = obj.g(t); 51 p = obj.g(t);
56 h = line(p(1,:),p(2,:));
57 end
58
59 % Plots the derivative gp(t) for 0<t<1, using n points. Returns a handle h to the plotted curve.
60 % h = plot_curve(gp,n)
61 function h = plot_derivative(obj,n)
62 default_arg('n',100);
63
64 t = linspace(0,1,n);
65 p = obj.gp(t);
66 h = line(p(1,:),p(2,:)); 52 h = line(p(1,:),p(2,:));
67 end 53 end
68 54
69 function h= plot_normals(obj,l,n,m) 55 function h= plot_normals(obj,l,n,m)
70 default_arg('l',0.1); 56 default_arg('l',0.1);
95 81
96 obj.plot(); 82 obj.plot();
97 end 83 end
98 % Shows curve with name and arrow for direction. 84 % Shows curve with name and arrow for direction.
99 85
86 % Length of arc from parameter t0 to t1 (which may be vectors).
87 % Computed using derivative.
88 function L = arcLength(obj,t0,t1)
89 assert(~isempty(obj.gp),'Curve has no derivative!');
90 speed = @(t) sqrt(sum(obj.gp(t).^2,1));
91 L = integral_vec(speed,t0,t1);
92 end
100 93
101 % Arc length parameterization. 94 % Arc length parameterization.
102 function curve = arcLengthStretch(obj,N) 95 function curve = arcLengthParametrization(obj,N)
103 default_arg('N',100); 96 default_arg('N',100);
104 assert(~isempty(obj.gp),'This curve does not yet have a derivative added.')
105 97
106 % Construct arcLength function using splines 98 % Construct arcLength function using splines
107 tvec = linspace(0,1,N); 99 tvec = linspace(0,1,N);
108 arcVec = grid.Curve.arcLength(obj.gp,0,tvec); 100 arcVec = obj.arcLength(0,tvec);
109 arcLength = grid.Curve.spline(tvec,arcVec); 101 arcLength = spline(tvec,arcVec);
110 102
111 % Stretch the parameter, construct function with splines 103 % Stretch the parameter, construct function with splines
112 arcParVec = util.fzero_vec(@(s)arcLength(s) - tvec*arcLength(1),[0-10*eps,1+10*eps]); 104 arcParVec = fzero_vec(@(s)arcLength(s) - tvec*arcLength(1),[0-10*eps,1+10*eps]);
113 arcPar = grid.Curve.spline(tvec,arcParVec); 105 arcPar = spline(tvec,arcParVec);
114 106
115 % New function and derivative 107 % New function and derivative
116 g_new = @(t)obj.g(arcPar(t)); 108 g_new = @(t)obj.g(arcPar(t));
117 gp_old = obj.gp; 109 gp_new = @(t) normalize(obj.gp(arcPar(t)));
118 gp_new = @(t) normalize(gp_old(arcPar(t)));
119
120 curve = grid.Curve(g_new,gp_new); 110 curve = grid.Curve(g_new,gp_new);
121 111
122 end 112 end
123 113
124 % how to make it work for methods without returns 114 % how to make it work for methods without returns
257 end 247 end
258 end 248 end
259 249
260 methods (Static) 250 methods (Static)
261 251
262 % Length of arc from parameter t0 to t1 (which may be vectors).
263 % Computed using derivative.
264 function L = arcLength(deriv,t0,t1)
265 speed = @(t) sp(deriv(t));
266
267 function s = sp(deriv)
268 s = sqrt(sum(deriv.^2,1));
269 end
270 L = util.integral_vec(speed,t0,t1);
271 end
272
273 function gp_out = numerical_derivative(g,D1) 252 function gp_out = numerical_derivative(g,D1)
274 m = length(D1); L = 1; % Assume curve parameter from 0 to 1. 253 m = length(D1); L = 1; % Assume curve parameter from 0 to 1.
275 t = linspace(0,L,m); 254 t = linspace(0,L,m);
276 g = g(t)'; 255 gVec = g(t)';
277 gp = (D1*g)'; 256 gpVec = (D1*gVec)';
278 257
279 gp1_fun = grid.Curve.spline(t,gp(1,:)); 258 gp1_fun = spline(t,gpVec(1,:));
280 gp2_fun = grid.Curve.spline(t,gp(2,:)); 259 gp2_fun = spline(t,gpVec(2,:));
281 gp_out = @(t) [gp1_fun(t);gp2_fun(t)]; 260 gp_out = @(t) [gp1_fun(t);gp2_fun(t)];
282 261
283 end
284
285 % Returns a function handle to the spline.
286 function f = spline(tval,fval,spline_order)
287 default_arg('spline_order',4);
288 [m,~] = size(tval);
289 assert(m==1,'Need row vectors.');
290
291 % make vectors longer to be safe slightly beyond edges.
292 dt0 = tval(2)-tval(1); dt1 = tval(end)-tval(end-1);
293 df0 = fval(2)-fval(1); df1 = fval(end)-fval(end-1);
294 tval = [tval(1)-dt0, tval, tval(end)+dt1];
295 fval = [fval(1)-df0, fval, fval(end)+df1];
296
297 f_spline = spapi( optknt(tval,spline_order), tval, fval );
298 f = @(t) fnval(f_spline,t);
299 end 262 end
300 263
301 function obj = line(p1, p2) 264 function obj = line(p1, p2)
302 265
303 function v = g_fun(t) 266 function v = g_fun(t)
378 341
379 end 342 end
380 end 343 end
381 344
382 function g_norm = normalize(g0) 345 function g_norm = normalize(g0)
383 g1 = g0(1,:); g2 = g0(2,:); 346 g1 = g0(1,:);
347 g2 = g0(2,:);
384 normalization = sqrt(sum(g0.^2,1)); 348 normalization = sqrt(sum(g0.^2,1));
385 g_norm = [g1./normalization; g2./normalization]; 349 g_norm = [g1./normalization; g2./normalization];
386 end 350 end
387 351
352 function I = integral_vec(f,a,b)
353 % Wrapper around the built-in function integral that
354 % handles multiple limits.
355
356 Na = length(a);
357 Nb = length(b);
358 assert(Na == 1 || Nb == 1 || Na==Nb,...
359 'a and b must have same length, unless one is a scalar.');
360
361 if(Na>Nb);
362 I = zeros(size(a));
363 for i = 1:Na
364 I(i) = integral(f,a(i),b);
365 end
366 elseif(Nb>Na)
367 I = zeros(size(b));
368 for i = 1:Nb
369 I(i) = integral(f,a,b(i));
370 end
371 else
372 I = zeros(size(b));
373 for i = 1:Nb
374 I(i) = integral(f,a(i),b(i));
375 end
376 end
377 end
378
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.
395 function f = spline(tval,fval,spline_order)
396 default_arg('spline_order',4);
397 [m,~] = size(tval);
398 assert(m==1,'Need row vectors.');
399
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 );
407 f = @(t) fnval(f_spline,t);
408 end
409
388 410
389 411