Mercurial > repos > public > sbplib

comparison +grid/Curve.m @ 98:b30f3d8845f4 feature/arclen-param

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Removed arclength property. Constructor compatible again. Translation and rotation are fast now, but ti is slow.
author Martin Almquist <martin.almquist@it.uu.se>
date Sun, 06 Dec 2015 13:47:12 +0100
parents 0a29a60e0b21
children ecf77a50d4fe
comparison
equal deleted inserted replaced
97:33dba20b4b9d 98:b30f3d8845f4
1 classdef Curve 1 classdef Curve
2 properties 2 properties
3 g 3 g
4 gp 4 gp
5 transformation 5 transformation
6 arc_length_fun 6
7
8 % arc_length_fun:
9 % function handle. arc_length(s)=distance
10 % between t=0 and t=s.
11 end 7 end
12 methods 8 methods
13 %TODO: 9 %TODO:
14 % Errors or FD if there is no derivative function added. 10 % Errors or FD if there is no derivative function added.
15 % -semi-done 11 % -semi-done
16 12
17 % Concatenation of curves 13 % Concatenation of curves
18 % Subsections of curves 14 % Subsections of curves
19 % Stretching of curve paramter - semi-done 15 % Stretching of curve parameter - semi-done
20 % Curve to cell array of linesegments 16 % Curve to cell array of linesegments
21 17
22 % Should supply either derivative or a difference operator D1 18 % Can supply either derivative or a difference operator D1
23 function obj = Curve(g,gp,transformation,D1,arc_length_fun) 19 function obj = Curve(g,gp,transformation,D1)
24 default_arg('gp',[]); 20 default_arg('gp',[]);
25 default_arg('transformation',[]); 21 default_arg('transformation',[]);
26 default_arg('D1',[]); 22 default_arg('D1',[]);
27 default_arg('arc_length_fun',[]);
28 p_test = g(0); 23 p_test = g(0);
29 assert(all(size(p_test) == [2,1]), 'A curve parametrization must return a 2x1 vector.'); 24 assert(all(size(p_test) == [2,1]), 'A curve parametrization must return a 2x1 vector.');
30 25
31 if(isempty(gp) && isempty(D1)); 26 if(isempty(gp) && ~isempty(D1))
32 % Should be error instead of disp.
33 disp(['You really should supply either the exact derivative ',...
34 'or a suitable difference operator to compute an approximation.']);
35 elseif(isempty(gp))
36 gp = grid.Curve.numerical_derivative(g,D1); 27 gp = grid.Curve.numerical_derivative(g,D1);
37 end 28 end
38 29
39 if ~isempty(transformation) 30 if ~isempty(transformation)
40 transformation.base_g = g; 31 transformation.base_g = g;
41 transformation.base_gp = gp; 32 transformation.base_gp = gp;
42 [g,gp] = grid.Curve.transform_g(g,gp,transformation); 33 [g,gp] = grid.Curve.transform_g(g,gp,transformation);
43 end 34 end
44
45 if(isempty(arc_length_fun))
46 if(~isempty(D1)); N = length(D1); % Same accuracy as for deriv.
47 else N = 101; % Best way to let user choose?
48 end
49 tvec = linspace(0,1,N);
50 arc_vec = grid.Curve.arc_length(gp,0,tvec);
51 arc_length_fun = grid.Curve.spline(tvec,arc_vec);
52 end
53 35
54 obj.g = g; 36 obj.g = g;
55 obj.gp = gp; 37 obj.gp = gp;
56 obj.transformation = transformation; 38 obj.transformation = transformation;
57 obj.arc_length_fun = arc_length_fun; 39
58
59
60 end
61
62 % Made up length calculation!! Science this before actual use!!
63 % Calculates the length of the curve. Makes sure the longet segment used is shorter than h_max.
64 function [L,t] = curve_length(C,h_max)
65 default_arg('h_max',0.001);
66 g = C.g;
67 h = 0.1;
68 m = 1/h+1;
69 t = linspace(0,1,m);
70
71 [p,d] = get_d(t,g);
72
73 while any(d>h_max)
74 I = find(d>h_max);
75
76 % plot(p(1,:),p(2,:),'.')
77 % waitforbuttonpress
78
79 new_t = [];
80 for i = I
81 new_t(end +1) = (t(i)+t(i+1))/2;
82 end
83 t = [t new_t];
84 t = sort(t);
85
86 [p,d] = get_d(t,g);
87 end
88
89 L = sum(d);
90
91 function [p,d] = get_d(t,g)
92 n = length(t);
93
94 p = zeros(2,n);
95 for i = 1:n
96 p(:,i) = g(t(i));
97 end
98
99 d = zeros(1,n-1);
100 for i = 1:n-1
101 d(i) = norm(p(:,i) - p(:,i+1));
102 end
103 end
104 end 40 end
105 41
106 function n = normal(obj,t) 42 function n = normal(obj,t)
107 deriv = obj.gp(t); 43 deriv = obj.gp(t);
108 normalization = sqrt(sum(deriv.^2,1)); 44 normalization = sqrt(sum(deriv.^2,1));
110 end 46 end
111 47
112 48
113 % Plots a curve g(t) for 0<t<1, using n points. Returns a handle h to the plotted curve. 49 % Plots a curve g(t) for 0<t<1, using n points. Returns a handle h to the plotted curve.
114 % h = plot_curve(g,n) 50 % h = plot_curve(g,n)
115 function h = plot(obj,n,marker) 51 function h = plot(obj,n)
116 default_arg('n',100); 52 default_arg('n',100);
117 default_arg('marker','line')
118 53
119 t = linspace(0,1,n); 54 t = linspace(0,1,n);
120
121 p = obj.g(t); 55 p = obj.g(t);
122 56 h = line(p(1,:),p(2,:));
123 switch marker
124 case 'line'
125 h = line(p(1,:),p(2,:));
126 otherwise
127 h = plot(p(1,:),p(2,:),marker);
128 end
129 end 57 end
130 58
131 % Plots the derivative gp(t) for 0<t<1, using n points. Returns a handle h to the plotted curve. 59 % Plots the derivative gp(t) for 0<t<1, using n points. Returns a handle h to the plotted curve.
132 % h = plot_curve(gp,n) 60 % h = plot_curve(gp,n)
133 function h = plot_derivative(obj,n,marker) 61 function h = plot_derivative(obj,n)
134 default_arg('n',100); 62 default_arg('n',100);
135 default_arg('marker','line')
136 63
137 t = linspace(0,1,n); 64 t = linspace(0,1,n);
138
139 p = obj.gp(t); 65 p = obj.gp(t);
140 66 h = line(p(1,:),p(2,:));
141 switch marker
142 case 'line'
143 h = line(p(1,:),p(2,:));
144 otherwise
145 h = plot(p(1,:),p(2,:),marker);
146 end
147 end 67 end
148 68
149 function h= plot_normals(obj,l,n,m) 69 function h= plot_normals(obj,l,n,m)
150 default_arg('l',0.1); 70 default_arg('l',0.1);
151 default_arg('n',10); 71 default_arg('n',10);
176 obj.plot(); 96 obj.plot();
177 end 97 end
178 % Shows curve with name and arrow for direction. 98 % Shows curve with name and arrow for direction.
179 99
180 100
181 function curve = stretch_parameter(obj,type) 101 % Arc length parameterization.
182 default_arg('type','arc_length'); 102 function curve = arcLengthStretch(obj,N)
183 switch type 103 default_arg('N',100);
184 % Arc length parameterization. 104 assert(~isempty(obj.gp),'This curve does not yet have a derivative added.')
185 case 'arc_length' 105
186 arcLength = obj.arc_length_fun; 106 % Construct arcLength function using splines
187 arcPar = @(t) util.fzero_vec(@(s)arcLength(s) - t*arcLength(1),[0-10*eps,1+10*eps]); 107 tvec = linspace(0,1,N);
188 g_new = @(t)obj.g(arcPar(t)); 108 arcVec = grid.Curve.arcLength(obj.gp,0,tvec);
189 gp_old = obj.gp; 109 arcLength = grid.Curve.spline(tvec,arcVec);
190 gp_new = @(t) normalize(gp_old(arcPar(t))); 110
191 111 % Stretch the parameter
192 arc_len_new = @(t) t; 112 arcPar = @(t) util.fzero_vec(@(s)arcLength(s) - t*arcLength(1),[0-10*eps,1+10*eps]);
193 curve = grid.Curve(g_new,gp_new,[],[],arc_len_new); 113
194 114 % New function and derivative
195 otherwise 115 g_new = @(t)obj.g(arcPar(t));
196 error('That stretching is not implemented.'); 116 gp_old = obj.gp;
197 end 117 gp_new = @(t) normalize(gp_old(arcPar(t)));
198 118
119 curve = grid.Curve(g_new,gp_new);
120
199 end 121 end
200 122
201 % how to make it work for methods without returns 123 % how to make it work for methods without returns
202 function p = subsref(obj,S) 124 function p = subsref(obj,S)
203 %Should i add error checking here? 125 %Should i add error checking here?
211 otherwise 133 otherwise
212 p = builtin('subsref',obj,S); 134 p = builtin('subsref',obj,S);
213 % error() 135 % error()
214 end 136 end
215 end 137 end
216
217
218 138
219 139
220 %% TRANSFORMATION OF A CURVE 140 %% TRANSFORMATION OF A CURVE
221 function D = reverse(C) 141 function D = reverse(C)
222 % g = C.g; 142 % g = C.g;
338 258
339 methods (Static) 259 methods (Static)
340 260
341 % Length of arc from parameter t0 to t1 (which may be vectors). 261 % Length of arc from parameter t0 to t1 (which may be vectors).
342 % Computed using derivative. 262 % Computed using derivative.
343 function L = arc_length(deriv,t0,t1) 263 function L = arcLength(deriv,t0,t1)
344 speed = @(t) sp(deriv(t)); 264 speed = @(t) sp(deriv(t));
345 265
346 function s = sp(deriv) 266 function s = sp(deriv)
347 s = sqrt(sum(deriv.^2,1)); 267 s = sqrt(sum(deriv.^2,1));
348 end 268 end