# HG changeset patch
# User Nicholas R. Jankowski <jankowskin@asme.org>
# Date 1628521299 14400
#      Mon Aug 09 11:01:39 2021 -0400
# Node ID a6fda424fe5c927c44a8827f94e4dafce753c54e
# Parent  133abbad39d221d2648f93f6188ee56aeb14b033
v6

diff -r 133abbad39d2 -r a6fda424fe5c scripts/geometry/delaunayn.m
--- a/scripts/geometry/delaunayn.m	Sun Aug 08 21:17:35 2021 +0200
+++ b/scripts/geometry/delaunayn.m	Mon Aug 09 11:01:39 2021 -0400
@@ -91,37 +91,107 @@
   ## (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;
+  ## orthogonal simplex with (ndim-1) edge lengths equal to the edge lengths of
+  ## the original simplex.  If the relative volume is smaller than 1e3*eps, the
+  ## simplex is rejected.  Note division of the two volumes means that the
+  ## factor factorial(ndim+1) is dropped from the 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, ...)
+  ## stored in 3D array such that
+  ## rows = nt simplexes, cols = coordinates, pages = simplex edges
+  edg_vec =  permute(reshape(pts(T(:,2:nd),:).',[dim, nt, dim]),[2 1 3]) ...
+               - pts(T(:,1), :, ones (1,1,dim));
+
+  ## Calculate simplex volumes according to dimensionality of problem
+  if  (dim <= 6)
+    ## faster vector-product code paths for up to 3D cases
+
+    if (dim == 2)
+      ## 2D Use simple component cross product to calculate 2D triangle area
+      vol = abs (edg_vec(:,1,1) .* edg_vec(:,2,2)...
+               - edg_vec(:,1,2) .* edg_vec(:,2,1));
+    else
+      ## >=3D: Use scalar triple product to calculate 3D tetrahedron volume, and
+      ## explicit laplace expansion for higher dimension determinants 
+      vol = abs (laplace_det (dim, edg_vec));
+    endif
+
   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
+    ## higher dimensions
+
+    ## Place edge vectors into a block diagonal matrix
+    eqs = sparse (dim * nt, dim * nt);
+    eqs(logical (kron (speye (nt, nt), true (dim)))) = ...
+          permute (edg_vec, [2,3,1])(:);
+    ## Extract diagonal of LU factorization of that block diagonal matrix
+    [~, eqs, p, ~] = lu (eqs, "vector");
+    R = abs (diag (eqs));
+
+    ## extract simplex volumes as product of diagonal elements of u
+    ## preserving order relative to simplexes from delaunay triangulation
+    [~, rev_sort]= sort (kron (1:nt, ones (1, dim))(p));
+    vol = prod (reshape (R(rev_sort), dim, nt), 1).';
+
   endif
+
+  ## Check for small volumes: compare simplex to volume with orthogonal edges
+  idx = (vol ./ prod (sqrt (sumsq (edg_vec, 2)), 3)) < tol;
+
+  ##Remove trivially small simplexes 
   T(idx,:) = [];
 
 endfunction
 
+function retval = laplace_det (dim, v)
+  ## using scalar triple product for 3D determinant, and for >3D use laplace
+  ## expansion recursively down to 3D.  
+  ## input vector shape: rows - simplex, column - coordinate, page - edgevectors
+  switch dim
+    case 3
+##      retval = dot (v(:,:,1), cross (v(:,:,2), v(:,:,3), 2), 2);
+      retval = v(:,1,1) .* (v(:,2,2) .* v(:,3,3) - v(:,3,2) .* v(:,2,3))...
+        - v(:,2,1) .* (v(:,1,2) .* v(:,3,3) - v(:,3,2) .* v(:,1,3))...
+        + v(:,3,1) .* (v(:,1,2) .* v(:,2,3) - v(:,2,2) .* v(:,1,3));
+
+    case 4
+      retval = v(:,1,1) .* laplace_det (3, v(:, 2:4, 2:4))...
+        - v(:,2,1) .* laplace_det (3, v(:, [1,3,4], 2:4))...
+        + v(:,3,1) .* laplace_det (3, v(:, [1,2,4], 2:4))...
+        - v(:,4,1) .* laplace_det (3, v(:, 1:3, 2:4));
+
+    case 5
+      retval = v(:,1,1) .* laplace_det (4, v(:, 2:5, 2:5))...
+        - v(:,2,1) .* laplace_det (4, v(:, [1,3:5], 2:5))...
+        + v(:,3,1) .* laplace_det (4, v(:, [1,2,4,5], 2:5))...
+        - v(:,4,1) .* laplace_det (4, v(:, [1:3,5], 2:5))...
+        + v(:,5,1) .* laplace_det (4, v(:, 1:4, 2:5));
+
+    case 6
+      retval = v(:,1,1) .* laplace_det (5, v(:,2:6, 2:6))...
+        - v(:,2,1) .* laplace_det (5, v(:, [1,3:6], 2:6))...
+        + v(:,3,1) .* laplace_det (5, v(:, [1,2,4:6], 2:6))...
+        - v(:,4,1) .* laplace_det (5, v(:, [1:3,5,6], 2:6))...
+        + v(:,5,1) .* laplace_det (5, v(:, [1:4,6], 2:6))...
+        - v(:,6,1) .* laplace_det (5, v(:, 1:5, 2:6));
+
+    case 7
+      retval = v(:,1,1) .* laplace_det (6, v(:, 2:7, 2:7))...
+        - v(:,2,1) .* laplace_det (6, v(:, [1,3:7], 2:7))...
+        + v(:,3,1) .* laplace_det (6, v(:, [1,2,4:7], 2:7))...
+        - v(:,4,1) .* laplace_det (6, v(:, [1:3,5:7], 2:7))...
+        + v(:,5,1) .* laplace_det (6, v(:, [1:4,6,7], 2:7))...
+        - v(:,6,1) .* laplace_det (6, v(:, [1:5,7], 2:7))...
+        + v(:,7,1) .* laplace_det (6, v(:, 1:6, 2:7));
+
+    otherwise
+      error ("delaunayn: DIM too high for laplace expansion deternimant");
+  endswitch
+
+endfunction
+
 
 %!testif HAVE_QHULL
 %! x = [-1, 0; 0, 1; 1, 0; 0, -1; 0, 0];
