1######################################################################## 2## 3## Copyright (C) 2000-2021 The Octave Project Developers 4## 5## See the file COPYRIGHT.md in the top-level directory of this 6## distribution or <https://octave.org/copyright/>. 7## 8## This file is part of Octave. 9## 10## Octave is free software: you can redistribute it and/or modify it 11## under the terms of the GNU General Public License as published by 12## the Free Software Foundation, either version 3 of the License, or 13## (at your option) any later version. 14## 15## Octave is distributed in the hope that it will be useful, but 16## WITHOUT ANY WARRANTY; without even the implied warranty of 17## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 18## GNU General Public License for more details. 19## 20## You should have received a copy of the GNU General Public License 21## along with Octave; see the file COPYING. If not, see 22## <https://www.gnu.org/licenses/>. 23## 24######################################################################## 25 26## -*- texinfo -*- 27## @deftypefn {} {} voronoi (@var{x}, @var{y}) 28## @deftypefnx {} {} voronoi (@var{x}, @var{y}, @var{options}) 29## @deftypefnx {} {} voronoi (@dots{}, "linespec") 30## @deftypefnx {} {} voronoi (@var{hax}, @dots{}) 31## @deftypefnx {} {@var{h} =} voronoi (@dots{}) 32## @deftypefnx {} {[@var{vx}, @var{vy}] =} voronoi (@dots{}) 33## Plot the Voronoi diagram of points @code{(@var{x}, @var{y})}. 34## 35## The Voronoi facets with points at infinity are not drawn. 36## 37## The @var{options} argument, which must be a string or cell array of strings, 38## contains options passed to the underlying qhull command. 39## See the documentation for the Qhull library for details 40## @url{http://www.qhull.org/html/qh-quick.htm#options}. 41## 42## If @qcode{"linespec"} is given it is used to set the color and line style of 43## the plot. 44## 45## If an axes graphics handle @var{hax} is supplied then the Voronoi diagram is 46## drawn on the specified axes rather than in a new figure. 47## 48## If a single output argument is requested then the Voronoi diagram will be 49## plotted and a graphics handle @var{h} to the plot is returned. 50## 51## [@var{vx}, @var{vy}] = voronoi (@dots{}) returns the Voronoi vertices 52## instead of plotting the diagram. 53## 54## @example 55## @group 56## x = rand (10, 1); 57## y = rand (size (x)); 58## h = convhull (x, y); 59## [vx, vy] = voronoi (x, y); 60## plot (vx, vy, "-b", x, y, "o", x(h), y(h), "-g"); 61## legend ("", "points", "hull"); 62## @end group 63## @end example 64## 65## @seealso{voronoin, delaunay, convhull} 66## @end deftypefn 67 68function [vx, vy] = voronoi (varargin) 69 70 if (nargin < 1) 71 print_usage (); 72 endif 73 74 narg = 1; 75 hax = NaN; 76 if (isscalar (varargin{1}) && ishghandle (varargin{1})) 77 hax = varargin{1}; 78 if (! isaxes (hax)) 79 error ("voronoi: HAX argument must be an axes object"); 80 endif 81 narg += 1; 82 endif 83 84 if (nargin < 1 + narg || nargin > 3 + narg) 85 print_usage (); 86 endif 87 88 x = varargin{narg++}; 89 y = varargin{narg++}; 90 91 opts = {}; 92 if (narg <= nargin) 93 if (iscell (varargin{narg})) 94 opts = varargin(narg++); 95 elseif (isnumeric (varargin{narg})) 96 ## Accept, but ignore, the triangulation 97 narg += 1; 98 endif 99 endif 100 101 linespec = {"b"}; 102 if (narg <= nargin && ischar (varargin{narg})) 103 linespec = varargin(narg); 104 endif 105 106 if (! isvector (x) || ! isvector (y) || numel (x) != numel (y)) 107 error ("voronoi: X and Y must be vectors of the same length"); 108 elseif (numel (x) < 2) 109 error ("voronoi: minimum of 2 points required"); 110 endif 111 x = x(:); 112 y = y(:); 113 114 ## Add box to approximate rays to infinity. For Voronoi diagrams the 115 ## box should be close to the points themselves. To make the job of 116 ## finding the exterior edges easier it should be bigger than the area 117 ## enclosed by the points themselves. 118 ## NOTE: Octave uses a factor of 2 although we don't have an mathematical 119 ## justification for that. 120 121 xmin = min (x); 122 xmax = max (x); 123 ymin = min (y); 124 ymax = max (y); 125 ## Factor for size of bounding box 126 scale = 2; 127 xdelta = xmax - xmin; 128 ydelta = ymax - ymin; 129 xbox = [xmin - scale * xdelta; xmin - scale * xdelta; 130 xmax + scale * xdelta; xmax + scale * xdelta]; 131 ybox = [ymin - scale * ydelta; ymax + scale * ydelta; 132 ymax + scale * ydelta; ymin - scale * ydelta]; 133 134 [p, c, infi] = __voronoi__ ("voronoi", [[x; xbox], [y; ybox]], opts{:}); 135 136 ## Build list of edges from points in facet. 137 c = c(! infi).'; 138 edges = zeros (2, 0); 139 for i = 1:numel (c) 140 facet = c{i}; 141 if (isempty (facet)) 142 continue; 143 endif 144 edges = [edges, [facet; [facet(end), facet(1:end-1)]]]; 145 endfor 146 147 ## Keep only the unique edges of the Voronoi diagram 148 edges = sortrows (sort (edges).').'; 149 edges = edges(:, [any(diff(edges, 1, 2)), true]); 150 151 if (numel (x) > 2) 152 ## Eliminate the edges of the diagram representing the box. 153 ## Exclude points outside a certain radius from the center of distribution. 154 ## FIXME: Factor should be at least 1.0. Octave uses 1.1 for margin. 155 ## There is no theoretical justification for this choice. 156 ctr = [(xmax + xmin)/2 , (ymax + ymin)/2]; 157 radius = 1.1 * sumsq ([xmin, ymin] - ctr); 158 dist = sumsq (p - ctr, 2); 159 160 p_inside = (1:rows (p))(dist < radius); 161 edge_inside = any (ismember (edges, p_inside)); 162 edges = edges(:, edge_inside); 163 else 164 ## look for the edge between the two given points 165 for edge = edges 166 if (det ([[[1;1],p(edge,1:2)];1,x(1),y(1)]) 167 * det ([[[1;1],p(edge,1:2)];1,x(2),y(2)]) < 0) 168 edges = edge; 169 break; 170 endif 171 endfor 172 ## Use larger plot limits to make it more likely single bisector is shown. 173 xdelta = ydelta = max (xdelta, ydelta); 174 endif 175 176 ## Get points of the diagram 177 Vvx = reshape (p(edges, 1), size (edges)); 178 Vvy = reshape (p(edges, 2), size (edges)); 179 180 if (nargout < 2) 181 if (isnan (hax)) 182 hax = gca (); 183 endif 184 h = plot (hax, Vvx, Vvy, linespec{:}, x, y, '+'); 185 lim = [xmin, xmax, ymin, ymax]; 186 axis (lim + 0.1 * [[-1, 1] * xdelta, [-1, 1] * ydelta]); 187 if (nargout == 1) 188 vx = h; 189 endif 190 else 191 vx = Vvx; 192 vy = Vvy; 193 endif 194 195endfunction 196 197 198%!demo 199%! voronoi (rand (10,1), rand (10,1)); 200 201%!testif HAVE_QHULL 202%! phi = linspace (-pi, 3/4*pi, 8); 203%! [x,y] = pol2cart (phi, 1); 204%! [vx,vy] = voronoi (x,y); 205%! assert (vx(2,:), zeros (1, columns (vx)), eps); 206%! assert (vy(2,:), zeros (1, columns (vy)), eps); 207 208%!testif HAVE_QHULL <*40996> 209%! ## Special case of just 2 points 210%! x = [0 1]; y = [1 0]; 211%! [vx, vy] = voronoi (x,y); 212%! assert (vx, [-0.7; 1.7], eps); 213%! assert (vy, [-0.7; 1.7], eps); 214 215%!testif HAVE_QHULL <*38295> 216%! x = [1,2,3]; y = [2,3,1]; 217%! [vx, vy] = voronoi (x,y); 218%! assert (columns (vx), 3); 219 220%!testif HAVE_QHULL <*37270> 221%! ## Duplicate points can cause an internal error 222%! x = [1,2,3, 3]; y = [2,3,1, 1]; 223%! [vx, vy] = voronoi (x,y); 224 225 226## Input validation tests 227%!error voronoi () 228%!error voronoi (ones (3,1)) 229%!error voronoi (ones (3,1), ones (3,1), "invalid1", "invalid2", "invalid3") 230%!error <HAX argument must be an axes object> voronoi (0, ones (3,1), ones (3,1)) 231%!error <X and Y must be vectors of the same length> voronoi (ones (3,1), ones (4,1)) 232%!error <minimum of 2 points required> voronoi (2.5, 3.5) 233