GCC Code Coverage Report

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