--- a/scripts/ChangeLog	Sun Dec 28 17:39:20 2008 +0100
+++ b/scripts/ChangeLog	Mon Dec 29 23:41:18 2008 +0100
@@ -1,3 +1,7 @@
+2008-12-29  David Bateman  <dbateman@free.fr>
+
+	* goemetry/voronoi.m: Speed up and handle dense grids.
+
 2008-12-28  Jaroslav Hajek <highegg@gmail.com>
 
 	* miscellaneous/delete.m: Allow filename globs. Display warnings if

--- a/scripts/geometry/voronoi.m	Sun Dec 28 17:39:20 2008 +0100
+++ b/scripts/geometry/voronoi.m	Mon Dec 29 23:41:18 2008 +0100
@@ -106,44 +106,31 @@
     error ("voronoi: arguments must be vectors of same length");
   endif
 
-  ## Add big box to approximate rays to infinity
+  ## Add box to approximate rays to infinity. For Voronoi diagrams the
+  ## box can (and should) be close to the points themselves. To make the
+  ## job of finding the exterior edges it should be at least two times the
+  ## delta below however
   xmax = max (x(:));
   xmin = min (x(:));
   ymax = max (y(:));
   ymin = min (y(:));
   xdelta = xmax - xmin;
   ydelta = ymax - ymin;
-  scale = 10e4;
+  scale = 2;
 
   xbox = [xmin - scale * xdelta; xmin - scale * xdelta; ...
 	  xmax + scale * xdelta; xmax + scale * xdelta];
   ybox = [xmin - scale * xdelta; xmax + scale * xdelta; ...
 	  xmax + scale * xdelta; xmin - scale * xdelta];
-  x = [x(:) ; xbox(:)];
-  y = [y(:) ; ybox(:)];
 
-  [p, c, infi] = __voronoi__ ([x(:), y(:)], opts{:});
+  [p, c, infi] = __voronoi__ ([[x(:) ; xbox(:)], [y(:); ybox(:)]], opts{:});
 
   idx = find (!infi);
-
   ll = length (idx);
-  k = 0;r = 1;
-
-  for i = 1 : ll
-    k += length (c{idx(i)});
-  endfor
-
-  edges = zeros (2, k);
-
-  for i = 1 : ll
-    fac = c{idx(i)};
-    lf = length (fac);
-    fac = [fac, fac(1)];
-    fst = fac (1 : length(fac) - 1);
-    sec = fac(2 : length(fac));
-    edges (:, r : r + lf - 1) = [fst; sec];
-    r += lf;
-  endfor
+  c = c(idx).';
+  k = sum (cellfun ('length', 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).').';