diff -r b7db401e1a99 scripts/geometry/delaunayn.m --- a/scripts/geometry/delaunayn.m Tue Jun 19 09:18:44 2018 -0700 +++ b/scripts/geometry/delaunayn.m Wed Jun 20 09:01:25 2018 +0200 @@ -75,6 +75,7 @@ ## is rejected. Note division of the two volumes means that the factor ## factorial(ndim+1) is dropped. [nt, nd] = size (T); + legacy = 0; if (nd == 3) ## 2-D case np = rows (pts); @@ -87,7 +88,7 @@ det = cross (p12, p23, 2); idx = abs (det(:,3) ./ sqrt (sumsq (p12, 2))) < tol & ... abs (det(:,3) ./ sqrt (sumsq (p23, 2))) < tol; - else + elseif legacy == 1 ## FIXME: Vectorize this for loop or convert delaunayn to .oct function idx = []; for i = 1:nt @@ -96,6 +97,35 @@ idx(end+1) = i; endif endfor + else + dim = nd - 1; + % collect all nt matrices X below each other in Xall + % Xall has size(nt*dim, dim) + Xall = pts(T(:,1:end-1),:) - pts(T(:,2:end),:); + Xall = reshape(Xall, nt, dim, dim); + Xall = reshape(permute(Xall, [2, 1, 3]), [], dim); + + % Arrange the nt matrices X in an nt*nt sparse block-diagonal + % matrix, each block size(dim,dim) + i = kron(ones(dim, 1), (1 : nt*dim)'); + [~, km] = ndgrid(1:dim, 1:dim); + j1 = kron(km, ones(nt, 1)); + j2 = kron((0 : nt-1)', ones(dim) * dim); + j = j1 + j2; + block_diag_mat = sparse(i(:), j(:), Xall(:), nt*dim, nt*dim); + + % sparse LU factorization of block_diag_mat + [L, U, P, Q] = lu(block_diag_mat); + + % the determinants are the products of the entries on the + % diagonal of the LU factors + % NOTE: dr_L guaranteed to have only "1" on the diagonal? + dr_L = full(reshape(diag(L), dim, [])); + dr_U = full(reshape(diag(U), dim, [])); + detX = prod(dr_L, 1) .* prod(dr_U, 1); + % veclen: size(nt,dim) + veclen = reshape(sqrt(sumsq(Xall, 2)), dim, nt)'; + idx = all(abs(detX(:)) ./veclen < tol, 2); endif T(idx,:) = [];