GCC Code Coverage Report

Directory: ./
File: tests/test_RVC/main.cpp
Date: 2026-09-07 02:36:43
Exec Total Coverage
Lines: 0 130 0.0%
Functions: 0 6 0.0%
Branches: 0 276 0.0%
Line Branch Exec Source
1 /*
2 * Copyright (c) 2000-2022 Inria
3 * All rights reserved.
4 *
5 * Redistribution and use in source and binary forms, with or without
6 * modification, are permitted provided that the following conditions are met:
7 *
8 * * Redistributions of source code must retain the above copyright notice,
9 * this list of conditions and the following disclaimer.
10 * * Redistributions in binary form must reproduce the above copyright notice,
11 * this list of conditions and the following disclaimer in the documentation
12 * and/or other materials provided with the distribution.
13 * * Neither the name of the ALICE Project-Team nor the names of its
14 * contributors may be used to endorse or promote products derived from this
15 * software without specific prior written permission.
16 *
17 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
18 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
19 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
20 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
21 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
22 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
23 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
24 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
25 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
26 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
27 * POSSIBILITY OF SUCH DAMAGE.
28 *
29 * Contact: Bruno Levy
30 *
31 * https://www.inria.fr/fr/bruno-levy
32 *
33 * Inria,
34 * Domaine de Voluceau,
35 * 78150 Le Chesnay - Rocquencourt
36 * FRANCE
37 *
38 */
39
40 #include <geogram/mesh/mesh.h>
41 #include <geogram/mesh/mesh_io.h>
42 #include <geogram/mesh/mesh_geometry.h>
43 #include <geogram/mesh/mesh_repair.h>
44 #include <geogram/delaunay/delaunay.h>
45 #include <geogram/voronoi/RVD.h>
46 #include <geogram/numerics/predicates.h>
47 #include <geogram/basic/logger.h>
48 #include <geogram/basic/command_line.h>
49 #include <geogram/basic/command_line_args.h>
50 #include <geogram/basic/file_system.h>
51 #include <geogram/basic/progress.h>
52 #include <stdarg.h>
53
54
55 namespace {
56
57 using namespace GEO;
58
59 /**
60 * \brief Creates a surfacic mesh from a cube.
61 * \param[out] M the resulting mesh
62 */
63 ✗ void initialize_mesh_with_box(Mesh& M) {
64 ✗ M.clear();
65 ✗ M.vertices.set_dimension(3);
66
67 ✗ const double d = 1.0;
68
69 ✗ M.vertices.create_vertex(vec3(-d, -d, -d).data());
70 ✗ M.vertices.create_vertex(vec3(-d, -d, d).data());
71 ✗ M.vertices.create_vertex(vec3(-d, d, -d).data());
72 ✗ M.vertices.create_vertex(vec3(-d, d, d).data());
73 ✗ M.vertices.create_vertex(vec3(d, -d, -d).data());
74 ✗ M.vertices.create_vertex(vec3(d, -d, d).data());
75 ✗ M.vertices.create_vertex(vec3(d, d, -d).data());
76 ✗ M.vertices.create_vertex(vec3(d, d, d).data());
77
78 ✗ M.facets.create_quad(7,6,2,3);
79 ✗ M.facets.create_quad(1,3,2,0);
80 ✗ M.facets.create_quad(5,7,3,1);
81 ✗ M.facets.create_quad(4,6,7,5);
82 ✗ M.facets.create_quad(4,5,1,0);
83 ✗ M.facets.create_quad(6,4,0,2);
84
85 ✗ M.facets.connect();
86 ✗ }
87
88 ✗ void center_scale_mesh(Mesh& M, vec3 center, double radius) {
89 double xyz_min[3];
90 double xyz_max[3];
91 ✗ get_bbox(M, xyz_min, xyz_max);
92 ✗ double scale = 2.0*radius / Geom::distance(xyz_min, xyz_max, 3);
93 ✗ vec3 g = 0.5*(vec3(xyz_min) + vec3(xyz_max));
94 ✗ for(index_t i=0; i<M.vertices.nb(); ++i) {
95 ✗ M.vertices.point(i) -= g;
96 ✗ M.vertices.point(i) *= scale;
97 ✗ M.vertices.point(i) += center;
98 }
99 ✗ }
100
101 /**
102 * \brief Shrinks a mesh.
103 * \param[in,out] M the mesh to be shrunk
104 * \param[in] factor the shrinking factor (1.0 means
105 * no shrinking, 0.5 means average shrinking).
106 */
107 ✗ void shrink_mesh(Mesh& M, double factor) {
108 double xyz_min[3];
109 double xyz_max[3];
110 ✗ get_bbox(M, xyz_min, xyz_max);
111 ✗ vec3 g = 0.5*(vec3(xyz_min) + vec3(xyz_max));
112 ✗ for(index_t i=0; i<M.vertices.nb(); ++i) {
113 ✗ M.vertices.point(i) -= g;
114 ✗ M.vertices.point(i) *= factor;
115 ✗ M.vertices.point(i) += g;
116 }
117 ✗ }
118
119
120 /**
121 * \brief Tests whether the facets of a mesh are exactly planar.
122 * \retval true if all the facets are exactly planar
123 * \retval false otherwise
124 */
125 ✗ bool mesh_facets_are_planar(const Mesh& M) {
126 ✗ for(index_t f=0; f<M.facets.nb(); ++f) {
127 ✗ for(index_t c=M.facets.corners_begin(f); c+3<M.facets.corners_end(f); ++c) {
128 ✗ index_t v1 = M.facet_corners.vertex(c);
129 ✗ index_t v2 = M.facet_corners.vertex(c+1);
130 ✗ index_t v3 = M.facet_corners.vertex(c+2);
131 ✗ index_t v4 = M.facet_corners.vertex(c+3);
132 ✗ if(
133 ✗ PCK::orient_3d(
134 M.vertices.point_ptr(v1),
135 M.vertices.point_ptr(v2),
136 M.vertices.point_ptr(v3),
137 M.vertices.point_ptr(v4)
138 ✗ ) != ZERO
139 ) {
140 ✗ return false;
141 }
142 }
143 }
144 ✗ return true;
145 }
146
147
148 /**
149 * \brief Tests whether all the vertices of a mesh are of degree 3.
150 * \retval true if all the vertices are of degree 3
151 * \retval false otherwise
152 */
153 ✗ bool mesh_vertices_are_degree_3(const Mesh& M) {
154 ✗ vector<int> degree(M.vertices.nb(),0);
155 ✗ for(index_t f=0; f<M.facets.nb(); ++f) {
156 ✗ for(index_t c=M.facets.corners_begin(f); c<M.facets.corners_end(f); ++c) {
157 ✗ ++degree[M.facet_corners.vertex(c)];
158 }
159 }
160 ✗ for(index_t v=0; v<degree.size(); ++v) {
161 ✗ if(degree[v] != 3) {
162 ✗ return false;
163 }
164 }
165 ✗ return true;
166 ✗ }
167 }
168
169
170 ✗ int main(int argc, char** argv) {
171
172 ✗ GEO::initialize(GEO::GEOGRAM_INSTALL_ALL);
173 ✗ GEO::Logger::instance()->set_quiet(false);
174 ✗ GEO::CmdLine::import_arg_group("standard");
175 ✗ GEO::CmdLine::import_arg_group("algo");
176 ✗ GEO::CmdLine::declare_arg_percent(
177 ✗ "size", 10.0, "elements size, in bbox diagonal percent"
178 );
179 ✗ GEO::CmdLine::declare_arg("shrink", 0.9, "cells shrink");
180 ✗ GEO::CmdLine::declare_arg(
181 ✗ "border_only", false, "output only RVC facets on the border"
182 );
183
184 ✗ std::vector<std::string> filenames;
185 ✗ if(
186 ✗ !GEO::CmdLine::parse(
187 ✗ argc, argv, filenames, "points_filename <cell_filename>"
188 )
189 ) {
190 ✗ return 1;
191 }
192
193 ✗ if(filenames.size() != 1 && filenames.size() != 2) {
194 ✗ return 1;
195 }
196
197 ✗ GEO::Mesh points;
198 ✗ GEO::MeshIOFlags flags;
199 ✗ flags.reset_element(GEO::MESH_FACETS);
200 ✗ flags.reset_element(GEO::MESH_CELLS);
201 ✗ GEO::mesh_load(filenames[0], points, flags);
202 ✗ GEO::mesh_repair(points);
203
204 ✗ double diag = GEO::bbox_diagonal(points);
205 ✗ double size = GEO::CmdLine::get_arg_percent("size",diag);
206 ✗ double shrink = GEO::CmdLine::get_arg_double("shrink");
207 ✗ bool border_only = GEO::CmdLine::get_arg_bool("border_only");
208
209 // Since we compute restricted Voronoi cells one cell at a
210 // time, the mesh argument of the restricted Voronoi diagram
211 // is not used.
212 ✗ GEO::Mesh dummy_mesh;
213
214 // Create a Delaunay API that encapsulates a Kd-tree
215 ✗ GEO::Delaunay_var delaunay = Delaunay::create(3,"NN");
216 ✗ delaunay->set_vertices(points.vertices.nb(), points.vertices.point_ptr(0));
217
218 GEO::RestrictedVoronoiDiagram_var RVD =
219 ✗ GEO::RestrictedVoronoiDiagram::create(delaunay, &dummy_mesh);
220
221 ✗ GEO::Mesh cell;
222 ✗ GEO::Mesh clipped;
223 ✗ GEO::Attribute<signed_index_t> facet_id;
224 ✗ if(border_only) {
225 ✗ facet_id.bind(clipped.facets.attributes(),"id");
226 }
227
228 ✗ if(filenames.size() == 2) {
229 ✗ mesh_load(filenames[1],cell);
230 } else {
231 ✗ initialize_mesh_with_box(cell);
232 }
233
234 ✗ if(!mesh_vertices_are_degree_3(cell)) {
235 ✗ Logger::err("RVC") << "Mesh vertices are not all of degree 3"
236 ✗ << std::endl;
237 ✗ exit(-1);
238 }
239
240 ✗ if(mesh_facets_are_planar(cell)) {
241 ✗ Logger::out("RVC") << "Mesh facets are planar (good)" << std::endl;
242 } else {
243 ✗ Logger::warn("RVC") << "Mesh facets are not planar" << std::endl;
244 }
245
246 ✗ std::ofstream out("RVC.obj");
247 ✗ index_t offset = 1;
248
249 index_t progress_divider =
250 ✗ (points.vertices.nb() > 10000) ? 100 : 1;
251
252 ✗ GEO::ProgressTask task("RVC.obj",points.vertices.nb()/progress_divider);
253 // For each point, create a cube centered on the point
254 // and clip it with the Voronoi cell of the point.
255 ✗ for(GEO::index_t i=0; i<points.vertices.nb(); ++i) {
256
257 ✗ if(!(i%progress_divider)) {
258 ✗ task.progress(i/progress_divider);
259 }
260
261 ✗ center_scale_mesh(cell, points.vertices.point(i), size);
262 ✗ RVD->compute_RVC(i,cell,clipped,facet_id.is_bound());
263
264 ✗ if(shrink != 1.0) {
265 ✗ shrink_mesh(clipped, shrink);
266 }
267
268 // Append the generated mesh to the output mesh.
269 ✗ for(index_t j=0; j<clipped.vertices.nb(); ++j) {
270 ✗ out << "v " << clipped.vertices.point(j) << std::endl;
271 }
272 ✗ for(index_t f=0; f<clipped.facets.nb(); ++f) {
273 ✗ if(border_only && facet_id[f] >= 0) {
274 ✗ continue;
275 }
276 ✗ out << "f ";
277 ✗ for(
278 ✗ index_t c=clipped.facets.corners_begin(f);
279 ✗ c<clipped.facets.corners_end(f); ++c
280 ) {
281 ✗ out << clipped.facet_corners.vertex(c) + offset << " ";
282 }
283 ✗ out << std::endl;
284 }
285 ✗ offset += clipped.vertices.nb();
286 }
287
288 ✗ return 0;
289 ✗ }
290