selector TensorCurvature (non terminé)

include/Algo/Decimation/edgeSelector.h

@@ -31,6 +31,8 @@
 #include "Container/fakeAttribute.h"
 #include "Utils/qem.h"
 #include "Utils/quadricRGBfunctions.h"
+#include "Algo/Geometry/normal.h"
+#include "Algo/Selection/collector.h"
 #include "Algo/Geometry/curvature.h"
 #include "Algo/Geometry/area.h"
 
@@ -438,7 +440,7 @@ private:
 		bool valid ;
 		static std::string CGoGNnameOfType() { return "CurvatureTensorEdgeInfo" ; }
 	} CurvatureTensorEdgeInfo ;
-	typedef NoMathIOAttribute<CurvatureEdgeInfo> EdgeInfo ;
+	typedef NoMathIOAttribute<CurvatureTensorEdgeInfo> EdgeInfo ;
 
 	EdgeAttribute<EdgeInfo> edgeInfo ;
 	EdgeAttribute<REAL> edgeangle ;
@@ -449,7 +451,7 @@ private:
 	Approximator<PFP, VEC3,EDGE>* m_positionApproximator ;
 
 	void initEdgeInfo(Dart d) ;
-	void updateEdgeInfo(Dart d, bool recompute) ;
+	void updateEdgeInfo(Dart d) ; // TODO : usually has a 2nd arg (, bool recompute) : why ??
 	void computeEdgeInfo(Dart d, EdgeInfo& einfo) ;
 
 public:

include/Algo/Decimation/edgeSelector.hpp

@@ -1225,7 +1225,7 @@ void EdgeSelector_Curvature<PFP>::computeEdgeInfo(Dart d, EdgeInfo& einfo)
  ************************************************************************************/
 
 template <typename PFP>
-bool CurvatureTensor<PFP>::init()
+bool EdgeSelector_CurvatureTensor<PFP>::init()
 {
 	MAP& m = this->m_map ;
 
@@ -1247,7 +1247,7 @@ bool CurvatureTensor<PFP>::init()
 
 	edges.clear() ;
 
-//	TraversorE<MAP> travE(m);
+	TraversorE<MAP> travE(m);
 //	for(Dart dit = travE.begin() ; dit != travE.end() ; dit = travE.next())
 //	{
 //		computeEdgeMatrix(dit);
@@ -1264,7 +1264,7 @@ bool CurvatureTensor<PFP>::init()
 }
 
 template <typename PFP>
-bool CurvatureTensor<PFP>::nextEdge(Dart& d)
+bool EdgeSelector_CurvatureTensor<PFP>::nextEdge(Dart& d)
 {
 	if(cur == edges.end() || edges.empty())
 		return false ;
@@ -1273,7 +1273,7 @@ bool CurvatureTensor<PFP>::nextEdge(Dart& d)
 }
 
 template <typename PFP>
-void CurvatureTensor<PFP>::updateBeforeCollapse(Dart d)
+void EdgeSelector_CurvatureTensor<PFP>::updateBeforeCollapse(Dart d)
 {
 	MAP& m = this->m_map ;
 	assert(!m.isBoundaryEdge(d));
@@ -1318,25 +1318,31 @@ void CurvatureTensor<PFP>::updateBeforeCollapse(Dart d)
 
 
 template <typename PFP>
-void CurvatureTensor<PFP>::updateAfterCollapse(Dart d2, Dart dd2)
+void EdgeSelector_CurvatureTensor<PFP>::updateAfterCollapse(Dart d2, Dart dd2)
 {
 	MAP& m = this->m_map ;
+	CellMarkerStore<EDGE> eMark (m);
 
-	// TODO : here update edge angles : should look like :
-//	edgeangle[d] = computeAngleBetweenNormalsOnEdge<PFP>(map, d, this->m_position) ;
-
+	// update edge angles
+	Traversor2VF<MAP> tf (m,d2);
+	for(Dart dit = tf.begin(); dit != tf.end(); dit = tf.next())
+	{
+		Traversor2FE<MAP> te (m,dit);
+		for(Dart dit2 = te.begin(); dit2 != te.end(); dit2=te.next())
+		{
+			if (!eMark.isMarked(dit2))
+			{
+				edgeangle[dit2] = Algo::Geometry::computeAngleBetweenNormalsOnEdge<PFP>(m, dit2, this->m_position) ;
+				eMark.mark(dit2);
+			}
+		}
+	}
 
-	// update the edge matrices
-//	Traversor2VE<MAP> te (m,d2);
-//	for(Dart dit = te.begin() ; dit != te.end() ; dit = te.next())
-//	{
-//		computeEdgeMatrix(dit);
-//	}
 
 	// update the multimap
 
 	Traversor2VVaE<MAP> tv (m,d2);
-	CellMarker<EDGE> eMark (m);
+	eMark.unmarkAll();
 
 	for(Dart dit = tv.begin() ; dit != tv.end() ; dit = tv.next())
 	{
@@ -1355,7 +1361,7 @@ void CurvatureTensor<PFP>::updateAfterCollapse(Dart d2, Dart dd2)
 }
 
 template <typename PFP>
-void CurvatureTensor<PFP>::initEdgeInfo(Dart d)
+void EdgeSelector_CurvatureTensor<PFP>::initEdgeInfo(Dart d)
 {
 	MAP& m = this->m_map ;
 	EdgeInfo einfo ;
@@ -1367,7 +1373,7 @@ void CurvatureTensor<PFP>::initEdgeInfo(Dart d)
 }
 
 template <typename PFP>
-void CurvatureTensor<PFP>::updateEdgeInfo(Dart d)
+void EdgeSelector_CurvatureTensor<PFP>::updateEdgeInfo(Dart d)
 {
 	MAP& m = this->m_map ;
 	EdgeInfo& einfo = edgeInfo[d] ;
@@ -1384,46 +1390,44 @@ void CurvatureTensor<PFP>::updateEdgeInfo(Dart d)
 template <typename PFP>
 void EdgeSelector_CurvatureTensor<PFP>::computeEdgeInfo(Dart d, EdgeInfo& einfo)
 {
-
-	// TODO : code this function : everything complex is here
-
 	MAP& m = this->m_map ;
-	Dart dd = m.phi2(d);
-	Geom::Matrix33f M1; // init zero included
-	Geom::Matrix33f M2; // init zero included
+	Dart dd = m.phi2(d) ;
 
-	assert(! m.isBoundaryEdge(d));
+	unsigned int v1 = m.template getEmbedding<VERTEX>(d) ;
+	unsigned int v2 = m.template getEmbedding<VERTEX>(dd) ;
 
-	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();
-	}
+	m_positionApproximator->approximate(d) ;
 
-	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();
-	}
+	// compute tensor before collapse
+	typename PFP::MATRIX33 tens1;
+	Algo::Selection::Collector_OneRing_AroundEdge<PFP> col1 (m);
+	col1.collectAll(d);
+	col1.computeNormalCyclesTensor(this->m_position,edgeangle,tens1); // edgeangle is up to date
+	Algo::Geometry::normalCycles_SortTensor<PFP>(tens1);
 
-	m_positionApproximator->approximate(d) ;
-	const VEC3& a = m_positionApproximator->getApprox(d) ;
+	// temporary edge collapse
+	Dart d2 = m.phi2(m.phi_1(d)) ;
+	Dart dd2 = m.phi2(m.phi_1(dd)) ;
+	m.extractTrianglePair(d) ;
+	const unsigned int newV = m.template setOrbitEmbeddingOnNewCell<VERTEX>(d2) ;
+	this->m_position[newV] = m_positionApproximator->getApprox(d) ;
 
-	const VEC3 av1 = a - this->m_position[d] ;
-	const VEC3 av2 = a - this->m_position[dd] ;
+	// compute tensor after collapse
+	typename PFP::MATRIX33 tens2;
+	Algo::Selection::Collector_OneRing<PFP> col2 (m);
+	col2.collectAll(d);
+	col2.computeNormalCyclesTensor(this->m_position,tens2); // edgeangle is not up to date
+	Algo::Geometry::normalCycles_SortTensor<PFP>(tens2);
 
-	REAL err = av1 * (M1 * av1) + av2 * (M2 * av2);
+	// vertex split to reset the initial connectivity and embeddings
+	m.insertTrianglePair(d, d2, dd2) ;
+	m.template setOrbitEmbedding<VERTEX>(d, v1) ;
+	m.template setOrbitEmbedding<VERTEX>(dd, v2) ;
+
+	// TODO : compute err from the tensors
+	REAL err = 0 ;
 
+	// update the priority queue and edgeinfo
 	einfo.it = edges.insert(std::make_pair(err, d)) ;
 	einfo.valid = true ;
 }

include/Algo/Selection/collector.h

@@ -118,6 +118,7 @@ public:
 		return 0.0;
 	}
 	virtual void computeNormalCyclesTensor (const VertexAttribute<VEC3>& pos, const EdgeAttribute<REAL>& edgeangle, typename PFP::MATRIX33&) {assert(!"Warning: Collector<PFP>::computeNormalCyclesTensor() should be overloaded in non-virtual derived classes"); }
+	virtual void computeNormalCyclesTensor (const VertexAttribute<VEC3>& pos, typename PFP::MATRIX33&) {assert(!"Warning: Collector<PFP>::computeNormalCyclesTensor() should be overloaded in non-virtual derived classes"); }
 
 };
 
@@ -147,6 +148,7 @@ public:
 
 	REAL computeArea(const VertexAttribute<VEC3>& pos);
 	void computeNormalCyclesTensor (const VertexAttribute<VEC3>& pos, const EdgeAttribute<REAL>&edgeangle, typename PFP::MATRIX33&);
+	void computeNormalCyclesTensor (const VertexAttribute<VEC3>& pos, typename PFP::MATRIX33&);
 };
 
 /*********************************************************
@@ -177,6 +179,7 @@ public:
 
 	REAL computeArea(const VertexAttribute<VEC3>& pos);
 	void computeNormalCyclesTensor (const VertexAttribute<VEC3>& pos, const EdgeAttribute<REAL>&edgeangle, typename PFP::MATRIX33&);
+	void computeNormalCyclesTensor (const VertexAttribute<VEC3>& pos, typename PFP::MATRIX33&);
 };
 
 
@@ -215,6 +218,7 @@ public:
 
 	REAL computeArea(const VertexAttribute<VEC3>& pos);
 	void computeNormalCyclesTensor (const VertexAttribute<VEC3>& pos, const EdgeAttribute<REAL>& edgeangle, typename PFP::MATRIX33&);
+	void computeNormalCyclesTensor (const VertexAttribute<VEC3>& pos, typename PFP::MATRIX33&);
 };
 
 /*********************************************************

include/Algo/Selection/collector.hpp

@@ -168,6 +168,23 @@ void Collector_OneRing<PFP>::computeNormalCyclesTensor (const VertexAttribute<VE
 	tensor /= computeArea(pos) ;
 }
 
+template <typename PFP>
+void Collector_OneRing<PFP>::computeNormalCyclesTensor (const VertexAttribute<VEC3>& pos, typename PFP::MATRIX33& tensor){
+	assert(this->isInsideCollected || !"computeNormalCyclesTensor: inside cells have not been collected.") ;
+
+	tensor.zero() ;
+
+	// collect edges inside the neighborhood
+	for (std::vector<Dart>::const_iterator it = this->insideEdges.begin(); it != this->insideEdges.end(); ++it)
+	{
+		const VEC3 e = Algo::Geometry::vectorOutOfDart<PFP>(this->map, *it, pos) ;
+		const REAL edgeangle = Algo::Geometry::computeAngleBetweenNormalsOnEdge<PFP>(this->map, *it, pos) ;
+		tensor += Geom::transposed_vectors_mult(e,e) * edgeangle * (1 / e.norm()) ;
+	}
+
+	tensor /= computeArea(pos) ;
+}
+
 /*********************************************************
  * Collector One Ring aroung edge
  *********************************************************/
@@ -257,6 +274,24 @@ void Collector_OneRing_AroundEdge<PFP>::computeNormalCyclesTensor (const VertexA
 	tensor /= computeArea(pos) ;
 }
 
+template <typename PFP>
+void Collector_OneRing_AroundEdge<PFP>::computeNormalCyclesTensor (const VertexAttribute<VEC3>& pos, typename PFP::MATRIX33& tensor){
+	assert(this->isInsideCollected || !"computeNormalCyclesTensor: inside cells have not been collected.") ;
+
+	tensor.zero() ;
+
+	// collect edges inside the neighborhood
+	for (std::vector<Dart>::const_iterator it = this->insideEdges.begin(); it != this->insideEdges.end(); ++it)
+	{
+		const VEC3 e = Algo::Geometry::vectorOutOfDart<PFP>(this->map, *it, pos) ;
+		const REAL edgeangle = Algo::Geometry::computeAngleBetweenNormalsOnEdge<PFP>(this->map, *it, pos) ;
+		tensor += Geom::transposed_vectors_mult(e,e) * edgeangle * (1 / e.norm()) ;
+	}
+
+	tensor /= computeArea(pos) ;
+}
+
+
 /*********************************************************
  * Collector Within Sphere
  *********************************************************/
@@ -434,6 +469,33 @@ void Collector_WithinSphere<PFP>::computeNormalCyclesTensor (const VertexAttribu
 }
 
 
+template <typename PFP>
+void Collector_WithinSphere<PFP>::computeNormalCyclesTensor (const VertexAttribute<VEC3>& pos, typename PFP::MATRIX33& tensor){
+	assert(this->isInsideCollected || !"computeNormalCyclesTensor: inside cells have not been collected.") ;
+
+	VEC3 centerPosition = pos[this->centerDart];
+	tensor.zero() ;
+
+	// collect edges inside the neighborhood
+	for (std::vector<Dart>::const_iterator it = this->insideEdges.begin(); it != this->insideEdges.end(); ++it)
+	{
+		const VEC3 e = Algo::Geometry::vectorOutOfDart<PFP>(this->map, *it, pos) ;
+		const REAL edgeangle = Algo::Geometry::computeAngleBetweenNormalsOnEdge<PFP>(this->map, *it, pos) ;
+		tensor += Geom::transposed_vectors_mult(e,e) * edgeangle * (1 / e.norm()) ;
+	}
+	// collect edges crossing the neighborhood's border
+	for (std::vector<Dart>::const_iterator it = this->border.begin(); it != this->border.end(); ++it)
+	{
+		const VEC3 e = Algo::Geometry::vectorOutOfDart<PFP>(this->map, *it, pos) ;
+		REAL alpha ;
+		Algo::Geometry::intersectionSphereEdge<PFP>(this->map, centerPosition, radius, *it, pos, alpha) ;
+		const REAL edgeangle = Algo::Geometry::computeAngleBetweenNormalsOnEdge<PFP>(this->map, *it, pos) ;
+		tensor += Geom::transposed_vectors_mult(e,e) * edgeangle * (1 / e.norm()) * alpha ;
+	}
+
+	tensor /= computeArea(pos) ;
+}
+
 /*********************************************************
  * Collector Normal Angle (Vertices)
  *********************************************************/