X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=octave_packages%2Fplot-1.1.0%2Fplotdecimate.m;fp=octave_packages%2Fplot-1.1.0%2Fplotdecimate.m;h=2bffbbc270fc58331e2e55891b9fb7706db6e854;hb=f5f7a74bd8a4900f0b797da6783be80e11a68d86;hp=0000000000000000000000000000000000000000;hpb=1705066eceaaea976f010f669ce8e972f3734b05;p=CreaPhase.git

diff --git a/octave_packages/plot-1.1.0/plotdecimate.m b/octave_packages/plot-1.1.0/plotdecimate.m
new file mode 100644
index 0000000..2bffbbc
--- /dev/null
+++ b/octave_packages/plot-1.1.0/plotdecimate.m
@@ -0,0 +1,128 @@
+## Copyright (C) 2006 Francesco Potortì <potorti@isti.cnr.it>
+##
+## This program is free software; you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 2 of the License, or
+## (at your option) any later version.
+##
+## This program is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with this program; if not, write to the Free Software Foundation,
+## Inc., 51 Franklin Street, 5th Floor, Boston, MA 02110-1301, USA.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {} plotdecimate (@var{P})
+## @deftypefnx {Function File} {} plotdecimate (@var{P}, @var{so})
+## @deftypefnx {Function File} {} plotdecimate (@var{P}, @var{so}, @var{res})
+## Optimise plot data by removing redundant points and segments
+##
+## The first parameter @var{P} is a two-column matrix to be plotted as X and
+## Y coordinates.   The second optional argument @var{so} disables segment
+## optimisation when set to @var{false} (default is @var{true}). The third
+## optional argument @var{res} is the size of the largest error on the plot:
+## if it is a scalar, it is meant relative to the range of X and Y values
+## (default 1e-3); if it is a 2x1 array, it contains the absolute errors for
+## X and Y.  Returns a two-column matrix containing a subset of the rows of
+## @var{P}. A line plot of @var{P} has the same appearance as a line plot of
+## the output, with errors smaller than @var{res}.  When creating point
+## plots, set @var{so} to @var{false}.
+## @end deftypefn
+
+function C = plotdecimate (P, so, res)
+
+  if (!ismatrix(P) || columns(P) != 2)
+    error("P must be a matrix with two columns");
+  endif
+  if (nargin < 2)
+    so = true;                  # do segment optimisation
+  endif
+  if (nargin < 3)
+    res = 1e-3;                 # default resolution is 1000 dots/axis
+  endif
+
+  ## Slack: admissible error on coordinates on the output plot
+  if (isscalar(res))
+    if (res <= 0)
+      error("res must be positive");
+    endif
+    E = range(P)' * res;        # build error vector using range of data
+  elseif (ismatrix(res))
+    if (!all(size(res) == [2 1]) || any(res <= 0))
+      error("res must be a 2x1 matrix with positive values");
+    endif
+    E = res;                    # take error vector as it is
+  else
+    error("res should be a scalar or matrix");
+  endif
+
+  if (rows(P) < 3)
+    C = P;
+    return;                     # nothing to do
+  endif
+  P ./= repmat(E',rows(P),1);   # normalize P
+  rot = [0,-1;1,0];             # rotate a vector pi/4 anticlockwise
+
+  ## Iteratively remove points too near to the previous point
+  while (1)
+    V = [true; sumsq(diff(P),2) > 1]; # points far from the previous ones
+    if (all(V)) break; endif
+    V = [true; diff(V) >= 0];   # identify the sequence leaders
+    P = P(V,:);                 # remove them
+  endwhile
+
+  ## Remove points laying near to a segment: for each segment R->S, build a
+  ## unitary-lenght projection vector D perpendicular to R->S, and project
+  ## R->T over D to compute the distance ot T from R->S.
+  if (so)                       # segment optimisation
+    ## For each segment, r and s are its extremes
+    r = 1; R = P(1,:)';         # start of segment
+    s = 2; S = P(2,:)';         # end of the segment
+    rebuild = true;             # build first projection vector
+
+    for t = 3:rows(P)
+      if (rebuild)              # build projection vector
+        D = rot*(S-R)/sqrt(sumsq(S-R)); # projection vector for distance
+        rebuild = false;        # keep current projection vector
+      endif
+
+      T = P(t,:)';              # next point
+
+      if (abs(sum((T-R).*D)) < 1 # T is aligned
+          && sum((T-R).*(S-R)) > 0) # going forward
+        V(s) = false;           # do not plot s
+      else                      # set a new segment
+        r = s; R = S;           # new start of segment
+        rebuild = true;         # rebuild projection vector
+      endif
+      s = t; S = T;             # new end of segment
+    endfor
+  endif
+
+  C = P(V,:) .* repmat(E',sum(V),1); # denormalize P
+endfunction
+
+%!test
+%! x = [ 0 1 2 3 4 8 8 8 8 8 9 ]';
+%! y = [ 0 1 1 1 1 1 1 2 3 4 5 ]';
+%!
+%! x1 = [0 1 8 8 9]';
+%! y1 = [0 1 1 4 5]';
+%!   # optimised for segment plot
+%!
+%! x2 = [ 0 1 2 3 4 8 8 8 8 9 ]';
+%! y2 = [ 0 1 1 1 1 1 2 3 4 5 ]';
+%!   # double points removed
+%!
+%! P = [x,y];
+%!   # Original
+%! P1 = [x1, y1];
+%!   # optimised segments
+%! P2 = [x2, y2];
+%!   # double points removed
+%!
+%! assert(plotdecimate(P), P1);
+%! assert(plotdecimate(P, false), P2);