GCC Code Coverage Report

Directory: ./
File: lib/geogram/delaunay/delaunay_3d.cpp
Date: 2026-09-07 02:25:23
Exec Total Coverage
Lines: 224 362 61.9%
Functions: 11 19 57.9%
Branches: 161 374 43.0%
Line Branch Exec Source
1 /*
2 * Copyright (c) 2000-2022 Inria
3 * All rights reserved.
4 *
5 * Redistribution and use in source and binary forms, with or without
6 * modification, are permitted provided that the following conditions are met:
7 *
8 * * Redistributions of source code must retain the above copyright notice,
9 * this list of conditions and the following disclaimer.
10 * * Redistributions in binary form must reproduce the above copyright notice,
11 * this list of conditions and the following disclaimer in the documentation
12 * and/or other materials provided with the distribution.
13 * * Neither the name of the ALICE Project-Team nor the names of its
14 * contributors may be used to endorse or promote products derived from this
15 * software without specific prior written permission.
16 *
17 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
18 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
19 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
20 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
21 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
22 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
23 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
24 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
25 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
26 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
27 * POSSIBILITY OF SUCH DAMAGE.
28 *
29 * Contact: Bruno Levy
30 *
31 * https://www.inria.fr/fr/bruno-levy
32 *
33 * Inria,
34 * Domaine de Voluceau,
35 * 78150 Le Chesnay - Rocquencourt
36 * FRANCE
37 *
38 */
39
40 #include <geogram/delaunay/delaunay_3d.h>
41 #include <geogram/basic/logger.h>
42 #include <geogram/basic/geometry_nd.h>
43 #include <geogram/basic/process.h>
44 #include <geogram/basic/command_line.h>
45 #include <geogram/basic/stopwatch.h>
46 #include <geogram/basic/matrix.h>
47 #include <geogram/basic/permutation.h>
48 #include <geogram/basic/algorithm.h>
49 #include <geogram/mesh/mesh_reorder.h>
50 #include <geogram/bibliography/bibliography.h>
51 #include <stack>
52
53 // TODO: optimizations:
54 // - convex hull traversal for nearest_vertex()
55
56 namespace GEO {
57
58 char Delaunay3d::halfedge_facet_[4][4] = {
59 {4, 2, 3, 1},
60 {3, 4, 0, 2},
61 {1, 3, 4, 0},
62 {2, 0, 1, 4}
63 };
64
65 // tet facet vertex is such that the tetrahedron
66 // formed with:
67 // vertex lv
68 // tet_facet_vertex[lv][0]
69 // tet_facet_vertex[lv][1]
70 // tet_facet_vertex[lv][2]
71 // has the same orientation as the original tetrahedron for
72 // any vertex lv.
73
74 char Delaunay3d::tet_facet_vertex_[4][3] = {
75 {1, 2, 3},
76 {0, 3, 2},
77 {3, 0, 1},
78 {1, 0, 2}
79 };
80
81 2 Delaunay3d::Delaunay3d(coord_index_t dimension) :
82
1/2
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
 2 Delaunay(dimension)
83 {
84
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 geo_cite_with_info(
85 "DBLP:journals/cj/Bowyer81",
86 "One of the two initial references to the algorithm, "
87 "discovered independently and simultaneously by Bowyer and Watson."
88 );
89
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 geo_cite_with_info(
90 "journals/cj/Watson81",
91 "One of the two initial references to the algorithm, "
92 "discovered independently and simultaneously by Bowyer and Watson."
93 );
94
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 geo_cite_with_info(
95 "DBLP:conf/compgeom/AmentaCR03",
96 "Using spatial sorting has a dramatic impact on the performances."
97 );
98
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 geo_cite_with_info(
99 "DBLP:journals/comgeo/FunkeMN05",
100 "Initializing \\verb|locate()| with a non-exact version "
101 " (structural filtering) gains (a bit of) performance."
102 );
103
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 geo_cite_with_info(
104 "DBLP:journals/comgeo/BoissonnatDPTY02",
105 "The idea of traversing the cavity from inside "
106 " used in GEOGRAM is inspired by the implementation of "
107 " \\verb|Delaunay_triangulation_3| in CGAL."
108 );
109
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 geo_cite_with_info(
110 "DBLP:conf/imr/Si06",
111 "The triangulation data structure used in GEOGRAM is inspired "
112 "by Tetgen."
113 );
114
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 geo_cite_with_info(
115 "DBLP:journals/ijfcs/DevillersPT02",
116 "Analysis of the different versions of the line walk algorithm "
117 " used by \\verb|locate()|."
118 );
119
120
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 if(dimension != 3 && dimension != 4) {
121 ✗ throw InvalidDimension(dimension, "Delaunay3d", "3 or 4");
122 }
123 2 first_free_ = END_OF_LIST;
124 2 weighted_ = (dimension == 4);
125 // In weighted mode, vertices are 4d but combinatorics is 3d.
126
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 if(weighted_) {
127 ✗ cell_size_ = 4;
128 ✗ cell_v_stride_ = 4;
129 ✗ cell_neigh_stride_ = 4;
130 }
131 2 cur_stamp_ = 0;
132
3/6
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
 2 debug_mode_ = CmdLine::get_arg_bool("dbg:delaunay");
133
3/6
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
 2 verbose_debug_mode_ = CmdLine::get_arg_bool("dbg:delaunay_verbose");
134
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 debug_mode_ = (debug_mode_ || verbose_debug_mode_);
135
3/8
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
 2 benchmark_mode_ = CmdLine::get_arg_bool("dbg:delaunay_benchmark");
136 2 }
137
138
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1 times.
 8 Delaunay3d::~Delaunay3d() {
139 8 }
140
141 2 void Delaunay3d::set_vertices(index_t nb_vertices, const double* vertices) {
142 Stopwatch* W = nullptr;
143
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 if(benchmark_mode_) {
144 ✗ W = new Stopwatch("DelInternal");
145 }
146 2 cur_stamp_ = 0;
147
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 if(weighted_) {
148 ✗ heights_.resize(nb_vertices);
149 ✗ for(index_t i = 0; i < nb_vertices; ++i) {
150 // Client code uses 4d embedding with ti = sqrt(W - wi)
151 // where W = max(wi)
152 // We recompute the standard "shifted" lifting on
153 // the paraboloid from it.
154 // (we use wi - W, everything is shifted by W, but
155 // we do not care since the power diagram is invariant
156 // by a translation of all weights).
157 ✗ double w = -geo_sqr(vertices[4 * i + 3]);
158 ✗ heights_[i] = -w +
159 ✗ geo_sqr(vertices[4 * i]) +
160 ✗ geo_sqr(vertices[4 * i + 1]) +
161 ✗ geo_sqr(vertices[4 * i + 2]);
162 }
163 }
164
165 2 Delaunay::set_vertices(nb_vertices, vertices);
166
167 2 index_t expected_tetra = nb_vertices * 7;
168
169 2 cell_to_v_store_.reserve(expected_tetra * 4);
170 2 cell_to_cell_store_.reserve(expected_tetra * 4);
171 2 cell_next_.reserve(expected_tetra);
172
173 2 cell_to_v_store_.resize(0);
174 2 cell_to_cell_store_.resize(0);
175 2 cell_next_.resize(0);
176 2 first_free_ = END_OF_LIST;
177
178 // Sort the vertices spatially. This makes localisation
179 // faster.
180
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 if(do_reorder_) {
181 2 compute_BRIO_order(
182 2 nb_vertices, vertex_ptr(0), reorder_, 3, dimension_
183 );
184 } else {
185 ✗ reorder_.resize(nb_vertices);
186 ✗ for(index_t i = 0; i < nb_vertices; ++i) {
187 ✗ reorder_[i] = i;
188 }
189 }
190
191 double sorting_time = 0;
192
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 if(benchmark_mode_) {
193 ✗ sorting_time = W->elapsed_time();
194 ✗ Logger::out("DelInternal1") << "BRIO sorting:"
195 << sorting_time
196 << std::endl;
197 }
198
199 // The indices of the vertices of the first tetrahedron.
200 index_t v0, v1, v2, v3;
201
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
 2 if(!create_first_tetrahedron(v0, v1, v2, v3)) {
202 ✗ Logger::warn("Delaunay3d") << "All the points are coplanar"
203 << std::endl;
204 ✗ return;
205 }
206
207 index_t hint = NO_TETRAHEDRON;
208 // Insert all the vertices incrementally.
209
2/2
✓ Branch 0 taken 229 times.
✓ Branch 1 taken 2 times.
 231 for(index_t i = 0; i < nb_vertices; ++i) {
210 229 index_t v = reorder_[i];
211 // Do not re-insert the first four vertices.
212
8/8
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 227 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 225 times.
✓ Branch 4 taken 2 times.
✓ Branch 5 taken 223 times.
✓ Branch 6 taken 2 times.
✓ Branch 7 taken 221 times.
 229 if(v != v0 && v != v1 && v != v2 && v != v3) {
213 221 index_t new_hint = insert(v, hint);
214
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 221 times.
 221 if(new_hint != NO_TETRAHEDRON) {
215 hint = new_hint;
216 }
217 }
218 }
219
220
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 if(benchmark_mode_) {
221 ✗ Logger::out("DelInternal2") << "Core insertion algo:"
222 ✗ << W->elapsed_time() - sorting_time
223 << std::endl;
224 }
225
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 delete W;
226
227
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 if(debug_mode_) {
228 ✗ check_combinatorics(verbose_debug_mode_);
229 ✗ check_geometry(verbose_debug_mode_);
230 }
231
232 // Compress cell_to_v_store_ and cell_to_cell_store_
233 // (remove free and virtual tetrahedra).
234 // Since cell_next_ is not used at this point,
235 // we reuse it for storing the conversion array that
236 // maps old tet indices to new tet indices
237 // Note: tet_is_real() uses the previous value of
238 // cell_next(), but we are processing indices
239 // in increasing order and since old2new[t] is always
240 // smaller or equal to t, we never overwrite a value
241 // before needing it.
242
243 vector<index_t>& old2new = cell_next_;
244 index_t nb_tets = 0;
245 index_t nb_tets_to_delete = 0;
246
247 {
248
2/2
✓ Branch 0 taken 1186 times.
✓ Branch 1 taken 2 times.
 2374 for(index_t t = 0; t < max_t(); ++t) {
249 if(
250
3/6
✗ Branch 0 not taken.
✓ Branch 1 taken 1186 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 779 times.
✓ Branch 5 taken 366 times.
 2331 (keep_infinite_ && !tet_is_free(t)) ||
251 tet_is_real(t)
252 ) {
253
1/2
✓ Branch 0 taken 779 times.
✗ Branch 1 not taken.
 779 if(t != nb_tets) {
254 Memory::copy(
255 &cell_to_v_store_[nb_tets * 4],
256 &cell_to_v_store_[t * 4],
257 4 * sizeof(index_t)
258 );
259 Memory::copy(
260 &cell_to_cell_store_[nb_tets * 4],
261 &cell_to_cell_store_[t * 4],
262 4 * sizeof(index_t)
263 );
264 }
265 779 old2new[t] = nb_tets;
266 779 ++nb_tets;
267 } else {
268 407 old2new[t] = NO_INDEX;
269 407 ++nb_tets_to_delete;
270 }
271 }
272 2 cell_to_v_store_.resize(4 * nb_tets);
273 2 cell_to_cell_store_.resize(4 * nb_tets);
274
2/2
✓ Branch 0 taken 3116 times.
✓ Branch 1 taken 2 times.
 3118 for(index_t i = 0; i < 4 * nb_tets; ++i) {
275 3116 index_t t = cell_to_cell_store_[i];
276 geo_debug_assert(t != NO_INDEX);
277 3116 t = old2new[t];
278 // Note: t can be equal to -1 when a real tet is
279 // adjacent to a virtual one (and this is how the
280 // rest of Vorpaline expects to see tets on the
281 // border).
282 3116 cell_to_cell_store_[i] = t;
283 }
284 }
285
286 // In "keep_infinite" mode, we reorder the cells in such
287 // a way that finite cells have indices [0..nb_finite_cells_-1]
288 // and infinite cells have indices [nb_finite_cells_ .. nb_cells_-1]
289
290
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 if(keep_infinite_) {
291 ✗ nb_finite_cells_ = 0;
292 index_t finite_ptr = 0;
293 ✗ index_t infinite_ptr = nb_tets - 1;
294 for(;;) {
295 ✗ while(tet_is_finite(finite_ptr)) {
296 ✗ old2new[finite_ptr] = finite_ptr;
297 ✗ ++finite_ptr;
298 ✗ ++nb_finite_cells_;
299 }
300 ✗ while(!tet_is_finite(infinite_ptr)) {
301 ✗ old2new[infinite_ptr] = infinite_ptr;
302 ✗ --infinite_ptr;
303 }
304 ✗ if(finite_ptr > infinite_ptr) {
305 break;
306 }
307 ✗ old2new[finite_ptr] = infinite_ptr;
308 ✗ old2new[infinite_ptr] = finite_ptr;
309 ✗ ++nb_finite_cells_;
310 ✗ for(index_t lf=0; lf<4; ++lf) {
311 ✗ std::swap(
312 ✗ cell_to_cell_store_[4*finite_ptr + lf],
313 ✗ cell_to_cell_store_[4*infinite_ptr + lf]
314 );
315 }
316 ✗ for(index_t lv=0; lv<4; ++lv) {
317 ✗ std::swap(
318 ✗ cell_to_v_store_[4*finite_ptr + lv],
319 ✗ cell_to_v_store_[4*infinite_ptr + lv]
320 );
321 }
322 ✗ ++finite_ptr;
323 ✗ --infinite_ptr;
324 ✗ }
325 ✗ for(index_t i = 0; i < 4 * nb_tets; ++i) {
326 ✗ index_t t = cell_to_cell_store_[i];
327 geo_debug_assert(t != NO_INDEX);
328 ✗ t = old2new[t];
329 geo_debug_assert(t != NO_INDEX);
330 ✗ cell_to_cell_store_[i] = t;
331 }
332 }
333
334
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 if(benchmark_mode_) {
335 ✗ if(keep_infinite_) {
336 ✗ Logger::out("DelCompress")
337 << "Removed " << nb_tets_to_delete
338 << " tets (free list)" << std::endl;
339 } else {
340 ✗ Logger::out("DelCompress")
341 << "Removed " << nb_tets_to_delete
342 << " tets (free list and infinite)" << std::endl;
343 }
344 }
345
346 2 set_arrays(
347 nb_tets,
348 cell_to_v_store_.data(), cell_to_cell_store_.data()
349 );
350 }
351
352 ✗ index_t Delaunay3d::nearest_vertex(const double* p) const {
353
354 // TODO: For the moment, we fallback to the (unefficient)
355 // baseclass implementation when in weighted mode.
356 ✗ if(weighted_) {
357 ✗ return Delaunay::nearest_vertex(p);
358 }
359
360 // Find a tetrahedron (real or virtual) that contains p
361 ✗ index_t t = locate(p, NO_TETRAHEDRON, thread_safe());
362
363 // If p is outside the convex hull of the inserted points,
364 // a special traversal is required (not implemented yet).
365 // TODO: implement convex hull boundary traversal
366 // (for now we fallback to linear search implemented
367 // in baseclass)
368 ✗ if(t == NO_TETRAHEDRON || tet_is_virtual(t)) {
369 ✗ return Delaunay::nearest_vertex(p);
370 }
371
372 double sq_dist = 1e30;
373 index_t result = NO_TETRAHEDRON;
374
375 // Find the nearest vertex among t's vertices
376 ✗ for(index_t lv = 0; lv < 4; ++lv) {
377 index_t v = tet_vertex(t, lv);
378 // If the tetrahedron is virtual, then the first vertex
379 // is the vertex at infinity and is skipped.
380 ✗ if(v == NO_INDEX) {
381 ✗ continue;
382 }
383 double cur_sq_dist = Geom::distance2(p, vertex_ptr(v), 3);
384 ✗ if(cur_sq_dist < sq_dist) {
385 sq_dist = cur_sq_dist;
386 result = v;
387 }
388 }
389 return result;
390 }
391
392
393
394 221 index_t Delaunay3d::locate_inexact(
395 const double* p, index_t hint, index_t max_iter
396 ) const {
397
398 // If no hint specified, find a tetrahedron randomly
399
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 219 times.
 221 while(hint == NO_TETRAHEDRON) {
400
2/4
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 2 times.
 2 hint = index_t(Numeric::random_int32()) % max_t();
401
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 if(tet_is_free(hint)) {
402 hint = NO_TETRAHEDRON;
403 }
404 }
405
406 // Always start from a real tet. If the tet is virtual,
407 // find its real neighbor (always opposite to the
408 // infinite vertex)
409
2/2
✓ Branch 0 taken 152 times.
✓ Branch 1 taken 69 times.
 221 if(tet_is_virtual(hint)) {
410
1/2
✓ Branch 0 taken 205 times.
✗ Branch 1 not taken.
 205 for(index_t lf = 0; lf < 4; ++lf) {
411
2/2
✓ Branch 0 taken 69 times.
✓ Branch 1 taken 136 times.
 205 if(tet_vertex(hint, lf) == VERTEX_AT_INFINITY) {
412 hint = tet_adjacent(hint, lf);
413 geo_debug_assert(hint != NO_TETRAHEDRON);
414 69 break;
415 }
416 }
417 }
418
419 index_t t = hint;
420 index_t t_pred = NO_TETRAHEDRON;
421
422 221 still_walking:
423 {
424 const double* pv[4];
425 1542 pv[0] = vertex_ptr(finite_tet_vertex(t,0));
426 1542 pv[1] = vertex_ptr(finite_tet_vertex(t,1));
427 1542 pv[2] = vertex_ptr(finite_tet_vertex(t,2));
428 1542 pv[3] = vertex_ptr(finite_tet_vertex(t,3));
429
430
2/2
✓ Branch 0 taken 3777 times.
✓ Branch 1 taken 116 times.
 3893 for(index_t f = 0; f < 4; ++f) {
431
432 index_t t_next = tet_adjacent(t,f);
433
434 // If the opposite tet is -1, then it means that
435 // we are trying to locate() (e.g. called from
436 // nearest_vertex) within a tetrahedralization
437 // from which the infinite tets were removed.
438
1/2
✓ Branch 0 taken 3777 times.
✗ Branch 1 not taken.
 3777 if(t_next == NO_INDEX) {
439 105 return NO_TETRAHEDRON;
440 }
441
442 // If the candidate next tetrahedron is the
443 // one we came from, then we know already that
444 // the orientation is positive, thus we examine
445 // the next candidate (or exit the loop if they
446 // are exhausted).
447
2/2
✓ Branch 0 taken 637 times.
✓ Branch 1 taken 3140 times.
 3777 if(t_next == t_pred) {
448 637 continue ;
449 }
450
451 // To test the orientation of p w.r.t. the facet f of
452 // t, we replace vertex number f with p in t (same
453 // convention as in CGAL).
454 3140 const double* pv_bkp = pv[f];
455 3140 pv[f] = p;
456 3140 Sign ori = PCK::orient_3d_inexact(pv[0], pv[1], pv[2], pv[3]);
457
458 // If the orientation is not negative, then we cannot
459 // walk towards t_next, and examine the next candidate
460 // (or exit the loop if they are exhausted).
461
2/2
✓ Branch 0 taken 1714 times.
✓ Branch 1 taken 1426 times.
 3140 if(ori != NEGATIVE) {
462 1714 pv[f] = pv_bkp;
463 1714 continue;
464 }
465
466 // If the opposite tet is a virtual tet, then
467 // the point has a positive orientation relative
468 // to the facet on the border of the convex hull,
469 // thus t_next is a tet in conflict and we are
470 // done.
471
2/2
✓ Branch 0 taken 1321 times.
✓ Branch 1 taken 105 times.
 1426 if(tet_is_virtual(t_next)) {
472 return t_next;
473 }
474
475 // If we reach this point, then t_next is a valid
476 // successor, thus we are still walking.
477 t_pred = t;
478 t = t_next;
479
1/2
✓ Branch 0 taken 1321 times.
✗ Branch 1 not taken.
 1321 if(--max_iter != 0) {
480 1321 goto still_walking;
481 }
482 }
483 }
484
485 // If we reach this point, we did not find a valid successor
486 // for walking (a face for which p has negative orientation),
487 // thus we reached the tet for which p has all positive
488 // face orientations (i.e. the tet that contains p).
489
490 116 return t;
491 }
492
493
494 221 index_t Delaunay3d::locate(
495 const double* p, index_t hint, bool thread_safe,
496 Sign* orient
497 ) const {
498
499 // Try improving the hint by using the
500 // inexact locate function. This gains
501 // (a little bit) performance (a few
502 // percent in total Delaunay computation
503 // time), but it is better than nothing...
504 // Note: there is a maximum number of tets
505 // traversed by locate_inexact() (2500)
506 // since there exists configurations in which
507 // locate_inexact() loops forever !
508 221 hint = locate_inexact(p, hint, 2500);
509
510 static Process::spinlock locate_lock = GEOGRAM_SPINLOCK_INIT;
511
512 // We need to have this spinlock because
513 // of random() that is not thread-safe
514 // (TODO: implement a random() function with
515 // thread local storage)
516
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 221 times.
 221 if(thread_safe) {
517 Process::acquire_spinlock(locate_lock);
518 }
519
520 // If no hint specified, find a tetrahedron randomly
521
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 221 times.
 221 while(hint == NO_TETRAHEDRON) {
522 ✗ hint = index_t(Numeric::random_int32()) % max_t();
523 ✗ if(tet_is_free(hint)) {
524 hint = NO_TETRAHEDRON;
525 }
526 }
527
528 // Always start from a real tet. If the tet is virtual,
529 // find its real neighbor (always opposite to the
530 // infinite vertex)
531
2/2
✓ Branch 0 taken 116 times.
✓ Branch 1 taken 105 times.
 221 if(tet_is_virtual(hint)) {
532
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 296 times.
 296 for(index_t lf = 0; lf < 4; ++lf) {
533
2/2
✓ Branch 0 taken 105 times.
✓ Branch 1 taken 191 times.
 296 if(tet_vertex(hint, lf) == VERTEX_AT_INFINITY) {
534 hint = tet_adjacent(hint, lf);
535 geo_debug_assert(hint != NO_TETRAHEDRON);
536 105 break;
537 }
538 }
539 }
540
541 index_t t = hint;
542 index_t t_pred = NO_TETRAHEDRON;
543 Sign orient_local[4];
544
1/2
✓ Branch 0 taken 221 times.
✗ Branch 1 not taken.
 221 if(orient == nullptr) {
545 orient = orient_local;
546 }
547
548
549 221 still_walking:
550 {
551 const double* pv[4];
552 226 pv[0] = vertex_ptr(finite_tet_vertex(t,0));
553 226 pv[1] = vertex_ptr(finite_tet_vertex(t,1));
554 226 pv[2] = vertex_ptr(finite_tet_vertex(t,2));
555 226 pv[3] = vertex_ptr(finite_tet_vertex(t,3));
556
557 // Start from a random facet
558 226 index_t f0 = index_t(Numeric::random_int32()) % 4;
559
2/2
✓ Branch 0 taken 711 times.
✓ Branch 1 taken 116 times.
 827 for(index_t df = 0; df < 4; ++df) {
560
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 711 times.
 711 index_t f = (f0 + df) % 4;
561
562 index_t t_next = tet_adjacent(t,f);
563
564 // If the opposite tet is -1, then it means that
565 // we are trying to locate() (e.g. called from
566 // nearest_vertex) within a tetrahedralization
567 // from which the infinite tets were removed.
568
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 711 times.
 711 if(t_next == NO_INDEX) {
569 ✗ if(thread_safe) {
570 Process::release_spinlock(locate_lock);
571 }
572 105 return NO_TETRAHEDRON;
573 }
574
575 // If the candidate next tetrahedron is the
576 // one we came from, then we know already that
577 // the orientation is positive, thus we examine
578 // the next candidate (or exit the loop if they
579 // are exhausted).
580
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 710 times.
 711 if(t_next == t_pred) {
581 1 orient[f] = POSITIVE ;
582 1 continue ;
583 }
584
585 // To test the orientation of p w.r.t. the facet f of
586 // t, we replace vertex number f with p in t (same
587 // convention as in CGAL).
588 // This is equivalent to tet_facet_point_orient3d(t,f,p)
589 // (but less costly, saves a couple of lookups)
590 710 const double* pv_bkp = pv[f];
591 710 pv[f] = p;
592 710 orient[f] = PCK::orient_3d(pv[0], pv[1], pv[2], pv[3]);
593
594 // If the orientation is not negative, then we cannot
595 // walk towards t_next, and examine the next candidate
596 // (or exit the loop if they are exhausted).
597
2/2
✓ Branch 0 taken 600 times.
✓ Branch 1 taken 110 times.
 710 if(orient[f] != NEGATIVE) {
598 600 pv[f] = pv_bkp;
599 600 continue;
600 }
601
602 // If the opposite tet is a virtual tet, then
603 // the point has a positive orientation relative
604 // to the facet on the border of the convex hull,
605 // thus t_next is a tet in conflict and we are
606 // done.
607
2/2
✓ Branch 0 taken 105 times.
✓ Branch 1 taken 5 times.
 110 if(tet_is_virtual(t_next)) {
608
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 105 times.
 105 if(thread_safe) {
609 Process::release_spinlock(locate_lock);
610 }
611
2/2
✓ Branch 0 taken 420 times.
✓ Branch 1 taken 105 times.
 525 for(index_t lf = 0; lf < 4; ++lf) {
612 420 orient[lf] = POSITIVE;
613 }
614 return t_next;
615 }
616
617 // If we reach this point, then t_next is a valid
618 // successor, thus we are still walking.
619 t_pred = t;
620 t = t_next;
621 5 goto still_walking;
622 }
623 }
624
625 // If we reach this point, we did not find a valid successor
626 // for walking (a face for which p has negative orientation),
627 // thus we reached the tet for which p has all positive
628 // face orientations (i.e. the tet that contains p).
629
630
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 116 times.
 116 if(thread_safe) {
631 Process::release_spinlock(locate_lock);
632 }
633 return t;
634 }
635
636
1/2
✓ Branch 0 taken 221 times.
✗ Branch 1 not taken.
 221 void Delaunay3d::find_conflict_zone(
637 index_t v,
638 index_t t, const Sign* orient,
639 index_t& t_bndry, index_t& f_bndry,
640 index_t& first, index_t& last
641 ) {
642 cavity_.clear();
643
644
1/2
✓ Branch 0 taken 221 times.
✗ Branch 1 not taken.
 221 first = last = END_OF_LIST;
645
646 // Generate a unique stamp from current vertex index,
647 // used for marking tetrahedra.
648 set_tet_mark_stamp(v);
649
650 // Pointer to the coordinates of the point to be inserted
651 const double* p = vertex_ptr(v);
652
653 geo_debug_assert(t != NO_TETRAHEDRON);
654
655 // Test whether the point already exists in
656 // the triangulation. The point already exists
657 // if it's located on three faces of the
658 // tetrahedron returned by locate().
659 221 int nb_zero =
660 221 (orient[0] == ZERO) +
661 221 (orient[1] == ZERO) +
662 221 (orient[2] == ZERO) +
663 221 (orient[3] == ZERO) ;
664
665
1/2
✓ Branch 0 taken 221 times.
✗ Branch 1 not taken.
 221 if(nb_zero >= 3) {
666 return;
667 }
668
669 // Weighted triangulations can have dangling
670 // vertices. Such vertices p are characterized by
671 // the fact that p is not in conflict with the
672 // tetrahedron returned by locate().
673
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 221 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 221 if(weighted_ && !tet_is_conflict(t, p)) {
674 return;
675 }
676
677 // Note: points on edges and on facets are
678 // handled by the way tet_is_in_conflict()
679 // is implemented, that naturally inserts
680 // the correct tetrahedra in the conflict list.
681
682
683 // Mark t as conflict
684 add_tet_to_list(t, first, last);
685
686 // A small optimization: if the point to be inserted
687 // is on some faces of the located tetrahedron, insert
688 // the neighbors accros those faces in the conflict list.
689 // It saves a couple of calls to the predicates in this
690 // specific case (combinatorics are in general less
691 // expensive than the predicates).
692
3/4
✓ Branch 0 taken 221 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 95 times.
✓ Branch 3 taken 126 times.
 221 if(!weighted_ && nb_zero != 0) {
693
2/2
✓ Branch 0 taken 380 times.
✓ Branch 1 taken 95 times.
 475 for(index_t lf = 0; lf < 4; ++lf) {
694
2/2
✓ Branch 0 taken 96 times.
✓ Branch 1 taken 284 times.
 380 if(orient[lf] == ZERO) {
695 index_t t2 = tet_adjacent(t, lf);
696 add_tet_to_list(t2, first, last);
697 }
698 }
699
2/2
✓ Branch 0 taken 380 times.
✓ Branch 1 taken 95 times.
 475 for(index_t lf = 0; lf < 4; ++lf) {
700
2/2
✓ Branch 0 taken 96 times.
✓ Branch 1 taken 284 times.
 380 if(orient[lf] == ZERO) {
701 index_t t2 = tet_adjacent(t, lf);
702 96 find_conflict_zone_iterative(
703 p,t2,t_bndry,f_bndry,first,last
704 );
705 }
706 }
707 }
708
709 // Determine the conflict list by greedy propagation from t.
710 221 find_conflict_zone_iterative(p,t,t_bndry,f_bndry,first,last);
711 }
712
713
1/2
✓ Branch 0 taken 317 times.
✗ Branch 1 not taken.
 317 void Delaunay3d::find_conflict_zone_iterative(
714 const double* p, index_t t_in,
715 index_t& t_bndry, index_t& f_bndry,
716 index_t& first, index_t& last
717 ) {
718
719 //std::stack<index_t> S;
720 S_.push(t_in);
721
722
2/2
✓ Branch 0 taken 3313 times.
✓ Branch 1 taken 317 times.
 3630 while(!S_.empty()) {
723
724
2/2
✓ Branch 0 taken 3311 times.
✓ Branch 1 taken 2 times.
 3313 index_t t = S_.top();
725 S_.pop();
726
727
2/2
✓ Branch 0 taken 13252 times.
✓ Branch 1 taken 3313 times.
 16565 for(index_t lf = 0; lf < 4; ++lf) {
728
2/2
✓ Branch 0 taken 5808 times.
✓ Branch 1 taken 7444 times.
 13252 index_t t2 = tet_adjacent(t, lf);
729
730
2/2
✓ Branch 0 taken 5808 times.
✓ Branch 1 taken 7444 times.
 13252 if(
731 tet_is_in_list(t2) // known as conflict
732 ) {
733 9300 continue;
734 }
735
736
2/2
✓ Branch 0 taken 496 times.
✓ Branch 1 taken 6948 times.
 7444 if(
737 tet_is_marked(t2) // known as non-conflict
738 ) {
739 496 cavity_.new_facet(
740 t, lf,
741 tet_vertex(t, tet_facet_vertex(lf,0)),
742 tet_vertex(t, tet_facet_vertex(lf,1)),
743 tet_vertex(t, tet_facet_vertex(lf,2))
744 );
745 496 continue;
746 }
747
748
749
2/2
✓ Branch 1 taken 2996 times.
✓ Branch 2 taken 3952 times.
 6948 if(tet_is_conflict(t2, p)) {
750 // Chain t2 in conflict list
751 add_tet_to_list(t2, first, last);
752 S_.push(t2);
753 2996 continue;
754 }
755
756 // At this point, t is in conflict
757 // and t2 is not in conflict.
758 // We keep a reference to a tet on the boundary
759 3952 t_bndry = t;
760 3952 f_bndry = lf;
761 // Mark t2 as visited (but not conflict)
762 mark_tet(t2);
763
764 3952 cavity_.new_facet(
765 t, lf,
766 tet_vertex(t, tet_facet_vertex(lf,0)),
767 tet_vertex(t, tet_facet_vertex(lf,1)),
768 tet_vertex(t, tet_facet_vertex(lf,2))
769 );
770
771 }
772 }
773 317 }
774
775
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
 3 index_t Delaunay3d::stellate_conflict_zone_iterative(
776 index_t v, index_t t1, index_t t1fbord, index_t t1fprev
777 ) {
778 // This function is de-recursified because some degenerate
779 // inputs can cause stack overflow (system stack is limited to
780 // a few megs). For instance, it can happen when a large number
781 // of points are on the same sphere exactly.
782
783 // To de-recursify, it uses class StellateConflictStack
784 // that emulates system's stack for storing functions's
785 // parameters and local variables in all the nested stack
786 // frames.
787
788 S2_.push(t1, t1fbord, t1fprev);
789
790 index_t new_t; // the newly created tetrahedron.
791
792 index_t t1ft2; // traverses the 4 facets of t1.
793
794 index_t t2; // the tetrahedron on the border of
795 // the conflict zone that shares an
796 // edge with t1 along t1ft2.
797
798 index_t t2fbord; // the facet of t2 on the border of
799 // the conflict zone.
800
801 index_t t2ft1; // the facet of t2 that is incident to t1.
802
803 500 entry_point:
804 S2_.get_parameters(t1, t1fbord, t1fprev);
805
806 geo_debug_assert(tet_is_in_list(t1));
807 geo_debug_assert(tet_adjacent(t1,t1fbord) != NO_INDEX);
808 geo_debug_assert(!tet_is_in_list(tet_adjacent(t1,t1fbord)));
809
810 // Create new tetrahedron with same vertices as t_bndry
811 500 new_t = new_tetrahedron(
812 tet_vertex(t1,0),
813 tet_vertex(t1,1),
814 tet_vertex(t1,2),
815 tet_vertex(t1,3)
816 );
817
818 // Replace in new_t the vertex opposite to t1fbord with v
819 set_tet_vertex(new_t, t1fbord, v);
820
821 // Connect new_t with t1's neighbor accros t1fbord
822 {
823 index_t tbord = tet_adjacent(t1,t1fbord);
824 set_tet_adjacent(new_t, t1fbord, tbord);
825 set_tet_adjacent(tbord, find_tet_adjacent(tbord,t1), new_t);
826 }
827
828 // Lookup new_t's neighbors accros its three other
829 // facets and connect them
830
2/2
✓ Branch 0 taken 2000 times.
✓ Branch 1 taken 500 times.
 2500 for(t1ft2=0; t1ft2<4; ++t1ft2) {
831
832
4/4
✓ Branch 0 taken 1503 times.
✓ Branch 1 taken 497 times.
✓ Branch 2 taken 753 times.
✓ Branch 3 taken 750 times.
 2000 if(t1ft2 == t1fprev || tet_adjacent(new_t,t1ft2) != NO_INDEX) {
833 1250 continue;
834 }
835
836 // Get t1's neighbor along the border of the conflict zone
837
2/2
✓ Branch 1 taken 497 times.
✓ Branch 2 taken 253 times.
 750 if(!get_neighbor_along_conflict_zone_border(
838 t1,t1fbord,t1ft2, t2,t2fbord,t2ft1
839 )) {
840 // If t1's neighbor is not a new tetrahedron,
841 // create a new tetrahedron through a recursive call.
842
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 497 times.
 497 S2_.save_locals(new_t, t1ft2, t2ft1);
843
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 497 times.
 497 S2_.push(t2, t2fbord, t2ft1);
844 497 goto entry_point;
845
846 return_point:
847 // This is the return value of the called function.
848 index_t result = new_t;
849 S2_.pop();
850
851 // Special case: we were in the outermost frame,
852 // then we (truly) return from the function.
853 500 if(S2_.empty()) {
854 3 return result;
855 }
856
857 S2_.get_parameters(t1, t1fbord, t1fprev);
858 S2_.get_locals(new_t, t1ft2, t2ft1);
859 497 t2 = result;
860 }
861
862 750 set_tet_adjacent(t2, t2ft1, new_t);
863 set_tet_adjacent(new_t, t1ft2, t2);
864 }
865
866 // Except for the initial call (see "Special case" above),
867 // the nested calls all come from the same location,
868 // thus there is only one possible return point
869 // (no need to push any return address).
870 500 goto return_point;
871 }
872
873 218 index_t Delaunay3d::stellate_cavity(index_t v) {
874
875 index_t new_tet = NO_INDEX;
876
877
2/2
✓ Branch 0 taken 3948 times.
✓ Branch 1 taken 218 times.
 4166 for(index_t f=0; f<cavity_.nb_facets(); ++f) {
878 index_t old_tet = cavity_.facet_tet(f);
879 index_t lf = cavity_.facet_facet(f);
880 index_t t_neigh = tet_adjacent(old_tet, lf);
881 index_t v1 = cavity_.facet_vertex(f,0);
882 index_t v2 = cavity_.facet_vertex(f,1);
883 index_t v3 = cavity_.facet_vertex(f,2);
884 3948 new_tet = new_tetrahedron(v, v1, v2, v3);
885 set_tet_adjacent(new_tet, 0, t_neigh);
886 set_tet_adjacent(
887 t_neigh, find_tet_adjacent(t_neigh,old_tet), new_tet
888 );
889 cavity_.set_facet_tet(f, new_tet);
890 }
891
892
2/2
✓ Branch 0 taken 3948 times.
✓ Branch 1 taken 218 times.
 4166 for(index_t f=0; f<cavity_.nb_facets(); ++f) {
893 new_tet = cavity_.facet_tet(f);
894 index_t neigh1, neigh2, neigh3;
895 3948 cavity_.get_facet_neighbor_tets(f, neigh1, neigh2, neigh3);
896 3948 set_tet_adjacent(new_tet, 1, neigh1);
897 3948 set_tet_adjacent(new_tet, 2, neigh2);
898 3948 set_tet_adjacent(new_tet, 3, neigh3);
899 }
900
901 218 return new_tet;
902 }
903
904 221 index_t Delaunay3d::insert(index_t v, index_t hint) {
905 221 index_t t_bndry = NO_TETRAHEDRON;
906 221 index_t f_bndry = NO_INDEX;
907 221 index_t first_conflict = NO_TETRAHEDRON;
908 221 index_t last_conflict = NO_TETRAHEDRON;
909
910 const double* p = vertex_ptr(v);
911
912 Sign orient[4];
913 221 index_t t = locate(p, hint, false, orient);
914 221 find_conflict_zone(
915 v,t,orient,t_bndry,f_bndry,first_conflict,last_conflict
916 );
917
918 // The conflict list can be empty if:
919 // - Vertex v already exists in the triangulation
920 // - The triangulation is weighted and v is not visible
921
1/2
✓ Branch 0 taken 221 times.
✗ Branch 1 not taken.
 221 if(first_conflict == END_OF_LIST) {
922 return NO_TETRAHEDRON;
923 }
924
925 index_t new_tet = NO_INDEX;
926
2/2
✓ Branch 0 taken 218 times.
✓ Branch 1 taken 3 times.
 221 if(cavity_.OK()) {
927 218 new_tet = stellate_cavity(v);
928 } else {
929 3 new_tet = stellate_conflict_zone_iterative(v,t_bndry,f_bndry);
930 }
931
932 // Recycle the tetrahedra of the conflict zone.
933 221 cell_next_[last_conflict] = first_free_;
934 221 first_free_ = first_conflict;
935
936 // Return one of the newly created tets
937 221 return new_tet;
938 }
939
940
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 bool Delaunay3d::create_first_tetrahedron(
941 index_t& iv0, index_t& iv1, index_t& iv2, index_t& iv3
942 ) {
943
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 if(nb_vertices() < 4) {
944 return false;
945 }
946
947 2 iv0 = 0;
948
949 2 iv1 = 1;
950 2 while(
951
2/4
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
 4 iv1 < nb_vertices() &&
952 2 PCK::points_are_identical_3d(
953 vertex_ptr(iv0), vertex_ptr(iv1)
954 )
955 ) {
956 ✗ ++iv1;
957 }
958
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 if(iv1 == nb_vertices()) {
959 return false;
960 }
961
962 2 iv2 = iv1 + 1;
963 2 while(
964
2/4
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
 4 iv2 < nb_vertices() &&
965 2 PCK::points_are_colinear_3d(
966 vertex_ptr(iv0), vertex_ptr(iv1), vertex_ptr(iv2)
967 )
968 ) {
969 ✗ ++iv2;
970 }
971
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 if(iv2 == nb_vertices()) {
972 return false;
973 }
974
975 2 iv3 = iv2 + 1;
976 Sign s = ZERO;
977 2 while(
978
3/4
✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 2 times.
 8 iv3 < nb_vertices() &&
979 4 (s = PCK::orient_3d(
980 vertex_ptr(iv0), vertex_ptr(iv1),
981 vertex_ptr(iv2), vertex_ptr(iv3)
982 )) == ZERO
983 ) {
984 2 ++iv3;
985 }
986
987
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 if(iv3 == nb_vertices()) {
988 return false;
989 }
990
991 geo_debug_assert(s != ZERO);
992
993
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1 times.
 2 if(s == NEGATIVE) {
994 std::swap(iv2, iv3);
995 }
996
997 // Create the first tetrahedron
998 2 index_t t0 = new_tetrahedron(iv0, iv1, iv2, iv3);
999
1000 // Create the first four virtual tetrahedra surrounding it
1001 index_t t[4];
1002
2/2
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 2 times.
 10 for(index_t f = 0; f < 4; ++f) {
1003 // In reverse order since it is an adjacent tetrahedron
1004 index_t v1 = tet_vertex(t0, tet_facet_vertex(f,2));
1005 index_t v2 = tet_vertex(t0, tet_facet_vertex(f,1));
1006 index_t v3 = tet_vertex(t0, tet_facet_vertex(f,0));
1007 8 t[f] = new_tetrahedron(VERTEX_AT_INFINITY, v1, v2, v3);
1008 }
1009
1010 // Connect the virtual tetrahedra to the real one
1011
2/2
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 2 times.
 10 for(index_t f=0; f<4; ++f) {
1012 8 set_tet_adjacent(t[f], 0, t0);
1013 set_tet_adjacent(t0, f, t[f]);
1014 }
1015
1016 // Interconnect the four virtual tetrahedra along their common
1017 // faces
1018
2/2
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 2 times.
 10 for(index_t f = 0; f < 4; ++f) {
1019 // In reverse order since it is an adjacent tetrahedron
1020 index_t lv1 = tet_facet_vertex(f,2);
1021 index_t lv2 = tet_facet_vertex(f,1);
1022 index_t lv3 = tet_facet_vertex(f,0);
1023 8 set_tet_adjacent(t[f], 1, t[lv1]);
1024 8 set_tet_adjacent(t[f], 2, t[lv2]);
1025 8 set_tet_adjacent(t[f], 3, t[lv3]);
1026 }
1027
1028 return true;
1029 }
1030
1031 /************************************************************************/
1032
1033 ✗ void Delaunay3d::show_tet(index_t t) const {
1034 std::cerr << "tet"
1035 ✗ << (tet_is_in_list(t) ? '*' : ' ')
1036 << t
1037 << ", v=["
1038 << tet_vertex(t, 0)
1039 ✗ << ' '
1040 << tet_vertex(t, 1)
1041 ✗ << ' '
1042 << tet_vertex(t, 2)
1043 ✗ << ' '
1044 << tet_vertex(t, 3)
1045 ✗ << "] adj=[";
1046 ✗ show_tet_adjacent(t, 0);
1047 ✗ show_tet_adjacent(t, 1);
1048 ✗ show_tet_adjacent(t, 2);
1049 ✗ show_tet_adjacent(t, 3);
1050 ✗ std::cerr << "] ";
1051
1052 ✗ for(index_t f = 0; f < 4; ++f) {
1053 ✗ std::cerr << 'f' << f << ':';
1054 ✗ for(index_t v = 0; v < 3; ++v) {
1055 std::cerr << tet_vertex(t, tet_facet_vertex(f,v))
1056 ✗ << ',';
1057 }
1058 ✗ std::cerr << ' ';
1059 }
1060 std::cerr << std::endl;
1061 ✗ }
1062
1063 ✗ void Delaunay3d::show_tet_adjacent(index_t t, index_t lf) const {
1064 index_t adj = tet_adjacent(t, lf);
1065 ✗ if(adj != NO_INDEX) {
1066 ✗ std::cerr << (tet_is_in_list(adj) ? '*' : ' ');
1067 }
1068 std::cerr << adj;
1069 ✗ std::cerr << ' ';
1070 ✗ }
1071
1072 ✗ void Delaunay3d::show_list(
1073 index_t first, const std::string& list_name
1074 ) const {
1075 index_t t = first;
1076 std::cerr << "tet list: " << list_name << std::endl;
1077 ✗ while(t != END_OF_LIST) {
1078 ✗ show_tet(t);
1079 t = tet_next(t);
1080 }
1081 std::cerr << "-------------" << std::endl;
1082 ✗ }
1083
1084 ✗ void Delaunay3d::check_combinatorics(bool verbose) const {
1085 ✗ if(verbose) {
1086 std::cerr << std::endl;
1087 }
1088 bool ok = true;
1089 ✗ std::vector<bool> v_has_tet(nb_vertices(), false);
1090 ✗ for(index_t t = 0; t < max_t(); ++t) {
1091 ✗ if(tet_is_free(t)) {
1092 /*
1093 if(verbose) {
1094 std::cerr << "-Deleted tet: ";
1095 show_tet(t);
1096 }
1097 */
1098 } else {
1099 /*
1100 if(verbose) {
1101 std::cerr << "Checking tet: ";
1102 show_tet(t);
1103 }
1104 */
1105 ✗ for(index_t lf = 0; lf < 4; ++lf) {
1106 ✗ if(tet_adjacent(t, lf) == NO_INDEX) {
1107 std::cerr << lf << ":Missing adjacent tet"
1108 << std::endl;
1109 ok = false;
1110 ✗ } else if(tet_adjacent(t, lf) == t) {
1111 std::cerr << lf << ":Tet is adjacent to itself"
1112 << std::endl;
1113 ok = false;
1114 } else {
1115 index_t t2 = tet_adjacent(t, lf);
1116 bool found = false;
1117 ✗ for(index_t lf2 = 0; lf2 < 4; ++lf2) {
1118 ✗ if(tet_adjacent(t2, lf2) == t) {
1119 found = true;
1120 }
1121 }
1122 ✗ if(!found) {
1123 std::cerr
1124 << lf << ":Adjacent link is not bidirectional"
1125 << std::endl;
1126 ok = false;
1127 }
1128 }
1129 }
1130 index_t nb_infinite = 0;
1131 ✗ for(index_t lv = 0; lv < 4; ++lv) {
1132 ✗ if(tet_vertex(t, lv) == NO_INDEX) {
1133 ✗ ++nb_infinite;
1134 }
1135 }
1136 ✗ if(nb_infinite > 1) {
1137 ok = false;
1138 std::cerr << "More than one infinite vertex"
1139 << std::endl;
1140 }
1141 }
1142 ✗ for(index_t lv = 0; lv < 4; ++lv) {
1143 index_t v = tet_vertex(t, lv);
1144 ✗ if(v != NO_INDEX) {
1145 v_has_tet[v] = true;
1146 }
1147 }
1148 }
1149 ✗ for(index_t v = 0; v < nb_vertices(); ++v) {
1150 ✗ if(!v_has_tet[v]) {
1151 ✗ if(verbose) {
1152 std::cerr << "Vertex " << v
1153 << " is isolated (duplicated ?)" << std::endl;
1154 }
1155 }
1156 }
1157 ✗ geo_assert(ok);
1158 ✗ if(verbose) {
1159 std::cerr << std::endl;
1160 }
1161 std::cerr << std::endl << "Delaunay Combi OK" << std::endl;
1162 ✗ }
1163
1164 ✗ void Delaunay3d::check_geometry(bool verbose) const {
1165 bool ok = true;
1166 ✗ for(index_t t = 0; t < max_t(); ++t) {
1167 ✗ if(!tet_is_free(t)) {
1168 index_t v0 = tet_vertex(t, 0);
1169 index_t v1 = tet_vertex(t, 1);
1170 index_t v2 = tet_vertex(t, 2);
1171 index_t v3 = tet_vertex(t, 3);
1172 ✗ for(index_t v = 0; v < nb_vertices(); ++v) {
1173 ✗ if(v == v0 || v == v1 || v == v2 || v == v3) {
1174 ✗ continue;
1175 }
1176 ✗ if(tet_is_conflict(t, vertex_ptr(v))) {
1177 ok = false;
1178 ✗ if(verbose) {
1179 std::cerr << "Tet " << t <<
1180 " is in conflict with vertex " << v
1181 << std::endl;
1182
1183 ✗ std::cerr << " offending tet: ";
1184 ✗ show_tet(t);
1185 }
1186 }
1187 }
1188 }
1189 }
1190 ✗ geo_assert(ok);
1191 std::cerr << std::endl << "Delaunay Geo OK" << std::endl;
1192 ✗ }
1193
1194 /************************************************************************/
1195
1196 ✗ RegularWeightedDelaunay3d::RegularWeightedDelaunay3d(
1197 coord_index_t dimension
1198 ✗ ) :
1199 ✗ Delaunay3d(4)
1200 {
1201 ✗ if(dimension != 4) {
1202 ✗ throw InvalidDimension(dimension, "RegularWeightedDelaunay3d", "4");
1203 }
1204 ✗ }
1205
1206 ✗ RegularWeightedDelaunay3d::~RegularWeightedDelaunay3d() {
1207 ✗ }
1208 }
1209