diff -r ae0518976e2b scripts/geometry/delaunayn.m --- a/scripts/geometry/delaunayn.m Sat Apr 21 22:27:17 2018 +0200 +++ b/scripts/geometry/delaunayn.m Mon Apr 23 21:30:54 2018 +0200 @@ -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),:)))/prod(1:np_simplex) ## (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); + ## prod(1:np_simplex) is dropped. + [nt, np_simplex] = 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 + if np_simplex == 3 + ## 2-D case + 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; + det = cross(p12, p23, 2); + idx = abs(det(:,3) ./ sqrt(sumsq(p12, 2))) < tol & ... + abs(det(:,3) ./ sqrt(sumsq(p23, 2))) < tol; + else + 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