xy_simplify

 XY_SIMPLIFY -- Simplify (reduce the number of points) of a polyline
  [xout, yout] = xy_simplify(x, y, nmax, tolerance) simplifies the (x, y)
   polyline, returning nmax points if (nmax > 0) or those
   points that are more significant than the given tolerance.
   The ls_sort() routine is used, in which deviation, not
   curvature, is the key principle.  NaNs are ignored.
  xy_simplify(...) returns [xout yout] in a single array.
  xy_simplify(nPoints) demonstrates itself with nPoints
   (random); default = 24.

0001 function [xout, yout] = xy_simplify(x, y, nmax, tolerance)
0002 % XY_SIMPLIFY -- Simplify (reduce the number of points) of a polyline
0003 %  [xout, yout] = xy_simplify(x, y, nmax, tolerance) simplifies the (x, y)
0004 %   polyline, returning nmax points if (nmax > 0) or those
0005 %   points that are more significant than the given tolerance.
0006 %   The ls_sort() routine is used, in which deviation, not
0007 %   curvature, is the key principle.  NaNs are ignored.
0008 %  xy_simplify(...) returns [xout yout] in a single array.
0009 %  xy_simplify(nPoints) demonstrates itself with nPoints
0010 %   (random); default = 24.
0011  
0012 % Copyright (C) 1997 Dr. Charles R. Denham, ZYDECO.
0013 %  All Rights Reserved.
0014 %   Disclosure without explicit written consent from the
0015 %    copyright owner does not constitute publication.
0016  
0017 % Version of 06-Jul-1998 06:40:26.
0018 % Updated    21-Sep-1998 17:13:27.
0019 
0020 % Reference: Douglas, D.H. and T.K. Peucker, Algorithms for
0021 %  the reduction of the number of points required to represent
0022 %  a digitized line of (sic?) its caricature, Can. Cartogr.,
0023 %  10, 112-122, 1973.  We have adopted a similar scheme for
0024 %  selecting "significant" points.  Unlike the D-P algorithm,
0025 %  we sort all the points first, then isolate those that satisfy
0026 %  our criteria for tolerance or total number of points retained.
0027 
0028 if nargin < 1, help(mfilename), x = 'demo'; end
0029 
0030 if isequal(x, 'demo'), x = 24; end
0031 
0032 if ischar(x), x = eval(x); end
0033 
0034 if length(x) == 1
0035         n = x;
0036         x = 1:n;
0037         y = cumsum(rand(size(x))-0.5);
0038         y(ceil(n/2)) = NaN;
0039         if (1), x(n) = x(1); y(n) = y(1); end
0040         nmax = fix((n+2) ./ 3);
0041         tolerance = 0;
0042         tic
0043         [xx, yy] = xy_simplify(x, y, nmax, tolerance);
0044         toc
0045         plot(xx, yy, '-*', x, y, 'o')
0046         xlabel('x'), ylabel('y')
0047         title('xy\_simplify')
0048         set(gcf, 'Name', 'xy_simplify')
0049         figure(gcf)
0050         return
0051 end
0052 
0053 if nargin < 3, nmax = 0; end
0054 if nargin < 4
0055         if nmax == 0, nmax = length(lon); end
0056         tolerance = 0;
0057 end
0058 
0059 if nmax == 0, nmax = length(x); end
0060 
0061 if nmax > 0 & nmax <= 1
0062         nmax = round(nmax*length(x));
0063 end
0064 
0065 % Save NaNs.
0066 
0067 xx = x(:);
0068 yy = y(:);
0069 
0070 i = find(isfinite(xx) & isfinite(yy));
0071 j = find(~isfinite(xx) | ~isfinite(yy));
0072 if any(i), xx = xx(i); yy = yy(i); end
0073 
0074 % Sort.
0075 
0076 [ind, dev] = ls_sort(xx, yy);
0077 
0078 % Restore NaNs: We need to decide how important
0079 %  the NaNs are -- most or least?  If least,
0080 %  then we do not need to include them.
0081 
0082 ind = i(ind);
0083 ind = [ind; j(:)];
0084 dev = [dev; NaN*j];
0085 
0086 if nmax > 0
0087         nmax = max(3, min(nmax, length(x)));
0088         ind = ind(1:nmax);
0089         dev = dev(1:nmax);
0090 else
0091         f = find(dev >= tolerance);
0092         if any(f)
0093                 ind = ind(f);
0094                 dev = dev(f);
0095         end
0096 end
0097 
0098 [ind, s] = sort(ind);
0099 dev = dev(s);
0100 
0101 xx = x(ind);
0102 yy = y(ind);
0103 
0104 if nargout == 1
0105         xout = [xx yy];
0106 elseif nargout > 1
0107         xout = xx;
0108         yout = yy;
0109 end
0110 
0111 function [theIndices, theDeviations] = ls_sort(x, y)
0112 
0113 % ls_sort -- Sort (x, y) data for line-simplification.
0114 %  [theIndices, theDeviations] = ls_sort(x, y) returns
0115 %   indices and deviations of the given (x, y) data,
0116 %   in order of decreasing importance for
0117 %   line-simplification.
0118  
0119 % Copyright (C) 1997 Dr. Charles R. Denham, ZYDECO.
0120 %  All Rights Reserved.
0121 %   Disclosure without explicit written consent from the
0122 %    copyright owner does not constitute publication.
0123  
0124 % Version of 02-Jul-1998 09:36:15.
0125 
0126 if nargin < 1, help(mfilename), x = 'demo'; end
0127 
0128 if isequal(x, 'demo'), x = 24; end
0129 
0130 if ischar(x), x = eval(x); end
0131 
0132 if length(x) == 1
0133         n = x;
0134         x = 1:n;
0135         y = cumsum(rand(size(x))-0.5);
0136         tic
0137         [i, d] = ls_sort(x, y);
0138         toc
0139         reduce = 3;   % Reduction-factor.
0140         k = fix((n+reduce-1)./reduce);
0141         ii = i(1:k);
0142         dd = d(1:k);
0143         xx = x(ii);
0144         yy = y(ii);
0145         [ii, s] = sort(ii);
0146         dd = dd(s);
0147         xx = xx(s);
0148         yy = yy(s);
0149         subplot(2, 1, 1)
0150         plot(x, y, 'bo', xx, yy, 'g-')
0151         title('Data and Subset')
0152         xlabel('x'), ylabel('y')
0153         subplot(2, 1, 2)
0154         plot(i, d, '.', xx, dd, 'o')
0155         title('Deviation')
0156         xlabel('x'), ylabel('y')
0157         set(gcf, 'name', 'ls_sort demo')
0158         figure(gcf)
0159         return
0160 end
0161 
0162 ind = zeros(size(x));
0163 dev = zeros(size(x));
0164 
0165 nstack = ceil(1 + 2*log(length(x))/log(2));
0166 stack = zeros(nstack, 2);
0167 nstack = 1;
0168 stack(nstack, :) = [1 length(x)];
0169 
0170 ind(1) = 1;
0171 ind(2) = length(x);
0172 k = 2;
0173 
0174 % Call ls_pivot() repeatedly and stack the
0175 %  longer interval before the shorter.
0176 
0177 while nstack > 0
0178         a = stack(nstack, 1);
0179         b = stack(nstack, 2);
0180         nstack = nstack - 1;
0181         if b-a > 1
0182                 [p, d] = ls_pivot(x([a:b]), y([a:b]));
0183                 if any(p)
0184                         p = a+p-1;
0185                         k = k+1;
0186                         ind(k) = p;
0187                         dev(k) = d;
0188                         if b-p >= p-a
0189                                 if b-p > 1
0190                                         nstack = nstack + 1;
0191                                         stack(nstack, :) = [p b];
0192                                 end
0193                         end
0194                         if p-a > 1
0195                                 nstack = nstack + 1;
0196                                 stack(nstack, :) = [a p];
0197                         end
0198                         if b-p < p-a
0199                                 if b-p > 1
0200                                         nstack = nstack + 1;
0201                                         stack(nstack, :) = [p b];
0202                                 end
0203                         end
0204                 end
0205         end
0206 end
0207 
0208 % Sort in descending order of importance.
0209 
0210 dev(1:2) = max(dev);
0211 [ignore, s] = sort(-dev);
0212 ind = ind(s);
0213 dev = dev(s);
0214 
0215 % Output.
0216 
0217 if nargout > 0
0218         theIndices = ind;
0219         theDeviations = dev;
0220 else
0221         ind, dev
0222 end
0223 
0224 function [thePivot, theDeviation] = ls_pivot(x, y)
0225 
0226 % ls_pivot -- Line-simplification pivot.
0227 %  ls_pivot(x, y) returns the index of the (x, y)
0228 %   point that deviates most from the straight-line
0229 %   connection between the first and last of the
0230 %   given points.
0231  
0232 % Copyright (C) 1997 Dr. Charles R. Denham, ZYDECO.
0233 %  All Rights Reserved.
0234 %   Disclosure without explicit written consent from the
0235 %    copyright owner does not constitute publication.
0236  
0237 % Version of 02-Jul-1998 08:32:40.
0238 
0239 if nargin < 1, help(mfilename), x = 'demo'; end
0240 
0241 if isequal(x, 'demo'), x = 10; end
0242 
0243 if ischar(x), x = eval(x); end
0244 
0245 if length(x) == 1
0246         n = x;
0247         x = 1:n;
0248         y = cumsum(rand(size(x))-0.5);
0249         tic
0250         [piv, dev] = ls_pivot(x, y);
0251         toc
0252         plot(x, y, 'b-', x([1 n]), y([1 n]), 'g-', x(piv), y(piv), '*')
0253         title(['ls\_pivot(' int2str(n) ')'])
0254         xlabel('x'), ylabel('y')
0255         figure(gcf)
0256         return
0257 end
0258 
0259 if nargin < 2, y = x; x = 1:length(y); end
0260 
0261 thePivot = 0;
0262 theDeviation = -inf;
0263 
0264 n = length(x);
0265 if n < 3, return, end
0266 
0267 % Translate to origin.
0268 
0269 x = x(:) - x(1);
0270 y = y(:) - y(1);
0271 
0272 % Rotate onto x-axis about origin, then
0273 %  select the furthest point in y.
0274 %  If end-points are the same, we instead pick
0275 %  the point most distant from point #1.
0276 
0277 if any([x(1) y(1)] ~= [x(n) y(n)])
0278         r = norm([x(n) y(n)]);
0279         c = x(n) ./ r;
0280         s = y(n) ./ r;
0281         z = [x y] * [c -s; s c];
0282 else
0283         z = x + sqrt(-1)*y;
0284         z = [z z];
0285 end
0286 
0287 % Find the pivot.
0288 
0289 temp = abs(z(2:n-1, 2));
0290 f = find(temp == max(temp));
0291 piv = f(1) + 1;
0292 dev = abs(z(piv));
0293 
0294 % Output.
0295 
0296 if nargout > 0
0297         thePivot = piv;
0298         theDeviation = dev;
0299 else
0300         piv, dev
0301 end