Mercurial > repos > public > sbplib

comparison +grid/Curve.m @ 0:48b6fb693025

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Initial commit.
author Jonatan Werpers <jonatan@werpers.com>
date Thu, 17 Sep 2015 10:12:50 +0200
parents
children 3c39dd714fb6
comparison
equal deleted inserted replaced
-1:000000000000 0:48b6fb693025
1 classdef Curve
2 properties
3 g
4 gp
5 transformation
6 end
7 methods
8 %TODO:
9 % Errors or FD if there is no derivative function added.
10
11 % Concatenation of curves
12 % Subsections of curves
13 % Stretching of curve paramter
14 % Curve to cell array of linesegments
15
16 % Should accept derivative of curve.
17 function obj = Curve(g,gp,transformation)
18 default_arg('gp',[]);
19 default_arg('transformation',[]);
20 p_test = g(0);
21 assert(all(size(p_test) == [2,1]), 'A curve parametrization must return a 2x1 vector.');
22
23 if ~isempty(transformation)
24 transformation.base_g = g;
25 transformation.base_gp = gp;
26 [g,gp] = grid.Curve.transform_g(g,gp,transformation);
27 end
28
29 obj.g = g;
30 obj.gp = gp;
31 obj.transformation = transformation;
32 end
33
34 % Made up length calculation!! Science this before actual use!!
35 % Calculates the length of the curve. Makes sure the longet segment used is shorter than h_max.
36 function [L,t] = curve_length(C,h_max)
37 default_arg('h_max',0.001);
38 g = C.g;
39 h = 0.1;
40 m = 1/h+1;
41 t = linspace(0,1,m);
42
43 [p,d] = get_d(t,g);
44
45 while any(d>h_max)
46 I = find(d>h_max);
47
48 % plot(p(1,:),p(2,:),'.')
49 % waitforbuttonpress
50
51 new_t = [];
52 for i = I
53 new_t(end +1) = (t(i)+t(i+1))/2;
54 end
55 t = [t new_t];
56 t = sort(t);
57
58 [p,d] = get_d(t,g);
59 end
60
61 L = sum(d);
62
63 function [p,d] = get_d(t,g)
64 n = length(t);
65
66 p = zeros(2,n);
67 for i = 1:n
68 p(:,i) = g(t(i));
69 end
70
71 d = zeros(1,n-1);
72 for i = 1:n-1
73 d(i) = norm(p(:,i) - p(:,i+1));
74 end
75 end
76 end
77
78 function n = normal(obj,t)
79 deriv = obj.gp(t);
80 normalization = sqrt(sum(deriv.^2,1));
81 n = [-deriv(2,:)./normalization; deriv(1,:)./normalization];
82 end
83
84
85 % Plots a curve g(t) for 0<t<1, using n points. Retruns a handle h to the plotted curve.
86 % h = plot_curve(g,n)
87 function h = plot(obj,n)
88 default_arg('n',100);
89
90 t = linspace(0,1,n);
91
92 p = obj.g(t);
93
94 h = line(p(1,:),p(2,:));
95 end
96
97 function h= plot_normals(obj,l,n,m)
98 default_arg('l',0.1);
99 default_arg('n',10);
100 default_arg('m',100);
101 t_n = linspace(0,1,n);
102
103 normals = obj.normal(t_n)*l;
104
105 n0 = obj.g(t_n);
106 n1 = n0 + normals;
107
108 h = line([n0(1,:); n1(1,:)],[n0(2,:); n1(2,:)]);
109 set(h,'Color',Color.red);
110 obj.plot(m);
111 end
112
113 function h= show(obj,name)
114 p = obj.g(1/2);
115 n = obj.normal(1/2);
116 p = p + n*0.1;
117
118 % Add arrow
119
120 h = text(p(1),p(2),name);
121 h.HorizontalAlignment = 'center';
122 h.VerticalAlignment = 'middle';
123
124 obj.plot();
125 end
126 % Shows curve with name and arrow for direction.
127
128
129 % how to make it work for methods without returns
130 function p = subsref(obj,S)
131 %Should i add error checking here?
132 %Maybe if you want performance you fetch obj.g and then use that
133 switch S(1).type
134 case '()'
135 p = obj.g(S.subs{1});
136 % case '.'
137
138 % p = obj.(S.subs);
139 otherwise
140 p = builtin('subsref',obj,S);
141 % error()
142 end
143 end
144
145
146
147
148 %% TRANSFORMATION OF A CURVE
149 function D = reverse(C)
150 % g = C.g;
151 % gp = C.gp;
152 % D = grid.Curve(@(t)g(1-t),@(t)-gp(1-t));
153 D = C.transform([],[],-1);
154 end
155
156 function D = transform(C,A,b,flip)
157 default_arg('A',[1 0; 0 1]);
158 default_arg('b',[0; 0]);
159 default_arg('flip',1);
160 if isempty(C.transformation)
161 g = C.g;
162 gp = C.gp;
163 transformation.A = A;
164 transformation.b = b;
165 transformation.flip = flip;
166 else
167 g = C.transformation.base_g;
168 gp = C.transformation.base_gp;
169 A_old = C.transformation.A;
170 b_old = C.transformation.b;
171 flip_old = C.transformation.flip;
172
173 transformation.A = A*A_old;
174 transformation.b = A*b_old + b;
175 transformation.flip = flip*flip_old;
176 end
177
178 D = grid.Curve(g,gp,transformation);
179
180 end
181
182 function D = translate(C,a)
183 g = C.g;
184 gp = C.gp;
185
186 % function v = g_fun(t)
187 % x = g(t);
188 % v(1,:) = x(1,:)+a(1);
189 % v(2,:) = x(2,:)+a(2);
190 % end
191
192 % D = grid.Curve(@g_fun,gp);
193
194 D = C.transform([],a);
195 end
196
197 function D = mirror(C, a, b)
198 assert_size(a,[2,1]);
199 assert_size(b,[2,1]);
200
201 g = C.g;
202 gp = C.gp;
203
204 l = b-a;
205 lx = l(1);
206 ly = l(2);
207
208
209 % fprintf('Singular?\n')
210
211 A = [lx^2-ly^2 2*lx*ly; 2*lx*ly ly^2-lx^2]/(l'*l);
212
213 % function v = g_fun(t)
214 % % v = a + A*(g(t)-a)
215 % x = g(t);
216
217 % ax1 = x(1,:)-a(1);
218 % ax2 = x(2,:)-a(2);
219 % v(1,:) = a(1)+A(1,:)*[ax1;ax2];
220 % v(2,:) = a(2)+A(2,:)*[ax1;ax2];
221 % end
222
223 % function v = gp_fun(t)
224 % v = A*gp(t);
225 % end
226
227 % D = grid.Curve(@g_fun,@gp_fun);
228
229 % g = A(g-a)+a = Ag - Aa + a;
230 b = - A*a + a;
231 D = C.transform(A,b);
232
233 end
234
235 function D = rotate(C,a,rad)
236 assert_size(a, [2,1]);
237 assert_size(rad, [1,1]);
238 g = C.g;
239 gp = C.gp;
240
241
242 A = [cos(rad) -sin(rad); sin(rad) cos(rad)];
243
244 % function v = g_fun(t)
245 % % v = a + A*(g(t)-a)
246 % x = g(t);
247
248 % ax1 = x(1,:)-a(1);
249 % ax2 = x(2,:)-a(2);
250 % v(1,:) = a(1)+A(1,:)*[ax1;ax2];
251 % v(2,:) = a(2)+A(2,:)*[ax1;ax2];
252 % end
253
254 % function v = gp_fun(t)
255 % v = A*gp(t);
256 % end
257
258 % D = grid.Curve(@g_fun,@gp_fun);
259
260
261 % g = A(g-a)+a = Ag - Aa + a;
262 b = - A*a + a;
263 D = C.transform(A,b);
264 end
265 end
266
267 methods (Static)
268 function obj = line(p1, p2)
269
270 function v = g_fun(t)
271 v(1,:) = p1(1) + t.*(p2(1)-p1(1));
272 v(2,:) = p1(2) + t.*(p2(2)-p1(2));
273 end
274 g = @g_fun;
275
276 obj = grid.Curve(g);
277 end
278
279 function obj = circle(c,r,phi)
280 default_arg('phi',[0; 2*pi])
281 default_arg('c',[0; 0])
282 default_arg('r',1)
283
284 function v = g_fun(t)
285 w = phi(1)+t*(phi(2)-phi(1));
286 v(1,:) = c(1) + r*cos(w);
287 v(2,:) = c(2) + r*sin(w);
288 end
289
290 function v = g_fun_deriv(t)
291 w = phi(1)+t*(phi(2)-phi(1));
292 v(1,:) = -(phi(2)-phi(1))*r*sin(w);
293 v(2,:) = (phi(2)-phi(1))*r*cos(w);
294 end
295
296 obj = grid.Curve(@g_fun,@g_fun_deriv);
297 end
298
299 function obj = bezier(p0, p1, p2, p3)
300 function v = g_fun(t)
301 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);
302 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);
303 end
304
305 function v = g_fun_deriv(t)
306 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));
307 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));
308 end
309
310 obj = grid.Curve(@g_fun,@g_fun_deriv);
311 end
312
313
314 function [g_out,gp_out] = transform_g(g,gp,tr)
315 A = tr.A;
316 b = tr.b;
317 flip = tr.flip;
318
319 function v = g_fun_noflip(t)
320 % v = A*g + b
321 x = g(t);
322
323 v(1,:) = A(1,:)*x+b(1);
324 v(2,:) = A(2,:)*x+b(2);
325 end
326
327 function v = g_fun_flip(t)
328 % v = A*g + b
329 x = g(1-t);
330
331 v(1,:) = A(1,:)*x+b(1);
332 v(2,:) = A(2,:)*x+b(2);
333 end
334
335
336 switch flip
337 case 1
338 g_out = @g_fun_noflip;
339 gp_out = @(t)A*gp(t);
340 case -1
341 g_out = @g_fun_flip;
342 gp_out = @(t)-A*gp(1-t);
343 end
344 end
345
346 end
347 end