# HG changeset patch
# User Nicholas R. Jankowski <jankowski.nicholas@gmail.com>
# Date 1637388642 18000
#      Sat Nov 20 01:10:42 2021 -0500
# Node ID d5585e5f375164a38a206e11dca34b6641a38f6f
# Parent  b568c29bb5be7869cfa4e4fbdacaad0270111687
delaunayn.m: Overhaul trivial simplex check for consistency and speed (bug #60818)

* scripts/geometry/delaunay.m: Rewrites trivial simplex check so it uses the
same volume check for all input dimensions. Vectorizes 3D codepath and large
portion of nD for loop for perforance improvement. Rewrites 2D codepath to align
with 3D and nD approach.

diff -r b568c29bb5be -r d5585e5f3751 scripts/geometry/delaunayn.m
--- a/scripts/geometry/delaunayn.m	Fri Nov 19 20:57:09 2021 -0800
+++ b/scripts/geometry/delaunayn.m	Sat Nov 20 01:10:42 2021 -0500
@@ -80,49 +80,80 @@
     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,etc)
+    ## 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 case - area = cross product of triange edge vectors
+      vol = edge_vecs(:,1,1) .* edge_vecs(:,2,2)...
+                  - edge_vecs(:,1,2) .* edge_vecs(:,2,1);
+
+    elseif (nd == 4)
+      ## 3-D case - 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.
+
+      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
 
-
 %!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]);
@@ -132,6 +163,25 @@
 %! x = [-1, -1, 1, 0, -1]; y = [-1, 1, 1, 0, -1]; z = [0, 0, 0, 1, 1];
 %! assert (sortrows (sort (delaunayn ([x(:) y(:) z(:)]), 2)), [1,2,3,4;1,2,4,5]);
 
+## 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]);
+
 ## FIXME: Need tests for delaunayn
 
 ## Input validation tests
