diff -r 4adeabc1bbfe scripts/geometry/delaunayn.m --- a/scripts/geometry/delaunayn.m Tue Apr 17 21:40:31 2018 -0400 +++ b/scripts/geometry/delaunayn.m Wed Apr 18 09:29:00 2018 +0200 @@ -74,15 +74,33 @@ ## 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 + dim = size(pts, 2); + if dim == 2 + np = size(pts, 1); + ptsz = [pts, zeros(np, 1)]; + p1 = ptsz(T(:,1), :); + p2 = ptsz(T(:,2), :); + p3 = ptsz(T(:,3), :); + p12 = p1 - p2; + p23 = p2 - p3; + a = cross(p12, p23, 2); + idx = abs(a(:,3) ./ sqrt(sumsq([p12, p23], 2))) < tol; + else + idx = []; + for i = 1:nt + X = pts(T(i,1:end-1),:) - pts(T(i,2:end),:); + if i == 1 + size(X) + det(X) + sqrt(sumsq(X,2)) + end + if (abs (det (X)) / sqrt (sumsq (X, 2)) < tol) + idx(end+1) = i; + endif + endfor + endif T(idx,:) = []; endfunction