Mercurial > repos > public > sbplib

comparison +grid/old/curve_discretise.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
comparison
equal deleted inserted replaced
-1:000000000000 0:48b6fb693025
1 % Discretises the curve g with the smallest number of points such that all segments
2 % are shorter than h. If do_plot is true the points of the discretisation and
3 % the normals of the curve in those points are plotted.
4 %
5 % [t,p,d] = curve_discretise(g,h,do_plot)
6 %
7 % t is a vector of input values to g.
8 % p is a cector of points.
9 % d are the length of the segments.
10 function [t,p,d] = curve_discretise(g,h,do_plot)
11 default_arg('do_plot',false)
12
13 n = 10;
14
15 [t,p,d] = curve_discretise_n(g,n);
16
17 % ni = 0;
18 while any(d>h)
19 [t,p,d] = curve_discretise_n(g,n);
20 n = ceil(n*d(1)/h);
21 % ni = ni+1;
22 end
23
24 % nj = 0;
25 while all(d<h)
26 [t,p,d] = curve_discretise_n(g,n);
27 n = n-1;
28 % nj = nj+1;
29 end
30 [t,p,d] = curve_discretise_n(g,n+1);
31
32 % fprintf('ni = %d, nj = %d\n',ni,nj);
33
34 if do_plot
35 fprintf('n:%d max: %f min: %f\n', n, max(d),min(d));
36 p = grid.map_curve(g,t);
37 figure
38 show(g,t,h);
39 end
40
41 end
42
43 function [t,p,d] = curve_discretise_n(g,n)
44 t = linspace(0,1,n);
45 t = equalize_d(g,t);
46 d = D(g,t);
47 p = grid.map_curve(g,t);
48 end
49
50 function d = D(g,t)
51 p = grid.map_curve(g,t);
52
53 d = zeros(1,length(t)-1);
54 for i = 1:length(d)
55 d(i) = norm(p(:,i) - p(:,i+1));
56 end
57 end
58
59 function t = equalize_d(g,t)
60 d = D(g,t);
61 v = d-mean(d);
62 while any(abs(v)>0.01*mean(d))
63 dt = t(2:end)-t(1:end-1);
64 t(2:end) = t(2:end) - cumsum(dt.*v./d);
65
66 t = t/t(end);
67 d = D(g,t);
68 v = d-mean(d);
69 end
70 end
71
72
73 function show(g,t,hh)
74 p = grid.map_curve(g,t);
75
76
77
78 h = grid.plot_curve(g);
79 h.LineWidth = 2;
80 axis equal
81 hold on
82 h = plot(p(1,:),p(2,:),'.');
83 h.Color = [0.8500 0.3250 0.0980];
84 h.MarkerSize = 24;
85 hold off
86
87 n = grid.curve_normals(g,t);
88 hold on
89 for i = 1:length(t)
90 p0 = p(:,i);
91 p1 = p0 + hh*n(:,i);
92 l = [p0, p1];
93 h = plot(l(1,:),l(2,:));
94 h.Color = [0.8500 0.3250 0.0980];
95 end
96
97 end