Mercurial > octave

comparison scripts/geometry/delaunayn.m @ 25328:3b96348d5ccd

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
delaunayn.m: Vectorize test for 2-D zero volume simplexes (bug #53689) * delaunayn.m: Test for 2-D input case and use vectorized fast path for calculation of volume of simplexes. Otherwise, use slow for loop.
author Rik <rik@octave.org>
date Sat, 28 Apr 2018 21:08:18 -0700
parents 6652d3823428
children 00f796120a6d
comparison
equal deleted inserted replaced
25327:abc5095d58c2 25328:3b96348d5ccd
64 else 64 else
65 tol = 1e3 * eps; 65 tol = 1e3 * eps;
66 endif 66 endif
67 67
68 ## Try to remove the zero volume simplices. The volume of the i-th simplex is 68 ## Try to remove the zero volume simplices. The volume of the i-th simplex is
69 ## given by abs(det(pts(T(i,1:end-1),:)-pts(T(i,2:end),:)))/prod(1:n) 69 ## given by abs(det(pts(T(i,1:end-1),:)-pts(T(i,2:end),:)))/factorial(ndim+1)
70 ## (reference http://en.wikipedia.org/wiki/Simplex). Any simplex with a 70 ## (reference http://en.wikipedia.org/wiki/Simplex). Any simplex with a
71 ## relative volume less than some arbitrary criteria is rejected. The 71 ## relative volume less than some arbitrary criteria is rejected. The
72 ## criteria we use is the volume of the simplex corresponding to an 72 ## criteria we use is the volume of the simplex corresponding to an
73 ## orthogonal simplex is equal edge length all equal to the edge length of 73 ## orthogonal simplex is equal edge length all equal to the edge length of
74 ## the original simplex. If the relative volume is 1e3*eps then the simplex 74 ## the original simplex. If the relative volume is 1e3*eps then the simplex
75 ## is rejected. Note division of the two volumes means that the factor 75 ## is rejected. Note division of the two volumes means that the factor
76 ## prod(1:n) is dropped. 76 ## factorial(ndim+1) is dropped.
77 idx = []; 77 [nt, nd] = size (T);
78 [nt, n] = size (T); 78 if (nd == 3)
79 ## FIXME: Vectorize this for loop or convert delaunayn to .oct function 79 ## 2-D case
80 for i = 1:nt 80 np = rows (pts);
81 X = pts(T(i,1:end-1),:) - pts(T(i,2:end),:); 81 ptsz = [pts, zeros(np, 1)];
82 if (abs (det (X)) / sqrt (sumsq (X, 2)) < tol) 82 p1 = ptsz(T(:,1), :);
83 idx(end+1) = i; 83 p2 = ptsz(T(:,2), :);
84 endif 84 p3 = ptsz(T(:,3), :);
85 endfor 85 p12 = p1 - p2;
86 p23 = p2 - p3;
87 det = cross (p12, p23, 2);
88 idx = abs (det(:,3) ./ sqrt (sumsq (p12, 2))) < tol & ...
89 abs (det(:,3) ./ sqrt (sumsq (p23, 2))) < tol;
90 else
91 ## FIXME: Vectorize this for loop or convert delaunayn to .oct function
92 idx = [];
93 for i = 1:nt
94 X = pts(T(i,1:end-1),:) - pts(T(i,2:end),:);
95 if (abs (det (X)) / sqrt (sumsq (X, 2)) < tol)
96 idx(end+1) = i;
97 endif
98 endfor
99 endif
86 T(idx,:) = []; 100 T(idx,:) = [];
87 101
88 endfunction 102 endfunction
89 103
90 104