Mercurial > octave

comparison scripts/geometry/delaunayn.m @ 31240:bf8f33249e86

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
delaunayn simplex check consistency and performance improvement (bug #60818) * delaunayn.m: Apply consistent volume calculation across all trivial simplex removal code paths. Vectorize 3D simplex removal code path and minimize function calls within >3D loop for performance improvement. Update FIXME note for future performance improvement. Add input type validation checks. Add BISTs for dimensions other than 2D, simplex removal, and input validation. * etc/News.8.md: Describe function improvements under General Improvements.
author Nicholas R. Jankowski <jankowski.nicholas@gmail.com>
date Wed, 28 Sep 2022 14:35:30 -0400
parents 66456820ff59
children 8b75954a4670
comparison
equal deleted inserted replaced
31239:dd6b37f67db2 31240:bf8f33249e86
66 66
67 if (nargin < 1) 67 if (nargin < 1)
68 print_usage (); 68 print_usage ();
69 endif 69 endif
70 70
71 ## NOTE: varargin options input validation is performed in __delaunayn__
72 if ((! isnumeric (pts)) || (ndims (pts) > 2))
73 error ("delaunayn: input points must be a two dimensional numeric array.");
74 endif
75
76 ## Perform delaunay calculation using either default or specified options
71 if (isempty (varargin) || isempty (varargin{1})) 77 if (isempty (varargin) || isempty (varargin{1}))
72 try 78 try
73 T = __delaunayn__ (pts); 79 T = __delaunayn__ (pts);
74 catch err 80 catch err
75 if (columns (pts) <= 2) 81 if (columns (pts) <= 2)
80 end_try_catch 86 end_try_catch
81 else 87 else
82 T = __delaunayn__ (pts, varargin{:}); 88 T = __delaunayn__ (pts, varargin{:});
83 endif 89 endif
84 90
85 if (isa (pts, "single")) 91 ## Begin check for and removal of trivial simplices
86 tol = 1e3 * eps ("single"); 92 if (! isequal (T, 0)) # skip trivial simplex check if no simplexes
87 else 93
88 tol = 1e3 * eps; 94 if (isa (pts, "single"))
89 endif 95 tol = 1e3 * eps ("single");
90 96 else
91 ## Try to remove the zero volume simplices. The volume of the i-th simplex is 97 tol = 1e3 * eps;
92 ## given by abs(det(pts(T(i,1:end-1),:)-pts(T(i,2:end),:)))/factorial(ndim+1) 98 endif
93 ## (reference http://en.wikipedia.org/wiki/Simplex). Any simplex with a 99
94 ## relative volume less than some arbitrary criteria is rejected. The 100 ## Try to remove the ~zero volume simplices. The volume of the i-th simplex
95 ## criteria we use is the volume of the simplex corresponding to an 101 ## is given by abs(det(pts(T(i,2:end),:)-pts(T(i,1),:)))/factorial(ndim+1)
96 ## orthogonal simplex is equal edge length all equal to the edge length of 102 ## (reference http://en.wikipedia.org/wiki/Simplex). Any simplex with a
97 ## the original simplex. If the relative volume is 1e3*eps then the simplex 103 ## relative volume less than some arbitrary criteria is rejected. The
98 ## is rejected. Note division of the two volumes means that the factor 104 ## criteria we use is the volume of a simplex corresponding to an
99 ## factorial(ndim+1) is dropped. 105 ## orthogonal simplex (rectangle, rectangular prism, etc.) with edge lengths
100 [nt, nd] = size (T); 106 ## equal to the common-origin edge lengths of the original simplex. If the
101 if (nd == 3) 107 ## relative volume is 1e3*eps then the simplex is rejected. Note division
102 ## 2-D case 108 ## of the two volumes means that the factor factorial(ndim+1) is dropped
103 np = rows (pts); 109 ## from volume calculations.
104 ptsz = [pts, zeros(np, 1)]; 110
105 p1 = ptsz(T(:,1), :); 111 [nt, nd] = size (T); # nt = simplex count, nd = # of simplex points
106 p2 = ptsz(T(:,2), :); 112 dim = nd - 1;
107 p3 = ptsz(T(:,3), :); 113
108 p12 = p1 - p2; 114 ## calculate common origin edge vectors for each simplex (p2-p1,p3-p1,...)
109 p23 = p2 - p3; 115 ## store in 3D array such that:
110 det = cross (p12, p23, 2); 116 ## rows = nt simplexes, cols = coordinates, pages = simplex edges
111 idx = abs (det (:,3) ./ sqrt (sumsq (p12, 2))) < tol & ... 117 edge_vecs = permute (reshape (pts(T(:, 2:nd), :).', [dim, nt, dim]), ...
112 abs (det (:,3) ./ sqrt (sumsq (p23, 2))) < tol; 118 [2, 1, 3]) - pts(T(:, 1), :, ones (1, 1, dim));
113 else 119
114 ## FIXME: Vectorize this for loop or convert delaunayn to .oct function 120 ## calculate orthogonal simplex volumes for comparison
115 idx = []; 121 orthog_simplex_vols = sqrt (prod (sumsq (edge_vecs, 2), 3));
116 for i = 1:nt 122
117 X = pts(T(i,1:end-1),:) - pts(T(i,2:end),:); 123 ## calculate simplex volumes according to problem dimension
118 if (abs (det (X)) / sqrt (sumsq (X, 2)) < tol) 124 if (nd == 3)
119 idx(end+1) = i; 125 ## 2-D: area = cross product of triangle edge vectors
120 endif 126 vol = edge_vecs(:,1,1) .* edge_vecs(:,2,2)...
121 endfor 127 - edge_vecs(:,1,2) .* edge_vecs(:,2,1);
122 endif 128
123 T(idx,:) = []; 129 elseif (nd == 4)
124 130 ## 3-D: vol = scalar triple product [a.(b x c)]
131 vol = edge_vecs(:,1,1) .* ...
132 (edge_vecs(:,2,2) .* edge_vecs(:,3,3) - ...
133 edge_vecs(:,3,2) .* edge_vecs(:,2,3)) ...
134 - edge_vecs(:,2,1) .* ...
135 (edge_vecs(:,1,2) .* edge_vecs(:,3,3) - ...
136 edge_vecs(:,3,2) .* edge_vecs(:,1,3)) ...
137 + edge_vecs(:,3,1) .* ...
138 (edge_vecs(:,1,2) .* edge_vecs(:,2,3) - ...
139 edge_vecs(:,2,2) .* edge_vecs(:,1,3));
140
141 else
142 ## >= 4-D: simplex 'volume' proportional to det|edge_vecs|
143
144 ## FIXME: Vectorize this for nD inputs without excessive memory impact
145 ## over __delaunayn__ itself. Perhaps with a paged determinant function.
146 ## see bug #60818 for speed/memory improvement attempts and concerns.
147
148 vol = zeros (nt, 1);
149
150 ##reshape so det can operate in dim1&2
151 edge_vecs = permute (edge_vecs, [3, 2, 1]);
152
153 ## calculate determinant for arbitrary problem dimension
154 for ii = 1:nt
155 vol(ii) = det (edge_vecs(:, :, ii));
156 endfor
157 endif
158
159 ## mark simplices with relative volume < tol for removal
160 idx = (abs ((vol) ./ orthog_simplex_vols)) < tol;
161
162 ##Remove trivially small simplexes from T
163 T(idx, :) = [];
164 endif
125 endfunction 165 endfunction
126 166
127 167 ## Test 1-D input
168 %!testif HAVE_QHULL
169 %! assert (sortrows (sort (delaunayn ([1;2]), 2)), [1, 2]);
170 %! assert (sortrows (sort (delaunayn ([1;2;3]), 2)), [1, 2; 2, 3]);
171
172 ## Test 2-D input
128 %!testif HAVE_QHULL 173 %!testif HAVE_QHULL
129 %! x = [-1, 0; 0, 1; 1, 0; 0, -1; 0, 0]; 174 %! x = [-1, 0; 0, 1; 1, 0; 0, -1; 0, 0];
130 %! assert (sortrows (sort (delaunayn (x), 2)), [1,2,5;1,4,5;2,3,5;3,4,5]); 175 %! assert (sortrows (sort (delaunayn (x), 2)), [1,2,5;1,4,5;2,3,5;3,4,5]);
131 176
132 ## Test 3-D input 177 ## Test 3-D input
133 %!testif HAVE_QHULL 178 %!testif HAVE_QHULL
134 %! x = [-1, -1, 1, 0, -1]; y = [-1, 1, 1, 0, -1]; z = [0, 0, 0, 1, 1]; 179 %! x = [-1, -1, 1, 0, -1]; y = [-1, 1, 1, 0, -1]; z = [0, 0, 0, 1, 1];
135 %! assert (sortrows (sort (delaunayn ([x(:) y(:) z(:)]), 2)), 180 %! assert (sortrows (sort (delaunayn ([x(:) y(:) z(:)]), 2)),
136 %! [1,2,3,4;1,2,4,5]); 181 %! [1,2,3,4;1,2,4,5]);
137 182
138 ## FIXME: Need tests for delaunayn 183 ## 3D test with trivial simplex removal
184 %!testif HAVE_QHULL
185 %! x = [0 0 0; 0 0 1; 0 1 0; 1 0 0; 0 1 1; 1 0 1; 1 1 0; 1 1 1; 0.5 0.5 0.5];
186 %! T = sortrows (sort (delaunayn (x), 2));
187 %! assert (rows (T), 12);
188
189 ## 4D single simplex test
190 %!testif HAVE_QHULL
191 %! x = [0 0 0 0; 1 0 0 0; 1 1 0 0; 0 0 1 0; 0 0 0 1];
192 %! T = sort (delaunayn (x), 2);
193 %! assert (T, [1 2 3 4 5]);
194
195 ## 4D two simplices test
196 %!testif HAVE_QHULL
197 %! x = [0 0 0 0; 1 0 0 0; 1 1 0 0; 0 0 1 0; 0 0 0 1; 0 0 0 2];
198 %! T = sortrows (sort (delaunayn (x), 2));
199 %! assert (rows (T), 2);
200 %! assert (T, [1 2 3 4 5; 2 3 4 5 6]);
139 201
140 ## Input validation tests 202 ## Input validation tests
141 %!error <Invalid call> delaunayn () 203 %!error <Invalid call> delaunayn ()
204 %!error <input points must be> delaunayn ("abc")
205 %!error <input points must be> delaunayn ({1})
206 %!error <input points must be> delaunayn (true)
207 %!error <input points must be> delaunayn (cat (3, [1 2 3], [4 5 6]))
208