Decimation for lightfield data v1

include/Algo/Decimation/edgeSelector.h

@@ -544,12 +544,14 @@ private:
 	EdgeAttribute<EdgeInfo> edgeInfo ;
 
 	VertexAttribute<VEC3> m_pos, m_frameT, m_frameB, m_frameN ;
-	//VertexAttribute<VEC3> *m_HF ;
+	std::vector<VertexAttribute<VEC3> > m_HF ;
 	int m_approxindex_pos, m_attrindex_pos ;
 	int m_approxindex_FN, m_attrindex_FN ;
+	std::vector<unsigned int> m_approxindex_HF, m_attrindex_HF ;
+	unsigned int m_K ;
 
 	VertexAttribute<Quadric<REAL> > m_quadricGeom ;
-	VertexAttribute<QuadricHF<REAL> > m_quadricHF ;
+	EdgeAttribute<QuadricHF<REAL> > m_quadricHF ;
 
 	std::vector<Approximator<PFP, typename PFP::VEC3>* > m_approx ;
 
@@ -567,11 +569,12 @@ public:
 		m_approxindex_pos(-1),
 		m_attrindex_pos(-1),
 		m_approxindex_FN(-1),
-		m_attrindex_FN(-1)
+		m_attrindex_FN(-1),
+		m_K(0)
 	{
 		edgeInfo = m.template addAttribute<EdgeInfo, EDGE>("edgeInfo") ;
 		m_quadricGeom = m.template addAttribute<Quadric<REAL>, VERTEX>("QEMquadric") ;
-		m_quadricHF = m.template addAttribute<QuadricHF<REAL>, VERTEX>("HFquadric") ;
+		m_quadricHF = m.template getAttribute<QuadricHF<REAL>, EDGE>("HFquadric") ;
 	}
 	~EdgeSelector_Lightfield()
 	{

include/Algo/Decimation/edgeSelector.hpp

@@ -1686,6 +1686,7 @@ bool EdgeSelector_Lightfield<PFP>::init()
 	unsigned int ok = 0 ;
 	for (unsigned int approxindex = 0 ; approxindex < this->m_approximators.size() ; ++approxindex)
 	{
+		unsigned int k = 0 ;
 		bool saved = false ;
 		for (unsigned int attrindex = 0 ; attrindex < this->m_approximators[approxindex]->getNbApproximated() ; ++ attrindex)
 		{
@@ -1741,10 +1742,31 @@ bool EdgeSelector_Lightfield<PFP>::init()
 					saved = true ;
 				}
 			}
+			else
+			{
+				std::stringstream s ;
+				s << "PBcoefs" << k ;
+				if(ok > 3 && this->m_approximators[approxindex]->getApproximatedAttributeName(attrindex) == s.str().c_str())
+				{
+					++ok ;
+					m_HF.push_back(m.template getAttribute<typename PFP::VEC3, VERTEX>(s.str().c_str())) ;
+					if (m_HF[k++].isValid())
+					{
+						m_approxindex_HF.push_back(approxindex) ;
+						m_attrindex_HF.push_back(attrindex) ;
+						if (!saved)
+						{
+							m_approx.push_back(reinterpret_cast<Approximator<PFP, VEC3>* >(this->m_approximators[approxindex])) ;
+							saved = true ;
+						}
+					}
+				}
+			}
 		}
+		m_K = k ;
 	}
 
-	if(ok != 4)
+	if(ok < 5)
 		return false ;
 
 	TraversorV<MAP> travV(m);
@@ -1829,7 +1851,7 @@ void EdgeSelector_Lightfield<PFP>::recomputeQuadric(const Dart d, const bool rec
 {
 	Dart dFront,dBack ;
 	Dart dInit = d ;
-//todo
+
 	// Init Front
 	dFront = dInit ;
 
@@ -1855,7 +1877,7 @@ template <typename PFP>
 void EdgeSelector_Lightfield<PFP>::updateAfterCollapse(Dart d2, Dart dd2)
 {
 	MAP& m = this->m_map ;
-	//todo
+
 	recomputeQuadric(d2, true) ;
 
 	Dart vit = d2 ;
@@ -1952,11 +1974,10 @@ void EdgeSelector_Lightfield<PFP>::computeEdgeInfo(Dart d, EdgeInfo& einfo)
 	const VEC3& newFN = this->m_approx[m_approxindex_FN]->getApprox(d,m_attrindex_FN) ; // get new frameN
 
 	// New function
-	// todo const VEC3& hfcoefs0 = this->m_approximators[m_approxindex_hf]->getApprox(d,m_attrindex_hf[0]) ; // get newHFcoefs0
-	// todo const VEC3& hfcoefs1 = this->m_approximators[m_approxindex_hf]->getApprox(d,m_attrindex_hf[1]) ; // get newHFcoefs1
-	// todo const VEC3& hfcoefs2 = this->m_approximators[m_approxindex_hf]->getApprox(d,m_attrindex_hf[2]) ; // get newHFcoefs2
-	// ...
-	// todo const VEC3& hfcoefsK = this->m_approximators[m_approxindex_hf]->getApprox(d,m_attrindex_hf[K]) ; // get newHFcoefsK
+	std::vector<VEC3> newHF ;
+	newHF.resize(m_K) ;
+	for (unsigned int k = 0 ; k < m_K ; ++k)
+		newHF[k] = this->m_approx[m_approxindex_HF[k]]->getApprox(d,m_attrindex_HF[k]) ; // get newHFcoefsK
 
 	// Compute errors
 	// Position
@@ -1965,21 +1986,26 @@ void EdgeSelector_Lightfield<PFP>::computeEdgeInfo(Dart d, EdgeInfo& einfo)
 	quadGeom += m_quadricGeom[dd] ;	// two vertices quadrics
 
 	// hemisphere difference error
-	double scal1 = abs(double(m_frameN[d] * newFN)) ;
-	scal1 = std::min(scal1, double(1)) ; // for epsilon normalization of newFN errors
+	double scal1 = double(m_frameN[d] * newFN) ;
+	double alpha1 = acos(std::max(std::min(scal1, double(1)),double(-1))) ; // for epsilon normalization of newFN errors
 	// angle 2
-	double scal2 = abs(double(m_frameN[dd] * newFN)) ;
-	scal2 = std::min(scal2, double(1)) ;
-	// Important: sum of abs values works only if newFN is in-between the two old ones (interpolated).
-	// This avoids computation of the sign of alpha1 and alpha2.
-	double alpha = acos(scal1 + scal2) ;
+	double scal2 = double(m_frameN[dd] * newFN) ;
+	double alpha2 = acos(std::max(std::min(scal2, double(1)),double(-1))) ; // for epsilon normalization of newFN errors
+
+	double alpha = alpha1 + alpha2 ;
+
+	if (isnan(alpha))
+		std::cerr << "Nan" << std::endl ;
+
+	assert(m_quadricHF.isValid() | !"EdgeSelector_Lightfield<PFP>::computeEdgeInfo: quadricHF is not valid") ;
+	QuadricHF<REAL> quadHF = m_quadricHF[d] ;
 
-	std::cout << quadGeom(newPos) << " vs " << alpha/M_PI << std::endl ;
+	//std::cout << quadGeom(newPos) << " vs " << alpha/M_PI << " vs " << quadHF(newHF) << std::endl ;
 	// sum of QEM metric and frame orientation difference
 	const REAL& err =
 			quadGeom(newPos) // geom
 			+ alpha / M_PI // frame
-			// todo+ quadHF(newHF) // function coefficients
+			+ quadHF(newHF) // function coefficients
 			 ;
 
 	// Check if errated values appear

include/Algo/Decimation/lightfieldApproximator.h

@@ -115,7 +115,7 @@ public:
 
 	std::vector<VertexAttribute<VEC3>* > m_coefs ;
 
-	VertexAttribute<QuadricHF<REAL> > m_quadricHF ;
+	EdgeAttribute<QuadricHF<REAL> > m_quadricHF ;
 
 public:
 	Approximator_HemiFuncCoefs(MAP& m, std::vector<VertexAttribute<VEC3>* >& attr, Predictor<PFP, VEC3>* pred = NULL) :
@@ -123,12 +123,8 @@ public:
 		m_nbCoefs(0),
 		m_HFtype(0) // SH = 0
 	{
-		// check name of number 0
-		if (this->m_attrV[0]->name().find("SH") != std::string::npos)
-			m_HFtype = 1 ;
-
 		unsigned int i ;
-		for (i = 1 ; i < 200 ; ++i)
+		for (i = 0 ; i < 200 ; ++i)
 		{
 			// check if number i is present
 			if ((this->m_attrV.size() <= i) || this->m_attrV[i]->name().find("coefs") == std::string::npos)
@@ -136,7 +132,15 @@ public:
 
 			m_coefs.push_back(this->m_attrV[i]) ;
 		}
-		m_nbCoefs = i - 1 ;
+
+		// check name of last valid
+		if (this->m_attrV[i-1]->name().find("PB") != std::string::npos)
+			m_HFtype = 1 ;
+
+		m_nbCoefs = i ;
+
+		// set quadric
+		m_quadricHF = this->m_map.template addAttribute<QuadricHF<REAL>, EDGE>("HFquadric") ;
 	}
 	~Approximator_HemiFuncCoefs()
 	{}

include/Algo/Decimation/lightfieldApproximator.hpp

@@ -45,11 +45,13 @@ void Approximator_FrameInterpolation<PFP>::approximate(Dart d)
 	{
 		for (unsigned int i = 0 ; i < 3 ; ++i)
 			this->m_approx[i][d] = this->m_attrV[i]->operator[](d) ;
+		//std::cout << "fallback to p1" << std::endl ;
 	}
 	else if (this->m_approxposition[d] == this->m_position[dd]) // new Position is position of vertex dd
 	{
 		for (unsigned int i = 0 ; i < 3 ; ++i)
 			this->m_approx[i][d] = this->m_attrV[i]->operator[](dd) ;
+		//std::cout << "fallback to p2" << std::endl ;
 	}
 	else
 	{
@@ -64,11 +66,15 @@ void Approximator_FrameInterpolation<PFP>::approximate(Dart d)
 		REAL t = (v0v1 * v0v) / v0v1.norm() ;
 		t = std::max (std::min (t , REAL(1)) , REAL(0) ) ; // clamp it to [0,1]
 
-		const VEC3& normal1 = this->m_attrV[2]->operator[](d) ;
-		const VEC3& normal2 = this->m_attrV[2]->operator[](dd) ;
+		VEC3& normal1 = this->m_attrV[2]->operator[](d) ;
+		VEC3& normal2 = this->m_attrV[2]->operator[](dd) ;
+
+		//Geom::Vector<3,double> n1 = normal1 ;
+		//Geom::Vector<3,double> n2 = normal2 ;
 
 		VEC3 newN = slerp(normal1,normal2,t) ; // spherical interpolation
 		newN.normalize() ;
+		assert (!isnan(newN[0])) ;
 
 		VEC3 newT = normal2 ^ normal1 ; // i is perpendicular to newNormal
 		newT.normalize() ;
@@ -104,10 +110,6 @@ bool Approximator_HemiFuncCoefs<PFP>::init()
 	assert((m_newFrameT.isValid() && m_newFrameB.isValid() && m_frameN.isValid())
 			|| !"New frame embeddings are not set") ;
 
-	// get quadric
-	m_quadricHF = this->m_map.template getAttribute<QuadricHF<REAL>, VERTEX>("HFquadric") ;
-	// Does not require to be valid (if it is not, altenatives will be used).
-
 	if(this->m_predictor)
 	{
 		return false ;
@@ -121,59 +123,78 @@ void Approximator_HemiFuncCoefs<PFP>::approximate(Dart d)
 {
 	const Dart& dd = this->m_map.phi1(d) ;
 
-	QuadricHF<REAL> q ;
-	if (m_quadricHF.isValid()) // if the selector does provide a quadricHF
+	// get coefs of v1 and v2
+	std::vector<VEC3> coefs1, coefs2 ;
+	coefs1.resize(m_nbCoefs) ; coefs2.resize(m_nbCoefs) ;
+	for (unsigned int i = 0 ; i < m_nbCoefs ; ++i)
 	{
-		q = m_quadricHF[d] ;
-		q += m_quadricHF[dd] ;
+		coefs1[i] = m_coefs[i]->operator[](d) ;
+		coefs2[i] = m_coefs[i]->operator[](dd) ;
 	}
-	else // the selector does not provide a quadricHF
+
+	const VEC3& T = m_newFrameT[d] ;
+	const VEC3& T1 = m_frameT[d] ;
+	const VEC3& T2 = m_frameT[dd] ;
+	const VEC3& B1 = m_frameB[d] ;
+	const VEC3& B2 = m_frameB[dd] ;
+	const VEC3& N = m_newFrameN[d] ;
+	const VEC3& N1 = m_frameN[d] ;
+	const VEC3& N2 = m_frameN[dd] ;
+
+	//assert(abs(T * N1) < 1e-6 || !"Approximator_FrameInterpolation<PFP>::approximate: T is not located in the first tangent plane") ;
+	//assert(abs(T * N2) < 1e-6 || !"Approximator_FrameInterpolation<PFP>::approximate: T is not located in the second tangent plane") ;
+
+	// Compute D1' and D2'
+	VEC3 B1prime = N1 ^ T ;
+	B1prime.normalize() ;
+	VEC3 B2prime = N2 ^ T ;
+	B2prime.normalize() ;
+
+	// Rotation dans sens trigo dans le plan tangent autour de N (T1 --> T)
+	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
+
+	assert (m_quadricHF.isValid() || !"LightfieldApproximator:approximate quadricHF is not valid") ;
+
+	if (T == T1)
 	{
-		// get coefs of v1 and v2
-		std::vector<VEC3> coefs1, coefs2 ;
-		coefs1.resize(m_nbCoefs) ; coefs2.resize(m_nbCoefs) ;
+		// std::cout << "Fallback to p1" << std::endl ;
+
+		m_quadricHF[d] = QuadricHF<REAL>(coefs1, gamma1, 0) ;
+		m_quadricHF[d] += QuadricHF<REAL>(coefs2, gamma2, alpha2 - alpha1) ;
+
 		for (unsigned int i = 0 ; i < m_nbCoefs ; ++i)
-		{
-			coefs1[i] = m_coefs[i]->operator[](d) ;
-			coefs2[i] = m_coefs[i]->operator[](dd) ;
-		}
-
-		const VEC3& T = m_newFrameT[d] ;
-		const VEC3& T1 = m_frameT[d] ;
-		const VEC3& T2 = m_frameT[dd] ;
-		const VEC3& B1 = m_newFrameB[d] ;
-		const VEC3& B2 = m_newFrameB[d] ;
-		const VEC3& N = m_frameN[d] ;
-		const VEC3& N1 = m_newFrameN[d] ;
-		const VEC3& N2 = m_newFrameN[dd] ;
-
-		assert(T * (N1 ^ N2) > 0.99 || !"T is not an intersection of the two tangent planes") ;
-
-		// Compute D1' and D2'
-		VEC3 B1prime = N1 ^ T ;
-		B1prime.normalize() ;
-		VEC3 B2prime = N2 ^ T ;
-		B2prime.normalize() ;
-
-		// Rotation dans sens trigo dans le plan tangent autour de N (T1 --> T)
-		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
+			this->m_approx[i][d] = m_coefs[i]->operator[](d) ;
+	}
+	else if (T == T2)
+	{
+		// std::cout << "Fallback to p2" << std::endl ;
 
-		// Create and sum quadrics
-		q = QuadricHF<REAL>(coefs1, gamma1, alpha1) ;
-		q += QuadricHF<REAL>(coefs2, gamma2, alpha2) ;
+		m_quadricHF[d] = QuadricHF<REAL>(coefs1, gamma1, alpha1 - alpha2) ;
+		m_quadricHF[d] += QuadricHF<REAL>(coefs2, gamma2, 0) ;
+
+		for (unsigned int i = 0 ; i < m_nbCoefs ; ++i)
+			this->m_approx[i][d] = m_coefs[i]->operator[](dd) ;
 	}
+	else
+	{
+		// Create and sum quadrics
+		m_quadricHF[d] = QuadricHF<REAL>(coefs1, gamma1, alpha1) ;
+		m_quadricHF[d] += QuadricHF<REAL>(coefs2, gamma2, alpha2) ;
 
-	std::vector<VEC3> coefs ;
-	// Compute new function
-	bool opt = q.findOptimizedCoefs(coefs) ;
-	// copy result
-	for (unsigned int i = 0 ; i < m_nbCoefs ; ++i)
-		this->m_approx[i][d] = opt ? (m_coefs[i]->operator[](d) + m_coefs[i]->operator[](dd))/2 : m_coefs[i]->operator[](d) ; // if fail take first one
-	assert(false || !"TODO: what todo when fail") ;
+		std::vector<VEC3> coefs ;
+		// Compute new function
+		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 take first one
+
+		if (!opt)
+			std::cerr << "LightfieldApproximator::Inversion failed : not treated" << std::endl ;
+	}
 }
 
 

include/Geometry/vector_gen.h

@@ -183,7 +183,7 @@ T tripleProduct(const Vector<DIM,T>& v1, const Vector<DIM,T>& v2, const Vector<D
 
 // returns a spherical interpolation of two vectors considering parameter t ((0 <= t <= 1) => result between v1 and v2)
 template <unsigned int DIM, typename T>
-Vector<DIM,T> slerp(const Vector<DIM,T> &v1, const Vector<DIM,T> &v2, const T &t) ;
+Vector<DIM,T> slerp(Vector<DIM,T> v1, Vector<DIM,T> v2, const T &t) ;
 
 /**********************************************/
 /*           SOME USEFUL TYPEDEFS             */

include/Geometry/vector_gen.hpp

@@ -336,8 +336,6 @@ std::istream& operator>>(std::istream& in, Vector<DIM,T>& v)
 	return in ;
 }
 
-
-
 template <unsigned int DIM, typename T>
 inline Vector<DIM,T> operator*(T a, const Vector<DIM,T>& v)
 {
@@ -357,11 +355,19 @@ inline T tripleProduct(const Vector<DIM,T>& v1, const Vector<DIM,T>& v2, const V
 }
 
 template <unsigned int DIM, typename T>
-inline Vector<DIM,T> slerp(const Vector<DIM,T>& v1, const Vector<DIM,T>& v2, const T& t)
+inline Vector<DIM,T> slerp(Vector<DIM,T> v1, Vector<DIM,T> v2, const T& t)
 {
 	Vector<DIM,T> res ;
 
-	T omega = acos(v1 * v2) ;	// assume v1,v2 normalized
+	T scal = v1 * v2 ;
+
+	// Prevention for floating point errors
+	if (1 < scal && scal < 1 + 1e-6)
+		scal = T(1) ;
+	if (-1. - 1e-6 < scal && scal < -1)
+		scal = -T(1) ;
+
+	T omega = acos(scal) ;
 	T den = sin(omega) ;
 
 	if (-1e-8 < den && den < 1e-8)