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/parameterization/mesh_segmentation.h>
41			#include <geogram/mesh/mesh.h>
42			#include <geogram/mesh/mesh_geometry.h>
43			#include <geogram/mesh/mesh_io.h>
44			#include <geogram/mesh/mesh_manifold_harmonics.h>
45			#include <geogram/voronoi/CVT.h>
46			#include <geogram/voronoi/RVD.h>
47			#include <geogram/voronoi/RVD_callback.h>
48			#include <geogram/voronoi/generic_RVD_polygon.h>
49			#include <geogram/points/principal_axes.h>
50			#include <geogram/numerics/matrix_util.h>
51			#include <geogram/basic/numeric.h>
52
53			#include <deque>
54			#include <stack>
55
56			/**************************************************************************
57			**** SMOOTH PARTITION ****
58			**************************************************************************/
59
60			namespace {
61			using namespace GEO;
62
63
64			/**
65			* \brief Iterator that traverses all facets incident to
66			* a given internal vertex.
67			* \param[in] M a reference to the mesh
68			* \param[in] f a facet incident to the vertex
69			* \param[in] lv the local index of the vertex in \p f
70			* \param[in] CB the function or lambda to be called for
71			* all facets incident to the vertex
72			*/
73		✗	inline void for_each_facet_around_internal_vertex(
74			const Mesh& M, index_t f, index_t lv,
75			std::function<void(index_t, index_t)> CB
76			) {
77			index_t v = M.facets.vertex(f,lv);
78			index_t cur_f = f;
79			index_t cur_lv = lv;
80			index_t count = 0;
81			do {
82		✗	CB(cur_f, cur_lv);
83			cur_f = M.facets.adjacent(cur_f,cur_lv);
84		✗	geo_assert(cur_f != NO_INDEX);
85		✗	cur_lv = M.facets.find_vertex(cur_f, v);
86		✗	geo_assert(cur_lv != NO_INDEX);
87		✗	++count;
88		✗	geo_assert(count < 10000); // sanity check (are we looping forever?)
89		✗	} while(cur_f != f);
90		✗	}
91
92			/**
93			* \brief Utility class for mesh partition smoothing
94			* \details Determines whether chart indices in the facets incident to
95			* a given vertex could be changed to reduce chart border length
96			*/
97			class SmoothVertex {
98			public:
99
100		✗	SmoothVertex() {
101			}
102
103			/**
104			* \brief SmoothVertex constructor
105			* \param[in] M a reference to the mesh
106			* \param[in] chart a reference to the partition facet attribute
107			* \param[in] f a facet incident to the vertex
108			* \param[in] lv the local index of the vertex in \p f
109			*/
110		✗	SmoothVertex(
111			const Mesh& M,
112			Attribute<index_t>& chart,
113			index_t f,
114			index_t lv
115		✗	) {
116		✗	f_ = f;
117		✗	lv_ = lv;
118		✗	v_ = M.facets.vertex(f_,lv_);
119		✗	is_valid_ = true;
120		✗	chart_id_ = NO_INDEX;
121
122			// Get chart1 and chart2 Ids
123		✗	index_t chart1 = NO_INDEX;
124		✗	index_t chart2 = NO_INDEX;
125		✗	index_t prev_chart = NO_INDEX;
126		✗	for_each_facet_around_internal_vertex(
127			M,f,lv,
128		✗	[&](index_t cur_f, index_t cur_lv) {
129			geo_argused(cur_lv);
130		✗	if(chart1 == NO_INDEX) {
131		✗	chart1 = chart[cur_f];
132		✗	} else if(chart[cur_f] != chart1 && chart2 == NO_INDEX) {
133		✗	chart2 = chart[cur_f];
134			}
135		✗	prev_chart = chart[cur_f];
136		✗	}
137			);
138
139			// Test that vertex is incident to at most
140			// two charts and that chart id changes at most
141			// twice when turning around the vertex
142		✗	index_t nb_change=0;
143		✗	index_t nb_chart1=0;
144		✗	index_t nb_chart2=0;
145		✗	for_each_facet_around_internal_vertex(
146			M,f,lv,
147		✗	[&](index_t cur_f, index_t cur_lv) {
148			geo_argused(cur_lv);
149		✗	if(chart[cur_f] != prev_chart) {
150		✗	++nb_change;
151			}
152		✗	prev_chart = chart[cur_f];
153		✗	if(chart[cur_f] == chart1) {
154		✗	++nb_chart1;
155		✗	} else if(chart[cur_f] == chart2) {
156		✗	++nb_chart2;
157			} else {
158		✗	is_valid_ = false;
159			}
160		✗	}
161			);
162		✗	is_valid_ = is_valid_ &&
163		✗	(chart1 != NO_INDEX) &&
164		✗	(chart2 != NO_INDEX) ;
165		✗	is_valid_ = is_valid_ && (nb_change <= 2);
166		✗	if(!is_valid_) {
167		✗	return;
168			}
169		✗	chart_id_ = (nb_chart1 > nb_chart2) ? chart1 : chart2;
170
171			// Compute delta len
172		✗	delta_len_ = 0.0;
173		✗	for_each_facet_around_internal_vertex(
174			M,f,lv,
175		✗	[&](index_t cur_f, index_t cur_lv) {
176		✗	index_t N = M.facets.nb_vertices(cur_f);
177		✗	index_t prev_lv = (cur_lv == 0) ? (N-1) : cur_lv - 1;
178		✗	index_t next_lv = (cur_lv == N-1) ? 0 : cur_lv + 1;
179			index_t prev_v = M.facets.vertex(cur_f, prev_lv);
180			index_t v = M.facets.vertex(cur_f, cur_lv);
181			index_t next_v = M.facets.vertex(cur_f, next_lv);
182		✗	vec3 prev_p = M.vertices.point(prev_v);
183		✗	vec3 p = M.vertices.point(v);
184		✗	vec3 next_p = M.vertices.point(next_v);
185		✗	if(chart[cur_f] !=
186			chart[M.facets.adjacent(cur_f, cur_lv)]) {
187		✗	delta_len_ += Geom::distance(p, next_p);
188			}
189		✗	if(chart[cur_f] != chart_id_) {
190		✗	delta_len_ -= Geom::distance(prev_p, p);
191			}
192		✗	}
193			);
194		✗	is_valid_ = is_valid_ && delta_len_ > 0;
195			}
196
197			/**
198			* \brief used to sort a vector of SmoothVertex and
199			* smooth them in order of priority.
200			*/
201			bool operator<(const SmoothVertex& rhs) const {
202		✗	return (delta_len_ < rhs.delta_len_) ;
203			}
204
205			/**
206			* \brief Changes chart ids in the facets incident to
207			* this SmoothVertex.
208			* \param[in] M a reference to the mesh
209			* \param[in,out] chart a reference to the segmentation facet attribute
210			* \param[in,out] v_is_locked a vector of booleans that forbids
211			* smoothing the neighbors of a vertex that was already smoothed
212			*/
213		✗	bool apply(
214			Mesh& M,
215			Attribute<index_t>& chart,
216			std::vector<bool>& v_is_locked
217			) {
218		✗	if(v_is_locked[M.facets.vertex(f_, lv_)]) {
219			return false;
220			}
221		✗	for_each_facet_around_internal_vertex(
222			M,f_,lv_,
223		✗	[&](index_t cur_f, index_t cur_lv) {
224		✗	chart[cur_f] = chart_id_;
225		✗	index_t N = M.facets.nb_vertices(cur_f);
226		✗	index_t next_lv = (cur_lv == N-1) ? 0 : cur_lv + 1;
227		✗	v_is_locked[M.facets.vertex(cur_f,next_lv)] = true;
228		✗	}
229			);
230		✗	return true;
231			}
232
233			/**
234			* \brief Tests whether this SmoothVertex can be applied.
235			* \details A SmoothVertex cannot be applied if it is incident
236			* to more than two charts or if one of its neighbors are incident
237			* to more than two charts, or if chart id changes more than twice
238			* when turning around it.
239			* \retval true if it can be applied, false otherwise
240			*/
241			bool is_valid() const {
242		✗	return is_valid_ ;
243			}
244
245			public:
246			index_t f_;
247			index_t lv_;
248			index_t v_;
249			index_t chart_id_;
250			double delta_len_;
251			bool is_valid_;
252			};
253
254			/*****************************************************/
255
256			/**
257			* \brief Smoothes a mesh segmentation by changing chart ids in order
258			* to reduce total chart border length
259			* \param[in,out] M a reference to a mesh
260			* \param[in] nb_iter number of iterations of partition smoothing
261			*/
262		✗	void mesh_smooth_segmentation(Mesh& M, index_t nb_iter=10) {
263
264			// For each vertex, store one facet incident to that vertex
265			vector<index_t> v_to_f(M.vertices.nb(), NO_INDEX);
266		✗	for(index_t c: M.facet_corners) {
267		✗	v_to_f[M.facet_corners.vertex(c)] =
268			M.facet_corners.adjacent_facet(c) ;
269			}
270
271		✗	vector<bool> v_on_border(M.vertices.nb(), false);
272		✗	for(index_t c: M.facet_corners) {
273			index_t v = M.facet_corners.vertex(c);
274		✗	if(M.facet_corners.adjacent_facet(c) == NO_INDEX) {
275		✗	v_on_border[v] = true;
276			}
277			}
278
279			// Remove vertices on border and vertices adjacent to a vertex
280			// on border
281		✗	for(index_t f: M.facets) {
282		✗	for(index_t c1: M.facets.corners(f)) {
283			index_t v1 = M.facet_corners.vertex(c1);
284			index_t c2 = M.facets.next_corner_around_facet(f,c1);
285			index_t v2 = M.facet_corners.vertex(c2);
286			if(
287		✗	M.facet_corners.adjacent_facet(c1) == NO_INDEX ||
288		✗	v_on_border[v2]
289			) {
290		✗	v_to_f[v1] = NO_INDEX;
291			}
292			}
293			}
294
295		✗	Attribute<index_t> chart(M.facets.attributes(),"chart");
296			vector<bool> v_is_locked;
297			vector<SmoothVertex> smooth_vertices;
298
299		✗	for(index_t i=0; i<nb_iter; ++i) {
300		✗	smooth_vertices.resize(0);
301		✗	v_is_locked.assign(M.vertices.nb(),false);
302		✗	for(index_t v: M.vertices) {
303			// skip vertices on border
304			// and vertices adjacent to vertices on border
305			// and isolated vertices
306		✗	if(v_to_f[v] == NO_INDEX) {
307		✗	continue;
308			}
309			SmoothVertex sv(
310			M, chart, v_to_f[v], M.facets.find_vertex(v_to_f[v],v)
311		✗	);
312		✗	if(sv.is_valid()) {
313			smooth_vertices.push_back(sv);
314			}
315			}
316		✗	std::sort(smooth_vertices.begin(), smooth_vertices.end()) ;
317			bool changed = false;
318		✗	for(SmoothVertex& sv: smooth_vertices) {
319		✗	if(sv.apply(M, chart, v_is_locked)) {
320			changed = true;
321			}
322			}
323		✗	if(!changed) {
324			break;
325			}
326			}
327		✗	}
328
329
330			/**
331			* \brief Makes sure that each chart of the segmentation is
332			* connected.
333			* \details Segmentation is stored in the "chart" facet attribute.
334			* Generates a new chart id for each connected component of
335			* the input charts.
336			* \return number of charts
337			*/
338		✗	index_t mesh_postprocess_segmentation(Mesh& M, bool verbose=false) {
339			Attribute<index_t> chart;
340		✗	chart.bind_if_is_defined(M.facets.attributes(),"chart");
341		✗	geo_assert(chart.is_bound());
342
343			// Mark facets as non-visited by negating chart id
344		✗	for(index_t f: M.facets) {
345		✗	signed_index_t id = -signed_index_t(chart[f])-1;
346		✗	chart[f] = index_t(id);
347			}
348
349			std::stack<index_t> S;
350			index_t cur_chart = 0;
351		✗	for(index_t f: M.facets) {
352		✗	index_t f_chart = chart[f];
353		✗	if(signed_index_t(f_chart) < 0) {
354		✗	chart[f] = cur_chart;
355			S.push(f);
356		✗	while(!S.empty()) {
357		✗	index_t cur_f = S.top();
358			S.pop();
359		✗	for(index_t e=0; e<M.facets.nb_vertices(cur_f); ++e) {
360		✗	index_t neigh_f = M.facets.adjacent(cur_f,e);
361			if(
362		✗	neigh_f != NO_INDEX &&
363		✗	chart[neigh_f] == f_chart
364			) {
365		✗	chart[neigh_f] = cur_chart;
366			S.push(neigh_f);
367			}
368			}
369			}
370		✗	++cur_chart;
371			}
372			}
373		✗	if(verbose) {
374		✗	Logger::out("Segmentation") << cur_chart << " charts" << std::endl;
375			}
376		✗	return cur_chart;
377			}
378			}
379
380			/***************************************************************************
381			***** MESH_SEGMENT CVT *****
382			***************************************************************************/
383
384			namespace {
385			using namespace GEO;
386
387			/**
388			* \brief Helper class for mesh_segment()
389			* \details In a restricted Voronoi diagram, finds for each facet
390			* of the mesh the id of the Voronoi cell that has the largest
391			* intersection with the mesh.
392			*/
393			class PartitionCB : public RVDPolygonCallback {
394			public:
395		✗	PartitionCB(const Mesh* mesh) : mesh_(mesh) {
396			}
397
398		✗	void begin() override {
399		✗	facet_seed_.assign(mesh_->facets.nb(), NO_INDEX);
400		✗	facet_RVD_area_.assign(mesh_->facets.nb(), 0.0);
401		✗	}
402
403		✗	void end() override {
404		✗	Attribute<index_t> chart(mesh_->facets.attributes(), "chart");
405		✗	for(index_t f:mesh_->facets) {
406		✗	chart[f] = facet_seed_[f];
407			}
408		✗	}
409
410		✗	void operator() (
411			index_t v,
412			index_t t,
413			const GEOGen::Polygon& C
414			) const override {
415		✗	double A = area(C);
416		✗	if(facet_seed_[t] == NO_INDEX || A > facet_RVD_area_[t]) {
417		✗	facet_seed_[t] = v;
418		✗	facet_RVD_area_[t] = A;
419			}
420		✗	}
421
422		✗	double area(const GEOGen::Polygon& C) const {
423			double result = 0.0;
424			vec3 p0(C.vertex(0).point());
425		✗	for(index_t i=1; i<C.nb_vertices()-1; ++i) {
426			vec3 pi(C.vertex(i).point());
427		✗	vec3 pj(C.vertex(i+1).point());
428		✗	result += Geom::triangle_area(p0,pi,pj);
429			}
430		✗	return result;
431			}
432
433			private:
434			const Mesh* mesh_;
435			mutable vector<index_t> facet_seed_;
436			mutable vector<double> facet_RVD_area_;
437			};
438			}
439
440			/***************************************************************************
441			***** MESH_SEGMENT PPAL AXIS *****
442			***************************************************************************/
443
444			namespace {
445			using namespace GEO;
446
447			/**
448			* \brief Splits a chart along one of its principal axis
449			* \details Greedily grow two charts from the two facets that
450			* are furthest away along the specified axis
451			* \param[in] M a reference to the Mesh
452			* \param[in] axis one of 0,1,2
453			* \retval true if chart boundary touches a border
454			* \retval false otherwise
455			*/
456		✗	bool split_chart_along_principal_axis(Mesh & M, index_t axis) {
457		✗	Attribute<index_t> chart(M.facets.attributes(), "chart");
458
459		✗	PrincipalAxes3d axes ;
460		✗	axes.begin() ;
461		✗	for(index_t f: M.facets) {
462		✗	for(index_t lv=0; lv<M.facets.nb_vertices(f); ++lv) {
463			index_t v = M.facets.vertex(f,lv);
464		✗	axes.add_point(M.vertices.point(v));
465			}
466			}
467		✗	axes.end() ;
468			vec3 center = axes.center() ;
469		✗	vec3 X = axes.axis(axis) ;
470
471			vector<double> X_coord(M.facets.nb());
472
473		✗	for(index_t f: M.facets) {
474		✗	X_coord[f] = dot(X,(Geom::mesh_facet_center(M,f) - center));
475			}
476
477			vector<double> axis_coord = X_coord;
478		✗	std::sort(axis_coord.begin(), axis_coord.end());
479		✗	double X_cutoff = axis_coord[axis_coord.size()/2];
480
481		✗	for(index_t f: M.facets) {
482		✗	chart[f] = (X_coord[f] > X_cutoff);;
483			}
484
485			// Test whether chart boundary touches mesh border
486			// (which is what we want if we split a cylindroid
487			// or a sockoid).
488		✗	for(index_t f: M.facets) {
489			index_t N = M.facets.nb_vertices(f);
490		✗	for(index_t e1=0; e1<N; ++e1) {
491		✗	index_t e2 = (e1+1)%N;
492
493			index_t adj1 = M.facets.adjacent(f,e1);
494			index_t adj2 = M.facets.adjacent(f,e2);
495
496			if(
497		✗	adj1 == NO_INDEX &&
498		✗	adj2 != NO_INDEX &&
499		✗	chart[adj2] != chart[f]
500			) {
501			return true;
502			}
503
504			if(
505		✗	adj2 == NO_INDEX &&
506		✗	adj1 != NO_INDEX &&
507		✗	chart[adj1] != chart[f]
508			) {
509			return true;
510			}
511
512			}
513
514			}
515
516			return false;
517			}
518			}
519
520			namespace GEO {
521
522		✗	index_t mesh_segment(
523			Mesh& M, MeshSegmenter segmenter, index_t nb_segments, bool verbose
524			) {
525		✗	geo_assert(M.facets.are_simplices());
526
527			double anisotropy = 0.0;
528		✗	if(segmenter == SEGMENT_GEOMETRIC_VSA_L12) {
529		✗	anisotropy = bbox_diagonal(M) * 100.0;
530			}
531			index_t dimension = 0;
532			index_t nb_manifold_harmonics=0;
533
534			switch(segmenter) {
535			case SEGMENT_GEOMETRIC_VSA_L2: break;
536			case SEGMENT_GEOMETRIC_VSA_L12: break;
537			case SEGMENT_INERTIA_AXIS: break;
538			case SEGMENT_SPECTRAL_8: dimension=8; break;
539			case SEGMENT_SPECTRAL_20: dimension=20; break;
540			case SEGMENT_SPECTRAL_100: dimension=100; break;
541			}
542
543			Attribute<double> geom_bkp;
544
545		✗	if(segmenter == SEGMENT_INERTIA_AXIS) {
546			// Pick the axis such that the segmentation obtained
547			// by splitting along it has a chart boundary that
548			// touches the mesh boundary.
549		✗	for(index_t axis=0; axis<3; ++axis) {
550		✗	if(split_chart_along_principal_axis(M, 2-axis)) {
551			break;
552			}
553			}
554		✗	mesh_smooth_segmentation(M);
555		✗	return mesh_postprocess_segmentation(M,verbose);
556			}
557
558		✗	if(dimension != 0) {
559		✗	nb_manifold_harmonics = dimension+20;
560		✗	geom_bkp.create_vector_attribute(
561		✗	M.vertices.attributes(), "bkp", 3
562			);
563		✗	for(index_t v: M.vertices) {
564		✗	geom_bkp[3*v] = M.vertices.point_ptr(v)[0];
565		✗	geom_bkp[3*v+1] = M.vertices.point_ptr(v)[1];
566		✗	geom_bkp[3*v+2] = M.vertices.point_ptr(v)[2];
567			}
568		✗	mesh_compute_manifold_harmonics(
569			M, nb_manifold_harmonics,
570		✗	FEM_P1_LUMPED, "eigen", 0.0, true
571			);
572			Attribute<double> eigen(
573		✗	M.vertices.attributes(), "eigen"
574			);
575		✗	M.vertices.set_dimension(dimension);
576		✗	for(index_t v: M.vertices) {
577		✗	for(index_t mh=0; mh<dimension; ++mh) {
578		✗	M.vertices.point_ptr(v)[mh] =
579		✗	eigen[nb_manifold_harmonics*v + mh + 1];
580			}
581			}
582		✗	eigen.destroy();
583		✗	} else if(anisotropy != 0.0) {
584		✗	compute_normals(M);
585			// smooth normals --------------.
586			// v
587		✗	simple_Laplacian_smooth(M, 3, true);
588		✗	set_anisotropy(M,anisotropy*0.02);
589			}
590
591		✗	CentroidalVoronoiTesselation CVT(&M);
592		✗	CVT.compute_initial_sampling(nb_segments);
593		✗	if(verbose) {
594		✗	Logger::out("RVD") << "Optimizing CVT" << std::endl;
595			}
596		✗	CVT.Lloyd_iterations(30);
597		✗	CVT.Newton_iterations(10);
598		✗	PartitionCB CB(&M);
599		✗	CVT.RVD()->for_each_polygon(CB);
600
601		✗	if(nb_manifold_harmonics != 0) {
602		✗	M.vertices.set_dimension(3);
603		✗	for(index_t v: M.vertices) {
604		✗	M.vertices.point_ptr(v)[0] = geom_bkp[3*v];
605		✗	M.vertices.point_ptr(v)[1] = geom_bkp[3*v+1];
606		✗	M.vertices.point_ptr(v)[2] = geom_bkp[3*v+2];
607			}
608		✗	geom_bkp.destroy();
609		✗	} else if(anisotropy != 0.0) {
610		✗	M.vertices.set_dimension(3);
611			}
612
613		✗	mesh_smooth_segmentation(M);
614		✗	return mesh_postprocess_segmentation(M,verbose);
615			}
616			}
617
618			/***************************************************************************/
619