# HG changeset patch
# User Rik <rik@octave.org>
# Date 1347737524 25200
# Node ID b99c52303d0bf22ca253adff218830d939973c4e
# Parent  c7fd43f5a89dc4dc916934fe6e63de0c75df3e7d
voronoi.m: Fix bug when there are multiple identical input points (bug #37270)
Also clean up coding style to reflect core Octave style.

* voronoi.m: Eliminate empty facets from output of __voronoi__ call before
proceeding with the rest of processing.

diff --git a/scripts/geometry/voronoi.m b/scripts/geometry/voronoi.m
--- a/scripts/geometry/voronoi.m
+++ b/scripts/geometry/voronoi.m
@@ -105,10 +105,7 @@
     linespec = varargin(narg);
   endif
 
-  lx = length (x);
-  ly = length (y);
-
-  if (lx != ly)
+  if (length (x) != length (y))
     error ("voronoi: X and Y must be vectors of the same length");
   endif
 
@@ -133,35 +130,33 @@
                               [[x(:) ; xbox(:)], [y(:); ybox(:)]],
                               opts{:});
 
-  idx = find (! infi);
-  ll = length (idx);
-  c = c(idx).';
-  k = sum (cellfun ("length", c));
+  c = c(! infi).';
+  ## Delete null entries which cause problems in next cellfun function
+  c(cellfun ("isempty", c)) = [];
   edges = cell2mat (cellfun (@(x) [x ; [x(end), x(1:end-1)]], c,
                              "uniformoutput", false));
 
   ## Identify the unique edges of the Voronoi diagram
   edges = sortrows (sort (edges).').';
-  edges = edges (:, [(edges(1, 1: end - 1) != edges(1, 2 : end) | ...
-                      edges(2, 1 :end - 1) != edges(2, 2 : end)), true]);
+  edges = edges(:, [(edges(1, 1 :end - 1) != edges(1, 2 : end) | ...
+                     edges(2, 1 :end - 1) != edges(2, 2 : end)), true]);
 
   ## Eliminate the edges of the diagram representing the box
-  poutside = (1 : rows (p)) ...
-      (p (:, 1) < xmin - xdelta | p (:, 1) > xmax + xdelta | ...
-       p (:, 2) < ymin - ydelta | p (:, 2) > ymax + ydelta);
-  edgeoutside = ismember (edges (1, :), poutside) & ...
-      ismember (edges (2, :), poutside);
-  edges (:, edgeoutside) = [];
+  poutside = (1:rows (p)) ...
+      (p(:, 1) < xmin - xdelta | p(:, 1) > xmax + xdelta | ...
+       p(:, 2) < ymin - ydelta | p(:, 2) > ymax + ydelta);
+  edgeoutside = ismember (edges(1, :), poutside) & ...
+                ismember (edges(2, :), poutside);
+  edges(:, edgeoutside) = [];
 
   ## Get points of the diagram
   Vvx = reshape (p(edges, 1), size (edges));
   Vvy = reshape (p(edges, 2), size (edges));
 
   if (nargout < 2)
+    h = plot (handl, Vvx, Vvy, linespec{:}, x, y, '+');
     lim = [xmin, xmax, ymin, ymax];
-    h = plot (handl, Vvx, Vvy, linespec{:}, x, y, '+');
-    axis (lim + 0.1 * [[-1, 1] * (lim (2) - lim (1)), ...
-                       [-1, 1] * (lim (4) - lim (3))]);
+    axis (lim + 0.1 * [[-1, 1] * xdelta, [-1, 1] * ydelta]);
     if (nargout == 1)
       vx = h;
     endif