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/NL/nl.h>
41
42			#include <algorithm>
43			#include <vector>
44			#include <set>
45			#include <string>
46			#include <iostream>
47			#include <sstream>
48			#include <fstream>
49			#include <cstdlib>
50			#include <cstring>
51			#include <cmath>
52			#include <cassert>
53
54			#ifdef __EMSCRIPTEN__
55			#include <emscripten.h>
56			#endif
57
58			/******************************************************************************/
59			/* Basic geometric types */
60
61			/**
62			* \brief A 2D vector
63			*/
64			class vec2 {
65			public:
66			/**
67			* \brief Constructs a 2D vector from its coordinates.
68			* \param x_in , y_in the coordinates.
69			*/
70		✗	vec2(double x_in, double y_in) :
71		✗	x(x_in), y(y_in) {
72		✗	}
73
74			/**
75			* \brief Constructs the zero vector.
76			*/
77		✗	vec2() : x(0), y(0) {
78		✗	}
79
80			double x;
81			double y;
82			};
83
84			/**
85			* \brief A 3D vector
86			*/
87			class vec3 {
88			public:
89
90			/**
91			* \brief Constructs a 3D vector from its coordinates.
92			* \param x_in , y_in , z_in the coordinates.
93			*/
94		✗	vec3(double x_in, double y_in, double z_in) :
95		✗	x(x_in), y(y_in), z(z_in) {
96		✗	}
97
98			/**
99			* \brief Constructs the zero vector.
100			*/
101		✗	vec3() : x(0), y(0), z(0) {
102		✗	}
103
104			/**
105			* \brief Gets the length of this vector.
106			* \return the length of this vector.
107			*/
108		✗	double length() const {
109		✗	return sqrt(x*x + y*y + z*z);
110			}
111
112			/**
113			* \brief Normalizes this vector.
114			* \details This makes the norm equal to 1.0
115			*/
116		✗	void normalize() {
117		✗	double l = length();
118		✗	x /= l; y /= l; z /= l;
119		✗	}
120
121			double x;
122			double y;
123			double z;
124			};
125
126			/**
127			* \brief Outputs a 2D vector to a stream.
128			* \param[out] out a reference to the stream
129			* \param[in] v a const reference to the vector
130			* \return the new state of the stream
131			* \relates vec2
132			*/
133		✗	inline std::ostream& operator<<(std::ostream& out, const vec2& v) {
134		✗	return out << v.x << " " << v.y;
135			}
136
137			/**
138			* \brief Outputs a 3D vector to a stream.
139			* \param[out] out a reference to the stream
140			* \param[in] v a const reference to the vector
141			* \return the new state of the stream
142			* \relates vec3
143			*/
144		✗	inline std::ostream& operator<<(std::ostream& out, const vec3& v) {
145		✗	return out << v.x << " " << v.y << " " << v.z;
146			}
147
148			/**
149			* \brief Reads a 2D vector from a stream
150			* \param[in] in a reference to the stream
151			* \param[out] v a reference to the vector
152			* \return the new state of the stream
153			* \relates vec2
154			*/
155			inline std::istream& operator>>(std::istream& in, vec2& v) {
156			return in >> v.x >> v.y;
157			}
158
159			/**
160			* \brief Reads a 3D vector from a stream
161			* \param[in] in a reference to the stream
162			* \param[out] v a reference to the vector
163			* \return the new state of the stream
164			* \relates vec3
165			*/
166		✗	inline std::istream& operator>>(std::istream& in, vec3& v) {
167		✗	return in >> v.x >> v.y >> v.z;
168			}
169
170			/**
171			* \brief Computes the dot product between two vectors
172			* \param[in] v1 , v2 const references to the two vectors
173			* \return the dot product between \p v1 and \p v2
174			* \relates vec3
175			*/
176		✗	inline double dot(const vec3& v1, const vec3& v2) {
177		✗	return v1.x * v2.x + v1.y * v2.y + v1.z * v2.z;
178			}
179
180			/**
181			* \brief Computes the cross product between two vectors
182			* \param[in] v1 , v2 const references to the two vectors
183			* \return the cross product between \p v1 and \p v2
184			* \relates vec3
185			*/
186		✗	inline vec3 cross(const vec3& v1, const vec3& v2) {
187			return vec3(
188		✗	v1.y*v2.z - v2.y*v1.z,
189		✗	v1.z*v2.x - v2.z*v1.x,
190		✗	v1.x*v2.y - v2.x*v1.y
191		✗	);
192			}
193
194			/**
195			* \brief Computes the sum of two vectors
196			* \param[in] v1 , v2 const references to the two vectors
197			* \return the sum of \p v1 and \p v2
198			* \relates vec3
199			*/
200			inline vec3 operator+(const vec3& v1, const vec3& v2) {
201			return vec3(
202			v1.x + v2.x,
203			v1.y + v2.y,
204			v1.z + v2.z
205			);
206			}
207
208			/**
209			* \brief Computes the difference between two vectors
210			* \param[in] v1 , v2 const references to the two vectors
211			* \return the difference between \p v1 and \p v2
212			* \relates vec3
213			*/
214		✗	inline vec3 operator-(const vec3& v1, const vec3& v2) {
215			return vec3(
216		✗	v1.x - v2.x,
217		✗	v1.y - v2.y,
218		✗	v1.z - v2.z
219		✗	);
220			}
221
222			/**
223			* \brief Computes the sum of two vectors
224			* \param[in] v1 , v2 const references to the two vectors
225			* \return the sum of \p v1 and \p v2
226			* \relates vec2
227			*/
228			inline vec2 operator+(const vec2& v1, const vec2& v2) {
229			return vec2(
230			v1.x + v2.x,
231			v1.y + v2.y
232			);
233			}
234
235			/**
236			* \brief Computes the difference between two vectors
237			* \param[in] v1 , v2 const references to the two vectors
238			* \return the difference between \p v1 and \p v2
239			* \relates vec2
240			*/
241		✗	inline vec2 operator-(const vec2& v1, const vec2& v2) {
242		✗	return vec2(
243		✗	v1.x - v2.x,
244		✗	v1.y - v2.y
245		✗	);
246			}
247
248			/******************************************************************************/
249			/* Mesh class */
250
251			/**
252			* \brief A vertex in an IndexedMesh
253			* \relates IndexedMesh
254			*/
255			class Vertex {
256			public:
257			/**
258			* \brief Vertex constructor.
259			*/
260			Vertex() : locked(false) {
261			}
262
263			/**
264			* \brief Vertex constructor from 3D point and texture
265			* coordinates.
266			* \param[in] p the 3D coordinates of the vertex
267			* \param[in] t the texture coordinates associated with the vertex
268			*/
269		✗	Vertex(
270			const vec3& p, const vec2& t
271		✗	) : point(p), tex_coord(t), locked(false) {
272		✗	}
273
274			/**
275			* \brief The 3D coordinates.
276			*/
277			vec3 point;
278
279			/**
280			* \brief The texture coordinates (2D).
281			*/
282			vec2 tex_coord;
283
284			/**
285			* \brief A boolean flag that indicates whether the vertex is
286			* locked, i.e. considered as constant in the optimizations.
287			*/
288			bool locked;
289			};
290
291
292			/**
293			* \brief A minimum mesh class.
294			* \details It does not have facet adjacency information
295			* (we do not need it for LSCM), it just stores for each facet the indices
296			* of its vertices. It has load() and save() functions that use the
297			* Alias Wavefront .obj file format.
298			*/
299			class IndexedMesh {
300			public:
301
302			/**
303			* \brief IndexedMesh constructor
304			*/
305		✗	IndexedMesh() : in_facet(false) {
306		✗	facet_ptr.push_back(0);
307		✗	}
308
309			/**
310			* \brief Gets the number of vertices.
311			* \return the number of vertices in this mesh.
312			*/
313		✗	NLuint nb_vertices() const {
314		✗	return NLuint(vertex.size());
315			}
316
317			/**
318			* \brief Gets the number of facets.
319			* \return the number of facets in this mesh.
320			*/
321		✗	NLuint nb_facets() const {
322		✗	return NLuint(facet_ptr.size()-1);
323			}
324
325			/**
326			* \brief Gets the number of vertices in a facet.
327			* \param[in] f the facet, in 0..nb_facets()-1
328			* \return the number of vertices in facet \p f
329			* \pre f < nb_facets()
330			*/
331		✗	NLuint facet_nb_vertices(NLuint f) {
332		✗	assert(f < nb_facets());
333		✗	return facet_ptr[f+1]-facet_ptr[f];
334			}
335
336			/**
337			* \brief Gets a facet vertex by facet index and
338			* local vertex index in facet.
339			* \param[in] f the facet, in 0..nb_facets()-1
340			* \param[in] lv the local vertex index in the facet,
341			* in 0..facet_nb_vertices(f)-1
342			* \return the global vertex index, in 0..nb_vertices()-1
343			* \pre f<nb_facets() && lv < facet_nb_vertices(f)
344			*/
345		✗	NLuint facet_vertex(NLuint f, NLuint lv) {
346		✗	assert(f < nb_facets());
347		✗	assert(lv < facet_nb_vertices(f));
348		✗	return corner[facet_ptr[f] + lv];
349			}
350
351			/**
352			* \brief Adds a new vertex to the mesh.
353			*/
354			void add_vertex() {
355			vertex.push_back(Vertex());
356			}
357
358			/**
359			* \brief Adds a new vertex to the mesh.
360			* \param[in] p the 3D coordinates of the vertex
361			* \param[in] t the texture coordinates of the vertex
362			*/
363		✗	void add_vertex(const vec3& p, const vec2& t) {
364		✗	vertex.push_back(Vertex(p,t));
365		✗	}
366
367			/**
368			* \brief Stats a new facet.
369			*/
370		✗	void begin_facet() {
371		✗	assert(!in_facet);
372		✗	in_facet = true;
373		✗	}
374
375			/**
376			* \brief Terminates the current facet.
377			*/
378		✗	void end_facet() {
379		✗	assert(in_facet);
380		✗	in_facet = false;
381		✗	facet_ptr.push_back(NLuint(corner.size()));
382		✗	}
383
384			/**
385			* \brief Adds a vertex to the current facet.
386			* \param[in] v the index of the vertex
387			* \pre v < vertex.size()
388			*/
389		✗	void add_vertex_to_facet(NLuint v) {
390		✗	assert(in_facet);
391		✗	assert(v < vertex.size());
392		✗	corner.push_back(v);
393		✗	}
394
395			/**
396			* \brief Removes all vertices and all facets from
397			* this mesh.
398			*/
399		✗	void clear() {
400		✗	vertex.clear();
401		✗	corner.clear();
402		✗	facet_ptr.clear();
403		✗	facet_ptr.push_back(0);
404		✗	}
405
406			/**
407			* \brief Loads a file in Alias Wavefront OFF format.
408			* \param[in] file_name the name of the file.
409			*/
410		✗	void load(const std::string& file_name) {
411		✗	std::ifstream input(file_name.c_str());
412		✗	clear();
413		✗	while(input) {
414		✗	std::string line;
415		✗	std::getline(input, line);
416		✗	std::stringstream line_input(line);
417		✗	std::string keyword;
418		✗	line_input >> keyword;
419		✗	if(keyword == "v") {
420		✗	vec3 p;
421		✗	line_input >> p;
422		✗	add_vertex(p,vec2(0.0,0.0));
423		✗	} else if(keyword == "vt") {
424			// Ignore tex vertices
425		✗	} else if(keyword == "f") {
426		✗	begin_facet();
427		✗	while(line_input) {
428		✗	std::string s;
429		✗	line_input >> s;
430		✗	if(s.length() > 0) {
431		✗	std::stringstream v_input(s.c_str());
432			NLuint index;
433		✗	v_input >> index;
434		✗	add_vertex_to_facet(index - 1);
435			char c;
436		✗	v_input >> c;
437		✗	if(c == '/') {
438		✗	v_input >> index;
439			// Ignore tex vertex index
440			}
441		✗	}
442		✗	}
443		✗	end_facet();
444			}
445		✗	}
446		✗	std::cout << "Loaded " << vertex.size() << " vertices and "
447		✗	<< nb_facets() << " facets" << std::endl;
448		✗	}
449
450			/**
451			* \brief Saves a file in Alias Wavefront OFF format.
452			* \param[in] file_name the name of the file.
453			*/
454		✗	void save(const std::string& file_name) {
455		✗	std::ofstream out(file_name.c_str());
456		✗	for(NLuint v=0; v<nb_vertices(); ++v) {
457		✗	out << "v " << vertex[v].point << std::endl;
458			}
459		✗	for(NLuint v=0; v<nb_vertices(); ++v) {
460		✗	out << "vt " << vertex[v].tex_coord << std::endl;
461			}
462		✗	for(NLuint f=0; f<nb_facets(); ++f) {
463		✗	NLuint nv = facet_nb_vertices(f);
464		✗	out << "f ";
465		✗	for(NLuint lv=0; lv<nv; ++lv) {
466		✗	NLuint v = facet_vertex(f,lv);
467		✗	out << (v + 1) << "/" << (v + 1) << " ";
468			}
469		✗	out << std::endl;
470			}
471		✗	for(NLuint v=0; v<nb_vertices(); ++v) {
472		✗	if(vertex[v].locked) {
473		✗	out << "# anchor " << v+1 << std::endl;
474			}
475			}
476		✗	}
477
478			std::vector<Vertex> vertex;
479			bool in_facet;
480
481			/**
482			* \brief All the vertices associated with the facet corners.
483			*/
484			std::vector<NLuint> corner;
485
486			/**
487			* \brief Indicates where facets start and end within the corner
488			* array (facet indicence matrix is stored in the compressed row
489			* storage format).
490			* \details The corners associated with facet f are in the range
491			* facet_ptr[f] ... facet_ptr[f+1]-1
492			*/
493			std::vector<NLuint> facet_ptr;
494			};
495
496			/**
497			* \brief Computes Least Squares Conformal Maps in least squares or
498			* spectral mode.
499			* \details The method is described in the following references:
500			* - Least Squares Conformal Maps, Levy, Petitjean, Ray, Maillot, ACM
501			* SIGGRAPH, 2002
502			* - Spectral Conformal Parameterization, Mullen, Tong, Alliez, Desbrun,
503			* Computer Graphics Forum (SGP conf. proc.), 2008
504			*/
505			class LSCM {
506			public:
507
508			/**
509			* \brief LSCM constructor
510			* \param[in] M a reference to the mesh. It needs to correspond to a
511			* topological disk (open surface with one border and no handle).
512			*/
513		✗	LSCM(IndexedMesh& M) : mesh_(&M) {
514		✗	spectral_ = false;
515		✗	}
516
517			/**
518			* \brief Sets whether spectral mode is used.
519			* \details In default mode, the trivial solution (all vertices to zero)
520			* is avoided by locking two vertices (that are as "extremal" as possible).
521			* In spectral mode, the trivial solution is avoided by finding the first
522			* minimizer that is orthogonal to it (more elegant, but more costly).
523			*/
524		✗	void set_spectral(bool x) {
525		✗	spectral_ = x;
526		✗	}
527
528			/**
529			* \brief Computes the least squares conformal map and stores it in
530			* the texture coordinates of the mesh.
531			* \details Outline of the algorithm (steps 1,2,3 are not used
532			* in spetral mode):
533			* - 1) Find an initial solution by projecting on a plane
534			* - 2) Lock two vertices of the mesh
535			* - 3) Copy the initial u,v coordinates to OpenNL
536			* - 4) Construct the LSCM equation with OpenNL
537			* - 5) Solve the equation with OpenNL
538			* - 6) Copy OpenNL solution to the u,v coordinates
539			*/
540
541		✗	void apply() {
542		✗	const int nb_eigens = 10;
543		✗	nlNewContext();
544		✗	if(spectral_) {
545		✗	if(nlInitExtension("ARPACK")) {
546		✗	std::cout << "ARPACK extension initialized"
547		✗	<< std::endl;
548			} else {
549		✗	std::cout << "Could not initialize ARPACK extension"
550		✗	<< std::endl;
551		✗	exit(-1);
552			}
553		✗	nlEigenSolverParameteri(NL_EIGEN_SOLVER, NL_ARPACK_EXT);
554		✗	nlEigenSolverParameteri(NL_NB_EIGENS, nb_eigens);
555		✗	nlEnable(NL_VERBOSE);
556			}
557		✗	NLuint nb_vertices = NLuint(mesh_->vertex.size());
558		✗	if(!spectral_) {
559		✗	project();
560			}
561		✗	nlSolverParameteri(NL_NB_VARIABLES, NLint(2*nb_vertices));
562		✗	nlSolverParameteri(NL_LEAST_SQUARES, NL_TRUE);
563		✗	nlSolverParameteri(NL_MAX_ITERATIONS, NLint(5*nb_vertices));
564		✗	if(spectral_) {
565		✗	nlSolverParameterd(NL_THRESHOLD, 0.0);
566			} else {
567		✗	nlSolverParameterd(NL_THRESHOLD, 1e-6);
568			}
569		✗	nlBegin(NL_SYSTEM);
570		✗	mesh_to_solver();
571		✗	nlBegin(NL_MATRIX);
572		✗	setup_lscm();
573		✗	nlEnd(NL_MATRIX);
574		✗	nlEnd(NL_SYSTEM);
575		✗	std::cout << "Solving ..." << std::endl;
576
577		✗	if(spectral_) {
578		✗	nlEigenSolve();
579		✗	for(NLuint i=0; i<nb_eigens; ++i) {
580		✗	std::cerr << "[" << i << "] "
581		✗	<< nlGetEigenValue(i) << std::endl;
582			}
583
584			// Find first "non-zero" eigenvalue
585		✗	double small_eigen = ::fabs(nlGetEigenValue(0)) ;
586		✗	eigen_ = 1;
587		✗	for(NLuint i=1; i<nb_eigens; ++i) {
588		✗	if(::fabs(nlGetEigenValue(i)) / small_eigen > 1e3) {
589		✗	eigen_ = i ;
590		✗	break ;
591			}
592			}
593			} else{
594		✗	nlSolve();
595			}
596
597		✗	solver_to_mesh();
598		✗	normalize_uv();
599
600		✗	if(!spectral_) {
601			double time;
602			NLint iterations;
603		✗	nlGetDoublev(NL_ELAPSED_TIME, &time);
604		✗	nlGetIntegerv(NL_USED_ITERATIONS, &iterations);
605		✗	std::cout << "Solver time: " << time << std::endl;
606		✗	std::cout << "Used iterations: " << iterations << std::endl;
607			}
608
609		✗	nlDeleteContext(nlGetCurrent());
610		✗	}
611
612			protected:
613
614			/**
615			* \brief Creates the LSCM equations in OpenNL.
616			*/
617		✗	void setup_lscm() {
618		✗	for(NLuint f=0; f<mesh_->nb_facets(); ++f) {
619		✗	setup_lscm(f);
620			}
621		✗	}
622
623			/**
624			* \brief Creates the LSCM equations in OpenNL, related
625			* with a given facet.
626			* \param[in] f the index of the facet.
627			* \details no-need to triangulate the facet,
628			* we do that "virtually", by creating triangles
629			* radiating around vertex 0 of the facet.
630			* (however, this may be invalid for concave facets)
631			*/
632		✗	void setup_lscm(NLuint f) {
633		✗	NLuint nv = mesh_->facet_nb_vertices(f);
634		✗	for(NLuint i=1; i<nv-1; ++i) {
635		✗	setup_conformal_map_relations(
636		✗	mesh_->facet_vertex(f,0),
637		✗	mesh_->facet_vertex(f,i),
638		✗	mesh_->facet_vertex(f,i+1)
639			);
640			}
641		✗	}
642
643			/**
644			* \brief Computes the coordinates of the vertices of a triangle
645			* in a local 2D orthonormal basis of the triangle's plane.
646			* \param[in] p0 , p1 , p2 the 3D coordinates of the vertices of
647			* the triangle
648			* \param[out] z0 , z1 , z2 the 2D coordinates of the vertices of
649			* the triangle
650			*/
651		✗	static void project_triangle(
652			const vec3& p0,
653			const vec3& p1,
654			const vec3& p2,
655			vec2& z0,
656			vec2& z1,
657			vec2& z2
658			) {
659		✗	vec3 X = p1 - p0;
660		✗	X.normalize();
661		✗	vec3 Z = cross(X,(p2 - p0));
662		✗	Z.normalize();
663		✗	vec3 Y = cross(Z,X);
664		✗	const vec3& O = p0;
665
666		✗	double x0 = 0;
667		✗	double y0 = 0;
668		✗	double x1 = (p1 - O).length();
669		✗	double y1 = 0;
670		✗	double x2 = dot((p2 - O),X);
671		✗	double y2 = dot((p2 - O),Y);
672
673		✗	z0 = vec2(x0,y0);
674		✗	z1 = vec2(x1,y1);
675		✗	z2 = vec2(x2,y2);
676		✗	}
677
678			/**
679			* \brief Creates the LSCM equation in OpenNL, related with
680			* a given triangle, specified by vertex indices.
681			* \param[in] v0 , v1 , v2 the indices of the three vertices of
682			* the triangle.
683			* \details Uses the geometric form of LSCM equation:
684			* (Z1 - Z0)(U2 - U0) = (Z2 - Z0)(U1 - U0)
685			* Where Uk = uk + i.vk is the complex number
686			* corresponding to (u,v) coords
687			* Zk = xk + i.yk is the complex number
688			* corresponding to local (x,y) coords
689			* There is no divide with this expression,
690			* this makes it more numerically stable in
691			* the presence of degenerate triangles.
692			*/
693		✗	void setup_conformal_map_relations(
694			NLuint v0, NLuint v1, NLuint v2
695			) {
696
697		✗	const vec3& p0 = mesh_->vertex[v0].point;
698		✗	const vec3& p1 = mesh_->vertex[v1].point;
699		✗	const vec3& p2 = mesh_->vertex[v2].point;
700
701		✗	vec2 z0,z1,z2;
702		✗	project_triangle(p0,p1,p2,z0,z1,z2);
703		✗	vec2 z01 = z1 - z0;
704		✗	vec2 z02 = z2 - z0;
705		✗	double a = z01.x;
706		✗	double b = z01.y;
707		✗	double c = z02.x;
708		✗	double d = z02.y;
709		✗	assert(b == 0.0);
710
711			// Note : 2*id + 0 --> u
712			// 2*id + 1 --> v
713		✗	NLuint u0_id = 2*v0 ;
714		✗	NLuint v0_id = 2*v0 + 1;
715		✗	NLuint u1_id = 2*v1 ;
716		✗	NLuint v1_id = 2*v1 + 1;
717		✗	NLuint u2_id = 2*v2 ;
718		✗	NLuint v2_id = 2*v2 + 1;
719
720			// Note : b = 0
721
722			// Real part
723		✗	nlBegin(NL_ROW);
724		✗	nlCoefficient(u0_id, -a+c) ;
725		✗	nlCoefficient(v0_id, b-d) ;
726		✗	nlCoefficient(u1_id, -c) ;
727		✗	nlCoefficient(v1_id, d) ;
728		✗	nlCoefficient(u2_id, a);
729		✗	nlEnd(NL_ROW);
730
731			// Imaginary part
732		✗	nlBegin(NL_ROW);
733		✗	nlCoefficient(u0_id, -b+d);
734		✗	nlCoefficient(v0_id, -a+c);
735		✗	nlCoefficient(u1_id, -d);
736		✗	nlCoefficient(v1_id, -c);
737		✗	nlCoefficient(v2_id, a);
738		✗	nlEnd(NL_ROW);
739		✗	}
740
741			/**
742			* \brief Copies u,v coordinates from OpenNL solver to the mesh.
743			*/
744		✗	void solver_to_mesh() {
745		✗	for(NLuint i=0; i<mesh_->vertex.size(); ++i) {
746		✗	Vertex& it = mesh_->vertex[i];
747		✗	double u = spectral_ ? nlMultiGetVariable(2 * i ,eigen_)
748		✗	: nlGetVariable(2 * i );
749		✗	double v = spectral_ ? nlMultiGetVariable(2 * i + 1,eigen_)
750		✗	: nlGetVariable(2 * i + 1);
751		✗	it.tex_coord = vec2(u,v);
752			}
753		✗	}
754
755			/**
756			* \brief Translates and scales tex coords in such a way that they fit
757			* within the unit square.
758			*/
759		✗	void normalize_uv() {
760		✗	double u_min=1e30, v_min=1e30, u_max=-1e30, v_max=-1e30;
761		✗	for(NLuint i=0; i<mesh_->vertex.size(); ++i) {
762		✗	u_min = std::min(u_min, mesh_->vertex[i].tex_coord.x);
763		✗	v_min = std::min(v_min, mesh_->vertex[i].tex_coord.y);
764		✗	u_max = std::max(u_max, mesh_->vertex[i].tex_coord.x);
765		✗	v_max = std::max(v_max, mesh_->vertex[i].tex_coord.y);
766			}
767		✗	double l = std::max(u_max-u_min,v_max-v_min);
768		✗	for(NLuint i=0; i<mesh_->vertex.size(); ++i) {
769		✗	mesh_->vertex[i].tex_coord.x -= u_min;
770		✗	mesh_->vertex[i].tex_coord.x /= l;
771		✗	mesh_->vertex[i].tex_coord.y -= v_min;
772		✗	mesh_->vertex[i].tex_coord.y /= l;
773			}
774		✗	}
775
776			/**
777			* \brief Copies u,v coordinates from the mesh to OpenNL solver.
778			*/
779		✗	void mesh_to_solver() {
780		✗	for(NLuint i=0; i<mesh_->vertex.size(); ++i) {
781		✗	Vertex& it = mesh_->vertex[i];
782		✗	double u = it.tex_coord.x;
783		✗	double v = it.tex_coord.y;
784		✗	nlSetVariable(2 * i , u);
785		✗	nlSetVariable(2 * i + 1, v);
786		✗	if(!spectral_ && it.locked) {
787		✗	nlLockVariable(2 * i );
788		✗	nlLockVariable(2 * i + 1);
789			}
790			}
791		✗	}
792
793			/**
794			* \brief Chooses an initial solution, and locks two vertices.
795			*/
796		✗	void project() {
797			// Get bbox
798			unsigned int i;
799
800		✗	double xmin = 1e30;
801		✗	double ymin = 1e30;
802		✗	double zmin = 1e30;
803		✗	double xmax = -1e30;
804		✗	double ymax = -1e30;
805		✗	double zmax = -1e30;
806
807		✗	for(i=0; i<mesh_->vertex.size(); i++) {
808		✗	const Vertex& v = mesh_->vertex[i];
809		✗	xmin = std::min(v.point.x, xmin);
810		✗	ymin = std::min(v.point.y, ymin);
811		✗	zmin = std::min(v.point.z, zmin);
812
813		✗	xmax = std::max(v.point.x, xmax);
814		✗	ymax = std::max(v.point.y, ymax);
815		✗	zmax = std::max(v.point.z, zmax);
816			}
817
818		✗	double dx = xmax - xmin;
819		✗	double dy = ymax - ymin;
820		✗	double dz = zmax - zmin;
821
822		✗	vec3 V1,V2;
823
824			// Find shortest bbox axis
825		✗	if(dx <= dy && dx <= dz) {
826		✗	if(dy > dz) {
827		✗	V1 = vec3(0,1,0);
828		✗	V2 = vec3(0,0,1);
829			} else {
830		✗	V2 = vec3(0,1,0);
831		✗	V1 = vec3(0,0,1);
832			}
833		✗	} else if(dy <= dx && dy <= dz) {
834		✗	if(dx > dz) {
835		✗	V1 = vec3(1,0,0);
836		✗	V2 = vec3(0,0,1);
837			} else {
838		✗	V2 = vec3(1,0,0);
839		✗	V1 = vec3(0,0,1);
840			}
841		✗	} else if(dz <= dx && dz <= dy) {
842		✗	if(dx > dy) {
843		✗	V1 = vec3(1,0,0);
844		✗	V2 = vec3(0,1,0);
845			} else {
846		✗	V2 = vec3(1,0,0);
847		✗	V1 = vec3(0,1,0);
848			}
849			}
850
851			// Project onto shortest bbox axis,
852			// and lock extrema vertices
853
854		✗	Vertex* vxmin = nullptr;
855		✗	double umin = 1e30;
856		✗	Vertex* vxmax = nullptr;
857		✗	double umax = -1e30;
858
859		✗	for(i=0; i<mesh_->vertex.size(); i++) {
860		✗	Vertex& V = mesh_->vertex[i];
861		✗	double u = dot(V.point,V1);
862		✗	double v = dot(V.point,V2);
863		✗	V.tex_coord = vec2(u,v);
864		✗	if(u < umin) {
865		✗	vxmin = &V;
866		✗	umin = u;
867			}
868		✗	if(u > umax) {
869		✗	vxmax = &V;
870		✗	umax = u;
871			}
872			}
873
874		✗	vxmin->locked = true;
875		✗	vxmax->locked = true;
876		✗	}
877
878			IndexedMesh* mesh_;
879
880			/**
881			* \brief true if spectral mode is used,
882			* false if locked least squares mode is used.
883			*/
884			bool spectral_;
885
886			/**
887			* \brief In spectral mode, the index of the first
888			* non-zero eigenvalue.
889			*/
890			NLuint eigen_;
891			};
892
893
894		✗	int main(int argc, char** argv) {
895		✗	bool spectral = false;
896		✗	bool OK = true;
897		✗	std::vector<std::string> filenames;
898
899		✗	nlInitialize(argc, argv);
900
901		✗	for(int i=1; i<argc; ++i) {
902		✗	if(!strcmp(argv[i],"spectral=true")) {
903		✗	spectral = true;
904		✗	} else if(!strcmp(argv[i],"spectral=false")) {
905		✗	spectral = false;
906		✗	} else if(strchr(argv[i],'=') == nullptr) {
907		✗	filenames.push_back(argv[i]);
908			}
909			}
910
911		✗	OK = OK && (filenames.size() >= 1) && (filenames.size() <= 2);
912
913		✗	if(!OK) {
914			std::cerr << "usage: " << argv[0]
915		✗	<< " infile.obj <outfile.obj> <spectral=true|false>"
916		✗	<< std::endl;
917		✗	return -1;
918			}
919
920		✗	if(filenames.size() == 1) {
921		✗	filenames.push_back("out.obj");
922			}
923
924		✗	IndexedMesh mesh;
925		✗	std::cout << "Loading " << filenames[0] << " ..." << std::endl;
926		✗	mesh.load(filenames[0]);
927
928		✗	LSCM lscm(mesh);
929		✗	lscm.set_spectral(spectral);
930		✗	lscm.apply();
931
932		✗	std::cout << "Saving " << filenames[1] << " ..." << std::endl;
933		✗	mesh.save(filenames[1]);
934		✗	}
935