Mercurial > octave

comparison scripts/geometry/delaunayn.m @ 31248:8b75954a4670

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
delaunayn: adjust node ordering for positive outward normal vectors (bug #53397) * delaunayn.m: Check sign of simplex volume, flip node order for negative volumes to ensure positive (outward-pointing) normal vectors. Add BISTs to check for positive volumes. * etc/News.8.md: Append function improvement note to delaunayn change paragraph under General Improvements.
author Nicholas R. Jankowski <jankowski.nicholas@gmail.com>
date Thu, 29 Sep 2022 23:09:05 -0400
parents bf8f33249e86
children a40c0b7aa376
comparison
equal deleted inserted replaced
31247:3dae836c598c 31248:8b75954a4670
137 + edge_vecs(:,3,1) .* ... 137 + edge_vecs(:,3,1) .* ...
138 (edge_vecs(:,1,2) .* edge_vecs(:,2,3) - ... 138 (edge_vecs(:,1,2) .* edge_vecs(:,2,3) - ...
139 edge_vecs(:,2,2) .* edge_vecs(:,1,3)); 139 edge_vecs(:,2,2) .* edge_vecs(:,1,3));
140 140
141 else 141 else
142 ## >= 4-D: simplex 'volume' proportional to det|edge_vecs| 142 ## 1D and >= 4-D: simplex 'volume' proportional to det|edge_vecs|
143 143
144 ## FIXME: Vectorize this for nD inputs without excessive memory impact 144 ## FIXME: Vectorize this for nD inputs without excessive memory impact
145 ## over __delaunayn__ itself. Perhaps with a paged determinant function. 145 ## over __delaunayn__ itself, or move simplex checking into __delaunayn__.
146 ## Perhaps with an optimized page-wise determinant.
146 ## see bug #60818 for speed/memory improvement attempts and concerns. 147 ## see bug #60818 for speed/memory improvement attempts and concerns.
147 148
148 vol = zeros (nt, 1); 149 vol = zeros (nt, 1);
149 150
150 ##reshape so det can operate in dim1&2 151 ##reshape so det can operate in dim1&2
159 ## mark simplices with relative volume < tol for removal 160 ## mark simplices with relative volume < tol for removal
160 idx = (abs ((vol) ./ orthog_simplex_vols)) < tol; 161 idx = (abs ((vol) ./ orthog_simplex_vols)) < tol;
161 162
162 ##Remove trivially small simplexes from T 163 ##Remove trivially small simplexes from T
163 T(idx, :) = []; 164 T(idx, :) = [];
165
166 ## ensure ccw node order for consistent outward normal (bug #53397)
167 ## simplest method of maintaining positive unit normal direction is to
168 ## reverse order of two nodes. this preserves 'nice' monotonic descending
169 ## node 1 ordering. Currently ignores 1D cases for compatibility.
170
171 if ((dim > 1) && any (negvol = (vol (!idx) < 0)))
172 T(negvol, [2, 3]) = T(negvol, [3, 2]);
173 endif
174
164 endif 175 endif
165 endfunction 176 endfunction
166 177
167 ## Test 1-D input 178 ## Test 1-D input
168 %!testif HAVE_QHULL 179 %!testif HAVE_QHULL
178 %!testif HAVE_QHULL 189 %!testif HAVE_QHULL
179 %! x = [-1, -1, 1, 0, -1]; y = [-1, 1, 1, 0, -1]; z = [0, 0, 0, 1, 1]; 190 %! x = [-1, -1, 1, 0, -1]; y = [-1, 1, 1, 0, -1]; z = [0, 0, 0, 1, 1];
180 %! assert (sortrows (sort (delaunayn ([x(:) y(:) z(:)]), 2)), 191 %! assert (sortrows (sort (delaunayn ([x(:) y(:) z(:)]), 2)),
181 %! [1,2,3,4;1,2,4,5]); 192 %! [1,2,3,4;1,2,4,5]);
182 193
183 ## 3D test with trivial simplex removal 194 ## 3-D test with trivial simplex removal
184 %!testif HAVE_QHULL 195 %!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]; 196 %! 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)); 197 %! T = sortrows (sort (delaunayn (x), 2));
187 %! assert (rows (T), 12); 198 %! assert (rows (T), 12);
188 199
189 ## 4D single simplex test 200 ## 4-D single simplex test
190 %!testif HAVE_QHULL 201 %!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]; 202 %! 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); 203 %! T = sort (delaunayn (x), 2);
193 %! assert (T, [1 2 3 4 5]); 204 %! assert (T, [1 2 3 4 5]);
194 205
195 ## 4D two simplices test 206 ## 4-D two simplices test
196 %!testif HAVE_QHULL 207 %!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]; 208 %! 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)); 209 %! T = sortrows (sort (delaunayn (x), 2));
199 %! assert (rows (T), 2); 210 %! assert (rows (T), 2);
200 %! assert (T, [1 2 3 4 5; 2 3 4 5 6]); 211 %! assert (T, [1 2 3 4 5; 2 3 4 5 6]);
212
213 ## Test negative simplex produce positive normals
214 ## 2-D test
215 %!testif HAVE_QHULL <*53397>
216 %! x = [-1, 0; 0, 1; 1, 0; 0, -1; 0, 0];
217 %! y = delaunayn (x);
218 %! edges = permute (reshape (x(y(:, 2:end), :).', [2, 4, 2]), [2, 1, 3]) - ...
219 %! x(y(:, 1), :, ones (1, 1, 2));
220 %! vol = edges(:,1,1) .* edges(:,2,2) - edges(:,1,2) .* edges(:,2,1);
221 %! assert (all (vol >= 0));
222
223 ## 3-D test
224 %!testif HAVE_QHULL <*53397>
225 %! x = [[-1, -1, 1, 0, -1]',[-1, 1, 1, 0, -1]',[0, 0, 0, 1, 1]'];
226 %! y = delaunayn (x);
227 %! edges = permute (reshape (x(y(:, 2:end), :).', [3, 2, 3]), [2, 1, 3]) - ...
228 %! x(y(:, 1), :, ones (1, 1, 3));
229 %! vol = edges(:,1,1) .* ...
230 %! (edges(:,2,2) .* edges(:,3,3) - edges(:,3,2) .* edges(:,2,3)) ...
231 %! - edges(:,2,1) .* ...
232 %! (edges(:,1,2) .* edges(:,3,3) - edges(:,3,2) .* edges(:,1,3)) ...
233 %! + edges(:,3,1) .* ...
234 %! (edges(:,1,2) .* edges(:,2,3) - edges(:,2,2) .* edges(:,1,3));
235 %! assert (all (vol >= 0));
201 236
202 ## Input validation tests 237 ## Input validation tests
203 %!error <Invalid call> delaunayn () 238 %!error <Invalid call> delaunayn ()
204 %!error <input points must be> delaunayn ("abc") 239 %!error <input points must be> delaunayn ("abc")
205 %!error <input points must be> delaunayn ({1}) 240 %!error <input points must be> delaunayn ({1})