diff --git a/include/Algo/Import/importMesh.hpp b/include/Algo/Import/importMesh.hpp index 8dca5e6a31ef461cff31044f1cc505bc55bb1c59..f71c2a1b9c82bf0db5cf8372d3d36ffa20b7767a 100644 --- a/include/Algo/Import/importMesh.hpp +++ b/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; diff --git a/include/Algo/Import/importNodeEle.hpp b/include/Algo/Import/importNodeEle.hpp new file mode 100644 index 0000000000000000000000000000000000000000..8608287c761d52cf31e63f4647bab250815ef639 --- /dev/null +++ b/include/Algo/Import/importNodeEle.hpp @@ -0,0 +1,248 @@ +/******************************************************************************* +* 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 + +namespace CGoGN +{ + +namespace Algo +{ + +namespace Import +{ + +template +bool importNodeWithELERegions(typename PFP::MAP& map, const std::string& filenameNode, const std::string& filenameELE, std::vector& attrNames) +{ + typedef typename PFP::VEC3 VEC3; + + AttributeHandler position = map.template addAttribute(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 > > 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 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 > 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(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 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 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& vec = vecDartsPerVertex[map.phi1(d)]; + + Dart good_dart = NIL; + for(typename std::vector::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 diff --git a/include/Topology/map/map2MR/map2MR_PM.h b/include/Topology/map/map2MR/map2MR_PM.h index 7e70c3baffc08874b70bdc29374ac354cc2ff290..26c32d0318045413962f9865f220699eac490750 100644 --- a/include/Topology/map/map2MR/map2MR_PM.h +++ b/include/Topology/map/map2MR/map2MR_PM.h @@ -72,7 +72,7 @@ public: void analysis() ; void synthesis() ; } ; -} ; + } // namespace CGoGN diff --git a/include/Topology/map/map3MR/filters_Primal.h b/include/Topology/map/map3MR/filters_Primal.h index 5b2d014c1b0c695e1cf7a28c3c621cff782525bd..c1e2c1be3dd543cdbf5cf0c09ea60b10264f8c84 100644 --- a/include/Topology/map/map3MR/filters_Primal.h +++ b/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 -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 trav(m_map) ; - for (Dart d = trav.begin(); d != trav.end(); d = trav.next()) + float a = 0.5; + + TraversorW 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(m_map, d, m_position); + + unsigned int count = 0; + typename PFP::VEC3 ec(0); + Traversor3WE 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 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 travF(m_map) ; + for (Dart d = travF.begin(); d != travF.end(); d = travF.next()) + { + typename PFP::VEC3 vf = Algo::Geometry::faceCentroid(m_map, d, m_position); + + typename PFP::VEC3 ef(0); + unsigned int count = 0; + Traversor3FE 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 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 +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 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 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 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 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(m_map, m_position, db) ; - m_position[db] += p ; + + typename PFP::VEC3 fv(0); + unsigned int count = 0; + Traversor2VF 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 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 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 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 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 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 -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 trav(m_map) ; - for (Dart d = trav.begin(); d != trav.end(); d = trav.next()) + float a = 0.5; + + TraversorV 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; + } + } - unsigned int degree = m_map.vertexDegreeOnBoundary(db) ; - float n = 3.0/8.0 + 1.0/4.0 * cos(2.0 * M_PI / degree) ; - n = 8.0/5.0 * (n * n) ; + TraversorE 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); - m_position[db] *= n ; + m_map.incCurrentLevel() ; + Dart midE = m_map.phi2(db); + m_position[midE] *= a ; + m_map.decCurrentLevel() ; + } + else + { + m_map.incCurrentLevel() ; + Dart midE = m_map.phi2(d); + m_position[midE] *= a * a; + m_map.decCurrentLevel() ; + } + } + + TraversorF travF(m_map) ; + for (Dart d = travF.begin(); d != travF.end(); d = travF.next()) + { + if(!m_map.isBoundaryFace(d)) + { + m_map.incCurrentLevel() ; + Dart midF = m_map.phi2(m_map.phi1(d)); + m_position[midF] *= a ; + m_map.decCurrentLevel() ; } } } } ; +/* Loop on Boundary Vertices and SHW04 on Insides Vertices + *********************************************************************************/ template class LoopOddSynthesisFilter : public MRFilter { @@ -453,6 +743,62 @@ public: } } ; +template +class LoopNormalisationSynthesisFilter : 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) + {} + + void operator() () + { + TraversorV 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); + + unsigned int degree = m_map.vertexDegreeOnBoundary(db) ; + float n = 3.0/8.0 + 1.0/4.0 * cos(2.0 * M_PI / degree) ; + n = 8.0/5.0 * (n * n) ; + + m_position[db] *= n ; + } + } + } +} ; + +template +class LoopEvenSynthesisFilter : 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) + {} + + void operator() () + { + TraversorV 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 p = loopEvenVertex(m_map, m_position, db) ; + m_position[db] += p ; + } + } + } +} ; + template class LoopVolumeSynthesisFilter : public MRFilter {