Mercurial > octave

changeset 25328:3b96348d5ccd

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
delaunayn.m: Vectorize test for 2-D zero volume simplexes (bug #53689) * delaunayn.m: Test for 2-D input case and use vectorized fast path for calculation of volume of simplexes. Otherwise, use slow for loop.
author Rik <rik@octave.org>
date Sat, 28 Apr 2018 21:08:18 -0700
parents abc5095d58c2
children 55c1b0e5b07a
files scripts/geometry/delaunayn.m
diffstat 1 files changed, 25 insertions(+), 11 deletions(-) [+]
line wrap: on
line diff
--- a/scripts/geometry/delaunayn.m	Fri Apr 27 14:57:16 2018 -0400
+++ b/scripts/geometry/delaunayn.m	Sat Apr 28 21:08:18 2018 -0700
@@ -66,23 +66,37 @@
   endif
 
   ## 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)
+  ## 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
-  ## prod(1:n) is dropped.
-  idx = [];
-  [nt, n] = size (T);
-  ## FIXME: Vectorize this for loop or convert delaunayn to .oct function
-  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
+  ## 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
+  endif
   T(idx,:) = [];
 
 endfunction