GCC Code Coverage Report

Directory: ./
File: examples/geogram/compute_delaunay/main.cpp
Date: 2026-09-07 02:25:23
Exec Total Coverage
Lines: 93 98 94.9%
Functions: 3 3 100.0%
Branches: 123 250 49.2%
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/basic/common.h>
41 #include <geogram/basic/logger.h>
42 #include <geogram/basic/command_line.h>
43 #include <geogram/basic/command_line_args.h>
44 #include <geogram/basic/stopwatch.h>
45 #include <geogram/basic/file_system.h>
46 #include <geogram/mesh/mesh.h>
47 #include <geogram/mesh/mesh_io.h>
48 #include <geogram/mesh/mesh_reorder.h>
49 #include <geogram/delaunay/delaunay.h>
50 #include <algorithm>
51
52 namespace {
53 using namespace GEO;
54
55 /**
56 * \brief Loads points from a file.
57 * \param[in] points_filename the name of the file with the points.
58 * -If the example was compiled with the Geogram library, then any
59 * mesh file handled by Geogram can be used.
60 * -if the example was compiled with Delaunay_psm (single file), then
61 * the file should be ASCII, with one point per line.
62 * \param[in] dimension number of coordinates of the points.
63 * \param[out] points the loaded points, in a single vector of coordinates.
64 * In the end, the number of loaded points is points.size()/dimension.
65 */
66 9 bool load_points(
67 const std::string& points_filename,
68 index_t dimension,
69 vector<double>& points
70 ) {
71 #ifdef GEOGRAM_PSM
72 // Simple data input: one point per line, coordinates in ASCII
73 LineInput input(points_filename);
74 if(!input.OK()) {
75 return false;
76 }
77 while(!input.eof() && input.get_line()) {
78 input.get_fields();
79 if(input.nb_fields() == dimension) {
80 for(index_t c=0; c<dimension; ++c) {
81 points.push_back(input.field_as_double(c));
82 }
83 }
84 }
85 #else
86 // Using Geogram mesh I/O
87 9 Mesh M;
88
1/2
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
 9 MeshIOFlags flags;
89 flags.reset_element(MESH_FACETS);
90 flags.reset_element(MESH_CELLS);
91
2/4
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 9 times.
✗ Branch 4 not taken.
 9 if(!mesh_load(points_filename, M, flags)) {
92 return false;
93 }
94
1/2
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
 9 M.vertices.set_dimension(dimension);
95 index_t nb_points = M.vertices.nb();
96
1/2
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
 9 points.resize(nb_points * dimension);
97 9 Memory::copy(
98 points.data(),
99 M.vertices.point_ptr(0),
100 9 M.vertices.nb()*dimension*sizeof(double)
101 );
102 #endif
103 9 return true;
104 9 }
105
106 /**
107 * \brief Saves a Delaunay triangulation to a file.
108 * \param[in] delaunay a pointer to the Delaunay triangulation.
109 * \param[in] filename the name of the file to be saved.
110 * -If the example was compiled with the Geogram library, then any
111 * mesh file handled by Geogram can be used.
112 * if the example was compiled with Delaunay_psm (single file), then
113 * the points and vertices of the triangulation are output in ASCII.
114 * \param[in] convex_hull_only if true, then only the triangles on the
115 * convex hull are output.
116 */
117
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 7 times.
 9 void save_Delaunay(
118 Delaunay* delaunay, const std::string& filename,
119 bool convex_hull_only = false
120 ) {
121 vector<index_t> tri2v;
122
123
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 7 times.
 9 if(convex_hull_only) {
124
125 // The convex hull can be efficiently traversed only if infinite
126 // tetrahedra are kept.
127
1/8
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
 2 geo_assert(delaunay->keeps_infinite());
128
129 // The convex hull can be retrieved as the finite facets
130 // of the infinite cells (note: it would be also possible to
131 // throw away the infinite cells and get the convex hull as
132 // the facets adjacent to no cell). Here we use the infinite
133 // cells to show an example with them.
134
135
136 // This block is just a sanity check
137 {
138
2/2
✓ Branch 0 taken 779 times.
✓ Branch 1 taken 2 times.
 781 for(index_t t=0; t < delaunay->nb_finite_cells(); ++t) {
139 geo_debug_assert(delaunay->cell_is_finite(t));
140 }
141
142 for(index_t t=delaunay->nb_finite_cells();
143 t < delaunay->nb_cells(); ++t) {
144 geo_debug_assert(delaunay->cell_is_infinite(t));
145 }
146 }
147
148 // This iterates on the infinite cells
149 366 for(
150 index_t t = delaunay->nb_finite_cells();
151
2/2
✓ Branch 0 taken 366 times.
✓ Branch 1 taken 2 times.
 368 t < delaunay->nb_cells(); ++t
152 ) {
153
2/2
✓ Branch 0 taken 1464 times.
✓ Branch 1 taken 366 times.
 1830 for(index_t lv=0; lv<4; ++lv) {
154 index_t v = delaunay->cell_vertex(t,lv);
155
2/2
✓ Branch 0 taken 1098 times.
✓ Branch 1 taken 366 times.
 1464 if(v != NO_INDEX) {
156
1/4
✓ Branch 1 taken 1098 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 1098 tri2v.push_back(index_t(v));
157 }
158 }
159 }
160 }
161
162 #ifdef GEOGRAM_PSM
163 // Simple data output: output vertices and simplices
164
165 Logger::out("Delaunay") << "Saving output to " << filename << std::endl;
166 std::ofstream out(filename.c_str());
167
168 out << delaunay->nb_vertices() << " vertices" << std::endl;
169 for(index_t v=0; v < delaunay->nb_vertices(); ++v) {
170 for(index_t c=0; c < delaunay->dimension(); ++c) {
171 out << delaunay->vertex_ptr(v)[c] << " ";
172 }
173 out << std::endl;
174 }
175 if(convex_hull_only) {
176 out << tri2v.size()/3 << " simplices" << std::endl;
177 for(index_t t=0; t<tri2v.size()/3; ++t) {
178 out << tri2v[3*t] << " "
179 << tri2v[3*t+1] << " "
180 << tri2v[3*t+2] << std::endl;
181 }
182 } else {
183 out << delaunay->nb_cells() << " simplices" << std::endl;
184 for(index_t t=0; t<delaunay->nb_cells(); ++t) {
185 for(index_t lv=0; lv<delaunay->cell_size(); ++lv) {
186 out << delaunay->cell_vertex(t,lv) << " ";
187 }
188 out << std::endl;
189 }
190 }
191
192 #else
193 // Using Geogram mesh I/O: copy Delaunay into a Geogram
194 // mesh and save it to disk.
195
196
1/2
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
 9 Mesh M_out;
197
1/2
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
 9 vector<double> pts(delaunay->nb_vertices() * 3);
198
2/2
✓ Branch 0 taken 1137 times.
✓ Branch 1 taken 9 times.
 1146 for(index_t v = 0; v < delaunay->nb_vertices(); ++v) {
199
2/2
✓ Branch 0 taken 687 times.
✓ Branch 1 taken 450 times.
 1137 pts[3 * v] = delaunay->vertex_ptr(v)[0];
200
2/2
✓ Branch 0 taken 687 times.
✓ Branch 1 taken 450 times.
 1137 pts[3 * v + 1] = delaunay->vertex_ptr(v)[1];
201 1137 pts[3 * v + 2] =
202
2/2
✓ Branch 0 taken 687 times.
✓ Branch 1 taken 450 times.
 1137 (delaunay->dimension() >= 3) ? delaunay->vertex_ptr(v)[2] : 0.0;
203 }
204
205
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 7 times.
 9 if(convex_hull_only) {
206
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 M_out.facets.assign_triangle_mesh(3, pts, tri2v, true);
207
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 3 times.
 7 } else if(delaunay->dimension() == 3) {
208
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 vector<index_t> tet2v(delaunay->nb_cells() * 4);
209
2/2
✓ Branch 0 taken 1558 times.
✓ Branch 1 taken 4 times.
 1562 for(index_t t = 0; t < delaunay->nb_cells(); ++t) {
210 1558 tet2v[4 * t] = index_t(delaunay->cell_vertex(t, 0));
211 1558 tet2v[4 * t + 1] = index_t(delaunay->cell_vertex(t, 1));
212 1558 tet2v[4 * t + 2] = index_t(delaunay->cell_vertex(t, 2));
213 1558 tet2v[4 * t + 3] = index_t(delaunay->cell_vertex(t, 3));
214 }
215
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 M_out.cells.assign_tet_mesh(3, pts, tet2v, true);
216
1/2
✓ Branch 0 taken 3 times.
✗ Branch 1 not taken.
 3 } else if(delaunay->dimension() == 2) {
217
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 3 tri2v.resize(delaunay->nb_cells() * 3);
218
2/2
✓ Branch 0 taken 153 times.
✓ Branch 1 taken 3 times.
 156 for(index_t t = 0; t < delaunay->nb_cells(); ++t) {
219 153 tri2v[3 * t] = index_t(delaunay->cell_vertex(t, 0));
220 153 tri2v[3 * t + 1] = index_t(delaunay->cell_vertex(t, 1));
221 153 tri2v[3 * t + 2] = index_t(delaunay->cell_vertex(t, 2));
222 }
223
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 3 M_out.facets.assign_triangle_mesh(3, pts, tri2v, true);
224 }
225
1/2
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
 9 M_out.show_stats();
226
227
2/6
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 9 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 9 Logger::div("Saving the result");
228
1/2
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
 9 MeshIOFlags flags;
229 flags.set_element(MESH_FACETS);
230 flags.set_element(MESH_CELLS);
231
1/2
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
 9 mesh_save(M_out, filename, flags);
232 #endif
233 18 }
234
235 }
236
237 9 int main(int argc, char** argv) {
238 using namespace GEO;
239
240 // Needs to be called once.
241 9 GEO::initialize(GEO::GEOGRAM_INSTALL_ALL);
242
243 try {
244
245
3/6
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 9 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 9 times.
✗ Branch 8 not taken.
 18 Stopwatch Wtot("Total time");
246
247 std::vector<std::string> filenames;
248
249
2/4
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 9 times.
✗ Branch 5 not taken.
 9 CmdLine::import_arg_group("standard");
250
2/4
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 9 times.
✗ Branch 5 not taken.
 9 CmdLine::import_arg_group("algo");
251
252
1/2
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
 9 CmdLine::declare_arg(
253
2/6
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 9 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
 18 "convex_hull", false,
254
1/2
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
 9 "compute just the convex hull of the points"
255 );
256
257
1/2
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
 9 CmdLine::declare_arg(
258
3/8
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 9 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 9 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
 18 "dimension", 3, "3 for 3D, 2 for 2D"
259 );
260
261
2/4
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 9 times.
✗ Branch 5 not taken.
 9 CmdLine::set_arg("algo:delaunay","default");
262
263 9 if(
264
1/2
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
 9 !CmdLine::parse(
265
2/4
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 9 times.
✗ Branch 4 not taken.
 18 argc, argv, filenames, "pointsfile <outputfile|none>"
266 )
267 ) {
268 return 1;
269 }
270
271
272 std::string points_filename = filenames[0];
273
274 std::string output_filename =
275
3/8
✗ Branch 0 not taken.
✓ Branch 1 taken 9 times.
✓ Branch 3 taken 9 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 9 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
 9 filenames.size() >= 2 ? filenames[1] : std::string("out.mesh");
276
277 bool output = (output_filename != "none");
278
279
2/4
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 9 times.
✗ Branch 5 not taken.
 9 Logger::div("Data I/O");
280
281
2/4
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 9 times.
✗ Branch 5 not taken.
 9 Logger::out("I/O") << "Output = " << output_filename << std::endl;
282
283
2/4
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 9 times.
✗ Branch 5 not taken.
 9 bool convex_hull_only = CmdLine::get_arg_bool("convex_hull");
284
3/6
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 9 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 9 times.
✗ Branch 7 not taken.
 9 index_t dimension = index_t(CmdLine::get_arg_int("dimension"));
285
286
2/6
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 9 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 9 std::string del = CmdLine::get_arg("algo:delaunay");
287
2/2
✓ Branch 0 taken 6 times.
✓ Branch 1 taken 3 times.
 9 if(del == "default") {
288
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 2 times.
 6 if(dimension == 3) {
289
3/6
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 4 times.
✗ Branch 7 not taken.
 8 if(DelaunayFactory::has_creator("PDEL")) {
290 // PDEL = Parallel 3D Delaunay
291
2/4
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
 8 CmdLine::set_arg("algo:delaunay", "PDEL");
292 } else {
293 // BDEL = Sequential 3D Delaunay
294 ✗ CmdLine::set_arg("algo:delaunay", "BDEL");
295 }
296
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 } else if(dimension == 2) {
297 // BDEL2d = Sequential 2D Delaunay
298
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 4 CmdLine::set_arg("algo:delaunay", "BDEL2d");
299 }
300 }
301
302
3/8
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 9 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 9 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
 18 Logger::out("Delaunay")
303
3/8
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 9 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 9 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
 18 << "Using " << CmdLine::get_arg("algo:delaunay") << std::endl;
304
305 // Note: To create a parallel Delaunay 3D, one can use directly:
306 // Delaunay_var delaunay = Delaunay::create(3,"PDEL") instead
307 // of the line below (that uses the command line to select the
308 // implementation of Delaunay).
309
1/2
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
 9 Delaunay_var delaunay = Delaunay::create(coord_index_t(dimension));
310
311 // If we want the convex hull, we keep the infinite facets,
312 // because the convex hull can be retreived as the finite facets
313 // of the infinite cells (note: it would be also possible to
314 // throw away the infinite cells and get the convex hull as
315 // the facets adjacent to no cell).
316
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 7 times.
 9 if(convex_hull_only) {
317
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 delaunay->set_keeps_infinite(true);
318 }
319
320 vector<double> points;
321
322
2/4
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 9 times.
 9 if(!load_points(points_filename, dimension, points)) {
323 ✗ Logger::err("Delaunay") << "Could not load points" << std::endl;
324 return 1;
325 }
326
327 9 index_t nb_points = points.size() / dimension;
328
329
2/4
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 9 times.
✗ Branch 5 not taken.
 9 Logger::out("Delaunay")
330 << "Loaded " << nb_points << " points" << std::endl;
331
332 double time = 0.0;
333 {
334
2/4
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 9 times.
✗ Branch 5 not taken.
 9 Stopwatch Wdel("Delaunay");
335 // Note: this does not transfer ownership of memory, caller
336 // is still responsible of the memory of the points (here the
337 // vector<double>). No memory is copied, Delaunay just keeps
338 // a pointer.
339
2/4
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 9 times.
✗ Branch 5 not taken.
 9 delaunay->set_vertices(nb_points, points.data());
340
1/2
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
 9 time = Wdel.elapsed_time();
341 9 }
342
343
3/6
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 9 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 9 times.
✗ Branch 8 not taken.
 9 Logger::out("Delaunay") << delaunay->nb_cells() << " tetrahedra"
344 << std::endl;
345
346
4/10
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 9 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 9 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 9 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
 9 Logger::out("Delaunay") << double(delaunay->nb_cells()) / time
347 << " tetrahedra / second"
348 << std::endl;
349
350
1/2
✓ Branch 0 taken 9 times.
✗ Branch 1 not taken.
 9 if(output) {
351
1/2
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
 9 save_Delaunay(delaunay, output_filename, convex_hull_only);
352 }
353 9 }
354 ✗ catch(const std::exception& e) {
355 ✗ std::cerr << "Received an exception: " << e.what() << std::endl;
356 return 1;
357 ✗ }
358
359
1/2
✓ Branch 2 taken 9 times.
✗ Branch 3 not taken.
 9 Logger::out("") << "Everything OK, Returning status 0" << std::endl;
360 9 return 0;
361 }
362