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_param_packer.h>
41			#include <geogram/parameterization/mesh_atlas_maker.h>
42			#include <geogram/parameterization/mesh_segmentation.h>
43			#include <geogram/mesh/mesh.h>
44			#include <geogram/mesh/mesh_geometry.h>
45			#include <geogram/points/principal_axes.h>
46			#include <geogram/basic/logger.h>
47			#include <geogram/numerics/matrix_util.h>
48			#include <geogram/third_party/xatlas/xatlas.h>
49			#include <algorithm>
50			#include <stack>
51			#include <math.h>
52
53			/**************************************************************
54			**** XATLAS ****
55			**************************************************************/
56
57			namespace {
58			using namespace GEO;
59
60			/**
61			* \brief Interface between geogram and xatlas to represent
62			* charts.
63			*/
64			class XAtlasChart {
65			public:
66			// TODO: create all charts simultaneously (because here
67			// we are O(m*n), not good, there can be many charts...)
68			// ... (seems to be OK though, even on large meshes)
69		✗	void initialize(const Mesh& M, index_t chart_id) {
70		✗	Attribute<index_t> chart;
71		✗	chart.bind_if_is_defined(M.facets.attributes(), "chart");
72		✗	geo_assert(chart.is_bound());
73		✗	Attribute<double> tex_coord;
74		✗	tex_coord.bind_if_is_defined(
75		✗	M.facet_corners.attributes(), "tex_coord"
76			);
77		✗	geo_assert(tex_coord.is_bound());
78			Attribute<index_t> xatlas_vertex_index(
79		✗	M.vertices.attributes(), "xatlas_index"
80		✗	);
81			Attribute<index_t> facet_corner_xatlas_vertex_index(
82		✗	M.facet_corners.attributes(), "xatlas_index"
83		✗	);
84
85
86		✗	for(index_t v: M.vertices) {
87		✗	xatlas_vertex_index[v] = NO_INDEX;
88			}
89		✗	index_t cur_xatlas_vertex = 0;
90		✗	for(index_t f: M.facets) {
91		✗	if(chart[f] == chart_id) {
92		✗	for(index_t c: M.facets.corners(f)) {
93		✗	index_t v = M.facet_corners.vertex(c);
94		✗	if(xatlas_vertex_index[v] == NO_INDEX) {
95		✗	xatlas_vertex_index[v] = cur_xatlas_vertex;
96		✗	const double* xyz = M.vertices.point_ptr(v);
97		✗	xyz_.push_back(float(xyz[0]));
98		✗	xyz_.push_back(float(xyz[1]));
99		✗	xyz_.push_back(float(xyz[2]));
100		✗	const double* uv = &tex_coord[2*c];
101		✗	uv_.push_back(float(uv[0]));
102		✗	uv_.push_back(float(uv[1]));
103		✗	++cur_xatlas_vertex;
104			}
105			}
106			}
107			}
108		✗	for(index_t f: M.facets) {
109		✗	if(chart[f] != chart_id) {
110		✗	continue;
111			}
112
113		✗	index_t c1 = M.facets.corners_begin(f);
114		✗	for(index_t c2 = c1+1;
115		✗	c2+1 < M.facets.corners_end(f); ++c2
116			) {
117		✗	index_t v1 = M.facet_corners.vertex(c1);
118		✗	index_t v2 = M.facet_corners.vertex(c2);
119		✗	index_t v3 = M.facet_corners.vertex(c2+1);
120
121		✗	triangles_.push_back(xatlas_vertex_index[v1]);
122		✗	triangles_.push_back(xatlas_vertex_index[v2]);
123		✗	triangles_.push_back(xatlas_vertex_index[v3]);
124			}
125			}
126
127		✗	for(index_t f: M.facets) {
128		✗	if(chart[f] != chart_id) {
129		✗	continue;
130			}
131		✗	for(index_t c: M.facets.corners(f)) {
132		✗	index_t v = M.facet_corners.vertex(c);
133		✗	facet_corner_xatlas_vertex_index[c] =
134		✗	xatlas_vertex_index[v];
135			}
136			}
137
138		✗	mesh_decl_.vertexCount = uint32_t(cur_xatlas_vertex);
139		✗	mesh_decl_.vertexPositionData = xyz_.data();
140		✗	mesh_decl_.vertexPositionStride = sizeof(float)*3;
141		✗	mesh_decl_.vertexUvData = uv_.data();
142		✗	mesh_decl_.vertexUvStride = sizeof(float)*2;
143		✗	mesh_decl_.indexCount = uint32_t(triangles_.size());
144		✗	mesh_decl_.indexData = triangles_.data();
145		✗	mesh_decl_.indexOffset = 0;
146		✗	mesh_decl_.indexFormat = xatlas::IndexFormat::UInt32;
147
148		✗	uv_mesh_decl_.vertexCount = uint32_t(cur_xatlas_vertex);
149		✗	uv_mesh_decl_.vertexStride = sizeof(float)*2;
150		✗	uv_mesh_decl_.vertexUvData = uv_.data();
151		✗	uv_mesh_decl_.indexCount = uint32_t(triangles_.size());
152		✗	uv_mesh_decl_.indexData = triangles_.data();
153		✗	uv_mesh_decl_.indexOffset = 0;
154		✗	uv_mesh_decl_.indexFormat = xatlas::IndexFormat::UInt32;
155		✗	uv_mesh_decl_.faceMaterialData = nullptr;
156		✗	uv_mesh_decl_.rotateCharts = false;
157		✗	}
158
159			xatlas::MeshDecl mesh_decl_;
160			xatlas::UvMeshDecl uv_mesh_decl_;
161			vector<float> xyz_;
162			vector<float> uv_;
163			vector<index_t> triangles_;
164			};
165
166
167			/**
168			* \brief Allocator function for geogram
169			* \details If we call realloc() with size == 0 to free memory
170			* (as xatlas does), then valgrind thinks that there is a memory leak,
171			* so I re-interpret the initial intent of the caller here,
172			* to make the memory debugger happy.
173			*/
174		✗	void* my_realloc(void* ptr, size_t size) {
175		✗	if(size == 0) {
176		✗	if(ptr != nullptr) {
177		✗	free(ptr);
178			}
179		✗	return nullptr;
180			}
181		✗	if(ptr == nullptr) {
182		✗	return malloc(size);
183			}
184		✗	return realloc(ptr, size);
185			}
186			}
187
188			namespace GEO {
189
190		✗	void pack_atlas_using_xatlas(Mesh& mesh) {
191
192		✗	Attribute<double> tex_coord;
193
194		✗	tex_coord.bind_if_is_defined(
195		✗	mesh.facet_corners.attributes(), "tex_coord"
196			);
197
198		✗	if(!tex_coord.is_bound()) {
199		✗	Logger::err("Atlas") << "Mesh has no texture coordinates"
200		✗	<< std::endl;
201		✗	return;
202			}
203
204
205			// Step 1: compute chart indices
206		✗	Attribute<index_t> chart(mesh.facets.attributes(), "chart");
207		✗	index_t nb_charts = mesh_get_charts(mesh);
208
209		✗	xatlas::Atlas* atlas = nullptr;
210
211			// Step 2: copy to xatlas and pack
212			{
213		✗	std::vector<XAtlasChart> charts;
214		✗	charts.resize(nb_charts);
215		✗	for(index_t c=0; c<nb_charts; ++c) {
216		✗	charts[c].initialize(mesh, c);
217			}
218
219		✗	atlas = xatlas::Create();
220		✗	for(index_t c=0; c<nb_charts; ++c) {
221		✗	xatlas::AddMeshError::Enum error = xatlas::AddUvMesh(
222		✗	atlas, charts[c].uv_mesh_decl_
223			);
224		✗	if (error != xatlas::AddMeshError::Success) {
225		✗	Logger::err("Packer")
226			<< "XAtlas error: "
227		✗	<< xatlas::StringForEnum(error)
228		✗	<< std::endl;
229			}
230			}
231			Attribute<index_t> xatlas_vertex_index(
232		✗	mesh.vertices.attributes(), "xatlas_index"
233		✗	);
234		✗	xatlas_vertex_index.destroy();
235		✗	}
236
237			// Step 3: pack
238			{
239		✗	xatlas::SetAlloc(my_realloc);
240		✗	xatlas::PackOptions packerOptions;
241		✗	packerOptions.padding = 1;
242		✗	xatlas::PackCharts(atlas, packerOptions);
243
244		✗	double u_min = Numeric::max_float64();
245		✗	double v_min = Numeric::max_float64();
246		✗	double u_max = Numeric::min_float64();
247		✗	double v_max = Numeric::min_float64();
248
249			// Get uv bbox
250			{
251			Attribute<index_t> facet_corner_xatlas_vertex_index(
252		✗	mesh.facet_corners.attributes(), "xatlas_index"
253		✗	);
254		✗	for(index_t f: mesh.facets) {
255		✗	index_t chart_id = chart[f];
256		✗	for(index_t c: mesh.facets.corners(f)) {
257		✗	index_t v = facet_corner_xatlas_vertex_index[c];
258		✗	const xatlas::Vertex& vertex =
259		✗	atlas->meshes[chart_id].vertexArray[v];
260		✗	double U = double(vertex.uv[0]);
261		✗	double V = double(vertex.uv[1]);
262		✗	u_min = std::min(u_min,U);
263		✗	v_min = std::min(v_min,V);
264		✗	u_max = std::max(u_max,U);
265		✗	v_max = std::max(v_max,V);
266			}
267			}
268		✗	}
269
270			// Scale texture coordinates to [0,1]
271			{
272		✗	double scale = 1.0 / std::max(u_max-u_min,v_max-v_min);
273			Attribute<index_t> facet_corner_xatlas_vertex_index(
274		✗	mesh.facet_corners.attributes(), "xatlas_index"
275		✗	);
276		✗	for(index_t f: mesh.facets) {
277		✗	index_t chart_id = chart[f];
278		✗	for(index_t c: mesh.facets.corners(f)) {
279		✗	index_t v = facet_corner_xatlas_vertex_index[c];
280		✗	const xatlas::Vertex& vertex =
281		✗	atlas->meshes[chart_id].vertexArray[v];
282		✗	tex_coord[2*c] = scale*(double(vertex.uv[0])-u_min);
283		✗	tex_coord[2*c+1] = scale*(double(vertex.uv[1])-v_min);
284			}
285			}
286		✗	facet_corner_xatlas_vertex_index.destroy();
287		✗	chart.destroy();
288		✗	}
289		✗	xatlas::Destroy(atlas);
290			}
291		✗	}
292
293			}
294
295			/**************************************************************
296			**** TETRIS PACKER ****
297			**************************************************************/
298
299			namespace GEO {
300			/**
301			* \brief A piece of a mesh.
302			* \details Stores a list of facet indices. The mesh it belongs
303			* to is supposed to have an Attribute<index_t> attached to the
304			* facets and indicating the id of the chart for each facet.
305			*/
306			struct Chart {
307
308			/**
309			* \brief Chart constructor.
310			* \param[in] mesh_in a reference to a surface mesh.
311			* \param[in] id_in the id of the chart.
312			*/
313		✗	Chart(Mesh& mesh_in, index_t id_in) :
314		✗	mesh(mesh_in), id(id_in) {
315		✗	}
316
317			/**
318			* \brief Chart copy constructor.
319			* \param[in] rhs a const reference to the Chart to be copied.
320			*/
321		✗	Chart(const Chart& rhs) :
322		✗	mesh(rhs.mesh), facets(rhs.facets), id(rhs.id) {
323		✗	}
324
325			/**
326			* \brief Chart affectation.
327			* \param[in] rhs a const reference to the Chart to be copied.
328			* \return a reference to this Chart after copy.
329			* \pre rhs.mesh == mesh
330			*/
331			Chart& operator=(const Chart& rhs) {
332			if(&rhs != this) {
333			geo_debug_assert(&rhs.mesh == &mesh);
334			facets = rhs.facets;
335			id = rhs.id;
336			}
337			return *this;
338			}
339
340			/**
341			* \brief Gets the number of edges on the border of this chart.
342			* \details An edge is on the border of a chart if it is on the
343			* border of the surface or if the adjacent facet is on a different
344			* chart.
345			*/
346			index_t nb_edges_on_border() const {
347			index_t result = 0;
348			Attribute<index_t> chart(mesh.facets.attributes(),"chart");
349			for(index_t f1: facets) {
350			for(index_t le = 0; le < mesh.facets.nb_vertices(f1); ++le) {
351			index_t f2 = mesh.facets.adjacent(f1,le);
352			if(f2 == NO_INDEX || chart[f2] != id) {
353			++result;
354			}
355			}
356			}
357			return result;
358			}
359
360			/**
361			* \brief Tests whether a chart is shaped like a sock.
362			* \details A chart is shaped like a sock if the area of
363			* the holes is smaller than a certain threshold relative
364			* to the surface area.
365			*/
366			bool is_sock(double min_area_ratio = 1.0/6.0) const {
367			Attribute<index_t> chart(mesh.facets.attributes(),"chart");
368
369			vec3 border_bary(0.0, 0.0, 0.0);
370			index_t bary_N = 0;
371
372			for(index_t f1: facets) {
373			for(index_t le = 0; le < mesh.facets.nb_vertices(f1); ++le) {
374			index_t f2 = mesh.facets.adjacent(f1,le);
375			if(f2 == NO_INDEX || chart[f2] != id) {
376			index_t v = mesh.facets.vertex(f1,le);
377			border_bary += mesh.vertices.point(v);
378			++bary_N;
379			}
380			}
381			}
382
383			border_bary *= (1.0 / double(bary_N));
384
385			double total_area = 0.0;
386			double border_area = 0.0;
387
388			for(index_t f1: facets) {
389			total_area += Geom::mesh_facet_area(mesh,f1);
390			index_t N = mesh.facets.nb_vertices(f1);
391			for(index_t le = 0; le < N; ++le) {
392			index_t f2 = mesh.facets.adjacent(f1,le);
393			if(f2 == NO_INDEX || chart[f2] != id) {
394			index_t v1 = mesh.facets.vertex(f1,le);
395			index_t v2 = mesh.facets.vertex(f1,(le+1) % N);
396			vec3 p1 = mesh.vertices.point(v1);
397			vec3 p2 = mesh.vertices.point(v2);
398			border_area += Geom::triangle_area(border_bary, p1, p2);
399			}
400			}
401			}
402			return ((border_area / total_area) < min_area_ratio);
403			}
404
405			/**
406			* \brief A reference to the mesh.
407			*/
408			Mesh& mesh;
409
410			/**
411			* \brief The list of facet indices of this chart.
412			* \details The mesh is supposed to have an Attribute<index_t>
413			* named "chart" attached to the facets, and the list of facets
414			* has all the facets f for which chart[f] == id.
415			*/
416			vector<index_t> facets;
417
418			/**
419			* \brief The id of this chart.
420			* \details The mesh is supposed to have an Attribute<index_t>
421			* named "chart" attached to the facets, and the list of facets
422			* has all the facets f for which chart[f] == id.
423			*/
424			index_t id;
425			};
426
427			namespace Geom {
428
429			/**
430			* \brief Computes the 2D bounding box of a parameterized mesh.
431			* \param[in] mesh a const reference to a parameterized surface mesh.
432			* \param[in] tex_coord the texture coordinates as a 2d vector
433			* attribute attached to the facet corners.
434			* \param[out] xmin , ymin , xmax , ymax references to the
435			* extremum coordinates of the parameterization.
436			*/
437		✗	static void get_mesh_bbox_2d(
438			const Mesh& mesh, Attribute<double>& tex_coord,
439			double& xmin, double& ymin,
440			double& xmax, double& ymax
441			) {
442		✗	xmin = Numeric::max_float64();
443		✗	ymin = Numeric::max_float64();
444		✗	xmax = Numeric::min_float64();
445		✗	ymax = Numeric::min_float64();
446		✗	for(index_t c: mesh.facet_corners) {
447		✗	xmin = std::min(xmin, tex_coord[2*c]);
448		✗	ymin = std::min(ymin, tex_coord[2*c+1]);
449		✗	xmax = std::max(xmax, tex_coord[2*c]);
450		✗	ymax = std::max(ymax, tex_coord[2*c+1]);
451			}
452		✗	}
453
454			/**
455			* \brief Computes the 2D bounding box of a parameterized chart.
456			* \param[in] chart a const reference to a parameterized surface chart.
457			* \param[in] tex_coord the texture coordinates as a
458			* 2d vector attribute attached to the facet corners.
459			* \param[out] xmin , ymin , xmax , ymax references to the
460			* extremum coordinates of the parameterization.
461			*/
462		✗	static void get_chart_bbox_2d(
463			const Chart& chart, Attribute<double>& tex_coord,
464			double& xmin, double& ymin,
465			double& xmax, double& ymax
466			) {
467		✗	xmin = Numeric::max_float64();
468		✗	ymin = Numeric::max_float64();
469		✗	xmax = Numeric::min_float64();
470		✗	ymax = Numeric::min_float64();
471		✗	for(index_t ff=0; ff<chart.facets.size(); ++ff) {
472		✗	index_t f = chart.facets[ff];
473		✗	for(index_t c: chart.mesh.facets.corners(f)) {
474		✗	xmin = std::min(xmin, tex_coord[2*c]);
475		✗	ymin = std::min(ymin, tex_coord[2*c+1]);
476		✗	xmax = std::max(xmax, tex_coord[2*c]);
477		✗	ymax = std::max(ymax, tex_coord[2*c+1]);
478			}
479			}
480		✗	}
481
482			/**
483			* \brief Computes the 2D parameter-space area of a facet in
484			* a parameterized mesh.
485			* \param[in] mesh a const reference to the mesh.
486			* \param[in] f a facet of \p mesh.
487			* \param[in] tex_coord the texture coordinates as a 2d vector
488			* attribute attached to the facet corners.
489			* \return the area of the facet in parameter space.
490			*/
491		✗	static double mesh_facet_area_2d(
492			const Mesh& mesh, index_t f, Attribute<double>& tex_coord
493			) {
494		✗	double result = 0.0;
495		✗	index_t c1 = mesh.facets.corners_begin(f);
496		✗	vec2 p1(tex_coord[2*c1], tex_coord[2*c1+1]);
497		✗	for(
498		✗	index_t c2 = c1+1;
499		✗	c2+1 < mesh.facets.corners_end(f); ++c2
500			) {
501		✗	index_t c3 = c2+1;
502		✗	vec2 p2(tex_coord[2*c2], tex_coord[2*c2+1]);
503		✗	vec2 p3(tex_coord[2*c3], tex_coord[2*c3+1]);
504		✗	result += Geom::triangle_area(p1,p2,p3);
505			}
506		✗	return result;
507			}
508
509			/**
510			* \brief Computes the 2D parameter-space area of a parameterized mesh.
511			* \param[in] mesh a const reference to a parameterized surface mesh.
512			* \param[in] tex_coord the texture coordinates as a
513			* 2d vector attribute attached to the facet corners.
514			* \return the area of the parameter space.
515			*/
516		✗	static double mesh_area_2d(
517			const Mesh& mesh, Attribute<double>& tex_coord
518			) {
519		✗	double result = 0.0;
520		✗	for(index_t f: mesh.facets) {
521		✗	result += mesh_facet_area_2d(mesh, f, tex_coord);
522			}
523		✗	return result;
524			}
525
526
527			/**
528			* \brief Computes the 2D parameter-space area of a parameterized chart.
529			* \param[in] chart a const reference to a parameterized surface chart.
530			* \param[in] tex_coord the texture coordinates as a 2d vector
531			* attribute attached to the facet corners.
532			* \return the area of the parameter space.
533			*/
534		✗	static double chart_area_2d(
535			const Chart& chart, Attribute<double>& tex_coord
536			) {
537		✗	double result = 0.0;
538		✗	for(index_t ff=0; ff<chart.facets.size(); ++ff) {
539		✗	index_t f = chart.facets[ff];
540		✗	result += mesh_facet_area_2d(chart.mesh, f, tex_coord);
541			}
542		✗	return result;
543			}
544
545			/**
546			* \brief Computes the area of a chart.
547			* \param[in] chart a const reference to a chart.
548			* \return the area of the chart in 3D.
549			*/
550		✗	static double chart_area(const Chart& chart) {
551		✗	double result = 0.0;
552		✗	for(index_t ff=0; ff<chart.facets.size(); ++ff) {
553		✗	index_t f = chart.facets[ff];
554		✗	result += Geom::mesh_facet_area(chart.mesh, f);
555			}
556		✗	return result;
557			}
558
559
560			}
561			}
562
563			namespace {
564			using namespace GEO;
565
566			class AverageDirection2d {
567			public:
568
569		✗	AverageDirection2d() {
570		✗	result_is_valid_ = false;
571		✗	}
572
573		✗	void begin() {
574		✗	for(int i=0; i<3; i++) {
575		✗	M_[i] = 0.0 ;
576			}
577		✗	result_is_valid_ = false;
578		✗	}
579
580		✗	void add_vector(const vec2& e) {
581		✗	M_[0] += e.x * e.x ;
582		✗	M_[1] += e.x * e.y ;
583		✗	M_[2] += e.y * e.y ;
584		✗	}
585
586		✗	void end() {
587			double eigen_vectors[4] ;
588			double eigen_values[2] ;
589		✗	MatrixUtil::semi_definite_symmetric_eigen(
590		✗	M_, 2, eigen_vectors, eigen_values
591			);
592		✗	int k = 0 ;
593		✗	result_ = vec2(
594		✗	eigen_vectors[2*k],
595		✗	eigen_vectors[2*k+1]
596		✗	);
597		✗	result_is_valid_ = true ;
598		✗	}
599
600		✗	const vec2& average_direction() const {
601		✗	geo_assert(result_is_valid_) ;
602		✗	return result_ ;
603			}
604
605			private:
606			double M_[3] ;
607			vec2 result_ ;
608			bool result_is_valid_ ;
609			} ;
610			}
611
612
613			namespace GEO {
614
615			/**
616			* \brief Stores a reference to a chart, its bounding box
617			* in parameter space, and its upper and lower horizons.
618			*/
619			class ChartBBox {
620			public:
621
622			/**
623			* \brief ChartBBox constructor.
624			* \param[in] chart a pointer to the chart.
625			* \param[in] min , max lower-left and upper-right
626			* corners of the parameter-space bounding box of
627			* the chart.
628			*/
629		✗	ChartBBox(
630			Chart* chart, const vec2& min, const vec2& max
631		✗	) :
632		✗	min_(min), max_(max), chart_(chart),
633		✗	min_func_(nullptr), max_func_(nullptr),
634		✗	nb_steps_(0) {
635		✗	}
636
637			/**
638			* \brief ChartBBox destructor.
639			*/
640		✗	~ChartBBox(){
641		✗	free();
642		✗	chart_ = nullptr;
643		✗	}
644
645			/**
646			* \brief ChartBBox copy-constructor.
647			* \param[in] rhs a const reference to the ChartBBox to be copied.
648			*/
649		✗	ChartBBox(const ChartBBox& rhs) {
650		✗	copy(rhs);
651		✗	}
652
653			/**
654			* \brief ChartBBox affectation.
655			* \param[in] rhs a const reference to the ChartBBox to be copied.
656			* \return A reference to this ChartBBox after affectation.
657			*/
658		✗	ChartBBox& operator=(const ChartBBox& rhs) {
659		✗	if(&rhs != this) {
660		✗	free();
661		✗	copy(rhs);
662			}
663		✗	return *this;
664			}
665
666			/**
667			* \brief Initializes the upper and lower horizons.
668			* \param[in] step pixel-size in parameter space
669			* \param[in] margin margin size in parameter space
670			* \param[in] margin_width_in_pixels margin size in pixels
671			* \param[in] tex_coord a 2d vector attribute attached to
672			* the facet corners of the mesh with the texture coordinates.
673			* \param[in] chart_attr an attribute attached to the facets of
674			* the mesh with the id of the chart the facet belongs to.
675			*/
676		✗	void init_max_and_min_func(
677			double step, double margin, index_t margin_width_in_pixels,
678			Attribute<double>& tex_coord, Attribute<index_t>& chart_attr
679			) {
680
681		✗	if(min_func_ != nullptr) {
682		✗	delete[] min_func_;
683			}
684		✗	if(max_func_ != nullptr) {
685		✗	delete[] max_func_;
686			}
687
688		✗	nb_steps_ = int(1.0 + (max_.x - min_.x) / step );
689		✗	min_func_ = new double[size_t(nb_steps_)];
690		✗	max_func_ = new double[size_t(nb_steps_)];
691		✗	for (int i=0;i<nb_steps_;i++){
692		✗	min_func(i) = max_.y - min_.y;
693		✗	max_func(i) = 0;
694			}
695
696		✗	for(index_t ff=0; ff<chart_->facets.size(); ++ff) {
697		✗	index_t f = chart_->facets[ff];
698		✗	for(
699		✗	index_t c1: chart_->mesh.facets.corners(f)) {
700			index_t neighf =
701		✗	chart_->mesh.facet_corners.adjacent_facet(c1);
702
703		✗	if(neighf != NO_FACET && chart_attr[neighf] == chart_->id) {
704		✗	continue;
705			}
706
707			index_t c2 =
708		✗	chart_->mesh.facets.next_corner_around_facet(f,c1);
709		✗	vec2 p1(tex_coord[2*c1], tex_coord[2*c1+1]);
710		✗	vec2 p2(tex_coord[2*c2], tex_coord[2*c2+1]);
711
712		✗	if(p2.x == p1.x) {
713		✗	continue;
714			}
715
716		✗	if(p2.x < p1.x) {
717		✗	std::swap(p1,p2);
718			}
719
720		✗	p1.x -= double(margin_width_in_pixels) * step;
721		✗	p2.x += double(margin_width_in_pixels) * step;
722
723		✗	double a = (p2.y - p1.y) / (p2.x - p1.x);
724		✗	for(double x = p1.x; x<p2.x; x+=step) {
725		✗	double y = p1.y + a * (x - p1.x);
726		✗	int cX = int((x - min_.x) / step);
727		✗	if (cX>=0 && cX<nb_steps_) {
728		✗	min_func(cX) = std::min(min_func(cX), y - margin);
729		✗	max_func(cX) = std::max(max_func(cX), y + margin);
730			}
731			}
732			}
733			}
734		✗	}
735
736		✗	double& max_func(int i) {
737		✗	geo_assert(i>=0 && i<nb_steps_);
738		✗	return max_func_[i];
739			}
740
741		✗	double& min_func(int i) {
742		✗	geo_assert(i>=0 && i<nb_steps_);
743		✗	return min_func_[i];
744			}
745
746		✗	vec2 size() const {
747		✗	return max_-min_;
748			}
749
750			/**
751			* \brief Gets the area of the bounding box.
752			* \return the area of the bounding box in parameter space.
753			*/
754		✗	double area() const {
755		✗	vec2 s = size();
756		✗	return s.x*s.y;
757			}
758
759			/**
760			* \brief Applies a translation vector to the texture coordinates.
761			* \param[in] v the translation vector.
762			* \param[in,out] tex_coord a dimension 2 vector attribute attached to
763			* the facet corners of the mesh with the texture coordinates.
764			*/
765		✗	void translate(const vec2& v, Attribute<double>& tex_coord){
766		✗	min_=min_+v;
767		✗	max_=max_+v;
768		✗	for(index_t ff=0; ff<chart_->facets.size(); ++ff) {
769		✗	index_t f = chart_->facets[ff];
770		✗	for(index_t c: chart_->mesh.facets.corners(f)) {
771		✗	tex_coord[2*c] += v.x;
772		✗	tex_coord[2*c+1] += v.y;
773			}
774			}
775		✗	}
776
777			/**
778			* \brief Gets the lower-left corner of the bounding box.
779			* \return A const reference to the parametric-space
780			* coordinates of the lower-left corner of the bounding box.
781			*/
782			const vec2& min() const {
783			return min_;
784			}
785
786			/**
787			* \brief Gets the upper-right corner of the bounding box.
788			* \return A const reference to the parametric-space
789			* coordinates of the upper-right corner of the bounding box.
790			*/
791			const vec2& max() const {
792			return max_;
793			}
794
795			/**
796			* \brief Gets the lower-left corner of the bounding box.
797			* \return A modifiable reference to the parametric-space
798			* coordinates of the lower-left corner of the bounding box.
799			*/
800		✗	vec2& min() {
801		✗	return min_;
802			}
803
804			/**
805			* \brief Gets the upper-right corner of the bounding box.
806			* \return A modifiable reference to the parametric-space
807			* coordinates of the upper-right corner of the bounding box.
808			*/
809		✗	vec2& max() {
810		✗	return max_;
811			}
812
813			/**
814			* \brief Gets the chart.
815			* \return a const pointer to the chart.
816			*/
817			const Chart* chart() const {
818			return chart_;
819			}
820
821
822			protected:
823
824			/**
825			* \brief Copies a ChartBBox.
826			* \param[in] rhs a const reference to the ChartBBox
827			* to be copied.
828			*/
829		✗	void copy(const ChartBBox& rhs) {
830		✗	chart_ = rhs.chart_;
831		✗	nb_steps_ = rhs.nb_steps_;
832
833		✗	if(rhs.min_func_ != nullptr) {
834		✗	min_func_ = new double[size_t(nb_steps_)];
835		✗	max_func_ = new double[size_t(nb_steps_)];
836		✗	for(int i=0; i<nb_steps_; i++) {
837		✗	min_func_[i] = rhs.min_func_[i];
838		✗	max_func_[i] = rhs.max_func_[i];
839			}
840			} else {
841		✗	geo_assert(rhs.max_func_ == nullptr);
842		✗	min_func_ = nullptr;
843		✗	max_func_ = nullptr;
844			}
845		✗	min_ = rhs.min_;
846		✗	max_ = rhs.max_;
847		✗	}
848
849			/**
850			* \brief Releases the memory allocated by
851			* this ChartBBox.
852			*/
853		✗	void free() {
854		✗	delete[] min_func_;
855		✗	delete[] max_func_;
856		✗	min_func_ = nullptr;
857		✗	max_func_ = nullptr;
858		✗	nb_steps_ = 0;
859		✗	}
860
861			private:
862			vec2 min_, max_;
863			Chart* chart_;
864			double* min_func_;
865			double* max_func_;
866			int nb_steps_;
867			};
868
869			/***********************************************************/
870
871			/**
872			* \brief Internal implementation of the packing algorithm.
873			* \details By Nicolas Ray. The algorithm is inspired by
874			* the Tetris game, and iteratively places the charts from
875			* bottom to top.
876			*/
877			class TetrisPacker {
878			public :
879
880			/**
881			* \brief TetrisPacker constructor.
882			* \param[in] mesh a reference the surface mesh to be packed.
883			* It needs to have a vector attribute of dimension 2 attached
884			* to the facet corners and named "tex_coord", as well as a
885			* facet attribute named "chart".
886			*/
887		✗	TetrisPacker(Mesh& mesh) {
888		✗	nb_xpos_ = 1024;
889		✗	height_ = new double[size_t(nb_xpos_)];
890		✗	image_size_in_pixels_ = 1024;
891		✗	margin_width_in_pixels_ = 4;
892		✗	tex_coord_.bind_if_is_defined(
893		✗	mesh.facet_corners.attributes(), "tex_coord"
894			);
895		✗	geo_assert(tex_coord_.is_bound() && tex_coord_.dimension() == 2);
896		✗	chart_attr_.bind_if_is_defined(
897		✗	mesh.facets.attributes(), "chart"
898			);
899		✗	if(!chart_attr_.is_bound()) {
900		✗	mesh_get_charts(mesh);
901		✗	chart_attr_.bind_if_is_defined(
902		✗	mesh.facets.attributes(), "chart"
903			);
904		✗	geo_assert(chart_attr_.is_bound());
905			}
906		✗	}
907
908			/**
909			* \brief TetrisPacker destructor.
910			*/
911		✗	~TetrisPacker() {
912		✗	delete[] height_;
913		✗	height_ = nullptr;
914		✗	}
915
916		✗	void set_image_size_in_pixels(index_t size) {
917		✗	image_size_in_pixels_ = size;
918		✗	}
919
920			index_t margin_width_in_pixels() const {
921			return margin_width_in_pixels_;
922			}
923
924		✗	void set_margin_width_in_pixels(index_t width) {
925		✗	margin_width_in_pixels_ = width;
926		✗	}
927
928		✗	void add(const ChartBBox& r){
929		✗	data_.push_back(r);
930		✗	}
931
932			/**
933			* \brief Compares two ChartBBox objects.
934			* \details Used by the packing algorithm to sort
935			* the boxes.
936			* \param[in] b0 , b1 const references to the
937			* two ChartBBox objects to be compared.
938			* \retval true if \p b0 is taller than \p b1.
939			* \retval false otherwise.
940			*/
941		✗	static bool compare(
942			const ChartBBox& b0, const ChartBBox& b1
943			) {
944		✗	return (b0.size().y > b1.size().y);
945			}
946
947		✗	void add_margin(double margin_size) {
948		✗	for (unsigned int i=0;i<data_.size();i++){
949		✗	data_[i].translate(
950		✗	vec2(
951		✗	-data_[i].min().x + margin_size,
952		✗	-data_[i].min().y + margin_size
953		✗	),
954		✗	tex_coord_
955			);
956		✗	data_[i].max() = data_[i].max() +
957		✗	vec2( 2.0 * margin_size, 2.0 * margin_size);
958		✗	geo_assert(data_[i].max().x>0);
959		✗	geo_assert(data_[i].max().y>0);
960		✗	data_[i].min() = vec2(0,0);
961			}
962		✗	}
963
964		✗	double max_height() {
965		✗	double result = 0;
966		✗	for (unsigned int i=0;i<data_.size();i++) {
967		✗	result = std::max(result, data_[i].max().y);
968			}
969		✗	return result;
970			}
971
972		✗	void recursive_apply() {
973			// compute margin
974		✗	double area = 0;
975		✗	for (unsigned int numrect = 0; numrect <data_.size();numrect++ ) {
976		✗	area += data_[numrect].area();
977			}
978			double margin =
979		✗	(::sqrt(area) / double(image_size_in_pixels_)) *
980		✗	double(margin_width_in_pixels_);
981		✗	add_margin(margin);
982
983			// find a first solution
984		✗	apply(margin);
985		✗	double scoreSup = max_height();
986		✗	double borneSup = width_;
987		✗	double decalborne = 0.5 * ::sqrt(scoreSup * width_);
988
989			// dichotomy
990		✗	for (int i=0;i<10;i++ ){
991		✗	double new_borne = borneSup - decalborne;
992		✗	apply(margin,new_borne);
993		✗	double max = max_height();
994		✗	if (max < new_borne){
995		✗	borneSup = borneSup - decalborne;
996			}
997		✗	decalborne/=2;
998			}
999		✗	apply(margin, borneSup);
1000		✗	}
1001
1002
1003		✗	void apply(double margin, double width =-1) {
1004		✗	width_ = width;
1005		✗	for (unsigned int i=0; i<data_.size(); i++) {
1006		✗	data_[i].translate(
1007		✗	vec2(
1008		✗	-data_[i].min().x,
1009		✗	-data_[i].min().y
1010		✗	),
1011		✗	tex_coord_
1012			);
1013			}
1014
1015			// sort ChartBBoxes by their heights ...
1016		✗	std::sort(data_.begin(),data_.end(),compare);
1017
1018			{ //find the best width
1019		✗	double max_bbox_width = 0;
1020		✗	for (unsigned int i=0; i<data_.size(); i++) {
1021		✗	max_bbox_width= std::max(
1022		✗	max_bbox_width, data_[i].size().x
1023			);
1024			}
1025
1026		✗	if ( width == -1){
1027			// try to have a square : width_ = sqrt(area);
1028		✗	double area =0;
1029
1030		✗	for (unsigned int i=0; i<data_.size(); i++){
1031		✗	area += data_[i].area();
1032			}
1033
1034		✗	width_ = ::sqrt(area) *1.1;
1035
1036
1037			// resize if a square seems to be too bad
1038			// (the first piece is too high)
1039
1040		✗	if (data_[0].size().y > width_) {
1041		✗	width_ = area / data_[0].size().y;
1042			}
1043			}
1044
1045			// be sure all surface can fit in the width
1046		✗	width_ = std::max(
1047		✗	max_bbox_width * (double(nb_xpos_ + 2) / double(nb_xpos_)),
1048		✗	width_
1049			);
1050			}
1051			// set the step depending on the width and the discretisation
1052		✗	step_ = width_ / nb_xpos_;
1053
1054
1055			// init local min and max height functions
1056		✗	for (unsigned int numrect = 0; numrect <data_.size(); numrect++) {
1057		✗	data_[numrect].init_max_and_min_func(
1058			step_, margin, margin_width_in_pixels_,
1059		✗	tex_coord_, chart_attr_
1060			);
1061			}
1062
1063			// init global height function
1064		✗	for (int x=0; x<nb_xpos_; x++) {
1065		✗	height(x)=0;
1066			}
1067
1068		✗	for (unsigned int i=0;i<data_.size();i++) {
1069		✗	place(data_[i]);
1070			}
1071		✗	}
1072
1073
1074			private:
1075
1076		✗	double& height(int x) {
1077		✗	geo_assert(x >= 0 && x < nb_xpos_);
1078		✗	return height_[x];
1079			}
1080
1081			const double& height(int x) const {
1082			geo_assert(x >= 0 && x < nb_xpos_);
1083			return height_[x];
1084			}
1085
1086
1087			/**
1088			* \brief Inserts a new chart.
1089			* \details Follows the "tetris" strategy.
1090			* \param[in] rect the chart to be inserted.
1091			*/
1092		✗	void place(ChartBBox& rect) {
1093
1094		✗	const int width_in_pas = int ( (rect.size().x / step_) + 1 );
1095
1096			// find the best position
1097		✗	int bestXPos=0;
1098		✗	double bestHeight = Numeric::max_float64();
1099		✗	double bestFreeArea = Numeric::max_float64();
1100
1101		✗	for (int x=0;x<nb_xpos_ - width_in_pas;x++) {
1102
1103		✗	double localHeight = max_height(x,width_in_pas,rect);
1104			double localFreeArea =
1105		✗	free_area(x,width_in_pas,localHeight)
1106		✗	+ localHeight * rect.size().x;
1107
1108			// the best position criterion is the position that
1109			// minimize area lost
1110			// area could be lost under and upper the bbox
1111		✗	if (localFreeArea < bestFreeArea){
1112		✗	bestXPos = x;
1113		✗	bestHeight = localHeight;
1114		✗	bestFreeArea = localFreeArea;
1115			}
1116			}
1117
1118			// be sure we have a solution
1119		✗	geo_assert(bestHeight != Numeric::max_float64());
1120
1121			// place the rectangle
1122		✗	rect.translate(vec2(bestXPos*step_,bestHeight),tex_coord_);
1123		✗	for (int i=bestXPos;i<bestXPos + width_in_pas;i++){
1124		✗	height(i) = rect.min().y + rect.max_func(i-bestXPos);
1125			}
1126		✗	}
1127
1128			// find the max height in the range [xpos, xpos + width]
1129		✗	double max_height(int xpos, int width, ChartBBox& rect){
1130		✗	double result = 0;
1131		✗	for (int i=xpos; i<xpos+width; i++) {
1132		✗	result = std::max( result, height(i)-rect.min_func(i-xpos) );
1133			}
1134		✗	return result;
1135			}
1136
1137			// find the free area in the range under height_max in range
1138			// [xpos, xpos + width]
1139
1140		✗	double free_area(int xpos, int width, double height_max) {
1141		✗	double result =0;
1142		✗	for (int i=xpos; i<xpos+width; i++) {
1143		✗	result += step_ * (height_max - height(i));
1144			}
1145		✗	return result;
1146			}
1147
1148			private:
1149			index_t image_size_in_pixels_;
1150			index_t margin_width_in_pixels_;
1151
1152			int nb_xpos_;
1153			double width_;
1154			double step_;
1155			double *height_;
1156
1157			std::vector<ChartBBox> data_;
1158
1159			Attribute<double> tex_coord_;
1160			Attribute<index_t> chart_attr_;
1161			};
1162
1163			}
1164
1165			/***************************************************************
1166			**** Tetris Packer ****
1167			***************************************************************/
1168
1169			namespace {
1170			using namespace GEO;
1171
1172			// There where some problems with nan (Not a number),
1173			// this code is used to track such problems (seems to
1174			// be ok now)
1175		✗	static bool chart_is_ok(Chart& chart, Attribute<double>& tex_coord) {
1176		✗	for(index_t ff=0; ff<chart.facets.size(); ++ff) {
1177		✗	index_t f = chart.facets[ff];
1178		✗	for(index_t c: chart.mesh.facets.corners(f)) {
1179		✗	if(Numeric::is_nan(tex_coord[2*c])) {
1180		✗	return false;
1181			}
1182		✗	if(Numeric::is_nan(tex_coord[2*c+1])) {
1183		✗	return false;
1184			}
1185			}
1186			}
1187		✗	return true;
1188			}
1189
1190			/**
1191			* \brief Packs a set of surfaces in parameter (texture) space.
1192			* \details This organizes the charts of a parameterization in
1193			* a way that tries to minimize unused texture space.
1194			* The algorithm, by Nicolas Ray, is described in the following
1195			* reference:
1196			* - least squares mode: Least Squares Conformal Maps,
1197			* Levy, Petitjean, Ray, Maillot, ACM SIGGRAPH, 2002
1198			*/
1199			class GEOGRAM_API Packer {
1200			public:
1201
1202			/**
1203			* \brief Packer constructor.
1204			*/
1205		✗	Packer() {
1206		✗	image_size_in_pixels_ = 1024;
1207		✗	margin_width_in_pixels_ = 4;
1208		✗	}
1209
1210			/**
1211			* \brief Forbids copy.
1212			*/
1213			Packer(const Packer& rhs) = delete;
1214
1215			/**
1216			* \brief Forbids copy.
1217			*/
1218			Packer& operator=(const Packer& rhs) = delete;
1219
1220			/**
1221			* \brief Packs a texture atlas.
1222			* \param[in,out] mesh a surface mesh.
1223			* \details The mesh is supposed to have texture coordinates
1224			* stored in a 2d vector attribute attached to the facet
1225			* corners and named "tex_coord".
1226			*/
1227		✗	void pack_surface(Mesh& mesh, bool normalize_only=false) {
1228		✗	tex_coord_.bind_if_is_defined(
1229		✗	mesh.facet_corners.attributes(), "tex_coord"
1230			);
1231		✗	geo_assert(tex_coord_.is_bound() && tex_coord_.dimension() == 2);
1232		✗	chart_attr_.bind_if_is_defined(
1233		✗	mesh.facets.attributes(), "chart"
1234			);
1235		✗	if(!chart_attr_.is_bound()) {
1236		✗	mesh_get_charts(mesh);
1237		✗	chart_attr_.bind_if_is_defined(
1238		✗	mesh.facets.attributes(), "chart"
1239			);
1240		✗	geo_assert(chart_attr_.is_bound());
1241			}
1242
1243			// Get the charts
1244		✗	vector<Chart> charts;
1245		✗	index_t nb_charts=0;
1246		✗	for(index_t f: mesh.facets) {
1247		✗	nb_charts = std::max(nb_charts, chart_attr_[f]);
1248			}
1249		✗	++nb_charts;
1250
1251		✗	for(index_t i=0; i<nb_charts; ++i) {
1252		✗	charts.push_back(Chart(mesh, i));
1253			}
1254
1255		✗	for(index_t f: mesh.facets) {
1256		✗	charts[chart_attr_[f]].facets.push_back(f);
1257			}
1258
1259		✗	Logger::out("Packer") << "Packing "
1260		✗	<< charts.size() << " charts" << std::endl;
1261
1262			// Sanity check
1263		✗	for(index_t i=0; i<nb_charts; ++i) {
1264		✗	if(!chart_is_ok(charts[i], tex_coord_)) {
1265		✗	for(index_t ff=0; ff<charts[i].facets.size(); ++ff) {
1266		✗	index_t f = charts[i].facets[ff];
1267		✗	for(index_t c: mesh.facets.corners(f)) {
1268		✗	tex_coord_[2*c] = 0.0;
1269		✗	tex_coord_[2*c+1] = 0.0;
1270			}
1271			}
1272			}
1273			}
1274
1275		✗	pack_charts(charts, normalize_only);
1276
1277			// Normalize tex coords in [0,1] x [0,1]
1278			{
1279			double u_min, v_min, u_max, v_max;
1280		✗	Geom::get_mesh_bbox_2d(
1281		✗	mesh, tex_coord_, u_min, v_min, u_max, v_max
1282			);
1283		✗	double l = std::max(u_max - u_min, v_max - v_min);
1284		✗	if(l > 1e-6) {
1285		✗	for(index_t c: mesh.facet_corners) {
1286		✗	tex_coord_[2*c] = (tex_coord_[2*c] - u_min) / l;
1287		✗	tex_coord_[2*c+1] = (tex_coord_[2*c+1] - v_min) / l;
1288			}
1289			}
1290			}
1291		✗	tex_coord_.unbind();
1292		✗	chart_attr_.unbind();
1293		✗	}
1294
1295			/**
1296			* \brief Gets the size of the target texture image in
1297			* pixels.
1298			* \return the size of the target texture image.
1299			*/
1300		✗	index_t image_size_in_pixels() const {
1301		✗	return image_size_in_pixels_;
1302			}
1303
1304			/**
1305			* \brief Gets the size of the margin (or "gutter") around the charts.
1306			* \details This may be required to avoid undesirable blends due to
1307			* mip-mapping.
1308			* \return the number of empty pixels to be preserved around each chart.
1309			*/
1310		✗	index_t margin_width_in_pixels() const {
1311		✗	return margin_width_in_pixels_;
1312			}
1313
1314			protected:
1315			/**
1316			* \brief Packs a set of charts.
1317			* \param[in,out] charts a const reference to a vector with the
1318			* charts to be packed.
1319			* \param[in] normalize_only if set, just normalize texture coordinates
1320			* and do not pack the charts.
1321			* \details All the charts are supposed to be attached to the same mesh.
1322			* Texture coordinates are stored in a 2d vector attribute attached to
1323			* the facet corners of the mesh and called "tex_coord".
1324			*/
1325		✗	void pack_charts(vector<Chart>& charts, bool normalize_only = false) {
1326
1327		✗	if(charts.size() == 0) {
1328		✗	return;
1329			}
1330
1331		✗	Mesh& mesh = charts[0].mesh;
1332
1333		✗	for(index_t i=0; i<charts.size(); ++i) {
1334		✗	geo_assert(chart_is_ok(charts[i], tex_coord_));
1335		✗	normalize_chart(charts[i]);
1336		✗	geo_assert(chart_is_ok(charts[i], tex_coord_));
1337			}
1338
1339		✗	if(normalize_only) {
1340		✗	return;
1341			}
1342
1343			// use the tetris packer (more efficient for large dataset)
1344			// set some application dependent const
1345		✗	TetrisPacker pack(mesh);
1346		✗	pack.set_image_size_in_pixels(image_size_in_pixels());
1347		✗	pack.set_margin_width_in_pixels(margin_width_in_pixels());
1348
1349		✗	for(index_t i=0; i<charts.size(); ++i) {
1350		✗	geo_assert(chart_is_ok(charts[i], tex_coord_));
1351			double u_min, v_min, u_max, v_max;
1352		✗	Geom::get_chart_bbox_2d(
1353		✗	charts[i], tex_coord_, u_min, v_min, u_max, v_max
1354			);
1355
1356		✗	geo_assert(!Numeric::is_nan(u_min));
1357		✗	geo_assert(!Numeric::is_nan(v_min));
1358		✗	geo_assert(!Numeric::is_nan(u_max));
1359		✗	geo_assert(!Numeric::is_nan(v_max));
1360		✗	geo_assert(u_max >= u_min);
1361		✗	geo_assert(v_max >= v_min);
1362
1363			ChartBBox r(
1364		✗	&charts[i],
1365		✗	vec2(u_min,v_min) ,
1366		✗	vec2(u_max,v_max)
1367		✗	);
1368		✗	pack.add(r);
1369		✗	}
1370
1371
1372		✗	pack.recursive_apply();
1373
1374			{
1375		✗	double total_area = Geom::mesh_area_2d(mesh, tex_coord_);
1376			double u_min, v_min, u_max, v_max;
1377		✗	Geom::get_mesh_bbox_2d(
1378		✗	mesh, tex_coord_, u_min, v_min, u_max, v_max
1379			);
1380		✗	double bbox_area = (u_max - u_min)*(v_max-v_min);
1381		✗	double filling_ratio = total_area / bbox_area;
1382		✗	Logger::out("Packer") << "BBox area:" << bbox_area << std::endl;
1383		✗	Logger::out("Packer") << "Filling ratio:"
1384		✗	<< filling_ratio << std::endl;
1385			}
1386		✗	}
1387
1388			/**
1389			* \brief Normalizes the parameterization of a chart.
1390			* \param[in,out] chart a reference to the chart to be normalized.
1391			* \details This rescales texture coordinates in such a way that the
1392			* chart has the same area in 3D and in texture space.
1393			* This also applies a rotation to the texture coordinates such that
1394			* the area of the bounding rectangle is minimized.
1395			* Texture coordinates are stored in a 2d vector attribute attached
1396			* to the facet corners of the mesh and called "tex_coord".
1397			*/
1398		✗	void normalize_chart(Chart& chart) {
1399			// TODO: choose best axes, rotate chart.
1400		✗	vec2 U(1.0, 0.0);
1401		✗	vec2 V(0.0, 1.0);
1402		✗	AverageDirection2d dir;
1403		✗	dir.begin();
1404		✗	for(index_t ff=0; ff<chart.facets.size(); ++ff) {
1405		✗	index_t f = chart.facets[ff];
1406		✗	for(index_t c1: chart.mesh.facets.corners(f)) {
1407		✗	index_t adj_f = chart.mesh.facet_corners.adjacent_facet(c1);
1408		✗	if(adj_f == NO_FACET || chart_attr_[adj_f] != chart.id) {
1409			index_t c2 =
1410		✗	chart.mesh.facets.next_corner_around_facet(f,c1);
1411		✗	vec2 uv1(tex_coord_[2*c1], tex_coord_[2*c1+1]);
1412		✗	vec2 uv2(tex_coord_[2*c2], tex_coord_[2*c2+1]);
1413		✗	dir.add_vector(uv2-uv1);
1414			}
1415			}
1416			}
1417		✗	dir.end();
1418		✗	vec2 W=dir.average_direction();
1419		✗	if(!Geom::has_nan(W)) {
1420		✗	W = normalize(W);
1421		✗	V = W;
1422		✗	U = vec2(-V.y, V.x);
1423			}
1424
1425		✗	for(index_t ff=0; ff<chart.facets.size(); ++ff) {
1426		✗	index_t f = chart.facets[ff];
1427		✗	for(index_t c: chart.mesh.facets.corners(f)) {
1428		✗	vec2 uv(tex_coord_[2*c], tex_coord_[2*c+1]);
1429		✗	tex_coord_[2*c] = dot(uv,U);
1430		✗	tex_coord_[2*c+1] = dot(uv,V);
1431			}
1432			}
1433
1434		✗	double area3d = Geom::chart_area(chart);
1435		✗	double area2d = Geom::chart_area_2d(chart, tex_coord_);
1436		✗	double factor = 1.0;
1437		✗	if(::fabs(area2d) > 1e-30) {
1438		✗	factor = ::sqrt(area3d/area2d);
1439			} else {
1440		✗	factor = 0.0;
1441			}
1442
1443		✗	geo_assert(!Numeric::is_nan(area2d));
1444		✗	geo_assert(!Numeric::is_nan(area3d));
1445		✗	geo_assert(!Numeric::is_nan(factor));
1446
1447		✗	for(index_t ff=0; ff<chart.facets.size(); ++ff) {
1448		✗	index_t f = chart.facets[ff];
1449		✗	for(index_t c: chart.mesh.facets.corners(f)) {
1450		✗	tex_coord_[2*c] *= factor;
1451		✗	tex_coord_[2*c+1] *= factor;
1452			}
1453			}
1454		✗	}
1455
1456			private:
1457			index_t image_size_in_pixels_ ;
1458			index_t margin_width_in_pixels_ ;
1459			Attribute<double> tex_coord_;
1460			Attribute<index_t> chart_attr_;
1461			} ;
1462			}
1463
1464			namespace GEO {
1465
1466		✗	void pack_atlas_only_normalize_charts(Mesh& mesh) {
1467		✗	Packer tetris;
1468		✗	bool normalize_tex_coords_only = true;
1469		✗	tetris.pack_surface(mesh, normalize_tex_coords_only);
1470		✗	}
1471
1472		✗	void pack_atlas_using_tetris_packer(Mesh& mesh) {
1473		✗	Packer tetris;
1474		✗	tetris.pack_surface(mesh);
1475		✗	}
1476			}
1477