GCC Code Coverage Report

Line	Branch	Exec	Source
1			/*
2			* Copyright (c) 2000-2022 Inria
3			* All rights reserved.
4			*
5			* Redistribution and use in source and binary forms, with or without
6			* modification, are permitted provided that the following conditions are met:
7			*
8			* * Redistributions of source code must retain the above copyright notice,
9			* this list of conditions and the following disclaimer.
10			* * Redistributions in binary form must reproduce the above copyright notice,
11			* this list of conditions and the following disclaimer in the documentation
12			* and/or other materials provided with the distribution.
13			* * Neither the name of the ALICE Project-Team nor the names of its
14			* contributors may be used to endorse or promote products derived from this
15			* software without specific prior written permission.
16			*
17			* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
18			* AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
19			* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
20			* ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
21			* LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
22			* CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
23			* SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
24			* INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
25			* CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
26			* ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
27			* POSSIBILITY OF SUCH DAMAGE.
28			*
29			* Contact: Bruno Levy
30			*
31			* https://www.inria.fr/fr/bruno-levy
32			*
33			* Inria,
34			* Domaine de Voluceau,
35			* 78150 Le Chesnay - Rocquencourt
36			* FRANCE
37			*
38			*/
39
40			#include <geogram/mesh/mesh_frame_field.h>
41			#include <geogram/mesh/mesh.h>
42			#include <geogram/mesh/mesh_geometry.h>
43			#include <geogram/numerics/matrix_util.h>
44			#include <geogram/basic/logger.h>
45			#include <geogram/basic/line_stream.h>
46			#include <geogram/basic/progress.h>
47			#include <geogram/bibliography/bibliography.h>
48			#include <geogram/NL/nl.h>
49
50			// Some member functions of NormalCycle are not used here.
51			// note: NormalCycle will be exported sometime, so for now
52			// we deactivate the warning.
53
54			#ifdef __clang__
55			#pragma GCC diagnostic ignored "-Wunused-member-function"
56			#endif
57
58			#ifdef __ICC
59			#pragma warning disable 177
60			#endif
61
62			namespace {
63			using namespace GEO;
64
65			/**
66			* \brief We manipulate 4-symmetry direction fields [Ray et.al]
67			*/
68			const double symd = 4.0;
69
70			/**
71			* \brief Represents a local orthonormal basis
72			* of a mesh facet.
73			*/
74			class MeshFacetBasis {
75			public:
76
77			/**
78			* \brief Constructs a new MeshFacetBasis.
79			* \param[in] M a const reference to the mesh
80			* \param[in] f the index of the facet in \p M
81			*/
82		✗	MeshFacetBasis(
83			const Mesh& M, index_t f
84		✗	) {
85		✗	X = normalize(
86		✗	Geom::mesh_corner_vector(M, M.facets.corners_begin(f))
87			);
88		✗	N = normalize(Geom::mesh_facet_normal(M,f));
89		✗	Y = cross(N,X);
90		✗	}
91
92			/**
93			* \brief Transforms a 3d vector into the local
94			* 2d basis.
95			* \param[in] v the input 3d vector
96			* \return the representation of \p v in the local
97			* 2d basis.
98			*/
99			vec2 project(const vec3& v) const {
100			return vec2(dot(v,X),dot(v,Y));
101			}
102
103			/**
104			* \brief Transforms a local 3d vector into the
105			* global 3d basis.
106			* \param[in] v the input 2d vector in the local basis
107			* \return the representation of \p v in the global
108			* 3d basis.
109			*/
110			vec3 unproject(const vec2& v) const {
111			return v.x*X + v.y*Y;
112			}
113
114			/**
115			* \brief Computes the angles between a 3d vector and
116			* the first axis of the local basis.
117			* \param[in] v the input vector in the global 3d basis
118			* \return the angle between \p v and the first axis of
119			* the local basis
120			*/
121		✗	double angle(const vec3& v) const {
122			vec2 v2=project(v);
123		✗	return atan2(v2.y,v2.x);
124			}
125
126			/**
127			* \brief Computes the rotation angle between the two reference vectors
128			* of two facets that share an edge.
129			* \param[in] M a reference to the mesh
130			* \param[in] c a corner index in \p M
131			* \return the angle between the reference frames of the
132			* two facets sharing the edge originating at \p c1
133			*/
134		✗	static double reference_rotation_accross_edge(
135			const Mesh& M, index_t c
136			) {
137			geo_debug_assert(M.facets.are_simplices());
138		✗	index_t f1 = c/3;
139			index_t f2 = M.facet_corners.adjacent_facet(c);
140			geo_debug_assert(f2 != NO_FACET);
141		✗	vec3 ref = Geom::mesh_corner_vector(M,c);
142		✗	MeshFacetBasis B1(M,f1);
143		✗	MeshFacetBasis B2(M,f2);
144		✗	return B2.angle(ref) - B1.angle(ref);
145			}
146
147			private:
148			vec3 X;
149			vec3 Y;
150			vec3 N;
151			};
152
153
154			/**
155			* \brief Solves for Periodic Global Parameterization variables.
156			* \details This function computes for each facet of the mesh \p M
157			* two variables, that correspond to the cosine and sine of an angle
158			* interpolated over the mesh. This angle is relative to the first
159			* edge of the facet, as defined in the MeshFacetBasis class.
160			* \param[in] M a const reference to the surface mesh
161			* \param[in,out] sincos_alpha a vector of 2*M.facets.nb() variables,
162			* that correspond to the interpolated variables. The initial value
163			* is taken into account in the fitting term if \p fitting is non-zero
164			* \param[in] locked a vector of M.facets.nb() booleans, indicating
165			* whether each facet is locked. The variables that correspond to
166			* a locked facet are unchanged.
167			* \param[in] global_fitting importance of the fitting term with respect to
168			* the initial value of \p sincos_alpha. If zero, no fitting term is
169			* installed
170			* \param[in] local_fitting an optional vector of M.facets.nb() doubles
171			* that specifies for each facet an individual factor that scales
172			* global_fitting
173			*/
174		✗	void solve_PGP(
175			const Mesh& M,
176			vector<double>& sincos_alpha,
177			const vector<bool>& locked,
178			double global_fitting,
179			const vector<double>& local_fitting = vector<double>()
180			) {
181			// Step 0: normalize variables
182		✗	for(index_t f: M.facets) {
183		✗	double c = sincos_alpha[2*f];
184		✗	double s = sincos_alpha[2*f+1];
185		✗	double scale = sqrt(s*s+c*c);
186		✗	if(scale > 1e-30) {
187		✗	sincos_alpha[2*f] = c/scale;
188		✗	sincos_alpha[2*f+1] = s/scale;
189			}
190			}
191
192			// Step 1: Setup the OpenNL solver
193		✗	nlNewContext();
194		✗	nlSolverParameteri(NL_NB_VARIABLES, NLint(2*M.facets.nb()));
195		✗	nlSolverParameteri(NL_LEAST_SQUARES, NL_TRUE);
196			#ifdef GEO_DEBUG
197			nlEnable(NL_VERBOSE);
198			#endif
199		✗	nlEnable(NL_NORMALIZE_ROWS);
200
201			// Step 2: setup the variables
202		✗	nlBegin(NL_SYSTEM);
203		✗	for(index_t f: M.facets) {
204		✗	nlSetVariable(2*f, sincos_alpha[2*f]);
205		✗	nlSetVariable(2*f+1, sincos_alpha[2*f+1]);
206		✗	if(locked.size() != 0 && locked[f]) {
207		✗	nlLockVariable(2*f);
208		✗	nlLockVariable(2*f+1);
209			}
210			}
211
212		✗	nlBegin(NL_MATRIX);
213
214			// Step 3: setup the PGP smoothness term
215		✗	for(index_t f1: M.facets) {
216		✗	for(index_t c1: M.facets.corners(f1)) {
217			index_t f2 = M.facet_corners.adjacent_facet(c1);
218		✗	if(f2 == NO_FACET || f1 < f2) {
219		✗	continue;
220			}
221
222			double angle = -symd*
223		✗	MeshFacetBasis::reference_rotation_accross_edge(
224			M,c1
225		✗	);
226
227		✗	double c = cos(angle);
228		✗	double s = sin(angle);
229
230		✗	nlBegin(NL_ROW);
231		✗	nlCoefficient(2*f1,c);
232		✗	nlCoefficient(2*f1+1,s);
233		✗	nlCoefficient(2*f2,-1.0);
234		✗	nlEnd(NL_ROW);
235
236		✗	nlBegin(NL_ROW);
237		✗	nlCoefficient(2*f1,-s);
238		✗	nlCoefficient(2*f1+1,c);
239		✗	nlCoefficient(2*f2+1,-1.0);
240		✗	nlEnd(NL_ROW);
241			}
242			}
243
244			// Step 4: setup the data fitting term
245		✗	if(global_fitting != 0) {
246		✗	for(index_t f: M.facets) {
247
248			double fitting = global_fitting;
249		✗	if(local_fitting.size() != 0) {
250		✗	fitting *= local_fitting[f];
251			}
252
253		✗	if(fitting == 0.0) {
254		✗	continue;
255			}
256
257		✗	nlRowScaling(fitting);
258		✗	nlBegin(NL_ROW);
259		✗	nlCoefficient(2*f,1.0);
260		✗	nlRightHandSide(sincos_alpha[2*f]);
261		✗	nlEnd(NL_ROW);
262
263		✗	nlRowScaling(fitting);
264		✗	nlBegin(NL_ROW);
265		✗	nlCoefficient(2*f+1,1.0);
266		✗	nlRightHandSide(sincos_alpha[2*f+1]);
267		✗	nlEnd(NL_ROW);
268			}
269			}
270
271		✗	nlEnd(NL_MATRIX);
272		✗	nlEnd(NL_SYSTEM);
273
274			// Step 5: solve the linear system
275		✗	nlSolve() ;
276
277			// Step 6: read the new values of the variables
278		✗	for(index_t f: M.facets) {
279		✗	sincos_alpha[2*f] = nlGetVariable(2*f);
280		✗	sincos_alpha[2*f+1] = nlGetVariable(2*f+1);
281			}
282
283			// Step 7: cleanup memory allocated by OpenNL
284		✗	nlDeleteContext(nlGetCurrent());
285		✗	}
286
287
288			/**
289			* \brief Estimates the curvature tensor using a set
290			* of samples. Each sample is a vector and a dihedral angle.
291			* \details The algorithm is detailed in the following reference:
292			* Restricted Delaunay Triangulation and Normal Cycle,
293			* D. Cohen-Steiner and J.M. Morvan,
294			* SOCG 2003
295			*/
296			class NormalCycle {
297			public:
298			/**
299			* \brief Constructs a new NormalCycle.
300			*/
301		✗	NormalCycle() {
302			clear();
303			}
304
305			/**
306			* \brief Clears the currently accumulated matrix.
307			*/
308			void clear() {
309		✗	for(index_t i=0; i<6; ++i) {
310		✗	M_[i] = 0.0;
311			}
312			}
313
314			/**
315			* \brief Computes the eigenvalues
316			* and eigenvectors of the accumulated tensor.
317			*/
318		✗	void compute() {
319		✗	double trace = M_[0] + M_[2] + M_[5] ;
320		✗	double s = 1e-6 * trace ;
321		✗	if (trace==0.0) {
322			s = 1e-6;
323			}
324		✗	M_[0] += s ;
325		✗	M_[2] += s ;
326		✗	M_[5] += s ;
327
328
329			double eigen_vectors[9] ;
330		✗	MatrixUtil::semi_definite_symmetric_eigen(
331		✗	M_, 3, eigen_vectors, eigen_value_
332			) ;
333
334		✗	axis_[0] = vec3(
335			eigen_vectors[0], eigen_vectors[1], eigen_vectors[2]
336			);
337
338		✗	axis_[1] = vec3(
339			eigen_vectors[3], eigen_vectors[4], eigen_vectors[5]
340			);
341
342		✗	axis_[2] = vec3(
343			eigen_vectors[6], eigen_vectors[7], eigen_vectors[8]
344			);
345
346			// Normalize the eigen vectors
347
348		✗	for(index_t i=0; i<3; ++i) {
349		✗	axis_[i] = normalize(axis_[i]) ;
350			}
351
352			// Sort the eigen vectors
353		✗	i_[0] = 0 ;
354		✗	i_[1] = 1 ;
355		✗	i_[2] = 2 ;
356
357		✗	double l0 = ::fabs(eigen_value_[0]) ;
358		✗	double l1 = ::fabs(eigen_value_[1]) ;
359		✗	double l2 = ::fabs(eigen_value_[2]) ;
360
361		✗	if(l1 > l0) {
362			std::swap(l0 , l1 ) ;
363			std::swap(i_[0], i_[1]) ;
364			}
365
366		✗	if(l2 > l1) {
367			std::swap(l1 , l2 ) ;
368			std::swap(i_[1], i_[2]) ;
369			}
370		✗	if(l1 > l0) {
371			std::swap(l0 , l1 ) ;
372			std::swap(i_[0],i_[1]) ;
373			}
374		✗	}
375
376			/**
377			* \brief Accumulates a dihedral angle to the current
378			* tensor.
379			* \details This function needs to be called between
380			* a begin() \ end() pair. If a geometric clipping
381			* neighborhood is used, the specified edge vector
382			* needs to be clipped by it.
383			* \param[in] edge the supporting edge of the dihedron
384			* \param[in] angle the angle of the dihedron
385			* \param[in] neigh_area the area of the clipped
386			* neighborhood
387			* to store the accumulated tensor.
388			*/
389		✗	void accumulate_dihedral_angle(
390			const vec3& edge, double angle, double neigh_area=1.0
391			) {
392		✗	vec3 e = normalize(edge) ;
393		✗	double s = length(edge) * angle * neigh_area ;
394		✗	M_[0] += s * e.x * e.x;
395		✗	M_[1] += s * e.x * e.y;
396		✗	M_[2] += s * e.y * e.y;
397		✗	M_[3] += s * e.x * e.z;
398		✗	M_[4] += s * e.y * e.z;
399		✗	M_[5] += s * e.z * e.z;
400		✗	}
401
402			/**
403			* \brief Gets an eigenvector by index
404			* \param[in] i the index of the eigenvector (0,1 or 2).
405			* The eigenvectors are sorted by increasing eigenvalue
406			* magnitude.
407			* \return the \p i%-th eigenvector
408			*/
409			const vec3& eigen_vector(int i) const {
410		✗	return axis_[i_[i]];
411			}
412
413			/**
414			* \brief Gets an eigenvalue by index
415			* \param[in] i the index of the eigenvalue (0,1 or 2).
416			* The eigenvalues are sorted by increasing eigenvalue
417			* magnitude.
418			* \return the \p i%-th eigenvalue
419			*/
420			double eigen_value(int i) const {
421		✗	return eigen_value_[i_[i]];
422			}
423
424			/**
425			* \brief Gets the estimated normal vector.
426			* \return the estimated normal vector
427			*/
428			const vec3& N() const {
429			return eigen_vector(2);
430			}
431
432			/**
433			* \brief Gets the estimated direction of maximum curvature.
434			* \return the estimated direction of maximum curvature
435			*/
436			const vec3& Kmax() const {
437			return eigen_vector(0);
438			}
439
440			/**
441			* \brief Gets the estimated direction of minimum curvature.
442			* \return the estimated direction of minimum curvature
443			*/
444			const vec3& Kmin() const {
445			return eigen_vector(1);
446			}
447
448			/**
449			* \brief Gets the estimated maximum curvature.
450			* \return the estimated maximum curvature
451			*/
452			double kmax() const {
453			return eigen_value(0);
454			}
455
456			/**
457			* \brief Gets the estimated minimum curvature.
458			* \return the estimated minimum curvature
459			*/
460			double kmin() const {
461			return eigen_value(1);
462			}
463
464			/**
465			* \brief Adds the currently accumulated tensor to a
466			* matrix.
467			* \param[out] M an array of 6 doubles that represents
468			* the tensor to which the current tensor will be added
469			*/
470			void add_to_matrix(double* M) const {
471		✗	for(index_t i=0; i<6; ++i) {
472		✗	M[i] += M_[i];
473			}
474			}
475
476			/**
477			* \brief Adds a matrix to the currently accumulated tensor.
478			* \param[in] M an array of 6 doubles that represents the
479			* matrix that should be added to the currently accumulated
480			* tensor.
481			*/
482			void add_matrix(const double* M) {
483		✗	for(index_t i=0; i<6; ++i) {
484		✗	M_[i] += M[i];
485			}
486			}
487
488			private:
489			vec3 axis_[3] ;
490			double eigen_value_[3] ;
491			double M_[6] ;
492			int i_[3] ;
493			};
494
495			/**
496			* \brief Estimates the direction of the maximum principal curvature.
497			* \details The direction of maximum principal curvature is encoded
498			* as the cosine and sine of the angle it makes relative to the first
499			* edge of teach triangle, as defined in the MeshFacetBasis class.
500			* \param[in] M a const reference to the surface mesh
501			* \param[in,out] sincos_alpha a vector of 2*M.facets.nb() doubles,
502			* that contains the cosines and sines of the angle between the estimated
503			* directions and the first edge of each facet.
504			* \param[in] locked a vector of M.facets.nb() booleans, that indicates
505			* for each facet whether it is locked. Directions of locked facets are
506			* kept unchanged
507			* \param[out] magnitude a vector of M.facets.nb() doubles that indicates
508			* for each facet the magnitude of the principal direction of curvature
509			*/
510		✗	void estimate_max_curvature_direction(
511			const Mesh& M, vector<double>& sincos_alpha, const vector<bool>& locked,
512			vector<double>& magnitude
513			) {
514			NormalCycle NC;
515		✗	vector<double> matrices(M.vertices.nb()*6,0.0);
516
517			// Compute tensors of vertex neighborhoods
518		✗	for(index_t f1: M.facets) {
519		✗	for(index_t c: M.facets.corners(f1)) {
520			index_t f2 = M.facet_corners.adjacent_facet(c);
521		✗	if(f2 == NO_FACET || f2 < f1) {
522		✗	continue;
523			}
524
525			index_t v1 = M.facet_corners.vertex(c);
526		✗	index_t v2 = M.facet_corners.vertex(c+1);
527
528		✗	vec3 e = Geom::mesh_corner_vector(M,c);
529		✗	double alpha = Geom::mesh_normal_angle(M,c);
530
531			NC.clear();
532		✗	NC.accumulate_dihedral_angle(e,alpha);
533		✗	NC.add_to_matrix(&matrices[6*v1]);
534		✗	NC.add_to_matrix(&matrices[6*v2]);
535			}
536			}
537
538
539			// For each facet, accumulate the tensors of all its
540			// vertices.
541		✗	for(index_t f: M.facets) {
542		✗	if(locked.size() != 0 && locked[f]) {
543			continue;
544			}
545			NC.clear();
546		✗	for(index_t c: M.facets.corners(f)) {
547			index_t v=M.facet_corners.vertex(c);
548		✗	NC.add_matrix(&matrices[6*v]);
549			}
550		✗	NC.compute();
551
552			// Compute eigenvector and encode it as PGP
553			// representation (cosine and sine of the
554			// angle relative to the first edge of the triangle).
555		✗	vec2 K = MeshFacetBasis(M,f).project(
556			NC.Kmax()
557			);
558		✗	double angle = atan2(K.y,K.x)*symd;
559		✗	sincos_alpha[2*f] = cos(angle);
560		✗	sincos_alpha[2*f+1] = sin(angle);
561
562		✗	magnitude[f] = fabs(NC.kmax());
563			}
564		✗	}
565
566
567			}
568
569			namespace GEO {
570
571		✗	bool FrameField::load(
572			const Mesh& M, bool volumetric, const std::string& filename
573			) {
574		✗	Logger::out("Frames") << "Loading frames from "
575			<< filename << std::endl;
576			frames_.clear();
577		✗	frames_.reserve(M.cells.nb() * 9);
578			centers_.clear();
579			bool result = true;
580			bool with_centers = false;
581			try {
582		✗	LineInput in(filename);
583		✗	if(!in.OK()) {
584			return false;
585			}
586			bool first_line = true;
587		✗	while(!in.eof() && in.get_line()) {
588		✗	in.get_fields();
589		✗	if(first_line) {
590		✗	if(in.nb_fields() == 12) {
591			with_centers = true;
592			}
593			first_line=false;
594			}
595		✗	if(!with_centers && in.nb_fields() != 3) {
596		✗	Logger::err("I/O")
597			<< "Line " << in.line_number()
598			<< ": invalid number of fields (expected 3, got"
599			<< in.nb_fields()
600			<< ")"
601			<< std::endl;
602			result = false;
603		✗	break;
604			}
605		✗	if(with_centers && in.nb_fields() != 12) {
606		✗	Logger::err("I/O")
607			<< "Line " << in.line_number()
608			<< ": invalid number of fields (expected 12, got"
609			<< in.nb_fields()
610			<< ")"
611			<< std::endl;
612			result = false;
613		✗	break;
614			}
615		✗	if(with_centers) {
616		✗	for(index_t i=0; i<3; ++i) {
617		✗	centers_.push_back(in.field_as_double(i));
618			}
619		✗	for(index_t i=3; i<12; ++i) {
620		✗	frames_.push_back(in.field_as_double(i));
621			}
622			} else {
623		✗	frames_.push_back(in.field_as_double(0));
624		✗	frames_.push_back(in.field_as_double(1));
625		✗	frames_.push_back(in.field_as_double(2));
626			}
627			}
628		✗	}
629		✗	catch(const std::exception& ex) {
630		✗	Logger::err("I/O") << ex.what() << std::endl;
631			result = false;
632		✗	}
633		✗	if(!result) {
634		✗	Logger::err("I/O")
635			<< "Could not load file: " << filename
636			<< std::endl;
637		✗	return false;
638			}
639		✗	if(!with_centers) {
640		✗	if(volumetric) {
641		✗	if(frames_.size() != M.cells.nb() * 9) {
642		✗	Logger::err("I/O")
643			<< "Invalid number of elements in frame for: "
644			<< filename << std::endl;
645		✗	return false;
646			}
647		✗	centers_.resize(M.cells.nb() * 3);
648		✗	for(index_t t : M.cells) {
649		✗	vec3 g = Geom::mesh_tet_center(M, t);
650		✗	centers_[3 * t] = g.x;
651		✗	centers_[3 * t + 1] = g.y;
652		✗	centers_[3 * t + 2] = g.z;
653			}
654			} else {
655		✗	if(frames_.size() != M.facets.nb() * 9) {
656		✗	Logger::err("I/O")
657			<< "Invalid number of elements in frame for: "
658			<< filename << std::endl;
659		✗	return false;
660			}
661		✗	centers_.resize(M.facets.nb() * 3);
662		✗	for(index_t f: M.facets) {
663		✗	vec3 g = Geom::mesh_facet_center(M, f);
664		✗	centers_[3 * f] = g.x;
665		✗	centers_[3 * f + 1] = g.y;
666		✗	centers_[3 * f + 2] = g.z;
667			}
668			}
669			}
670		✗	geo_assert((centers_.size()/3)*3 == centers_.size());
671		✗	geo_assert((frames_.size()/9)*9 == frames_.size());
672
673		✗	index_t nb_vectors = frames_.size()/3;
674		✗	for(index_t i=0; i<nb_vectors; ++i) {
675			double s = 0.0;
676		✗	for(index_t c=0; c<3; ++c) {
677		✗	s += frames_[3*i+c]*frames_[3*i+c];
678			}
679		✗	s = ::sqrt(s);
680		✗	if(s == 0.0) {
681		✗	Logger::warn("Frames")
682			<< "Zero-length vector in frame" << std::endl;
683			} else {
684		✗	for(index_t c=0; c<3; ++c) {
685		✗	frames_[3*i+c] /= s;
686			}
687			}
688			}
689
690
691		✗	Logger::out("Frames") << "Loaded " << centers_.size()/3
692			<< " frames" << std::endl;
693		✗	Logger::out("Frames") << "Creating NN search" << std::endl;
694		✗	NN_ = NearestNeighborSearch::create(3, "default");
695		✗	NN_->set_points(centers_.size()/3, centers_.data());
696		✗	return true;
697			}
698
699		✗	void FrameField::create_from_surface_mesh(
700			const Mesh& M, bool volumetric, double sharp_angle_threshold
701			) {
702
703		✗	geo_cite("DBLP:journals/tog/RayVLL08");
704		✗	geo_cite("DBLP:journals/tog/RayVAL09");
705
706		✗	sharp_angle_threshold *= M_PI/180.0 ;
707
708		✗	vector<double> alpha_sincos(2*M.facets.nb(),0.0);
709			vector<bool> locked(M.facets.nb());
710
711			// Step 1: setup the fixed variables
712			index_t nb_constrained = 0;
713
714		✗	for(index_t f1: M.facets) {
715		✗	for(index_t c1: M.facets.corners(f1)) {
716			index_t f2 = M.facet_corners.adjacent_facet(c1);
717			if(
718		✗	f2 == NO_FACET || (
719		✗	::fabs(Geom::mesh_normal_angle(M,c1)) >
720			sharp_angle_threshold
721			)
722			) {
723		✗	vec2 v = MeshFacetBasis(M,f1).project(
724		✗	Geom::mesh_corner_vector(M,c1)
725			);
726		✗	double angle = atan2(v.y,v.x)*symd;
727
728		✗	locked[f1]=true;
729		✗	alpha_sincos[2*f1] = cos(angle);
730		✗	alpha_sincos[2*f1+1] = sin(angle);
731
732		✗	++nb_constrained;
733			}
734			}
735			}
736
737		✗	Logger::out("Frames")
738			<< nb_constrained << " constrained edges" << std::endl;
739
740			vector<double> certainty(M.facets.nb());
741		✗	estimate_max_curvature_direction(M,alpha_sincos,locked,certainty);
742		✗	double max_certainty = 0.0;
743		✗	for(index_t f: M.facets) {
744		✗	max_certainty = std::max(max_certainty,certainty[f]);
745			}
746		✗	for(index_t f: M.facets) {
747		✗	certainty[f] /= max_certainty;
748		✗	if(Numeric::is_nan(certainty[f])) {
749		✗	certainty[f] = 0.0;
750			}
751			}
752
753			// Step 2: solve for sines and cosines
754			// (Periodic Global Parameterization)
755			try {
756		✗	ProgressTask progress("Frames Smth.",4);
757		✗	for(index_t k=0; k<5; ++k) {
758		✗	solve_PGP(M,alpha_sincos,locked,1.0,certainty);
759		✗	progress.progress(k);
760			}
761		✗	} catch(const TaskCanceled&) {
762		✗	}
763
764
765			// Step 3: deduce the frame field from
766			// the solution of the linear system
767		✗	frames_.resize(M.facets.nb()*9);
768		✗	centers_.resize(M.facets.nb()*3);
769		✗	for(index_t f: M.facets) {
770		✗	double angle = atan2(
771		✗	alpha_sincos[2*f+1],
772		✗	alpha_sincos[2*f]
773		✗	) / symd;
774		✗	vec3 U = MeshFacetBasis(M,f).unproject(vec2(cos(angle),sin(angle)));
775		✗	vec3 W = normalize(Geom::mesh_facet_normal(M,f));
776			vec3 V = cross(W,U);
777
778		✗	frames_[9*f+0] = U.x;
779		✗	frames_[9*f+1] = U.y;
780		✗	frames_[9*f+2] = U.z;
781		✗	frames_[9*f+3] = V.x;
782		✗	frames_[9*f+4] = V.y;
783		✗	frames_[9*f+5] = V.z;
784		✗	frames_[9*f+6] = W.x;
785		✗	frames_[9*f+7] = W.y;
786		✗	frames_[9*f+8] = W.z;
787
788		✗	vec3 g = Geom::mesh_facet_center(M,f);
789		✗	centers_[3*f+0] = g.x;
790		✗	centers_[3*f+1] = g.y;
791		✗	centers_[3*f+2] = g.z;
792			}
793
794		✗	if(use_NN_ || volumetric) {
795		✗	NN_ = NearestNeighborSearch::create(3, "default");
796		✗	NN_->set_points(centers_.size()/3, centers_.data());
797			}
798
799			// Step 4: In volumetric mode, for each tet we find the nearest
800			// facet and lookup the frame field from it.
801		✗	if(volumetric) {
802		✗	vector<double> new_frames(9*M.cells.nb());
803		✗	vector<double> new_centers(3*M.cells.nb());
804		✗	for(index_t t: M.cells) {
805		✗	vec3 g = Geom::mesh_tet_center(M,t);
806		✗	get_nearest_frame(g.data(), &new_frames[9*t]);
807		✗	new_centers[3*t+0] = g.x;
808		✗	new_centers[3*t+1] = g.y;
809		✗	new_centers[3*t+2] = g.z;
810			}
811			frames_.swap(new_frames);
812			centers_.swap(new_centers);
813		✗	NN_->set_points(centers_.size()/3, centers_.data());
814			}
815		✗	}
816
817
818
819
820		✗	void FrameField::scale_frame_vector(
821			double* frame, const vec3& N, double s
822			) {
823			index_t max_index = 0;
824			double max_prod = -1e30;
825		✗	for(index_t i = 0; i < 3; ++i) {
826		✗	double cur_prod =
827			::fabs(
828		✗	N.x * frame[3 * i] +
829		✗	N.y * frame[3 * i + 1] +
830		✗	N.z * frame[3 * i + 2]
831			);
832		✗	if(cur_prod > max_prod) {
833			max_prod = cur_prod;
834			max_index = i;
835			}
836			}
837		✗	vec3 W = s*normalize(vec3(&frame[3*max_index]));
838		✗	frame[3*max_index ]=W.x;
839		✗	frame[3*max_index+1]=W.y;
840		✗	frame[3*max_index+2]=W.z;
841		✗	}
842
843		✗	void FrameField::fix_frame(double* frame, const vec3& N) {
844			index_t w_index=0;
845			double max_prod = -1e30;
846		✗	for(index_t i=0; i<3; ++i) {
847		✗	double cur_prod =
848		✗	::fabs(N.x*frame[3*i]+N.y*frame[3*i+1]+N.z*frame[3*i+2]);
849		✗	if(cur_prod > max_prod) {
850			max_prod = cur_prod;
851			w_index = i;
852			}
853			}
854		✗	index_t u_index = (w_index + 1)%3;
855		✗	index_t v_index = (u_index + 1)%3;
856		✗	vec3 U = normalize(vec3(frame+3*u_index));
857		✗	vec3 V = normalize(vec3(frame+3*v_index));
858		✗	vec3 W = normalize(N);
859
860		✗	if(dot(cross(U,V),W) < 0.0) {
861		✗	U = -U;
862			}
863
864			U -= dot(U,W)*W;
865			V -= dot(V,W)*W;
866
867		✗	U = normalize(U);
868		✗	V = normalize(V);
869
870
871		✗	vec3 X = normalize(U+V);
872		✗	vec3 Y = normalize(U-V);
873		✗	U = normalize(X+Y);
874		✗	V = normalize(X-Y);
875
876		✗	frame[0] = U.x;
877		✗	frame[1] = U.y;
878		✗	frame[2] = U.z;
879		✗	frame[3] = V.x;
880		✗	frame[4] = V.y;
881		✗	frame[5] = V.z;
882		✗	frame[6] = W.x;
883		✗	frame[7] = W.y;
884		✗	frame[8] = W.z;
885		✗	}
886
887
888			}
889