lots of modification in MR

CMakeLists.txt

@@ -128,7 +128,7 @@ IF(WIN32)
 	set(CMAKE_CONFIGURATION_TYPES "${CMAKE_CONFIGURATION_TYPES}" CACHE STRING "Only Release or Debug" FORCE)
 #	set(CMAKE_CONFIGURATION_TYPES "Release Debug" CACHE STRING "Only Release or Debug" FORCE)
 ELSE(WIN32)
-	SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall")
+	SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -fPIC")
 	add_subdirectory(Release)
 	add_subdirectory(Debug)
 	add_subdirectory(${CGoGN_ROOT_DIR}/Apps Apps)

include/Algo/Multiresolution/map2MR/map2MR_PM.h

@@ -76,6 +76,8 @@ public:
 
 	void addNewLevel(unsigned int percentWantedVertices);
 
+	void collapseEdge(Dart d);
+
 	//coarsen the mesh -> analysis
 	void coarsen() ;
 

include/Algo/Multiresolution/map2MR/map2MR_PM.hpp

@@ -35,8 +35,7 @@ namespace Multiresolution
 template <typename PFP>
 Map2MR_PM<PFP>::Map2MR_PM(typename PFP::MAP& map, VertexAttribute<typename PFP::VEC3>& position) : m_map(map), m_position(position)
 {
-	//TODO
-	//map.enableMR();
+
 }
 
 template <typename PFP>
@@ -48,15 +47,11 @@ Map2MR_PM<PFP>::~Map2MR_PM()
 		delete (*it) ;
 	for(typename std::vector<Algo::Decimation::PredictorGen<PFP>*>::iterator it = m_predictors.begin(); it != m_predictors.end(); ++it)
 		delete (*it) ;
-
-	//TODO
-	//map.disableMR()
 }
 
 template <typename PFP>
-void Map2MR_PM<PFP>::createPM(Algo::Decimation::SelectorType s, Algo::Decimation::ApproximatorType a, const FunctorSelect& select = allDarts)
+void Map2MR_PM<PFP>::createPM(Algo::Decimation::SelectorType s, Algo::Decimation::ApproximatorType a, const FunctorSelect& select)
 {
-
 	CGoGNout << "  creating approximator and predictor.." << CGoGNflush ;
 	switch(a)
 	{
@@ -143,8 +138,9 @@ void Map2MR_PM<PFP>::addNewLevel(unsigned int percentWantedVertices)
 {
 	// level handling
 	m_map.pushLevel() ;
-	m_map.addLevel();
-	m_map.setCurrentLevel(m_map.getMaxLevel()) ;
+	m_map.addLevelBack();
+	m_map.duplicateDarts(m_map.getMaxLevel());
+	m_map.setCurrentLevel(m_map.getMaxLevel());
 
 	unsigned int nbVertices = m_map.template getNbOrbits<VERTEX>() ;
 	unsigned int nbWantedVertices = nbVertices * percentWantedVertices / 100 ;
@@ -161,8 +157,6 @@ void Map2MR_PM<PFP>::addNewLevel(unsigned int percentWantedVertices)
 		Dart d2 = m_map.phi2(m_map.phi_1(d)) ;
 		Dart dd2 = m_map.phi2(m_map.phi_1(m_map.phi2(d))) ;
 
-
-
 		for(typename std::vector<Algo::Decimation::ApproximatorGen<PFP>*>::iterator it = m_approximators.begin(); it != m_approximators.end(); ++it)
 		{
 			(*it)->approximate(d) ;					// compute approximated attributes with its associated detail
@@ -171,11 +165,7 @@ void Map2MR_PM<PFP>::addNewLevel(unsigned int percentWantedVertices)
 
 		m_selector->updateBeforeCollapse(d) ;		// update selector
 
-		//m_map.collapseEdge(d);
-
-		inactiveMarker.markOrbit<FACE>(d) ;
-		inactiveMarker.markOrbit<FACE>(m_map.phi2(d)) ;
-		m_map.extractTrianglePair(d);
+		collapseEdge(d);
 
 		unsigned int newV = m_map.template embedNewCell<VERTEX>(d2) ;
 		unsigned int newE1 = m_map.template embedNewCell<EDGE>(d2) ;
@@ -198,6 +188,29 @@ void Map2MR_PM<PFP>::addNewLevel(unsigned int percentWantedVertices)
 	CGoGNout << "..done (" << nbVertices << " vertices)" << CGoGNendl ;
 }
 
+
+template <typename PFP>
+void Map2MR_PM<PFP>::collapseEdge(Dart d)
+{
+	//duplication :
+	m_map.duplicateMRDart(m_map.phi2(m_map.phi1(d)));
+	m_map.duplicateMRDart(m_map.phi2(m_map.phi_1(d)));
+	m_map.duplicateMRDart(m_map.phi2(m_map.phi1(m_map.phi2(d))));
+	m_map.duplicateMRDart(m_map.phi2(m_map.phi_1(m_map.phi2(d))));
+
+	//effacer :
+//	m_map.unrefMRDart(m_map.phi1(m_map.phi2(d)));
+//	m_map.unrefMRDart(m_map.phi_1(m_map.phi2(d)));
+//	m_map.unrefMRDart(m_map.phi2(d));
+//
+//	m_map.unrefMRDart(m_map.phi1(d));
+//	m_map.unrefMRDart(m_map.phi_1(d));
+//	m_map.unrefMRDart(d);
+
+	m_map.collapseEdge(d);
+}
+
+
 template <typename PFP>
 void Map2MR_PM<PFP>::coarsen()
 {

include/Topology/map/map2MR/map2MR_PrimalRegular.h → include/Algo/Multiresolution/map2MR/map2MR_PrimalRegular.h

@@ -34,23 +34,42 @@
 namespace CGoGN
 {
 
-class Map2MR_PrimalRegular : public EmbeddedMap2
+namespace Algo
 {
+
+namespace MR
+{
+
+namespace Primal
+{
+
+namespace Regular
+{
+
+template <typename PFP>
+class Map2MR
+{
+
+public:
+	typedef typename PFP::MAP MAP ;
+	typedef typename PFP::VEC3 VEC3 ;
+	typedef typename PFP::REAL REAL ;
+
 protected:
+	MAP& m_map;
 	bool shareVertexEmbeddings ;
 
-	std::vector<Multiresolution::MRFilter*> synthesisFilters ;
-	std::vector<Multiresolution::MRFilter*> analysisFilters ;
+	std::vector<CGoGN::Multiresolution::MRFilter*> synthesisFilters ;
+	std::vector<CGoGN::Multiresolution::MRFilter*> analysisFilters ;
 
 public:
-	Map2MR_PrimalRegular() ;
+	Map2MR(MAP& map) ;
 
-	virtual std::string mapTypeName() const { return "Map2MR_PrimalRegular" ; }
 
 	void addNewLevel(bool embedNewVertices = true) ;
 
-	void addSynthesisFilter(Multiresolution::MRFilter* f) { synthesisFilters.push_back(f) ; }
-	void addAnalysisFilter(Multiresolution::MRFilter* f) { analysisFilters.push_back(f) ; }
+	void addSynthesisFilter(CGoGN::Multiresolution::MRFilter* f) { synthesisFilters.push_back(f) ; }
+	void addAnalysisFilter(CGoGN::Multiresolution::MRFilter* f) { analysisFilters.push_back(f) ; }
 
 	void clearSynthesisFilters() { synthesisFilters.clear() ; }
 	void clearAnalysisFilters() { analysisFilters.clear() ; }
@@ -59,6 +78,16 @@ public:
 	void synthesis() ;
 } ;
 
+} // namespace Regular
+
+} // namespace Primal
+
+} // namespace MR
+
+} // namespace Algo
+
 } // namespace CGoGN
 
+#include "Algo/Multiresolution/map2MR/map2MR_PrimalRegular.hpp"
+
 #endif

src/Topology/map/map2MR/map2MR_PrimalRegular.cpp → include/Algo/Multiresolution/map2MR/map2MR_PrimalRegular.hpp

@@ -22,23 +22,37 @@
 *                                                                              *
 *******************************************************************************/
 
-#include "Topology/map/map2MR/map2MR_PrimalRegular.h"
-
 namespace CGoGN
 {
 
-Map2MR_PrimalRegular::Map2MR_PrimalRegular() :
+namespace Algo
+{
+
+namespace MR
+{
+
+namespace Primal
+{
+
+namespace Regular
+{
+
+template <typename PFP>
+Map2MR<PFP>::Map2MR(typename PFP::MAP& map) :
+	m_map(map),
 	shareVertexEmbeddings(true)
 {
-	initMR() ;
+
 }
 
-void Map2MR_PrimalRegular::addNewLevel(bool embedNewVertices)
+template <typename PFP>
+void Map2MR<PFP>::addNewLevel(bool embedNewVertices)
 {
-	pushLevel() ;
+	m_map.pushLevel() ;
 
-	addLevel() ;
-	setCurrentLevel(getMaxLevel()) ;
+	m_map.addLevelBack() ;
+	m_map.duplicateDarts(m_map.getMaxLevel());
+	m_map.setCurrentLevel(m_map.getMaxLevel()) ;
 
 //	for(unsigned int i = m_mrattribs.begin(); i != m_mrattribs.end(); m_mrattribs.next(i))
 //	{
@@ -49,102 +63,113 @@ void Map2MR_PrimalRegular::addNewLevel(bool embedNewVertices)
 //	}
 
 	// cut edges
-	TraversorE<Map2MR_PrimalRegular> travE(*this) ;
+	TraversorE<typename PFP::MAP> travE(m_map) ;
 	for (Dart d = travE.begin(); d != travE.end(); d = travE.next())
 	{
 		if(!shareVertexEmbeddings && embedNewVertices)
 		{
-			if(getEmbedding<VERTEX>(d) == EMBNULL)
-				embedNewCell<VERTEX>(d) ;
-			if(getEmbedding<VERTEX>(phi1(d)) == EMBNULL)
-				embedNewCell<VERTEX>(d) ;
+			if(m_map.template getEmbedding<VERTEX>(d) == EMBNULL)
+				m_map.template embedNewCell<VERTEX>(d) ;
+			if(m_map.template getEmbedding<VERTEX>(m_map.phi1(d)) == EMBNULL)
+				m_map.template embedNewCell<VERTEX>(d) ;
 		}
 
-		cutEdge(d) ;
+		m_map.cutEdge(d) ;
 		travE.skip(d) ;
-		travE.skip(phi1(d)) ;
+		travE.skip(m_map.phi1(d)) ;
 
 		if(embedNewVertices)
-			embedNewCell<VERTEX>(phi1(d)) ;
+			m_map.template embedNewCell<VERTEX>(m_map.phi1(d)) ;
 	}
 
 	// split faces
-	TraversorF<Map2MR_PrimalRegular> travF(*this) ;
+	TraversorF<typename PFP::MAP> travF(m_map) ;
 	for (Dart d = travF.begin(); d != travF.end(); d = travF.next())
 	{
 		Dart old = d ;
-		if(getDartLevel(old) == getMaxLevel())
-			old = phi1(old) ;
+		if(m_map.getDartLevel(old) == m_map.getMaxLevel())
+			old = m_map.phi1(old) ;
 
-		decCurrentLevel() ;
-		unsigned int degree = faceDegree(old) ;
-		incCurrentLevel() ;
+		m_map.decCurrentLevel() ;
+		unsigned int degree = m_map.faceDegree(old) ;
+		m_map.incCurrentLevel() ;
 
 		if(degree == 3)					// if subdividing a triangle
 		{
-			Dart dd = phi1(old) ;
-			Dart e = phi1(phi1(dd)) ;
-			splitFace(dd, e) ;
+			Dart dd = m_map.phi1(old) ;
+			Dart e = m_map.phi1(m_map.phi1(dd)) ;
+			m_map.splitFace(dd, e) ;
 			travF.skip(dd) ;
 
 			dd = e ;
-			e = phi1(phi1(dd)) ;
-			splitFace(dd, e) ;
+			e = m_map.phi1(m_map.phi1(dd)) ;
+			m_map.splitFace(dd, e) ;
 			travF.skip(dd) ;
 
 			dd = e ;
-			e = phi1(phi1(dd)) ;
-			splitFace(dd, e) ;
+			e = m_map.phi1(m_map.phi1(dd)) ;
+			m_map.splitFace(dd, e) ;
 			travF.skip(dd) ;
 
 			travF.skip(e) ;
 		}
 		else							// if subdividing a polygonal face
 		{
-			Dart dd = phi1(old) ;
-			Dart next = phi1(phi1(dd)) ;
-			splitFace(dd, next) ;		// insert a first edge
+			Dart dd = m_map.phi1(old) ;
+			Dart next = m_map.phi1(m_map.phi1(dd)) ;
+			m_map.splitFace(dd, next) ;		// insert a first edge
 
-			Dart ne = alpha1(dd) ;
-			cutEdge(ne) ;				// cut the new edge to insert the central vertex
+			Dart ne = m_map.alpha1(dd) ;
+			m_map.cutEdge(ne) ;				// cut the new edge to insert the central vertex
 			travF.skip(dd) ;
 
 			if(embedNewVertices)
-				embedNewCell<VERTEX>(phi1(ne)) ;
+				m_map.template embedNewCell<VERTEX>(m_map.phi1(ne)) ;
 
-			dd = phi1(phi1(next)) ;
+			dd = m_map.phi1(m_map.phi1(next)) ;
 			while(dd != ne)				// turn around the face and insert new edges
 			{							// linked to the central vertex
-				Dart tmp = phi1(ne) ;
-				splitFace(tmp, dd) ;
+				Dart tmp = m_map.phi1(ne) ;
+				m_map.splitFace(tmp, dd) ;
 				travF.skip(tmp) ;
-				dd = phi1(phi1(dd)) ;
+				dd = m_map.phi1(m_map.phi1(dd)) ;
 			}
 			travF.skip(ne) ;
 		}
 	}
 
-	popLevel() ;
+	m_map.popLevel() ;
 }
 
-void Map2MR_PrimalRegular::analysis()
+template <typename PFP>
+void Map2MR<PFP>::analysis()
 {
-	assert(getCurrentLevel() > 0 || !"analysis : called on level 0") ;
+	assert(m_map.getCurrentLevel() > 0 || !"analysis : called on level 0") ;
 
-	decCurrentLevel() ;
+	m_map.decCurrentLevel() ;
 
 	for(unsigned int i = 0; i < analysisFilters.size(); ++i)
 		(*analysisFilters[i])() ;
 }
 
-void Map2MR_PrimalRegular::synthesis()
+template <typename PFP>
+void Map2MR<PFP>::synthesis()
 {
-	assert(getCurrentLevel() < getMaxLevel() || !"synthesis : called on max level") ;
+	assert(m_map.getCurrentLevel() < m_map.getMaxLevel() || !"synthesis : called on max level") ;
 
 	for(unsigned int i = 0; i < synthesisFilters.size(); ++i)
 		(*synthesisFilters[i])() ;
 
-	incCurrentLevel() ;
+	m_map.incCurrentLevel() ;
 }
 
+
+} // namespace Regular
+
+} // namespace Primal
+
+} // namespace MR
+
+} // namespace Algo
+
 } // namespace CGoGN

include/Algo/Multiresolution/map3MR/filters_Primal.h

@@ -872,6 +872,298 @@ public:
 	}
 } ;
 
+/* Catmull-clark on Boundary Vertices and MJ96 on Insides Vertices
+ *********************************************************************************/
+template <typename PFP>
+class MJ96VertexSubdivision : public MRFilter
+{
+protected:
+	typename PFP::MAP& m_map ;
+	VertexAttribute<typename PFP::VEC3>& m_position ;
+
+public:
+	MJ96VertexSubdivision(typename PFP::MAP& m, VertexAttribute<typename PFP::VEC3>& p) : m_map(m), m_position(p)
+	{}
+
+	void operator() ()
+	{
+		TraversorV<typename PFP::MAP> trav(m_map) ;
+		for (Dart d = trav.begin(); d != trav.end(); d = trav.next())
+		{
+			if(m_map.isBoundaryVertex(d))
+			{
+				Dart db = m_map.findBoundaryFaceOfVertex(d);
+
+				typename PFP::VEC3 np1(0) ;
+				typename PFP::VEC3 np2(0) ;
+				unsigned int degree1 = 0 ;
+				unsigned int degree2 = 0 ;
+				Dart it = db ;
+				do
+				{
+					++degree1 ;
+					Dart dd = m_map.phi1(it) ;
+					np1 += m_position[dd] ;
+					Dart end = m_map.phi_1(it) ;
+					dd = m_map.phi1(dd) ;
+					do
+					{
+						++degree2 ;
+						np2 += m_position[dd] ;
+						dd = m_map.phi1(dd) ;
+					} while(dd != end) ;
+					it = m_map.phi2(m_map.phi_1(it)) ;
+				} while(it != db) ;
+
+				float beta = 3.0 / (2.0 * degree1) ;
+				float gamma = 1.0 / (4.0 * degree2) ;
+				np1 *= beta / degree1 ;
+				np2 *= gamma / degree2 ;
+
+				typename PFP::VEC3 vp = m_position[db] ;
+				vp *= 1.0 - beta - gamma ;
+
+				m_map.incCurrentLevel() ;
+
+				m_position[d] = np1 + np2 + vp ;
+
+				m_map.decCurrentLevel() ;
+
+			}
+			else
+			{
+				typename PFP::VEC3 P = m_position[d];
+
+				//vertex points
+				typename PFP::VEC3 Cavg = typename PFP::VEC3(0);
+				unsigned int degree = 0;
+				Traversor3VW<typename PFP::MAP> travVW(m_map, d);
+				for(Dart dit = travVW.begin() ; dit != travVW.end() ; dit = travVW.next())
+				{
+					Cavg += Algo::Geometry::volumeCentroid<PFP>(m_map, dit, m_position);
+					++degree;
+				}
+				Cavg /= degree;
+
+				typename PFP::VEC3 Aavg = typename PFP::VEC3(0);
+				degree = 0;
+				Traversor3VF<typename PFP::MAP> travVF(m_map, d);
+				for(Dart dit = travVF.begin() ; dit != travVF.end() ; dit = travVF.next())
+				{
+					Aavg += Algo::Geometry::faceCentroid<PFP>(m_map, dit, m_position);
+					++degree;
+				}
+				Aavg /= degree;
+
+				typename PFP::VEC3 Mavg = typename PFP::VEC3(0);
+				degree = 0;
+				Traversor3VE<typename PFP::MAP> travVE(m_map, d);
+				for(Dart dit = travVE.begin() ; dit != travVE.end() ; dit = travVE.next())
+				{
+					Dart d2 = m_map.phi2(dit);
+					Aavg += (m_position[dit] + m_position[d2]) * typename PFP::REAL(0.5);
+					++degree;
+				}
+				Aavg /= degree;
+
+				typename PFP::VEC3 vp = Cavg + Aavg * 3 + Mavg * 3 + P;
+				vp /= 8;
+
+				m_map.incCurrentLevel() ;
+
+				m_position[d] = P;//vp;
+
+				m_map.decCurrentLevel() ;
+			}
+		}
+	}
+};
+
+template <typename PFP>
+class MJ96EdgeSubdivision : public MRFilter
+{
+protected:
+	typename PFP::MAP& m_map ;
+	VertexAttribute<typename PFP::VEC3>& m_position ;
+
+public:
+	MJ96EdgeSubdivision(typename PFP::MAP& m, VertexAttribute<typename PFP::VEC3>& p) : m_map(m), m_position(p)
+	{}
+
+	void operator() ()
+	{
+		TraversorE<typename PFP::MAP> trav(m_map) ;
+		for (Dart d = trav.begin(); d != trav.end(); d = trav.next())
+		{
+			if(m_map.isBoundaryEdge(d))
+			{
+				Dart db = m_map.findBoundaryFaceOfEdge(d);
+
+				Dart d1 = m_map.phi2(db) ;
+				Dart d2 = m_map.phi2(d1) ;
+				Dart d3 = m_map.phi_1(d1) ;
+				Dart d4 = m_map.phi_1(d2) ;
+				Dart d5 = m_map.phi1(m_map.phi1(d1)) ;
+				Dart d6 = m_map.phi1(m_map.phi1(d2)) ;
+
+				typename PFP::VEC3 p1 = m_position[d1] ;
+				typename PFP::VEC3 p2 = m_position[d2] ;
+				typename PFP::VEC3 p3 = m_position[d3] ;
+				typename PFP::VEC3 p4 = m_position[d4] ;
+				typename PFP::VEC3 p5 = m_position[d5] ;
+				typename PFP::VEC3 p6 = m_position[d6] ;
+
+				p1 *= 3.0 / 8.0 ;
+				p2 *= 3.0 / 8.0 ;
+				p3 *= 1.0 / 16.0 ;
+				p4 *= 1.0 / 16.0 ;
+				p5 *= 1.0 / 16.0 ;
+				p6 *= 1.0 / 16.0 ;
+
+				m_map.incCurrentLevel() ;
+
+				Dart midV = m_map.phi2(d);
+
+				m_position[midV] = p1 + p2 + p3 + p4 + p5 + p6 ;
+
+				m_map.decCurrentLevel() ;
+			}
+			else
+			{
+				//edge points
+				typename PFP::VEC3 Cavg = typename PFP::VEC3(0);
+				unsigned int degree = 0;
+				Traversor3EW<typename PFP::MAP> travEW(m_map, d);
+				for(Dart dit = travEW.begin() ; dit != travEW.end() ; dit = travEW.next())
+				{
+					Cavg += Algo::Geometry::volumeCentroid<PFP>(m_map, dit, m_position);
+					++degree;
+				}
+				Cavg /= degree;
+
+				typename PFP::VEC3 Aavg = typename PFP::VEC3(0);
+				degree = 0;
+				Traversor3EF<typename PFP::MAP> travEF(m_map, d);
+				for(Dart dit = travEF.begin() ; dit != travEF.end() ; dit = travEF.next())
+				{
+					Aavg += Algo::Geometry::faceCentroid<PFP>(m_map, dit, m_position);
+					++degree;
+				}
+				Aavg /= degree;
+
+				Dart d2 = m_map.phi2(d);
+				typename PFP::VEC3 M = (m_position[d] + m_position[d2]) * typename PFP::REAL(0.5);
+
+				typename PFP::VEC3 ep = Cavg + Aavg * 2 + M * (degree - 3);
+				ep /= degree;
+
+				m_map.incCurrentLevel() ;
+
+				Dart midV = m_map.phi2(d);
+
+				m_position[midV] = ep;
+
+				m_map.decCurrentLevel() ;
+			}
+		}
+	}
+};
+
+template <typename PFP>
+class MJ96FaceSubdivision : public MRFilter
+{
+protected:
+	typename PFP::MAP& m_map ;
+	VertexAttribute<typename PFP::VEC3>& m_position ;
+
+public:
+	MJ96FaceSubdivision(typename PFP::MAP& m, VertexAttribute<typename PFP::VEC3>& p) : m_map(m), m_position(p)
+	{}
+
+	void operator() ()
+	{
+		TraversorF<typename PFP::MAP> trav(m_map) ;
+		for (Dart d = trav.begin(); d != trav.end(); d = trav.next())
+		{
+			if(m_map.isBoundaryFace(d))
+			{
+				Dart db = m_map.phi3(d);
+
+				typename PFP::VEC3 p(0) ;
+				unsigned int degree = 0 ;
+				Traversor2FV<typename PFP::MAP> trav(m_map, db) ;
+				for(Dart it = trav.begin(); it != trav.end(); it = trav.next())
+				{
+					++degree ;
+					p += m_position[it] ;
+				}
+				p /= degree ;
+
+				m_map.incCurrentLevel() ;
+
+				Dart df = m_map.phi1(m_map.phi1(d)) ;
+
+				m_position[df] = p ;
+
+				m_map.decCurrentLevel() ;
+			}
+			else
+			{
+				//face points
+				typename PFP::VEC3 C0 = Algo::Geometry::volumeCentroid<PFP>(m_map, d, m_position);
+				typename PFP::VEC3 C1 = Algo::Geometry::volumeCentroid<PFP>(m_map, m_map.phi3(d), m_position);
+
+				typename PFP::VEC3 A = Algo::Geometry::faceCentroid<PFP>(m_map, m_map.phi3(d), m_position);
+
+				typename PFP::VEC3 fp = C0 + A * 2 + C1;
+				fp /= 4;
+
+				m_map.incCurrentLevel() ;
+
+				Dart df = m_map.phi1(m_map.phi1(d)) ;
+				m_position[df] = fp;
+
+				m_map.decCurrentLevel() ;
+			}
+		}
+	}
+};
+
+template <typename PFP>
+class MJ96VolumeSubdivision : public MRFilter
+{
+protected:
+	typename PFP::MAP& m_map ;
+	VertexAttribute<typename PFP::VEC3>& m_position ;
+
+public:
+	MJ96VolumeSubdivision(typename PFP::MAP& m, VertexAttribute<typename PFP::VEC3>& p) : m_map(m), m_position(p)
+	{}
+
+	void operator() ()
+	{
+		TraversorW<typename PFP::MAP> trav(m_map) ;
+		for (Dart d = trav.begin(); d != trav.end(); d = trav.next())
+		{
+			//cell points : these points are the average of the
+			//vertices of the lattice that bound the cell
+			typename PFP::VEC3 p = Algo::Geometry::volumeCentroid<PFP>(m_map, d, m_position);
+
+			m_map.incCurrentLevel() ;
+
+			if(!Algo::Modelisation::Tetrahedralization::isTetrahedron<PFP>(m_map,d))
+			{
+				Dart midV = m_map.phi_1(m_map.phi2(m_map.phi1(d)));
+				m_position[midV] = p ;
+			}
+
+			m_map.decCurrentLevel() ;
+
+		}
+	}
+};
+
+
 
 } // namespace Multiresolution