# HG changeset patch
# User Nicholas R. Jankowski <jankowski.nicholas@gmail.com>
# Date 1664507345 14400
# Node ID 8b75954a467046412fc0d0ad971bd1b09dbababd
# Parent  3dae836c598c8b91dc9604aa470788cc11fa3853
delaunayn: adjust node ordering for positive outward normal vectors (bug #53397)

* delaunayn.m: Check sign of simplex volume, flip node order for negative
volumes to ensure positive (outward-pointing) normal vectors. Add BISTs to
check for positive volumes.
* etc/News.8.md: Append function improvement note to delaunayn change
paragraph under General Improvements.

diff -r 3dae836c598c -r 8b75954a4670 etc/NEWS.8.md
--- a/etc/NEWS.8.md	Fri Sep 30 06:38:59 2022 -0400
+++ b/etc/NEWS.8.md	Thu Sep 29 23:09:05 2022 -0400
@@ -24,8 +24,9 @@
 
 - `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.
+checking performance has been improved.  Simplexes points are now ordered so
+they will all have positive outward normal vectors. Input type checking has
+also been added for improved error handling.
 
 ### Graphical User Interface
 
diff -r 3dae836c598c -r 8b75954a4670 scripts/geometry/delaunayn.m
--- a/scripts/geometry/delaunayn.m	Fri Sep 30 06:38:59 2022 -0400
+++ b/scripts/geometry/delaunayn.m	Thu Sep 29 23:09:05 2022 -0400
@@ -139,10 +139,11 @@
                 edge_vecs(:,2,2) .* edge_vecs(:,1,3));
 
     else
-      ## >= 4-D: simplex 'volume' proportional to det|edge_vecs|
+      ## 1D and >= 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.
+      ## over __delaunayn__ itself, or move simplex checking into __delaunayn__.
+      ## Perhaps with an optimized page-wise determinant.
       ## see bug #60818 for speed/memory improvement attempts and concerns.
 
       vol = zeros (nt, 1);
@@ -161,6 +162,16 @@
 
     ##Remove trivially small simplexes from T
     T(idx, :) = [];
+
+    ## ensure ccw node order for consistent outward normal (bug #53397)
+    ## simplest method of maintaining positive unit normal direction is to
+    ## reverse order of two nodes. this preserves 'nice' monotonic descending
+    ## node 1 ordering.  Currently ignores 1D cases for compatibility.
+
+    if ((dim > 1) && any (negvol = (vol (!idx) < 0)))
+      T(negvol, [2, 3]) = T(negvol, [3, 2]);
+    endif
+
   endif
 endfunction
 
@@ -180,25 +191,49 @@
 %! assert (sortrows (sort (delaunayn ([x(:) y(:) z(:)]), 2)),
 %!         [1,2,3,4;1,2,4,5]);
 
-## 3D test with trivial simplex removal
+## 3-D 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
+## 4-D 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
+## 4-D 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]);
 
+## Test negative simplex produce positive normals
+## 2-D test
+%!testif HAVE_QHULL <*53397>
+%! x = [-1, 0; 0, 1; 1, 0; 0, -1; 0, 0];
+%! y = delaunayn (x);
+%! edges = permute (reshape (x(y(:, 2:end), :).', [2, 4, 2]), [2, 1, 3]) - ...
+%!                   x(y(:, 1), :, ones (1, 1, 2));
+%! vol = edges(:,1,1) .* edges(:,2,2) - edges(:,1,2) .* edges(:,2,1);
+%! assert (all (vol >= 0));
+
+## 3-D test
+%!testif HAVE_QHULL <*53397>
+%! x = [[-1, -1, 1, 0, -1]',[-1, 1, 1, 0, -1]',[0, 0, 0, 1, 1]'];
+%! y = delaunayn (x);
+%! edges = permute (reshape (x(y(:, 2:end), :).', [3, 2, 3]), [2, 1, 3]) - ...
+%!                   x(y(:, 1), :, ones (1, 1, 3));
+%! vol = edges(:,1,1) .* ...
+%!            (edges(:,2,2) .* edges(:,3,3) - edges(:,3,2) .* edges(:,2,3)) ...
+%!       - edges(:,2,1) .* ...
+%!            (edges(:,1,2) .* edges(:,3,3) - edges(:,3,2) .* edges(:,1,3)) ...
+%!       + edges(:,3,1) .* ...
+%!            (edges(:,1,2) .* edges(:,2,3) - edges(:,2,2) .* edges(:,1,3));
+%! assert (all (vol >= 0));
+
 ## Input validation tests
 %!error <Invalid call> delaunayn ()
 %!error <input points must be> delaunayn ("abc")