GCC Code Coverage Report

Directory: ./
File: lib/geogram/voronoi/CVT.cpp
Date: 2026-09-07 02:37:58
Exec Total Coverage
Lines: 130 208 62.5%
Functions: 11 16 68.8%
Branches: 82 268 30.6%
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/CVT.h>
41 #include <geogram/voronoi/RVD.h>
42 #include <geogram/mesh/mesh_repair.h>
43 #include <geogram/mesh/mesh_geometry.h>
44 #include <geogram/numerics/optimizer.h>
45 #include <geogram/basic/progress.h>
46 #include <geogram/basic/argused.h>
47 #include <geogram/bibliography/bibliography.h>
48
49 /****************************************************************************/
50
51 namespace GEO {
52
53 CentroidalVoronoiTesselation*
54 CentroidalVoronoiTesselation::instance_ = nullptr;
55
56 10 CentroidalVoronoiTesselation::CentroidalVoronoiTesselation(
57 Mesh* mesh, coord_index_t dim, const std::string& delaunay
58 10 ) {
59 10 use_RVC_centroids_ = true;
60 10 show_iterations_ = false;
61 10 constrained_cvt_ = false;
62 10 dimension_ =
63
3/4
✓ Branch 0 taken 6 times.
✓ Branch 1 taken 4 times.
✓ Branch 3 taken 6 times.
✗ Branch 4 not taken.
 10 (dim != 0) ? dim : coord_index_t(mesh->vertices.dimension());
64
2/8
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 10 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
 10 geo_assert(index_t(dimension_) <= mesh->vertices.dimension());
65 10 is_projection_ = true;
66
2/4
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 10 times.
✗ Branch 5 not taken.
 10 delaunay_ = Delaunay::create(dimension_, delaunay);
67
2/4
✓ Branch 2 taken 10 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 10 times.
✗ Branch 6 not taken.
 10 RVD_ = RestrictedVoronoiDiagram::create(delaunay_, mesh);
68 10 mesh_ = mesh;
69
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 10 geo_assert(instance_ == nullptr);
70 10 instance_ = this;
71 10 progress_ = nullptr;
72
1/2
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
 10 geo_cite("Lloyd82leastsquares");
73
1/2
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
 10 geo_cite("Du:1999:CVT:340312.340319");
74
1/2
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
 10 geo_cite("DBLP:journals/tog/LiuWLSYLY09");
75 10 }
76
77 ✗ CentroidalVoronoiTesselation::CentroidalVoronoiTesselation(
78 Mesh* mesh, const vector<vec3>& R3_embedding, coord_index_t dim,
79 const std::string& delaunay
80 ✗ ) {
81 ✗ use_RVC_centroids_ = true;
82 ✗ show_iterations_ = false;
83 ✗ constrained_cvt_ = false;
84 ✗ dimension_ =
85 ✗ (dim != 0) ? dim : coord_index_t(mesh->vertices.dimension());
86 ✗ geo_assert(index_t(dimension_) <= mesh->vertices.dimension());
87 ✗ is_projection_ = (R3_embedding.size() == 0);
88 ✗ delaunay_ = Delaunay::create(dimension_, delaunay);
89 ✗ if(is_projection_) {
90 ✗ RVD_ = RestrictedVoronoiDiagram::create(delaunay_, mesh);
91 } else {
92 RVD_ = RestrictedVoronoiDiagram::create(
93 delaunay_, mesh, R3_embedding
94 ✗ );
95 }
96 ✗ mesh_ = mesh;
97 ✗ geo_assert(instance_ == nullptr);
98 ✗ instance_ = this;
99 ✗ progress_ = nullptr;
100 ✗ geo_cite("Lloyd82leastsquares");
101 ✗ geo_cite("Du:1999:CVT:340312.340319");
102 ✗ geo_cite("DBLP:journals/tog/LiuWLSYLY09");
103 ✗ }
104
105 20 CentroidalVoronoiTesselation::~CentroidalVoronoiTesselation() {
106 20 instance_ = nullptr;
107 20 }
108
109 10 bool CentroidalVoronoiTesselation::compute_initial_sampling(
110 index_t nb_samples, bool verbose
111 ) {
112 10 points_.resize(dimension_ * nb_samples);
113 10 return RVD_->compute_initial_sampling(
114 points_.data(), nb_samples, verbose
115 10 );
116 }
117
118 ✗ void CentroidalVoronoiTesselation::set_points(
119 index_t nb_points, const double* points
120 ) {
121 ✗ points_.resize(dimension_ * nb_points);
122 ✗ for(index_t i = 0; i < points_.size(); i++) {
123 ✗ points_[i] = points[i];
124 }
125 ✗ }
126
127 ✗ void CentroidalVoronoiTesselation::resize_points(
128 index_t nb_points
129 ) {
130 ✗ points_.resize(dimension_ * nb_points);
131 ✗ }
132
133 10 void CentroidalVoronoiTesselation::Lloyd_iterations(index_t nb_iter) {
134 10 index_t nb_points = index_t(points_.size() / dimension_);
135
136 10 vector<double> mg;
137 10 vector<double> m;
138
139
2/4
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 10 times.
✗ Branch 5 not taken.
 10 RVD_->set_check_SR(false);
140
141
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
 10 if(progress_ != nullptr) {
142
1/2
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
 10 progress_->reset(nb_iter);
143 }
144
145 10 cur_iter_ = 0;
146 10 nb_iter_ = nb_iter;
147
148
2/2
✓ Branch 0 taken 50 times.
✓ Branch 1 taken 10 times.
 60 for(index_t i = 0; i < nb_iter; i++) {
149
1/2
✓ Branch 1 taken 50 times.
✗ Branch 2 not taken.
 50 mg.assign(nb_points * dimension_, 0.0);
150
1/2
✓ Branch 1 taken 50 times.
✗ Branch 2 not taken.
 50 m.assign(nb_points, 0.0);
151
2/4
✓ Branch 1 taken 50 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 50 times.
✗ Branch 6 not taken.
 50 delaunay_->set_vertices(nb_points, points_.data());
152
2/4
✓ Branch 1 taken 50 times.
✗ Branch 2 not taken.
✓ Branch 6 taken 50 times.
✗ Branch 7 not taken.
 50 RVD_->compute_centroids(mg.data(), m.data());
153 50 index_t cur = 0;
154
2/2
✓ Branch 0 taken 156000 times.
✓ Branch 1 taken 50 times.
 156050 for(index_t j = 0; j < nb_points; j++) {
155
7/10
✓ Branch 1 taken 156000 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 155966 times.
✓ Branch 4 taken 34 times.
✓ Branch 6 taken 155966 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 155966 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 155966 times.
✓ Branch 11 taken 34 times.
 156000 if(m[j] > 1e-30 && !point_is_locked(j)) {
156
1/2
✓ Branch 1 taken 155966 times.
✗ Branch 2 not taken.
 155966 double s = 1.0 / m[j];
157
2/2
✓ Branch 0 taken 781296 times.
✓ Branch 1 taken 155966 times.
 937262 for(index_t coord = 0; coord < dimension_; coord++) {
158
2/4
✓ Branch 1 taken 781296 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 781296 times.
✗ Branch 5 not taken.
 781296 points_[cur + coord] = s * mg[cur + coord];
159 }
160 }
161 156000 cur += dimension_;
162 }
163
1/2
✓ Branch 1 taken 50 times.
✗ Branch 2 not taken.
 50 newiteration();
164 }
165
166 10 progress_ = nullptr;
167 10 }
168
169 6 void CentroidalVoronoiTesselation::compute_surface(
170 Mesh* mesh, bool multinerve
171 ) {
172 6 index_t nb_points = index_t(points_.size() / dimension_);
173
2/4
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 6 times.
✗ Branch 6 not taken.
 6 delaunay_->set_vertices(nb_points, points_.data());
174
175 6 vector<index_t> triangles;
176 6 vector<double> vertices;
177 6 vector<double> vertices_R3;
178
179 6 RestrictedVoronoiDiagram::RDTMode mode =
180 RestrictedVoronoiDiagram::RDTMode(0);
181
182
1/2
✓ Branch 0 taken 6 times.
✗ Branch 1 not taken.
 6 if(multinerve) {
183 6 mode = RestrictedVoronoiDiagram::RDTMode(
184 6 mode | RestrictedVoronoiDiagram::RDT_MULTINERVE
185 );
186 }
187
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 2 times.
 6 if(use_RVC_centroids_) {
188 4 mode = RestrictedVoronoiDiagram::RDTMode(
189 4 mode | RestrictedVoronoiDiagram::RDT_RVC_CENTROIDS
190 4 | RestrictedVoronoiDiagram::RDT_PREFER_SEEDS
191 );
192 }
193
194
2/4
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6 times.
✗ Branch 5 not taken.
 6 RVD_->set_check_SR(true);
195
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 RVD_->compute_RDT(
196
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 triangles, vertices, mode, point_is_locked_
197 );
198
199 // TODO: projection is not good when the embedding is not one-to-one,
200 // for instance with Gauss map. We should use barycentric coordinates
201 // instead.
202 6 index_t nb_vertices = index_t(vertices.size() / dimension_);
203
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 vertices_R3.resize(nb_vertices * 3);
204
1/2
✓ Branch 0 taken 6 times.
✗ Branch 1 not taken.
 6 if(is_projection_) {
205 6 double* cur = vertices.data();
206
2/2
✓ Branch 0 taken 30000 times.
✓ Branch 1 taken 6 times.
 30006 for(index_t v = 0; v < nb_vertices; v++) {
207
1/2
✓ Branch 1 taken 30000 times.
✗ Branch 2 not taken.
 30000 vertices_R3[3 * v] = cur[0];
208
1/2
✓ Branch 1 taken 30000 times.
✗ Branch 2 not taken.
 30000 vertices_R3[3 * v + 1] = cur[1];
209
1/2
✓ Branch 1 taken 30000 times.
✗ Branch 2 not taken.
 30000 vertices_R3[3 * v + 2] = cur[2];
210 30000 cur += dimension_;
211 }
212 } else {
213 ✗ RVD_->project_points_on_surface(
214 nb_vertices, vertices.data(), vertices_R3
215 );
216 }
217
218
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 mesh->clear();
219
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 mesh->facets.assign_triangle_mesh(3, vertices_R3, triangles, true);
220
221
1/2
✓ Branch 0 taken 6 times.
✗ Branch 1 not taken.
 6 if(multinerve) {
222
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 mesh_postprocess_RDT(*mesh);
223 } else {
224 // The 'repair' phase is needed to reconstruct the
225 // facet-facet links, that are not initialized by
226 // Mesh::assign_triangle_mesh()
227 ✗ double radius = bbox_diagonal(*mesh);
228 ✗ mesh_repair(*mesh, MESH_REPAIR_DEFAULT, 1e-6 * radius);
229 // TODO: check: is it really good to have some tolerance here,
230 // not sure, may cause some Moebius configs sometimes.
231 }
232 6 }
233
234 ✗ void CentroidalVoronoiTesselation::compute_volume(
235 Mesh* mesh
236 ) {
237 ✗ geo_assert(volumetric());
238 ✗ index_t nb_points = index_t(points_.size() / dimension_);
239 ✗ delaunay_->set_vertices(nb_points, points_.data());
240
241 ✗ vector<index_t> tets;
242 ✗ vector<double> vertices;
243 ✗ vector<double> vertices_R3;
244
245 ✗ RVD_->set_check_SR(true);
246 ✗ RVD_->compute_RDT(
247 tets, vertices
248 );
249
250 // TODO: projection is not good when the embedding is not one-to-one,
251 // for instance with Gauss map. We should use barycentric coordinates
252 // instead.
253 ✗ index_t nb_vertices = index_t(vertices.size() / dimension_);
254 ✗ vertices_R3.resize(nb_vertices * 3);
255 ✗ if(is_projection_) {
256 ✗ double* cur = vertices.data();
257 ✗ for(index_t v = 0; v < nb_vertices; v++) {
258 ✗ vertices_R3[3 * v] = cur[0];
259 ✗ vertices_R3[3 * v + 1] = cur[1];
260 ✗ vertices_R3[3 * v + 2] = cur[2];
261 ✗ cur += dimension_;
262 }
263 } else {
264 // TODO: map vertices from embedding space to 3D space
265 // when we are not in projection mode
266 ✗ geo_assert_not_reached;
267 }
268 ✗ mesh->clear();
269 ✗ mesh->cells.assign_tet_mesh(3, vertices_R3, tets, true);
270 ✗ }
271
272 10 void CentroidalVoronoiTesselation::Newton_iterations(
273 index_t nb_iter, index_t m
274 ) {
275
2/4
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 10 times.
✗ Branch 5 not taken.
 10 Optimizer_var optimizer = Optimizer::create("HLBFGS");
276
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 10 times.
 10 if(optimizer.is_null()) {
277 ✗ Logger::warn("CVT") << "This geogram was not compiled with HLBFGS"
278 ✗ << " (falling back to Lloyd iterations)"
279 ✗ << std::endl;
280 ✗ Lloyd_iterations(nb_iter);
281 ✗ return;
282 }
283
284 10 index_t n = index_t(points_.size());
285
286
2/4
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 10 times.
✗ Branch 5 not taken.
 10 RVD_->set_check_SR(true);
287
288
1/2
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
 10 if(progress_ != nullptr) {
289
1/2
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
 10 progress_->reset(nb_iter);
290 }
291
292 10 cur_iter_ = 0;
293 10 nb_iter_ = nb_iter;
294
295
1/2
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
 10 optimizer->set_epsg(0.0);
296
1/2
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
 10 optimizer->set_epsf(0.0);
297
1/2
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
 10 optimizer->set_epsx(0.0);
298
1/2
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
 10 optimizer->set_newiteration_callback(newiteration_CB);
299
1/2
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
 10 optimizer->set_funcgrad_callback(funcgrad_CB);
300
1/2
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
 10 optimizer->set_N(n);
301
1/2
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
 10 optimizer->set_M(m);
302
1/2
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
 10 optimizer->set_max_iter(nb_iter);
303
2/4
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 10 times.
✗ Branch 6 not taken.
 10 optimizer->optimize(points_.data());
304
305
1/2
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
 10 simplex_func_.reset();
306 10 progress_ = nullptr;
307
1/2
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
 10 }
308
309 335 void CentroidalVoronoiTesselation::constrain_points(double* g) const {
310
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 335 times.
 335 if(point_is_locked_.size() != 0) {
311 ✗ double* cur_g = g;
312 ✗ for(index_t i = 0; i < nb_points(); ++i) {
313 ✗ if(point_is_locked_[i]) {
314 ✗ for(index_t c = 0; c < dimension_; c++) {
315 ✗ cur_g[c] = 0.0;
316 }
317 }
318 ✗ cur_g += dimension_;
319 }
320 }
321 335 }
322
323 335 void CentroidalVoronoiTesselation::funcgrad(
324 index_t n, double* x, double& f, double* g
325 ) {
326 335 index_t nb_points = n / dimension_;
327 335 delaunay_->set_vertices(nb_points, x);
328 335 Memory::clear(g, n * sizeof(double));
329 335 f = 0.0;
330
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 335 times.
 335 if(!simplex_func_.is_null()) {
331 ✗ RVD_->compute_integration_simplex_func_grad(
332 f,g,simplex_func_
333 );
334 } else {
335 335 RVD_->compute_CVT_func_grad(f, g);
336 }
337 335 constrain_points(g);
338 335 }
339
340 360 void CentroidalVoronoiTesselation::newiteration() {
341
1/2
✓ Branch 0 taken 360 times.
✗ Branch 1 not taken.
 360 if(progress_ != nullptr) {
342 360 progress_->next();
343 }
344 360 cur_iter_++;
345 360 }
346
347 335 void CentroidalVoronoiTesselation::funcgrad_CB(
348 index_t n, double* x, double& f, double* g
349 ) {
350 335 instance_->funcgrad(n, x, f, g);
351 335 }
352
353 310 void CentroidalVoronoiTesselation::newiteration_CB(
354 index_t n, const double* x, double f, const double* g, double gnorm
355 ) {
356 310 geo_argused(n);
357 310 geo_argused(x);
358 310 geo_argused(f);
359 310 geo_argused(g);
360 310 geo_argused(gnorm);
361 310 instance_->newiteration();
362 310 }
363
364 ✗ void CentroidalVoronoiTesselation::compute_R3_embedding() {
365 ✗ index_t nb_points = index_t(points_.size() / dimension_);
366 ✗ points_R3_.resize(nb_points);
367 ✗ if(is_projection_ && !constrained_cvt_) {
368 ✗ double* cur = points_.data();
369 ✗ for(index_t p = 0; p < nb_points; p++) {
370 ✗ points_R3_[p] = vec3(cur[0], cur[1], cur[2]);
371 ✗ cur += dimension_;
372 }
373 ✗ } else {
374 ✗ RVD_->project_points_on_surface(
375 ✗ nb_points, points_.data(), points_R3_, constrained_cvt_
376 );
377 }
378 ✗ }
379 }
380