Mercurial > octave
view scripts/geometry/voronoi.m @ 33577:2506c2d30b32 bytecode-interpreter tip
maint: Merge default to bytecode-interpreter
author | Arun Giridhar <arungiridhar@gmail.com> |
---|---|
date | Sat, 11 May 2024 18:49:01 -0400 |
parents | 2e484f9f1f18 |
children |
line wrap: on
line source
######################################################################## ## ## Copyright (C) 2000-2024 The Octave Project Developers ## ## See the file COPYRIGHT.md in the top-level directory of this ## distribution or <https://octave.org/copyright/>. ## ## This file is part of Octave. ## ## Octave 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 3 of the License, or ## (at your option) any later version. ## ## Octave 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 Octave; see the file COPYING. If not, see ## <https://www.gnu.org/licenses/>. ## ######################################################################## ## -*- texinfo -*- ## @deftypefn {} {} voronoi (@var{x}, @var{y}) ## @deftypefnx {} {} voronoi (@var{x}, @var{y}, @var{options}) ## @deftypefnx {} {} voronoi (@dots{}, "linespec") ## @deftypefnx {} {} voronoi (@var{hax}, @dots{}) ## @deftypefnx {} {@var{h} =} voronoi (@dots{}) ## @deftypefnx {} {[@var{vx}, @var{vy}] =} voronoi (@dots{}) ## Plot the Voronoi diagram of points @code{(@var{x}, @var{y})}. ## ## The Voronoi facets with points at infinity are not drawn. ## ## The @var{options} argument, which must be a string or cell array of strings, ## contains options passed to the underlying qhull command. ## See the documentation for the Qhull library for details ## @url{http://www.qhull.org/html/qh-quick.htm#options}. ## ## If @qcode{"linespec"} is given it is used to set the color and line style of ## the plot. ## ## If an axes graphics handle @var{hax} is supplied then the Voronoi diagram is ## drawn on the specified axes rather than in a new figure. ## ## If a single output argument is requested then the Voronoi diagram will be ## plotted and a graphics handle @var{h} to the plot is returned. ## ## [@var{vx}, @var{vy}] = voronoi (@dots{}) returns the Voronoi vertices ## instead of plotting the diagram. ## ## @example ## @group ## x = rand (10, 1); ## y = rand (size (x)); ## h = convhull (x, y); ## [vx, vy] = voronoi (x, y); ## plot (vx, vy, "-b", x, y, "o", x(h), y(h), "-g"); ## legend ("", "points", "hull"); ## @end group ## @end example ## ## @seealso{voronoin, delaunay, convhull} ## @end deftypefn function [vx, vy] = voronoi (varargin) if (nargin < 1) print_usage (); endif narg = 1; hax = NaN; if (isscalar (varargin{1}) && ishghandle (varargin{1})) hax = varargin{1}; if (! isaxes (hax)) error ("voronoi: HAX argument must be an axes object"); endif narg += 1; endif if (nargin < 1 + narg || nargin > 3 + narg) print_usage (); endif x = varargin{narg++}; y = varargin{narg++}; opts = {}; if (narg <= nargin) if (iscell (varargin{narg})) opts = varargin(narg++); elseif (isnumeric (varargin{narg})) ## Accept, but ignore, the triangulation narg += 1; endif endif linespec = {"b"}; if (narg <= nargin && ischar (varargin{narg})) linespec = varargin(narg); endif if (! isvector (x) || ! isvector (y) || numel (x) != numel (y)) error ("voronoi: X and Y must be vectors of the same length"); elseif (numel (x) < 2) error ("voronoi: minimum of 2 points required"); endif x = x(:); y = y(:); ## Add box to approximate rays to infinity. For Voronoi diagrams the ## box should be close to the points themselves. To make the job of ## finding the exterior edges easier it should be bigger than the area ## enclosed by the points themselves. ## NOTE: Octave uses a factor of 2 although we don't have an mathematical ## justification for that. xmin = min (x); xmax = max (x); ymin = min (y); ymax = max (y); ## Factor for size of bounding box scale = 2; xdelta = xmax - xmin; ydelta = ymax - ymin; xbox = [xmin - scale * xdelta; xmin - scale * xdelta; xmax + scale * xdelta; xmax + scale * xdelta]; ybox = [ymin - scale * ydelta; ymax + scale * ydelta; ymax + scale * ydelta; ymin - scale * ydelta]; [p, c, infi] = __voronoi__ ("voronoi", [[x; xbox], [y; ybox]], opts{:}); ## Build list of edges from points in facet. c = c(! infi).'; edges = zeros (2, 0); for i = 1:numel (c) facet = c{i}; if (isempty (facet)) continue; endif edges = [edges, [facet; [facet(end), facet(1:end-1)]]]; endfor ## Keep only the unique edges of the Voronoi diagram edges = sortrows (sort (edges).').'; edges = edges(:, [any(diff(edges, 1, 2)), true]); if (numel (x) > 2) ## Eliminate the edges of the diagram representing the box. ## Exclude points outside a certain radius from the center of distribution. ## FIXME: Factor should be at least 1.0. Octave uses 1.1 for margin. ## There is no theoretical justification for this choice. ctr = [(xmax + xmin)/2 , (ymax + ymin)/2]; radius = 1.1 * sumsq ([xmin, ymin] - ctr); dist = sumsq (p - ctr, 2); p_inside = (1:rows (p))(dist < radius); edge_inside = any (ismember (edges, p_inside)); edges = edges(:, edge_inside); else ## look for the edge between the two given points for edge = edges if (det ([[[1;1],p(edge,1:2)];1,x(1),y(1)]) * det ([[[1;1],p(edge,1:2)];1,x(2),y(2)]) < 0) edges = edge; break; endif endfor ## Use larger plot limits to make it more likely single bisector is shown. xdelta = ydelta = max (xdelta, ydelta); endif ## Get points of the diagram Vvx = reshape (p(edges, 1), size (edges)); Vvy = reshape (p(edges, 2), size (edges)); if (nargout < 2) if (isnan (hax)) hax = gca (); endif h = plot (hax, Vvx, Vvy, linespec{:}, x, y, '+'); lim = [xmin, xmax, ymin, ymax]; axis (lim + 0.1 * [[-1, 1] * xdelta, [-1, 1] * ydelta]); if (nargout == 1) vx = h; endif else vx = Vvx; vy = Vvy; endif endfunction %!demo %! voronoi (rand (10,1), rand (10,1)); %!testif HAVE_QHULL %! phi = linspace (-pi, 3/4*pi, 8); %! [x,y] = pol2cart (phi, 1); %! [vx,vy] = voronoi (x,y); %! assert (vx(2,:), zeros (1, columns (vx)), eps); %! assert (vy(2,:), zeros (1, columns (vy)), eps); %!testif HAVE_QHULL <*40996> %! ## Special case of just 2 points %! x = [0 1]; y = [1 0]; %! [vx, vy] = voronoi (x,y); %! assert (vx, [-0.7; 1.7], eps); %! assert (vy, [-0.7; 1.7], eps); %!testif HAVE_QHULL <*38295> %! x = [1,2,3]; y = [2,3,1]; %! [vx, vy] = voronoi (x,y); %! assert (columns (vx), 3); %!testif HAVE_QHULL <*37270> %! ## Duplicate points can cause an internal error %! x = [1,2,3, 3]; y = [2,3,1, 1]; %! [vx, vy] = voronoi (x,y); ## Input validation tests %!error <Invalid call> voronoi () %!error voronoi (ones (3,1)) %!error voronoi (ones (3,1), ones (3,1), "invalid1", "invalid2", "invalid3") %!error <HAX argument must be an axes object> voronoi (0, ones (3,1), ones (3,1)) %!error <X and Y must be vectors of the same length> voronoi (ones (3,1), ones (4,1)) %!error <minimum of 2 points required> voronoi (2.5, 3.5)