# HG changeset patch # User Nicholas R. Jankowski # Date 1664390130 14400 # Node ID bf8f33249e8666b9110585b1232cccd0d29aa28d # Parent dd6b37f67db2939f39cd99717a69b64e917cbc7e delaunayn simplex check consistency and performance improvement (bug #60818) * delaunayn.m: Apply consistent volume calculation across all trivial simplex removal code paths. Vectorize 3D simplex removal code path and minimize function calls within >3D loop for performance improvement. Update FIXME note for future performance improvement. Add input type validation checks. Add BISTs for dimensions other than 2D, simplex removal, and input validation. * etc/News.8.md: Describe function improvements under General Improvements. diff -r dd6b37f67db2 -r bf8f33249e86 etc/NEWS.8.md --- a/etc/NEWS.8.md Sun Sep 25 06:22:25 2022 -0400 +++ b/etc/NEWS.8.md Wed Sep 28 14:35:30 2022 -0400 @@ -22,6 +22,11 @@ - `quadgk` can now accept the `ArrayValued` input parameter to handle array-valued input functions. +- `delaunayn` now has consistent trivial simplex checking and removal for all +input dimensions, simplex checking 3D inputs is now vectorized, and >3D simplex +checking performance has been improved. Input type checking has also been +added for improved error handling. + ### Graphical User Interface diff -r dd6b37f67db2 -r bf8f33249e86 scripts/geometry/delaunayn.m --- a/scripts/geometry/delaunayn.m Sun Sep 25 06:22:25 2022 -0400 +++ b/scripts/geometry/delaunayn.m Wed Sep 28 14:35:30 2022 -0400 @@ -68,6 +68,12 @@ print_usage (); endif + ## NOTE: varargin options input validation is performed in __delaunayn__ + if ((! isnumeric (pts)) || (ndims (pts) > 2)) + error ("delaunayn: input points must be a two dimensional numeric array."); + endif + + ## Perform delaunay calculation using either default or specified options if (isempty (varargin) || isempty (varargin{1})) try T = __delaunayn__ (pts); @@ -82,49 +88,88 @@ T = __delaunayn__ (pts, varargin{:}); endif - if (isa (pts, "single")) - tol = 1e3 * eps ("single"); - else - tol = 1e3 * eps; - endif + ## Begin check for and removal of trivial simplices + if (! isequal (T, 0)) # skip trivial simplex check if no simplexes + + if (isa (pts, "single")) + tol = 1e3 * eps ("single"); + else + tol = 1e3 * eps; + endif + + ## Try to remove the ~zero volume simplices. The volume of the i-th simplex + ## is given by abs(det(pts(T(i,2:end),:)-pts(T(i,1),:)))/factorial(ndim+1) + ## (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 a simplex corresponding to an + ## orthogonal simplex (rectangle, rectangular prism, etc.) with edge lengths + ## equal to the common-origin edge lengths 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 factorial(ndim+1) is dropped + ## from volume calculations. + + [nt, nd] = size (T); # nt = simplex count, nd = # of simplex points + dim = nd - 1; + + ## calculate common origin edge vectors for each simplex (p2-p1,p3-p1,...) + ## store in 3D array such that: + ## rows = nt simplexes, cols = coordinates, pages = simplex edges + edge_vecs = permute (reshape (pts(T(:, 2:nd), :).', [dim, nt, dim]), ... + [2, 1, 3]) - pts(T(:, 1), :, ones (1, 1, dim)); + + ## calculate orthogonal simplex volumes for comparison + orthog_simplex_vols = sqrt (prod (sumsq (edge_vecs, 2), 3)); - ## 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),:)))/factorial(ndim+1) - ## (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 - ## factorial(ndim+1) is dropped. - [nt, nd] = size (T); - if (nd == 3) - ## 2-D case - np = rows (pts); - ptsz = [pts, zeros(np, 1)]; - p1 = ptsz(T(:,1), :); - p2 = ptsz(T(:,2), :); - p3 = ptsz(T(:,3), :); - p12 = p1 - p2; - p23 = p2 - p3; - det = cross (p12, p23, 2); - idx = abs (det (:,3) ./ sqrt (sumsq (p12, 2))) < tol & ... - abs (det (:,3) ./ sqrt (sumsq (p23, 2))) < tol; - else - ## FIXME: Vectorize this for loop or convert delaunayn to .oct function - idx = []; - for i = 1:nt - X = pts(T(i,1:end-1),:) - pts(T(i,2:end),:); - if (abs (det (X)) / sqrt (sumsq (X, 2)) < tol) - idx(end+1) = i; - endif - endfor + ## calculate simplex volumes according to problem dimension + if (nd == 3) + ## 2-D: area = cross product of triangle edge vectors + vol = edge_vecs(:,1,1) .* edge_vecs(:,2,2)... + - edge_vecs(:,1,2) .* edge_vecs(:,2,1); + + elseif (nd == 4) + ## 3-D: vol = scalar triple product [a.(b x c)] + vol = edge_vecs(:,1,1) .* ... + (edge_vecs(:,2,2) .* edge_vecs(:,3,3) - ... + edge_vecs(:,3,2) .* edge_vecs(:,2,3)) ... + - edge_vecs(:,2,1) .* ... + (edge_vecs(:,1,2) .* edge_vecs(:,3,3) - ... + edge_vecs(:,3,2) .* edge_vecs(:,1,3)) ... + + edge_vecs(:,3,1) .* ... + (edge_vecs(:,1,2) .* edge_vecs(:,2,3) - ... + edge_vecs(:,2,2) .* edge_vecs(:,1,3)); + + else + ## >= 4-D: simplex 'volume' proportional to det|edge_vecs| + + ## FIXME: Vectorize this for nD inputs without excessive memory impact + ## over __delaunayn__ itself. Perhaps with a paged determinant function. + ## see bug #60818 for speed/memory improvement attempts and concerns. + + vol = zeros (nt, 1); + + ##reshape so det can operate in dim1&2 + edge_vecs = permute (edge_vecs, [3, 2, 1]); + + ## calculate determinant for arbitrary problem dimension + for ii = 1:nt + vol(ii) = det (edge_vecs(:, :, ii)); + endfor + endif + + ## mark simplices with relative volume < tol for removal + idx = (abs ((vol) ./ orthog_simplex_vols)) < tol; + + ##Remove trivially small simplexes from T + T(idx, :) = []; endif - T(idx,:) = []; - endfunction +## Test 1-D input +%!testif HAVE_QHULL +%! assert (sortrows (sort (delaunayn ([1;2]), 2)), [1, 2]); +%! assert (sortrows (sort (delaunayn ([1;2;3]), 2)), [1, 2; 2, 3]); +## Test 2-D input %!testif HAVE_QHULL %! x = [-1, 0; 0, 1; 1, 0; 0, -1; 0, 0]; %! assert (sortrows (sort (delaunayn (x), 2)), [1,2,5;1,4,5;2,3,5;3,4,5]); @@ -135,7 +180,29 @@ %! assert (sortrows (sort (delaunayn ([x(:) y(:) z(:)]), 2)), %! [1,2,3,4;1,2,4,5]); -## FIXME: Need tests for delaunayn +## 3D test with trivial simplex removal +%!testif HAVE_QHULL +%! x = [0 0 0; 0 0 1; 0 1 0; 1 0 0; 0 1 1; 1 0 1; 1 1 0; 1 1 1; 0.5 0.5 0.5]; +%! T = sortrows (sort (delaunayn (x), 2)); +%! assert (rows (T), 12); + +## 4D single simplex test +%!testif HAVE_QHULL +%! x = [0 0 0 0; 1 0 0 0; 1 1 0 0; 0 0 1 0; 0 0 0 1]; +%! T = sort (delaunayn (x), 2); +%! assert (T, [1 2 3 4 5]); + +## 4D two simplices test +%!testif HAVE_QHULL +%! x = [0 0 0 0; 1 0 0 0; 1 1 0 0; 0 0 1 0; 0 0 0 1; 0 0 0 2]; +%! T = sortrows (sort (delaunayn (x), 2)); +%! assert (rows (T), 2); +%! assert (T, [1 2 3 4 5; 2 3 4 5 6]); ## Input validation tests %!error delaunayn () +%!error delaunayn ("abc") +%!error delaunayn ({1}) +%!error delaunayn (true) +%!error delaunayn (cat (3, [1 2 3], [4 5 6])) +