octave_packages/plot-1.1.0/plotdecimate.m

   1 ## Copyright (C) 2006 Francesco Potortì <potorti@isti.cnr.it>
   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 2 of the License, or
   6 ## (at your option) 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, write to the Free Software Foundation,
  15 ## Inc., 51 Franklin Street, 5th Floor, Boston, MA 02110-1301, USA.
  16
  17 ## -*- texinfo -*-
  18 ## @deftypefn {Function File} {} plotdecimate (@var{P})
  19 ## @deftypefnx {Function File} {} plotdecimate (@var{P}, @var{so})
  20 ## @deftypefnx {Function File} {} plotdecimate (@var{P}, @var{so}, @var{res})
  21 ## Optimise plot data by removing redundant points and segments
  22 ##
  23 ## The first parameter @var{P} is a two-column matrix to be plotted as X and
  24 ## Y coordinates.   The second optional argument @var{so} disables segment
  25 ## optimisation when set to @var{false} (default is @var{true}). The third
  26 ## optional argument @var{res} is the size of the largest error on the plot:
  27 ## if it is a scalar, it is meant relative to the range of X and Y values
  28 ## (default 1e-3); if it is a 2x1 array, it contains the absolute errors for
  29 ## X and Y.  Returns a two-column matrix containing a subset of the rows of
  30 ## @var{P}. A line plot of @var{P} has the same appearance as a line plot of
  31 ## the output, with errors smaller than @var{res}.  When creating point
  32 ## plots, set @var{so} to @var{false}.
  33 ## @end deftypefn
  34
  35 function C = plotdecimate (P, so, res)
  36
  37   if (!ismatrix(P) || columns(P) != 2)
  38     error("P must be a matrix with two columns");
  39   endif
  40   if (nargin < 2)
  41     so = true;                  # do segment optimisation
  42   endif
  43   if (nargin < 3)
  44     res = 1e-3;                 # default resolution is 1000 dots/axis
  45   endif
  46
  47   ## Slack: admissible error on coordinates on the output plot
  48   if (isscalar(res))
  49     if (res <= 0)
  50       error("res must be positive");
  51     endif
  52     E = range(P)' * res;        # build error vector using range of data
  53   elseif (ismatrix(res))
  54     if (!all(size(res) == [2 1]) || any(res <= 0))
  55       error("res must be a 2x1 matrix with positive values");
  56     endif
  57     E = res;                    # take error vector as it is
  58   else
  59     error("res should be a scalar or matrix");
  60   endif
  61
  62   if (rows(P) < 3)
  63     C = P;
  64     return;                     # nothing to do
  65   endif
  66   P ./= repmat(E',rows(P),1);   # normalize P
  67   rot = [0,-1;1,0];             # rotate a vector pi/4 anticlockwise
  68
  69   ## Iteratively remove points too near to the previous point
  70   while (1)
  71     V = [true; sumsq(diff(P),2) > 1]; # points far from the previous ones
  72     if (all(V)) break; endif
  73     V = [true; diff(V) >= 0];   # identify the sequence leaders
  74     P = P(V,:);                 # remove them
  75   endwhile
  76
  77   ## Remove points laying near to a segment: for each segment R->S, build a
  78   ## unitary-lenght projection vector D perpendicular to R->S, and project
  79   ## R->T over D to compute the distance ot T from R->S.
  80   if (so)                       # segment optimisation
  81     ## For each segment, r and s are its extremes
  82     r = 1; R = P(1,:)';         # start of segment
  83     s = 2; S = P(2,:)';         # end of the segment
  84     rebuild = true;             # build first projection vector
  85
  86     for t = 3:rows(P)
  87       if (rebuild)              # build projection vector
  88         D = rot*(S-R)/sqrt(sumsq(S-R)); # projection vector for distance
  89         rebuild = false;        # keep current projection vector
  90       endif
  91
  92       T = P(t,:)';              # next point
  93
  94       if (abs(sum((T-R).*D)) < 1 # T is aligned
  95           && sum((T-R).*(S-R)) > 0) # going forward
  96         V(s) = false;           # do not plot s
  97       else                      # set a new segment
  98         r = s; R = S;           # new start of segment
  99         rebuild = true;         # rebuild projection vector
 100       endif
 101       s = t; S = T;             # new end of segment
 102     endfor
 103   endif
 104
 105   C = P(V,:) .* repmat(E',sum(V),1); # denormalize P
 106 endfunction
 107
 108 %!test
 109 %! x = [ 0 1 2 3 4 8 8 8 8 8 9 ]';
 110 %! y = [ 0 1 1 1 1 1 1 2 3 4 5 ]';
 111 %!
 112 %! x1 = [0 1 8 8 9]';
 113 %! y1 = [0 1 1 4 5]';
 114 %!   # optimised for segment plot
 115 %!
 116 %! x2 = [ 0 1 2 3 4 8 8 8 8 9 ]';
 117 %! y2 = [ 0 1 1 1 1 1 2 3 4 5 ]';
 118 %!   # double points removed
 119 %!
 120 %! P = [x,y];
 121 %!   # Original
 122 %! P1 = [x1, y1];
 123 %!   # optimised segments
 124 %! P2 = [x2, y2];
 125 %!   # double points removed
 126 %!
 127 %! assert(plotdecimate(P), P1);
 128 %! assert(plotdecimate(P, false), P2);