|
0
|
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
|
