Mercurial > repos > public > sbplib

comparison +parametrization/Curve.m @ 151:3a3cf386bb7e feature/grids

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Moved Curves and Tis from package grid to package parametrization.
author Jonatan Werpers <jonatan@werpers.com>
date Mon, 21 Dec 2015 15:09:32 +0100
parents +grid/Curve.m@06c3034966b7
children 81e0ead29431
comparison
equal deleted inserted replaced
150:d5a794c734bc 151:3a3cf386bb7e
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] = grid.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 p = obj.g(1/2);
72 n = obj.normal(1/2);
73 p = p + n*0.1;
74
75 % Add arrow
76
77 h = text(p(1),p(2),name);
78 h.HorizontalAlignment = 'center';
79 h.VerticalAlignment = 'middle';
80
81 obj.plot();
82 end
83 % Shows curve with name and arrow for direction.
84
85 % Length of arc from parameter t0 to t1 (which may be vectors).
86 % Computed using derivative.
87 function L = arcLength(obj,t0,t1)
88 assert(~isempty(obj.gp),'Curve has no derivative!');
89 speed = @(t) sqrt(sum(obj.gp(t).^2,1));
90 L = integral_vec(speed,t0,t1);
91 end
92
93 % Creates the arc length parameterization of a curve.
94 % N -- number of points used to approximate the arclength function
95 function curve = arcLengthParametrization(obj,N)
96 default_arg('N',100);
97 assert(~isempty(obj.gp),'Curve has no derivative!');
98
99 % Construct arcLength function using splines
100 tvec = linspace(0,1,N);
101 arcVec = obj.arcLength(0,tvec);
102 tFunc = spline(arcVec,tvec); % t as a function of arcLength
103 L = obj.arcLength(0,1);
104 arcPar = @(s) tFunc(s*L);
105
106 % New function and derivative
107 g_new = @(t)obj.g(arcPar(t));
108 gp_new = @(t) normalize(obj.gp(arcPar(t)));
109 curve = grid.Curve(g_new,gp_new);
110
111 end
112
113 % how to make it work for methods without returns
114 function p = subsref(obj,S)
115 %Should i add error checking here?
116 %Maybe if you want performance you fetch obj.g and then use that
117 switch S(1).type
118 case '()'
119 p = obj.g(S.subs{1});
120 % case '.'
121
122 % p = obj.(S.subs);
123 otherwise
124 p = builtin('subsref',obj,S);
125 % error()
126 end
127 end
128
129
130 %% TRANSFORMATION OF A CURVE
131 function D = reverse(C)
132 % g = C.g;
133 % gp = C.gp;
134 % D = grid.Curve(@(t)g(1-t),@(t)-gp(1-t));
135 D = C.transform([],[],-1);
136 end
137
138 function D = transform(C,A,b,flip)
139 default_arg('A',[1 0; 0 1]);
140 default_arg('b',[0; 0]);
141 default_arg('flip',1);
142 if isempty(C.transformation)
143 g = C.g;
144 gp = C.gp;
145 transformation.A = A;
146 transformation.b = b;
147 transformation.flip = flip;
148 else
149 g = C.transformation.base_g;
150 gp = C.transformation.base_gp;
151 A_old = C.transformation.A;
152 b_old = C.transformation.b;
153 flip_old = C.transformation.flip;
154
155 transformation.A = A*A_old;
156 transformation.b = A*b_old + b;
157 transformation.flip = flip*flip_old;
158 end
159
160 D = grid.Curve(g,gp,transformation);
161
162 end
163
164 function D = translate(C,a)
165 g = C.g;
166 gp = C.gp;
167
168 % function v = g_fun(t)
169 % x = g(t);
170 % v(1,:) = x(1,:)+a(1);
171 % v(2,:) = x(2,:)+a(2);
172 % end
173
174 % D = grid.Curve(@g_fun,gp);
175
176 D = C.transform([],a);
177 end
178
179 function D = mirror(C, a, b)
180 assert_size(a,[2,1]);
181 assert_size(b,[2,1]);
182
183 g = C.g;
184 gp = C.gp;
185
186 l = b-a;
187 lx = l(1);
188 ly = l(2);
189
190
191 % fprintf('Singular?\n')
192
193 A = [lx^2-ly^2 2*lx*ly; 2*lx*ly ly^2-lx^2]/(l'*l);
194
195 % function v = g_fun(t)
196 % % v = a + A*(g(t)-a)
197 % x = g(t);
198
199 % ax1 = x(1,:)-a(1);
200 % ax2 = x(2,:)-a(2);
201 % v(1,:) = a(1)+A(1,:)*[ax1;ax2];
202 % v(2,:) = a(2)+A(2,:)*[ax1;ax2];
203 % end
204
205 % function v = gp_fun(t)
206 % v = A*gp(t);
207 % end
208
209 % D = grid.Curve(@g_fun,@gp_fun);
210
211 % g = A(g-a)+a = Ag - Aa + a;
212 b = - A*a + a;
213 D = C.transform(A,b);
214
215 end
216
217 function D = rotate(C,a,rad)
218 assert_size(a, [2,1]);
219 assert_size(rad, [1,1]);
220 g = C.g;
221 gp = C.gp;
222
223
224 A = [cos(rad) -sin(rad); sin(rad) cos(rad)];
225
226 % function v = g_fun(t)
227 % % v = a + A*(g(t)-a)
228 % x = g(t);
229
230 % ax1 = x(1,:)-a(1);
231 % ax2 = x(2,:)-a(2);
232 % v(1,:) = a(1)+A(1,:)*[ax1;ax2];
233 % v(2,:) = a(2)+A(2,:)*[ax1;ax2];
234 % end
235
236 % function v = gp_fun(t)
237 % v = A*gp(t);
238 % end
239
240 % D = grid.Curve(@g_fun,@gp_fun);
241
242
243 % g = A(g-a)+a = Ag - Aa + a;
244 b = - A*a + a;
245 D = C.transform(A,b);
246 end
247 end
248
249 methods (Static)
250
251 % Computes the derivative of g: R -> R^2 using an operator D1
252 function gp_out = numericalDerivative(g,D1)
253 m = length(D1);
254 t = linspace(0,1,m);
255 gVec = g(t)';
256 gpVec = (D1*gVec)';
257
258 gp1_fun = spline(t,gpVec(1,:));
259 gp2_fun = spline(t,gpVec(2,:));
260 gp_out = @(t) [gp1_fun(t);gp2_fun(t)];
261 end
262
263 function obj = line(p1, p2)
264
265 function v = g_fun(t)
266 v(1,:) = p1(1) + t.*(p2(1)-p1(1));
267 v(2,:) = p1(2) + t.*(p2(2)-p1(2));
268 end
269 g = @g_fun;
270
271 obj = grid.Curve(g);
272 end
273
274 function obj = circle(c,r,phi)
275 default_arg('phi',[0; 2*pi])
276 default_arg('c',[0; 0])
277 default_arg('r',1)
278
279 function v = g_fun(t)
280 w = phi(1)+t*(phi(2)-phi(1));
281 v(1,:) = c(1) + r*cos(w);
282 v(2,:) = c(2) + r*sin(w);
283 end
284
285 function v = g_fun_deriv(t)
286 w = phi(1)+t*(phi(2)-phi(1));
287 v(1,:) = -(phi(2)-phi(1))*r*sin(w);
288 v(2,:) = (phi(2)-phi(1))*r*cos(w);
289 end
290
291 obj = grid.Curve(@g_fun,@g_fun_deriv);
292 end
293
294 function obj = bezier(p0, p1, p2, p3)
295 function v = g_fun(t)
296 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);
297 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);
298 end
299
300 function v = g_fun_deriv(t)
301 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));
302 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));
303 end
304
305 obj = grid.Curve(@g_fun,@g_fun_deriv);
306 end
307
308
309 function [g_out,gp_out] = transform_g(g,gp,tr)
310 A = tr.A;
311 b = tr.b;
312 flip = tr.flip;
313
314 function v = g_fun_noflip(t)
315 % v = A*g + b
316 x = g(t);
317
318 v(1,:) = A(1,:)*x+b(1);
319 v(2,:) = A(2,:)*x+b(2);
320 end
321
322 function v = g_fun_flip(t)
323 % v = A*g + b
324 x = g(1-t);
325
326 v(1,:) = A(1,:)*x+b(1);
327 v(2,:) = A(2,:)*x+b(2);
328 end
329
330
331 switch flip
332 case 1
333 g_out = @g_fun_noflip;
334 gp_out = @(t)A*gp(t);
335 case -1
336 g_out = @g_fun_flip;
337 gp_out = @(t)-A*gp(1-t);
338 end
339 end
340
341 end
342 end
343
344
345
346 function g_norm = normalize(g0)
347 g1 = g0(1,:);
348 g2 = g0(2,:);
349 normalization = sqrt(sum(g0.^2,1));
350 g_norm = [g1./normalization; g2./normalization];
351 end
352
353 function I = integral_vec(f,a,b)
354 % Wrapper around the built-in function integral that
355 % handles multiple limits.
356
357 Na = length(a);
358 Nb = length(b);
359 assert(Na == 1 || Nb == 1 || Na==Nb,...
360 'a and b must have same length, unless one is a scalar.');
361
362 if(Na>Nb);
363 I = zeros(size(a));
364 for i = 1:Na
365 I(i) = integral(f,a(i),b);
366 end
367 elseif(Nb>Na)
368 I = zeros(size(b));
369 for i = 1:Nb
370 I(i) = integral(f,a,b(i));
371 end
372 else
373 I = zeros(size(b));
374 for i = 1:Nb
375 I(i) = integral(f,a(i),b(i));
376 end
377 end
378 end
379
380 % Returns a function handle to the spline.
381 function f = spline(tval,fval,spline_order)
382 default_arg('spline_order',4);
383 [m,~] = size(tval);
384 assert(m==1,'Need row vectors.');
385
386 f_spline = spapi( optknt(tval,spline_order), tval, fval );
387 f = @(t) fnval(f_spline,t);
388 end