octave_packages/geometry-1.5.0/polygons2d/simplifypolyline.m

   1 %% Copyright (c) 2012 Juan Pablo Carbajal <carbajal@ifi.uzh.ch>
   2 %%
   3 %%    This program is free software: you can redistribute it and/or modify
   4 %%    it under the terms of the GNU General Public License as published by
   5 %%    the Free Software Foundation, either version 3 of the License, or
   6 %%    any later version.
   7 %%
   8 %%    This program is distributed in the hope that it will be useful,
   9 %%    but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 %%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  11 %%    GNU General Public License for more details.
  12 %%
  13 %%    You should have received a copy of the GNU General Public License
  14 %%    along with this program. If not, see <http://www.gnu.org/licenses/>.
  15
  16 %% -*- texinfo -*-
  17 %% @deftypefn {Function File} {[@var{pline2} @var{idx}] = } simplifypolyline (@var{pline})
  18 %% @deftypefnx {Function File} {@dots{} = } simplifypolyline (@dots{},@var{property},@var{value},@dots{})
  19 %% Simplify or subsample a polyline using the Ramer-Douglas-Peucker algorithm,
  20 %% a.k.a. the iterative end-point fit algorithm or the split-and-merge algorithm.
  21 %%
  22 %% The @var{pline} as a N-by-2 matrix. Rows correspond to the
  23 %% verices (compatible with @code{polygons2d}). The vector @var{idx} constains
  24 %% the indexes on vetices in @var{pline} that generates @var{pline2}, i.e.
  25 %% @code{pline2 = pline(idx,:)}.
  26 %%
  27 %% @strong{Parameters}
  28 %% @table @samp
  29 %% @item 'Nmax'
  30 %% Maximum number of vertices. Default value @code{1e3}.
  31 %% @item 'Tol'
  32 %% Tolerance for the error criteria. Default value @code{1e-4}.
  33 %% @item 'MaxIter'
  34 %% Maximum number of iterations. Default value @code{10}.
  35 %% @item 'Method'
  36 %% Not implemented.
  37 %% @end table
  38 %%
  39 %% Run @code{demo simplifypolyline} to see an example.
  40 %%
  41 %% @seealso{curve2polyline, curveval}
  42 %% @end deftypefn
  43
  44 function [pline idx] = simplifypolyline (pline_o, varargin)
  45 %% TODO do not print warnings if user provided Nmax or MaxIter.
  46
  47   # --- Parse arguments --- #
  48   parser = inputParser ();
  49   parser.FunctionName = "simplifypolyline";
  50   parser = addParamValue (parser,'Nmax', 1e3, @(x)x>0);
  51   toldef = 1e-4;%max(sumsq(diff(pline_o),2))*2;
  52   parser = addParamValue (parser,'Tol', toldef, @(x)x>0);
  53   parser = addParamValue (parser,'MaxIter', 100, @(x)x>0);
  54   parser = parse(parser,varargin{:});
  55
  56   Nmax      = parser.Results.Nmax;
  57   tol       = parser.Results.Tol;
  58   MaxIter   = parser.Results.MaxIter;
  59
  60   clear parser toldef
  61   msg = ["simplifypolyline: Maximum number of points reached with maximum error %g." ...
  62        " Increase %s if the result is not satisfactory."];
  63   # ------ #
  64
  65   [N dim] = size(pline_o);
  66   idx = [1 N];
  67
  68   for iter = 1:MaxIter
  69     % Find the point with the maximum distance.
  70     [dist ii] = maxdistance (pline_o, idx);
  71
  72     tf = dist > tol;
  73     n = sum(tf);
  74     if all(!tf);
  75       break;
  76     end
  77
  78     idx(end+1:end+n) = ii(tf);
  79     idx = sort(idx);
  80
  81     if length(idx) >= Nmax
  82       %% TODO remove extra points
  83       warning('geometry:MayBeWrongOutput', sprintf(msg,max(dist),'Nmax'));
  84       break;
  85     end
  86
  87   end
  88   if iter == MaxIter
  89     warning('geometry:MayBeWrongOutput', sprintf(msg,max(dist),'MaxIter'));
  90   end
  91
  92   pline = pline_o(idx,:);
  93 endfunction
  94
  95 function [dist ii] = maxdistance (p, idx)
  96
  97   %% Separate the groups of points according to the edge they can divide.
  98   func = @(x,y) x:y;
  99   idxc   = arrayfun (func, idx(1:end-1), idx(2:end), "UniformOutput",false);
 100   points = cellfun (@(x)p(x,:), idxc, "UniformOutput",false);
 101
 102   %% Build the edges
 103   edges = [p(idx(1:end-1),:) p(idx(2:end),:)];
 104   edges = mat2cell (edges, ones(1,size(edges,1)), 4)';
 105
 106   %% Calculate distance between the points and the corresponding edge
 107   [dist ii] = cellfun(@dd, points,edges,idxc);
 108
 109 endfunction
 110
 111 function [dist ii] = dd (p,e,idx)
 112   [d pos] = distancePointEdge(p,e);
 113   [dist ii] = max(d);
 114   ii = idx(ii);
 115 endfunction
 116
 117 %!demo
 118 %! t     = linspace(0,1,100).';
 119 %! y     = polyval([1 -1.5 0.5 0],t);
 120 %! pline = [t y];
 121 %!
 122 %! figure(1)
 123 %! clf
 124 %! plot (t,y,'-r;Original;','linewidth',2);
 125 %! hold on
 126 %!
 127 %! tol    = [8 2  1 0.5]*1e-2;
 128 %! colors = jet(4);
 129 %!
 130 %! for i=1:4
 131 %!   pline_ = simplifypolyline(pline,'tol',tol(i));
 132 %!   msg = sprintf('-;%g;',tol(i));
 133 %!   h = plot (pline_(:,1),pline_(:,2),msg);
 134 %!   set(h,'color',colors(i,:),'linewidth',2,'markersize',4);
 135 %! end
 136 %! hold off
 137 %!
 138 %! % ---------------------------------------------------------
 139 %! % Four approximations of the initial polyline with decreasing tolerances.
 140
 141 %!demo
 142 %! P       = [0 0; 3 1; 3 4; 1 3; 2 2; 1 1];
 143 %! func    = @(x,y) linspace(x,y,5);
 144 %! P2(:,1) = cell2mat( ...
 145 %!             arrayfun (func, P(1:end-1,1),P(2:end,1), ...
 146 %!             'uniformoutput',false))'(:);
 147 %! P2(:,2) = cell2mat( ...
 148 %!             arrayfun (func, P(1:end-1,2),P(2:end,2), ...
 149 %!             'uniformoutput',false))'(:);
 150 %!
 151 %! P2s = simplifypolyline (P2);
 152 %!
 153 %! plot(P(:,1),P(:,2),'s',P2(:,1),P2(:,2),'o',P2s(:,1),P2s(:,2),'-ok');
 154 %!
 155 %! % ---------------------------------------------------------
 156 %! % Simplification of a polyline in the plane.