GCC Code Coverage Report


Directory: ./
File: examples/geogram/compute_RVD/main.cpp
Date: 2026-09-07 02:36:43
Exec Total Coverage
Lines: 123 311 39.5%
Functions: 2 13 15.4%
Branches: 194 766 25.3%

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/basic/process.h>
47 #include <geogram/mesh/mesh.h>
48 #include <geogram/mesh/mesh_geometry.h>
49 #include <geogram/mesh/mesh_topology.h>
50 #include <geogram/mesh/mesh_io.h>
51 #include <geogram/mesh/mesh_repair.h>
52 #include <geogram/mesh/mesh_fill_holes.h>
53 #include <geogram/mesh/mesh_preprocessing.h>
54 #include <geogram/mesh/mesh_degree3_vertices.h>
55 #include <geogram/mesh/mesh_tetrahedralize.h>
56 #include <geogram/delaunay/delaunay.h>
57 #include <geogram/delaunay/periodic_delaunay_3d.h>
58 #include <geogram/voronoi/RVD.h>
59 #include <geogram/voronoi/RVD_callback.h>
60 #include <geogram/voronoi/RVD_mesh_builder.h>
61 #include <geogram/voronoi/convex_cell.h>
62 #include <geogram/numerics/predicates.h>
63
64 namespace {
65
66 using namespace GEO;
67
68 /**
69 * \brief Removes zero area facets in a mesh
70 * \param[in] M the input mesh
71 */
72 36 void check_for_zero_area_facets(Mesh& M) {
73 36 vector<index_t> remove_f;
74 36 vec3 q1(0, 0, 0);
75 36 vec3 q2(0, 0, 1);
76 36 vec3 q3(0, 1, 0);
77 36 vec3 q4(1, 0, 0);
78
2/2
✓ Branch 1 taken 60502 times.
✓ Branch 2 taken 36 times.
60538 for(index_t f = 0; f < M.facets.nb(); ++f) {
79
1/2
✓ Branch 1 taken 60502 times.
✗ Branch 2 not taken.
60502 index_t c = M.facets.corners_begin(f);
80
1/2
✓ Branch 1 taken 60502 times.
✗ Branch 2 not taken.
60502 index_t v1 = M.facet_corners.vertex(c);
81
1/2
✓ Branch 1 taken 60502 times.
✗ Branch 2 not taken.
60502 index_t v2 = M.facet_corners.vertex(c + 1);
82
1/2
✓ Branch 1 taken 60502 times.
✗ Branch 2 not taken.
60502 index_t v3 = M.facet_corners.vertex(c + 2);
83
1/2
✓ Branch 1 taken 60502 times.
✗ Branch 2 not taken.
60502 const vec3& p1 = M.vertices.point(v1);
84
1/2
✓ Branch 1 taken 60502 times.
✗ Branch 2 not taken.
60502 const vec3& p2 = M.vertices.point(v2);
85
1/2
✓ Branch 1 taken 60502 times.
✗ Branch 2 not taken.
60502 const vec3& p3 = M.vertices.point(v3);
86
87 // Colinearity is tested by using four coplanarity
88 // tests with points q1,q2,q3,q4 that are
89 // not coplanar.
90 60502 if(
91
1/2
✓ Branch 1 taken 60502 times.
✗ Branch 2 not taken.
60502 PCK::orient_3d(p1, p2, p3, q1) == 0 &&
92
3/4
✓ Branch 1 taken 22272 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 128 times.
✓ Branch 4 taken 22144 times.
22272 PCK::orient_3d(p1, p2, p3, q2) == 0 &&
93
5/8
✓ Branch 0 taken 22272 times.
✓ Branch 1 taken 38230 times.
✓ Branch 3 taken 128 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✓ Branch 6 taken 128 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 60502 times.
82774 PCK::orient_3d(p1, p2, p3, q3) == 0 &&
94 PCK::orient_3d(p1, p2, p3, q4) == 0
95 ) {
96 Logger::warn("Validate") << "Found a zero-area facet"
97 << std::endl;
98 remove_f.resize(M.facets.nb(), 0);
99 remove_f[f] = 1;
100 }
101 }
102
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 36 times.
36 if(remove_f.size() != 0) {
103 Logger::warn("Validate") << "Removing zero-area facet(s)"
104 << std::endl;
105 M.facets.delete_elements(remove_f);
106 }
107 36 }
108
109 /**
110 * \brief The callback called for each RVD polyhedron. Constructs a
111 * mesh with the boundary of all cells.
112 * \details Its member functions are called for each RVD polyhedron,
113 * i.e. the intersections between the volumetric mesh tetrahedra and
114 * the Voronoi cells. Based on set_simplify_xxx(), a smaller number of
115 * polyhedra can be generated.
116 */
117 class SaveRVDCells : public RVDPolyhedronCallback {
118 public:
119
120 /**
121 * \brief SaveRVDCells constructor.
122 * \param[out] output_mesh a reference to the generated mesh
123 */
124 SaveRVDCells(Mesh& output_mesh) :
125 output_mesh_(output_mesh),
126 facet_seed_attr_(output_mesh.facets.attributes(), "region"),
127 current_seed_(NO_INDEX)
128 {
129 my_vertex_map_ = nullptr;
130
131 // If set, then only one polyhedron per (connected
132 // component of) restricted Voronoi cell is generated.
133 set_simplify_internal_tet_facets(
134 CmdLine::get_arg_bool("RVD_cells:simplify_tets")
135 );
136
137 // If set, then only one polygon per Voronoi facet is generated.
138 set_simplify_voronoi_facets(
139 CmdLine::get_arg_bool("RVD_cells:simplify_voronoi")
140 );
141
142 // If set, then the intersection between a Voronoi cell
143 // and the boundary surface is replaced with a single
144 // polygon whenever possible (i.e. when its topology is a
145 // disk and when it has at least 3 corners).
146 set_simplify_boundary_facets(
147 CmdLine::get_arg_bool("RVD_cells:simplify_boundary"),
148 CmdLine::get_arg_double(
149 "RVD_cells:simplify_boundary_angle_threshold"
150 )
151 );
152
153 // If set, then the intersections are available as Mesh
154 // objects through the function
155 // process_polyhedron_mesh(). Note that this is implied by
156 // simplify_voronoi_facets or simplify_boundary.
157 if(CmdLine::get_arg_double("RVD_cells:shrink") != 0.0) {
158 set_use_mesh(true);
159 }
160 }
161
162 ~SaveRVDCells() override {
163 delete my_vertex_map_;
164 my_vertex_map_ = nullptr;
165 }
166
167 /**
168 * \brief Called at the beginning of RVD traversal.
169 */
170 void begin() override {
171 RVDPolyhedronCallback::begin();
172 output_mesh_.clear();
173 output_mesh_.vertices.set_dimension(3);
174 }
175
176 /**
177 * \brief Called at the end of RVD traversal.
178 */
179 void end() override {
180 RVDPolyhedronCallback::end();
181 output_mesh_.facets.connect();
182 }
183
184 /**
185 * \brief Called at the beginning of each RVD polyhedron.
186 * \param[in] seed , tetrahedron the (seed,tetrahedron) pair that
187 * defines the RVD polyhedron, as the intersection between the Voronoi
188 * cell of the seed and the tetrahedron.
189 */
190 void begin_polyhedron(index_t seed, index_t tetrahedron) override {
191 geo_argused(tetrahedron);
192 current_seed_ = seed;
193
194 // The RVDVertexMap is used to map the symbolic
195 // representation of vertices to indices. Here we reset
196 // indexing for each new cell, so that vertices shared by
197 // the faces of two different cells will be duplicated. We
198 // do that because we construct the boundary of the cells
199 // in a surfacic mesh (for visualization purposes). Client
200 // code that has a data structure for polyhedral
201 // volumetric mesh will not want to reset indexing (and
202 // will comment-out the following three lines). It will
203 // also construct the RVDVertexMap in the constructor.
204
205 delete my_vertex_map_;
206 my_vertex_map_ = new RVDVertexMap;
207 my_vertex_map_->set_first_vertex_index(output_mesh_.vertices.nb());
208 }
209
210 /**
211 * \brief Called at the beginning of each RVD polyhedron.
212 * \param[in] facet_seed if the facet is on a Voronoi bisector,
213 * the index of the Voronoi seed on the other side of the bisector,
214 * else index_t(-1)
215 * \param[in] facet_tet if the facet is on a tethedral facet, then
216 * the index of the tetrahedron on the other side, else index_t(-1)
217 */
218 void begin_facet(index_t facet_seed, index_t facet_tet) override {
219 geo_argused(facet_seed);
220 geo_argused(facet_tet);
221 current_facet_.resize(0);
222 }
223
224 void vertex(
225 const double* geometry, const GEOGen::SymbolicVertex& symb
226 ) override {
227 // Find the index of the vertex associated with its
228 // symbolic representation.
229 index_t vid = my_vertex_map_->find_or_create_vertex(seed(), symb);
230
231 // If the vertex does not exist in the mesh, create it.
232 if(vid >= output_mesh_.vertices.nb()) {
233 output_mesh_.vertices.create_vertex(geometry);
234 }
235
236 // Memorize the current facet.
237 current_facet_.push_back(vid);
238 }
239
240 void end_facet() override {
241 // Create the facet from the memorized indices.
242 index_t f = output_mesh_.facets.nb();
243 output_mesh_.facets.create_polygon(current_facet_.size());
244 for(index_t i=0; i<current_facet_.size(); ++i) {
245 output_mesh_.facets.set_vertex(f,i,current_facet_[i]);
246 }
247 facet_seed_attr_[f] = current_seed_;
248 }
249
250 void end_polyhedron() override {
251 // Nothing to do.
252 }
253
254 void process_polyhedron_mesh() override {
255 // This function is called for each cell if
256 // set_use_mesh(true) was called. It is the case if
257 // simplify_voronoi_facets(true) or
258 // simplify_boundary_facets(true) was called. Note1: most
259 // users will not need to overload this function (advanded
260 // use only). Note2: mesh_ is managed internally by
261 // RVDPolyhedronCallback class, as an intermediary
262 // representation to store the cell before calling the
263 // callbacks. It is distinct from the output_mesh_
264 // constructed by the callbacks.
265
266 // The current cell represented by a Mesh can be
267 // filtered/modified/post-processed (member variable
268 // mesh_) here, before calling base class's
269 // implementation. As an example, we shrink the
270 // cells. More drastic modifications/ transformations of
271 // the mesh can be done (see base class's implementation
272 // in geogram/voronoi/RVD_polyhedron_callback.cpp).
273
274 double shrink = CmdLine::get_arg_double("RVD_cells:shrink");
275 if(shrink != 0.0 && mesh_.vertices.nb() != 0) {
276 vec3 center(0.0, 0.0, 0.0);
277 for(index_t v=0; v<mesh_.vertices.nb(); ++v) {
278 center += vec3(mesh_.vertices.point_ptr(v));
279 }
280 center = (1.0 / double(mesh_.vertices.nb())) * center;
281 for(index_t v=0; v<mesh_.vertices.nb(); ++v) {
282 vec3 p(mesh_.vertices.point_ptr(v));
283 p = shrink * center + (1.0 - shrink) * p;
284 mesh_.vertices.point_ptr(v)[0] = p.x;
285 mesh_.vertices.point_ptr(v)[1] = p.y;
286 mesh_.vertices.point_ptr(v)[2] = p.z;
287 }
288 }
289
290 // The default implementation simplifies Voronoi facets
291 // and boundary mesh facets based on the boolean flags
292 // defined by set_simplify_xxx(). Then it calls the callbacks
293 // for each mesh facet.
294 RVDPolyhedronCallback::process_polyhedron_mesh();
295
296 }
297
298 private:
299 vector<index_t> current_facet_;
300 Mesh& output_mesh_;
301 RVDVertexMap* my_vertex_map_;
302 Attribute<index_t> facet_seed_attr_;
303 index_t current_seed_;
304 };
305
306 void compute_RVD_cells(RestrictedVoronoiDiagram* RVD, Mesh& RVD_mesh) {
307 SaveRVDCells callback(RVD_mesh);
308 RVD->for_each_polyhedron(callback);
309 }
310
311 }
312
313 37 int main(int argc, char** argv) {
314 using namespace GEO;
315
316 37 GEO::initialize(GEO::GEOGRAM_INSTALL_ALL);
317
318 try {
319
320
2/4
✓ Branch 1 taken 37 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 37 times.
✗ Branch 5 not taken.
37 Stopwatch Wtot("Total time");
321
322 37 std::vector<std::string> filenames;
323
324
2/4
✓ Branch 1 taken 37 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 37 times.
✗ Branch 5 not taken.
74 CmdLine::import_arg_group("standard");
325
2/4
✓ Branch 1 taken 37 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 37 times.
✗ Branch 5 not taken.
74 CmdLine::import_arg_group("algo");
326
3/6
✓ Branch 1 taken 37 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 37 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 37 times.
✗ Branch 8 not taken.
148 CmdLine::declare_arg("volumetric", false, "volumetric or surfacic RVD");
327
1/2
✓ Branch 1 taken 37 times.
✗ Branch 2 not taken.
37 CmdLine::declare_arg(
328
2/4
✓ Branch 1 taken 37 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 37 times.
✗ Branch 5 not taken.
185 "cell_borders", false, "generate only cell borders"
329 );
330
1/2
✓ Branch 1 taken 37 times.
✗ Branch 2 not taken.
37 CmdLine::declare_arg(
331
1/2
✓ Branch 1 taken 37 times.
✗ Branch 2 not taken.
74 "integration_smplx", false,
332
1/2
✓ Branch 1 taken 37 times.
✗ Branch 2 not taken.
111 "in volumetric mode, generate integration simplices"
333 );
334
3/6
✓ Branch 1 taken 37 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 37 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 37 times.
✗ Branch 8 not taken.
148 CmdLine::declare_arg("RDT", false, "save RDT");
335
3/6
✓ Branch 1 taken 37 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 37 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 37 times.
✗ Branch 8 not taken.
148 CmdLine::declare_arg("RVD", true, "save RVD");
336
1/2
✓ Branch 1 taken 37 times.
✗ Branch 2 not taken.
37 CmdLine::declare_arg(
337
1/2
✓ Branch 1 taken 37 times.
✗ Branch 2 not taken.
74 "RVD_cells", false,
338
1/2
✓ Branch 1 taken 37 times.
✗ Branch 2 not taken.
111 "use new API for computing RVD cells (implies volumetric)"
339 );
340
1/2
✓ Branch 1 taken 37 times.
✗ Branch 2 not taken.
37 CmdLine::declare_arg_group(
341
2/4
✓ Branch 1 taken 37 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 37 times.
✗ Branch 5 not taken.
185 "RVD_cells", "RVD cells simplification flags"
342 );
343
1/2
✓ Branch 1 taken 37 times.
✗ Branch 2 not taken.
37 CmdLine::declare_arg(
344
2/4
✓ Branch 1 taken 37 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 37 times.
✗ Branch 5 not taken.
185 "RVD_cells:simplify_tets", true, "Simplify tets intersections"
345 );
346
1/2
✓ Branch 1 taken 37 times.
✗ Branch 2 not taken.
37 CmdLine::declare_arg(
347
2/4
✓ Branch 1 taken 37 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 37 times.
✗ Branch 5 not taken.
185 "RVD_cells:simplify_voronoi", true, "Simplify Voronoi facets"
348 );
349
1/2
✓ Branch 1 taken 37 times.
✗ Branch 2 not taken.
37 CmdLine::declare_arg(
350
2/4
✓ Branch 1 taken 37 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 37 times.
✗ Branch 5 not taken.
185 "RVD_cells:simplify_boundary", false, "Simplify boundary facets"
351 );
352
1/2
✓ Branch 1 taken 37 times.
✗ Branch 2 not taken.
37 CmdLine::declare_arg(
353
1/2
✓ Branch 1 taken 37 times.
✗ Branch 2 not taken.
74 "RVD_cells:simplify_boundary_angle_threshold", 45.0,
354
1/2
✓ Branch 1 taken 37 times.
✗ Branch 2 not taken.
111 "Angle below which boundary facets are simplified."
355 "Only applies if simplify_boundary is `true`."
356 );
357
1/2
✓ Branch 1 taken 37 times.
✗ Branch 2 not taken.
37 CmdLine::declare_arg(
358
2/4
✓ Branch 1 taken 37 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 37 times.
✗ Branch 5 not taken.
185 "RVD_cells:shrink", 0.0, "Shrink factor for computed cells"
359 );
360
361
362
1/2
✓ Branch 1 taken 37 times.
✗ Branch 2 not taken.
37 CmdLine::declare_arg_percent(
363
1/2
✓ Branch 1 taken 37 times.
✗ Branch 2 not taken.
74 "epsilon",0.001,
364
1/2
✓ Branch 1 taken 37 times.
✗ Branch 2 not taken.
111 "Tolerance for merging vertices relative to bbox diagonal"
365 );
366
3/6
✓ Branch 1 taken 37 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 37 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 37 times.
✗ Branch 8 not taken.
148 CmdLine::declare_arg("constrained", false, "constrained Delaunay");
367
1/2
✓ Branch 1 taken 37 times.
✗ Branch 2 not taken.
37 CmdLine::declare_arg(
368
1/2
✓ Branch 1 taken 37 times.
✗ Branch 2 not taken.
74 "prefer_seeds", false,
369
1/2
✓ Branch 1 taken 37 times.
✗ Branch 2 not taken.
111 "in constrained mode, use seeds whenever possible"
370 );
371
372 37 if(
373
1/2
✓ Branch 1 taken 37 times.
✗ Branch 2 not taken.
37 !CmdLine::parse(
374
3/4
✓ Branch 1 taken 37 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✓ Branch 5 taken 36 times.
111 argc, argv, filenames, "meshfile <pointsfile> <outputfile>"
375 )
376 ) {
377 1 return 1;
378 }
379
380
3/6
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 36 times.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✓ Branch 8 taken 36 times.
72 if(CmdLine::get_arg_bool("RVD_cells")) {
381 CmdLine::set_arg("volumetric",true);
382 }
383
384
1/2
✓ Branch 2 taken 36 times.
✗ Branch 3 not taken.
36 std::string mesh_filename = filenames[0];
385
1/2
✓ Branch 2 taken 36 times.
✗ Branch 3 not taken.
36 std::string points_filename = filenames[0];
386
2/2
✓ Branch 1 taken 19 times.
✓ Branch 2 taken 17 times.
36 if(filenames.size() >= 2) {
387
1/2
✓ Branch 2 taken 19 times.
✗ Branch 3 not taken.
19 points_filename = filenames[1];
388 }
389
390
2/4
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 36 times.
✗ Branch 5 not taken.
72 bool volumetric = CmdLine::get_arg_bool("volumetric");
391
2/4
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 36 times.
✗ Branch 5 not taken.
72 bool cell_borders = CmdLine::get_arg_bool("cell_borders");
392
2/4
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 36 times.
✗ Branch 5 not taken.
36 bool integ_smplx = CmdLine::get_arg_bool("integration_smplx");
393 36 std::string output_filename;
394
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 36 times.
36 if(filenames.size() >= 3) {
395 output_filename = filenames[2];
396 } else {
397
7/18
✓ Branch 0 taken 36 times.
✗ Branch 1 not taken.
✓ Branch 3 taken 36 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 36 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 36 times.
✓ Branch 10 taken 36 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 36 times.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✓ Branch 16 taken 36 times.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
108 if(volumetric || CmdLine::get_arg_bool("constrained")) {
398 if(CmdLine::get_arg_bool("RVD_cells")) {
399 output_filename = "out.obj";
400 } else {
401 output_filename = "out.meshb";
402 }
403 } else {
404
1/2
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
36 output_filename = "out.eobj";
405 }
406 }
407
408
5/10
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 36 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 36 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 36 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 36 times.
✗ Branch 14 not taken.
72 Logger::out("I/O") << "Output = " << output_filename << std::endl;
409
410
2/4
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 36 times.
✗ Branch 5 not taken.
36 Logger::div("Loading data");
411
412
2/4
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 36 times.
✗ Branch 5 not taken.
36 Mesh M_in, points_in;
413
1/2
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
36 Mesh M_out;
414
1/2
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
36 bool cube = (mesh_filename == "cube");
415
416
1/2
✓ Branch 0 taken 36 times.
✗ Branch 1 not taken.
36 if(!cube) {
417
1/2
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
36 MeshIOFlags flags;
418
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 36 times.
36 if(volumetric) {
419 flags.set_element(MESH_CELLS);
420 }
421
2/4
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 36 times.
36 if(!mesh_load(mesh_filename, M_in, flags)) {
422 return 1;
423 }
424
1/2
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
36 }
425
426
1/2
✓ Branch 0 taken 36 times.
✗ Branch 1 not taken.
36 if(!volumetric) {
427
1/2
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
36 mesh_repair(M_in);
428
1/2
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
36 check_for_zero_area_facets(M_in);
429 } else {
430 if(M_in.cells.nb() == 0) {
431 Logger::out("RVD")
432 << "Mesh does not have tetrahedra, tetrahedralizing"
433 << std::endl;
434 mesh_tetrahedralize(M_in);
435 }
436 }
437
438
3/6
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 36 times.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✓ Branch 8 taken 36 times.
36 if(!mesh_load(points_filename, points_in)) {
439 return 1;
440 }
441
442
2/4
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 36 times.
✗ Branch 5 not taken.
72 double epsilon = CmdLine::get_arg_percent(
443
1/2
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
72 "epsilon",bbox_diagonal(points_in)
444 );
445
1/2
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
36 points_in.facets.clear();
446
1/2
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
36 mesh_repair(points_in, MESH_REPAIR_COLOCATE, epsilon);
447
448
2/8
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 36 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
36 geo_assert(points_in.vertices.dimension() == 3);
449
450
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 36 times.
36 if(cube) {
451 double shrink = CmdLine::get_arg_double("RVD_cells:shrink");
452 SmartPointer<PeriodicDelaunay3d> delaunay =
453 new PeriodicDelaunay3d(false);
454 delaunay->set_keeps_infinite(true);
455 delaunay->set_vertices(
456 points_in.vertices.nb(), points_in.vertices.point_ptr(0)
457 );
458 delaunay->compute();
459 ConvexCell C;
460 PeriodicDelaunay3d::IncidentTetrahedra W;
461 index_t cur_v_index = 1;
462 if(FileSystem::extension(output_filename) != "obj") {
463 Logger::err("RVD")
464 << "cube mode only available in .obj file format"
465 << std::endl;
466 exit(-1);
467 }
468
469
470 std::ofstream out(output_filename);
471 for(index_t v=0; v<delaunay->nb_vertices(); ++v) {
472 delaunay->copy_Laguerre_cell_from_Delaunay(v, C, W);
473 C.clip_by_plane(vec4( 1.0, 0.0, 0.0, 0.0));
474 C.clip_by_plane(vec4(-1.0, 0.0, 0.0, 1.0));
475 C.clip_by_plane(vec4( 0.0, 1.0, 0.0, 0.0));
476 C.clip_by_plane(vec4( 0.0,-1.0, 0.0, 1.0));
477 C.clip_by_plane(vec4( 0.0, 0.0, 1.0, 0.0));
478 C.clip_by_plane(vec4( 0.0, 0.0,-1.0, 1.0));
479 out << "# CELL " << v << std::endl;
480 cur_v_index += C.save(out, cur_v_index, shrink);
481 }
482 exit(0);
483
3/6
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 36 times.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✓ Branch 8 taken 36 times.
72 } else if(CmdLine::get_arg_bool("constrained")) {
484
485 Mesh surface;
486 vector<double> inner_points;
487
488
489 {
490 Logger::div("Computing the surface");
491 Delaunay_var delaunay = Delaunay::create(3);
492 RestrictedVoronoiDiagram_var RVD =
493 RestrictedVoronoiDiagram::create(delaunay,&M_in);
494 delaunay->set_vertices(
495 points_in.vertices.nb(), points_in.vertices.point_ptr(0)
496 );
497
498
499 RestrictedVoronoiDiagram::RDTMode mode =
500 RestrictedVoronoiDiagram::RDTMode(
501 RestrictedVoronoiDiagram::RDT_MULTINERVE |
502 RestrictedVoronoiDiagram::RDT_RVC_CENTROIDS
503 );
504
505 if(CmdLine::get_arg_bool("prefer_seeds")) {
506 mode = RestrictedVoronoiDiagram::RDTMode(
507 mode | RestrictedVoronoiDiagram::RDT_PREFER_SEEDS
508 );
509 }
510
511
512 RVD->compute_RDT(surface, mode);
513
514 mesh_repair(surface);
515 remove_small_connected_components(surface,0.0,100);
516 fill_holes(surface, 1e30);
517 double radius = bbox_diagonal(surface);
518 remove_degree3_vertices(surface, 0.01*radius);
519 mesh_save(surface,"surface.meshb");
520
521 vector<double> m(points_in.vertices.nb());
522 vector<double> mg(points_in.vertices.nb()*3);
523 RVD->compute_centroids_on_surface(&mg[0], &m[0]);
524 for(index_t v=0; v<points_in.vertices.nb(); ++v) {
525 if(m[v] == 0.0) {
526 inner_points.push_back(
527 points_in.vertices.point_ptr(v)[0]
528 );
529 inner_points.push_back(
530 points_in.vertices.point_ptr(v)[1]
531 );
532 inner_points.push_back(
533 points_in.vertices.point_ptr(v)[2]
534 );
535 }
536 }
537 }
538
539 Logger::div("Calling tetgen");
540 Delaunay_var delaunay = Delaunay::create(3,"tetgen");
541 delaunay->set_constraints(&surface);
542 delaunay->set_vertices(inner_points.size()/3, &inner_points[0]);
543
544 vector<double> pts(delaunay->nb_vertices() * 3);
545 vector<index_t> tet2v(delaunay->nb_cells() * 4);
546 for(index_t v = 0; v < delaunay->nb_vertices(); ++v) {
547 pts[3 * v] = delaunay->vertex_ptr(v)[0];
548 pts[3 * v + 1] = delaunay->vertex_ptr(v)[1];
549 pts[3 * v + 2] = delaunay->vertex_ptr(v)[2];
550 }
551 for(index_t t = 0; t < delaunay->nb_cells(); ++t) {
552 tet2v[4 * t] = index_t(delaunay->cell_vertex(t, 0));
553 tet2v[4 * t + 1] = index_t(delaunay->cell_vertex(t, 1));
554 tet2v[4 * t + 2] = index_t(delaunay->cell_vertex(t, 2));
555 tet2v[4 * t + 3] = index_t(delaunay->cell_vertex(t, 3));
556 }
557 M_out.cells.assign_tet_mesh(3, pts, tet2v, true);
558 M_out.show_stats();
559
560 Logger::div("Saving the result");
561 MeshIOFlags flags;
562 flags.set_element(MESH_CELLS);
563 mesh_save(M_out, output_filename, flags);
564 } else {
565
1/2
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
36 Delaunay_var delaunay = Delaunay::create(3);
566 RestrictedVoronoiDiagram_var RVD = RestrictedVoronoiDiagram::create(
567 delaunay, &M_in
568
1/2
✓ Branch 2 taken 36 times.
✗ Branch 3 not taken.
36 );
569 {
570
2/4
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 36 times.
✗ Branch 5 not taken.
36 Stopwatch W("Delaunay");
571
2/4
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 36 times.
✗ Branch 5 not taken.
72 delaunay->set_vertices(
572
1/2
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
36 points_in.vertices.nb(), points_in.vertices.point_ptr(0)
573 );
574 36 }
575
576
2/4
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 36 times.
✗ Branch 5 not taken.
36 RVD->set_volumetric(volumetric);
577
578
3/6
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 36 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 36 times.
✗ Branch 8 not taken.
72 if(CmdLine::get_arg_bool("RVD")) {
579
2/4
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 36 times.
✗ Branch 5 not taken.
72 Logger::div("Restricted Voronoi Diagram");
580 {
581
2/4
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 36 times.
✗ Branch 5 not taken.
72 Stopwatch W("RVD");
582
3/6
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 36 times.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✓ Branch 8 taken 36 times.
72 if(CmdLine::get_arg_bool("RVD_cells")) {
583 compute_RVD_cells(RVD, M_out);
584 } else {
585
2/4
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 36 times.
✗ Branch 5 not taken.
36 RVD->compute_RVD(M_out, 0, cell_borders, integ_smplx);
586
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
36 if(integ_smplx && volumetric) {
587 M_out.cells.connect();
588 M_out.cells.compute_borders();
589 }
590 }
591 36 }
592
2/4
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 36 times.
✗ Branch 5 not taken.
36 Logger::div("Result");
593
594
1/2
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
36 MeshIOFlags flags;
595 36 flags.set_attribute(MESH_FACET_REGION);
596 36 flags.set_attribute(MESH_CELL_REGION);
597 36 flags.set_element(MESH_CELLS);
598
1/2
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
36 mesh_save(M_out, output_filename, flags);
599 36 }
600
601
3/6
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 36 times.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✓ Branch 8 taken 36 times.
72 if(CmdLine::get_arg_bool("RDT")) {
602 Logger::out("RDT") << "Computing RDT..." << std::endl;
603 Mesh RDT ;
604 RVD->compute_RDT(RDT);
605 MeshIOFlags flags;
606 if(volumetric) {
607 flags.set_elements(MESH_CELLS);
608 }
609 mesh_save(RDT, "RDT.meshb", flags);
610 }
611
612
6/8
✓ Branch 0 taken 36 times.
✗ Branch 1 not taken.
✓ Branch 3 taken 36 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
✓ Branch 6 taken 35 times.
✓ Branch 7 taken 1 times.
✓ Branch 8 taken 35 times.
36 if(!volumetric && !meshes_have_same_topology(M_in, M_out, true)) {
613
4/8
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
1 Logger::out("") << "Returning error code (2)" << std::endl;
614 1 return 2;
615 }
616
4/4
✓ Branch 1 taken 35 times.
✓ Branch 2 taken 1 times.
✓ Branch 4 taken 35 times.
✓ Branch 5 taken 1 times.
37 }
617
16/16
✓ Branch 1 taken 35 times.
✓ Branch 2 taken 1 times.
✓ Branch 4 taken 35 times.
✓ Branch 5 taken 1 times.
✓ Branch 7 taken 35 times.
✓ Branch 8 taken 1 times.
✓ Branch 10 taken 35 times.
✓ Branch 11 taken 1 times.
✓ Branch 13 taken 35 times.
✓ Branch 14 taken 1 times.
✓ Branch 16 taken 35 times.
✓ Branch 17 taken 1 times.
✓ Branch 19 taken 35 times.
✓ Branch 20 taken 2 times.
✓ Branch 22 taken 35 times.
✓ Branch 23 taken 2 times.
45 }
618 catch(const std::exception& e) {
619 std::cerr << "Received an exception: " << e.what() << std::endl;
620 return 1;
621 }
622
623
4/8
✓ Branch 1 taken 35 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 35 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 35 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 35 times.
✗ Branch 11 not taken.
35 Logger::out("") << "Everything OK, Returning status 0" << std::endl;
624 35 return 0;
625 }
626