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