Mercurial > repos > public > sbplib

comparison +grid/Curve.m @ 107:06c3034966b7 feature/arclen-param

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Added and changed some comments. Also my text editor removed a bunch of whitespace.
author Jonatan Werpers <jonatan@werpers.com>
date Tue, 08 Dec 2015 09:50:39 +0100
parents 1df4f3704b76
children
comparison
equal deleted inserted replaced
106:eb7f592b9512 107:06c3034966b7
1 classdef Curve 1 classdef Curve
2 properties 2 properties
3 g 3 g
4 gp 4 gp
5 transformation 5 transformation
6
7 end 6 end
7
8 methods 8 methods
9 %TODO: 9 %TODO:
10 % Concatenation of curves 10 % Concatenation of curves
11 % Subsections of curves 11 % Subsections of curves
12 % Stretching of curve parameter - done for arc length. 12 % Stretching of curve parameter - done for arc length.
13 % Curve to cell array of linesegments 13 % Curve to cell array of linesegments
14 14
15 % The curve parameter t must be in [0,1]. 15 % Returns a curve object.
16 % Can supply derivative. 16 % g -- curve parametrization for parameter between 0 and 1
17 % gp -- parametrization of curve derivative
17 function obj = Curve(g,gp,transformation) 18 function obj = Curve(g,gp,transformation)
18 default_arg('gp',[]); 19 default_arg('gp',[]);
19 default_arg('transformation',[]); 20 default_arg('transformation',[]);
20 p_test = g(0); 21 p_test = g(0);
21 assert(all(size(p_test) == [2,1]), 'A curve parametrization must return a 2x1 vector.'); 22 assert(all(size(p_test) == [2,1]), 'A curve parametrization must return a 2x1 vector.');
22 23
23 if ~isempty(transformation) 24 if ~isempty(transformation)
24 transformation.base_g = g; 25 transformation.base_g = g;
25 transformation.base_gp = gp; 26 transformation.base_gp = gp;
26 [g,gp] = grid.Curve.transform_g(g,gp,transformation); 27 [g,gp] = grid.Curve.transform_g(g,gp,transformation);
27 end 28 end
28 29
29 obj.g = g; 30 obj.g = g;
30 obj.gp = gp; 31 obj.gp = gp;
31 obj.transformation = transformation; 32 obj.transformation = transformation;
32 33
33 end 34 end
34 35
35 function n = normal(obj,t) 36 function n = normal(obj,t)
36 assert(~isempty(obj.gp),'Curve has no derivative!'); 37 assert(~isempty(obj.gp),'Curve has no derivative!');
37 deriv = obj.gp(t); 38 deriv = obj.gp(t);
86 function L = arcLength(obj,t0,t1) 87 function L = arcLength(obj,t0,t1)
87 assert(~isempty(obj.gp),'Curve has no derivative!'); 88 assert(~isempty(obj.gp),'Curve has no derivative!');
88 speed = @(t) sqrt(sum(obj.gp(t).^2,1)); 89 speed = @(t) sqrt(sum(obj.gp(t).^2,1));
89 L = integral_vec(speed,t0,t1); 90 L = integral_vec(speed,t0,t1);
90 end 91 end
91 92
92 % Arc length parameterization. 93 % Creates the arc length parameterization of a curve.
94 % N -- number of points used to approximate the arclength function
93 function curve = arcLengthParametrization(obj,N) 95 function curve = arcLengthParametrization(obj,N)
94 default_arg('N',100); 96 default_arg('N',100);
95 assert(~isempty(obj.gp),'Curve has no derivative!'); 97 assert(~isempty(obj.gp),'Curve has no derivative!');
96 98
97 % Construct arcLength function using splines 99 % Construct arcLength function using splines
98 tvec = linspace(0,1,N); 100 tvec = linspace(0,1,N);
99 arcVec = obj.arcLength(0,tvec); 101 arcVec = obj.arcLength(0,tvec);
100 tFunc = spline(arcVec,tvec); % t as a function of arcLength 102 tFunc = spline(arcVec,tvec); % t as a function of arcLength
101 L = obj.arcLength(0,1); 103 L = obj.arcLength(0,1);
102 arcPar = @(s) tFunc(s*L); 104 arcPar = @(s) tFunc(s*L);
103 105
104 % New function and derivative 106 % New function and derivative
105 g_new = @(t)obj.g(arcPar(t)); 107 g_new = @(t)obj.g(arcPar(t));
106 gp_new = @(t) normalize(obj.gp(arcPar(t))); 108 gp_new = @(t) normalize(obj.gp(arcPar(t)));
107 curve = grid.Curve(g_new,gp_new); 109 curve = grid.Curve(g_new,gp_new);
108 110
109 end 111 end
110 112
111 % how to make it work for methods without returns 113 % how to make it work for methods without returns
112 function p = subsref(obj,S) 114 function p = subsref(obj,S)
113 %Should i add error checking here? 115 %Should i add error checking here?
114 %Maybe if you want performance you fetch obj.g and then use that 116 %Maybe if you want performance you fetch obj.g and then use that
115 switch S(1).type 117 switch S(1).type
243 D = C.transform(A,b); 245 D = C.transform(A,b);
244 end 246 end
245 end 247 end
246 248
247 methods (Static) 249 methods (Static)
248 250
249 251 % Computes the derivative of g: R -> R^2 using an operator D1
250 function gp_out = numericalDerivative(g,D1) 252 function gp_out = numericalDerivative(g,D1)
251 m = length(D1); L = 1; % Assume curve parameter from 0 to 1. 253 m = length(D1);
252 t = linspace(0,L,m); 254 t = linspace(0,1,m);
253 gVec = g(t)'; 255 gVec = g(t)';
254 gpVec = (D1*gVec)'; 256 gpVec = (D1*gVec)';
255 257
256 gp1_fun = spline(t,gpVec(1,:)); 258 gp1_fun = spline(t,gpVec(1,:));
257 gp2_fun = spline(t,gpVec(2,:)); 259 gp2_fun = spline(t,gpVec(2,:));
258 gp_out = @(t) [gp1_fun(t);gp2_fun(t)]; 260 gp_out = @(t) [gp1_fun(t);gp2_fun(t)];
259 261 end
260 end 262
261
262 function obj = line(p1, p2) 263 function obj = line(p1, p2)
263 264
264 function v = g_fun(t) 265 function v = g_fun(t)
265 v(1,:) = p1(1) + t.*(p2(1)-p1(1)); 266 v(1,:) = p1(1) + t.*(p2(1)-p1(1));
266 v(2,:) = p1(2) + t.*(p2(2)-p1(2)); 267 v(2,:) = p1(2) + t.*(p2(2)-p1(2));
335 g_out = @g_fun_flip; 336 g_out = @g_fun_flip;
336 gp_out = @(t)-A*gp(1-t); 337 gp_out = @(t)-A*gp(1-t);
337 end 338 end
338 end 339 end
339 340
340 end 341 end
341 end 342 end
342 343
344
345
343 function g_norm = normalize(g0) 346 function g_norm = normalize(g0)
344 g1 = g0(1,:); 347 g1 = g0(1,:);
345 g2 = g0(2,:); 348 g2 = g0(2,:);
346 normalization = sqrt(sum(g0.^2,1)); 349 normalization = sqrt(sum(g0.^2,1));
347 g_norm = [g1./normalization; g2./normalization]; 350 g_norm = [g1./normalization; g2./normalization];
348 end 351 end
349 352
353 356
354 Na = length(a); 357 Na = length(a);
355 Nb = length(b); 358 Nb = length(b);
356 assert(Na == 1 || Nb == 1 || Na==Nb,... 359 assert(Na == 1 || Nb == 1 || Na==Nb,...
357 'a and b must have same length, unless one is a scalar.'); 360 'a and b must have same length, unless one is a scalar.');
358 361
359 if(Na>Nb); 362 if(Na>Nb);
360 I = zeros(size(a)); 363 I = zeros(size(a));
361 for i = 1:Na 364 for i = 1:Na
362 I(i) = integral(f,a(i),b); 365 I(i) = integral(f,a(i),b);
363 end 366 end
381 assert(m==1,'Need row vectors.'); 384 assert(m==1,'Need row vectors.');
382 385
383 f_spline = spapi( optknt(tval,spline_order), tval, fval ); 386 f_spline = spapi( optknt(tval,spline_order), tval, fval );
384 f = @(t) fnval(f_spline,t); 387 f = @(t) fnval(f_spline,t);
385 end 388 end
386
387
388