GCC Code Coverage Report

Directory: ./
File: voronoi/RVD.cpp
Date: 2026-09-27 03:24:14
Exec Total Coverage
Lines: 595 1082 55.0%
Functions: 127 579 21.9%
Branches: 261 938 27.8%
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/voronoi/RVD.h>
41 #include <geogram/voronoi/generic_RVD.h>
42 #include <geogram/voronoi/RVD_mesh_builder.h>
43 #include <geogram/voronoi/integration_simplex.h>
44 #include <geogram/voronoi/RVD_callback.h>
45 #include <geogram/mesh/mesh_partition.h>
46 #include <geogram/mesh/mesh_sampling.h>
47 #include <geogram/mesh/mesh_repair.h>
48 #include <geogram/mesh/mesh_AABB.h>
49 #include <geogram/delaunay/delaunay.h>
50 #include <geogram/basic/geometry_nd.h>
51 #include <geogram/basic/process.h>
52 #include <geogram/basic/command_line.h>
53 #include <geogram/basic/argused.h>
54 #include <geogram/basic/algorithm.h>
55 #include <geogram/bibliography/bibliography.h>
56
57 /*
58 * There are three levels of implementation:
59 * Level 1: RestrictedVoronoiDiagram is the abstract API seen from client code
60 * Level 2: RVD_Nd_Impl<DIM> implements RestrictedVoronoiDiagram
61 * Level 3: RVD_Nd_Impl<DIM>::GenRestrictedVoronoiDiagram is
62 * an instantiation of GEOGen::RestrictedVoronoiDiagram (from generic_RVD.h)
63 *
64 * Warning: there are approx. 1000 lines of boring code ahead.
65 */
66
67 namespace {
68
69 using namespace GEO;
70
71 /**
72 * \brief Generic implementation of RestrictedVoronoiDiagram.
73 * \tparam DIM dimension
74 */
75 template <unsigned int DIM>
76 class RVD_Nd_Impl : public GEO::RestrictedVoronoiDiagram {
77
78 /** \brief This class type */
79 typedef RVD_Nd_Impl<DIM> thisclass;
80
81 /** \brief The base class of this class */
82 typedef RestrictedVoronoiDiagram baseclass;
83
84 public:
85 /** \brief Implementation based on the generic version. */
86 typedef GEOGen::RestrictedVoronoiDiagram<DIM>
87 GenRestrictedVoronoiDiagram;
88
89 /** \brief Representation of points. */
90 typedef vecng<DIM, double> Point;
91
92 /** \brief Representation of vectors. */
93 typedef vecng<DIM, double> Vector;
94
95 /** \brief Represents a point and its symbolic information. */
96 typedef typename GenRestrictedVoronoiDiagram::Vertex Vertex;
97
98 /**
99 * \brief Specifies the computation done by the threads.
100 */
101 enum ThreadMode {
102 MT_NONE, /**< uninitialized */
103 MT_LLOYD, /**< Lloyd iteration */
104 MT_NEWTON, /**< Newton optimization */
105 MT_INT_SMPLX, /**< Newton with integration simplex */
106 MT_POLYG, /**< Polygon callback */
107 MT_POLYH /**< Polyhedron callback */
108 };
109
110 /**
111 * \brief Gets a mesh vertex from its index.
112 * \param[in] v index of the vertex
113 * \return a const reference to a Point
114 */
115 ✗ const Point& mesh_vertex(index_t v) {
116 ✗ return *(const Point*) mesh_->vertices.point_ptr(v);
117 }
118
119
120 /**
121 * \brief Creates a RVD_Nd_Impl.
122 *
123 * \details The dimension is determined by \p mesh->dimension().
124 * \param[in] delaunay the Delaunay triangulation
125 * \param[in] mesh the input mesh
126 * \param[in] R3_embedding gives for each vertex
127 * its mapping in 3D space.
128 * \param[in] R3_embedding_stride gives the stride between
129 * two consecutive vertices in R3_embedding
130 */
131 96 RVD_Nd_Impl(
132 Delaunay* delaunay, Mesh* mesh,
133 const double* R3_embedding, index_t R3_embedding_stride
134 ) :
135 RestrictedVoronoiDiagram(
136 delaunay, mesh, R3_embedding, R3_embedding_stride
137 ),
138
1/2
✓ Branch 2 taken 48 times.
✗ Branch 3 not taken.
 96 RVD_(delaunay, mesh) {
139 96 use_exact_projection_ = false;
140 96 is_slave_ = false;
141 96 master_ = nullptr;
142 96 has_weights_ = false;
143
4/6
✓ Branch 2 taken 48 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 48 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 1 times.
✓ Branch 9 taken 47 times.
 288 if(mesh->vertices.attributes().is_defined("weight")) {
144
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
 4 vertex_weight_.bind(mesh->vertices.attributes(), "weight");
145 2 has_weights_ = true;
146 }
147 96 parts_ = nullptr;
148 96 nb_parts_ = 0;
149 96 funcval_ = 0.0;
150 96 simplex_func_ = nullptr;
151 96 polygon_callback_ = nullptr;
152 96 polyhedron_callback_ = nullptr;
153 96 arg_vectors_ = nullptr;
154 96 arg_scalars_ = nullptr;
155 96 thread_mode_ = MT_NONE;
156 96 nb_triangles_ = 0;
157 96 }
158
159 /**
160 * \brief Constructor for parts, used in multithreading mode.
161 */
162 96 RVD_Nd_Impl() :
163 RestrictedVoronoiDiagram(nullptr, nullptr, nullptr, 0),
164
1/2
✓ Branch 2 taken 48 times.
✗ Branch 3 not taken.
 96 RVD_(nullptr, nullptr) {
165 96 use_exact_projection_ = false;
166 96 is_slave_ = true;
167 96 master_ = nullptr;
168 96 mesh_ = nullptr;
169 96 parts_ = nullptr;
170 96 nb_parts_ = 0;
171 96 facets_begin_ = NO_INDEX;
172 96 facets_end_ = NO_INDEX;
173 96 funcval_ = 0.0;
174 96 simplex_func_ = nullptr;
175 96 polygon_callback_ = nullptr;
176 96 polyhedron_callback_ = nullptr;
177 96 arg_vectors_ = nullptr;
178 96 arg_scalars_ = nullptr;
179 96 thread_mode_ = MT_NONE;
180 96 nb_triangles_ = 0;
181 96 }
182
183 96 void set_delaunay(Delaunay* delaunay) override {
184 96 baseclass::set_delaunay(delaunay);
185 96 RVD_.set_delaunay(delaunay);
186
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 48 times.
 96 for(index_t p = 0; p < nb_parts_; ++p) {
187 ✗ parts_[p].set_delaunay(delaunay);
188 }
189 96 }
190
191 350 void set_check_SR(bool x) override {
192 350 RVD_.set_check_SR(x);
193
2/2
✓ Branch 0 taken 96 times.
✓ Branch 1 taken 175 times.
 542 for(index_t p = 0; p < nb_parts_; ++p) {
194 192 parts_[p].set_check_SR(x);
195 }
196 350 }
197
198 206 void set_exact_predicates(bool x) override {
199 206 RVD_.set_exact_predicates(x);
200
2/2
✓ Branch 0 taken 16 times.
✓ Branch 1 taken 103 times.
 238 for(index_t p = 0; p < nb_parts_; ++p) {
201 32 parts_[p].set_exact_predicates(x);
202 }
203 206 }
204
205 ✗ bool exact_predicates() const override {
206 ✗ return RVD_.exact_predicates();
207 }
208
209 /********************************************************************/
210
211 /**
212 * \brief Place holder, "no locking" policy.
213 * \details NoLocks is used by algorithms templated
214 * by locking policy, for the single-threaded instances
215 * that do not need synchronization. The multi-threaded
216 * instances are parameterized by SpinLockArray.
217 */
218 class NoLocks {
219 public:
220 /**
221 * \brief Acquires a spinlock.
222 * \details Does nothing in this version
223 * \param[in] i index of the spinlock to acquire
224 */
225 ✗ void acquire_spinlock(index_t i) {
226 ✗ geo_argused(i);
227 ✗ }
228
229 /**
230 * \brief Releases a spinlock.
231 * \details Does nothing in this version
232 * \param[in] i index of the spinlock to release
233 */
234 ✗ void release_spinlock(index_t i) {
235 ✗ geo_argused(i);
236 ✗ }
237 };
238
239 /******************************************************************/
240
241 /**
242 * \brief Implementation class for surfacic Lloyd relaxation.
243 * \details To be used as a template argument
244 * to RVD::for_each_triangle().
245 * This version ignores the weights.
246 *
247 * Computes for each RVD cell:
248 * - mg[v] (v's Voronoi cell's total area times centroid)
249 * - m[v] (v's total area)
250 * \tparam LOCKS locking policy
251 * (can be one of Process::SpinLockArray, NoLocks)
252 */
253 template <class LOCKS>
254 class ComputeCentroids {
255 public:
256 /**
257 * \brief Constructs a ComputeCentroids.
258 * \param[out] mg where to store the centroids
259 * \param[out] m where to store the masses
260 * \param[in] locks the array of locks
261 * (or NoLocks in single thread mode)
262 */
263 280 ComputeCentroids(
264 double* mg,
265 double* m,
266 LOCKS& locks
267 ) :
268 280 mg_(mg),
269 280 m_(m),
270 280 locks_(locks) {
271 280 }
272
273 /**
274 * \brief The callback called for each integration simplex.
275 * \param[in] v index of current center vertex
276 * \param[in] p1 first vertex of current integration simplex
277 * \param[in] p2 second vertex of current integration simplex
278 * \param[in] p3 third vertex of current integration simplex
279 */
280 4523982 void operator() (
281 index_t v,
282 const double* p1,
283 const double* p2,
284 const double* p3
285 ) const {
286 4523982 double cur_m = Geom::triangle_area(p1, p2, p3, DIM);
287 4523982 double s = cur_m / 3.0;
288 4523982 locks_.acquire_spinlock(v);
289 4523982 m_[v] += cur_m;
290 4523982 double* cur_mg_out = mg_ + v * DIM;
291
2/2
✓ Branch 0 taken 8987061 times.
✓ Branch 1 taken 2261991 times.
 22498104 for(coord_index_t coord = 0; coord < DIM; coord++) {
292 17974122 cur_mg_out[coord] +=
293 17974122 s * (p1[coord] + p2[coord] + p3[coord]);
294 }
295 4523982 locks_.release_spinlock(v);
296 4523982 }
297
298 private:
299 double* mg_;
300 double* m_;
301 LOCKS& locks_;
302 };
303
304 /**
305 * \brief Implementation class for surfacic Lloyd relaxation.
306 * \details To be used as a template
307 * argument to RVD::for_each_triangle().
308 * This version takes the weights into account.
309 *
310 * Computes for each RVD cell:
311 * - mg[v] (v's Voronoi cell's total area times centroid)
312 * - m[v] (v's total area)
313 * \tparam LOCKS locking policy
314 * (can be one of Process::SpinLockArray, NoLocks)
315 */
316 template <class LOCKS>
317 class ComputeCentroidsWeighted {
318 public:
319 /**
320 * \brief Constructs a ComputeCentroidsWeighted.
321 * \param[out] mg where to store the centroids
322 * \param[out] m where to store the masses
323 * \param[in] locks the array of locks
324 * (or NoLocks in single thread mode)
325 */
326 40 ComputeCentroidsWeighted(
327 double* mg,
328 double* m,
329 LOCKS& locks
330 ) :
331 40 mg_(mg),
332 40 m_(m),
333 40 locks_(locks) {
334 40 }
335
336 /**
337 * \brief The callback called for each integration simplex.
338 * \param[in] v index of current center vertex
339 * \param[in] v1 first vertex of current integration simplex
340 * \param[in] v2 second vertex of current integration simplex
341 * \param[in] v3 third vertex of current integration simplex
342 */
343 1578712 void operator() (
344 index_t v,
345 const Vertex& v1,
346 const Vertex& v2,
347 const Vertex& v3
348 ) const {
349 double cur_m;
350 double cur_Vg[DIM];
351
1/2
✓ Branch 7 taken 789356 times.
✗ Branch 8 not taken.
 1578712 Geom::triangle_centroid(
352 v1.point(), v2.point(), v3.point(),
353 v1.weight(), v2.weight(), v3.weight(),
354 cur_Vg, cur_m, DIM
355 );
356
1/2
✓ Branch 1 taken 789356 times.
✗ Branch 2 not taken.
 1578712 locks_.acquire_spinlock(v);
357 1578712 m_[v] += cur_m;
358 1578712 double* cur_mg_out = mg_ + v * DIM;
359
2/2
✓ Branch 0 taken 4736136 times.
✓ Branch 1 taken 789356 times.
 11050984 for(coord_index_t coord = 0; coord < DIM; coord++) {
360 9472272 cur_mg_out[coord] += cur_Vg[coord];
361 }
362
1/2
✓ Branch 1 taken 789356 times.
✗ Branch 2 not taken.
 1578712 locks_.release_spinlock(v);
363 1578712 }
364
365 private:
366 double* mg_;
367 double* m_;
368 LOCKS& locks_;
369 };
370
371 400 void compute_centroids_on_surface(double* mg, double* m) override {
372 400 create_threads();
373
2/2
✓ Branch 1 taken 160 times.
✓ Branch 2 taken 40 times.
 400 if(nb_parts() == 0) {
374
1/2
✓ Branch 0 taken 160 times.
✗ Branch 1 not taken.
 320 if(master_ != nullptr) {
375
2/2
✓ Branch 0 taken 20 times.
✓ Branch 1 taken 140 times.
 320 if(has_weights_) {
376
1/2
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
 40 RVD_.for_each_triangle(
377 40 ComputeCentroidsWeighted<Process::SpinLockArray>(
378 40 mg, m, master_->spinlocks_
379 )
380 );
381 } else {
382
1/2
✓ Branch 1 taken 140 times.
✗ Branch 2 not taken.
 280 RVD_.for_each_triangle(
383 280 ComputeCentroids<Process::SpinLockArray>(
384 280 mg, m, master_->spinlocks_
385 )
386 );
387 }
388 } else {
389 NoLocks nolocks;
390 ✗ if(has_weights_) {
391 ✗ RVD_.for_each_triangle(
392 ✗ ComputeCentroidsWeighted<NoLocks>(
393 mg, m, nolocks
394 )
395 );
396 } else {
397 ✗ RVD_.for_each_triangle(
398 ✗ ComputeCentroids<NoLocks>(mg, m, nolocks)
399 );
400 }
401 }
402 } else {
403 80 thread_mode_ = MT_LLOYD;
404 80 arg_vectors_ = mg;
405 80 arg_scalars_ = m;
406 80 spinlocks_.resize(delaunay_->nb_vertices());
407
1/2
✓ Branch 2 taken 40 times.
✗ Branch 3 not taken.
 80 parallel_for(
408 0, nb_parts(),
409 320 [this](index_t i) { run_thread(i); }
410 );
411 }
412 400 }
413
414 /********************************************************************/
415
416 /**
417 * \brief Implementation class for surfacic Lloyd relaxation.
418 * \details To be used as a template argument
419 * to RVD::for_each_volumetric_integration_simplex().
420 * This version ignores the weights.
421 *
422 * Computes for each RVD cell:
423 * - mg[v] (v's Voronoi cell's total area times centroid)
424 * - m[v] (v's total area)
425 * \tparam LOCKS locking policy
426 * (can be one of Process::SpinLockArray, NoLocks)
427 */
428 template <class LOCKS>
429 class ComputeCentroidsVolumetric {
430 public:
431 /**
432 * \brief Constructs a ComputeCentroidsVolumetric.
433 * \param[out] mg where to store the centroids
434 * \param[out] m where to store the masses
435 * \param[in] delaunay the Delaunay triangulation
436 * \param[in] locks the array of locks
437 * (or NoLocks in single thread mode)
438 */
439 160 ComputeCentroidsVolumetric(
440 double* mg,
441 double* m,
442 const Delaunay* delaunay,
443 LOCKS& locks
444 ) :
445 160 mg_(mg),
446 160 m_(m),
447 160 delaunay_(delaunay),
448 160 locks_(locks) {
449 160 }
450
451 /**
452 * \brief The callback called for each integration simplex.
453 * \param[in] v index of current center vertex
454 * \param[in] v_adj (unused here) is the index of the Voronoi cell
455 * adjacent to t accros facet (\p v1, \p v2, \p v3) or
456 * NO_INDEX if it does not exists
457 * \param[in] t (unused here) is the index of the current
458 * tetrahedron
459 * \param[in] t_adj (unused here) is the index of the
460 * tetrahedron adjacent to t accros facet (\p v1, \p v2, \p v3)
461 * or NO_INDEX if it does not exists
462 * \param[in] p0 first vertex of current integration simplex
463 * \param[in] p1 second vertex of current integration simplex
464 * \param[in] p2 third vertex of current integration simplex
465 * \param[in] p3 fourth vertex of current integration simplex
466 */
467 1506142 void operator() (
468 index_t v, index_t v_adj,
469 index_t t, index_t t_adj,
470 const double* p0,
471 const double* p1,
472 const double* p2,
473 const double* p3
474 ) const {
475 1506142 geo_argused(v_adj);
476 1506142 geo_argused(t);
477 1506142 geo_argused(t_adj);
478 1506142 double cur_m = Geom::tetra_volume<DIM>(
479 p0, p1, p2, p3
480 );
481 1506142 double s = cur_m / 4.0;
482 1506142 locks_.acquire_spinlock(v);
483 1506142 m_[v] += cur_m;
484 1506142 double* cur_mg_out = mg_ + v * DIM;
485
2/2
✓ Branch 0 taken 4158741 times.
✓ Branch 1 taken 753071 times.
 9823624 for(coord_index_t coord = 0; coord < DIM; coord++) {
486 8317482 cur_mg_out[coord] += s * (
487 8317482 p0[coord] + p1[coord] + p2[coord] + p3[coord]
488 );
489 }
490 1506142 locks_.release_spinlock(v);
491 1506142 }
492
493 private:
494 double* mg_;
495 double* m_;
496 const Delaunay* delaunay_;
497 LOCKS& locks_;
498 };
499
500 200 void compute_centroids_in_volume(double* mg, double* m) override {
501 200 create_threads();
502
2/2
✓ Branch 1 taken 80 times.
✓ Branch 2 taken 20 times.
 200 if(nb_parts() == 0) {
503
1/2
✓ Branch 0 taken 80 times.
✗ Branch 1 not taken.
 160 if(master_ != nullptr) {
504
1/2
✓ Branch 1 taken 80 times.
✗ Branch 2 not taken.
 160 RVD_.for_each_tetrahedron(
505 320 ComputeCentroidsVolumetric<Process::SpinLockArray>(
506 160 mg, m, RVD_.delaunay(), master_->spinlocks_
507 )
508 );
509 } else {
510 NoLocks nolocks;
511 ✗ RVD_.for_each_tetrahedron(
512 ✗ ComputeCentroidsVolumetric<NoLocks>(
513 ✗ mg, m, RVD_.delaunay(), nolocks
514 )
515 );
516 }
517 } else {
518 40 thread_mode_ = MT_LLOYD;
519 40 arg_vectors_ = mg;
520 40 arg_scalars_ = m;
521 40 spinlocks_.resize(delaunay_->nb_vertices());
522
1/2
✓ Branch 2 taken 20 times.
✗ Branch 3 not taken.
 40 parallel_for(
523 0, nb_parts(),
524 160 [this](index_t i) { run_thread(i); }
525 );
526 }
527 200 }
528
529 /********************************************************************/
530
531 /**
532 * \brief Implementation class for Newton-based restricted CVT.
533 * \details To be used as a template argument
534 * to RVD::for_each_triangle().
535 * This version ignores the weights.
536 *
537 * Computes for each RVD cell:
538 * - g (gradient)
539 * - f (CVT energy)
540 * \tparam LOCKS locking policy
541 * (can be one of Process::SpinLockArray, NoLocks)
542 */
543 template <class LOCKS>
544 class ComputeCVTFuncGrad {
545 public:
546 /**
547 * \brief Constructs a ComputeCVTFuncGrad.
548 * \param[in] RVD the restricted Voronoi diagram
549 * \param[out] f the computed function value
550 * \param[out] g the computed gradient of f,
551 * allocated by caller, and managed
552 * by caller
553 * \param[in] locks the array of locks
554 * (or NoLocks in single thread mode)
555 */
556 1688 ComputeCVTFuncGrad(
557 const GenRestrictedVoronoiDiagram& RVD,
558 double& f,
559 double* g,
560 LOCKS& locks
561 ) :
562 1688 f_(f),
563 1688 g_(g),
564 1688 locks_(locks),
565 1688 RVD_(RVD) {
566 1688 }
567
568 /**
569 * \brief The callback called for each integration simplex.
570 * \param[in] v index of current center vertex
571 * \param[in] p1 first vertex of current integration simplex
572 * \param[in] p2 second vertex of current integration simplex
573 * \param[in] p3 third vertex of current integration simplex
574 */
575 19718710 void operator() (
576 index_t v,
577 const double* p1,
578 const double* p2,
579 const double* p3
580 ) const {
581
582 19718710 const double* p0 = RVD_.delaunay()->vertex_ptr(v);
583
584 19718710 double t_area = Geom::triangle_area(p1, p2, p3, DIM);
585
586 19718710 double cur_f = 0.0;
587
2/2
✓ Branch 0 taken 43964100 times.
✓ Branch 1 taken 9859355 times.
 107646910 for(index_t c = 0; c < DIM; c++) {
588 87928200 double u0 = p0[c] - p1[c];
589 87928200 double u1 = p0[c] - p2[c];
590 87928200 double u2 = p0[c] - p3[c];
591 87928200 cur_f += u0 * u0;
592 87928200 cur_f += u1 * (u0 + u1);
593 87928200 cur_f += u2 * (u0 + u1 + u2);
594 }
595
596 19718710 f_ += t_area * cur_f / 6.0;
597
598 19718710 locks_.acquire_spinlock(v);
599
2/2
✓ Branch 0 taken 43964100 times.
✓ Branch 1 taken 9859355 times.
 107646910 for(index_t c = 0; c < DIM; c++) {
600 87928200 double Gc = (1.0 / 3.0) * (p1[c] + p2[c] + p3[c]);
601 87928200 g_[DIM * v + c] += (2.0 * t_area) * (p0[c] - Gc);
602 }
603 19718710 locks_.release_spinlock(v);
604 19718710 }
605
606 double& f_;
607 double* g_;
608 LOCKS& locks_;
609 const GenRestrictedVoronoiDiagram& RVD_;
610 };
611
612 /**
613 * \brief Implementation class for Newton-based restricted CVT.
614 * \details To be used as a template argument
615 * to RVD::for_each_triangle().
616 * This version takes the weights into account.
617 *
618 * Computes for each RVD cell:
619 * - g (gradient)
620 * - f (CVT energy)
621 * \tparam LOCKS locking policy
622 * (can be one of Process::SpinLockArray, NoLocks)
623 */
624 template <class LOCKS>
625 class ComputeCVTFuncGradWeighted {
626 public:
627 /**
628 * \brief Constructs a ComputeCVTFuncGradWeighted.
629 * \param[in] RVD the restricted Voronoi diagram
630 * \param[out] f the computed function value
631 * \param[out] g the computed gradient of f,
632 * allocated by caller, and managed by caller
633 * \param[in] locks the array of locks
634 * (or NoLocks in single thread mode)
635 */
636 280 ComputeCVTFuncGradWeighted(
637 const GenRestrictedVoronoiDiagram& RVD,
638 double& f,
639 double* g,
640 LOCKS& locks
641 ) :
642 280 f_(f),
643 280 g_(g),
644 280 locks_(locks),
645 280 RVD_(RVD) {
646 280 }
647
648 /**
649 * \brief The callback called for each integration simplex.
650 * \param[in] v index of current center vertex
651 * \param[in] v1 first vertex of current integration simplex
652 * \param[in] v2 second vertex of current integration simplex
653 * \param[in] v3 third vertex of current integration simplex
654 */
655 11020392 void operator() (
656 index_t v,
657 const Vertex& v1,
658 const Vertex& v2,
659 const Vertex& v3
660 ) const {
661
662
1/2
✓ Branch 2 taken 5510196 times.
✗ Branch 3 not taken.
 11020392 const double* p0 = RVD_.delaunay()->vertex_ptr(v);
663
664 11020392 const double* p1 = v1.point();
665 11020392 const double* p2 = v2.point();
666 11020392 const double* p3 = v3.point();
667
668
1/2
✓ Branch 1 taken 5510196 times.
✗ Branch 2 not taken.
 11020392 double t_area = Geom::triangle_area(p1, p2, p3, DIM);
669
670 11020392 double Sp = v1.weight() + v2.weight() + v3.weight();
671 double rho[3], alpha[3];
672 11020392 rho[0] = v1.weight();
673 11020392 rho[1] = v2.weight();
674 11020392 rho[2] = v3.weight();
675 11020392 alpha[0] = Sp + rho[0];
676 11020392 alpha[1] = Sp + rho[1];
677 11020392 alpha[2] = Sp + rho[2];
678
679 11020392 double dotprod_00 = 0.0;
680 11020392 double dotprod_10 = 0.0;
681 11020392 double dotprod_11 = 0.0;
682 11020392 double dotprod_20 = 0.0;
683 11020392 double dotprod_21 = 0.0;
684 11020392 double dotprod_22 = 0.0;
685
2/2
✓ Branch 0 taken 33061176 times.
✓ Branch 1 taken 5510196 times.
 77142744 for(unsigned int c = 0; c < DIM; c++) {
686 66122352 double sp0 = p0[c] - p1[c];
687 66122352 double sp1 = p0[c] - p2[c];
688 66122352 double sp2 = p0[c] - p3[c];
689 66122352 dotprod_00 += sp0 * sp0;
690 66122352 dotprod_10 += sp1 * sp0;
691 66122352 dotprod_11 += sp1 * sp1;
692 66122352 dotprod_20 += sp2 * sp0;
693 66122352 dotprod_21 += sp2 * sp1;
694 66122352 dotprod_22 += sp2 * sp2;
695 }
696
697 11020392 double cur_f = 0.0;
698 11020392 cur_f += (alpha[0] + rho[0]) * dotprod_00; // 0 0
699 11020392 cur_f += (alpha[1] + rho[0]) * dotprod_10; // 1 0
700 11020392 cur_f += (alpha[1] + rho[1]) * dotprod_11; // 1 1
701 11020392 cur_f += (alpha[2] + rho[0]) * dotprod_20; // 2 0
702 11020392 cur_f += (alpha[2] + rho[1]) * dotprod_21; // 2 1
703 11020392 cur_f += (alpha[2] + rho[2]) * dotprod_22; // 2 2
704
705 11020392 f_ += t_area * cur_f / 30.0;
706 11020392 double* g_out = g_ + v * DIM;
707
1/2
✓ Branch 1 taken 5510196 times.
✗ Branch 2 not taken.
 11020392 locks_.acquire_spinlock(v);
708
2/2
✓ Branch 0 taken 33061176 times.
✓ Branch 1 taken 5510196 times.
 77142744 for(index_t c = 0; c < DIM; c++) {
709 66122352 g_out[c] += (t_area / 6.0) * (
710 66122352 4.0 * Sp * p0[c] - (
711 66122352 alpha[0] * p1[c] +
712 66122352 alpha[1] * p2[c] +
713 66122352 alpha[2] * p3[c]
714 )
715 );
716 }
717
1/2
✓ Branch 1 taken 5510196 times.
✗ Branch 2 not taken.
 11020392 locks_.release_spinlock(v);
718 11020392 }
719
720 double& f_;
721 double* g_;
722 LOCKS& locks_;
723 const GenRestrictedVoronoiDiagram& RVD_;
724 };
725
726 2460 void compute_CVT_func_grad_on_surface(double& f, double* g) override {
727 2460 create_threads();
728
2/2
✓ Branch 1 taken 984 times.
✓ Branch 2 taken 246 times.
 2460 if(nb_parts() == 0) {
729
1/2
✓ Branch 0 taken 984 times.
✗ Branch 1 not taken.
 1968 if(master_ != nullptr) {
730
2/2
✓ Branch 0 taken 140 times.
✓ Branch 1 taken 844 times.
 1968 if(has_weights_) {
731
1/2
✓ Branch 1 taken 140 times.
✗ Branch 2 not taken.
 280 RVD_.for_each_triangle(
732 280 ComputeCVTFuncGradWeighted<Process::SpinLockArray>(
733 280 RVD_, f, g, master_->spinlocks_
734 )
735 );
736 } else {
737
1/2
✓ Branch 1 taken 844 times.
✗ Branch 2 not taken.
 1688 RVD_.for_each_triangle(
738 1688 ComputeCVTFuncGrad<Process::SpinLockArray>(
739 1688 RVD_, f, g, master_->spinlocks_
740 )
741 );
742 }
743 } else {
744 NoLocks nolocks;
745 ✗ if(has_weights_) {
746 ✗ RVD_.for_each_triangle(
747 ✗ ComputeCVTFuncGradWeighted<NoLocks>(
748 ✗ RVD_, f, g, nolocks
749 )
750 );
751 } else {
752 ✗ RVD_.for_each_triangle(
753 ✗ ComputeCVTFuncGrad<NoLocks>(
754 ✗ RVD_, f, g, nolocks
755 )
756 );
757 }
758 }
759 } else {
760 492 thread_mode_ = MT_NEWTON;
761 492 arg_vectors_ = g;
762 492 spinlocks_.resize(delaunay_->nb_vertices());
763
2/2
✓ Branch 1 taken 984 times.
✓ Branch 2 taken 246 times.
 2460 for(index_t t = 0; t < nb_parts(); t++) {
764 1968 part(t).funcval_ = 0.0;
765 }
766
1/2
✓ Branch 2 taken 246 times.
✗ Branch 3 not taken.
 492 parallel_for(
767 0, nb_parts(),
768 1968 [this](index_t i) { run_thread(i); }
769 );
770
2/2
✓ Branch 1 taken 984 times.
✓ Branch 2 taken 246 times.
 2460 for(index_t t = 0; t < nb_parts(); t++) {
771 1968 f += part(t).funcval_;
772 }
773 }
774 2460 }
775
776 /********************************************************************/
777
778 /**
779 * \brief Implementation class for Newton-based restricted CVT
780 * in volume.
781 * \details To be used as a template argument
782 * to RVD::for_each_volumetric_integration_simplex().
783 * This version ignores the weights.
784 *
785 * Computes for each RVD cell:
786 * - g (gradient)
787 * - f (CVT energy)
788 * \tparam LOCKS locking policy
789 * (can be one of Process::SpinLockArray, NoLocks)
790 */
791 template <class LOCKS>
792 class ComputeCVTFuncGradVolumetric {
793 public:
794 /**
795 * \brief Constructs a ComputeCentroidsFuncGradVolumetric.
796 * \param[in] RVD the restricted Voronoi diagram
797 * \param[out] f the computed function value
798 * \param[out] g the computed gradient of f,
799 * allocated by caller, and managed
800 * by caller
801 * \param[in] locks the array of locks
802 * (or NoLocks in single thread mode)
803 */
804 1096 ComputeCVTFuncGradVolumetric(
805 const GenRestrictedVoronoiDiagram& RVD,
806 double& f,
807 double* g,
808 LOCKS& locks
809 ) :
810 1096 f_(f),
811 1096 g_(g),
812 1096 locks_(locks),
813 1096 RVD_(RVD) {
814 1096 }
815
816 /**
817 * \brief The callback called for each integration simplex.
818 * \param[in] v index of current center vertex
819 * \param[in] v_adj (unused here) is the index of the Voronoi cell
820 * adjacent to t accros facet (\p v1, \p v2, \p v3) or
821 * NO_INDEX if it does not exists
822 * \param[in] t (unused here) is the index of the current
823 * tetrahedron
824 * \param[in] t_adj (unused here) is the index of the
825 * tetrahedron adjacent to t accros facet (\p v1, \p v2, \p v3)
826 * or NO_INDEX if it does not exists
827 * \param[in] p1 first vertex of current integration simplex
828 * \param[in] p2 second vertex of current integration simplex
829 * \param[in] p3 third vertex of current integration simplex
830 */
831 10369872 void operator() (
832 index_t v,
833 index_t v_adj,
834 index_t t,
835 index_t t_adj,
836 const double* p1,
837 const double* p2,
838 const double* p3
839 ) const {
840 10369872 geo_argused(v_adj);
841 10369872 geo_argused(t);
842 10369872 geo_argused(t_adj);
843 10369872 const double* p0 = RVD_.delaunay()->vertex_ptr(v);
844
845 10369872 double mi = Geom::tetra_volume<DIM>(p0, p1, p2, p3);
846
847 // fi = (mi/10)*(U.U + V.V + W.W + U.V + V.W + W.U)
848 // where: U = p1-p0 ; V = p2-p0 and W=p3-p0
849 10369872 double fi = 0.0;
850
2/2
✓ Branch 0 taken 28237498 times.
✓ Branch 1 taken 5184936 times.
 66844868 for(coord_index_t c = 0; c < DIM; ++c) {
851 56474996 double Uc = p1[c] - p0[c];
852 56474996 double Vc = p2[c] - p0[c];
853 56474996 double Wc = p3[c] - p0[c];
854 56474996 fi += geo_sqr(Uc) + geo_sqr(Vc) + geo_sqr(Wc);
855 56474996 fi += (Uc * Vc + Vc * Wc + Wc * Uc);
856 }
857 10369872 fi *= (mi / 10.0);
858 10369872 f_ += fi;
859
860 // gi = 2*mi(p0 - 1/4(p0 + p1 + p2 + p3))
861 10369872 double* g_out = g_ + v * DIM;
862 10369872 locks_.acquire_spinlock(v);
863
2/2
✓ Branch 0 taken 28237498 times.
✓ Branch 1 taken 5184936 times.
 66844868 for(coord_index_t c = 0; c < DIM; ++c) {
864 56474996 g_out[c] += 2.0 * mi * (
865 56474996 0.75 * p0[c]
866 56474996 - 0.25 * p1[c] - 0.25 * p2[c] - 0.25 * p3[c]
867 );
868 }
869 10369872 locks_.release_spinlock(v);
870 10369872 }
871
872 double& f_;
873 double* g_;
874 LOCKS& locks_;
875 const GenRestrictedVoronoiDiagram& RVD_;
876 };
877
878 1370 void compute_CVT_func_grad_in_volume(double& f, double* g) override {
879 1370 create_threads();
880
2/2
✓ Branch 1 taken 548 times.
✓ Branch 2 taken 137 times.
 1370 if(nb_parts() == 0) {
881
1/2
✓ Branch 0 taken 548 times.
✗ Branch 1 not taken.
 1096 if(master_ != nullptr) {
882
1/2
✓ Branch 1 taken 548 times.
✗ Branch 2 not taken.
 1096 RVD_.for_each_volumetric_integration_simplex(
883 1096 ComputeCVTFuncGradVolumetric<Process::SpinLockArray>(
884 1096 RVD_, f, g, master_->spinlocks_
885 )
886 );
887 } else {
888 NoLocks nolocks;
889 ✗ RVD_.for_each_volumetric_integration_simplex(
890 ✗ ComputeCVTFuncGradVolumetric<NoLocks>(
891 ✗ RVD_, f, g, nolocks
892 )
893 );
894 }
895 } else {
896 274 thread_mode_ = MT_NEWTON;
897 274 arg_vectors_ = g;
898 274 spinlocks_.resize(delaunay_->nb_vertices());
899
2/2
✓ Branch 1 taken 548 times.
✓ Branch 2 taken 137 times.
 1370 for(index_t t = 0; t < nb_parts(); t++) {
900 1096 part(t).funcval_ = 0.0;
901 }
902
1/2
✓ Branch 2 taken 137 times.
✗ Branch 3 not taken.
 274 parallel_for(
903 0, nb_parts(),
904 1096 [this](index_t i) { run_thread(i); }
905 );
906
2/2
✓ Branch 1 taken 548 times.
✓ Branch 2 taken 137 times.
 1370 for(index_t t = 0; t < nb_parts(); t++) {
907 1096 f += part(t).funcval_;
908 }
909 }
910 1370 }
911
912 /********************************************************************/
913
914 /**
915 * \brief Implementation class for computing function integrals
916 * and gradients over integration simplices.
917 * \details To be used as a template argument
918 * to RVD::for_each_triangle() and
919 * RVD::for_each_volumetric_integration_simplex()
920 */
921 class ComputeCVTFuncGradIntegrationSimplex {
922 public:
923 /**
924 * \brief Constructs a ComputeCVTFuncGradIntegrationSimplex.
925 * \param[in] RVD the Restricted Voronoi Diagram
926 * \param[in] F the IntegrationSimplex
927 */
928 ✗ ComputeCVTFuncGradIntegrationSimplex(
929 const GenRestrictedVoronoiDiagram& RVD,
930 IntegrationSimplex* F
931 ) :
932 ✗ f_(0.0),
933 ✗ RVD_(RVD),
934 ✗ simplex_func_(F) {
935 ✗ simplex_func_->reset_thread_local_storage();
936 ✗ }
937
938 /**
939 * \brief The callback called for each surfacic integration simplex.
940 * \param[in] i index of current center vertex
941 * \param[in] v1 first vertex of current integration simplex
942 * \param[in] v2 second vertex of current integration simplex
943 * \param[in] v3 third vertex of current integration simplex
944 */
945 ✗ void operator() (
946 index_t i,
947 const Vertex& v1,
948 const Vertex& v2,
949 const Vertex& v3
950 ) {
951 ✗ f_ += simplex_func_->eval(
952 ✗ i,v1,v2,v3,RVD_.current_facet()
953 );
954 ✗ }
955
956 /**
957 * \brief The callback called for each volumetric
958 * integration simplex.
959 * \param[in] v index of current center vertex
960 * \param[in] v_adj index of the Voronoi cell adjacent to t accros
961 * facet (\p v1, \p v2, \p v3) or NO_INDEX if it does not exists
962 * \param[in] t index of the current tetrahedron
963 * \param[in] t_adj index of the tetrahedron adjacent to t accros
964 * facet (\p v1, \p v2, \p v3) or NO_INDEX if it does not exists
965 * \param[in] v1 first vertex of current integration simplex
966 * \param[in] v2 second vertex of current integration simplex
967 * \param[in] v3 third vertex of current integration simplex
968 */
969 ✗ void operator() (
970 index_t v,
971 index_t v_adj,
972 index_t t,
973 index_t t_adj,
974 const Vertex& v1,
975 const Vertex& v2,
976 const Vertex& v3
977 ) {
978 ✗ geo_argused(v_adj);
979 ✗ geo_argused(t_adj);
980 ✗ f_ += simplex_func_->eval(v,v1,v2,v3,t,t_adj,v_adj);
981 ✗ }
982
983 /**
984 * \brief Gets the function value.
985 * \return The function value accumulated so far.
986 */
987 ✗ double f() const {
988 ✗ return f_;
989 }
990
991 private:
992 double f_;
993 const GenRestrictedVoronoiDiagram& RVD_;
994 IntegrationSimplex* simplex_func_;
995 };
996
997 ✗ void compute_integration_simplex_func_grad(
998 double& f, double* g, IntegrationSimplex* F
999 ) override {
1000 ✗ create_threads();
1001 ✗ if(nb_parts() == 0) {
1002 ✗ if(master_ == nullptr) {
1003 ✗ F->set_points_and_gradient(
1004 ✗ delaunay()->dimension(),
1005 delaunay()->nb_vertices(),
1006 delaunay()->vertex_ptr(0),
1007 g
1008 );
1009 }
1010 ✗ ComputeCVTFuncGradIntegrationSimplex C(RVD_,F);
1011 ✗ bool sym = RVD_.symbolic();
1012 ✗ RVD_.set_symbolic(true);
1013 ✗ if(F->volumetric()) {
1014 ✗ RVD_.for_each_volumetric_integration_simplex(
1015 C,
1016 ✗ F->background_mesh_has_varying_attribute(),
1017 false /* Coherent triangles */
1018 );
1019 } else {
1020 ✗ RVD_.for_each_triangle(C);
1021 }
1022 ✗ RVD_.set_symbolic(sym);
1023 ✗ funcval_ = C.f();
1024 ✗ f = C.f();
1025 } else {
1026 ✗ spinlocks_.resize(delaunay_->nb_vertices());
1027 ✗ F->set_points_and_gradient(
1028 ✗ delaunay()->dimension(),
1029 delaunay()->nb_vertices(),
1030 delaunay()->vertex_ptr(0),
1031 g,
1032 &spinlocks_
1033 );
1034 ✗ thread_mode_ = MT_INT_SMPLX;
1035 ✗ arg_vectors_ = g;
1036 ✗ simplex_func_ = F;
1037 ✗ funcval_ = 0.0;
1038 ✗ for(index_t t = 0; t < nb_parts(); t++) {
1039 ✗ part(t).arg_vectors_ = g;
1040 ✗ part(t).simplex_func_ = F;
1041 ✗ part(t).funcval_ = 0.0;
1042 }
1043
1044 ✗ parallel_for(
1045 0, nb_parts(),
1046 ✗ [this](index_t i) { run_thread(i); }
1047 );
1048
1049 ✗ f = 0.0;
1050 ✗ for(index_t t = 0; t < nb_parts(); t++) {
1051 ✗ f += part(t).funcval_;
1052 }
1053 }
1054 ✗ }
1055
1056 /********************************************************************/
1057
1058 /**
1059 * \brief Adapter class used internally to implement for_each_polygon()
1060 * \details Gets the current triangle from the RVD and passes it back
1061 * to the callback. It is needed because GenericRVD::for_each_polygon()
1062 * does not pass the current triangle.
1063 */
1064 // TODO: pass it through all the callbacks, because it is ridiculous:
1065 // we pass it through the first levels, then throw it, then retrieve it
1066 // (see GenRVD)
1067 class PolygonCallbackAction {
1068 public:
1069 /**
1070 * \brief PolygonCallbackAction constructor
1071 * \param[in] RVD a pointer to the restricted Voronoi diagram
1072 * \param[in] callback a pointer to the PolygonCallback
1073 */
1074 ✗ PolygonCallbackAction(
1075 GenRestrictedVoronoiDiagram& RVD,
1076 GEO::RVDPolygonCallback& callback
1077 ) :
1078 ✗ RVD_(RVD),
1079 ✗ callback_(callback) {
1080 ✗ }
1081
1082 /**
1083 * \brief Callback called for each polygon.
1084 * \details Routes the callback to the wrapped user action class.
1085 * \param[in] v index of current Delaunay seed
1086 * \param[in] P intersection between current mesh facet
1087 * and the Voronoi cell of \p v
1088 */
1089 ✗ void operator() (
1090 index_t v,
1091 const GEOGen::Polygon& P
1092 ) const {
1093 ✗ callback_(v, RVD_.current_facet(), P);
1094 ✗ }
1095
1096 protected:
1097 GenRestrictedVoronoiDiagram& RVD_;
1098 GEO::RVDPolygonCallback& callback_;
1099 };
1100
1101
1102 ✗ virtual void compute_with_polygon_callback(
1103 GEO::RVDPolygonCallback& polygon_callback
1104 ) {
1105 ✗ create_threads();
1106 ✗ if(nb_parts() == 0) {
1107 ✗ PolygonCallbackAction action(RVD_,polygon_callback);
1108 ✗ RVD_.for_each_polygon(action);
1109 } else {
1110 ✗ for(index_t t = 0; t < nb_parts(); t++) {
1111 ✗ part(t).RVD_.set_symbolic(RVD_.symbolic());
1112 ✗ part(t).RVD_.set_connected_components_priority(
1113 ✗ RVD_.connected_components_priority()
1114 );
1115 }
1116 ✗ spinlocks_.resize(delaunay_->nb_vertices());
1117 ✗ thread_mode_ = MT_POLYG;
1118 ✗ polygon_callback_ = &polygon_callback;
1119 ✗ polygon_callback_->set_spinlocks(&spinlocks_);
1120 // Note: callback begin()/end() is called in for_each_polygon()
1121 ✗ parallel_for(
1122 0, nb_parts(),
1123 ✗ [this](index_t i) { run_thread(i); }
1124 );
1125 ✗ polygon_callback_->set_spinlocks(nullptr);
1126 }
1127 ✗ }
1128
1129 ✗ virtual void compute_with_polyhedron_callback(
1130 GEO::RVDPolyhedronCallback& polyhedron_callback
1131 ) {
1132 ✗ create_threads();
1133 ✗ if(nb_parts() == 0) {
1134 ✗ RVD_.for_each_polyhedron(polyhedron_callback);
1135 } else {
1136 ✗ for(index_t t = 0; t < nb_parts(); t++) {
1137 ✗ part(t).RVD_.set_symbolic(RVD_.symbolic());
1138 ✗ part(t).RVD_.set_connected_components_priority(
1139 ✗ RVD_.connected_components_priority()
1140 );
1141 }
1142 ✗ spinlocks_.resize(delaunay_->nb_vertices());
1143 ✗ thread_mode_ = MT_POLYH;
1144 ✗ polyhedron_callback_ = &polyhedron_callback;
1145 ✗ polyhedron_callback_->set_spinlocks(&spinlocks_);
1146 // Note: callback begin()/end() is
1147 // called in for_each_polyhedron()
1148 ✗ parallel_for(
1149 0, nb_parts(),
1150 ✗ [this](index_t i) { run_thread(i); }
1151 );
1152 ✗ polyhedron_callback_->set_spinlocks(nullptr);
1153 }
1154 ✗ }
1155
1156 /********************************************************************/
1157
1158 /**
1159 * \brief Implementation class for explicitly constructing
1160 * a surfacic mesh that corresponds to the surfacic
1161 * restricted Voronoi diagram.
1162 * \details To be used as a template argument
1163 * to RVD::for_each_polygon(). The current Vornoi cell is
1164 * reported in facet region.
1165 * \tparam BUILDER a class that implements iterative mesh building,
1166 * e.g., MeshBuilder.
1167 */
1168 template <class BUILDER>
1169 class BuildRVD {
1170 public:
1171 /**
1172 * \brief Constructs a new BuildRVD.
1173 * \param[in] RVD_in the restricted Voronoi diagram
1174 * \param[in] builder the lesh builder
1175 */
1176 72 BuildRVD(
1177 const GenRestrictedVoronoiDiagram& RVD_in,
1178 BUILDER& builder
1179 ) :
1180 72 RVD(RVD_in),
1181 72 builder_(builder),
1182 72 current_facet_(NO_INDEX) {
1183 72 builder_.begin_surface();
1184 72 }
1185
1186 /**
1187 * \brief The destructor
1188 * \details Terminates the current facet
1189 * and the current surface.
1190 */
1191 72 ~BuildRVD() {
1192
1/2
✓ Branch 0 taken 36 times.
✗ Branch 1 not taken.
 72 if(current_facet_ != NO_INDEX) {
1193 72 builder_.end_reference_facet();
1194 }
1195 72 builder_.end_surface();
1196 72 }
1197
1198 /**
1199 * \brief The callback called for each restricted Voronoi cell.
1200 * \param[in] v index of current center vertex
1201 * \param[in] P current restricted Voronoi cell
1202 */
1203 420818 void operator() (
1204 index_t v,
1205 const typename GenRestrictedVoronoiDiagram::Polygon& P
1206 ) {
1207 420818 index_t f = RVD.current_facet();
1208
2/2
✓ Branch 0 taken 60502 times.
✓ Branch 1 taken 149907 times.
 420818 if(f != current_facet_) {
1209
2/2
✓ Branch 0 taken 60466 times.
✓ Branch 1 taken 36 times.
 121004 if(current_facet_ != NO_INDEX) {
1210 120932 builder_.end_reference_facet();
1211 }
1212 121004 current_facet_ = f;
1213 121004 builder_.begin_reference_facet(f);
1214 }
1215 420818 builder_.begin_facet(v);
1216
2/2
✓ Branch 1 taken 841011 times.
✓ Branch 2 taken 210409 times.
 2102840 for(index_t i = 0; i < P.nb_vertices(); i++) {
1217 1682022 const Vertex& ve = P.vertex(i);
1218 1682022 builder_.add_vertex_to_facet(ve.point(), ve.sym());
1219 }
1220 420818 builder_.end_facet();
1221 420818 }
1222
1223 private:
1224 const GenRestrictedVoronoiDiagram& RVD;
1225 BUILDER& builder_;
1226 index_t current_facet_;
1227 };
1228
1229 /**
1230 * \brief Implementation class for explicitly constructing
1231 * a volumetric mesh that corresponds to the volumetric
1232 * restricted Voronoi diagram.
1233 * \details To be used as a template argument
1234 * to RVD::for_each_volumetric_integration_simplex().
1235 * The current Voronoi cell is reported in tetrahedron region.
1236 * \note For the moment, vertices are duplicated (will be fixed
1237 * in a future version).
1238 */
1239 class BuildVolumetricRVD {
1240 public:
1241 /**
1242 * Constructs a new BuildVolumetricRVD.
1243 * \param[in] RVD the volumetric restricted Voronoi diagram
1244 * \param[in] dim dimension of the points (can be smaller
1245 * than actual dimension of the RVD).
1246 * \param[out] vertices coordinates of the generated vertices
1247 * \param[out] triangle_vertex_indices generated triangles, as
1248 * vertex indices triplets
1249 * \param[out] tet_vertex_indices generated tetrahedra, as
1250 * vertex indices 4-uples
1251 * \param[out] triangle_regions each generated triangle has
1252 * a region index, that corresponds to the index of the
1253 * Voronoi cell the triangle belongs to
1254 * \param[out] tet_regions each generated tetrahedron has
1255 * a region index, that corresponds to the index of the
1256 * Voronoi cell the tetrahedron belongs to
1257 * \param[in] cell_borders_only if true, only the surfacic
1258 * borders of the volumetric cells are saved in the mesh, else
1259 * volumetric cells are tetrahedralized.
1260 * \pre dim <= delaunay->dimension()
1261 */
1262 ✗ BuildVolumetricRVD(
1263 GenRestrictedVoronoiDiagram& RVD,
1264 coord_index_t dim,
1265 vector<double>& vertices,
1266 vector<index_t>& triangle_vertex_indices,
1267 vector<index_t>& tet_vertex_indices,
1268 vector<index_t>& triangle_regions,
1269 vector<index_t>& tet_regions,
1270 bool cell_borders_only
1271 ) :
1272 ✗ delaunay_(RVD.delaunay()),
1273 ✗ mesh_(RVD.mesh()),
1274 ✗ dim_(dim),
1275 ✗ vertices_(vertices),
1276 ✗ triangle_vertex_indices_(triangle_vertex_indices),
1277 ✗ tet_vertex_indices_(tet_vertex_indices),
1278 ✗ triangle_regions_(triangle_regions),
1279 ✗ tet_regions_(tet_regions),
1280 ✗ cell_borders_only_(cell_borders_only)
1281 {
1282 ✗ vertices_.clear();
1283 ✗ triangle_vertex_indices_.clear();
1284 ✗ tet_vertex_indices_.clear();
1285 ✗ triangle_regions_.clear();
1286 ✗ tet_regions_.clear();
1287
1288 // The first vertices are copied from Delaunay,
1289 // the other ones will be created during the traversal
1290 ✗ nb_vertices_ = delaunay_->nb_vertices();
1291 ✗ vertices_.resize(nb_vertices_ * dim);
1292 ✗ for(index_t v = 0; v < delaunay_->nb_vertices(); ++v) {
1293 ✗ for(coord_index_t c = 0; c < dim; ++c) {
1294 ✗ vertices_[v * dim + c] = delaunay_->vertex_ptr(v)[c];
1295 }
1296 }
1297 ✗ vertex_map_.set_first_vertex_index(nb_vertices_);
1298 ✗ }
1299
1300 /**
1301 * \brief The callback called for each integration simplex.
1302 * \param[in] v index of current center vertex
1303 * \param[in] v_adj (unused here) is the index of the Voronoi cell
1304 * adjacent to t accros facet (\p v1, \p v2, \p v3) or
1305 * NO_INDEX if it does not exists
1306 * \param[in] t (unused here) is the index of the current
1307 * tetrahedron
1308 * \param[in] t_adj (unused here) is the index of the
1309 * tetrahedron adjacent to t accros facet (\p v1, \p v2, \p v3)
1310 * or NO_INDEX if it does not exists
1311 * \param[in] v1 first vertex of current integration simplex
1312 * \param[in] v2 second vertex of current integration simplex
1313 * \param[in] v3 third vertex of current integration simplex
1314 */
1315 ✗ void operator() (
1316 index_t v, index_t v_adj,
1317 index_t t, index_t t_adj,
1318 const Vertex& v1, const Vertex& v2, const Vertex& v3
1319 ) {
1320 ✗ geo_argused(v_adj);
1321 ✗ geo_argused(t);
1322
1323 ✗ if(cell_borders_only_) {
1324 ✗ if(v > v_adj) {
1325 ✗ index_t iv2 = find_or_create_vertex(v, v1);
1326 ✗ index_t iv3 = find_or_create_vertex(v, v2);
1327 ✗ index_t iv4 = find_or_create_vertex(v, v3);
1328 ✗ triangle_vertex_indices_.push_back(iv4);
1329 ✗ triangle_vertex_indices_.push_back(iv3);
1330 ✗ triangle_vertex_indices_.push_back(iv2);
1331 ✗ triangle_regions_.push_back(v);
1332 }
1333 } else {
1334 ✗ index_t iv1 = v;
1335 ✗ index_t iv2 = find_or_create_vertex(v, v1);
1336 ✗ index_t iv3 = find_or_create_vertex(v, v2);
1337 ✗ index_t iv4 = find_or_create_vertex(v, v3);
1338
1339 // Triangle v1,v2,v3 is on border if there is
1340 // no adjacent seed and no adjacent tet.
1341 ✗ if(v_adj == NO_INDEX && t_adj == NO_INDEX) {
1342 ✗ triangle_vertex_indices_.push_back(iv4);
1343 ✗ triangle_vertex_indices_.push_back(iv3);
1344 ✗ triangle_vertex_indices_.push_back(iv2);
1345 ✗ triangle_regions_.push_back(v);
1346 }
1347
1348 ✗ tet_vertex_indices_.push_back(iv1);
1349 ✗ tet_vertex_indices_.push_back(iv2);
1350 ✗ tet_vertex_indices_.push_back(iv3);
1351 ✗ tet_vertex_indices_.push_back(iv4);
1352 ✗ tet_regions_.push_back(v);
1353 }
1354 ✗ }
1355
1356 /**
1357 * \brief The callback called for each tetrahedron
1358 * \param[in] v index of current center vertex
1359 * \param[in] v_adj (unused here) is the index of the Voronoi cell
1360 * adjacent to t accros facet (\p v1, \p v2, \p v3) or
1361 * NO_INDEX if it does not exists
1362 * \param[in] t (unused here) is the index of the current
1363 * tetrahedron
1364 * \param[in] t_adj (unused here) is the index of the
1365 * tetrahedron adjacent to t accros facet (\p v1, \p v2, \p v3)
1366 * or NO_INDEX if it does not exists
1367 * \param[in] v1 first vertex of current tetrahedron
1368 * \param[in] v2 second vertex of current tetrahedron
1369 * \param[in] v3 third vertex of current tetrahedron
1370 * \param[in] v4 fourth vertex of current tetrahedron
1371 */
1372 ✗ void operator() (
1373 index_t v, index_t v_adj,
1374 index_t t, index_t t_adj,
1375 const Vertex& v1, const Vertex& v2,
1376 const Vertex& v3, const Vertex& v4
1377 ) {
1378 ✗ geo_argused(v_adj);
1379 ✗ geo_argused(t);
1380 ✗ geo_argused(t_adj);
1381 ✗ index_t iv1 = vertices_.size() / dim_;
1382 ✗ for(index_t c = 0; c < dim_; ++c) {
1383 ✗ vertices_.push_back(v1.point()[c]);
1384 }
1385 ✗ index_t iv2 = vertices_.size() / dim_;
1386 ✗ for(index_t c = 0; c < dim_; ++c) {
1387 ✗ vertices_.push_back(v2.point()[c]);
1388 }
1389 ✗ index_t iv3 = vertices_.size() / dim_;
1390 ✗ for(index_t c = 0; c < dim_; ++c) {
1391 ✗ vertices_.push_back(v3.point()[c]);
1392 }
1393 ✗ index_t iv4 = vertices_.size() / dim_;
1394 ✗ for(index_t c = 0; c < dim_; ++c) {
1395 ✗ vertices_.push_back(v4.point()[c]);
1396 }
1397 ✗ tet_vertex_indices_.push_back(iv1);
1398 ✗ tet_vertex_indices_.push_back(iv2);
1399 ✗ tet_vertex_indices_.push_back(iv3);
1400 ✗ tet_vertex_indices_.push_back(iv4);
1401 ✗ tet_regions_.push_back(v);
1402 ✗ }
1403
1404 protected:
1405 /**
1406 * \brief Retrieves the index of a vertex given its symbolic
1407 * representation.
1408 * \param[in] center_vertex_id index of current Voronoi seed
1409 * \param[in] v symbolic and geometric representation of the vertex
1410 * \return the index of the vertex
1411 */
1412 ✗ index_t find_or_create_vertex(
1413 index_t center_vertex_id, const Vertex& v
1414 ) {
1415 ✗ index_t result = vertex_map_.find_or_create_vertex(
1416 ✗ center_vertex_id, v.sym()
1417 );
1418 ✗ if(result >= nb_vertices_) {
1419 ✗ geo_assert(result == nb_vertices_);
1420 ✗ nb_vertices_ = result + 1;
1421 ✗ for(coord_index_t c = 0; c < dim_; ++c) {
1422 ✗ vertices_.push_back(v.point()[c]);
1423 }
1424 }
1425 ✗ return result;
1426 }
1427
1428 private:
1429 const Delaunay* delaunay_;
1430 const Mesh* mesh_;
1431 coord_index_t dim_;
1432 vector<double>& vertices_;
1433 vector<index_t>& triangle_vertex_indices_;
1434 vector<index_t>& tet_vertex_indices_;
1435 vector<index_t>& triangle_regions_;
1436 vector<index_t>& tet_regions_;
1437 RVDVertexMap vertex_map_;
1438 index_t nb_vertices_;
1439 bool cell_borders_only_;
1440 };
1441
1442 72 void compute_RVD(
1443 Mesh& M, coord_index_t dim, bool cell_borders_only,
1444 bool integration_simplices
1445 ) override {
1446 72 bool sym = RVD_.symbolic();
1447 72 RVD_.set_symbolic(true);
1448
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 36 times.
 72 if(volumetric_) {
1449 ✗ if(dim == 0) {
1450 ✗ dim = dimension();
1451 }
1452 ✗ vector<double> vertices;
1453 ✗ vector<index_t> triangle_vertices;
1454 ✗ vector<index_t> tet_vertices;
1455 ✗ vector<index_t> triangle_regions;
1456 ✗ vector<index_t> tet_regions;
1457 ✗ if(cell_borders_only) {
1458 ✗ RVD_.for_each_volumetric_integration_simplex(
1459 ✗ BuildVolumetricRVD(
1460 ✗ RVD_, dim,
1461 vertices,
1462 triangle_vertices,
1463 tet_vertices,
1464 triangle_regions,
1465 tet_regions,
1466 cell_borders_only
1467 ),
1468 false, // Do not visit inner tetrahedra.
1469 true // Ensure that polygonal facets are triangulated
1470 // coherently.
1471 );
1472 } else {
1473 ✗ if(integration_simplices) {
1474 ✗ RVD_.for_each_volumetric_integration_simplex(
1475 ✗ BuildVolumetricRVD(
1476 ✗ RVD_, dim,
1477 vertices,
1478 triangle_vertices,
1479 tet_vertices,
1480 triangle_regions,
1481 tet_regions,
1482 cell_borders_only
1483 ),
1484 false, // Do not visit inner tetrahedra.
1485 true // Ensure that polygonal facets are
1486 // triangulated coherently.
1487 );
1488 } else {
1489 ✗ RVD_.for_each_tetrahedron(
1490 ✗ BuildVolumetricRVD(
1491 ✗ RVD_, dim,
1492 vertices,
1493 triangle_vertices,
1494 tet_vertices,
1495 triangle_regions,
1496 tet_regions,
1497 cell_borders_only
1498 )
1499 );
1500 }
1501 }
1502
1503 ✗ M.clear(true); // keep attributes
1504
1505 ✗ M.vertices.assign_points(vertices,dim,true);
1506 ✗ M.facets.assign_triangle_mesh(triangle_vertices, true);
1507 ✗ M.cells.assign_tet_mesh(tet_vertices, true);
1508
1509 // TODO: use Attribute::assign(vector, steal_args)
1510 // when it is there...
1511
1512 ✗ if(M.facets.nb() != 0) {
1513 ✗ Attribute<index_t> facet_region_attr(
1514 ✗ M.facets.attributes(), "region"
1515 );
1516 ✗ for(index_t f=0; f<M.facets.nb(); ++f) {
1517 ✗ facet_region_attr[f] = triangle_regions[f];
1518 }
1519 ✗ }
1520
1521 ✗ if(M.cells.nb() != 0) {
1522 ✗ Attribute<index_t> cell_region_attr(
1523 ✗ M.cells.attributes(), "region"
1524 );
1525 ✗ for(index_t c=0; c<M.cells.nb(); ++c) {
1526 ✗ cell_region_attr[c] = tet_regions[c];
1527 }
1528 ✗ }
1529
1530 ✗ if(cell_borders_only) {
1531 ✗ mesh_repair(M, MESH_REPAIR_TOPOLOGY);
1532 } else {
1533 ✗ M.facets.connect();
1534 }
1535 ✗ } else {
1536
1/2
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
 72 RVDMeshBuilder builder(
1537 &M, mesh_, delaunay_
1538 );
1539
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 36 times.
 72 if(dim != 0) {
1540 ✗ builder.set_dimension(dim);
1541 }
1542
1/2
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
 72 RVD_.for_each_polygon(
1543
1/2
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
 144 BuildRVD<RVDMeshBuilder>(RVD_, builder)
1544 );
1545 72 }
1546 72 RVD_.set_symbolic(sym);
1547
2/4
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 36 times.
✗ Branch 5 not taken.
 72 M.show_stats("RVD");
1548 72 }
1549
1550 /********************************************************************/
1551
1552 ✗ void compute_RVC(
1553 index_t i,
1554 Mesh& M,
1555 Mesh& result,
1556 bool copy_symbolic_info
1557 ) override {
1558 ✗ Mesh* tmp_mesh = mesh_;
1559 ✗ mesh_ = &M;
1560 ✗ RVD_.set_mesh(&M);
1561 ✗ typename GenRestrictedVoronoiDiagram::Polyhedron Cell(dimension());
1562 ✗ Cell.initialize_from_surface_mesh(&M, RVD_.symbolic());
1563 ✗ RVD_.intersect_cell_cell(i, Cell);
1564 ✗ Cell.convert_to_mesh(&result, copy_symbolic_info);
1565 ✗ mesh_ = tmp_mesh;
1566 ✗ RVD_.set_mesh(tmp_mesh);
1567 ✗ }
1568
1569 /********************************************************************/
1570
1571 ✗ void for_each_polygon(
1572 GEO::RVDPolygonCallback& callback,
1573 bool symbolic,
1574 bool connected_comp_priority,
1575 bool parallel
1576 ) override {
1577 ✗ bool sym_backup = RVD_.symbolic();
1578 ✗ RVD_.set_symbolic(symbolic);
1579 ✗ RVD_.set_connected_components_priority(connected_comp_priority);
1580 ✗ callback.begin();
1581 ✗ if(parallel) {
1582 ✗ compute_with_polygon_callback(callback);
1583 } else {
1584 ✗ PolygonCallbackAction action(RVD_,callback);
1585 ✗ RVD_.for_each_polygon(action);
1586 }
1587 ✗ callback.end();
1588 ✗ RVD_.set_symbolic(sym_backup);
1589 ✗ RVD_.set_connected_components_priority(false);
1590 ✗ }
1591
1592 /********************************************************************/
1593
1594 8 void for_each_polyhedron(
1595 GEO::RVDPolyhedronCallback& callback,
1596 bool symbolic,
1597 bool connected_comp_priority,
1598 bool parallel
1599 ) override {
1600 8 bool sym_backup = RVD_.symbolic();
1601 8 RVD_.set_symbolic(symbolic);
1602 8 RVD_.set_connected_components_priority(connected_comp_priority);
1603 8 callback.set_dimension(RVD_.mesh()->vertices.dimension());
1604 8 callback.begin();
1605
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 8 if(parallel) {
1606 ✗ compute_with_polyhedron_callback(callback);
1607 } else {
1608 8 RVD_.for_each_polyhedron(callback);
1609 }
1610 8 callback.end();
1611 8 RVD_.set_symbolic(sym_backup);
1612 8 RVD_.set_connected_components_priority(false);
1613 8 }
1614
1615 /********************************************************************/
1616
1617 /**
1618 * \brief Does the actual computation for a specific part
1619 * in multithread mode.
1620 * \param[in] t the index of the part.
1621 * \pre \p t < nb_parts()
1622 */
1623 3544 void run_thread(index_t t) {
1624
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 1772 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 3544 geo_assert(t < nb_parts());
1625 3544 thisclass& T = part(t);
1626
2/7
✓ Branch 0 taken 240 times.
✓ Branch 1 taken 1532 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
 3544 switch(thread_mode_) {
1627 480 case MT_LLOYD:
1628 {
1629 480 T.compute_centroids(arg_vectors_, arg_scalars_);
1630 480 } break;
1631 3064 case MT_NEWTON:
1632 {
1633 3064 T.compute_CVT_func_grad(T.funcval_, arg_vectors_);
1634 3064 } break;
1635 ✗ case MT_INT_SMPLX:
1636 {
1637 ✗ T.compute_integration_simplex_func_grad(
1638 ✗ T.funcval_, arg_vectors_, simplex_func_
1639 );
1640 ✗ } break;
1641 ✗ case MT_POLYG:
1642 {
1643 ✗ T.compute_with_polygon_callback(
1644 ✗ *polygon_callback_
1645 );
1646 ✗ } break;
1647 ✗ case MT_POLYH:
1648 {
1649 ✗ T.compute_with_polyhedron_callback(
1650 ✗ *polyhedron_callback_
1651 );
1652 ✗ } break;
1653 ✗ case MT_NONE:
1654 ✗ geo_assert_not_reached;
1655 }
1656 3544 }
1657
1658 16 bool compute_initial_sampling_on_surface(
1659 double* p, index_t nb_points, bool verbose
1660 ) override {
1661
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 8 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 16 geo_assert(mesh_->facets.are_simplices());
1662
1663 // We do that here, since this triggers partitioning,
1664 // that improves data locality. Then data locality is
1665 // inherited by the generated points.
1666 16 create_threads();
1667
1668
4/6
✓ Branch 0 taken 7 times.
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 7 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 7 times.
✗ Branch 5 not taken.
 16 if(verbose && facets_begin_ == NO_INDEX && facets_end_ == NO_INDEX) {
1669
2/4
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 7 times.
✗ Branch 5 not taken.
 42 Logger::out("RVD")
1670
1/2
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
 14 << "Computing initial sampling on surface, using dimension="
1671
2/4
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 7 times.
✗ Branch 5 not taken.
 14 << index_t(dimension_) << std::endl;
1672 }
1673
1674 32 return mesh_generate_random_samples_on_surface<DIM>(
1675 16 *mesh_, p, nb_points, vertex_weight_, facets_begin_, facets_end_
1676 16 );
1677 }
1678
1679 8 bool compute_initial_sampling_in_volume(
1680 double* p, index_t nb_points, bool verbose
1681 ) override {
1682
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 4 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 8 geo_assert(mesh_->cells.nb() != 0);
1683
1684 // We do that here, since this triggers partitioning,
1685 // that improves data locality. Then data locality is
1686 // inherited by the generated points.
1687 8 create_threads();
1688
1689
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 8 if(verbose && tets_begin_ == NO_INDEX && tets_end_ == NO_INDEX) {
1690 ✗ Logger::out("RVD")
1691 ✗ << "Computing initial sampling in volume, using dimension="
1692 ✗ << index_t(dimension_) << std::endl;
1693 }
1694
1695 16 return mesh_generate_random_samples_in_volume<DIM>(
1696 8 *mesh_, p, nb_points, vertex_weight_, tets_begin_, tets_end_
1697 8 );
1698 }
1699
1700 /**
1701 * \brief Creates the data structures for fast projection.
1702 * \details It decomposes the surface into triangles, and
1703 * stores the vertices in a KdTree.
1704 */
1705 ✗ void prepare_projection() {
1706 ✗ if(!mesh_vertices_.is_null()) {
1707 ✗ return;
1708 }
1709
1710 // Step 1: get triangles
1711 ✗ for(index_t f = 0; f < mesh_->facets.nb(); f++) {
1712 ✗ index_t i = mesh_->facets.corners_begin(f);
1713 ✗ for(
1714 ✗ index_t j = i + 1;
1715 ✗ j + 1 < mesh_->facets.corners_end(f); j++
1716 ) {
1717 ✗ triangles_.push_back(mesh_->facet_corners.vertex(i));
1718 ✗ triangles_.push_back(mesh_->facet_corners.vertex(j));
1719 ✗ triangles_.push_back(mesh_->facet_corners.vertex(j + 1));
1720 }
1721 }
1722 ✗ nb_triangles_ = index_t(triangles_.size() / 3);
1723
1724 // Step 2: get vertices stars
1725 // Step 2.1: get one-ring neighborhood
1726 ✗ vector<vector<index_t> > stars2(mesh_->vertices.nb());
1727 ✗ for(index_t t = 0; t < nb_triangles_; t++) {
1728 ✗ stars2[triangles_[3 * t]].push_back(t);
1729 ✗ stars2[triangles_[3 * t + 1]].push_back(t);
1730 ✗ stars2[triangles_[3 * t + 2]].push_back(t);
1731 }
1732
1733 // Step 2.2: get two-ring neighborhood
1734 ✗ stars_.resize(mesh_->vertices.nb());
1735 ✗ for(index_t i = 0; i < stars2.size(); i++) {
1736 ✗ vector<index_t> Ni;
1737 ✗ for(index_t j = 0; j < stars2[i].size(); j++) {
1738 ✗ index_t t = stars2[i][j];
1739 ✗ for(index_t iv = 0; iv < 3; iv++) {
1740 ✗ index_t v = triangles_[3 * t + iv];
1741 ✗ if(v != i) {
1742 ✗ Ni.push_back(v);
1743 }
1744 }
1745 }
1746 ✗ sort_unique(Ni);
1747 ✗ for(index_t j = 0; j < Ni.size(); j++) {
1748 ✗ index_t k = Ni[j];
1749 ✗ stars_[i].insert(
1750 ✗ stars_[i].end(), stars2[k].begin(), stars2[k].end()
1751 );
1752 }
1753 ✗ sort_unique(stars_[i]);
1754 }
1755
1756 // Step 3: create search structure
1757 ✗ mesh_vertices_ = Delaunay::create(dimension_, "NN");
1758 ✗ index_t nb_vertices = mesh_->vertices.nb();
1759
1760 // TODO: BUG !! mesh_vertices_ keeps a ref. to mesh_vertices
1761 // that is destroyed when leaving this function.
1762 ✗ vector<double> mesh_vertices(nb_vertices * dimension_);
1763 ✗ for(index_t i = 0; i < nb_vertices; i++) {
1764 ✗ for(index_t coord = 0; coord < dimension_; coord++) {
1765 ✗ mesh_vertices[i * dimension_ + coord] =
1766 ✗ mesh_->vertices.point_ptr(i)[coord];
1767 }
1768 }
1769 ✗ mesh_vertices_->set_vertices(nb_vertices, mesh_vertices.data());
1770 ✗ }
1771
1772 ✗ void project_points_on_surface(
1773 index_t nb_points, double* points, vec3* nearest, bool do_project
1774 ) override {
1775
1776 ✗ prepare_projection();
1777
1778 ✗ if(use_exact_projection_) {
1779 ✗ for(index_t p = 0; p < nb_points; p++) {
1780 ✗ Point P(points + p * dimension_);
1781 ✗ double d2 = Numeric::max_float64();
1782 ✗ for(index_t t = 0; t < nb_triangles_; t++) {
1783 ✗ const Point& p1 = mesh_vertex(triangles_[3 * t]);
1784 ✗ const Point& p2 = mesh_vertex(triangles_[3 * t + 1]);
1785 ✗ const Point& p3 = mesh_vertex(triangles_[3 * t + 2]);
1786
1787 double l1, l2, l3;
1788 ✗ Point nearestP;
1789
1790 ✗ double cur_d2 = Geom::point_triangle_squared_distance(
1791 P, p1, p2, p3, nearestP, l1, l2, l3
1792 );
1793
1794 ✗ if(cur_d2 < d2) {
1795 ✗ d2 = cur_d2;
1796 const vec3& p1_R3 =
1797 ✗ R3_embedding(triangles_[3 * t]);
1798 const vec3& p2_R3 =
1799 ✗ R3_embedding(triangles_[3 * t + 1]);
1800 const vec3& p3_R3 =
1801 ✗ R3_embedding(triangles_[3 * t + 2]);
1802 ✗ nearest[p] = l1 * p1_R3 + l2 * p2_R3 + l3 * p3_R3;
1803 ✗ if(do_project) {
1804 ✗ for(coord_index_t
1805 ✗ coord = 0; coord < dimension_; coord++) {
1806 ✗ (points + p * dimension_)[coord] =
1807 ✗ nearestP[coord];
1808 }
1809 }
1810 }
1811 }
1812 }
1813 ✗ return;
1814 }
1815
1816 // find nearest point on surface in star of nearest vertex
1817 ✗ for(index_t p = 0; p < nb_points; p++) {
1818 ✗ Point P(points + p * dimension_);
1819 ✗ index_t v = mesh_vertices_->nearest_vertex(
1820 ✗ points + p * dimension_
1821 );
1822 ✗ double d2 = Numeric::max_float64();
1823 ✗ nearest[p] = R3_embedding(v);
1824 ✗ for(index_t i = 0; i < stars_[v].size(); i++) {
1825 ✗ index_t t = stars_[v][i];
1826 ✗ const Point& p1 = mesh_vertex(triangles_[3 * t]);
1827 ✗ const Point& p2 = mesh_vertex(triangles_[3 * t + 1]);
1828 ✗ const Point& p3 = mesh_vertex(triangles_[3 * t + 2]);
1829
1830 double l1, l2, l3;
1831 ✗ Point nearestP;
1832 ✗ double cur_d2 = Geom::point_triangle_squared_distance(
1833 P, p1, p2, p3, nearestP, l1, l2, l3
1834 );
1835 ✗ if(cur_d2 < d2) {
1836 ✗ d2 = cur_d2;
1837 ✗ const vec3& p1_R3 = R3_embedding(triangles_[3 * t]);
1838 ✗ const vec3& p2_R3 = R3_embedding(triangles_[3 * t + 1]);
1839 ✗ const vec3& p3_R3 = R3_embedding(triangles_[3 * t + 2]);
1840 ✗ nearest[p] = l1 * p1_R3 + l2 * p2_R3 + l3 * p3_R3;
1841 ✗ if(do_project) {
1842 ✗ for(coord_index_t coord = 0;
1843 ✗ coord < dimension_; coord++
1844 ) {
1845 ✗ (points + p * dimension_)[coord] =
1846 ✗ nearestP[coord];
1847 }
1848 }
1849 }
1850 }
1851 }
1852 }
1853
1854 /********************************************************************/
1855
1856 /**
1857 * \brief Implementation class for computing the restricted Delaunay
1858 * triangulation.
1859 * \details To be used as a template argument
1860 * to RVD::for_each_primal_triangle().
1861 */
1862 class GetPrimalTriangles {
1863 public:
1864 /**
1865 * \brief Creates a new GetPrimalTriangles.
1866 * \param[out] triangles where to store the triangles
1867 */
1868 ✗ GetPrimalTriangles(
1869 vector<index_t>& triangles
1870 ) :
1871 ✗ triangles_(triangles) {
1872 ✗ }
1873
1874 /**
1875 * \brief The callback called for each primal triangle.
1876 * \param[in] v1 index of the first vertex
1877 * \param[in] v2 index of the second vertex
1878 * \param[in] v3 index of the third vertex
1879 */
1880 ✗ void operator() (index_t v1, index_t v2, index_t v3) {
1881 ✗ triangles_.push_back(v1);
1882 ✗ triangles_.push_back(v2);
1883 ✗ triangles_.push_back(v3);
1884 ✗ }
1885
1886 private:
1887 vector<index_t>& triangles_;
1888 };
1889
1890 /**
1891 * \brief Implementation class for computing the restricted Delaunay
1892 * triangulation of the connected components.
1893 *
1894 * \details The difference with GetPrimalTriangles is that when a
1895 * restricted Voronoi cell has multiple connected components,
1896 * more triangles are generated to account for the topology.
1897 * To be used as a template argument to RVD::for_each_polygon().
1898 * The RestrictedVoronoiDiagram needs to be in connected-components
1899 * priority mode.
1900 */
1901 class GetConnectedComponentsPrimalTriangles {
1902 public:
1903 /** Internal representation of the polygons. */
1904 typedef typename GenRestrictedVoronoiDiagram::Polygon Polygon;
1905
1906 /** Internal representation of the vertices. */
1907 typedef typename GenRestrictedVoronoiDiagram::Vertex Vertex;
1908
1909 static constexpr index_t UNINITIALIZED = index_t(-1);
1910 static constexpr index_t MULTI_COMP = index_t(-2);
1911 static constexpr index_t ON_BORDER = index_t(-3);
1912
1913 /**
1914 * \brief Constructs a new GetConnectedComponentsPrimalTriangles.
1915 * \param[in] RVD the restricted Voronoi diagram
1916 * \param[out] triangles where to store the triangles
1917 * \param[out] vertices where to store the vertices
1918 * \param[in] dimension dimension of the restricted Voronoi diagram
1919 * \param[in] mode a combination of constants defined in RDTMode
1920 * \param[in] seed_is_locked specifies for each seed whether it
1921 * can be moved (to RVC centroid or projected on surface). If
1922 * left uninitialized, all the seeds can be moved.
1923 * \param[in] AABB an axis-aligned bounding box tree defined on
1924 * the input surface. It is used if one of (RDT_SELECT_NEAREST,
1925 * RDT_PROJECT_ON_SURFACE) is set in \p mode. If needed and not
1926 * specified, then a new one is created locally.
1927 */
1928 14 GetConnectedComponentsPrimalTriangles(
1929 const GenRestrictedVoronoiDiagram& RVD,
1930 vector<index_t>& triangles,
1931 vector<double>& vertices,
1932 coord_index_t dimension,
1933 RDTMode mode,
1934 const std::vector<bool>& seed_is_locked,
1935 MeshFacetsAABB* AABB = nullptr
1936 ) :
1937 14 RVD_(RVD),
1938 14 dimension_(dimension),
1939 14 triangles_(triangles),
1940 14 vertices_(vertices),
1941 14 m_(0.0),
1942 14 cur_seed_(NO_INDEX),
1943 14 cur_vertex_(0),
1944 14 use_RVC_centroids_((mode & RDT_RVC_CENTROIDS) != 0),
1945 14 select_nearest_((mode & RDT_SELECT_NEAREST) != 0),
1946 14 project_on_surface_((mode & RDT_PROJECT_ON_SURFACE) != 0),
1947 14 seed_is_locked_(seed_is_locked),
1948 14 prefer_seeds_((mode & RDT_PREFER_SEEDS) != 0),
1949 14 AABB_(AABB)
1950 {
1951
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 2 times.
 14 if(prefer_seeds_) {
1952
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 10 seed_to_vertex_.assign(
1953 10 RVD.delaunay()->nb_vertices(), UNINITIALIZED
1954 );
1955 }
1956 14 }
1957
1958 /**
1959 * \brief The callback called for each restricted Voronoi cell.
1960 * \param[in] s1 index of current center vertex
1961 * \param[in] P current restricted Voronoi cell
1962 */
1963 359096 void operator() (index_t s1, const Polygon& P) {
1964
2/2
✓ Branch 1 taken 35000 times.
✓ Branch 2 taken 144548 times.
 359096 if(RVD_.connected_component_changed()) {
1965
2/2
✓ Branch 0 taken 34993 times.
✓ Branch 1 taken 7 times.
 70000 if(cur_seed_ != NO_INDEX) {
1966 69986 end_connected_component();
1967 }
1968 70000 begin_connected_component(s1);
1969 }
1970
1971
4/4
✓ Branch 0 taken 160862 times.
✓ Branch 1 taken 18686 times.
✓ Branch 2 taken 152810 times.
✓ Branch 3 taken 8052 times.
 359096 if(prefer_seeds_ && !component_on_border_) {
1972 // NOTE: there is one vertex shift between
1973 // adjacent facet and adjacent seed, there
1974 // must be something wrong in the way the
1975 // combinatorial information is initialized,
1976 // to be checked !!! (should be the i index
1977 // for both)
1978
2/2
✓ Branch 1 taken 597533 times.
✓ Branch 2 taken 152547 times.
 1500160 for(index_t i=0; i<P.nb_vertices(); ++i) {
1979 1195066 index_t j = (i+1) % P.nb_vertices();
1980 1195066 if(
1981
6/6
✓ Branch 2 taken 278953 times.
✓ Branch 3 taken 318580 times.
✓ Branch 4 taken 263 times.
✓ Branch 5 taken 278690 times.
✓ Branch 6 taken 263 times.
✓ Branch 7 taken 597270 times.
 1752972 P.vertex(i).adjacent_facet() == -1 &&
1982 557906 P.vertex(j).adjacent_seed() == -1
1983 ) {
1984 526 component_on_border_ = true;
1985 526 break;
1986 }
1987 }
1988 }
1989
1990 // Accumulate mass and barycenter
1991 359096 index_t vbase = cur_vertex_ * dimension_;
1992
2/2
✓ Branch 1 taken 358564 times.
✓ Branch 2 taken 179548 times.
 1076224 for(index_t i = 1; i + 1 < P.nb_vertices(); ++i) {
1993 717128 double cur_m = Geom::triangle_area(
1994 717128 P.vertex(0).point(),
1995 717128 P.vertex(i).point(),
1996 717128 P.vertex(i + 1).point(), dimension_
1997 );
1998
2/2
✓ Branch 0 taken 1980642 times.
✓ Branch 1 taken 358564 times.
 4678412 for(coord_index_t c = 0; c < dimension_; ++c) {
1999 3961284 vertices_[vbase + c] += cur_m / 3.0 * (
2000 3961284 P.vertex(0).point()[c] +
2001 3961284 P.vertex(i).point()[c] +
2002 3961284 P.vertex(i + 1).point()[c]
2003 );
2004 }
2005 717128 m_ += cur_m;
2006 }
2007
2008 // Detect Voronoi vertices and generate
2009 // triangles.
2010 // Note: they can be generated several times,
2011 // since we cannot know in advance whether
2012 // the other instances of the Voronoi vertex
2013 // were finalized or not (i.e. have their
2014 // three vertices ready).
2015 // Duplicate triangles are then filtered-out
2016 // by client code.
2017
2/2
✓ Branch 1 taken 717660 times.
✓ Branch 2 taken 179548 times.
 1794416 for(index_t i = 0; i < P.nb_vertices(); ++i) {
2018 1435320 const Vertex& V = P.vertex(i);
2019
2/2
✓ Branch 2 taken 209199 times.
✓ Branch 3 taken 508461 times.
 1435320 if(V.sym().nb_bisectors() == 2) {
2020
1/2
✓ Branch 2 taken 209199 times.
✗ Branch 3 not taken.
 418398 index_t s2 = V.sym().bisector(0);
2021
1/2
✓ Branch 2 taken 209199 times.
✗ Branch 3 not taken.
 418398 index_t s3 = V.sym().bisector(1);
2022
1/2
✓ Branch 2 taken 209199 times.
✗ Branch 3 not taken.
 418398 index_t f = V.sym().boundary_facet(0);
2023
2024 418398 index_t v1 = RVD_.current_connected_component();
2025 418398 index_t v2 = index_t(
2026
1/2
✓ Branch 1 taken 209199 times.
✗ Branch 2 not taken.
 418398 RVD_.get_facet_seed_connected_component(f,s2)
2027 );
2028 418398 index_t v3 = index_t(
2029
1/2
✓ Branch 1 taken 209199 times.
✗ Branch 2 not taken.
 418398 RVD_.get_facet_seed_connected_component(f,s3)
2030 );
2031
2032
4/4
✓ Branch 0 taken 103793 times.
✓ Branch 1 taken 105406 times.
✓ Branch 2 taken 69733 times.
✓ Branch 3 taken 34060 times.
 418398 if(v2 != NO_INDEX && v3 != NO_INDEX) {
2033
1/2
✓ Branch 1 taken 69733 times.
✗ Branch 2 not taken.
 139466 triangles_.push_back(v1);
2034
1/2
✓ Branch 1 taken 69733 times.
✗ Branch 2 not taken.
 139466 triangles_.push_back(v2);
2035
1/2
✓ Branch 1 taken 69733 times.
✗ Branch 2 not taken.
 139466 triangles_.push_back(v3);
2036 }
2037 }
2038 }
2039 359096 }
2040
2041 /**
2042 * \brief The destructor
2043 */
2044 14 ~GetConnectedComponentsPrimalTriangles() {
2045
2046 14 bool owns_AABB = false;
2047
2/4
✓ Branch 0 taken 7 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 7 times.
 14 if(select_nearest_ || project_on_surface_) {
2048 ✗ if(AABB_ == nullptr) {
2049 // Construct an axis-aligned bounding box tree,
2050 // do not reorder the mesh (needs to be pre-reordered)
2051 ✗ AABB_ = new MeshFacetsAABB(
2052 ✗ *const_cast<Mesh*>(RVD_.mesh()), AABB_NOREORDER
2053 );
2054 ✗ owns_AABB = true;
2055 }
2056 }
2057
2058
1/2
✓ Branch 0 taken 7 times.
✗ Branch 1 not taken.
 14 if(cur_seed_ != NO_INDEX) {
2059 14 end_connected_component();
2060 }
2061
2062
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 2 times.
 14 if(prefer_seeds_) {
2063
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
 10 if(select_nearest_) {
2064 ✗ for(index_t s=0; s<seed_to_vertex_.size(); ++s) {
2065 ✗ if(
2066 ✗ seed_to_vertex_[s] != MULTI_COMP &&
2067 ✗ seed_to_vertex_[s] != UNINITIALIZED &&
2068 ✗ seed_to_vertex_[s] != ON_BORDER
2069 ) {
2070 ✗ index_t vbase = seed_to_vertex_[s] * dimension_;
2071
2072 const double* seed_ptr =
2073 ✗ RVD_.delaunay()->vertex_ptr(s);
2074
2075
2076 // At this step, vertex_ptr contains
2077 // the centroid of the connected component
2078 // of the RVC, we now determine whether it
2079 // should be replaced by the
2080 // seed (or by a projection onto the surface).
2081
2082 ✗ const double* vertex_ptr = &(vertices_[vbase]);
2083
2084 // If the seed is nearer to the surface
2085 // than the
2086 // centroid of the connected component of the
2087 // restricted Voronoi cell, then use the seed.
2088
2089 double seed_dist;
2090 ✗ vec3 seed_projection;
2091 double vertex_dist;
2092 ✗ vec3 vertex_projection;
2093
2094 ✗ AABB_->nearest_facet(
2095 ✗ vec3(seed_ptr), seed_projection, seed_dist
2096 );
2097 ✗ AABB_->nearest_facet(
2098 ✗ vec3(vertex_ptr),
2099 vertex_projection, vertex_dist
2100 );
2101
2102 ✗ if(seed_dist < vertex_dist) {
2103 ✗ if(project_on_surface_) {
2104 ✗ for(
2105 ✗ coord_index_t c = 0;
2106 ✗ c < dimension_; ++c
2107 ) {
2108 ✗ vertices_[vbase + c] =
2109 seed_projection[c];
2110 }
2111 } else {
2112 ✗ for(coord_index_t c = 0;
2113 ✗ c < dimension_; ++c
2114 ) {
2115 ✗ vertices_[vbase + c] = seed_ptr[c];
2116 }
2117 }
2118 } else {
2119 ✗ if(project_on_surface_) {
2120 ✗ for(
2121 ✗ coord_index_t c = 0;
2122 ✗ c < dimension_; ++c
2123 ) {
2124 ✗ vertices_[vbase + c] =
2125 vertex_projection[c];
2126 }
2127 }
2128 }
2129 }
2130 }
2131 } else {
2132
2133 // Current mode: prefer seeds and not select nearest
2134 // Replace all points with the seeds (provided that
2135 // they do not correspond to multiple
2136 // connected components).
2137
2138
2/2
✓ Branch 1 taken 25000 times.
✓ Branch 2 taken 5 times.
 50010 for(index_t s=0; s<seed_to_vertex_.size(); ++s) {
2139 50000 if(
2140 50000 seed_to_vertex_[s] != MULTI_COMP &&
2141
4/6
✓ Branch 0 taken 25000 times.
✗ Branch 1 not taken.
✓ Branch 3 taken 25000 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 24737 times.
✓ Branch 6 taken 263 times.
 100000 seed_to_vertex_[s] != UNINITIALIZED &&
2142
2/2
✓ Branch 1 taken 24737 times.
✓ Branch 2 taken 263 times.
 50000 seed_to_vertex_[s] != ON_BORDER
2143 ) {
2144 49474 index_t vbase = seed_to_vertex_[s] * dimension_;
2145
2146 const double* seed_ptr =
2147 49474 RVD_.delaunay()->vertex_ptr(s);
2148
2149
2/2
✓ Branch 0 taken 148422 times.
✓ Branch 1 taken 24737 times.
 346318 for(coord_index_t c = 0; c < dimension_; ++c) {
2150 296844 vertices_[vbase + c] = seed_ptr[c];
2151 }
2152 }
2153 }
2154 }
2155 }
2156
2157 14 if(
2158
4/6
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 5 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 7 times.
 14 (!prefer_seeds_ || !select_nearest_) && project_on_surface_
2159 ) {
2160 ✗ for(index_t v=0; v<vertices_.size()/3; ++v) {
2161 ✗ vec3 p(
2162 ✗ vertices_[3*v], vertices_[3*v+1], vertices_[3*v+2]
2163 );
2164 ✗ vec3 q;
2165 double sq_dist;
2166 ✗ AABB_->nearest_facet(p,q,sq_dist);
2167 ✗ vertices_[3*v ] = q.x;
2168 ✗ vertices_[3*v+1] = q.y;
2169 ✗ vertices_[3*v+2] = q.z;
2170 }
2171 }
2172
2173
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 7 times.
 14 if(owns_AABB) {
2174 ✗ delete AABB_;
2175 ✗ AABB_ = nullptr;
2176 }
2177 14 }
2178
2179 protected:
2180 /**
2181 * \brief Tests whether a given seed is locked.
2182 */
2183 50000 bool seed_is_locked(index_t s) {
2184 return
2185
1/4
✗ Branch 1 not taken.
✓ Branch 2 taken 25000 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 50000 seed_is_locked_.size() > 0 &&
2186 50000 seed_is_locked_[s]
2187 ;
2188 }
2189
2190 /**
2191 * \brief Starts a new connected component.
2192 * \param[in] s the seed the connected component
2193 * is associated with.
2194 */
2195 70000 void begin_connected_component(index_t s) {
2196 70000 cur_seed_ = s;
2197
2/2
✓ Branch 0 taken 180000 times.
✓ Branch 1 taken 35000 times.
 430000 for(coord_index_t c = 0; c < dimension_; ++c) {
2198
1/2
✓ Branch 1 taken 180000 times.
✗ Branch 2 not taken.
 360000 vertices_.push_back(0.0);
2199 }
2200 70000 m_ = 0.0;
2201 70000 component_on_border_ = false;
2202 70000 }
2203
2204 /**
2205 * \brief Terminates the current connected component.
2206 */
2207 70000 void end_connected_component() {
2208
2209 70000 if(
2210 120000 !use_RVC_centroids_ ||
2211
5/6
✓ Branch 0 taken 25000 times.
✓ Branch 1 taken 10000 times.
✓ Branch 3 taken 25000 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 10263 times.
✓ Branch 6 taken 24737 times.
 120000 seed_is_locked(cur_seed_) ||
2212
2/2
✓ Branch 0 taken 263 times.
✓ Branch 1 taken 24737 times.
 50000 component_on_border_
2213 ) {
2214 // Copy seed
2215 20526 index_t vbase = cur_vertex_ * dimension_;
2216 const double* seed_ptr =
2217 20526 RVD_.delaunay()->vertex_ptr(cur_seed_);
2218
2/2
✓ Branch 0 taken 31578 times.
✓ Branch 1 taken 10263 times.
 83682 for(coord_index_t c = 0; c < dimension_; ++c) {
2219 63156 vertices_[vbase + c] = seed_ptr[c];
2220 }
2221 } else {
2222 // Use restricted Voronoi
2223 // cell component's centroid.
2224
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 24737 times.
 49474 double scal = (m_ < 1e-30 ? 0.0 : 1.0 / m_);
2225 49474 index_t vbase = cur_vertex_ * dimension_;
2226
2/2
✓ Branch 0 taken 148422 times.
✓ Branch 1 taken 24737 times.
 346318 for(coord_index_t c = 0; c < dimension_; ++c) {
2227 296844 vertices_[vbase + c] *= scal;
2228 }
2229 }
2230
2/2
✓ Branch 0 taken 25000 times.
✓ Branch 1 taken 10000 times.
 70000 if(prefer_seeds_) {
2231
2/2
✓ Branch 0 taken 263 times.
✓ Branch 1 taken 24737 times.
 50000 if(component_on_border_) {
2232 526 seed_to_vertex_[cur_seed_] = ON_BORDER;
2233 }
2234
2/3
✓ Branch 1 taken 24737 times.
✓ Branch 2 taken 263 times.
✗ Branch 3 not taken.
 50000 switch(seed_to_vertex_[cur_seed_]) {
2235 49474 case UNINITIALIZED:
2236 49474 seed_to_vertex_[cur_seed_] = cur_vertex_;
2237 49474 break;
2238 526 case ON_BORDER:
2239 526 break;
2240 ✗ default:
2241 ✗ seed_to_vertex_[cur_seed_] = MULTI_COMP;
2242 ✗ break;
2243 }
2244 }
2245 70000 ++cur_vertex_;
2246 70000 }
2247
2248 private:
2249 const GenRestrictedVoronoiDiagram& RVD_;
2250 coord_index_t dimension_;
2251 vector<index_t>& triangles_;
2252 vector<double>& vertices_;
2253 double m_;
2254 index_t cur_seed_;
2255 index_t cur_vertex_;
2256 bool use_RVC_centroids_;
2257 bool select_nearest_;
2258 bool project_on_surface_;
2259 const std::vector<bool>& seed_is_locked_;
2260 bool prefer_seeds_;
2261 vector<index_t> seed_to_vertex_;
2262 bool component_on_border_;
2263 MeshFacetsAABB* AABB_;
2264 };
2265
2266 /**
2267 * \brief Implementation class for computing the restricted Delaunay
2268 * triangulation in volume mode.
2269 * \details To be used as a template argument
2270 * to RVD::for_each_primal_tetrahedron().
2271 */
2272 class GetPrimalTetrahedra {
2273 public:
2274 /**
2275 * \brief Creates a new GetPrimalTetrahedra.
2276 * \param[out] tetrahedra where to store the tetrahedra
2277 */
2278 ✗ GetPrimalTetrahedra(
2279 vector<index_t>& tetrahedra
2280 ) :
2281 ✗ tetrahedra_(tetrahedra) {
2282 ✗ }
2283
2284 /**
2285 * \brief The callback called for each primal tetrahedron.
2286 * \param[in] v1 index of the first vertex
2287 * \param[in] v2 index of the second vertex
2288 * \param[in] v3 index of the third vertex
2289 * \param[in] v4 index of the fourth vertex
2290 */
2291 ✗ void operator() (index_t v1, index_t v2, index_t v3, index_t v4) {
2292 ✗ tetrahedra_.push_back(v1);
2293 ✗ tetrahedra_.push_back(v2);
2294 ✗ tetrahedra_.push_back(v3);
2295 ✗ tetrahedra_.push_back(v4);
2296 ✗ }
2297
2298 private:
2299 vector<index_t>& tetrahedra_;
2300 };
2301
2302 14 void compute_RDT(
2303 vector<index_t>& simplices,
2304 vector<double>& embedding,
2305 RDTMode mode,
2306 const vector<bool>& seed_is_locked,
2307 MeshFacetsAABB* AABB
2308 ) override {
2309
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 7 times.
 14 if(volumetric_) {
2310 // For the moment, only simple mode is supported
2311 ✗ simplices.clear();
2312 ✗ RVD_.for_each_primal_tetrahedron(
2313 ✗ GetPrimalTetrahedra(simplices)
2314 );
2315 // Reorient the tetrahedra
2316 ✗ index_t nb_tetrahedra = simplices.size() / 4;
2317 ✗ for(index_t t = 0; t < nb_tetrahedra; ++t) {
2318 const double* p1 =
2319 ✗ delaunay()->vertex_ptr(simplices[4 * t]);
2320 const double* p2 =
2321 ✗ delaunay()->vertex_ptr(simplices[4 * t + 1]);
2322 const double*
2323 ✗ p3 = delaunay()->vertex_ptr(simplices[4 * t + 2]);
2324 const double*
2325 ✗ p4 = delaunay()->vertex_ptr(simplices[4 * t + 3]);
2326 ✗ if(PCK::orient_3d(p1, p2, p3, p4) < 0) {
2327 ✗ std::swap(simplices[4 * t], simplices[4 * t + 1]);
2328 }
2329 }
2330 ✗ embedding.clear();
2331 ✗ embedding.reserve(dimension_ * delaunay_->nb_vertices());
2332 ✗ for(index_t i = 0; i < delaunay_->nb_vertices(); i++) {
2333 ✗ for(coord_index_t coord = 0; coord < dimension_; coord++) {
2334 ✗ embedding.push_back(delaunay_->vertex_ptr(i)[coord]);
2335 }
2336 }
2337 } else {
2338
1/2
✓ Branch 0 taken 7 times.
✗ Branch 1 not taken.
 14 if((mode & RDT_MULTINERVE) != 0) {
2339 14 simplices.clear();
2340 14 embedding.clear();
2341 14 bool sym = RVD_.symbolic();
2342 14 RVD_.set_symbolic(true);
2343 14 RVD_.set_connected_components_priority(true);
2344
1/2
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
 14 RVD_.for_each_polygon(
2345 14 GetConnectedComponentsPrimalTriangles(
2346
1/2
✓ Branch 2 taken 7 times.
✗ Branch 3 not taken.
 14 RVD_, simplices, embedding, RVD_.dimension(),
2347 mode, seed_is_locked, AABB
2348 )
2349 );
2350 14 RVD_.set_symbolic(sym);
2351 14 RVD_.set_connected_components_priority(false);
2352 } else {
2353 // Simple mode: compute RDT, without any post-processing
2354 ✗ simplices.clear();
2355 ✗ RVD_.for_each_primal_triangle(
2356 ✗ GetPrimalTriangles(simplices)
2357 );
2358 ✗ embedding.clear();
2359 ✗ embedding.reserve(dimension_ * delaunay_->nb_vertices());
2360 ✗ for(index_t i = 0; i < delaunay_->nb_vertices(); i++) {
2361 ✗ for(
2362 ✗ coord_index_t coord = 0;
2363 ✗ coord < dimension_; ++coord
2364 ){
2365 ✗ embedding.push_back(
2366 ✗ delaunay_->vertex_ptr(i)[coord]
2367 );
2368 }
2369 }
2370 }
2371 }
2372 14 }
2373
2374 4454 void create_threads() override {
2375 // TODO: check if number of facets is not smaller than
2376 // number of threads
2377 // TODO: create parts even if facets range is specified
2378 // (and subdivide facets range)
2379 4454 if(
2380
2/2
✓ Branch 0 taken 455 times.
✓ Branch 1 taken 1772 times.
 4454 is_slave_ ||
2381
2/4
✓ Branch 0 taken 455 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 455 times.
 910 facets_begin_ != NO_INDEX || facets_end_ != NO_INDEX
2382 ) {
2383 3544 return;
2384 }
2385 910 index_t nb_parts_in = Process::maximum_concurrent_threads();
2386
2/2
✓ Branch 1 taken 12 times.
✓ Branch 2 taken 443 times.
 910 if(nb_parts() != nb_parts_in) {
2387
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 12 times.
 24 if(nb_parts_in == 1) {
2388 ✗ delete_threads();
2389 } else {
2390 24 vector<index_t> facet_ptr;
2391 24 vector<index_t> tet_ptr;
2392 24 mesh_partition(
2393
1/2
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
 24 *mesh_, MESH_PARTITION_HILBERT,
2394 facet_ptr, tet_ptr, nb_parts_in
2395 );
2396
1/2
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
 24 delete_threads();
2397
6/20
✓ Branch 0 taken 12 times.
✗ Branch 1 not taken.
✓ Branch 3 taken 12 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 48 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 48 times.
✓ Branch 9 taken 12 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 12 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
 120 parts_ = new thisclass[nb_parts_in];
2398 24 nb_parts_ = nb_parts_in;
2399
2/2
✓ Branch 1 taken 48 times.
✓ Branch 2 taken 12 times.
 120 for(index_t i = 0; i < nb_parts(); ++i) {
2400
1/2
✓ Branch 1 taken 48 times.
✗ Branch 2 not taken.
 96 part(i).mesh_ = mesh_;
2401
2/4
✓ Branch 1 taken 48 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 48 times.
✗ Branch 5 not taken.
 96 part(i).set_delaunay(delaunay_);
2402
1/2
✓ Branch 1 taken 48 times.
✗ Branch 2 not taken.
 96 part(i).R3_embedding_base_ = R3_embedding_base_;
2403
1/2
✓ Branch 1 taken 48 times.
✗ Branch 2 not taken.
 96 part(i).R3_embedding_stride_ = R3_embedding_stride_;
2404
1/2
✓ Branch 1 taken 48 times.
✗ Branch 2 not taken.
 96 part(i).has_weights_ = has_weights_;
2405
1/2
✓ Branch 1 taken 48 times.
✗ Branch 2 not taken.
 96 part(i).master_ = this;
2406
1/2
✓ Branch 1 taken 48 times.
✗ Branch 2 not taken.
 96 part(i).RVD_.set_mesh(mesh_);
2407
4/8
✓ Branch 1 taken 48 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 48 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 48 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 48 times.
✗ Branch 11 not taken.
 96 part(i).set_facets_range(
2408 facet_ptr[i], facet_ptr[i + 1]
2409 );
2410
2/4
✓ Branch 1 taken 48 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 48 times.
✗ Branch 6 not taken.
 96 part(i).set_exact_predicates(RVD_.exact_predicates());
2411
2/4
✓ Branch 1 taken 48 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 48 times.
✗ Branch 6 not taken.
 96 part(i).set_volumetric(volumetric());
2412
2/4
✓ Branch 1 taken 48 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 48 times.
✗ Branch 6 not taken.
 96 part(i).set_check_SR(RVD_.check_SR());
2413 }
2414
2/2
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 8 times.
 24 if(mesh_->cells.nb() != 0) {
2415
2/2
✓ Branch 1 taken 16 times.
✓ Branch 2 taken 4 times.
 40 for(index_t i = 0; i < nb_parts(); ++i) {
2416
4/8
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 16 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 16 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 16 times.
✗ Branch 11 not taken.
 32 part(i).set_tetrahedra_range(
2417 tet_ptr[i], tet_ptr[i + 1]
2418 );
2419 }
2420 }
2421
2/8
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 12 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
 24 geo_assert(!Process::is_running_threads());
2422 24 }
2423 }
2424 }
2425
2426 176 void set_volumetric(bool x) override {
2427 176 volumetric_ = x;
2428
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 88 times.
 176 for(index_t i = 0; i < nb_parts(); ++i) {
2429 ✗ part(i).set_volumetric(x);
2430 }
2431 176 }
2432
2433 96 void set_facets_range(
2434 index_t facets_begin, index_t facets_end
2435 ) override {
2436 96 RVD_.set_facets_range(facets_begin, facets_end);
2437 96 facets_begin_ = facets_begin;
2438 96 facets_end_ = facets_end;
2439 96 }
2440
2441 32 void set_tetrahedra_range(
2442 index_t tets_begin, index_t tets_end
2443 ) override {
2444 32 RVD_.set_tetrahedra_range(tets_begin, tets_end);
2445 32 tets_begin_ = tets_begin;
2446 32 tets_end_ = tets_end;
2447 32 }
2448
2449 230 void delete_threads() override {
2450
4/4
✓ Branch 0 taken 12 times.
✓ Branch 1 taken 103 times.
✓ Branch 2 taken 48 times.
✓ Branch 3 taken 12 times.
 326 delete[] parts_;
2451 230 parts_ = nullptr;
2452 230 nb_parts_ = 0;
2453 230 }
2454
2455 /**
2456 * \brief Gets the number of parts (or number of threads).
2457 */
2458 28526 index_t nb_parts() const {
2459 28526 return nb_parts_;
2460 }
2461
2462 /**
2463 * \brief Gets a given part from its index.
2464 * \param[in] i index of the part
2465 * \pre \p i < nb_parts()
2466 */
2467 10760 thisclass& part(index_t i) {
2468
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 5380 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 10760 geo_debug_assert(i < nb_parts());
2469 10760 return parts_[i];
2470 }
2471
2472 /**
2473 * \copydoc RestrictedVoronoiDiagram::point_allocator()
2474 */
2475 ✗ GEOGen::PointAllocator* point_allocator() override {
2476 ✗ return RVD_.point_allocator();
2477 }
2478
2479
2480 protected:
2481
2482 GenRestrictedVoronoiDiagram RVD_;
2483
2484 // For projection
2485 bool use_exact_projection_;
2486 index_t nb_triangles_;
2487 vector<index_t> triangles_;
2488 vector<vector<index_t> > stars_;
2489 Delaunay_var mesh_vertices_;
2490
2491 // One of MT_NONE, MT_LLOYD, MT_NEWTON
2492 ThreadMode thread_mode_;
2493
2494 bool is_slave_;
2495
2496 // Variables for 'master' in multithreading mode
2497 thisclass* parts_;
2498 index_t nb_parts_;
2499 Process::SpinLockArray spinlocks_;
2500
2501 // Newton mode with int. simplex
2502 IntegrationSimplex* simplex_func_;
2503
2504 // PolygonCallback mode.
2505 RVDPolygonCallback* polygon_callback_;
2506
2507 // PolyhedronCallback mode.
2508 RVDPolyhedronCallback* polyhedron_callback_;
2509
2510 // master stores argument for compute_centroids() and
2511 // compute_CVT_func_grad() to pass it to the parts.
2512 double* arg_vectors_;
2513 double* arg_scalars_;
2514
2515 // Variables for 'slaves' in multithreading mode
2516 thisclass* master_;
2517 double funcval_; // Newton mode: function value
2518
2519 protected:
2520 /**
2521 * \brief Destructor
2522 */
2523 288 ~RVD_Nd_Impl() override {
2524 192 delete_threads();
2525 288 }
2526
2527 private:
2528 /** \brief Forbids construction by copy. */
2529 RVD_Nd_Impl(const thisclass&);
2530
2531 /** \brief Forbids assignment. */
2532 thisclass& operator= (const thisclass&);
2533 };
2534 }
2535
2536 /****************************************************************************/
2537
2538 namespace GEO {
2539
2540 48 RestrictedVoronoiDiagram* RestrictedVoronoiDiagram::create(
2541 Delaunay* delaunay, Mesh* mesh,
2542 const double* R3_embedding, index_t R3_embedding_stride
2543 ) {
2544
2545 48 geo_cite("DBLP:journals/tog/EdelsbrunnerM90");
2546 48 geo_cite("DBLP:conf/compgeom/Shewchuk96");
2547 48 geo_cite("meyer:inria-00344297");
2548 48 geo_cite("DBLP:conf/gmp/YanWLL10");
2549 48 geo_cite("DBLP:journals/cad/YanWLL13");
2550 48 geo_cite("DBLP:journals/cad/Levy16");
2551
2552 48 delaunay->set_stores_neighbors(true);
2553 48 RestrictedVoronoiDiagram* result = nullptr;
2554
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 48 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 48 geo_assert(delaunay != nullptr);
2555 48 coord_index_t dim = delaunay->dimension();
2556
4/8
✗ Branch 0 not taken.
✓ Branch 1 taken 40 times.
✓ Branch 2 taken 1 times.
✓ Branch 3 taken 6 times.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 48 switch(dim) {
2557 ✗ case 2:
2558 ✗ result = new RVD_Nd_Impl<2>(
2559 delaunay, mesh, R3_embedding, R3_embedding_stride
2560 ✗ );
2561 ✗ break;
2562 40 case 3:
2563
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 40 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 40 result = new RVD_Nd_Impl<3>(
2564 delaunay, mesh, R3_embedding, R3_embedding_stride
2565
1/2
✓ Branch 2 taken 40 times.
✗ Branch 3 not taken.
 40 );
2566 40 break;
2567 1 case 4:
2568
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 1 result = new RVD_Nd_Impl<4>(
2569 delaunay, mesh, R3_embedding, R3_embedding_stride
2570
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 );
2571 1 break;
2572 6 case 6:
2573
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 6 result = new RVD_Nd_Impl<6>(
2574 delaunay, mesh, R3_embedding, R3_embedding_stride
2575
1/2
✓ Branch 2 taken 6 times.
✗ Branch 3 not taken.
 6 );
2576 6 break;
2577 1 case 8:
2578
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 1 result = new RVD_Nd_Impl<8>(
2579 delaunay, mesh, R3_embedding, R3_embedding_stride
2580
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 );
2581 1 break;
2582 ✗ case 20:
2583 ✗ result = new RVD_Nd_Impl<20>(
2584 delaunay, mesh, R3_embedding, R3_embedding_stride
2585 ✗ );
2586 ✗ break;
2587 ✗ case 100:
2588 ✗ result = new RVD_Nd_Impl<100>(
2589 delaunay, mesh, R3_embedding, R3_embedding_stride
2590 ✗ );
2591 ✗ break;
2592 ✗ default:
2593 ✗ geo_assert_not_reached;
2594 }
2595
5/8
✓ Branch 1 taken 48 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 48 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 48 times.
✗ Branch 8 not taken.
✓ Branch 11 taken 35 times.
✓ Branch 12 taken 13 times.
 96 if(CmdLine::get_arg("algo:predicates") == "exact") {
2596 35 result->set_exact_predicates(true);
2597 }
2598 48 return result;
2599 }
2600
2601 144 void RestrictedVoronoiDiagram::set_delaunay(Delaunay* delaunay) {
2602 144 delaunay_ = delaunay;
2603
2/2
✓ Branch 0 taken 96 times.
✓ Branch 1 taken 48 times.
 144 if(delaunay_ != nullptr) {
2604 96 dimension_ = delaunay->dimension();
2605 } else {
2606 48 dimension_ = 0;
2607 }
2608 144 }
2609
2610 192 RestrictedVoronoiDiagram::~RestrictedVoronoiDiagram() {
2611 192 }
2612
2613 96 RestrictedVoronoiDiagram::RestrictedVoronoiDiagram(
2614 Delaunay* delaunay, Mesh* mesh,
2615 const double* R3_embedding, index_t R3_embedding_stride
2616 96 ) :
2617 96 dimension_(0),
2618 96 mesh_(mesh),
2619 96 R3_embedding_base_(R3_embedding),
2620
1/2
✓ Branch 2 taken 96 times.
✗ Branch 3 not taken.
 96 R3_embedding_stride_(R3_embedding_stride) {
2621
1/2
✓ Branch 1 taken 96 times.
✗ Branch 2 not taken.
 96 set_delaunay(delaunay);
2622 96 has_weights_ = false;
2623 96 facets_begin_ = NO_INDEX;
2624 96 facets_end_ = NO_INDEX;
2625 96 tets_begin_ = NO_INDEX;
2626 96 tets_end_ = NO_INDEX;
2627 96 volumetric_ = false;
2628 96 }
2629
2630
2631 ✗ void RestrictedVoronoiDiagram::compute_RDT(
2632 Mesh& RDT,
2633 RDTMode mode,
2634 const vector<bool>& seed_is_locked,
2635 MeshFacetsAABB* AABB
2636 ) {
2637 ✗ vector<index_t> simplices;
2638 ✗ vector<double> embedding;
2639 ✗ compute_RDT(
2640 simplices, embedding,
2641 mode, seed_is_locked,
2642 AABB
2643 );
2644 ✗ if(volumetric()) {
2645 ✗ RDT.cells.assign_tet_mesh(dimension(),embedding,simplices,true);
2646 } else {
2647 ✗ RDT.facets.assign_triangle_mesh(
2648 ✗ dimension(),embedding,simplices,true
2649 );
2650 ✗ if((mode & RDT_DONT_REPAIR) == 0) {
2651 ✗ mesh_repair(RDT); // Needed to reorient triangles
2652 }
2653 }
2654 ✗ }
2655
2656
2657 }
2658