Mercurial > repos > public > sbplib

comparison +grid/Curve.m @ 109:142974097efc

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merged in feature/arclen-param (pull request #1) Feature/arclen param
author Martin Almquist <martin.almquist@it.uu.se>
date Tue, 08 Dec 2015 10:42:20 +0100
parents 06c3034966b7
children
comparison
equal deleted inserted replaced
101:9933169d2651 109:142974097efc
2 properties 2 properties
3 g 3 g
4 gp 4 gp
5 transformation 5 transformation
6 end 6 end
7
7 methods 8 methods
8 %TODO: 9 %TODO:
9 % Errors or FD if there is no derivative function added.
10
11 % Concatenation of curves 10 % Concatenation of curves
12 % Subsections of curves 11 % Subsections of curves
13 % Stretching of curve paramter 12 % Stretching of curve parameter - done for arc length.
14 % Curve to cell array of linesegments 13 % Curve to cell array of linesegments
15 14
16 % Should accept derivative of curve. 15 % Returns a curve object.
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.');
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 end 33
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 34 end
77 35
78 function n = normal(obj,t) 36 function n = normal(obj,t)
37 assert(~isempty(obj.gp),'Curve has no derivative!');
79 deriv = obj.gp(t); 38 deriv = obj.gp(t);
80 normalization = sqrt(sum(deriv.^2,1)); 39 normalization = sqrt(sum(deriv.^2,1));
81 n = [-deriv(2,:)./normalization; deriv(1,:)./normalization]; 40 n = [-deriv(2,:)./normalization; deriv(1,:)./normalization];
82 end 41 end
83 42
84 43
85 % Plots a curve g(t) for 0<t<1, using n points. Retruns a handle h to the plotted curve. 44 % Plots a curve g(t) for 0<t<1, using n points. Returns a handle h to the plotted curve.
86 % h = plot_curve(g,n) 45 % h = plot_curve(g,n)
87 function h = plot(obj,n) 46 function h = plot(obj,n)
88 default_arg('n',100); 47 default_arg('n',100);
89 48
90 t = linspace(0,1,n); 49 t = linspace(0,1,n);
91
92 p = obj.g(t); 50 p = obj.g(t);
93
94 h = line(p(1,:),p(2,:)); 51 h = line(p(1,:),p(2,:));
95 end 52 end
96 53
97 function h= plot_normals(obj,l,n,m) 54 function h= plot_normals(obj,l,n,m)
98 default_arg('l',0.1); 55 default_arg('l',0.1);
123 80
124 obj.plot(); 81 obj.plot();
125 end 82 end
126 % Shows curve with name and arrow for direction. 83 % Shows curve with name and arrow for direction.
127 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
128 112
129 % how to make it work for methods without returns 113 % how to make it work for methods without returns
130 function p = subsref(obj,S) 114 function p = subsref(obj,S)
131 %Should i add error checking here? 115 %Should i add error checking here?
132 %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
139 otherwise 123 otherwise
140 p = builtin('subsref',obj,S); 124 p = builtin('subsref',obj,S);
141 % error() 125 % error()
142 end 126 end
143 end 127 end
144
145
146 128
147 129
148 %% TRANSFORMATION OF A CURVE 130 %% TRANSFORMATION OF A CURVE
149 function D = reverse(C) 131 function D = reverse(C)
150 % g = C.g; 132 % g = C.g;
263 D = C.transform(A,b); 245 D = C.transform(A,b);
264 end 246 end
265 end 247 end
266 248
267 methods (Static) 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
268 function obj = line(p1, p2) 263 function obj = line(p1, p2)
269 264
270 function v = g_fun(t) 265 function v = g_fun(t)
271 v(1,:) = p1(1) + t.*(p2(1)-p1(1)); 266 v(1,:) = p1(1) + t.*(p2(1)-p1(1));
272 v(2,:) = p1(2) + t.*(p2(2)-p1(2)); 267 v(2,:) = p1(2) + t.*(p2(2)-p1(2));
343 end 338 end
344 end 339 end
345 340
346 end 341 end
347 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