Cleanup of debug outputs

include/Algo/Decimation/edgeSelector.hpp

@@ -2010,9 +2010,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] ;
 
@@ -2027,7 +2024,7 @@ void EdgeSelector_Lightfield<PFP>::computeEdgeInfo(Dart d, EdgeInfo& einfo)
 		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

@@ -207,7 +207,7 @@ void Approximator_HemiFuncCoefs<PFP>::approximate(Dart d)
 			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)
 
 		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/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 ;
-
-	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>
 Geom::Matrix33d
 QuadricHF<REAL>::buildRotateMatrix(const REAL& gamma)