# HG changeset patch # User Rik # Date 1319563043 25200 # Node ID 7ff0bdc3dc4c78fcf4d796b180721a955ef97fcd # Parent 3c3b74677fa0c4fe141540c0e2d96990c12f4038 Revamp geometry functions dependent on Qhull (Bug #34604, Bug #33346) * NEWS : Document new options being passed to Qhull * convhull.m, delaunay.m, delaunay3.m, delaunayn.m, voronoi.m, voronoin.m: Update docstrings. Put input validation first. Use same variable names as Matlab. Restore random state altered in demos. * __delaunayn__.cc: Use common syntax for parsing OPTIONS input. Add 'Qz' option to qhull command for 2D,3D data. Correctly free all Qhull memory and avoid segfault with non-simplicial facets. * __voronoi__.cc: Use common syntax for parsing OPTIONS input. Correctly free all Qhull memory. * convhulln.cc: Use common syntax for parsing OPTIONS input. Use Matlab-compatible options for qhull command. Correctly free all Qhull memory. Allow return of non-simplicial facets without causing a segfault. diff -r 3c3b74677fa0 -r 7ff0bdc3dc4c NEWS --- a/NEWS Mon Oct 24 16:16:50 2011 -0400 +++ b/NEWS Tue Oct 25 10:17:23 2011 -0700 @@ -46,6 +46,16 @@ in pattern and in string. Note that Matlab documents this behavior but the implementation does *not* always follow it. + ** Geometry functions derived from Qhull (convhull, delaunay, voronoi) + have been revamped. The options passed to the underlying qhull command + have been changed for better results or for Matlab compatibility. + + convhull : Default options are "Qt" for 2D, 3D, 4D inputs + Default options are "Qt Qx" for 5D and higher + delaunay : Default options are "Qt Qbb Qc Qz" for 2D and 3D inputs + Default options are "Qt Qbb Qc Qx" for 4D and higher + voronoi : No default arguments + ** Matlab-compatible preference functions: addpref getpref ispref rmpref setpref diff -r 3c3b74677fa0 -r 7ff0bdc3dc4c scripts/geometry/convhull.m --- a/scripts/geometry/convhull.m Mon Oct 24 16:16:50 2011 -0400 +++ b/scripts/geometry/convhull.m Tue Oct 25 10:17:23 2011 -0700 @@ -18,35 +18,44 @@ ## -*- texinfo -*- ## @deftypefn {Function File} {@var{H} =} convhull (@var{x}, @var{y}) -## @deftypefnx {Function File} {@var{H} =} convhull (@var{x}, @var{y}, @var{opt}) -## Return the index vector to the points of the enclosing convex hull. The -## data points are defined by the x and y vectors. +## @deftypefnx {Function File} {@var{H} =} convhull (@var{x}, @var{y}, @var{options}) +## Compute the convex hull of the set of points defined by the +## vectors @var{x} and @var{y}. The hull @var{H} is an index vector into +## the set of points and specifies which points form the enclosing hull. ## -## A third optional argument, which must be a string, contains extra options -## passed to the underlying qhull command. See the documentation for the -## Qhull library for details. +## An optional third 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}. +## The default option is @code{@{"Qt"@}}. ## -## @seealso{delaunay, convhulln} +## If @var{options} is not present or @code{[]} then the default arguments are +## used. Otherwise, @var{options} replaces the default argument list. +## To append user options to the defaults it is necessary to repeat the +## default arguments in @var{options}. Use a null string to pass no arguments. +## +## @seealso{convhulln, delaunay, voronoi} ## @end deftypefn ## Author: Kai Habel -function H = convhull (x, y, opt) +function H = convhull (x, y, options) if (nargin != 2 && nargin != 3) print_usage (); endif - if (isvector (x) && isvector (y) && length (x) == length (y)) - if (nargin == 2) - i = convhulln ([x(:), y(:)]); - elseif (ischar (opt) || iscell (opt)) - i = convhulln ([x(:), y(:)], opt); - else - error ("convhull: third argument must be a string or cell array of strings"); - endif + if (! (isvector (x) && isvector (y) && length (x) == length (y)) + && ! size_equal (x, y)) + error ("convhull: X and Y must be the same size"); + elseif (nargin == 3 && ! (ischar (options) || iscellstr (options))) + error ("convhull: OPTIONS must be a string or cell array of strings"); + endif + + if (nargin == 2) + i = convhulln ([x(:), y(:)]); else - error ("convhull: first two input arguments must be vectors of same size"); + i = convhulln ([x(:), y(:)], options); endif n = rows (i); @@ -71,16 +80,21 @@ endfor H(n + 1) = H(1); + endfunction -%!testif HAVE_QHULL -%! x = -3:0.5:3; -%! y = abs (sin (x)); -%! assert (convhull (x, y, {"s","Qci","Tcv","Pp"}), [1;7;13;12;11;10;4;3;2;1]) %!demo %! x = -3:0.05:3; %! y = abs (sin (x)); %! k = convhull (x, y); -%! plot (x(k),y(k),'r-',x,y,'b+'); +%! plot (x(k),y(k),"r-;convex hull;", x,y,"b+;points;"); %! axis ([-3.05, 3.05, -0.05, 1.05]); + +%!testif HAVE_QHULL +%! x = -3:0.5:3; +%! y = abs (sin (x)); +%! assert (convhull (x, y), [1;7;13;12;11;10;4;3;2;1]) + +%% FIXME: Need input validation tests + diff -r 3c3b74677fa0 -r 7ff0bdc3dc4c scripts/geometry/delaunay.m --- a/scripts/geometry/delaunay.m Mon Oct 24 16:16:50 2011 -0400 +++ b/scripts/geometry/delaunay.m Tue Oct 25 10:17:23 2011 -0700 @@ -17,74 +17,103 @@ ## . ## -*- texinfo -*- -## @deftypefn {Function File} {@var{tri} =} delaunay (@var{x}, @var{y}) -## @deftypefnx {Function File} {@var{tri} =} delaunay (@var{x}, @var{y}, @var{opt}) -## The return matrix of size [n, 3] contains a set triangles which are -## described by the indices to the data point x and y vector. -## The triangulation satisfies the Delaunay circum-circle criterion. -## No other data point is in the circum-circle of the defining triangle. +## @deftypefn {Function File} {} delaunay (@var{x}, @var{y}) +## @deftypefnx {Function File} {@var{tri} =} delaunay (@var{x}, @var{y}) +## @deftypefnx {Function File} {@var{tri} =} delaunay (@var{x}, @var{y}, @var{options}) +## Compute the Delaunay triangulation for a 2-D set of points. +## The return value @var{tri} is a set of triangles which satisfies the +## Delaunay circum-circle criterion, i.e., only a single data point from +## [@var{x}, @var{y}] is within the circum-circle of the defining triangle. +## +## The set of triangles @var{tri} is a matrix of size [n, 3]. Each +## row defines a triangle and the three columns are the three vertices +## of the triangle. The value of @code{@var{tri}(i,j)} is an index into +## @var{x} and @var{y} for the location of the j-th vertex of the i-th +## triangle. ## -## A third optional argument, which must be a string, contains extra options -## passed to the underlying qhull command. See the documentation for the -## Qhull library for details. +## An optional third 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}. +## The default options are @code{@{"Qt", "Qbb", "Qc", "Qz"@}}. +## +## If @var{options} is not present or @code{[]} then the default arguments are +## used. Otherwise, @var{options} replaces the default argument list. +## To append user options to the defaults it is necessary to repeat the +## default arguments in @var{options}. Use a null string to pass no arguments. +## +## If no output argument is specified the resulting Delaunay triangulation +## is plotted along with the original set of points. ## ## @example ## @group ## x = rand (1, 10); -## y = rand (size (x)); +## y = rand (1, 10); ## T = delaunay (x, y); -## X = [x(T(:,1)); x(T(:,2)); x(T(:,3)); x(T(:,1))]; -## Y = [y(T(:,1)); y(T(:,2)); y(T(:,3)); y(T(:,1))]; +## VX = [ x(T(:,1)); x(T(:,2)); x(T(:,3)); x(T(:,1)) ]; +## VY = [ y(T(:,1)); y(T(:,2)); y(T(:,3)); y(T(:,1)) ]; ## axis ([0,1,0,1]); -## plot (X, Y, "b", x, y, "r*"); +## plot (VX, VY, "b", x, y, "r*"); ## @end group ## @end example -## @seealso{voronoi, delaunay3, delaunayn} +## @seealso{delaunay3, delaunayn, convhull, voronoi} ## @end deftypefn ## Author: Kai Habel -function ret = delaunay (x, y, opt) +function tri = delaunay (x, y, options) if (nargin != 2 && nargin != 3) print_usage (); endif - if ((isvector (x) && isvector (y) && length (x) == length (y)) - || size_equal (x, y)) - if (nargin == 2) - tri = delaunayn ([x(:), y(:)]); - elseif (ischar (opt) || iscellstr (opt)) - tri = delaunayn ([x(:), y(:)], opt); - else - error ("delaunay: third argument must be a string"); - endif + if (! (isvector (x) && isvector (y) && length (x) == length (y)) + && ! size_equal (x, y)) + error ("delaunay: X and Y must be the same size"); + elseif (nargin == 3 && ! (ischar (options) || iscellstr (options))) + error ("delaunay: OPTIONS must be a string or cell array of strings"); + endif + + if (nargin == 2) + T = delaunayn ([x(:), y(:)]); else - error ("delaunay: first two input arguments must be matrices of same size"); + T = delaunayn ([x(:), y(:)], options); endif if (nargout == 0) x = x(:).'; y = y(:).'; - X = [x(tri(:,1)); x(tri(:,2)); x(tri(:,3)); x(tri(:,1))]; - Y = [y(tri(:,1)); y(tri(:,2)); y(tri(:,3)); y(tri(:,1))]; - plot(X, Y, 'b', x, y, 'r*'); + VX = [ x(T(:,1)); x(T(:,2)); x(T(:,3)); x(T(:,1)) ]; + VY = [ y(T(:,1)); y(T(:,2)); y(T(:,3)); y(T(:,1)) ]; + plot (VX, VY, "b", x, y, "r*"); else - ret = tri; + tri = T; endif + endfunction + +%!demo +%! old_state = rand ("state"); +%! restore_state = onCleanup (@() rand ("state", old_state)); +%! rand ("state", 1); +%! x = rand (1,10); +%! y = rand (1,10); +%! T = delaunay (x,y); +%! VX = [ x(T(:,1)); x(T(:,2)); x(T(:,3)); x(T(:,1)) ]; +%! VY = [ y(T(:,1)); y(T(:,2)); y(T(:,3)); y(T(:,1)) ]; +%! axis ([0,1,0,1]); +%! plot (VX,VY,"b", x,y,"r*"); + +%!testif HAVE_QHULL +%! x = [-1, 0, 1, 0]; +%! y = [0, 1, 0, -1]; +%! assert (sortrows (sort (delaunay (x, y), 2)), [1,2,4;2,3,4]); + %!testif HAVE_QHULL %! x = [-1, 0, 1, 0, 0]; %! y = [0, 1, 0, -1, 0]; -%! assert (sortrows (sort (delaunay (x, y), 2)), [1,2,5;1,4,5;2,3,5;3,4,5]) +%! assert (sortrows (sort (delaunay (x, y), 2)), [1,2,5;1,4,5;2,3,5;3,4,5]); -%!demo -%! rand ('state', 1); -%! x = rand(1,10); -%! y = rand(size(x)); -%! T = delaunay(x,y); -%! X = [ x(T(:,1)); x(T(:,2)); x(T(:,3)); x(T(:,1)) ]; -%! Y = [ y(T(:,1)); y(T(:,2)); y(T(:,3)); y(T(:,1)) ]; -%! axis([0,1,0,1]); -%! plot(X,Y,'b',x,y,'r*'); +%% FIXME: Need input validation tests + diff -r 3c3b74677fa0 -r 7ff0bdc3dc4c scripts/geometry/delaunay3.m --- a/scripts/geometry/delaunay3.m Mon Oct 24 16:16:50 2011 -0400 +++ b/scripts/geometry/delaunay3.m Tue Oct 25 10:17:23 2011 -0700 @@ -17,42 +17,61 @@ ## . ## -*- texinfo -*- -## @deftypefn {Function File} {@var{T} =} delaunay3 (@var{x}, @var{y}, @var{z}) -## @deftypefnx {Function File} {@var{T} =} delaunay3 (@var{x}, @var{y}, @var{z}, @var{opt}) -## A matrix of size [n, 4] is returned. Each row contains a -## set of tetrahedron which are -## described by the indices to the data point vectors (x,y,z). +## @deftypefn {Function File} {@var{tetr} =} delaunay3 (@var{x}, @var{y}, @var{z}) +## @deftypefnx {Function File} {@var{tetr} =} delaunay3 (@var{x}, @var{y}, @var{z}, @var{options}) +## Compute the Delaunay triangulation for a 3-D set of points. +## The return value @var{tetr} is a set of tetrahedrons which satisfies the +## Delaunay circum-circle criterion, i.e., only a single data point from +## [@var{x}, @var{y}, @var{z}] is within the circum-circle of the defining +## tetrahedron. ## -## A fourth optional argument, which must be a string or cell array of strings, -## contains extra options passed to the underlying qhull command. See the -## documentation for the Qhull library for details. -## @seealso{delaunay, delaunayn} +## The set of tetrahedrons @var{tetr} is a matrix of size [n, 4]. Each +## row defines a tetrahedron and the four columns are the four vertices +## of the tetrahedron. The value of @code{@var{tetr}(i,j)} is an index into +## @var{x}, @var{y}, @var{z} for the location of the j-th vertex of the i-th +## tetrahedron. +## +## An optional fourth 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}. +## The default options are @code{@{"Qt", "Qbb", "Qc", "Qz"@}}. +## +## If @var{options} is not present or @code{[]} then the default arguments are +## used. Otherwise, @var{options} replaces the default argument list. +## To append user options to the defaults it is necessary to repeat the +## default arguments in @var{options}. Use a null string to pass no arguments. +## +## @seealso{delaunay, delaunayn, convhull, voronoi} ## @end deftypefn ## Author: Kai Habel -function tetr = delaunay3 (x, y, z, opt) +function tetr = delaunay3 (x, y, z, options) - if (nargin != 3 && nargin != 4) + if (nargin < 3 || nargin > 4) print_usage (); endif - if (isvector (x) && isvector (y) &&isvector (z) - && length (x) == length (y) && length(x) == length (z)) - if (nargin == 3) - tetr = delaunayn ([x(:), y(:), z(:)]); - elseif (ischar (opt) || iscell (opt)) - tetr = delaunayn ([x(:), y(:), z(:)], opt); - else - error ("delaunay3: fourth argument must be a string or cell array of strings"); - endif + if (! (isvector (x) && isvector (y) && isvector (z) + && length (x) == length (y) && length(x) == length (z))) + error ("delaunay: X, Y, and Z must be the same size"); + elseif (nargin == 4 && ! (ischar (options) || iscellstr (options))) + error ("delaunay3: OPTIONS must be a string or cell array of strings"); + endif + + if (nargin == 3) + tetr = delaunayn ([x(:), y(:), z(:)]); else - error ("delaunay3: first three input arguments must be vectors of same size"); + tetr = delaunayn ([x(:), y(:), z(:)], options); endif endfunction + %!testif HAVE_QHULL %! x = [-1, -1, 1, 0, -1]; y = [-1, 1, 1, 0, -1]; z = [0, 0, 0, 1, 1]; %! assert (sortrows (sort (delaunay3 (x, y, z), 2)), [1,2,3,4;1,2,4,5]) +%% FIXME: Need input validation tests + diff -r 3c3b74677fa0 -r 7ff0bdc3dc4c scripts/geometry/delaunayn.m --- a/scripts/geometry/delaunayn.m Mon Oct 24 16:16:50 2011 -0400 +++ b/scripts/geometry/delaunayn.m Tue Oct 25 10:17:23 2011 -0700 @@ -17,66 +17,77 @@ ## . ## -*- texinfo -*- -## @deftypefn {Function File} {@var{t} =} delaunayn (@var{p}) -## @deftypefnx {Function File} {@var{t} =} delaunayn (@var{p}, @var{opt}) -## Form the Delaunay triangulation for a set of points. -## The Delaunay triangulation is a tessellation of the convex hull of the -## points such that no n-sphere defined by the n-triangles contains +## @deftypefn {Function File} {@var{T} =} delaunayn (@var{pts}) +## @deftypefnx {Function File} {@var{T} =} delaunayn (@var{pts}, @var{options}) +## Compute the Delaunay triangulation for an N-dimensional set of points. +## The Delaunay triangulation is a tessellation of the convex hull of a set +## of points such that no N-sphere defined by the N-triangles contains ## any other points from the set. -## The input matrix @var{p} of size @code{[n, dim]} contains @math{n} -## points in a space of dimension dim. The return matrix @var{t} has the -## size @code{[m, dim+1]}. It contains for each row a set of indices to -## the points, which describes a simplex of dimension dim. For example, -## a 2-D simplex is a triangle and 3-D simplex is a tetrahedron. +## +## The input matrix @var{pts} of size [n, dim] contains n points in a space of +## dimension dim. The return matrix @var{T} has size [m, dim+1]. Each row +## of @var{T} contains a set of indices back into the original set of points +## @var{pts} which describes a simplex of dimension dim. For example, a 2-D +## simplex is a triangle and 3-D simplex is a tetrahedron. ## -## Extra options for the underlying Qhull command can be specified by the -## second argument. This argument is a cell array of strings. The default -## options depend on the dimension of the input: +## An optional second 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}. +## The default options depend on the dimension of the input: ## ## @itemize -## @item 2D and 3D: @var{opt} = @code{@{"Qt", "Qbb", "Qc"@}} +## @item 2-D and 3-D: @var{options} = @code{@{"Qt", "Qbb", "Qc", "Qz"@}} ## -## @item 4D and higher: @var{opt} = @code{@{"Qt", "Qbb", "Qc", "Qz"@}} +## @item 4-D and higher: @var{options} = @code{@{"Qt", "Qbb", "Qc", "Qx"@}} ## @end itemize ## -## If @var{opt} is [], then the default arguments are used. If @var{opt} -## is @code{@{"@w{}"@}}, then none of the default arguments are used by Qhull. -## See the Qhull documentation for the available options. +## If @var{options} is not present or @code{[]} then the default arguments are +## used. Otherwise, @var{options} replaces the default argument list. +## To append user options to the defaults it is necessary to repeat the +## default arguments in @var{options}. Use a null string to pass no arguments. ## -## All options can also be specified as single string, for example -## @code{"Qt Qbb Qc Qz"}. -## +## @seealso{delaunay, delaunay3, convhulln, voronoin} ## @end deftypefn -function t = delaunayn (p, varargin) +function T = delaunayn (pts, varargin) + if (nargin < 1) print_usage (); endif - t = __delaunayn__ (p, varargin{:}); + T = __delaunayn__ (pts, varargin{:}); - if (isa (p, "single")) + if (isa (pts, "single")) myeps = eps ("single"); else myeps = eps; endif - ## Try to remove the zero volume simplices. The volume of the i-th simplex is - ## given by abs(det(p(t(i,1:end-1),:)-p(t(i,2:end),:)))/prod(1:n) - ## (reference http://en.wikipedia.org/wiki/Simplex). Any simplex with a - ## relative volume less than some arbitrary criteria is rejected. The + ## Try to remove the zero volume simplices. The volume of the i-th simplex is + ## given by abs(det(pts(T(i,1:end-1),:)-pts(T(i,2:end),:)))/prod(1:n) + ## (reference http://en.wikipedia.org/wiki/Simplex). Any simplex with a + ## relative volume less than some arbitrary criteria is rejected. The ## criteria we use is the volume of the simplex corresponding to an ## orthogonal simplex is equal edge length all equal to the edge length of - ## the original simplex. If the relative volume is 1e3*eps then the simplex - ## is rejected. Note division of the two volumes means that the factor + ## the original simplex. If the relative volume is 1e3*eps then the simplex + ## is rejected. Note division of the two volumes means that the factor ## prod(1:n) is dropped. idx = []; - [nt, n] = size (t); + [nt, n] = size (T); + ## FIXME: Vectorize this for loop or convert to delaunayn to .oct function for i = 1:nt - X = p(t(i,1:end-1),:) - p(t(i,2:end),:); - if (abs (det (X)) / sqrt (sum (X .^ 2, 2)) < 1e3 * myeps) - idx = [idx, i]; + X = pts(T(i,1:end-1),:) - pts(T(i,2:end),:); + if (abs (det (X)) / sqrt (sum (X .^ 2, 2)) < 1e3 * myeps) + idx(end+1) = i; endif endfor - t(idx,:) = []; + T(idx,:) = []; + endfunction + + +%% FIXME: Need tests for delaunayn + +%% FIXME: Need input validation tests + diff -r 3c3b74677fa0 -r 7ff0bdc3dc4c scripts/geometry/voronoi.m --- a/scripts/geometry/voronoi.m Mon Oct 24 16:16:50 2011 -0400 +++ b/scripts/geometry/voronoi.m Tue Oct 25 10:17:23 2011 -0700 @@ -18,18 +18,27 @@ ## -*- texinfo -*- ## @deftypefn {Function File} {} voronoi (@var{x}, @var{y}) -## @deftypefnx {Function File} {} voronoi (@var{x}, @var{y}, "plotstyle") -## @deftypefnx {Function File} {} voronoi (@var{x}, @var{y}, "plotstyle", @var{options}) +## @deftypefnx {Function File} {} voronoi (@var{x}, @var{y}, @var{options}) +## @deftypefnx {Function File} {} voronoi (@dots{}, "linespec") +## @deftypefnx {Function File} {} voronoi (@var{hax}, @dots{}) +## @deftypefnx {Function File} {@var{h} =} voronoi (@dots{}) ## @deftypefnx {Function File} {[@var{vx}, @var{vy}] =} voronoi (@dots{}) -## Plot Voronoi diagram of points @code{(@var{x}, @var{y})}. +## Plot the Voronoi diagram of points @code{(@var{x}, @var{y})}. ## The Voronoi facets with points at infinity are not drawn. -## [@var{vx}, @var{vy}] = voronoi(@dots{}) returns the vertices instead of -## plotting the -## diagram. plot (@var{vx}, @var{vy}) shows the Voronoi diagram. +## +## If "linespec" is given it is used to set the color and line style of the +## plot. If an axis graphics handle @var{hax} is supplied then the Voronoi +## diagram is drawn on the specified axis rather than in a new figure. ## -## A fourth optional argument, which must be a string, contains extra options -## passed to the underlying qhull command. See the documentation for the -## Qhull library for details. +## 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 a single output argument is requested then the Voronoi diagram will be +## plotted and a graphics handle to the plot is returned. +## [@var{vx}, @var{vy}] = voronoi(@dots{}) returns the Voronoi vertices +## instead of plotting the diagram. ## ## @example ## @group @@ -57,23 +66,23 @@ ## Added optional fourth argument to pass options to the underlying ## qhull command -function [vvx, vvy] = voronoi (varargin) +function [vx, vy] = voronoi (varargin) if (nargin < 1) print_usage (); + elseif (nargout > 2) + error ("voronoi: No more than two output arguments supported"); endif narg = 1; if (isscalar (varargin{1}) && ishandle (varargin{1})) handl = varargin{1}; - narg++; if (! strcmp (get (handl, "type"), "axes")) error ("voronoi: expecting first argument to be an axes object"); endif - else - if (nargout < 2) - handl = gca (); - endif + narg++; + elseif (nargout < 2) + handl = gca (); endif if (nargin < 1 + narg || nargin > 3 + narg) @@ -87,24 +96,22 @@ if (narg <= nargin) if (iscell (varargin{narg})) opts = varargin(narg++); - elseif (ismatrix (varargin{narg})) - ## Accept but ignore the triangulation + elseif (isnumeric (varargin{narg})) + ## Accept, but ignore, the triangulation narg++; endif endif linespec = {"b"}; - if (narg <= nargin) - if (ischar (varargin{narg})) - linespec = varargin(narg); - endif + if (narg <= nargin && ischar (varargin{narg})) + linespec = varargin(narg); endif lx = length (x); ly = length (y); if (lx != ly) - error ("voronoi: arguments must be vectors of same length"); + error ("voronoi: X and Y must be vectors of the same length"); endif ## Add box to approximate rays to infinity. For Voronoi diagrams the @@ -126,12 +133,12 @@ [p, c, infi] = __voronoi__ ([[x(:) ; xbox(:)], [y(:); ybox(:)]], opts{:}); - idx = find (!infi); + idx = find (! infi); ll = length (idx); c = c(idx).'; k = sum (cellfun ("length", c)); - edges = cell2mat(cellfun (@(x) [x ; [x(end), x(1:end-1)]], c, - "uniformoutput", false)); + 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).').'; @@ -147,32 +154,34 @@ edges (:, edgeoutside) = []; ## Get points of the diagram - vx = reshape (p (edges, 1), size(edges)); - vy = reshape (p (edges, 2), size(edges)); + Vvx = reshape (p(edges, 1), size (edges)); + Vvy = reshape (p(edges, 2), size (edges)); if (nargout < 2) lim = [xmin, xmax, ymin, ymax]; - h = plot (handl, vx, vy, linespec{:}, x, y, '+'); + h = plot (handl, Vvx, Vvy, linespec{:}, x, y, '+'); axis (lim + 0.1 * [[-1, 1] * (lim (2) - lim (1)), ... [-1, 1] * (lim (4) - lim (3))]); if (nargout == 1) - vxx = h; + vx = h; endif - elseif (nargout == 2) - vvx = vx; - vvy = vy; else - error ("voronoi: only two or zero output arguments supported"); + vx = Vvx; + vy = Vvy; endif endfunction -%!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,size(vx,2)),eps); -%! assert(vy(2,:),zeros(1,size(vy,2)),eps); %!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); + +%% FIXME: Need input validation tests + diff -r 3c3b74677fa0 -r 7ff0bdc3dc4c scripts/geometry/voronoin.m --- a/scripts/geometry/voronoin.m Mon Oct 24 16:16:50 2011 -0400 +++ b/scripts/geometry/voronoin.m Tue Oct 25 10:17:23 2011 -0700 @@ -20,14 +20,15 @@ ## @deftypefn {Function File} {[@var{C}, @var{F}] =} voronoin (@var{pts}) ## @deftypefnx {Function File} {[@var{C}, @var{F}] =} voronoin (@var{pts}, @var{options}) ## Compute N-dimensional Voronoi facets. The input matrix @var{pts} -## of size [n, dim] contains n points of dimension dim. +## of size [n, dim] contains n points in a space of dimension dim. ## @var{C} contains the points of the Voronoi facets. The list @var{F} -## contains for each facet the indices of the Voronoi points. +## contains, for each facet, the indices of the Voronoi points. ## -## A second optional argument, which must be a string, contains extra options -## passed to the underlying qhull command. See the documentation for the -## Qhull library for details. -## @seealso{voronoin, delaunay, convhull} +## An optional second 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}. +## @seealso{voronoi, convhulln, delaunayn} ## @end deftypefn ## Author: Kai Habel @@ -43,17 +44,24 @@ print_usage (); endif - [np, dims] = size (pts); - if (np > dims) - if (nargin == 1) - [C, F, infi] = __voronoi__ (pts); - elseif (ischar (options) || iscellstr (options)) - [C, F, infi] = __voronoi__ (pts, options); - else - error ("voronoin: second argument must be a string or cell array of strings"); - endif + [np, dim] = size (pts); + + if (np <= dim) + error ("voronoin: number of points must be greater than their dimension"); + elseif (nargin == 2 && ! (ischar (options) || iscellstr (options))) + error ("voronoin: OPTIONS argument must be a string or cell array of strings"); + endif + if (nargin == 1) + [C, F] = __voronoi__ (pts); else - error ("voronoin: number of points must be greater than their dimension"); + [C, F] = __voronoi__ (pts, options); endif + endfunction + + +%% FIXME: Need functional tests + +%% FIXME: Need input validation tests + diff -r 3c3b74677fa0 -r 7ff0bdc3dc4c src/DLD-FUNCTIONS/__delaunayn__.cc --- a/src/DLD-FUNCTIONS/__delaunayn__.cc Mon Oct 24 16:16:50 2011 -0400 +++ b/src/DLD-FUNCTIONS/__delaunayn__.cc Tue Oct 25 10:17:23 2011 -0700 @@ -62,8 +62,8 @@ DEFUN_DLD (__delaunayn__, args, , "-*- texinfo -*-\n\ -@deftypefn {Loadable Function} {@var{T} =} __delaunayn__ (@var{P})\n\ -@deftypefnx {Loadable Function} {@var{T} =} __delaunayn__ (@var{P}, @var{opt})\n\ +@deftypefn {Loadable Function} {@var{T} =} __delaunayn__ (@var{pts})\n\ +@deftypefnx {Loadable Function} {@var{T} =} __delaunayn__ (@var{pts}, @var{options})\n\ Undocumented internal function.\n\ @end deftypefn") @@ -73,7 +73,6 @@ #ifdef HAVE_QHULL retval(0) = 0.0; - std::string options = ""; int nargin = args.length (); if (nargin < 1 || nargin > 2) @@ -86,67 +85,52 @@ const octave_idx_type dim = p.columns (); const octave_idx_type n = p.rows (); - // default options + // Default options + std::string options; if (dim <= 3) - options = "Qt Qbb Qc"; + options = "Qt Qbb Qc Qz"; else options = "Qt Qbb Qc Qx"; - if (nargin == 2) { - if (args(1).is_empty ()) - { - // keep default options - } - else if (args(1).is_string ()) - { - // option string is directly provided + if (args(1).is_string ()) options = args(1).string_value (); - } - else if (args(1).is_cell ()) - { - options = ""; - - Cell c = args(1).cell_value (); - for (octave_idx_type i = 0; i < c.numel (); i++) - { + else if (args(1).is_empty ()) + ; // Use default options + else if (args(1).is_cellstr ()) + { + options = ""; + Array tmp = args(1).cellstr_value (); - if (! c.elem(i).is_string ()) - { - error ("__delaunayn__: all options must be strings"); - return retval; - } - - options = options + c.elem(i).string_value () + " "; - } - } - else - { - error ("__delaunayn__: OPT argument must be a string, cell array of strings, or empty"); - return retval; - } + for (octave_idx_type i = 0; i < tmp.numel (); i++) + options += tmp(i) + " "; + } + else + { + error ("__delaunayn__: OPTIONS argument must be a string, cell array of strings, or empty"); + return retval; + } } - //octave_stdout << "options " << options << std::endl; - if (n > dim + 1) { p = p.transpose (); double *pt_array = p.fortran_vec (); boolT ismalloc = false; - OCTAVE_LOCAL_BUFFER (char, flags, 250); + // Qhull flags argument is not const char* + OCTAVE_LOCAL_BUFFER (char, flags, 9 + options.length()); sprintf (flags, "qhull d %s", options.c_str ()); - // If you want some debugging information replace the 0 pointer - // with stdout or some other file open for writing. - + // Replace the 0 pointer with stdout for debugging information FILE *outfile = 0; FILE *errfile = stderr; - - if (! qh_new_qhull (dim, n, pt_array, ismalloc, flags, outfile, errfile)) + + int exitcode = qh_new_qhull (dim, n, pt_array, + ismalloc, flags, outfile, errfile); + if (! exitcode) { // triangulate non-simplicial facets qh_triangulate (); @@ -160,54 +144,48 @@ if (! facet->upperdelaunay) nf++; - // Double check + // Double check. Non-simplicial facets will cause segfault below if (! facet->simplicial) { error ("__delaunayn__: Qhull returned non-simplicial facets -- try delaunayn with different options"); + exitcode = 1; break; } } - Matrix simpl (nf, dim+1); - - FORALLfacets + if (! exitcode) { - if (! facet->upperdelaunay) - { - octave_idx_type j = 0; + Matrix simpl (nf, dim+1); - FOREACHvertex_ (facet->vertices) + FORALLfacets + { + if (! facet->upperdelaunay) { - // if delaunayn crashes, enable this check -#if 0 - if (j > dim) - { - error ("__delaunayn__: internal error. Qhull returned non-simplicial facets"); - return retval; - } -#endif + octave_idx_type j = 0; - simpl(i, j++) = 1 + qh_pointid(vertex->point); + FOREACHvertex_ (facet->vertices) + { + simpl(i, j++) = 1 + qh_pointid(vertex->point); + } + i++; } - i++; } - } - - retval(0) = simpl; - // free long memory - qh_freeqhull (! qh_ALL); - - // free short memory and memory allocator - int curlong, totlong; - qh_memfreeshort (&curlong, &totlong); - - if (curlong || totlong) - warning ("__delaunay__: did not free %d bytes of long memory (%d pieces)", - totlong, curlong); + retval(0) = simpl; + } } else error ("__delaunayn__: qhull failed"); + + // Free memory from Qhull + qh_freeqhull (! qh_ALL); + + int curlong, totlong; + qh_memfreeshort (&curlong, &totlong); + + if (curlong || totlong) + warning ("__delaunay__: did not free %d bytes of long memory (%d pieces)", + totlong, curlong); } else if (n == dim + 1) { diff -r 3c3b74677fa0 -r 7ff0bdc3dc4c src/DLD-FUNCTIONS/__voronoi__.cc --- a/src/DLD-FUNCTIONS/__voronoi__.cc Mon Oct 24 16:16:50 2011 -0400 +++ b/src/DLD-FUNCTIONS/__voronoi__.cc Tue Oct 25 10:17:23 2011 -0700 @@ -53,10 +53,11 @@ #endif #endif -DEFUN_DLD (__voronoi__, args, , +DEFUN_DLD (__voronoi__, args, nargout, "-*- texinfo -*-\n\ -@deftypefn {Loadable Function} {@var{tri} =} __voronoi__ (@var{point})\n\ -@deftypefnx {Loadable Function} {@var{tri} =} __voronoi__ (@var{point}, @var{options})\n\ +@deftypefn {Loadable Function} {@var{C}, @var{F} =} __voronoi__ (@var{pts})\n\ +@deftypefnx {Loadable Function} {@var{C}, @var{F} =} __voronoi__ (@var{pts}, @var{options})\n\ +@deftypefnx {Loadable Function} {@var{C}, @var{F}, @var{Inf_Pts} =} __voronoi__ (@dots{})\n\ Undocumented internal function.\n\ @end deftypefn") { @@ -73,57 +74,50 @@ return retval; } - std::string options = "qhull v Fv T0"; + std::string options = ""; if (nargin == 2) { - if (args(1).is_cellstr ()) + if (args(1).is_string ()) + options = args(1).string_value (); + else if (args(1).is_empty ()) + ; // Use default options + else if (args(1).is_cellstr ()) { + options = ""; Array tmp = args(1).cellstr_value (); for (octave_idx_type i = 0; i < tmp.numel (); i++) - options += " " + tmp(i); + options += tmp(i) + " "; } - else if (args(1).is_string ()) - options += " " + args(1).string_value (); else { - error ("__voronoi__: OPTIONS argument must be a string or cellstr"); + error ("__voronoi__: OPTIONS argument must be a string, cell array of strings, or empty"); return retval; } } Matrix p (args(0).matrix_value ()); - const octave_idx_type dim = p.columns (); const octave_idx_type np = p.rows (); + p = p.transpose (); - double *pt_array = p.fortran_vec (); - //double pt_array[dim * np]; - //for (int i = 0; i < np; i++) - // { - // for (int j = 0; j < dim; j++) - // { - // pt_array[j+i*dim] = p(i,j); - // } - // } - boolT ismalloc = false; - // If you want some debugging information replace the 0 pointer - // with stdout or some other file open for writing. - + // Replace the 0 pointer with stdout for debugging information FILE *outfile = 0; FILE *errfile = stderr; - // Qhull flags argument is not const char*... - OCTAVE_LOCAL_BUFFER (char, flags, options.length () + 1); + // Qhull flags argument is not const char* + OCTAVE_LOCAL_BUFFER (char, flags, 12 + options.length ()); - strcpy (flags, options.c_str ()); + sprintf (flags, "qhull v FV %s", options.c_str ()); - if (! qh_new_qhull (dim, np, pt_array, ismalloc, flags, outfile, errfile)) + int exitcode = qh_new_qhull (dim, np, pt_array, + ismalloc, flags, outfile, errfile); + if (! exitcode) { facetT *facet; vertexT *vertex; @@ -134,6 +128,7 @@ for (i = 0; i < np; i++) ni[i] = 0; qh_setvoronoi_all (); + bool infinity_seen = false; facetT *neighbor, **neighborp; coordT *voronoi_vertex; @@ -147,6 +142,7 @@ { if (qh hull_dim == 3) qh_order_vertexneighbors (vertex); + infinity_seen = false; FOREACHneighbor_ (vertex) @@ -173,9 +169,7 @@ for (octave_idx_type d = 0; d < dim; d++) v(0,d) = octave_Inf; - boolMatrix AtInf (np, 1); - for (i = 0; i < np; i++) - AtInf(i) = false; + boolMatrix AtInf (np, 1, false); octave_value_list F (np, octave_value ()); k = 0; i = 0; @@ -189,6 +183,7 @@ { if (qh hull_dim == 3) qh_order_vertexneighbors(vertex); + infinity_seen = false; RowVector facet_list (ni[k++]); m = 0; @@ -229,24 +224,27 @@ for (i = 0; i < r; i++) C.elem (i) = F(i); - retval(0) = v; + if (nargout == 3) + { + AtInf.resize (r, 1); + retval(2) = AtInf; + } retval(1) = C; - AtInf.resize (r, 1); - retval(2) = AtInf; - - // free long memory - qh_freeqhull (! qh_ALL); - - // free short memory and memory allocator - int curlong, totlong; - qh_memfreeshort (&curlong, &totlong); - - if (curlong || totlong) - warning ("__voronoi__: did not free %d bytes of long memory (%d pieces)", totlong, curlong); + retval(0) = v; } else error ("__voronoi__: qhull failed"); + // free memory from Qhull + qh_freeqhull (! qh_ALL); + + int curlong, totlong; + qh_memfreeshort (&curlong, &totlong); + + if (curlong || totlong) + warning ("__voronoi__: did not free %d bytes of long memory (%d pieces)", + totlong, curlong); + #else error ("__voronoi__: not available in this version of Octave"); #endif diff -r 3c3b74677fa0 -r 7ff0bdc3dc4c src/DLD-FUNCTIONS/convhulln.cc --- a/src/DLD-FUNCTIONS/convhulln.cc Mon Oct 24 16:16:50 2011 -0400 +++ b/src/DLD-FUNCTIONS/convhulln.cc Tue Oct 25 10:17:23 2011 -0700 @@ -53,24 +53,39 @@ DEFUN_DLD (convhulln, args, nargout, "-*- texinfo -*-\n\ -@deftypefn {Loadable Function} {@var{h} =} convhulln (@var{p})\n\ -@deftypefnx {Loadable Function} {@var{h} =} convhulln (@var{p}, @var{opt})\n\ +@deftypefn {Loadable Function} {@var{h} =} convhulln (@var{pts})\n\ +@deftypefnx {Loadable Function} {@var{h} =} convhulln (@var{pts}, @var{options})\n\ @deftypefnx {Loadable Function} {[@var{h}, @var{v}] =} convhulln (@dots{})\n\ -Return an index vector to the points of the enclosing convex hull.\n\ -The input matrix of size [n, dim] contains n points of dimension dim.\n\n\ -If a second optional argument is given, it must be a string or cell array\n\ -of strings containing options for the underlying qhull command. (See\n\ -the Qhull documentation for the available options.) The default options\n\ -are \"s Qci Tcv\".\n\ -If the second output @var{v} is requested the volume of the convex hull is\n\ -calculated.\n\n\ -@seealso{convhull, delaunayn}\n\ +Compute the convex hull of the set of points @var{pts} which is a matrix\n\ +of size [n, dim] containing n points in a space of dimension dim.\n\ +The hull @var{h} is an index vector into the set of points and specifies\n\ +which points form the enclosing hull.\n\ +\n\ +An optional second argument, which must be a string or cell array of strings,\n\ +contains options passed to the underlying qhull command.\n\ +See the documentation for the Qhull library for details\n\ +@url{http://www.qhull.org/html/qh-quick.htm#options}.\n\ +The default options depend on the dimension of the input:\n\ +\n\ +@itemize\n\ +@item 2D, 3D, 4D: @var{options} = @code{@{\"Qt\"@}}\n\ +\n\ +@item 5D and higher: @var{options} = @code{@{\"Qt\", \"Qx\"@}}\n\ +@end itemize\n\ +\n\ +If @var{options} is not present or @code{[]} then the default arguments are\n\ +used. Otherwise, @var{options} replaces the default argument list.\n\ +To append user options to the defaults it is necessary to repeat the\n\ +default arguments in @var{options}. Use a null string to pass no arguments.\n\ +\n\ +If the second output @var{v} is requested the volume of the enclosing\n\ +convex hull is calculated.\n\n\ +@seealso{convhull, delaunayn, voronoin}\n\ @end deftypefn") { octave_value_list retval; #ifdef HAVE_QHULL - std::string options; int nargin = args.length (); if (nargin < 1 || nargin > 2) @@ -79,85 +94,88 @@ return retval; } + Matrix p (args(0).matrix_value ()); + const octave_idx_type dim = p.columns (); + const octave_idx_type n = p.rows (); + + // Default options + std::string options; + if (dim <= 4) + options = "Qt"; + else + options = "Qt Qx"; + if (nargin == 2) { - if (args (1).is_string ()) + if (args(1).is_string ()) options = args(1).string_value (); - else if (args(1).is_cell ()) + else if (args(1).is_empty ()) + ; // Use default options + else if (args(1).is_cellstr ()) { - Cell c = args(1).cell_value (); options = ""; - for (octave_idx_type i = 0; i < c.numel (); i++) - { - if (! c.elem(i).is_string ()) - { - error ("convhulln: OPT must be a string or cell array of strings"); - return retval; - } + Array tmp = args(1).cellstr_value (); - options = options + c.elem(i).string_value() + " "; - } + for (octave_idx_type i = 0; i < tmp.numel (); i++) + options += tmp(i) + " "; } else { - error ("convhulln: OPT must be a string or cell array of strings"); + error ("convhulln: OPTIONS must be a string, cell array of strings, or empty"); return retval; } - } - else - // turn on some consistency checks - options = "s Qci Tcv"; + } - Matrix p (args(0).matrix_value ()); - const octave_idx_type dim = p.columns (); - const octave_idx_type n = p.rows (); p = p.transpose (); - double *pt_array = p.fortran_vec (); - - boolT ismalloc = False; - - std::ostringstream buf; + boolT ismalloc = false; - buf << "qhull QJ " << options; - - std::string buf_string = buf.str (); + // FIXME: we can't just pass options.c_str () to qh_new_qhull + // because the argument is not declared const. Ugh. Unless qh_new_qhull + // really needs to modify this argument, someone should fix QHULL. + OCTAVE_LOCAL_BUFFER (char, flags, 7 + options.length ()); - // FIXME -- we can't just pass buf_string.c_str () to qh_new_qhull - // because the argument is not declared const. Ugh. Unless - // qh_new_qhull really needs to modify this argument, someone should - // fix QHULL. + sprintf (flags, "qhull %s", options.c_str ()); - OCTAVE_LOCAL_BUFFER (char, flags, buf_string.length () + 1); - - strcpy (flags, buf_string.c_str ()); - - if (! qh_new_qhull (dim, n, pt_array, ismalloc, flags, 0, stderr)) + // Replace the 0 pointer with stdout for debugging information + FILE *outfile = 0; + FILE *errfile = stderr; + + int exitcode = qh_new_qhull (dim, n, pt_array, + ismalloc, flags, outfile, errfile); + if (! exitcode) { - // If you want some debugging information replace the NULL - // pointer with stdout - vertexT *vertex, **vertexp; facetT *facet; setT *vertices; - unsigned int nf = qh num_facets; + bool nonsimp_seen = false; + octave_idx_type nf = qh num_facets; - Matrix idx (nf, dim); + Matrix idx (nf, dim + 1); - octave_idx_type j, i = 0; + octave_idx_type i = 0, j; FORALLfacets { j = 0; - if (! facet->simplicial) - // should never happen with QJ - error ("convhulln: non-simplicial facet"); + + if (! nonsimp_seen && ! facet->simplicial) + { + nonsimp_seen = true; + if (options.find ("QJ") != std::string::npos) + { + // should never happen with QJ + error ("convhulln: qhull failed. Option 'QJ' returned non-simplicial facet"); + break; + } + } if (dim == 3) { vertices = qh_facet3vertex (facet); FOREACHvertex_ (vertices) idx(i, j++) = 1 + qh_pointid(vertex->point); + qh_settempfree (&vertices); } else @@ -174,11 +192,15 @@ } } if (j < dim) - // likewise but less fatal warning ("facet %d only has %d vertices", i, j); + i++; } + // Remove extra dimension if all facets were simplicial + if (! nonsimp_seen) + idx.resize (nf, dim, 0.0); + if (nargout == 2) // calculate volume of convex hull // taken from qhull src/geom2.c @@ -213,19 +235,21 @@ retval(1) = octave_value (qh totvol); } - retval(0) = octave_value (idx); + retval(0) = idx; } + else + error ("convhulln: qhull failed"); - // free long memory + // free memory from Qhull qh_freeqhull (! qh_ALL); - // free short memory and memory allocator int curlong, totlong; qh_memfreeshort (&curlong, &totlong); if (curlong || totlong) warning ("convhulln: did not free %d bytes of long memory (%d pieces)", - totlong, curlong); + totlong, curlong); + #else error ("convhulln: not available in this version of Octave"); #endif @@ -236,10 +260,23 @@ /* %!testif HAVE_QHULL %! cube = [0 0 0;1 0 0;1 1 0;0 1 0;0 0 1;1 0 1;1 1 1;0 1 1]; -%! [h, v] = convhulln(cube,'Pp'); +%! [h, v] = convhulln (cube); +%! assert (size (h), [6 4]); +%! h = sortrows (sort (h, 2), [1:4]); +%! assert (h, [1 2 3 4; 1 2 5 6; 1 4 5 8; 2 3 6 7; 3 4 7 8; 5 6 7 8]); +%! assert (v, 1, 10*eps); + +%!testif HAVE_QHULL +%! cube = [0 0 0;1 0 0;1 1 0;0 1 0;0 0 1;1 0 1;1 1 1;0 1 1]; +%! [h, v] = convhulln (cube, "QJ"); +%! assert (size (h), [12 3]); +%! assert (sortrows (sort (h, 2), [1:3]), [1 2 4; 1 2 5; 1 4 5; 2 3 4; 2 3 6; 2 5 6; 3 4 8; 3 6 7; 3 7 8; 4 5 8; 5 6 8; 6 7 8]); %! assert (v, 1.0, 1e6*eps); + %!testif HAVE_QHULL %! tetrahedron = [1 1 1;-1 -1 1;-1 1 -1;1 -1 -1]; -%! [h, v] = convhulln(tetrahedron,'Pp'); -%! assert (v, 8/3, 1e6*eps); +%! [h, v] = convhulln (tetrahedron); +%! h = sortrows (sort (h, 2), [1 2 3]); +%! assert (h, [1 2 3;1 2 4; 1 3 4; 2 3 4]); +%! assert (v, 8/3, 10*eps); */