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 <exploragram/optimal_transport/sampling.h>
41			#include <geogram/voronoi/CVT.h>
42			#include <geogram/mesh/mesh.h>
43			#include <geogram/mesh/mesh_AABB.h>
44			#include <geogram/mesh/mesh_tetrahedralize.h>
45			#include <geogram/mesh/mesh_repair.h>
46			#include <geogram/mesh/mesh_geometry.h>
47			#include <geogram/mesh/mesh_reorder.h>
48			#include <geogram/basic/command_line.h>
49			#include <geogram/basic/permutation.h>
50			#include <geogram/basic/progress.h>
51
52			#ifdef GEOGRAM_WITH_VORPALINE
53			#include <vorpalib/voronoi/LpCVT.h>
54			#endif
55
56			namespace {
57			using namespace GEO;
58
59
60			/**
61			* \brief Reorders the points in a Centroidal Voronoi Tesselation
62			* in such a way that continguous index ranges correspond to
63			* multiple resolutions.
64			* \param[in,out] CVT the CentroidalVoronoiTesselation
65			* \param[out] levels sample indices that correspond to level l are
66			* in the range levels[l] (included) ... levels[l+1] (excluded)
67			* \param[in] ratio number of samples ratio between two consecutive
68			* levels
69			* \param[in] threshold minimum number of samples in a level
70			*/
71		✗	void BRIO_reorder(
72			CentroidalVoronoiTesselation& CVT,
73			vector<index_t>& levels,
74			double ratio,
75			index_t threshold
76			) {
77		✗	vector<index_t> sorted_indices;
78		✗	compute_BRIO_order(
79		✗	CVT.nb_points(), CVT.embedding(0), sorted_indices,
80		✗	CVT.dimension(), CVT.dimension(), threshold, ratio, &levels
81			);
82		✗	Permutation::apply(
83		✗	CVT.embedding(0), sorted_indices,
84		✗	index_t(CVT.dimension() * sizeof(double))
85			);
86		✗	}
87
88			/**
89			* \brief Internal implementation function for
90			* compute_hierarchical_sampling().
91			* \param[in,out] CVT the CentroidalVoronoiTesselation, initialized
92			* with the volume to be sampled. On output, it stores the samples
93			* \param[in] nb_samples total number of samples to generate
94			* \param[out] levels sample indices that correspond to level l are
95			* in the range levels[l] (included) ... levels[l+1] (excluded)
96			* \param[in] ratio number of samples ratio between two consecutive
97			* levels
98			* \param[in] threshold minimum number of samples in a level
99			* \param[in] b first element of the level to be generated
100			* \param[in] e one position past the last element of the
101			* level to be generated
102			* \param[in,out] points work vector allocated by caller,
103			* of size 3*nb_samples
104			*/
105		✗	void compute_hierarchical_sampling_recursive(
106			CentroidalVoronoiTesselation& CVT,
107			index_t nb_samples,
108			vector<index_t>& levels,
109			double ratio,
110			index_t threshold,
111			index_t b, index_t e,
112			vector<double>& points
113			) {
114		✗	index_t m = b;
115
116			// Recurse in [b...m) range
117		✗	if(e - b > threshold) {
118		✗	m = b + index_t(double(e - b) * ratio);
119		✗	compute_hierarchical_sampling_recursive(
120			CVT, nb_samples, levels, ratio, threshold, b, m, points
121			);
122			}
123
124			// Initialize random points in [m...e) range
125		✗	CVT.RVD()->compute_initial_sampling(&points[3 * m], e - m);
126
127			// Set the points in [b...e) range
128		✗	CVT.set_points(e - b, &points[0]);
129
130			// Lock [b...m) range
131			{
132		✗	for(index_t i = b; i < e; ++i) {
133		✗	if(i < m) {
134		✗	CVT.lock_point(i);
135			} else {
136		✗	CVT.unlock_point(i);
137			}
138			}
139			}
140
141		✗	Logger::div(
142		✗	std::string("Generating level ") +
143		✗	String::to_string(levels.size())
144			);
145
146		✗	Logger::out("Sample") << " generating a level with " << e - m
147		✗	<< " samples" << std::endl;
148
149			try {
150		✗	ProgressTask progress("Lloyd", 100);
151		✗	CVT.set_progress_logger(&progress);
152		✗	CVT.Lloyd_iterations(CmdLine::get_arg_uint("opt:nb_Lloyd_iter"));
153		✗	}
154		✗	catch(const TaskCanceled&) {
155		✗	}
156
157			try {
158		✗	ProgressTask progress("Newton", 100);
159		✗	CVT.set_progress_logger(&progress);
160		✗	CVT.Newton_iterations(CmdLine::get_arg_uint("opt:nb_Newton_iter"));
161		✗	}
162		✗	catch(const TaskCanceled&) {
163		✗	}
164
165		✗	levels.push_back(e);
166		✗	}
167
168			/**
169			* \brief Computes a hierarchical sampling of a volume.
170			* \param[in,out] CVT the CentroidalVoronoiTesselation, initialized
171			* with the volume to be sampled. On output, it stores the samples
172			* \param[in] nb_samples total number of samples to generate
173			* \param[out] levels sample indices that correspond to level l are
174			* in the range levels[l] (included) ... levels[l+1] (excluded)
175			* \param[in] ratio number of samples ratio between two consecutive
176			* levels
177			* \param[in] threshold minimum number of samples in a level
178			*/
179		✗	void compute_hierarchical_sampling(
180			CentroidalVoronoiTesselation& CVT,
181			index_t nb_samples,
182			vector<index_t>& levels,
183			double ratio = 0.125,
184			index_t threshold = 300
185			) {
186		✗	levels.push_back(0);
187		✗	vector<double> points(nb_samples * 3);
188		✗	compute_hierarchical_sampling_recursive(
189			CVT, nb_samples, levels, ratio, threshold,
190			0, nb_samples,
191			points
192			);
193		✗	CVT.unlock_all_points();
194		✗	}
195
196			/**
197			* \brief Computes a sampling of a volume.
198			* \param[in,out] CVT the CentroidalVoronoiTesselation, initialized
199			* with the volume to be sampled. On output, it stores the samples
200			* \param[in] nb_samples total number of samples to generate
201			*/
202		✗	void compute_single_level_sampling(
203			CentroidalVoronoiTesselation& CVT,
204			index_t nb_samples
205			) {
206
207		✗	CVT.compute_initial_sampling(nb_samples);
208
209			try {
210		✗	ProgressTask progress("Lloyd", 100);
211		✗	CVT.set_progress_logger(&progress);
212		✗	CVT.Lloyd_iterations(CmdLine::get_arg_uint("opt:nb_Lloyd_iter"));
213		✗	}
214		✗	catch(const TaskCanceled&) {
215		✗	}
216
217			try {
218		✗	ProgressTask progress("Newton", 100);
219		✗	CVT.set_progress_logger(&progress);
220		✗	CVT.Newton_iterations(CmdLine::get_arg_uint("opt:nb_Newton_iter"));
221		✗	}
222		✗	catch(const TaskCanceled&) {
223		✗	}
224		✗	}
225
226			/**
227			* \brief Projects the points of a volumetric sampling
228			* onto the border of the volume.
229			*/
230		✗	void project_sampling_on_border(
231			CentroidalVoronoiTesselation& CVT
232			) {
233			try {
234		✗	ProgressTask progress("Surf. Lloyd", 100);
235		✗	CVT.set_progress_logger(&progress);
236		✗	CVT.set_volumetric(false);
237		✗	CVT.Lloyd_iterations(
238		✗	CmdLine::get_arg_uint("opt:nb_Lloyd_iter") * 2
239			);
240		✗	}
241		✗	catch(const TaskCanceled&) {
242		✗	}
243
244			{
245
246			#ifdef GEOGRAM_WITH_VORPALINExxx
247			CVT.done_current();
248			{
249			LpCentroidalVoronoiTesselation LpCVT(
250			CVT.mesh(), 0
251			);
252			LpCVT.set_points(CVT.nb_points(), CVT.embedding(0));
253			try {
254			ProgressTask progress("LpCVT", 100);
255			LpCVT.set_progress_logger(&progress);
256			LpCVT.set_normal_anisotropy(5.0);
257			LpCVT.Newton_iterations(30, 7);
258			}
259			catch(const TaskCanceled&) {
260			}
261			CVT.set_points(LpCVT.nb_points(), LpCVT.embedding(0));
262			}
263			CVT.make_current();
264			#endif
265			}
266
267		✗	vector<double> mg(3 * CVT.nb_points());
268		✗	vector<double> m(CVT.nb_points());
269		✗	CVT.RVD()->compute_centroids(&mg[0], &m[0]);
270		✗	for(index_t i = 0; i < CVT.nb_points(); ++i) {
271		✗	if(m[i] == 0.0) {
272		✗	CVT.unlock_point(i);
273			} else {
274		✗	CVT.lock_point(i);
275			}
276			}
277
278		✗	CVT.set_volumetric(true);
279
280			try {
281		✗	ProgressTask progress("Relax. vol.", 100);
282		✗	CVT.set_progress_logger(&progress);
283		✗	CVT.Lloyd_iterations(
284		✗	CmdLine::get_arg_uint("opt:nb_Lloyd_iter") * 2
285			);
286		✗	}
287		✗	catch(const TaskCanceled&) {
288		✗	}
289		✗	}
290
291			}
292
293			namespace GEO {
294
295		✗	void recenter_mesh(const Mesh& M1, Mesh& M2) {
296			double xyzmin1[3];
297			double xyzmax1[3];
298			double xyzmin2[3];
299			double xyzmax2[3];
300			double xlat[3];
301		✗	get_bbox(M1, xyzmin1, xyzmax1);
302		✗	get_bbox(M2, xyzmin2, xyzmax2);
303		✗	for(coord_index_t c=0; c<3; ++c) {
304		✗	xlat[c] = 0.5*
305		✗	((xyzmin1[c] + xyzmax1[c]) - (xyzmin2[c] + xyzmax2[c]));
306			}
307		✗	for(index_t v=0; v<M2.vertices.nb(); ++v) {
308		✗	for(coord_index_t c=0; c<3; ++c) {
309		✗	M2.vertices.point_ptr(v)[c] += xlat[c];
310			}
311			}
312		✗	}
313
314		✗	double mesh_tets_volume(const Mesh& M) {
315		✗	double result = 0.0;
316		✗	for(index_t t = 0; t < M.cells.nb(); ++t) {
317		✗	result += Geom::tetra_volume<3>(
318			M.vertices.point_ptr(M.cells.tet_vertex(t, 0)),
319			M.vertices.point_ptr(M.cells.tet_vertex(t, 1)),
320			M.vertices.point_ptr(M.cells.tet_vertex(t, 2)),
321			M.vertices.point_ptr(M.cells.tet_vertex(t, 3))
322			);
323			}
324		✗	return result;
325			}
326
327		✗	void rescale_mesh(const Mesh& M1, Mesh& M2) {
328			double xyzmin[3];
329			double xyzmax[3];
330		✗	get_bbox(M2, xyzmin, xyzmax);
331		✗	double s = pow(mesh_tets_volume(M1)/mesh_tets_volume(M2), 1.0/3.0);
332		✗	for(unsigned int v=0; v<M2.vertices.nb(); ++v) {
333		✗	for(index_t c=0; c<3; ++c) {
334		✗	double gc = 0.5*(xyzmin[c]+xyzmax[c]);
335		✗	M2.vertices.point_ptr(v)[c] =
336		✗	gc + s * (M2.vertices.point_ptr(v)[c] - gc);
337			}
338			}
339		✗	}
340
341			enum DensityFunction {
342			DENSITY_X=0,
343			DENSITY_Y=1,
344			DENSITY_Z=2,
345			DENSITY_R,
346			DENSITY_SIN,
347			DENSITY_DIST
348			};
349
350		✗	void set_density(
351			Mesh& M, double mass1, double mass2, const std::string& function_str_in,
352			Mesh* density_distance_reference
353			) {
354		✗	std::string function_str = function_str_in;
355
356		✗	if(mass1 == mass2) {
357		✗	return;
358			}
359
360		✗	bool minus = false;
361		✗	if(function_str.length() > 1 && function_str[0] == '-') {
362		✗	minus = true;
363		✗	function_str = function_str.substr(1,function_str.length()-1);
364			}
365		✗	double density_pow = 1.0;
366			{
367		✗	std::size_t found = function_str.find('^');
368		✗	if(found != std::string::npos) {
369			std::string pow_str =
370		✗	function_str.substr(found+1, function_str.length()-found-1);
371		✗	density_pow = String::to_double(pow_str);
372		✗	function_str = function_str.substr(0,found);
373		✗	}
374			}
375
376		✗	Logger::out("OTM")
377			<< "Using density: "
378			<< (minus ? "-" : "+")
379		✗	<< function_str << "^"
380		✗	<< density_pow
381		✗	<< " rescaled to ("
382		✗	<< mass1 << "," << mass2
383		✗	<< ")"
384		✗	<< std::endl;
385
386			DensityFunction function;
387		✗	if(function_str == "X") {
388		✗	function = DENSITY_X;
389		✗	} else if(function_str == "Y") {
390		✗	function = DENSITY_Y;
391		✗	} else if(function_str == "Z") {
392		✗	function = DENSITY_Z;
393		✗	} else if(function_str == "R") {
394		✗	function = DENSITY_R;
395		✗	} else if(function_str == "sin") {
396		✗	function = DENSITY_SIN;
397		✗	} else if(function_str == "dist") {
398		✗	function = DENSITY_DIST;
399			} else {
400		✗	Logger::err("OTM") << function_str << ": no such density function"
401		✗	<< std::endl;
402		✗	return;
403			}
404
405		✗	Attribute<double> mass(M.vertices.attributes(),"weight");
406
407		✗	switch(function) {
408		✗	case DENSITY_X:
409			case DENSITY_Y:
410			case DENSITY_Z: {
411		✗	for(index_t v=0; v<M.vertices.nb(); ++v) {
412		✗	mass[v] = M.vertices.point_ptr(v)[index_t(function)];
413			}
414		✗	} break;
415		✗	case DENSITY_R: {
416			double xyz_min[3];
417			double xyz_max[3];
418		✗	get_bbox(M, xyz_min, xyz_max);
419		✗	for(index_t v=0; v<M.vertices.nb(); ++v) {
420		✗	double r=0;
421		✗	const double* p = M.vertices.point_ptr(v);
422		✗	for(coord_index_t c=0; c<3; ++c) {
423		✗	r += geo_sqr(p[c] - 0.5*(xyz_min[c] + xyz_max[c]));
424			}
425		✗	r = ::sqrt(r);
426		✗	mass[v] = r;
427			}
428		✗	} break;
429		✗	case DENSITY_SIN: {
430			double xyz_min[3];
431			double xyz_max[3];
432		✗	get_bbox(M, xyz_min, xyz_max);
433		✗	for(index_t v=0; v<M.vertices.nb(); ++v) {
434		✗	double f = 1.0;
435		✗	const double* p = M.vertices.point_ptr(v);
436		✗	for(coord_index_t c=0; c<3; ++c) {
437		✗	double coord =
438		✗	(p[c] - xyz_min[c]) / (xyz_max[c] - xyz_min[c]);
439		✗	f *= sin(coord * M_PI * 2.0 * 2.0);
440			}
441		✗	mass[v] = f;
442			}
443		✗	} break;
444		✗	case DENSITY_DIST: {
445		✗	if(density_distance_reference != nullptr) {
446		✗	MeshFacetsAABB AABB(*density_distance_reference);
447		✗	for(index_t v=0; v<M.vertices.nb(); ++v) {
448		✗	mass[v] =
449		✗	::sqrt(AABB.squared_distance(
450		✗	vec3(M.vertices.point_ptr(v)))
451			);
452			}
453		✗	} else {
454		✗	MeshFacetsAABB AABB(M);
455		✗	for(index_t v=0; v<M.vertices.nb(); ++v) {
456		✗	mass[v] = ::sqrt(
457		✗	AABB.squared_distance(vec3(M.vertices.point_ptr(v)))
458			);
459			}
460		✗	}
461		✗	} break;
462			}
463
464			// Compute min and max mass
465		✗	double mass_min = Numeric::max_float64();
466		✗	double mass_max = Numeric::min_float64();
467		✗	for(index_t v=0; v<M.vertices.nb(); ++v) {
468		✗	mass_min = std::min(mass_min, mass[v]);
469		✗	mass_max = std::max(mass_max, mass[v]);
470			}
471
472			// Normalize mass, apply power, and rescale to (mass1 - mass2)
473		✗	for(index_t v=0; v<M.vertices.nb(); ++v) {
474		✗	double f = (mass[v] - mass_min) / (mass_max - mass_min);
475		✗	if(minus) {
476		✗	f = 1.0 - f;
477			}
478		✗	f = ::pow(f,density_pow);
479		✗	mass[v] = mass1 + f*(mass2 - mass1);
480			}
481		✗	}
482
483		✗	void sample(
484			CentroidalVoronoiTesselation& CVT,
485			index_t nb_points, bool project_on_border,
486			bool BRIO, bool multilevel, double ratio,
487			vector<index_t>* levels_out
488			) {
489		✗	vector<index_t> levels;
490		✗	multilevel = multilevel | BRIO;
491		✗	if(CmdLine::get_arg_bool("RVD_iter") && multilevel) {
492		✗	Logger::warn("OTM") << "Deactivating multilevel mode" << std::endl;
493		✗	Logger::warn("OTM") << "(because RVD_iter is set)" << std::endl;
494		✗	multilevel = false;
495			}
496		✗	if(multilevel) {
497		✗	if(BRIO) {
498		✗	compute_single_level_sampling(CVT, nb_points);
499		✗	BRIO_reorder(CVT, levels, ratio, 300);
500			} else {
501		✗	compute_hierarchical_sampling(
502			CVT, nb_points,levels,ratio
503			);
504			}
505			} else {
506		✗	compute_single_level_sampling(CVT, nb_points);
507			}
508		✗	if(levels_out != nullptr) {
509		✗	*levels_out = levels;
510			}
511		✗	if(project_on_border) {
512		✗	project_sampling_on_border(CVT);
513			}
514		✗	}
515
516
517			}
518