approximator NormalArea

include/Algo/Decimation/approximator.h

@@ -43,6 +43,7 @@ enum ApproximatorType
 	A_CornerCutting,
 	A_TangentPredict1,
 	A_TangentPredict2,
+	A_NormalArea,
 	A_ColorNaive,
 	A_ColorQEMext,
 	A_Lightfield,

include/Algo/Decimation/decimation.hpp

@@ -65,6 +65,9 @@ void decimate(
 		case A_TangentPredict2 :
 			approximators.push_back(new Approximator_MidEdge<PFP>(map, attribs)) ;
 			break ;
+		case A_NormalArea :
+			approximators.push_back(new Approximator_NormalArea<PFP>(map, attribs)) ;
+			break ;
 		case A_hHalfCollapse :
 			approximators.push_back(new Approximator_HalfCollapse<PFP>(map, attribs)) ;
 			break ;

include/Algo/Decimation/edgeSelector.hpp

@@ -781,7 +781,7 @@ bool EdgeSelector_NormalArea<PFP>::init()
 	{
 		if((*it)->getApproximatedAttributeName() == "position")
 		{
-			assert((*it)->getType() == A_MidEdge || !"Only MidEdge Approximator is valid") ;
+//			assert((*it)->getType() == A_MidEdge || (*it)->getType() == A_NormalArea || !"Only MidEdge and NormalArea Approximator are valid") ;
 			m_positionApproximator = reinterpret_cast<Approximator<PFP, VEC3,EDGE>* >(*it) ;
 			ok = true ;
 		}
@@ -991,7 +991,7 @@ void EdgeSelector_NormalArea<PFP>::computeEdgeMatrix(Dart d)
 	edgeMatrix[d] *= e.norm2();
 	edgeMatrix[d] -= Geom::transposed_vectors_mult(e,e) ;
 	// TODO : test : try to normalize by area
-//	edgeMatrix[d] /= e.norm2(); // pas d'amélioration significative par rapport à la version sans normalisation
+	edgeMatrix[d] /= e.norm2(); // pas d'amélioration significative par rapport à la version sans normalisation
 //	edgeMatrix[d] /= e.norm2() * e.norm2(); // étonnament bon : sur certains maillages équivalant à la QEMml
 }
 

include/Algo/Decimation/geometryApproximator.h

@@ -26,6 +26,7 @@
 #define __GEOMETRY_APPROXIMATOR_H__
 
 #include "Algo/Decimation/approximator.h"
+#include <Eigen/Dense>
 
 namespace CGoGN
 {
@@ -144,6 +145,32 @@ public:
 	void approximate(Dart d) ;
 } ;
 
+template <typename PFP>
+class Approximator_NormalArea : public Approximator<PFP, typename PFP::VEC3, EDGE>
+{
+public:
+	typedef typename PFP::MAP MAP ;
+	typedef typename PFP::VEC3 VEC3 ;
+	typedef typename PFP::REAL REAL ;
+
+protected:
+//	VertexAttribute<Utils::Quadric<REAL> > m_quadric ;
+	EdgeAttribute<Geom::Matrix<3,3,REAL> > edgeMatrix ;
+
+public:
+	Approximator_NormalArea(MAP& m, std::vector<VertexAttribute<VEC3>* > pos, Predictor<PFP, VEC3>* pred = NULL) :
+		Approximator<PFP, VEC3, EDGE>(m, pos, pred)
+	{
+		assert(pos.size() > 0 || !"Approximator_NormalArea: attribute vector is empty") ;
+	}
+	~Approximator_NormalArea()
+	{}
+	ApproximatorType getType() const { return A_NormalArea ; }
+	bool init() ;
+	void approximate(Dart d) ;
+} ;
+
+
 } //namespace Decimation
 
 } //namespace Algo

include/Algo/Decimation/geometryApproximator.hpp

@@ -368,6 +368,74 @@ void Approximator_CornerCutting<PFP>::approximate(Dart d)
 	}
 }
 
+/************************************************************************************
+ *                            NORMAL AREA METRIC                                    *
+ ************************************************************************************/
+template <typename PFP>
+bool Approximator_NormalArea<PFP>::init()
+{
+	edgeMatrix = this->m_map.template getAttribute<Geom::Matrix<3,3,REAL>, EDGE>("NormalAreaMatrix") ;
+	assert(edgeMatrix.isValid());
+
+//	m_quadric = this->m_map.template getAttribute<Utils::Quadric<REAL>, VERTEX>("QEMquadric") ;
+	// Does not require to be valid (if it is not, altenatives will be used).
+
+	if(this->m_predictor)
+	{
+		return false ;
+	}
+	return true ;
+}
+
+template <typename PFP>
+void Approximator_NormalArea<PFP>::approximate(Dart d)
+{
+	MAP& m = this->m_map ;
+	Dart dd = m.phi2(d);
+	Geom::Matrix33f M1; // init zero included
+	Geom::Matrix33f M2; // init zero included
+
+	assert(! m.isBoundaryEdge(d));
+
+	Traversor2VF<MAP> td (m,d);
+	Dart it = td.begin();
+	it = td.next();
+	Dart it2 = td.next();
+	while( it2 != td.end())
+	{
+		M1 += edgeMatrix[m.phi1(it)];
+		it = it2;
+		it2 = td.next();
+	}
+
+	Traversor2VF<MAP> tdd (m,dd);
+	it = tdd.begin();
+	it = tdd.next();
+	it2 = tdd.next();
+	while( it2 != tdd.end())
+	{
+		M2 += edgeMatrix[m.phi1(it)];
+		it = it2;
+		it2 = tdd.next();
+	}
+
+	const VEC3 & v1 = (*this->m_attrV[0])[d] ;
+	const VEC3 & v2 = (*this->m_attrV[0])[dd] ;
+
+
+	Eigen::Matrix3f A ;
+	A << M1(0,0)+M2(0,0) , M1(0,1)+M2(0,1) , M1(0,2)+M2(0,2) , M1(1,0)+M2(1,0) , M1(1,1)+M2(1,1) , M1(1,2)+M2(1,2) , M1(2,0)+M2(2,0) , M1(2,1)+M2(2,1) , M1(2,2)+M2(2,2) ;
+
+	VEC3 mb = M1*v1 + M2*v2 ;
+	Eigen::Vector3f b (mb[0],mb[1],mb[2]);
+
+	Eigen::LDLT<Eigen::Matrix3f> decompo (A);
+	Eigen::Vector3f x = decompo.solve(b);
+
+	this->m_approx[0][d] = VEC3 (x(0),x(1),x(2)) ;
+}
+
+
 } //namespace Decimation
 
 } //namespace Algo