Merge cgogn:~untereiner/CGoGN

include/Algo/Import/importMesh.hpp

@@ -271,6 +271,9 @@ bool importMeshV(typename PFP::MAP& map, const std::string& filename, std::vecto
 	if ((filename.rfind(".off") != std::string::npos) || (filename.rfind(".OFF") != std::string::npos))
 		kind = ImportVolumique::OFF;
 
+	if ((filename.rfind(".node") != std::string::npos) || (filename.rfind(".NODE") != std::string::npos))
+		kind = ImportVolumique::NODE;
+
 	if ((filename.rfind(".ts") != std::string::npos) || (filename.rfind(".TS") != std::string::npos))
 		kind = ImportVolumique::TS;
 

include/Algo/Import/importNodeEle.hpp

+/*******************************************************************************
+* CGoGN: Combinatorial and Geometric modeling with Generic N-dimensional Maps  *
+* version 0.1                                                                  *
+* Copyright (C) 2009-2012, IGG Team, LSIIT, University of Strasbourg           *
+*                                                                              *
+* This library is free software; you can redistribute it and/or modify it      *
+* under the terms of the GNU Lesser General Public License as published by the *
+* Free Software Foundation; either version 2.1 of the License, or (at your     *
+* option) any later version.                                                   *
+*                                                                              *
+* This library is distributed in the hope that it will be useful, but WITHOUT  *
+* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or        *
+* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License  *
+* for more details.                                                            *
+*                                                                              *
+* You should have received a copy of the GNU Lesser General Public License     *
+* along with this library; if not, write to the Free Software Foundation,      *
+* Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301 USA.           *
+*                                                                              *
+* Web site: http://cgogn.unistra.fr/                                           *
+* Contact information: cgogn@unistra.fr                                        *
+*                                                                              *
+*******************************************************************************/
+
+#include "Algo/Modelisation/polyhedron.h"
+#include <vector>
+
+namespace CGoGN
+{
+
+namespace Algo
+{
+
+namespace Import 
+{
+
+template <typename PFP>
+bool importNodeWithELERegions(typename PFP::MAP& map, const std::string& filenameNode, const std::string& filenameELE, std::vector<std::string>& attrNames)
+{
+	typedef typename PFP::VEC3 VEC3;
+
+	AttributeHandler<VEC3> position = map.template addAttribute<VEC3>(VERTEX, "position") ;
+	attrNames.push_back(position.name()) ;
+
+	AttributeContainer& container = map.getAttributeContainer(VERTEX) ;
+
+	unsigned int m_nbVertices = 0, m_nbFaces = 0, m_nbEdges = 0, m_nbVolumes = 0;
+	AutoAttributeHandler< NoMathIONameAttribute< std::vector<Dart> > > vecDartsPerVertex(map, VERTEX, "incidents");
+
+	//open file
+	std::ifstream fnode(filenameNode.c_str(), std::ios::in);
+	if (!fnode.good())
+	{
+		CGoGNerr << "Unable to open file " << filenameNode << CGoGNendl;
+		return false;
+	}
+
+	std::ifstream fele(filenameELE.c_str(), std::ios::in);
+	if (!fele.good())
+	{
+		CGoGNerr << "Unable to open file " << filenameELE << CGoGNendl;
+		return false;
+	}
+
+	std::string line;
+
+//	do
+//	{
+//		std::getline(fnode,line);
+//	}while(line.rfind("#") == 0);
+
+	//Reading NODE file
+	//First line: [# of points] [dimension (must be 3)] [# of attributes] [# of boundary markers (0 or 1)]
+	unsigned int nbe;
+	{
+		do
+		{
+			std::getline(fnode,line);
+		}while(line.size() == 0);
+
+		std::stringstream oss(line);
+		oss >> m_nbVertices;
+		oss >> nbe;
+		oss >> nbe;
+		oss >> nbe;
+	}
+
+	//Reading number of tetrahedra in ELE file
+	unsigned int nbv;
+	{
+		do
+		{
+			std::getline(fele,line);
+		}while(line.size() == 0);
+
+		std::stringstream oss(line);
+		oss >> m_nbVolumes;
+		oss >> nbv ; oss >> nbv;
+	}
+
+	CGoGNout << "nb points = " << m_nbVertices << " / nb faces = " << m_nbFaces << " / nb edges = " << m_nbEdges << " / nb tet = " << m_nbVolumes << CGoGNendl;
+
+	//Reading vertices
+	//Remaining lines: [point #] [x] [y] [z] [optional attributes] [optional boundary marker]
+	std::vector<unsigned int> verticesID;
+	verticesID.reserve(m_nbVertices);
+
+	for(unsigned int i = 0 ; i < m_nbVertices ; ++i)
+	{
+		do
+		{
+			std::getline(fnode,line);
+		}while(line.size() == 0);
+
+		std::stringstream oss(line);
+
+		int idv;
+		oss >> idv;
+
+		float x,y,z;
+		oss >> x;
+		oss >> y;
+		oss >> z;
+		//we can read colors informations if exists
+		VEC3 pos(x,y,z);
+
+		unsigned int id = container.insertLine();
+		position[id] = pos;
+
+		verticesID.push_back(id);
+	}
+
+	std::vector<std::vector<Dart> > vecDartPtrEmb;
+	vecDartPtrEmb.reserve(m_nbVertices);
+
+	DartMarkerNoUnmark m(map) ;
+
+	//Read and embed tetrahedra TODO
+	for(unsigned i = 0; i < m_nbVolumes ; ++i)
+	{
+		do
+		{
+			std::getline(fele,line);
+		} while(line.size() == 0);
+
+		std::stringstream oss(line);
+		oss >> nbe;
+
+		Dart d = Algo::Modelisation::createTetrahedron<PFP>(map);
+		Geom::Vec4ui pt;
+		oss >> pt[0];
+		--(pt[0]);
+		oss >> pt[1];
+		--(pt[1]);
+		oss >> pt[2];
+		--(pt[2]);
+		oss >> pt[3];
+		--(pt[3]);
+
+		//regions ?
+		//oss >> nbe;
+
+		// Embed three vertices
+		for(unsigned int j = 0 ; j < 3 ; ++j)
+		{
+			FunctorSetEmb<typename PFP::MAP> fsetemb(map, VERTEX, verticesID[pt[2-j]]);
+			map.foreach_dart_of_orbit(PFP::MAP::ORBIT_IN_PARENT(VERTEX), d, fsetemb);
+
+			//store darts per vertices to optimize reconstruction
+			Dart dd = d;
+			do
+			{
+				m.mark(dd) ;
+				vecDartsPerVertex[pt[2-j]].push_back(dd);
+				dd = map.phi1(map.phi2(dd));
+			} while(dd != d);
+
+			d = map.phi1(d);
+
+		}
+
+		//Embed the last vertex
+		d = map.phi_1(map.phi2(d));
+
+		FunctorSetEmb<typename PFP::MAP> fsetemb(map, VERTEX, verticesID[pt[3]]);
+		map.foreach_dart_of_orbit( PFP::MAP::ORBIT_IN_PARENT(VERTEX), d, fsetemb);
+
+		//store darts per vertices to optimize reconstruction
+		Dart dd = d;
+		do
+		{
+			m.mark(dd) ;
+			vecDartsPerVertex[pt[3]].push_back(dd);
+			dd = map.phi1(map.phi2(dd));
+		} while(dd != d);
+
+	}
+
+	fnode.close();
+	fele.close();
+
+	//Association des phi3
+	unsigned int nbBoundaryFaces = 0 ;
+	for (Dart d = map.begin(); d != map.end(); map.next(d))
+	{
+		if (m.isMarked(d))
+		{
+			std::vector<Dart>& vec = vecDartsPerVertex[map.phi1(d)];
+
+			Dart good_dart = NIL;
+			for(typename std::vector<Dart>::iterator it = vec.begin(); it != vec.end() && good_dart == NIL; ++it)
+			{
+				if(map.getEmbedding(VERTEX, map.phi1(*it)) == map.getEmbedding(VERTEX, d) &&
+				   map.getEmbedding(VERTEX, map.phi_1(*it)) == map.getEmbedding(VERTEX, map.phi_1(d)) /*&&
+				   map.getEmbedding(VERTEX, *it) == map.getEmbedding(VERTEX, map.phi1(d)) */)
+				{
+					good_dart = *it ;
+				}
+			}
+
+			if (good_dart != NIL)
+			{
+				map.sewVolumes(d, good_dart, false);
+				m.unmarkOrbit(FACE, d);
+			}
+			else
+			{
+				m.unmarkOrbit(PFP::MAP::ORBIT_IN_PARENT(FACE), d);
+				++nbBoundaryFaces;
+			}
+		}
+	}
+
+	if (nbBoundaryFaces > 0)
+	{
+		std::cout << "closing" << std::endl ;
+		map.closeMap();
+		CGoGNout << "Map closed (" << nbBoundaryFaces << " boundary faces)" << CGoGNendl;
+	}
+
+	return true;
+}
+
+} // namespace Import
+
+} // namespace Algo
+
+} // namespace CGoGN

include/Topology/map/map2MR/map2MR_PM.h

@@ -72,7 +72,7 @@ public:
 	void analysis() ;
 	void synthesis() ;
 } ;
-} ;
+
 
 } // namespace CGoGN
 

include/Topology/map/map3MR/filters_Primal.h

@@ -352,64 +352,354 @@ public:
 	}
 } ;
 
-/* Loop on Boundary Vertices and SHW04 on Insides Vertices
+/* Ber02
  *********************************************************************************/
+
+//w-lift(a)
 template <typename PFP>
-class LoopEvenSynthesisFilter : public MRFilter
+class Ber02OddSynthesisFilter : public MRFilter
 {
 protected:
 	typename PFP::MAP& m_map ;
 	typename PFP::TVEC3& m_position ;
 
 public:
-	LoopEvenSynthesisFilter(typename PFP::MAP& m, typename PFP::TVEC3& p) : m_map(m), m_position(p)
+	Ber02OddSynthesisFilter(typename PFP::MAP& m, typename PFP::TVEC3& 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())
+		float a = 0.5;
+
+		TraversorW<typename PFP::MAP> travW(m_map) ;
+		for (Dart d = travW.begin(); d != travW.end(); d = travW.next())
 		{
-			if(m_map.isBoundaryVertex(d))
+			typename PFP::VEC3 vc = Algo::Geometry::volumeCentroid<PFP>(m_map, d, m_position);
+
+			unsigned int count = 0;
+			typename PFP::VEC3 ec(0);
+			Traversor3WE<typename PFP::MAP> travWE(m_map, d);
+			for (Dart dit = travWE.begin(); dit != travWE.end(); dit = travWE.next())
+			{
+				m_map.incCurrentLevel();
+				ec += m_position[m_map.phi2(dit)];
+				m_map.decCurrentLevel();
+				++count;
+			}
+			ec /= count;
+
+			count = 0;
+			typename PFP::VEC3 fc(0);
+			Traversor3WF<typename PFP::MAP> travWF(m_map, d);
+			for (Dart dit = travWF.begin(); dit != travWF.end(); dit = travWF.next())
+			{
+				m_map.incCurrentLevel();
+				fc += m_position[m_map.phi2(m_map.phi1(dit))];
+				m_map.decCurrentLevel();
+				++count;
+			}
+
+			fc /= count;
+
+			m_map.incCurrentLevel() ;
+			Dart midV = m_map.phi_1(m_map.phi2(m_map.phi1(d)));
+			m_position[midV] += 8 * a * a * a * vc + 12 * a * a * ec + 6 * a * fc;
+			m_map.decCurrentLevel() ;
+		}
+
+		TraversorF<typename PFP::MAP> travF(m_map) ;
+		for (Dart d = travF.begin(); d != travF.end(); d = travF.next())
+		{
+			typename PFP::VEC3 vf = Algo::Geometry::faceCentroid<PFP>(m_map, d, m_position);
+
+			typename PFP::VEC3 ef(0);
+			unsigned int count = 0;
+			Traversor3FE<typename PFP::MAP> travFE(m_map, d);
+			for (Dart dit = travFE.begin(); dit != travFE.end(); dit = travFE.next())
+			{
+				m_map.incCurrentLevel();
+				ef += m_position[m_map.phi2(dit)];
+				m_map.decCurrentLevel();
+				++count;
+			}
+			ef /= count;
+
+			m_map.incCurrentLevel() ;
+			Dart midF = m_map.phi2(m_map.phi1(d));
+			m_position[midF] += vf * 4.0 * a * a + ef * 4.0 * a;
+			m_map.decCurrentLevel() ;
+		}
+
+		TraversorE<typename PFP::MAP> travE(m_map) ;
+		for (Dart d = travE.begin(); d != travE.end(); d = travE.next())
+		{
+			typename PFP::VEC3 ve = (m_position[d] + m_position[m_map.phi2(d)]) * typename PFP::REAL(0.5);
+
+			m_map.incCurrentLevel() ;
+			Dart midV = m_map.phi2(d) ;
+			m_position[midV] += ve * a * 2.0;
+			m_map.decCurrentLevel() ;
+		}
+
+	}
+} ;
+
+// s-lift(a)
+template <typename PFP>
+class Ber02EvenSynthesisFilter : public MRFilter
+{
+protected:
+	typename PFP::MAP& m_map ;
+	typename PFP::TVEC3& m_position ;
+
+public:
+	Ber02EvenSynthesisFilter(typename PFP::MAP& m, typename PFP::TVEC3& p) : m_map(m), m_position(p)
+	{}
+
+	void operator() ()
+	{
+		float a = 0.5;
+
+		TraversorV<typename PFP::MAP> travV(m_map);
+		for(Dart d = travV.begin() ; d != travV.end() ; d = travV.next())
+		{
+			if(!m_map.isBoundaryVertex(d))
+			{
+				typename PFP::VEC3 cv(0);
+				unsigned int count = 0;
+				Traversor3VW<typename PFP::MAP> travVW(m_map,d);
+				for(Dart dit = travVW.begin(); dit != travVW.end() ; dit = travVW.next())
+				{
+					m_map.incCurrentLevel() ;
+					Dart midV = m_map.phi_1(m_map.phi2(m_map.phi1(dit)));
+					cv += m_position[midV];
+					m_map.decCurrentLevel() ;
+					++count;
+				}
+				cv /= count;
+
+				typename PFP::VEC3 fv(0);
+				count = 0;
+				Traversor3VF<typename PFP::MAP> travVF(m_map,d);
+				for(Dart dit = travVF.begin(); dit != travVF.end() ; dit = travVF.next())
+				{
+					m_map.incCurrentLevel() ;
+					Dart midV = m_map.phi2(m_map.phi1(dit));
+					fv += m_position[midV];
+					m_map.decCurrentLevel() ;
+					++count;
+				}
+				fv /= count;
+
+				typename PFP::VEC3 ev(0);
+				count = 0;
+				Traversor3VE<typename PFP::MAP> travVE(m_map,d);
+				for(Dart dit = travVE.begin(); dit != travVE.end() ; dit = travVE.next())
+				{
+					m_map.incCurrentLevel() ;
+					Dart midV = m_map.phi2(dit);
+					ev += m_position[midV];
+					m_map.decCurrentLevel() ;
+					++count;
+				}
+				ev /= count;
+
+				m_position[d] += cv * 8 * a * a * a + fv * 12 * a * a + ev * 6 * a;
+			}
+			else
 			{
 				Dart db = m_map.findBoundaryFaceOfVertex(d);
-				typename PFP::VEC3 p = loopEvenVertex<PFP>(m_map, m_position, db) ;
-				m_position[db] += p ;
+
+				typename PFP::VEC3 fv(0);
+				unsigned int count = 0;
+				Traversor2VF<typename PFP::MAP> travVF(m_map,db);
+				for(Dart dit = travVF.begin(); dit != travVF.end() ; dit = travVF.next())
+				{
+					m_map.incCurrentLevel() ;
+					Dart midV = m_map.phi2(m_map.phi1(dit));
+					fv += m_position[midV];
+					m_map.decCurrentLevel() ;
+					++count;
+				}
+				fv /= count;
+
+				typename PFP::VEC3 ev(0);
+				count = 0;
+				Traversor2VE<typename PFP::MAP> travVE(m_map,db);
+				for(Dart dit = travVE.begin(); dit != travVE.end() ; dit = travVE.next())
+				{
+					m_map.incCurrentLevel() ;
+					Dart midV = m_map.phi2(dit);
+					ev += m_position[midV];
+					m_map.decCurrentLevel() ;
+					++count;
+				}
+				ev /= count;
+
+				m_position[db] += fv * 4 * a * a + ev * 4 * a;
+			}
+		}
+
+		TraversorE<typename PFP::MAP> travE(m_map);
+		for(Dart d = travE.begin() ; d != travE.end() ; d = travE.next())
+		{
+			if(m_map.isBoundaryEdge(d))
+			{
+				Dart db = m_map.findBoundaryFaceOfEdge(d);
+
+				typename PFP::VEC3 fe(0);
+
+				m_map.incCurrentLevel() ;
+				Dart midV = m_map.phi2(m_map.phi1(db));
+				fe += m_position[midV];
+				m_map.decCurrentLevel() ;
+
+				m_map.incCurrentLevel() ;
+				midV = m_map.phi2(m_map.phi1(m_map.phi2(db)));
+				fe += m_position[midV];
+				m_map.decCurrentLevel() ;
+
+				fe /= 2;
+
+				m_map.incCurrentLevel() ;
+				Dart midF = m_map.phi2(db);
+				m_position[midF] += fe * 2 * a;
+				m_map.decCurrentLevel() ;
+			}
+			else
+			{
+				typename PFP::VEC3 ce(0);
+				unsigned int count = 0;
+				Traversor3EW<typename PFP::MAP> travEW(m_map, d);
+				for(Dart dit = travEW.begin() ; dit != travEW.end() ; dit = travEW.next())
+				{
+					m_map.incCurrentLevel() ;
+					Dart midV = m_map.phi_1(m_map.phi2(m_map.phi1(dit)));
+					ce += m_position[midV];
+					m_map.decCurrentLevel() ;
+					++count;
+				}
+
+				ce /= count;
+
+				typename PFP::VEC3 fe(0);
+				count = 0;
+				Traversor3FW<typename PFP::MAP> travFW(m_map, d);
+				for(Dart dit = travFW.begin() ; dit != travFW.end() ; dit = travFW.next())
+				{
+					m_map.incCurrentLevel() ;
+					Dart midV = m_map.phi2(m_map.phi1(dit));
+					fe += m_position[midV];
+					m_map.decCurrentLevel() ;
+					++count;
+				}
+
+				fe /= count;
+
+				m_map.incCurrentLevel() ;
+				Dart midF = m_map.phi2(d);
+				m_position[midF] += ce * 4 * a * a + fe * 4 * a;
+				m_map.decCurrentLevel() ;
 			}
 		}
+
+		TraversorF<typename PFP::MAP> travF(m_map) ;
+		for (Dart d = travF.begin(); d != travF.end(); d = travF.next())
+		{
+			typename PFP::VEC3 cf(0);
+
+			m_map.incCurrentLevel();
+			Dart midV = m_map.phi_1(m_map.phi2(m_map.phi1(d)));
+			cf += m_position[midV];
+			m_map.decCurrentLevel();
+
+			if(!m_map.isBoundaryFace(d))
+			{
+				Dart d3 = m_map.phi3(d);
+				m_map.incCurrentLevel();
+				Dart midV = m_map.phi_1(m_map.phi2(m_map.phi1(d3)));
+				cf += m_position[midV];
+				m_map.decCurrentLevel();
+
+				cf /= 2;
+			}
+
+			m_map.incCurrentLevel() ;
+			Dart midF = m_map.phi2(m_map.phi1(d));
+			m_position[midF] += cf * 2 * a;
+			m_map.decCurrentLevel() ;
+		}
+
 	}
 } ;
 
+// s-scale(a)
 template <typename PFP>
-class LoopNormalisationSynthesisFilter : public MRFilter
+class Ber02ScaleSynthesisFilter : public MRFilter
 {
 protected:
 	typename PFP::MAP& m_map ;
 	typename PFP::TVEC3& m_position ;
 
 public:
-	LoopNormalisationSynthesisFilter(typename PFP::MAP& m, typename PFP::TVEC3& p) : m_map(m), m_position(p)
+	Ber02ScaleSynthesisFilter(typename PFP::MAP& m, typename PFP::TVEC3& 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())
+		float a = 0.5;
+
+		TraversorV<typename PFP::MAP> travV(m_map) ;
+		for (Dart d = travV.begin(); d != travV.end(); d = travV.next())
 		{
 			if(m_map.isBoundaryVertex(d))
 			{
 				Dart db = m_map.findBoundaryFaceOfVertex(d);
+				m_position[db] *= a * a;
+			}
+			else
+			{
+				m_position[d] *= a * a * a;
+			}
+		}