Mercurial > repos > public > sbplib

comparison +parametrization/Curve.m @ 886:8894e9c49e40 feature/timesteppers

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge with default for latest changes
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Thu, 15 Nov 2018 16:36:21 -0800
parents 0776fa4754ff
children acf19ecd1338 465fc81e12e8
comparison
equal deleted inserted replaced
816:b5e5b195da1e 886:8894e9c49e40
1 classdef Curve
2 properties
3 g
4 gp
5 transformation
6 end
7
8 methods
9 %TODO:
10 % Concatenation of curves
11 % Subsections of curves
12 % Stretching of curve parameter - done for arc length.
13 % Curve to cell array of linesegments
14
15 % Returns a curve object.
16 % g -- curve parametrization for parameter between 0 and 1
17 % gp -- parametrization of curve derivative
18 function obj = Curve(g,gp,transformation)
19 default_arg('gp',[]);
20 default_arg('transformation',[]);
21 p_test = g(0);
22 assert(all(size(p_test) == [2,1]), 'A curve parametrization must return a 2x1 vector.');
23
24 if ~isempty(transformation)
25 transformation.base_g = g;
26 transformation.base_gp = gp;
27 [g,gp] = parametrization.Curve.transform_g(g,gp,transformation);
28 end
29
30 obj.g = g;
31 obj.gp = gp;
32 obj.transformation = transformation;
33
34 end
35
36 function n = normal(obj,t)
37 assert(~isempty(obj.gp),'Curve has no derivative!');
38 deriv = obj.gp(t);
39 normalization = sqrt(sum(deriv.^2,1));
40 n = [-deriv(2,:)./normalization; deriv(1,:)./normalization];
41 end
42
43
44 % Plots a curve g(t) for 0<t<1, using n points. Returns a handle h to the plotted curve.
45 % h = plot_curve(g,n)
46 function h = plot(obj,n)
47 default_arg('n',100);
48
49 t = linspace(0,1,n);
50 p = obj.g(t);
51 h = line(p(1,:),p(2,:));
52 end
53
54 function h= plot_normals(obj,l,n,m)
55 default_arg('l',0.1);
56 default_arg('n',10);
57 default_arg('m',100);
58 t_n = linspace(0,1,n);
59
60 normals = obj.normal(t_n)*l;
61
62 n0 = obj.g(t_n);
63 n1 = n0 + normals;
64
65 h = line([n0(1,:); n1(1,:)],[n0(2,:); n1(2,:)]);
66 set(h,'Color',Color.red);
67 obj.plot(m);
68 end
69
70 function h= show(obj,name)
71 default_arg('name', '')
72 p = obj.g(1/2);
73 n = obj.normal(1/2);
74 p = p + n*0.1;
75
76 % Add arrow
77
78
79 if ~isempty(name)
80 h = text(p(1),p(2),name);
81 h.HorizontalAlignment = 'center';
82 h.VerticalAlignment = 'middle';
83 end
84
85 h = obj.plot();
86 end
87 % Shows curve with name and arrow for direction.
88
89 % Length of arc from parameter t0 to t1 (which may be vectors).
90 % Computed using derivative.
91 function L = arcLength(obj,t0,t1)
92 assert(~isempty(obj.gp),'Curve has no derivative!');
93 speed = @(t) sqrt(sum(obj.gp(t).^2,1));
94 L = integral_vec(speed,t0,t1);
95 end
96
97 % Creates the arc length parameterization of a curve.
98 % N -- number of points used to approximate the arclength function
99 function curve = arcLengthParametrization(obj,N)
100 default_arg('N',100);
101 assert(~isempty(obj.gp),'Curve has no derivative!');
102
103 % Construct arcLength function using splines
104 tvec = linspace(0,1,N);
105 arcVec = obj.arcLength(0,tvec);
106 tFunc = spline(arcVec,tvec); % t as a function of arcLength
107 L = obj.arcLength(0,1);
108 arcPar = @(s) tFunc(s*L);
109
110 % New function and derivative
111 g_new = @(t)obj.g(arcPar(t));
112 gp_new = @(t) normalize(obj.gp(arcPar(t)));
113 curve = parametrization.Curve(g_new,gp_new);
114
115 end
116
117 % how to make it work for methods without returns
118 function p = subsref(obj,S)
119 %Should i add error checking here?
120 %Maybe if you want performance you fetch obj.g and then use that
121 switch S(1).type
122 case '()'
123 p = obj.g(S.subs{1});
124 % case '.'
125
126 % p = obj.(S.subs);
127 otherwise
128 p = builtin('subsref',obj,S);
129 % error()
130 end
131 end
132
133
134 %% TRANSFORMATION OF A CURVE
135 function D = reverse(C)
136 % g = C.g;
137 % gp = C.gp;
138 % D = parametrization.Curve(@(t)g(1-t),@(t)-gp(1-t));
139 D = C.transform([],[],-1);
140 end
141
142 function D = transform(C,A,b,flip)
143 default_arg('A',[1 0; 0 1]);
144 default_arg('b',[0; 0]);
145 default_arg('flip',1);
146 if isempty(C.transformation)
147 g = C.g;
148 gp = C.gp;
149 transformation.A = A;
150 transformation.b = b;
151 transformation.flip = flip;
152 else
153 g = C.transformation.base_g;
154 gp = C.transformation.base_gp;
155 A_old = C.transformation.A;
156 b_old = C.transformation.b;
157 flip_old = C.transformation.flip;
158
159 transformation.A = A*A_old;
160 transformation.b = A*b_old + b;
161 transformation.flip = flip*flip_old;
162 end
163
164 D = parametrization.Curve(g,gp,transformation);
165
166 end
167
168 function D = translate(C,a)
169 g = C.g;
170 gp = C.gp;
171
172 % function v = g_fun(t)
173 % x = g(t);
174 % v(1,:) = x(1,:)+a(1);
175 % v(2,:) = x(2,:)+a(2);
176 % end
177
178 % D = parametrization.Curve(@g_fun,gp);
179
180 D = C.transform([],a);
181 end
182
183 function D = mirror(C, a, b)
184 assertSize(a,[2,1]);
185 assertSize(b,[2,1]);
186
187 g = C.g;
188 gp = C.gp;
189
190 l = b-a;
191 lx = l(1);
192 ly = l(2);
193
194
195 % fprintf('Singular?\n')
196
197 A = [lx^2-ly^2 2*lx*ly; 2*lx*ly ly^2-lx^2]/(l'*l);
198
199 % function v = g_fun(t)
200 % % v = a + A*(g(t)-a)
201 % x = g(t);
202
203 % ax1 = x(1,:)-a(1);
204 % ax2 = x(2,:)-a(2);
205 % v(1,:) = a(1)+A(1,:)*[ax1;ax2];
206 % v(2,:) = a(2)+A(2,:)*[ax1;ax2];
207 % end
208
209 % function v = gp_fun(t)
210 % v = A*gp(t);
211 % end
212
213 % D = parametrization.Curve(@g_fun,@gp_fun);
214
215 % g = A(g-a)+a = Ag - Aa + a;
216 b = - A*a + a;
217 D = C.transform(A,b);
218
219 end
220
221 function D = rotate(C,a,rad)
222 assertSize(a, [2,1]);
223 assertSize(rad, [1,1]);
224 g = C.g;
225 gp = C.gp;
226
227
228 A = [cos(rad) -sin(rad); sin(rad) cos(rad)];
229
230 % function v = g_fun(t)
231 % % v = a + A*(g(t)-a)
232 % x = g(t);
233
234 % ax1 = x(1,:)-a(1);
235 % ax2 = x(2,:)-a(2);
236 % v(1,:) = a(1)+A(1,:)*[ax1;ax2];
237 % v(2,:) = a(2)+A(2,:)*[ax1;ax2];
238 % end
239
240 % function v = gp_fun(t)
241 % v = A*gp(t);
242 % end
243
244 % D = parametrization.Curve(@g_fun,@gp_fun);
245
246
247 % g = A(g-a)+a = Ag - Aa + a;
248 b = - A*a + a;
249 D = C.transform(A,b);
250 end
251 end
252
253 methods (Static)
254
255 % Computes the derivative of g: R -> R^2 using an operator D1
256 function gp_out = numericalDerivative(g,D1)
257 m = length(D1);
258 t = linspace(0,1,m);
259 gVec = g(t)';
260 gpVec = (D1*gVec)';
261
262 gp1_fun = spline(t,gpVec(1,:));
263 gp2_fun = spline(t,gpVec(2,:));
264 gp_out = @(t) [gp1_fun(t);gp2_fun(t)];
265 end
266
267 function obj = line(p1, p2)
268
269 function v = g_fun(t)
270 v(1,:) = p1(1) + t.*(p2(1)-p1(1));
271 v(2,:) = p1(2) + t.*(p2(2)-p1(2));
272 end
273
274 function v = g_fun_deriv(t)
275 v(1,:) = t.*0 + (p2(1)-p1(1));
276 v(2,:) = t.*0 + (p2(2)-p1(2));
277 end
278
279 obj = parametrization.Curve(@g_fun, @g_fun_deriv);
280 end
281
282 function obj = circle(c,r,phi)
283 default_arg('phi',[0; 2*pi])
284 default_arg('c',[0; 0])
285 default_arg('r',1)
286
287 function v = g_fun(t)
288 w = phi(1)+t*(phi(2)-phi(1));
289 v(1,:) = c(1) + r*cos(w);
290 v(2,:) = c(2) + r*sin(w);
291 end
292
293 function v = g_fun_deriv(t)
294 w = phi(1)+t*(phi(2)-phi(1));
295 v(1,:) = -(phi(2)-phi(1))*r*sin(w);
296 v(2,:) = (phi(2)-phi(1))*r*cos(w);
297 end
298
299 obj = parametrization.Curve(@g_fun,@g_fun_deriv);
300 end
301
302 function obj = bezier(p0, p1, p2, p3)
303 function v = g_fun(t)
304 v(1,:) = (1-t).^3*p0(1) + 3*(1-t).^2.*t*p1(1) + 3*(1-t).*t.^2*p2(1) + t.^3*p3(1);
305 v(2,:) = (1-t).^3*p0(2) + 3*(1-t).^2.*t*p1(2) + 3*(1-t).*t.^2*p2(2) + t.^3*p3(2);
306 end
307
308 function v = g_fun_deriv(t)
309 v(1,:) = 3*(1-t).^2*(p1(1)-p0(1)) + 6*(1-t).*t*(p2(1)-p1(1)) + 3*t.^2*(p3(1)-p2(1));
310 v(2,:) = 3*(1-t).^2*(p1(2)-p0(2)) + 6*(1-t).*t*(p2(2)-p1(2)) + 3*t.^2*(p3(2)-p2(2));
311 end
312
313 obj = parametrization.Curve(@g_fun,@g_fun_deriv);
314 end
315
316
317 function [g_out,gp_out] = transform_g(g,gp,tr)
318 A = tr.A;
319 b = tr.b;
320 flip = tr.flip;
321
322 function v = g_fun_noflip(t)
323 % v = A*g + b
324 x = g(t);
325
326 v(1,:) = A(1,:)*x+b(1);
327 v(2,:) = A(2,:)*x+b(2);
328 end
329
330 function v = g_fun_flip(t)
331 % v = A*g + b
332 x = g(1-t);
333
334 v(1,:) = A(1,:)*x+b(1);
335 v(2,:) = A(2,:)*x+b(2);
336 end
337
338
339 switch flip
340 case 1
341 g_out = @g_fun_noflip;
342 gp_out = @(t)A*gp(t);
343 case -1
344 g_out = @g_fun_flip;
345 gp_out = @(t)-A*gp(1-t);
346 end
347 end
348
349 end
350 end
351
352
353
354 function g_norm = normalize(g0)
355 g1 = g0(1,:);
356 g2 = g0(2,:);
357 normalization = sqrt(sum(g0.^2,1));
358 g_norm = [g1./normalization; g2./normalization];
359 end
360
361 function I = integral_vec(f,a,b)
362 % Wrapper around the built-in function integral that
363 % handles multiple limits.
364
365 Na = length(a);
366 Nb = length(b);
367 assert(Na == 1 || Nb == 1 || Na==Nb,...
368 'a and b must have same length, unless one is a scalar.');
369
370 if(Na>Nb);
371 I = zeros(size(a));
372 for i = 1:Na
373 I(i) = integral(f,a(i),b);
374 end
375 elseif(Nb>Na)
376 I = zeros(size(b));
377 for i = 1:Nb
378 I(i) = integral(f,a,b(i));
379 end
380 else
381 I = zeros(size(b));
382 for i = 1:Nb
383 I(i) = integral(f,a(i),b(i));
384 end
385 end
386 end
387
388 % Returns a function handle to the spline.
389 function f = spline(tval,fval,spline_order)
390 default_arg('spline_order',4);
391 [m,~] = size(tval);
392 assert(m==1,'Need row vectors.');
393
394 f_spline = spapi( optknt(tval,spline_order), tval, fval );
395 f = @(t) fnval(f_spline,t);
396 end