noting arb's comment #5, this is now actually fairly trivial to implement since now that the bug #60818 patch made simplex checking for all dimensionalities consistently calculate the determinant(volume). can conditionally flip nodes 2 and 3 to change sign of determinant.  checking against matlab, for many cases where delaunay produces the same simplexes, the ordering is now identical.  the only one that seems to ignore this is the trivial 1D case where matlab allows negative oriented vectors in the output.

implemented with addition of the following to the end of the function:

if ((dim > 1) && any (negvol = (vol (!idx) < 0)))
  T(negvol, [2, 3]) = T(negvol, [3, 2]);
endif

(could add a dim==1 path in there where it switches 1 and 2, but i don't think it's necessary. easy enough to add later if anyone objects.)

added a couple BISTs for 2d and 3d.  updated the NEWS paragraph about delaunayn improvement with this bit.  pushed to default as
https://hg.savannah.gnu.org/hgweb/octave/rev/8b75954a4670
marking as Ready for Test.

The QHull FAQ has this question and answer:

When I compute a plane equation from a facet, I sometimes get an outward-pointing normal and sometimes an inward-pointing normal

    Qhull orients simplicial facets, and prints oriented output for 'i', 'Ft', and other options. The orientation depends on both the vertex order and the flag facet->toporient.

    Qhull does not orient non-simplicial facets. Instead it orients the facet's ridges. These are printed with the 'Qt' and 'Ft' option. The facet's hyperplane is oriented.

So, some combination of options to QHull + post-processing is probably required.

https://www.mathworks.com/matlabcentral/answers/321387-what-is-matlab-convention-for-the-direction-of-normal-vector-of-a-triangulated-mesh#answer_252306

In Octave the delaunay/delaunayn function calls returns inconsistent (both clockwise and counter-clockwise) orientations in 2d, while Matlab seems to return only counter-clockwise oriented triangulations:

p = [0 0;1 0;1 1;0 1;.5 .5]
t = delaunay( p );

d12 = p(t(:,2),:) - p(t(:,1),:);
d13 = p(t(:,3),:) - p(t(:,1),:);
a = ( d12(:,1).*d13(:,2) - d12(:,2).*d13(:,1) )/2;

is_clockwise = a<0

This can easily be fixed in post triangulation with

t(is_clockwise,2:3) =  t(is_clockwise,[3,2]);

but it could be good to make it consistent and built-in by default.