-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathdpsimplify.m
244 lines (203 loc) · 6.44 KB
/
dpsimplify.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
function [ps,ix] = dpsimplify(p,tol)
% Recursive Douglas-Peucker Polyline Simplification, Simplify
%
% [ps,ix] = dpsimplify(p,tol)
%
% dpsimplify uses the recursive Douglas-Peucker line simplification
% algorithm to reduce the number of vertices in a piecewise linear curve
% according to a specified tolerance. The algorithm is also know as
% Iterative Endpoint Fit. It works also for polylines and polygons
% in higher dimensions.
%
% In case of nans (missing vertex coordinates) dpsimplify assumes that
% nans separate polylines. As such, dpsimplify treats each line
% separately.
%
% For additional information on the algorithm follow this link
% http://en.wikipedia.org/wiki/Ramer-Douglas-Peucker_algorithm
%
% Input arguments
%
% p polyline n*d matrix with n vertices in d
% dimensions.
% tol tolerance (maximal euclidean distance allowed
% between the new line and a vertex)
%
% Output arguments
%
% ps simplified line
% ix linear index of the vertices retained in p (ps = p(ix))
%
% Examples
%
% 1. Simplify line
%
% tol = 1;
% x = 1:0.1:8*pi;
% y = sin(x) + randn(size(x))*0.1;
% p = [x' y'];
% ps = dpsimplify(p,tol);
%
% plot(p(:,1),p(:,2),'k')
% hold on
% plot(ps(:,1),ps(:,2),'r','LineWidth',2);
% legend('original polyline','simplified')
%
% 2. Reduce polyline so that only knickpoints remain by
% choosing a very low tolerance
%
% p = [(1:10)' [1 2 3 2 4 6 7 8 5 2]'];
% p2 = dpsimplify(p,eps);
% plot(p(:,1),p(:,2),'k+--')
% hold on
% plot(p2(:,1),p2(:,2),'ro','MarkerSize',10);
% legend('original line','knickpoints')
%
% 3. Simplify a 3d-curve
%
% x = sin(1:0.01:20)';
% y = cos(1:0.01:20)';
% z = x.*y.*(1:0.01:20)';
% ps = dpsimplify([x y z],0.1);
% plot3(x,y,z);
% hold on
% plot3(ps(:,1),ps(:,2),ps(:,3),'k*-');
%
%
%
% Author: Wolfgang Schwanghart, 13. July, 2010.
% w.schwanghart[at]unibas.ch
if nargin == 0
help dpsimplify
return
end
error(nargchk(2, 2, nargin))
% error checking
if ~isscalar(tol) || tol<0;
error('tol must be a positive scalar')
end
% nr of dimensions
nrvertices = size(p,1);
dims = size(p,2);
% anonymous function for starting point and end point comparision
% using a relative tolerance test
compare = @(a,b) abs(a-b)/max(abs(a),abs(b)) <= eps;
% what happens, when there are NaNs?
% NaNs divide polylines.
Inan = any(isnan(p),2);
% any NaN at all?
Inanp = any(Inan);
% if there is only one vertex
if nrvertices == 1 || isempty(p);
ps = p;
ix = 1;
% if there are two
elseif nrvertices == 2 && ~Inanp;
% when the line has no vertices (except end and start point of the
% line) check if the distance between both is less than the tolerance.
% If so, return the center.
if dims == 2;
d = hypot(p(1,1)-p(2,1),p(1,2)-p(2,2));
else
d = sqrt(sum((p(1,:)-p(2,:)).^2));
end
if d <= tol;
ps = sum(p,1)/2;
ix = 1;
else
ps = p;
ix = [1;2];
end
elseif Inanp;
% case: there are nans in the p array
% --> find start and end indices of contiguous non-nan data
Inan = ~Inan;
sIX = strfind(Inan',[0 1])' + 1;
eIX = strfind(Inan',[1 0])';
if Inan(end)==true;
eIX = [eIX;nrvertices];
end
if Inan(1);
sIX = [1;sIX];
end
% calculate length of non-nan components
lIX = eIX-sIX+1;
% put each component into a single cell
c = mat2cell(p(Inan,:),lIX,dims);
% now call dpsimplify again inside cellfun.
if nargout == 2;
[ps,ix] = cellfun(@(x) dpsimplify(x,tol),c,'uniformoutput',false);
ix = cellfun(@(x,six) x+six-1,ix,num2cell(sIX),'uniformoutput',false);
else
ps = cellfun(@(x) dpsimplify(x,tol),c,'uniformoutput',false);
end
% write the data from a cell array back to a matrix
ps = cellfun(@(x) [x;nan(1,dims)],ps,'uniformoutput',false);
ps = cell2mat(ps);
ps(end,:) = [];
% ix wanted? write ix to a matrix, too.
if nargout == 2;
ix = cell2mat(ix);
end
else
% if there are no nans than start the recursive algorithm
ixe = size(p,1);
ixs = 1;
% logical vector for the vertices to be retained
I = true(ixe,1);
% call recursive function
p = simplifyrec(p,tol,ixs,ixe);
ps = p(I,:);
% if desired return the index of retained vertices
if nargout == 2;
ix = find(I);
end
end
% _________________________________________________________
function p = simplifyrec(p,tol,ixs,ixe)
% check if startpoint and endpoint are the same
% better comparison needed which included a tolerance eps
c1 = num2cell(p(ixs,:));
c2 = num2cell(p(ixe,:));
% same start and endpoint with tolerance
sameSE = all(cell2mat(cellfun(compare,c1(:),c2(:),'UniformOutput',false)));
if sameSE;
% calculate the shortest distance of all vertices between ixs and
% ixe to ixs only
if dims == 2;
d = hypot(p(ixs,1)-p(ixs+1:ixe-1,1),p(ixs,2)-p(ixs+1:ixe-1,2));
else
d = sqrt(sum(bsxfun(@minus,p(ixs,:),p(ixs+1:ixe-1,:)).^2,2));
end
else
% calculate shortest distance of all points to the line from ixs to ixe
% subtract starting point from other locations
pt = bsxfun(@minus,p(ixs+1:ixe,:),p(ixs,:));
% end point
a = pt(end,:)';
beta = (a' * pt')./(a'*a);
b = pt-bsxfun(@times,beta,a)';
if dims == 2;
% if line in 2D use the numerical more robust hypot function
d = hypot(b(:,1),b(:,2));
else
d = sqrt(sum(b.^2,2));
end
end
% identify maximum distance and get the linear index of its location
[dmax,ixc] = max(d);
ixc = ixs + ixc;
% if the maximum distance is smaller than the tolerance remove vertices
% between ixs and ixe
if dmax <= tol;
if ixs ~= ixe-1;
I(ixs+1:ixe-1) = false;
end
% if not, call simplifyrec for the segments between ixs and ixc (ixc
% and ixe)
else
p = simplifyrec(p,tol,ixs,ixc);
p = simplifyrec(p,tol,ixc,ixe);
end
end
end