octave_packages/m/geometry/voronoi.m

   1 ## Copyright (C) 2000-2012 Kai Habel
   2 ##
   3 ## This file is part of Octave.
   4 ##
   5 ## Octave is free software; you can redistribute it and/or modify it
   6 ## under the terms of the GNU General Public License as published by
   7 ## the Free Software Foundation; either version 3 of the License, or (at
   8 ## your option) any later version.
   9 ##
  10 ## Octave is distributed in the hope that it will be useful, but
  11 ## WITHOUT ANY WARRANTY; without even the implied warranty of
  12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13 ## General Public License for more details.
  14 ##
  15 ## You should have received a copy of the GNU General Public License
  16 ## along with Octave; see the file COPYING.  If not, see
  17 ## <http://www.gnu.org/licenses/>.
  18
  19 ## -*- texinfo -*-
  20 ## @deftypefn  {Function File} {} voronoi (@var{x}, @var{y})
  21 ## @deftypefnx {Function File} {} voronoi (@var{x}, @var{y}, @var{options})
  22 ## @deftypefnx {Function File} {} voronoi (@dots{}, "linespec")
  23 ## @deftypefnx {Function File} {} voronoi (@var{hax}, @dots{})
  24 ## @deftypefnx {Function File} {@var{h} =} voronoi (@dots{})
  25 ## @deftypefnx {Function File} {[@var{vx}, @var{vy}] =} voronoi (@dots{})
  26 ## Plot the Voronoi diagram of points @code{(@var{x}, @var{y})}.
  27 ## The Voronoi facets with points at infinity are not drawn.
  28 ##
  29 ## If "linespec" is given it is used to set the color and line style of the
  30 ## plot.  If an axis graphics handle @var{hax} is supplied then the Voronoi
  31 ## diagram is drawn on the specified axis rather than in a new figure.
  32 ##
  33 ## The @var{options} argument, which must be a string or cell array of strings,
  34 ## contains options passed to the underlying qhull command.
  35 ## See the documentation for the Qhull library for details
  36 ## @url{http://www.qhull.org/html/qh-quick.htm#options}.
  37 ##
  38 ## If a single output argument is requested then the Voronoi diagram will be
  39 ## plotted and a graphics handle @var{h} to the plot is returned.
  40 ## [@var{vx}, @var{vy}] = voronoi(@dots{}) returns the Voronoi vertices
  41 ## instead of plotting the diagram.
  42 ##
  43 ## @example
  44 ## @group
  45 ## x = rand (10, 1);
  46 ## y = rand (size (x));
  47 ## h = convhull (x, y);
  48 ## [vx, vy] = voronoi (x, y);
  49 ## plot (vx, vy, "-b", x, y, "o", x(h), y(h), "-g");
  50 ## legend ("", "points", "hull");
  51 ## @end group
  52 ## @end example
  53 ##
  54 ## @seealso{voronoin, delaunay, convhull}
  55 ## @end deftypefn
  56
  57 ## Author: Kai Habel <kai.habel@gmx.de>
  58 ## First Release: 20/08/2000
  59
  60 ## 2002-01-04 Paul Kienzle <pkienzle@users.sf.net>
  61 ## * limit the default graph to the input points rather than the whole diagram
  62 ## * provide example
  63 ## * use unique(x,"rows") rather than __unique_rows__
  64
  65 ## 2003-12-14 Rafael Laboissiere <rafael@laboissiere.net>
  66 ## Added optional fourth argument to pass options to the underlying
  67 ## qhull command
  68
  69 function [vx, vy] = voronoi (varargin)
  70
  71   if (nargin < 1)
  72     print_usage ();
  73   endif
  74
  75   narg = 1;
  76   if (isscalar (varargin{1}) && ishandle (varargin{1}))
  77     handl = varargin{1};
  78     if (! strcmp (get (handl, "type"), "axes"))
  79       error ("voronoi: expecting first argument to be an axes object");
  80     endif
  81     narg++;
  82   elseif (nargout < 2)
  83     handl = gca ();
  84   endif
  85
  86   if (nargin < 1 + narg || nargin > 3 + narg)
  87     print_usage ();
  88   endif
  89
  90   x = varargin{narg++};
  91   y = varargin{narg++};
  92
  93   opts = {};
  94   if (narg <= nargin)
  95     if (iscell (varargin{narg}))
  96       opts = varargin(narg++);
  97     elseif (isnumeric (varargin{narg}))
  98       ## Accept, but ignore, the triangulation
  99       narg++;
 100     endif
 101   endif
 102
 103   linespec = {"b"};
 104   if (narg <= nargin && ischar (varargin{narg}))
 105     linespec = varargin(narg);
 106   endif
 107
 108   lx = length (x);
 109   ly = length (y);
 110
 111   if (lx != ly)
 112     error ("voronoi: X and Y must be vectors of the same length");
 113   endif
 114
 115   ## Add box to approximate rays to infinity. For Voronoi diagrams the
 116   ## box can (and should) be close to the points themselves. To make the
 117   ## job of finding the exterior edges it should be at least two times the
 118   ## delta below however
 119   xmax = max (x(:));
 120   xmin = min (x(:));
 121   ymax = max (y(:));
 122   ymin = min (y(:));
 123   xdelta = xmax - xmin;
 124   ydelta = ymax - ymin;
 125   scale = 2;
 126
 127   xbox = [xmin - scale * xdelta; xmin - scale * xdelta; ...
 128           xmax + scale * xdelta; xmax + scale * xdelta];
 129   ybox = [xmin - scale * xdelta; xmax + scale * xdelta; ...
 130           xmax + scale * xdelta; xmin - scale * xdelta];
 131
 132   [p, c, infi] = __voronoi__ ("voronoi",
 133                               [[x(:) ; xbox(:)], [y(:); ybox(:)]],
 134                               opts{:});
 135
 136   idx = find (! infi);
 137   ll = length (idx);
 138   c = c(idx).';
 139   k = sum (cellfun ("length", c));
 140   edges = cell2mat (cellfun (@(x) [x ; [x(end), x(1:end-1)]], c,
 141                              "uniformoutput", false));
 142
 143   ## Identify the unique edges of the Voronoi diagram
 144   edges = sortrows (sort (edges).').';
 145   edges = edges (:, [(edges(1, 1: end - 1) != edges(1, 2 : end) | ...
 146                       edges(2, 1 :end - 1) != edges(2, 2 : end)), true]);
 147
 148   ## Eliminate the edges of the diagram representing the box
 149   poutside = (1 : rows(p)) ...
 150       (p (:, 1) < xmin - xdelta | p (:, 1) > xmax + xdelta | ...
 151        p (:, 2) < ymin - ydelta | p (:, 2) > ymax + ydelta);
 152   edgeoutside = ismember (edges (1, :), poutside) & ...
 153       ismember (edges (2, :), poutside);
 154   edges (:, edgeoutside) = [];
 155
 156   ## Get points of the diagram
 157   Vvx = reshape (p(edges, 1), size (edges));
 158   Vvy = reshape (p(edges, 2), size (edges));
 159
 160   if (nargout < 2)
 161     lim = [xmin, xmax, ymin, ymax];
 162     h = plot (handl, Vvx, Vvy, linespec{:}, x, y, '+');
 163     axis (lim + 0.1 * [[-1, 1] * (lim (2) - lim (1)), ...
 164                        [-1, 1] * (lim (4) - lim (3))]);
 165     if (nargout == 1)
 166       vx = h;
 167     endif
 168   else
 169     vx = Vvx;
 170     vy = Vvy;
 171   endif
 172
 173 endfunction
 174
 175
 176 %!demo
 177 %! voronoi (rand(10,1), rand(10,1));
 178
 179 %!testif HAVE_QHULL
 180 %! phi = linspace (-pi, 3/4*pi, 8);
 181 %! [x,y] = pol2cart (phi, 1);
 182 %! [vx,vy] = voronoi (x,y);
 183 %! assert(vx(2,:), zeros (1, columns (vx)), eps);
 184 %! assert(vy(2,:), zeros (1, columns (vy)), eps);
 185
 186 %% FIXME: Need input validation tests
 187