Merge branch 'master' of cgogn:CGoGN

Conflicts:
	include/Geometry/vector_gen.h

include/Algo/Decimation/decimation.h

@@ -62,7 +62,7 @@ void decimate(
 	std::vector<VertexAttribute<typename PFP::VEC3> *>& position,
 	unsigned int nbWantedVertices,
 	const FunctorSelect& selected = allDarts,
-	VertexAttribute<typename PFP::VEC3> *edgeErrors = NULL,
+	EdgeAttribute<typename PFP::REAL> *edgeErrors = NULL,
 	void (*callback_wrapper)(void*, const void*) = NULL, void *callback_object = NULL
 ) ;
 

include/Algo/Decimation/edgeSelector.hpp

@@ -2070,9 +2070,6 @@ void EdgeSelector_Lightfield<PFP>::computeEdgeInfo(Dart d, EdgeInfo& einfo)
 
 	double alpha = alpha1 + alpha2 ;
 
-	if (isnan(alpha))
-		std::cerr << "Nan: " << m_frameN[d] << " ; " << m_frameN[dd] << " ; " << newFN << std::endl ;
-
 	assert(m_quadricHF.isValid() | !"EdgeSelector_Lightfield<PFP>::computeEdgeInfo: quadricHF is not valid") ;
 	Utils::QuadricHF<REAL> quadHF = m_quadricHF[d] ;
 
@@ -2083,11 +2080,11 @@ void EdgeSelector_Lightfield<PFP>::computeEdgeInfo(Dart d, EdgeInfo& einfo)
 	const REAL& errLF = quadHF(newHF) ; // function coefficients
 
 	// Check if errated values appear
-	if (errG < -1e-10 || errAngle < -1e-10 || errLF < -1e-10)
+	if (errG < -1e-6 || errAngle < -1e-6 || errLF < -1e-6)
 		einfo.valid = false ;
 	else
 	{
-		einfo.it = edges.insert(std::make_pair(std::max(errG + errAngle + errLF, REAL(0)), d)) ;
+		einfo.it = edges.insert(std::make_pair(std::max(/*errG +*/ errAngle + errLF, REAL(0)), d)) ;
 		einfo.valid = true ;
 	}
 }

include/Algo/Decimation/halfEdgeSelector.hpp

@@ -939,9 +939,6 @@ void HalfEdgeSelector_Lightfield<PFP>::computeHalfEdgeInfo(Dart d, HalfEdgeInfo&
 
 	double alpha = alpha1 + alpha2 ;
 
-	if (isnan(alpha))
-		std::cerr << "Nan: " << m_frameN[d] << " ; " << m_frameN[dd] << " ; " << newFN << std::endl ;
-
 	assert(m_quadricHF.isValid() | !"EdgeSelector_Lightfield<PFP>::computeEdgeInfo: quadricHF is not valid") ;
 	Utils::QuadricHF<REAL> quadHF = m_quadricHF[d] ;
 

include/Algo/Decimation/lightfieldApproximator.hpp

@@ -204,10 +204,10 @@ void Approximator_HemiFuncCoefs<PFP>::approximate(Dart d)
 		bool opt = m_quadricHF[d].findOptimizedCoefs(coefs) ;
 		// copy result
 		for (unsigned int i = 0 ; i < m_nbCoefs ; ++i)
-			this->m_approx[i][d] = opt ? coefs[i] : (m_coefs[i]->operator[](d) + m_coefs[i]->operator[](dd))/2 ; // if fail average coefs (TODO better)
+			this->m_approx[i][d] = coefs[i] ;
 
 		if (!opt)
-			std::cerr << "LightfieldApproximator::Inversion failed: not treated correctly" << std::endl ;
+			std::cerr << "LightfieldApproximator::Inversion failed: coefs are averaged" << std::endl ;
 	}
 }
 
@@ -337,8 +337,8 @@ void Approximator_HemiFuncCoefsHalfEdge<PFP>::approximate(Dart d)
 	const REAL gamma1 = ((B1 * T) > 0 ? 1 : -1) * acos( std::max(std::min(1.0f, T1 * T ), -1.0f)) ; // angle positif ssi
 	const REAL gamma2 = ((B2 * T) > 0 ? 1 : -1) * acos( std::max(std::min(1.0f, T2 * T ), -1.0f)) ; // -PI/2 < angle(i,j1) < PI/2  ssi i*j1 > 0
 	// Rotation dans le sens trigo autour de l'axe T (N1 --> N)
-	const REAL alpha1 = ((N * B1prime) > 0 ? -1 : 1) * acos( std::max(std::min(1.0f, N * N1), -1.0f) ) ; // angle positif ssi
-	const REAL alpha2 = ((N * B2prime) > 0 ? -1 : 1) * acos( std::max(std::min(1.0f, N * N2), -1.0f) ) ; // PI/2 < angle(j1',n) < -PI/2 ssi j1'*n < 0
+	const REAL alpha1 = (N == N1) ? 0 : ((N * B1prime) > 0 ? -1 : 1) * acos( std::max(std::min(1.0f, N * N1), -1.0f) ) ; // angle positif ssi
+	const REAL alpha2 = (N == N2) ? 0 : ((N * B2prime) > 0 ? -1 : 1) * acos( std::max(std::min(1.0f, N * N2), -1.0f) ) ; // PI/2 < angle(j1',n) < -PI/2 ssi j1'*n < 0
 
 	assert (m_quadricHF.isValid() || !"LightfieldApproximator:approximate quadricHF is not valid") ;
 
@@ -383,10 +383,13 @@ void Approximator_HemiFuncCoefsHalfEdge<PFP>::approximate(Dart d)
 		bool opt = m_quadricHF[d].findOptimizedCoefs(coefs) ;
 		// copy result
 		for (unsigned int i = 0 ; i < m_nbCoefs ; ++i)
-			this->m_approx[i][d] = opt ? coefs[i] : (m_coefs[i]->operator[](d) + m_coefs[i]->operator[](dd))/2 ; // if fail average coefs (TODO better)
+			this->m_approx[i][d] = coefs[i] ;
 
 		if (!opt)
-			std::cerr << "LightfieldApproximator::Inversion failed: not treated correctly" << std::endl ;
+		{
+			std::cerr << "LightfieldApproximatorHalfCollapse::Optimization failed (should never happen since no optim is done)" << std::endl ;
+			std::cout << alpha1 << std::endl ;
+		}
 
 		// Add second one for error evaluation
 		m_quadricHF[d] += Utils::QuadricHF<REAL>(coefs2, gamma2, alpha2) ;

include/Algo/Geometry/curvature.h

@@ -31,6 +31,9 @@
 #include "OpenNL/sparse_matrix.h"
 #include "OpenNL/full_vector.h"
 
+#include <Eigen/Core>
+#include <Eigen/Eigenvalues>
+
 namespace CGoGN
 {
 

include/Algo/Geometry/curvature.hpp

@@ -313,6 +313,7 @@ void computeCurvatureVertex_NormalCycles(
 	typedef typename PFP::VEC3 VEC3 ;
 	typedef typename PFP::REAL REAL ;
 
+	// collect the normal cycle tensor
 	Algo::Selection::Collector_WithinSphere<PFP> neigh(map, position, radius, thread) ;
 	neigh.collectAll(dart) ;
 	neigh.computeArea() ;
@@ -321,14 +322,14 @@ void computeCurvatureVertex_NormalCycles(
 
 	typename PFP::MATRIX33 tensor(0) ;
 
-	// inside
+	// collect edges inside the neighborhood
 	const std::vector<Dart>& vd1 = neigh.getInsideEdges() ;
 	for (std::vector<Dart>::const_iterator it = vd1.begin(); it != vd1.end(); ++it)
 	{
 		const VEC3 e = Algo::Geometry::vectorOutOfDart<PFP>(map, *it, position) ;
 		tensor += Geom::transposed_vectors_mult(e,e) * edgeangle[*it] * (1 / e.norm()) ;
 	}
-	// border
+	// collect edges crossing the neighborhood's border
 	const std::vector<Dart>& vd2 = neigh.getBorder() ;
 	for (std::vector<Dart>::const_iterator it = vd2.begin(); it != vd2.end(); ++it)
 	{
@@ -340,43 +341,40 @@ void computeCurvatureVertex_NormalCycles(
 
 	tensor /= neigh.getArea() ;
 
-	long int n = 3, lda = 3, info, lwork = 9 ;
-	char jobz='V', uplo = 'U' ;
-	float work[9] ;
-	float w[3] ;
-	float a[3*3] = {
-		tensor(0,0), 0.0f, 0.0f,
-		tensor(1,0), tensor(1,1), 0.0f,
-		tensor(2,0), tensor(2,1), tensor(2,2)
-	} ;
-	// Solve eigenproblem
-	ssyev_(&jobz, &uplo, (integer*)&n, a, (integer*)&lda, w, work, (integer*)&lwork, (integer*)&info) ;
-	// Check for convergence
-	if(info > 0)
-		std::cerr << "clapack ssyev_ failed to compute eigenvalues : exit with code " << info << std::endl ;
-	// sort eigen components : w[s[0]] has minimal absolute value ; kmin = w[s[1]] <= w[s[2]] = kmax
+	// solve eigen problem
+	Eigen::Matrix3f m3 ;
+	m3 << tensor(0,0) , tensor(0,1) , tensor(0,2) , tensor(1,0) , tensor(1,1) , tensor(1,2) , tensor(2,0) , tensor(2,1) , tensor(2,2) ;
+	Eigen::SelfAdjointEigenSolver<Eigen::Matrix3f> solver (m3);
+	Eigen::Vector3f ev = solver.eigenvalues();
+	Eigen::Matrix3f evec = solver.eigenvectors();
+
+	// sort eigen components : ev[s[0]] has minimal absolute value ; kmin = ev[s[1]] <= ev[s[2]] = kmax
 	int s[3] = {0, 1, 2} ;
 	int tmp ;
-	if (abs(w[s[2]]) < abs(w[s[1]])) { tmp = s[1] ; s[1] = s[2] ; s[2] = tmp ; }
-	if (abs(w[s[1]]) < abs(w[s[0]])) { tmp = s[0] ; s[0] = s[1] ; s[1] = tmp ; }
-	if (w[s[2]] < w[s[1]]) { tmp = s[1] ; s[1] = s[2] ; s[2] = tmp ; }
+	if (abs(ev[s[2]]) < abs(ev[s[1]])) { tmp = s[1] ; s[1] = s[2] ; s[2] = tmp ; }
+	if (abs(ev[s[1]]) < abs(ev[s[0]])) { tmp = s[0] ; s[0] = s[1] ; s[1] = tmp ; }
+	if (ev[s[2]] < ev[s[1]]) { tmp = s[1] ; s[1] = s[2] ; s[2] = tmp ; }
 
-	kmin[dart] = w[s[1]] ;
-	kmax[dart] = w[s[2]] ;
+	// set curvatures from sorted eigen components
+	// warning : Kmin and Kmax are switched w.r.t. kmin and kmax
+	// normal direction : minimal absolute eigen value
+	VEC3& dirNormal = Knormal[dart] ;
+	dirNormal[0] = evec(0,s[0]);
+	dirNormal[1] = evec(1,s[0]);
+	dirNormal[2] = evec(2,s[0]);
+	if (dirNormal * normal[dart] < 0) dirNormal *= -1; // change orientation
+	// min curvature
+	kmin[dart] = ev[s[1]] ;
 	VEC3& dirMin = Kmin[dart] ;
-	dirMin[0] = a[3*s[2]];
-	dirMin[1] = a[3*s[2]+1];
-	dirMin[2] = a[3*s[2]+2]; // warning : Kmin and Kmax are switched
+	dirMin[0] = evec(0,s[2]);
+	dirMin[1] = evec(1,s[2]);
+	dirMin[2] = evec(2,s[2]);
+	// max curvature
+	kmax[dart] = ev[s[2]] ;
 	VEC3& dirMax = Kmax[dart] ;
-	dirMax[0] = a[3*s[1]];
-	dirMax[1] = a[3*s[1]+1];
-	dirMax[2] = a[3*s[1]+2]; // warning : Kmin and Kmax are switched
-	VEC3& dirNormal = Knormal[dart] ;
-	dirNormal[0] = a[3*s[0]];
-	dirNormal[1] = a[3*s[0]+1];
-	dirNormal[2] = a[3*s[0]+2];
-	if (dirNormal * normal[dart] < 0)
-		dirNormal *= -1; // change orientation
+	dirMax[0] = evec(0,s[1]);
+	dirMax[1] = evec(1,s[1]);
+	dirMax[2] = evec(2,s[1]);
 }
 
 

include/Algo/Import/import2tablesSurface.hpp

@@ -1010,12 +1010,21 @@ bool MeshTablesSurface<PFP>::importPlySLFgenericBin(const std::string& filename,
 	attrNames.push_back(positions.name()) ;
 
 	VertexAttribute<typename PFP::VEC3> *frame = new VertexAttribute<typename PFP::VEC3>[3] ;
-	frame[0] = m_map.template addAttribute<typename PFP::VEC3, VERTEX>("frameT") ; // Tangent
-	frame[1] = m_map.template addAttribute<typename PFP::VEC3, VERTEX>("frameB") ; // Bitangent
-	frame[2] = m_map.template addAttribute<typename PFP::VEC3, VERTEX>("frameN") ; // Normal
-	attrNames.push_back(frame[0].name()) ;
-	attrNames.push_back(frame[1].name()) ;
-	attrNames.push_back(frame[2].name()) ;
+	if (tangent)
+	{
+		frame[0] = m_map.template addAttribute<typename PFP::VEC3, VERTEX>("frameT") ; // Tangent
+		attrNames.push_back(frame[0].name()) ;
+	}
+	if (binormal)
+	{
+		frame[1] = m_map.template addAttribute<typename PFP::VEC3, VERTEX>("frameB") ; // Bitangent
+		attrNames.push_back(frame[0].name()) ;
+	}
+	if (normal)
+	{
+		frame[2] = m_map.template addAttribute<typename PFP::VEC3, VERTEX>("frameN") ; // Normal
+		attrNames.push_back(frame[0].name()) ;
+	}
 
 	VertexAttribute<typename PFP::VEC3> *PBcoefs = NULL, *SHcoefs = NULL ;
 	if (PTM)
@@ -1064,14 +1073,15 @@ bool MeshTablesSurface<PFP>::importPlySLFgenericBin(const std::string& filename,
 
 		fp.read((char*)properties,nbProps * propSize) ;
 
-		positions[id] = VEC3(properties[0],properties[1],properties[2]) ; // position
-		for (unsigned int k = 0 ; k < 3 ; ++k) // frame
-			for (unsigned int l = 0 ; l < 3 ; ++l)
-				frame[k][id][l] = (typename PFP::REAL)(properties[3+(3*k+l)]) ;
-		/*
-		 	for (unsigned int k = 0 ; k < 3 ; ++k) // coefficients
-			for (unsigned int l = 0 ; l < nbCoefs ; ++l)
-		 */
+		// positions
+		if (nbProps > 2)
+			positions[id] = VEC3(properties[0],properties[1],properties[2]) ; // position
+
+		if (tangent && binormal && normal) // == if (nbprops > 11)
+			for (unsigned int k = 0 ; k < 3 ; ++k) // frame
+				for (unsigned int l = 0 ; l < 3 ; ++l)
+					frame[k][id][l] = (typename PFP::REAL)(properties[3+(3*k+l)]) ;
+
 		for (unsigned int l = 0 ; l < nbCoefs ; ++l) // coefficients
 			for (unsigned int k = 0 ; k < 3 ; ++k)
 			{
@@ -1080,6 +1090,7 @@ bool MeshTablesSurface<PFP>::importPlySLFgenericBin(const std::string& filename,
 				else /* if SH */
 					SHcoefs[l][id][k] = (typename PFP::REAL)(properties[12+(3*l+k)]) ;
 			}
+
 		unsigned int cur = 12+3*nbCoefs ;
 		for (unsigned int k = 0 ; k < nbRemainders ; ++k) // remaining data
 			remainders[k][id] = (typename PFP::REAL)(properties[cur + k]) ;

include/Algo/ProgressiveMesh/pmesh.h

@@ -65,7 +65,7 @@ private:
 	std::vector<VSplit<PFP>*> m_splits ;
 	unsigned int m_cur ;
 
-	Algo::Decimation::HalfEdgeApproximator<PFP, VEC3, EDGE>* m_positionApproximator ;
+	Algo::Decimation::Approximator<PFP, VEC3, EDGE>* m_positionApproximator ;
 
 	bool m_initOk ;
 

include/Algo/ProgressiveMesh/pmesh.hpp

@@ -108,7 +108,7 @@ ProgressiveMesh<PFP>::ProgressiveMesh(
 		if(! (*it)->init())
 			m_initOk = false ;
 		if((*it)->getApproximatedAttributeName() == "position")
-			m_positionApproximator = reinterpret_cast<Algo::Decimation::HalfEdgeApproximator<PFP, VEC3, EDGE>*>(*it) ;
+			m_positionApproximator = reinterpret_cast<Algo::Decimation::Approximator<PFP, VEC3, EDGE>*>(*it) ;
 	}
 	CGoGNout << "..done" << CGoGNendl ;
 

include/Utils/pointSprite.frag

@@ -6,42 +6,44 @@ uniform vec3 ambiant;
 uniform float size;
 #ifdef WITH_PLANE
 	uniform vec3 eyePos;
-	uniform vec3 eyeY;
-	VARYING_FRAG vec3 eyePosFrag;
+	VARYING_FRAG vec3 shiftedEye;
 #endif
 
-
 #ifdef WITH_COLOR_PER_VERTEX 
 	VARYING_FRAG vec3 colorsprite;
 #else
 	uniform vec3 colorsprite;
 #endif
 
-
-VARYING_FRAG vec2 texCoord;
+VARYING_FRAG vec2 spriteCoord;
 VARYING_FRAG vec3 sphereCenter;
 
-
-
 void main(void)
 {
-
-	vec3 billboard_frag_pos = sphereCenter + vec3(texCoord, 0.0) * size;
 #ifdef WITH_PLANE
-	vec3 ray_direction = normalize(billboard_frag_pos - eyePosFrag;);
+	vec3 billboard_frag_pos = vec3(spriteCoord, 0.0) * size;
+	vec3 ray_direction = normalize(billboard_frag_pos - shiftedEye);
+	float av = dot(shiftedEye,ray_direction);
+	float arg = av*av - dot(shiftedEye,shiftedEye) + size*size;
+	if (arg< 0.0)
+		discard;
+	float t = -av - sqrt(arg);
+	vec3 frag_position_eye = ray_direction * t + eyePos	;
 #else
+	vec3 billboard_frag_pos = sphereCenter + vec3(spriteCoord, 0.0) * size;
 	vec3 ray_direction = normalize(billboard_frag_pos);
-#endif
-
 	float TD = -dot(ray_direction,sphereCenter);
 	float c = dot(sphereCenter, sphereCenter) - size * size;
 	float arg = TD * TD - c;
-      
+    
 	if (arg < 0.0)
 		discard;
-  
+
 	float t = -c / (TD - sqrt(arg));
-	vec3 frag_position_eye = ray_direction * t;
+	vec3 frag_position_eye = ray_direction * t ;
+#endif	
+
+
 	vec4 pos = ProjectionMatrix * vec4(frag_position_eye, 1.0);
 	gl_FragDepth = (pos.z / pos.w + 1.0) / 2.0;
 

include/Utils/pointSprite.geom

@@ -5,12 +5,12 @@ uniform mat4 ModelViewMatrix;
 uniform mat4 ProjectionMatrix;
 #ifdef WITH_PLANE
 	uniform vec3 eyePos;
-	uniform vec3 eyeY;
-	VARYING_OUT vec4 eyePosFrag;
+	VARYING_OUT vec3 shiftedEye;
 #endif
 
-VARYING_OUT vec2 texCoord;
+VARYING_OUT vec2 spriteCoord;
 VARYING_OUT vec3 sphereCenter;
+
 #ifdef WITH_COLOR_PER_VERTEX 
 	VARYING_IN vec3 color[1];
 	VARYING_OUT vec3 colorsprite;
@@ -19,16 +19,16 @@ VARYING_OUT vec3 sphereCenter;
 #ifdef WITH_PLANE
 	void corner( vec4 center, vec3 planeX, vec3 planeY, float x, float y)
 	{
-		texCoord = vec2(1.4*x,1.4*y);
-		vec4 pos = center + size( x*vec4(planeX,0.0) + y*vec4(planeY,0.0)+ vec4(0.0,0.0,0.5,0.0)); 
+		spriteCoord = vec2(x,y);
+		vec4 pos = center + size*( x*vec4(planeX,0.0) + y*vec4(planeY,0.0)+ vec4(0.0,0.0,0.5,0.0)); 
 		gl_Position = ProjectionMatrix *  pos;
 		EmitVertex();
 	}
 #else
 	void corner( vec4 center, float x, float y)
 	{
-		texCoord = vec2(1.4*x,1.4*y);
-		vec4 pos = center + vec4( size*x, size*y, 0.0, 0.0);
+		spriteCoord = vec2(x,y);
+		vec4 pos = center + vec4(size*x, size*y, 0.0, 0.0);
 		gl_Position = ProjectionMatrix *  pos;
 		EmitVertex();
 	}
@@ -42,28 +42,27 @@ void main()
 #endif
 
 	vec4 posCenter = ModelViewMatrix * POSITION_IN(0);
-
 	sphereCenter = posCenter.xyz;
-
+	
 #ifdef WITH_PLANE
-	vec4 EPF = ModelViewMatrix * vec4(eyePos,1.0);
-	eyePosFrag = EPF.xyz;
+	shiftedEye = eyePos - sphereCenter;
+	vec3 V = -shiftedEye;
+	normalize(V);
 
-	vec3 V = sphereCenter-eyePosFrag;
-	V.normalize();
-	vec3 planeX = cross(V, eyeY);
-	vec3 planeY = cross(X,V);
-	corner(posCenter, planeX, planeY, -1.0, 1.0);
-	corner(posCenter, planeX, planeY, -1.0,-1.0);
-	corner(posCenter, planeX, planeY,  1.0, 1.0);
-	corner(posCenter, planeX, planeY,  1.0,-1.0);
+	vec3 planeX = vec3(-V[2],0.0,V[0]); //cross(V, vec3(0.0,1.0,0.0));
+	normalize(planeX);
+	vec3 planeY = cross(planeX,V);
+
+	corner(posCenter, planeX, planeY, -1.4, 1.4);
+	corner(posCenter, planeX, planeY, -1.4,-1.4);
+	corner(posCenter, planeX, planeY,  1.4, 1.4);
+	corner(posCenter, planeX, planeY,  1.4,-1.4);
 #else
-	corner(posCenter, -1.0, 1.0);
-	corner(posCenter, -1.0,-1.0);
-	corner(posCenter,  1.0, 1.0);
-	corner(posCenter,  1.0,-1.0);
+	corner(posCenter, -1.4, 1.4);
+	corner(posCenter, -1.4,-1.4);
+	corner(posCenter,  1.4, 1.4);
+	corner(posCenter,  1.4,-1.4);
 #endif	
-	
-	
 	EndPrimitive();
 }
+

include/Utils/pointSprite.h

@@ -64,7 +64,7 @@ protected:
 	
 	CGoGNGLuint m_uniform_EyePos;
 	
-	CGoGNGLuint m_uniform_EyeY;
+//	CGoGNGLuint m_uniform_EyeY;
 
 	CGoGNGLuint m_uniform_ambiant;
 
@@ -109,7 +109,7 @@ public:
 	/**
 	* set the plane of rendering for VR rendering
 	*/
-	void setEyePosition(const Geom::Vec3f& ox, const Geom::Vec3f& oy);
+	void setEyePosition(const Geom::Vec3f& ep);
 
 
 	/**

include/Utils/qem.h

@@ -583,6 +583,8 @@ private:
 	Eigen::MatrixXd m_A ; /*!< The first QuadricHF member matrix A */
 	Eigen::VectorXd m_b[3] ; /*!< The second QuadricHF member vector b */
 	double m_c[3] ; /*!< The third QuadricHF member scalar c */
+	std::vector<VEC3> m_coefs ; /*!< The coefficients in cas optim fails */
+	bool m_noAlphaRot ; /*!< If alpha = 0 then optim will fail */
 
 	/*!
 	 * \brief method to evaluate the error for a given lightfield function
@@ -593,18 +595,6 @@ private:
 	 */
 	REAL evaluate(const std::vector<VEC3>& coefs) const ;
 
-	/*!
-	 * \brief method to deduce an optimal coefficients in space
-	 * w.r.t. the current QuadricHF.
-	 *
-	 * \param coefs will contain the optimal result (if it can be computed)
-	 *
-	 * \return true if an optimal result was correctly computed
-	 */
-	bool optimize(std::vector<VEC3>& coefs) const ;
-
-//	Geom::Tensor3d rotate(const Geom::Tensor3d& T, const Geom::Matrix33d& R) ;
-
 	/*!
 	 * \brief method to build a rotate matrix (rotation in tangent plane)
 	 * given angle gamma

include/Utils/qem.hpp

@@ -300,11 +300,13 @@ QuadricNd<REAL,N>::optimize(VECN& v) const
 }
 
 template <typename REAL>
-QuadricHF<REAL>::QuadricHF()
+QuadricHF<REAL>::QuadricHF():
+m_noAlphaRot(false)
 {}
 
 template <typename REAL>
-QuadricHF<REAL>::QuadricHF(int i)
+QuadricHF<REAL>::QuadricHF(int i):
+m_noAlphaRot(false)
 {
 	m_A.resize(i,i) ;
 	for (unsigned int c = 0 ; c < 3 ; ++c)
@@ -323,18 +325,17 @@ QuadricHF<REAL>::QuadricHF(const std::vector<VEC3>& v, const REAL& gamma, const
 }
 
 template <typename REAL>
-QuadricHF<REAL>::QuadricHF(const Geom::Tensor3d* T, const REAL& gamma, const REAL& alpha)
+QuadricHF<REAL>::QuadricHF(const Geom::Tensor3d* T, const REAL& gamma, const REAL& alpha):
+m_noAlphaRot(fabs(alpha) < 1e-13)
 {
 	const unsigned int nbcoefs = ((T[0].order() + 1) * (T[0].order() + 2)) / 2. ;
 
-	//m_A.resize(nbcoefs, nbcoefs) ;
-
 	// 2D rotation
 	const Geom::Matrix33d R = buildRotateMatrix(gamma) ;
 	Geom::Tensor3d* Trot = new Geom::Tensor3d[3] ;
 	for (unsigned int c = 0 ; c < 3 ; ++c)
 		Trot[c] = rotate(T[c],R) ;
-	std::vector<VEC3> coefsR = coefsFromTensors(Trot) ;
+	m_coefs = coefsFromTensors(Trot) ;
 	delete[] Trot ;
 
 	// parameterized integral on intersection
@@ -348,7 +349,7 @@ QuadricHF<REAL>::QuadricHF(const Geom::Tensor3d* T, const REAL& gamma, const REA
 	{
 		Eigen::VectorXd v ;	v.resize(nbcoefs) ;
 		for (unsigned int i = 0 ; i < nbcoefs ; ++i) // copy into vector
-			v[i] = coefsR[i][c] ;
+			v[i] = m_coefs[i][c] ;
 
 		m_b[c] = integ_b * v ; // Vector b
 		m_c[c] = v.transpose() * (integ_c * v) ; // Constant c
@@ -369,6 +370,8 @@ QuadricHF<REAL>::zero()
 		m_b[c].setZero() ;
 		m_c[c] = 0 ;
 	}
+	m_coefs.clear() ;
+	m_noAlphaRot = false ;
 }
 
 template <typename REAL>
@@ -381,6 +384,8 @@ QuadricHF<REAL>::operator= (const QuadricHF<REAL>& q)
 		m_b[c] = q.m_b[c] ;
 		m_c[c] = q.m_c[c] ;
 	}
+	m_coefs = q.m_coefs ;
+	m_noAlphaRot = q.m_noAlphaRot ;
 
 	return *this ;
 }
@@ -399,6 +404,11 @@ QuadricHF<REAL>::operator+= (const QuadricHF<REAL>& q)
 		m_b[c] += q.m_b[c] ;
 		m_c[c] += q.m_c[c] ;
 	}
+	m_coefs.resize(m_coefs.size()) ;
+	for (unsigned int i = 0 ; i < m_coefs.size() ; ++i)
+		m_coefs[i] += q.m_coefs[i] ;
+
+	m_noAlphaRot &= q.m_noAlphaRot ;
 
 	return *this ;
 }
@@ -418,6 +428,12 @@ QuadricHF<REAL>::operator -= (const QuadricHF<REAL>& q)
 		m_c[c] -= q.m_c[c] ;
 	}
 
+	m_coefs.resize(m_coefs.size()) ;
+	for (unsigned int i = 0 ; i < m_coefs.size() ; ++i)
+		m_coefs[i] -= q.m_coefs[i] ;
+
+	m_noAlphaRot &= q.m_noAlphaRot ;
+
 	return *this ;
 }
 
@@ -425,6 +441,7 @@ template <typename REAL>
 QuadricHF<REAL>&
 QuadricHF<REAL>::operator *= (const REAL& v)
 {
+	std::cout << "Warning: QuadricHF<REAL>::operator *= should not be used !" << std::endl ;
 	m_A *= v ;
 	for (unsigned int c = 0 ; c < 3 ; ++c)
 	{
@@ -439,6 +456,7 @@ template <typename REAL>
 QuadricHF<REAL>&
 QuadricHF<REAL>::operator /= (const REAL& v)
 {
+	std::cout << "Warning: QuadricHF<REAL>::operator /= should not be used !" << std::endl ;
 	const REAL& inv = 1. / v ;
 
 	(*this) *= inv ;
@@ -457,7 +475,24 @@ template <typename REAL>
 bool
 QuadricHF<REAL>::findOptimizedCoefs(std::vector<VEC3>& coefs)
 {
-	return optimize(coefs) ;
+	coefs.resize(m_b[0].size()) ;
+
+	if (fabs(m_A.determinant()) < 1e-10 )
+	{
+		coefs = m_coefs ;
+		return m_noAlphaRot ; // if true inversion failed (normal!) and m_coefs forms a valid solution
+	}
+	Eigen::MatrixXd Ainv = m_A.inverse() ;
+
+	for (unsigned int c = 0 ; c < 3 ; ++c)
+	{
+		Eigen::VectorXd tmp(m_b[0].size()) ;
+		tmp = Ainv * m_b[c] ;
+		for (unsigned int i = 0 ; i < m_b[c].size() ; ++i)
+			coefs[i][c] = tmp[i] ;
+	}
+
+	return true ;
 }
 
 template <typename REAL>
@@ -480,28 +515,6 @@ QuadricHF<REAL>::evaluate(const std::vector<VEC3>& coefs) const
 	return (res[0] + res[1] + res[2]) / 3. ;
 }
 
-template <typename REAL>
-bool
-QuadricHF<REAL>::optimize(std::vector<VEC3>& coefs) const
-{
-	coefs.resize(m_b[0].size()) ;
-
-	if (fabs(m_A.determinant()) < 1e-10 )
-		return false ;