Re: plot save time scales with plot complexity

[Francesco Potortì]

>On 03/28/2012 05:43 PM, BobM wrote:
>
>> If I make the plot more complicated, but with the same number of points
>> t = (1:100000)
>> b = sin(2*pi*t/1000)
>> plot(t,b)
>> et to save ~ 19 sec
>> b = sin(2*pi*t/100) time to save as png goes to ~ 50 sec.
>
>For what it's worth, on my 2.3 GHz Linux box it takes 6 and 1.6 seconds 
>respectively as measured by tic(),print 'test.png' -dpng,toc().
>
>In any case, you are trying to plot 100,000 points in a space of 500 or 
>so pixels; I would suggest trimming your data

For cases like this, it is worthy to decimate data in nonuniform way,
depending on data density, so that the visual appearance of the plot is
not affected by variable data density.

The following code does it (optimally, I think).


===File /home/work/math/octavelib/plot/plotdecimate.m=======
## Copyright © 2006  Francesco Potortì
##
## 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

## Author: Francesco Potortì <address@hidden>
## $Revision: 2.7 $
## Usage: plotdecimate(P[, so[, res]])
## Description: Optimise plot data by removing redundant points and segments

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
============================================================

-- 
Francesco Potortì (ricercatore)        Voice:  +39.050.315.3058 (op.2111)
ISTI - Area della ricerca CNR          Mobile: +39.348.8283.107
via G. Moruzzi 1, I-56124 Pisa         Fax:    +39.050.315.2040  
(entrance 20, 1st floor, room C71)     Web:    http://fly.isti.cnr.it