diff --git a/include/Algo/Modelisation/tetrahedralization.hpp b/include/Algo/Modelisation/tetrahedralization.hpp index 3422a380548a271daf1129961ce4890f9b98420a..28043312d125026f92832698445a2c31da08f456 100644 --- a/include/Algo/Modelisation/tetrahedralization.hpp +++ b/include/Algo/Modelisation/tetrahedralization.hpp @@ -443,8 +443,10 @@ void swapGen3To2(typename PFP::MAP& map, Dart d) template void swapGen2To3(typename PFP::MAP& map, Dart d) { + unsigned int n = map.edgeDegree(d); + //- a single 2-3 swap, followed by n − 3 3-2 swaps, or -//– a single 4-4 swap, followed by n − 4 3-2 swaps. +//- a single 4-4 swap, followed by n − 4 3-2 swaps. } diff --git a/include/Algo/Multiresolution/Map2MR/Filters/bertram.h b/include/Algo/Multiresolution/Map2MR/Filters/bertram.h new file mode 100644 index 0000000000000000000000000000000000000000..569208856193e23cad3f30611521648400fa3796 --- /dev/null +++ b/include/Algo/Multiresolution/Map2MR/Filters/bertram.h @@ -0,0 +1,659 @@ +/******************************************************************************* +* 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 * +* * +*******************************************************************************/ + +#ifndef __3MR_BERTRAM_FILTER__ +#define __3MR_BERTRAM_FILTER__ + +#include +#include "Algo/Geometry/centroid.h" +#include "Algo/Modelisation/tetrahedralization.h" +#include "Algo/Multiresolution/filter.h" + +namespace CGoGN +{ + +namespace Algo +{ + +namespace MR +{ + +namespace Primal +{ + +namespace Filters +{ + +/******************************************************************************* + * Without features preserving + *******************************************************************************/ + +// +// Synthesis +// + +//w-lift(a) +template +class Ber02OddSynthesisFilter : public Filter +{ +protected: + typename PFP::MAP& m_map ; + VertexAttribute& m_position ; + typename PFP::VEC3::DATA_TYPE m_a; + + +public: + Ber02OddSynthesisFilter(typename PFP::MAP& m, VertexAttribute& p, typename PFP::VEC3::DATA_TYPE a) : m_map(m), m_position(p), m_a(a) + {} + + void operator() () + { + TraversorF travF(m_map) ; + for (Dart d = travF.begin(); d != travF.end(); d = travF.next()) + { + typename PFP::VEC3 vf(0.0); + typename PFP::VEC3 ef(0.0); + + unsigned int count = 0; + Traversor2FE travFE(m_map, d); + for (Dart dit = travFE.begin(); dit != travFE.end(); dit = travFE.next()) + { + vf += m_position[dit]; + m_map.incCurrentLevel(); + ef += m_position[m_map.phi1(dit)]; + m_map.decCurrentLevel(); + ++count; + } + ef /= count; + ef *= 4.0 * m_a; + + vf /= count; + vf *= 4.0 * m_a * m_a; + + m_map.incCurrentLevel() ; + Dart midF = m_map.phi1(m_map.phi1(d)); + m_position[midF] += vf + ef ; + 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.phi1(d)]) * typename PFP::REAL(0.5); + ve *= 2.0 * m_a; + + m_map.incCurrentLevel() ; + Dart midV = m_map.phi1(d) ; + m_position[midV] += ve; + m_map.decCurrentLevel() ; + } + } +} ; + +// s-lift(a) +template +class Ber02EvenSynthesisFilter : public Filter +{ +protected: + typename PFP::MAP& m_map ; + VertexAttribute& m_position ; + typename PFP::VEC3::DATA_TYPE m_a; + bool m_threshold; + unsigned int m_current; + unsigned int percentWanted; + +public: + Ber02EvenSynthesisFilter(typename PFP::MAP& m, VertexAttribute& p, typename PFP::VEC3::DATA_TYPE a) : m_map(m), m_position(p), m_a(a), m_threshold(false) + {} + +// void activateThreshold(bool b) { m_threshold = b; percentWanted = m_map.getMaxLevel() * 30 / 100 ; } + + void operator() () + { + + TraversorV travV(m_map); + for(Dart d = travV.begin() ; d != travV.end() ; d = travV.next()) + { +// if(m_threshold && m_current > percentWanted) +// { +// m_position[d] += typename PFP::VEC3(0.0); +// } +// else +// { + typename PFP::VEC3 ev(0.0); + typename PFP::VEC3 fv(0.0); + if(m_map.isBoundaryVertex(d)) + { + Dart db = m_map.findBoundaryEdgeOfVertex(d); + m_map.incCurrentLevel() ; + ev = (m_position[m_map.phi1(db)] + m_position[m_map.phi_1(db)]); + m_map.decCurrentLevel() ; + ev *= m_a; + + m_position[d] += ev; + } + else + { + unsigned int count = 0; + + Traversor2VF travVF(m_map,d); + for(Dart dit = travVF.begin(); dit != travVF.end() ; dit = travVF.next()) + { + m_map.incCurrentLevel() ; + + Dart midEdgeV = m_map.phi1(dit); + ev += m_position[midEdgeV]; + fv += m_position[m_map.phi1(midEdgeV)]; + + m_map.decCurrentLevel() ; + ++count; + } + fv /= count; + fv *= 4 * m_a * m_a; + + ev /= count; + ev *= 4 * m_a; + + m_position[d] += fv + ev; + } +// } + } + + TraversorE travE(m_map); + for(Dart d = travE.begin() ; d != travE.end() ; d = travE.next()) + { +// if(m_threshold && m_current > percentWanted) +// { +// m_map.incCurrentLevel() ; +// Dart midF = m_map.phi1(d); +// m_position[midF] += typename PFP::VEC3(0.0); +// m_map.decCurrentLevel() ; +// } +// else +// { + unsigned int count = 0; + + typename PFP::VEC3 fe(0.0); + Traversor2EF travEF(m_map, d); + for(Dart dit = travEF.begin() ; dit != travEF.end() ; dit = travEF.next()) + { + m_map.incCurrentLevel() ; + Dart midV = m_map.phi1(m_map.phi1(dit)); + fe += m_position[midV]; + m_map.decCurrentLevel() ; + ++count; + } + + fe /= count; + fe *= 2 * m_a; + + m_map.incCurrentLevel() ; + Dart midF = m_map.phi1(d); + m_position[midF] += fe; + m_map.decCurrentLevel() ; +// } + } + +// if(m_threshold && m_current < percentWanted) +// { +// ++m_current; +// } + } +} ; + +// s-scale(a) +template +class Ber02ScaleSynthesisFilter : public Filter +{ +protected: + typename PFP::MAP& m_map ; + VertexAttribute& m_position ; + typename PFP::VEC3::DATA_TYPE m_a; + +public: + Ber02ScaleSynthesisFilter(typename PFP::MAP& m, VertexAttribute& p, typename PFP::VEC3::DATA_TYPE a) : m_map(m), m_position(p), m_a(a) + {} + + void operator() () + { + TraversorV travV(m_map) ; + for (Dart d = travV.begin(); d != travV.end(); d = travV.next()) + { + if(m_map.isBoundaryVertex(d)) + m_position[d] *= m_a; + else + m_position[d] *= m_a * m_a; + } + + TraversorE travE(m_map) ; + for (Dart d = travE.begin(); d != travE.end(); d = travE.next()) + { + m_map.incCurrentLevel() ; + Dart midE = m_map.phi1(d); + if(!m_map.isBoundaryVertex(midE)) + m_position[midE] *= m_a ; + m_map.decCurrentLevel() ; + } + } +} ; + +// +// Analysis +// + +//w-lift(a) +template +class Ber02OddAnalysisFilter : public Filter +{ +protected: + typename PFP::MAP& m_map ; + VertexAttribute& m_position ; + typename PFP::VEC3::DATA_TYPE m_a; + +public: + Ber02OddAnalysisFilter(typename PFP::MAP& m, VertexAttribute& p, typename PFP::VEC3::DATA_TYPE a) : m_map(m), m_position(p), m_a(a) + {} + + void operator() () + { + 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.phi1(d)]) * typename PFP::REAL(0.5); + ve *= 2.0 * m_a; + + m_map.incCurrentLevel() ; + Dart midV = m_map.phi1(d) ; + m_position[midV] -= ve; + m_map.decCurrentLevel() ; + } + + + TraversorF travF(m_map) ; + for (Dart d = travF.begin(); d != travF.end(); d = travF.next()) + { + typename PFP::VEC3 vf(0.0); + typename PFP::VEC3 ef(0.0); + + unsigned int count = 0; + Traversor2FE travFE(m_map, d); + for (Dart dit = travFE.begin(); dit != travFE.end(); dit = travFE.next()) + { + vf += m_position[dit]; + m_map.incCurrentLevel(); + ef += m_position[m_map.phi1(dit)]; + m_map.decCurrentLevel(); + ++count; + } + ef /= count; + ef *= 4.0 * m_a; + + vf /= count; + vf *= 4.0 * m_a * m_a; + + m_map.incCurrentLevel() ; + Dart midF = m_map.phi1(m_map.phi1(d)); + m_position[midF] -= vf + ef ; + m_map.decCurrentLevel() ; + } + + } +}; + +// s-lift(a) +template +class Ber02EvenAnalysisFilter : public Filter +{ +protected: + typename PFP::MAP& m_map ; + VertexAttribute& m_position ; + typename PFP::VEC3::DATA_TYPE m_a; + +public: + Ber02EvenAnalysisFilter(typename PFP::MAP& m, VertexAttribute& p, typename PFP::VEC3::DATA_TYPE a) : m_map(m), m_position(p), m_a(a) + {} + + void operator() () + { + TraversorE travE(m_map); + for(Dart d = travE.begin() ; d != travE.end() ; d = travE.next()) + { + if(!m_map.isBoundaryEdge(d)) + { + unsigned int count = 0; + + typename PFP::VEC3 fe(0); + Traversor2EF travEF(m_map, d); + for(Dart dit = travEF.begin() ; dit != travEF.end() ; dit = travEF.next()) + { + m_map.incCurrentLevel() ; + Dart midV = m_map.phi1(m_map.phi1(dit)); + fe += m_position[midV]; + m_map.decCurrentLevel() ; + ++count; + } + + fe /= count; + fe *= 2 * m_a; + + m_map.incCurrentLevel() ; + Dart midF = m_map.phi1(d); + m_position[midF] -= fe; + m_map.decCurrentLevel() ; + } + } + + TraversorV travV(m_map); + for(Dart d = travV.begin() ; d != travV.end() ; d = travV.next()) + { + typename PFP::VEC3 ev(0.0); + typename PFP::VEC3 fv(0.0); + if(m_map.isBoundaryVertex(d)) + { + Dart db = m_map.findBoundaryEdgeOfVertex(d); + m_map.incCurrentLevel() ; + ev = (m_position[m_map.phi1(db)] + m_position[m_map.phi_1(db)]); + m_map.decCurrentLevel() ; + ev *= m_a; + + m_position[d] -= ev; + } + else + { + unsigned int count = 0; + + Traversor2VF travVF(m_map,d); + for(Dart dit = travVF.begin(); dit != travVF.end() ; dit = travVF.next()) + { + m_map.incCurrentLevel() ; + + Dart midEdgeV = m_map.phi1(dit); + ev += m_position[midEdgeV]; + fv += m_position[m_map.phi1(midEdgeV)]; + + m_map.decCurrentLevel() ; + ++count; + } + fv /= count; + fv *= 4 * m_a * m_a; + + ev /= count; + ev *= 4 * m_a; + + m_position[d] -= fv + ev; + } + } + } +}; + +// s-scale(a) +template +class Ber02ScaleAnalysisFilter : public Filter +{ +protected: + typename PFP::MAP& m_map ; + VertexAttribute& m_position ; + typename PFP::VEC3::DATA_TYPE m_a; + +public: + Ber02ScaleAnalysisFilter(typename PFP::MAP& m, VertexAttribute& p, typename PFP::VEC3::DATA_TYPE a) : m_map(m), m_position(p), m_a(a) + {} + + void operator() () + { + TraversorV travV(m_map) ; + for (Dart d = travV.begin(); d != travV.end(); d = travV.next()) + { + if(m_map.isBoundaryVertex(d)) + m_position[d] /= m_a; + else + m_position[d] /= m_a * m_a; + } + + TraversorE travE(m_map) ; + for (Dart d = travE.begin(); d != travE.end(); d = travE.next()) + { + m_map.incCurrentLevel() ; + Dart midE = m_map.phi1(d); + if(!m_map.isBoundaryVertex(midE)) + m_position[midE] /= m_a ; + m_map.decCurrentLevel() ; + } + } +}; + + + +/************************************************************************************* + * With features preserving + *************************************************************************************/ + +// +// Synthesis +// + +//w-lift(a) +template +class Ber02OddSynthesisFilterFeatures : public Ber02OddSynthesisFilter +{ +protected: + const VertexAttribute& m_Vfeature; + const EdgeAttribute& m_Efeature; + +public: + Ber02OddSynthesisFilterFeatures(typename PFP::MAP& m, VertexAttribute& p, + const VertexAttribute v, const EdgeAttribute e, typename PFP::VEC3::DATA_TYPE a) : + Ber02OddSynthesisFilter(m, p, a), m_Vfeature(v), m_Efeature(e) + {} +} ; + +// s-lift(a) +template +class Ber02EvenSynthesisFilterFeatures : public Filter +{ +protected: + typename PFP::MAP& m_map ; + VertexAttribute& m_position ; + const VertexAttribute& m_Vfeature; + const EdgeAttribute& m_Efeature; + typename PFP::VEC3::DATA_TYPE m_a; + +public: + Ber02EvenSynthesisFilterFeatures(typename PFP::MAP& m, VertexAttribute& p, + const VertexAttribute v, const EdgeAttribute e, typename PFP::VEC3::DATA_TYPE a) : + m_map(m), m_position(p), m_Vfeature(v), m_Efeature(e), m_a(a) + {} + + void operator() () + { + TraversorV travV(m_map); + for(Dart d = travV.begin() ; d != travV.end() ; d = travV.next()) + { + if(m_Vfeature[d] == 2) //is a part of a face-feature + { + unsigned int count = 0; + typename PFP::VEC3 ev(0.0); + typename PFP::VEC3 fv(0.0); + Traversor2VF travVF(m_map,d); + for(Dart dit = travVF.begin(); dit != travVF.end() ; dit = travVF.next()) + { + m_map.incCurrentLevel() ; + + Dart midEdgeV = m_map.phi1(dit); + ev += m_position[midEdgeV]; + fv += m_position[m_map.phi1(midEdgeV)]; + + m_map.decCurrentLevel() ; + ++count; + } + fv /= count; + fv *= 4 * m_a * m_a; + + ev /= count; + ev *= 4 * m_a; + + m_position[d] += fv + ev; + } + else if(m_Vfeature[d] == 1) //is a part of an edge-feature + { + unsigned int count = 0; + typename PFP::VEC3 ev(0.0); + + Traversor2VF travVF(m_map,d); + for(Dart dit = travVF.begin(); dit != travVF.end() ; dit = travVF.next()) + { + m_map.incCurrentLevel() ; + + Dart midEdgeV = m_map.phi1(dit); + + if(m_Vfeature[midEdgeV] == 1) + { + ev += m_position[midEdgeV]; + ++count; + } + + m_map.decCurrentLevel() ; + } + + ev /= count; + ev *= 2 * m_a; + + m_position[d] += ev; + } + else + { + m_position[d] += 0.0; + } + } + + TraversorE travE(m_map); + for(Dart d = travE.begin() ; d != travE.end() ; d = travE.next()) + { + if(m_Efeature[d] == 2) + { + unsigned int count = 0; + + typename PFP::VEC3 fe(0); + Traversor2EF travEF(m_map, d); + for(Dart dit = travEF.begin() ; dit != travEF.end() ; dit = travEF.next()) + { + m_map.incCurrentLevel() ; + Dart midV = m_map.phi1(m_map.phi1(dit)); + fe += m_position[midV]; + m_map.decCurrentLevel() ; + ++count; + } + + fe /= count; + fe *= 2 * m_a; + + m_map.incCurrentLevel() ; + Dart midF = m_map.phi1(d); + m_position[midF] += fe; + m_map.decCurrentLevel() ; + } + else + { + m_map.incCurrentLevel() ; + Dart midF = m_map.phi1(d); + m_position[midF] += 0.0; + m_map.decCurrentLevel() ; + } + } + + } +} ; + +// s-scale(a) +template +class Ber02ScaleFilterFeatures : public Filter +{ +protected: + typename PFP::MAP& m_map ; + VertexAttribute& m_position ; + const VertexAttribute& m_Vfeature; + const EdgeAttribute& m_Efeature; + typename PFP::VEC3::DATA_TYPE m_a; + +public: + Ber02ScaleFilterFeatures(typename PFP::MAP& m, VertexAttribute& p, + const VertexAttribute v, const EdgeAttribute e, typename PFP::VEC3::DATA_TYPE a) : + m_map(m), m_position(p), m_Vfeature(v), m_Efeature(e), m_a(a) + {} + + void operator() () + { + TraversorV travV(m_map) ; + for (Dart d = travV.begin(); d != travV.end(); d = travV.next()) + { + m_position[d] *= pow(m_a, m_Vfeature[d]); + } + + TraversorE travE(m_map) ; + for (Dart d = travE.begin(); d != travE.end(); d = travE.next()) + { + if(m_Efeature[d] == 2) + { + m_map.incCurrentLevel() ; + Dart midE = m_map.phi1(d); + m_position[midE] *= m_a ; + m_map.decCurrentLevel() ; + } + } + } +} ; + +// +// Analysis +// + +//w-lift(a) +template +class Ber02OddAnalysisFilterFeatures : public Ber02OddAnalysisFilter +{ +protected: + const VertexAttribute& m_Vfeature; + const EdgeAttribute& m_Efeature; + +public: + Ber02OddAnalysisFilterFeatures(typename PFP::MAP& m, VertexAttribute& p, + const VertexAttribute v, const EdgeAttribute e, typename PFP::VEC3::DATA_TYPE a) : + Ber02OddAnalysisFilter(m, p, a), m_Vfeature(v), m_Efeature(e) + {} +} ; + + +} // namespace Filters + +} // namespace Primal + +} // namespace MR + +} // namespace Algo + +} // namespace CGoGN + + +#endif /* __3MR_FILTERS_PRIMAL__ */ diff --git a/include/Algo/Multiresolution/Map2MR/Filters/sqrt3.h b/include/Algo/Multiresolution/Map2MR/Filters/sqrt3.h index 026206588c0f93efb54a24bf7cf5dbe952709444..f2a3f54b46f59bf371286ab7b2551a6a4ef087ec 100644 --- a/include/Algo/Multiresolution/Map2MR/Filters/sqrt3.h +++ b/include/Algo/Multiresolution/Map2MR/Filters/sqrt3.h @@ -99,34 +99,6 @@ inline double omega0(unsigned int n) * Lazy Wavelet *********************************************************************************/ -template -class Sqrt3InitVertexSynthesisFilter : public Filter -{ -protected: - typename PFP::MAP& m_map ; - VertexAttribute& m_position ; - bool first; - -public: - Sqrt3InitVertexSynthesisFilter(typename PFP::MAP& m, VertexAttribute& p) : m_map(m), m_position(p), first(true) - {} - - void operator() () - { - if(first) - { - TraversorV trav(m_map) ; - for (Dart d = trav.begin(); d != trav.end(); d = trav.next()) - { - m_map.incCurrentLevel() ; - m_position[d] = typename PFP::VEC3(0.0); // ou phi2(d) - m_map.decCurrentLevel() ; - } - first = false; - } - } -} ; - template class Sqrt3FaceSynthesisFilter : public Filter { @@ -154,12 +126,19 @@ public: p1 *= 1.0 / 3.0 ; p2 *= 1.0 / 3.0 ; - m_map.incCurrentLevel() ; - - m_position[m_map.phi2(d)] += p0 + p1 + p2 ; + if(m_map.isBoundaryFace(d)) + { + Dart df = m_map.findBoundaryEdgeOfFace(d); + m_map.incCurrentLevel() ; + m_position[m_map.phi_1(m_map.phi2(df))] += p0 + p1 + p2 ; + } + else + { + m_map.incCurrentLevel() ; + m_position[m_map.phi2(d)] += p0 + p1 + p2 ; + } m_map.decCurrentLevel() ; - } } } ; @@ -180,36 +159,70 @@ public: TraversorV trav(m_map) ; for (Dart d = trav.begin(); d != trav.end(); d = trav.next()) { - typename PFP::VEC3 nf(0) ; - unsigned int degree = 0 ; + if(m_map.isBoundaryVertex(d)) + { + Dart df = m_map.findBoundaryEdgeOfVertex(d); - m_map.incCurrentLevel() ; - Dart df = m_map.phi2(m_map.phi1(d)); + if((m_map.getCurrentLevel()%2 == 0)) + { + typename PFP::VEC3 np(0) ; + typename PFP::VEC3 nl(0) ; + typename PFP::VEC3 nr(0) ; - Traversor2VVaE trav(m_map, df) ; - for(Dart it = trav.begin(); it != trav.end(); it = trav.next()) - { - ++degree ; - nf += m_position[it] ; + typename PFP::VEC3 pi = m_position[df]; + typename PFP::VEC3 pi_1 = m_position[m_map.phi_1(df)]; + typename PFP::VEC3 pi1 = m_position[m_map.phi1(df)]; + + np += pi_1 * 4 + pi * 19 + pi1 * 4; + np /= 27; + + nl += pi_1 * 10 + pi * 16 + pi1; + nl /= 27; + + nr += pi_1 + pi * 16 + pi1 * 10; + nr /= 27; + m_map.incCurrentLevel() ; + + m_position[df] = np; + m_position[m_map.phi_1(df)] = nl; + m_position[m_map.phi1(df)] = nr; + + m_map.decCurrentLevel() ; + } } - m_map.decCurrentLevel() ; + else + { + typename PFP::VEC3 nf(0) ; + unsigned int degree = 0 ; - float alpha = 1.0/9.0 * ( 4.0 - 2.0 * cos(2.0 * M_PI / degree)); - float teta = 1 - (3 * alpha) / 2; - float sigma = (3 * alpha) / (2 * degree); + m_map.incCurrentLevel() ; + Dart df = m_map.phi2(m_map.phi1(d)); - nf *= sigma; + Traversor2VVaE trav(m_map, df) ; + for(Dart it = trav.begin(); it != trav.end(); it = trav.next()) + { + ++degree ; + nf += m_position[it] ; - typename PFP::VEC3 vp = m_position[d] ; - vp *= teta ; + } + m_map.decCurrentLevel() ; - m_map.incCurrentLevel() ; + float alpha = 1.0/9.0 * ( 4.0 - 2.0 * cos(2.0 * M_PI / degree)); + float teta = 1 - (3 * alpha) / 2; + float sigma = (3 * alpha) / (2 * degree); - m_position[df] = vp + nf; + nf *= sigma; - m_map.decCurrentLevel() ; + typename PFP::VEC3 vp = m_position[d] ; + vp *= teta ; + + m_map.incCurrentLevel() ; + m_position[df] = vp + nf; + + m_map.decCurrentLevel() ; + } } } } ; diff --git a/include/Algo/Multiresolution/Map2MR/Masks/sqrt3.h b/include/Algo/Multiresolution/Map2MR/Masks/sqrt3.h index f7268717789b542de98a60698d5613516c377952..66aa28360d37755ae3ba2e1ba1aa9c5b16e1a81a 100644 --- a/include/Algo/Multiresolution/Map2MR/Masks/sqrt3.h +++ b/include/Algo/Multiresolution/Map2MR/Masks/sqrt3.h @@ -56,6 +56,44 @@ public: bool operator() (Dart d) { + if(m_map.isBoundaryVertex(d)) + { + Dart df = m_map.findBoundaryEdgeOfVertex(d); + + if(m_map.getDartLevel(df) < m_map.getCurrentLevel()) + { + if((m_map.getCurrentLevel()%2 == 0)) + { + m_map.decCurrentLevel() ; + + typename PFP::VEC3 np(0) ; + typename PFP::VEC3 nl(0) ; + typename PFP::VEC3 nr(0) ; + + typename PFP::VEC3 pi = m_position[df]; + typename PFP::VEC3 pi_1 = m_position[m_map.phi_1(df)]; + typename PFP::VEC3 pi1 = m_position[m_map.phi1(df)]; + + np += pi_1 * 4 + pi * 19 + pi1 * 4; + np /= 27; + + nl += pi_1 * 10 + pi * 16 + pi1; + nl /= 27; + + nr += pi_1 + pi * 16 + pi1 * 10; + nr /= 27; + + m_map.incCurrentLevel() ; + + m_position[df] = np; + m_position[m_map.phi_1(df)] = nl; + m_position[m_map.phi1(df)] = nr; + + return false; + } + } + } + m_map.decCurrentLevel() ; typename PFP::VEC3 np(0) ; @@ -97,8 +135,8 @@ public: m_map.decCurrentLevel() ; - Dart d1 = m_map.phi1(d1) ; - Dart d2 = m_map.phi1(d2) ; + Dart d1 = m_map.phi1(d0) ; + Dart d2 = m_map.phi1(d1) ; typename PFP::VEC3 p0 = m_position[d0] ; typename PFP::VEC3 p1 = m_position[d1] ; diff --git a/include/Algo/Multiresolution/Map2MR/map2MR_PrimalAdapt.h b/include/Algo/Multiresolution/Map2MR/map2MR_PrimalAdapt.h index 1f740fee4e9c58f8a8dfc75085430fc603444cc6..61c0a03a786d89021134ec52a207438a32cf68ba 100644 --- a/include/Algo/Multiresolution/Map2MR/map2MR_PrimalAdapt.h +++ b/include/Algo/Multiresolution/Map2MR/map2MR_PrimalAdapt.h @@ -131,6 +131,11 @@ protected: */ void splitFace(Dart d, Dart e) ; + /** + * + */ + void flipBackEdge(Dart d); + /*************************************************** * SUBDIVISION * ***************************************************/ @@ -151,6 +156,11 @@ public: */ unsigned int subdivideFace(Dart d, bool triQuad = true, bool OneLevelDifference = true); + /** + * + */ + unsigned int subdivideFaceSqrt3(Dart d); + /** * coarsen the face of d from the next level */ diff --git a/include/Algo/Multiresolution/Map2MR/map2MR_PrimalAdapt.hpp b/include/Algo/Multiresolution/Map2MR/map2MR_PrimalAdapt.hpp index 3b57ca6bc74e4c23323e7e1b690d352e9d2b5067..a2b453705a55f2c8ae9dba8f1ae7f2af21eb4664 100644 --- a/include/Algo/Multiresolution/Map2MR/map2MR_PrimalAdapt.hpp +++ b/include/Algo/Multiresolution/Map2MR/map2MR_PrimalAdapt.hpp @@ -302,6 +302,27 @@ void Map2MR::splitFace(Dart d, Dart e) m_map.splitFace(d, e) ; } +template +void Map2MR::flipBackEdge(Dart d) +{ + Dart dprev = m_map.phi_1(d) ; + Dart dnext = m_map.phi_1(d) ; + + Dart d2 = m_map.phi2(d); + Dart d2prev = m_map.phi_1(d2) ; + Dart d2next = m_map.phi_1(d2) ; + + m_map.duplicateDart(d); + m_map.duplicateDart(dprev); + m_map.duplicateDart(dnext); + + m_map.duplicateDart(d2); + m_map.duplicateDart(d2prev); + m_map.duplicateDart(d2next); + + m_map.flipBackEdge(d) ; +} + template void Map2MR::subdivideEdge(Dart d) { @@ -369,7 +390,7 @@ unsigned int Map2MR::subdivideFace(Dart d, bool triQuad, bool OneLevelDiffe m_map.setCurrentLevel(fLevel + 1) ; // go to the next level to perform face subdivision - if(triQuad & degree == 3) // if subdividing a triangle + if(triQuad && degree == 3) // if subdividing a triangle { Dart dd = m_map.phi1(old) ; Dart e = m_map.phi1(dd) ; @@ -468,6 +489,76 @@ void Map2MR::coarsenFace(Dart d) m_map.removeLevelBack() ; } +template +unsigned int Map2MR::subdivideFaceSqrt3(Dart d) +{ + assert(m_map.getDartLevel(d) <= m_map.getCurrentLevel() || !"subdivideFace : called with a dart inserted after current level") ; + //assert(!faceIsSubdivided(d) || !"Trying to subdivide an already subdivided face") ; + + unsigned int fLevel = faceLevel(d) ; + Dart old = faceOldestDart(d) ; + + m_map.pushLevel() ; + m_map.setCurrentLevel(fLevel) ; // go to the level of the face to subdivide its edges + + if(m_map.getCurrentLevel() == m_map.getMaxLevel()) + m_map.addLevelBack(); + + m_map.setCurrentLevel(fLevel + 1) ; // go to the next level to perform face subdivision + + //if it is an even level (triadic refinement) and a boundary face + if((m_map.getCurrentLevel()%2 == 0) && m_map.isBoundaryFace(d)) + { +// //find the boundary edge +// Dart df = m_map.findBoundaryEdgeOfFace(d); +// //trisection of the boundary edge +// cutEdge(df) ; +// splitFace(m_map.phi2(df), m_map.phi1(m_map.phi1(m_map.phi2(df)))) ; +// (*vertexVertexFunctor)(m_map.phi1(df) ; +// +// df = m_map.phi1(df); +// cutEdge(df) ; +// splitFace(m_map.phi2(df), m_map.phi1(m_map.phi1(m_map.phi2(df)))) ; +// //(*vertexVertexFunctor)(m_map.phi1(df)) ; + } + else + { + Dart d1 = m_map.phi1(old); + (*vertexVertexFunctor)(old) ; + splitFace(old, d1) ; + (*vertexVertexFunctor)(d1) ; + cutEdge(m_map.phi_1(old)) ; + Dart x = m_map.phi2(m_map.phi_1(old)) ; + Dart dd = m_map.phi1(m_map.phi1(m_map.phi1((x)))); + while(dd != x) + { + (*vertexVertexFunctor)(dd) ; + Dart next = m_map.phi1(dd) ; + splitFace(dd, m_map.phi1(x)) ; + dd = next ; + } + + Dart cd = m_map.phi2(x); + (*faceVertexFunctor)(cd) ; + + Dart dit = cd; + do + { + Dart dit12 = m_map.phi2(m_map.phi1(dit)); + + dit = m_map.phi2(m_map.phi_1(dit)); + + if(faceLevel(dit12) == (fLevel + 1) && !m_map.isBoundaryEdge(dit12)) + flipBackEdge(dit12); + } + while(dit != cd); + } + + m_map.popLevel() ; + + return fLevel ; +} + } // namespace Adaptive } // namespace Primal diff --git a/include/Algo/Multiresolution/Map2MR/map2MR_PrimalRegular.h b/include/Algo/Multiresolution/Map2MR/map2MR_PrimalRegular.h index 0805d0620799e65d3ebf8253e82da049cbca1b17..3c1b215ec5a65e51ae8d07abde9e8b972e078bb0 100644 --- a/include/Algo/Multiresolution/Map2MR/map2MR_PrimalRegular.h +++ b/include/Algo/Multiresolution/Map2MR/map2MR_PrimalRegular.h @@ -79,6 +79,10 @@ public: void analysis() ; void synthesis() ; + + //threshold + + void filtering(); } ; } // namespace Regular diff --git a/include/Algo/Multiresolution/Map2MR/map2MR_PrimalRegular.hpp b/include/Algo/Multiresolution/Map2MR/map2MR_PrimalRegular.hpp index d00c6959a178466a6dc93f7e1d4e7772ad6f8550..67c3ce58f888272e968f601e230ce060e1adcc46 100644 --- a/include/Algo/Multiresolution/Map2MR/map2MR_PrimalRegular.hpp +++ b/include/Algo/Multiresolution/Map2MR/map2MR_PrimalRegular.hpp @@ -163,9 +163,14 @@ void Map2MR::addNewLevelSqrt3(bool embedNewVertices) if((m_map.getCurrentLevel()%2 == 0) && m_map.isBoundaryFace(dit)) { //find the boundary edge - //Dart df = m_map.findBoundaryEdgeOfFace(dit); - + Dart df = m_map.findBoundaryEdgeOfFace(dit); //trisection of the boundary edge + m_map.cutEdge(df) ; + m_map.splitFace(m_map.phi2(df), m_map.phi1(m_map.phi1(m_map.phi2(df)))) ; + + df = m_map.phi1(df); + m_map.cutEdge(df) ; + m_map.splitFace(m_map.phi2(df), m_map.phi1(m_map.phi1(m_map.phi2(df)))) ; } else { @@ -199,7 +204,7 @@ void Map2MR::addNewLevelSqrt3(bool embedNewVertices) TraversorE te(m_map) ; for (Dart dit = te.begin(); dit != te.end(); dit = te.next()) { - if(m_map.getDartLevel(dit) < m_map.getCurrentLevel() && !m_map.isBoundaryEdge(dit)) + if(!m_map.isBoundaryEdge(dit) && m_map.getDartLevel(dit) < m_map.getCurrentLevel()) m_map.flipEdge(dit); } diff --git a/include/Algo/Multiresolution/Map3MR/Filters/bertram.h b/include/Algo/Multiresolution/Map3MR/Filters/bertram.h index 62690d750664979cc75d722ae0cb24df59719cd4..a297d877fcda819e978ffc7b96dcaccc13d598cc 100644 --- a/include/Algo/Multiresolution/Map3MR/Filters/bertram.h +++ b/include/Algo/Multiresolution/Map3MR/Filters/bertram.h @@ -45,6 +45,14 @@ namespace Primal namespace Filters { +/******************************************************************************* + * Without features preserving + *******************************************************************************/ + +// +// Synthesis +// + //w-lift(a) template class Ber02OddSynthesisFilter : public Filter @@ -52,85 +60,91 @@ class Ber02OddSynthesisFilter : public Filter protected: typename PFP::MAP& m_map ; VertexAttribute& m_position ; + typename PFP::VEC3::DATA_TYPE m_a; public: - Ber02OddSynthesisFilter(typename PFP::MAP& m, VertexAttribute& p) : m_map(m), m_position(p) + Ber02OddSynthesisFilter(typename PFP::MAP& m, VertexAttribute& p, typename PFP::VEC3::DATA_TYPE a) : m_map(m), m_position(p), m_a(a) {} void operator() () { - float a = 0.5; - TraversorW travW(m_map) ; for (Dart d = travW.begin(); d != travW.end(); d = travW.next()) { typename PFP::VEC3 vc = Algo::Geometry::volumeCentroid(m_map, d, m_position); + vc *= 8 * m_a * m_a * m_a; unsigned int count = 0; - typename PFP::VEC3 ec(0); + typename PFP::VEC3 ec(0.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)]; + ec += m_position[m_map.phi1(dit)]; m_map.decCurrentLevel(); ++count; } ec /= count; + ec *= 12 * m_a * m_a; count = 0; - typename PFP::VEC3 fc(0); + typename PFP::VEC3 fc(0.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))]; + fc += m_position[m_map.phi1(m_map.phi1(dit))]; m_map.decCurrentLevel(); ++count; } - fc /= count; + fc *= 6 * m_a; 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_position[midV] += vc + ec + 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 vf(0.0); + typename PFP::VEC3 ef(0.0); - 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()) { + vf += m_position[dit]; m_map.incCurrentLevel(); - ef += m_position[m_map.phi2(dit)]; + ef += m_position[m_map.phi1(dit)]; m_map.decCurrentLevel(); ++count; } ef /= count; + ef *= 4.0 * m_a; + + vf /= count; + vf *= 4.0 * m_a * m_a; m_map.incCurrentLevel() ; - Dart midF = m_map.phi2(m_map.phi1(d)); - m_position[midF] += vf * 4.0 * a * a + ef * 4.0 * a; + Dart midF = m_map.phi1(m_map.phi1(d)); + m_position[midF] += vf + ef ; 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); + typename PFP::VEC3 ve = (m_position[d] + m_position[m_map.phi1(d)]) * typename PFP::REAL(0.5); + ve *= 2.0 * m_a; m_map.incCurrentLevel() ; - Dart midV = m_map.phi2(d) ; - m_position[midV] += ve * a * 2.0; + Dart midV = m_map.phi1(d) ; + m_position[midV] += ve; m_map.decCurrentLevel() ; } - } } ; @@ -141,22 +155,74 @@ class Ber02EvenSynthesisFilter : public Filter protected: typename PFP::MAP& m_map ; VertexAttribute& m_position ; + typename PFP::VEC3::DATA_TYPE m_a; public: - Ber02EvenSynthesisFilter(typename PFP::MAP& m, VertexAttribute& p) : m_map(m), m_position(p) + Ber02EvenSynthesisFilter(typename PFP::MAP& m, VertexAttribute& p, typename PFP::VEC3::DATA_TYPE a) : m_map(m), m_position(p), m_a(a) {} 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)) + if(m_map.isBoundaryVertex(d)) + { + Dart db = m_map.findBoundaryFaceOfVertex(d); + + unsigned int count = 0; + typename PFP::VEC3 ev(0.0); + typename PFP::VEC3 fv(0.0); + + std::cout << "db = " << db << std::endl; + + Dart dit = db; + do + { + std::cout << "dit = " << dit << std::endl; + m_map.incCurrentLevel() ; + + Dart midEdgeV = m_map.phi1(dit); + ev += m_position[midEdgeV]; + fv += m_position[m_map.phi1(midEdgeV)]; + + m_map.decCurrentLevel() ; + ++count; + + dit = m_map.phi2(m_map.phi_1(dit)); + + }while(dit != db); + +// Traversor2VF travVF(m_map,db); +// for(Dart dit = travVF.begin(); dit != travVF.end() ; dit = travVF.next()) +// { +// std::cout << "dit = " << dit << std::endl; +// m_map.incCurrentLevel() ; +// +// Dart midEdgeV = m_map.phi1(dit); +// ev += m_position[midEdgeV]; +// fv += m_position[m_map.phi1(midEdgeV)]; +// +// m_map.decCurrentLevel() ; +// ++count; +// } + + std::cout << std::endl; + + fv /= count; + fv *= 4 * m_a * m_a; + + ev /= count; + ev *= 4 * m_a; + + m_position[d] += fv + ev; + + } + else { - typename PFP::VEC3 cv(0); unsigned int count = 0; + + typename PFP::VEC3 cv(0.0); Traversor3VW travVW(m_map,d); for(Dart dit = travVW.begin(); dit != travVW.end() ; dit = travVW.next()) { @@ -167,66 +233,37 @@ public: ++count; } cv /= count; + cv *= 8 * m_a * m_a * m_a; - typename PFP::VEC3 fv(0); + typename PFP::VEC3 fv(0.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)); + Dart midV = m_map.phi1(m_map.phi1(dit)); fv += m_position[midV]; m_map.decCurrentLevel() ; ++count; } fv /= count; + fv *= 12 * m_a * m_a; - typename PFP::VEC3 ev(0); + typename PFP::VEC3 ev(0.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 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); + Dart midV = m_map.phi1(dit); ev += m_position[midV]; m_map.decCurrentLevel() ; ++count; } ev /= count; + ev *= 6 * m_a; - m_position[db] += fv * 4 * a * a + ev * 4 * a; + m_position[d] += cv + fv + ev; } } @@ -237,29 +274,40 @@ public: { Dart db = m_map.findBoundaryFaceOfEdge(d); - typename PFP::VEC3 fe(0); + unsigned int count = 0; - m_map.incCurrentLevel() ; - Dart midV = m_map.phi2(m_map.phi1(db)); - fe += m_position[midV]; - m_map.decCurrentLevel() ; + typename PFP::VEC3 fe(0.0); m_map.incCurrentLevel() ; - midV = m_map.phi2(m_map.phi1(m_map.phi2(db))); + Dart midV = m_map.phi1(m_map.phi1(db)); + fe += m_position[midV]; + midV = m_map.phi_1(m_map.phi2(db)); fe += m_position[midV]; m_map.decCurrentLevel() ; +// Traversor2EF travEF(m_map, db); +// for(Dart dit = travEF.begin() ; dit != travEF.end() ; dit = travEF.next()) +// { +// m_map.incCurrentLevel() ; +// Dart midV = m_map.phi1(m_map.phi1(dit)); +// fe += m_position[midV]; +// m_map.decCurrentLevel() ; +// ++count; +// } + fe /= 2; + fe *= 2 * m_a; m_map.incCurrentLevel() ; - Dart midF = m_map.phi2(db); - m_position[midF] += fe * 2 * a; + Dart midE = m_map.phi1(db); + m_position[midE] += fe; m_map.decCurrentLevel() ; } else { - typename PFP::VEC3 ce(0); unsigned int count = 0; + + typename PFP::VEC3 ce(0.0); Traversor3EW travEW(m_map, d); for(Dart dit = travEW.begin() ; dit != travEW.end() ; dit = travEW.next()) { @@ -269,26 +317,26 @@ public: m_map.decCurrentLevel() ; ++count; } - ce /= count; + ce *= 4 * m_a * m_a; - typename PFP::VEC3 fe(0); + typename PFP::VEC3 fe(0.0); count = 0; - Traversor3FW travFW(m_map, d); - for(Dart dit = travFW.begin() ; dit != travFW.end() ; dit = travFW.next()) + Traversor3EF travEF(m_map, d); + for(Dart dit = travEF.begin() ; dit != travEF.end() ; dit = travEF.next()) { m_map.incCurrentLevel() ; - Dart midV = m_map.phi2(m_map.phi1(dit)); + Dart midV = m_map.phi1(m_map.phi1(dit)); fe += m_position[midV]; m_map.decCurrentLevel() ; ++count; } - fe /= count; + fe *= 4 * m_a; m_map.incCurrentLevel() ; - Dart midF = m_map.phi2(d); - m_position[midF] += ce * 4 * a * a + fe * 4 * a; + Dart midE = m_map.phi1(d); + m_position[midE] += ce + fe; m_map.decCurrentLevel() ; } } @@ -296,30 +344,29 @@ public: 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(); + unsigned int count = 0; - cf /= 2; - } + typename PFP::VEC3 cf(0.0); + Traversor3FW travFW(m_map, d); + for(Dart dit = travFW.begin() ; dit != travFW.end() ; dit = travFW.next()) + { + m_map.incCurrentLevel(); + Dart midV = m_map.phi_1(m_map.phi2(m_map.phi1(dit))); + cf += m_position[midV]; + m_map.decCurrentLevel(); + ++count; + } + cf /= count; + cf *= 2 * m_a; - m_map.incCurrentLevel() ; - Dart midF = m_map.phi2(m_map.phi1(d)); - m_position[midF] += cf * 2 * a; - m_map.decCurrentLevel() ; + m_map.incCurrentLevel() ; + Dart midF = m_map.phi1(m_map.phi1(d)); + m_position[midF] += cf; + m_map.decCurrentLevel() ; + } } - } } ; @@ -330,60 +377,418 @@ class Ber02ScaleSynthesisFilter : public Filter protected: typename PFP::MAP& m_map ; VertexAttribute& m_position ; + typename PFP::VEC3::DATA_TYPE m_a; public: - Ber02ScaleSynthesisFilter(typename PFP::MAP& m, VertexAttribute& p) : m_map(m), m_position(p) + Ber02ScaleSynthesisFilter(typename PFP::MAP& m, VertexAttribute& p, typename PFP::VEC3::DATA_TYPE a) : m_map(m), m_position(p), m_a(a) {} 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)) + m_position[d] *= m_a * m_a; + else + m_position[d] *= m_a *m_a * m_a; + } + + TraversorE travE(m_map) ; + for (Dart d = travE.begin(); d != travE.end(); d = travE.next()) + { + + m_map.incCurrentLevel() ; + Dart midE = m_map.phi1(d); + if(m_map.isBoundaryVertex(midE)) + m_position[midE] *= m_a; + else + m_position[midE] *= m_a * m_a; + m_map.decCurrentLevel() ; + + } + + TraversorF travF(m_map) ; + for (Dart d = travF.begin(); d != travF.end(); d = travF.next()) + { + m_map.incCurrentLevel() ; + Dart midF = m_map.phi1(m_map.phi1(d)); + if(!m_map.isBoundaryVertex(midF)) + m_position[midF] *= m_a ; + m_map.decCurrentLevel() ; + + } + } +} ; + +// +// Analysis +// + +//w-lift(a) +template +class Ber02OddAnalysisFilter : public Filter +{ +protected: + typename PFP::MAP& m_map ; + VertexAttribute& m_position ; + typename PFP::VEC3::DATA_TYPE m_a; + +public: + Ber02OddAnalysisFilter(typename PFP::MAP& m, VertexAttribute& p, typename PFP::VEC3::DATA_TYPE a) : m_map(m), m_position(p), m_a(a) + {} + + void operator() () + { + 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.phi1(d)]) * typename PFP::REAL(0.5); + ve *= 2.0 * m_a; + + m_map.incCurrentLevel() ; + Dart midV = m_map.phi1(d) ; + m_position[midV] -= ve; + m_map.decCurrentLevel() ; + } + + TraversorF travF(m_map) ; + for (Dart d = travF.begin(); d != travF.end(); d = travF.next()) + { + typename PFP::VEC3 vf(0.0); + typename PFP::VEC3 ef(0.0); + + unsigned int count = 0; + Traversor3FE travFE(m_map, d); + for (Dart dit = travFE.begin(); dit != travFE.end(); dit = travFE.next()) { - Dart db = m_map.findBoundaryFaceOfVertex(d); - m_position[db] *= a * a; + vf += m_position[dit]; + m_map.incCurrentLevel(); + ef += m_position[m_map.phi1(dit)]; + m_map.decCurrentLevel(); + ++count; } - else + ef /= count; + ef *= 4.0 * m_a; + + vf /= count; + vf *= 4.0 * m_a * m_a; + + m_map.incCurrentLevel() ; + Dart midF = m_map.phi1(m_map.phi1(d)); + m_position[midF] -= vf + ef ; + m_map.decCurrentLevel() ; + } + + TraversorW travW(m_map) ; + for (Dart d = travW.begin(); d != travW.end(); d = travW.next()) + { + typename PFP::VEC3 vc = Algo::Geometry::volumeCentroid(m_map, d, m_position); + vc *= 8 * m_a * m_a * m_a; + + unsigned int count = 0; + typename PFP::VEC3 ec(0.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.phi1(dit)]; + m_map.decCurrentLevel(); + ++count; + } + ec /= count; + ec *= 12 * m_a * m_a; + + count = 0; + typename PFP::VEC3 fc(0.0); + Traversor3WF travWF(m_map, d); + for (Dart dit = travWF.begin(); dit != travWF.end(); dit = travWF.next()) { - m_position[d] *= a * a * a; + m_map.incCurrentLevel(); + fc += m_position[m_map.phi1(m_map.phi1(dit))]; + m_map.decCurrentLevel(); + ++count; } + fc /= count; + fc *= 6 * m_a; + + m_map.incCurrentLevel() ; + Dart midV = m_map.phi_1(m_map.phi2(m_map.phi1(d))); + m_position[midV] -= vc + ec + fc; + m_map.decCurrentLevel() ; } + } +} ; - TraversorE travE(m_map) ; - for (Dart d = travE.begin(); d != travE.end(); d = travE.next()) +// s-lift(a) +template +class Ber02EvenAnalysisFilter : public Filter +{ +protected: + typename PFP::MAP& m_map ; + VertexAttribute& m_position ; + typename PFP::VEC3::DATA_TYPE m_a; + +public: + Ber02EvenAnalysisFilter(typename PFP::MAP& m, VertexAttribute& p, typename PFP::VEC3::DATA_TYPE a) : m_map(m), m_position(p), m_a(a) + {} + + void operator() () + { + TraversorF travF(m_map) ; + for (Dart d = travF.begin(); d != travF.end(); d = travF.next()) + { + if(!m_map.isBoundaryFace(d)) + { + unsigned int count = 0; + + typename PFP::VEC3 cf(0.0); + Traversor3FW travFW(m_map, d); + for(Dart dit = travFW.begin() ; dit != travFW.end() ; dit = travFW.next()) + { + m_map.incCurrentLevel(); + Dart midV = m_map.phi_1(m_map.phi2(m_map.phi1(dit))); + cf += m_position[midV]; + m_map.decCurrentLevel(); + ++count; + } + cf /= count; + cf *= 2 * m_a; + + m_map.incCurrentLevel() ; + Dart midF = m_map.phi1(m_map.phi1(d)); + m_position[midF] -= cf; + m_map.decCurrentLevel() ; + } + } + + 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); + unsigned int count = 0; + + typename PFP::VEC3 fe(0.0); + + m_map.incCurrentLevel() ; + Dart midV = m_map.phi1(m_map.phi1(db)); + fe += m_position[midV]; + midV = m_map.phi_1(m_map.phi2(db)); + fe += m_position[midV]; + m_map.decCurrentLevel() ; + +// Traversor2EF travEF(m_map, db); +// for(Dart dit = travEF.begin() ; dit != travEF.end() ; dit = travEF.next()) +// { +// m_map.incCurrentLevel() ; +// Dart midV = m_map.phi1(m_map.phi1(dit)); +// fe += m_position[midV]; +// m_map.decCurrentLevel() ; +// ++count; +// } + + fe /= 2; + fe *= 2 * m_a; + m_map.incCurrentLevel() ; - Dart midE = m_map.phi2(db); - m_position[midE] *= a ; + Dart midE = m_map.phi1(db); + m_position[midE] -= fe; m_map.decCurrentLevel() ; } else { + unsigned int count = 0; + + typename PFP::VEC3 ce(0.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; + ce *= 4 * m_a * m_a; + + typename PFP::VEC3 fe(0.0); + count = 0; + Traversor3EF travEF(m_map, d); + for(Dart dit = travEF.begin() ; dit != travEF.end() ; dit = travEF.next()) + { + m_map.incCurrentLevel() ; + Dart midV = m_map.phi1(m_map.phi1(dit)); + fe += m_position[midV]; + m_map.decCurrentLevel() ; + ++count; + } + fe /= count; + fe *= 4 * m_a; + m_map.incCurrentLevel() ; - Dart midE = m_map.phi2(d); - m_position[midE] *= a * a; + Dart midE = m_map.phi1(d); + m_position[midE] -= ce + fe; m_map.decCurrentLevel() ; } } + 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); + + unsigned int count = 0; + typename PFP::VEC3 ev(0.0); + typename PFP::VEC3 fv(0.0); + + std::cout << "db = " << db << std::endl; + + Dart dit = db; + do + { + std::cout << "dit = " << dit << std::endl; + m_map.incCurrentLevel() ; + + Dart midEdgeV = m_map.phi1(dit); + ev += m_position[midEdgeV]; + fv += m_position[m_map.phi1(midEdgeV)]; + + m_map.decCurrentLevel() ; + ++count; + + dit = m_map.phi2(m_map.phi_1(dit)); + + }while(dit != db); + +// Traversor2VF travVF(m_map,db); +// for(Dart dit = travVF.begin(); dit != travVF.end() ; dit = travVF.next()) +// { +// std::cout << "dit = " << dit << std::endl; +// m_map.incCurrentLevel() ; +// +// Dart midEdgeV = m_map.phi1(dit); +// ev += m_position[midEdgeV]; +// fv += m_position[m_map.phi1(midEdgeV)]; +// +// m_map.decCurrentLevel() ; +// ++count; +// } + + std::cout << std::endl; + + fv /= count; + fv *= 4 * m_a * m_a; + + ev /= count; + ev *= 4 * m_a; + + m_position[d] -= fv + ev; + + } + else + { + unsigned int count = 0; + + typename PFP::VEC3 cv(0.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; + cv *= 8 * m_a * m_a * m_a; + + typename PFP::VEC3 fv(0.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.phi1(m_map.phi1(dit)); + fv += m_position[midV]; + m_map.decCurrentLevel() ; + ++count; + } + fv /= count; + fv *= 12 * m_a * m_a; + + typename PFP::VEC3 ev(0.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.phi1(dit); + ev += m_position[midV]; + m_map.decCurrentLevel() ; + ++count; + } + ev /= count; + ev *= 6 * m_a; + + m_position[d] -= cv + fv + ev; + } + } + } +} ; + +// s-scale(a) +template +class Ber02ScaleAnalysisFilter : public Filter +{ +protected: + typename PFP::MAP& m_map ; + VertexAttribute& m_position ; + typename PFP::VEC3::DATA_TYPE m_a; + +public: + Ber02ScaleAnalysisFilter(typename PFP::MAP& m, VertexAttribute& p, typename PFP::VEC3::DATA_TYPE a) : m_map(m), m_position(p), m_a(a) + {} + + void operator() () + { 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 ; + Dart midF = m_map.phi1(m_map.phi1(d)); + if(!m_map.isBoundaryVertex(midF)) + m_position[midF] /= m_a ; m_map.decCurrentLevel() ; - } + + } + + TraversorE travE(m_map) ; + for (Dart d = travE.begin(); d != travE.end(); d = travE.next()) + { + + m_map.incCurrentLevel() ; + Dart midE = m_map.phi1(d); + if(m_map.isBoundaryVertex(midE)) + m_position[midE] /= m_a; + else + m_position[midE] /= m_a * m_a; + m_map.decCurrentLevel() ; + + } + + TraversorV travV(m_map) ; + for (Dart d = travV.begin(); d != travV.end(); d = travV.next()) + { + if(m_map.isBoundaryVertex(d)) + m_position[d] /= m_a * m_a; + else + m_position[d] /= m_a *m_a * m_a; } } } ; diff --git a/include/Algo/Multiresolution/Map3MR/Filters/lerp.h b/include/Algo/Multiresolution/Map3MR/Filters/lerp.h index c57f1174a1ade975f2d27837b7758aa59d88fbe0..447dcdacd3a82a0992316e56f48da8fe4b2ee506 100644 --- a/include/Algo/Multiresolution/Map3MR/Filters/lerp.h +++ b/include/Algo/Multiresolution/Map3MR/Filters/lerp.h @@ -98,11 +98,10 @@ public: typename PFP::VEC3 p = Algo::Geometry::faceCentroid(m_map, d, m_position); m_map.incCurrentLevel() ; - if(m_map.faceDegree(d) != 3) - { - Dart midF = m_map.phi2(m_map.phi1(d)); - m_position[midF] = p ; - } + + Dart midF = m_map.phi2(m_map.phi1(d)); + m_position[midF] = p ; + m_map.decCurrentLevel() ; } @@ -135,11 +134,36 @@ public: Dart midV = m_map.phi_1(m_map.phi2(m_map.phi1(d))); m_position[midV] = p ; } - else - { - Dart midV = m_map.phi_1(m_map.phi2(d)); - m_position[midV] = p; - } + + m_map.decCurrentLevel() ; + + } + } +} ; + +template +class LerpSqrt3VolumeSynthesisFilter : public Filter +{ +protected: + typename PFP::MAP& m_map ; + VertexAttribute& m_position ; + +public: + LerpSqrt3VolumeSynthesisFilter(typename PFP::MAP& m, VertexAttribute& p) : m_map(m), m_position(p) + {} + + void operator() () + { + TraversorW trav(m_map) ; + for (Dart d = trav.begin(); d != trav.end(); d = trav.next()) + { + typename PFP::VEC3 p = Algo::Geometry::volumeCentroid(m_map, d, m_position); + + m_map.incCurrentLevel() ; + + Dart midV = m_map.phi_1(m_map.phi2(d)); + m_position[midV] = p; + m_map.decCurrentLevel() ; } diff --git a/include/Algo/Multiresolution/Map3MR/map3MR_PrimalAdapt.h b/include/Algo/Multiresolution/Map3MR/map3MR_PrimalAdapt.h index 36fbf9f6eb2e382d92d3a4c7a462a45a75b3dfe7..0aa6ba81026bbe4c33f8a8a41c99485ad26ea7b3 100644 --- a/include/Algo/Multiresolution/Map3MR/map3MR_PrimalAdapt.h +++ b/include/Algo/Multiresolution/Map3MR/map3MR_PrimalAdapt.h @@ -58,11 +58,6 @@ public: typedef typename PFP::REAL REAL ; protected: - enum SubdivideType - { - S_TRI, - S_QUAD - } ; MAP& m_map; bool shareVertexEmbeddings; @@ -74,33 +69,6 @@ protected: public: Map3MR(MAP& map); -private: - /*! @name Topological helping functions - * - *************************************************************************/ - //@{ - //! - /*! - */ - void swapEdges(Dart d, Dart e); - - void splitSurfaceInVolume(std::vector& vd, bool firstSideClosed = true, bool secondSideClosed = false); - - Dart cutEdgeInVolume(Dart d); - - void splitFaceInVolume(Dart d, Dart e); - - void splitVolume(std::vector& vd); - - void saveRelationsAroundVertex(Dart d, std::vector >& vd); - - void unsewAroundVertex(std::vector >& vd); - - Dart cutEdge(Dart d) ; - void splitFace(Dart d, Dart e) ; - //@} - -public: /*! @name Cells informations * *************************************************************************/ @@ -149,33 +117,58 @@ public: bool volumeIsSubdivided(Dart d); //@} - /*! @name Subdivision - * - *************************************************************************/ protected: - //@{ - //! Subdivide the edge of d to the next level - /*! @param d Dart from the edge - */ - void subdivideEdge(Dart d) ; - - //! Subdivide the edge of d to the next level - /*! @param d Dart frome the face - */ - void subdivideFace(Dart d, SubdivideType sType = S_QUAD) ; +// /*! @name Topological helping functions +// * +// *************************************************************************/ +// //@{ +// //! +// /*! +// */ +// void swapEdges(Dart d, Dart e); +// +// void splitSurfaceInVolume(std::vector& vd, bool firstSideClosed = true, bool secondSideClosed = false); +// +// Dart cutEdgeInVolume(Dart d); +// +// void splitFaceInVolume(Dart d, Dart e); +// +// void splitVolume(std::vector& vd); +// +// void saveRelationsAroundVertex(Dart d, std::vector >& vd); +// +// void unsewAroundVertex(std::vector >& vd); +// +// Dart cutEdge(Dart d) ; +// void splitFace(Dart d, Dart e) ; +// //@} +// +// /*! @name Subdivision +// * +// *************************************************************************/ +// //@{ +// //! Subdivide the edge of d to the next level +// /*! @param d Dart from the edge +// */ +// void subdivideEdge(Dart d) ; +// +// //! Subdivide the edge of d to the next level +// /*! @param d Dart frome the face +// */ +// void subdivideFace(Dart d, bool triQuad) ; public: - //! Subdivide the volume of d to hexahedral cells - /*! @param d Dart from the volume - */ - void subdivideVolume(Dart d) ; - - //! Subdivide the volume of d to hexahedral cells - /*! @param d Dart from the volume - */ - void subdivideVolumeTetOcta(Dart d) ; - - void subdivideVolumeTetOctaTemp(Dart d); +// //! Subdivide the volume of d to hexahedral cells +// /*! @param d Dart from the volume +// */ +// unsigned int subdivideVolume(Dart d, bool triQuad = true, bool OneLevelDifference = true); +// +// //! Subdivide the volume of d to hexahedral cells +// /*! @param d Dart from the volume +// */ +// void subdivideVolumeTetOcta(Dart d) ; +// +// void subdivideVolumeTetOctaTemp(Dart d); //@} /*! @name Vertices Attributes management diff --git a/include/Algo/Multiresolution/Map3MR/map3MR_PrimalAdapt.hpp b/include/Algo/Multiresolution/Map3MR/map3MR_PrimalAdapt.hpp index 9f6f4f094650aad48a665f5dcb1f6897b8e7b11f..5e936f8c9c1a87a778683f83f313815af5e40361 100644 --- a/include/Algo/Multiresolution/Map3MR/map3MR_PrimalAdapt.hpp +++ b/include/Algo/Multiresolution/Map3MR/map3MR_PrimalAdapt.hpp @@ -22,7 +22,7 @@ * * *******************************************************************************/ -#include "Algo/Multiresolution/map3MR/map3MR_PrimalAdapt.h" +#include "Algo/Multiresolution/Map3MR/map3MR_PrimalAdapt.h" namespace CGoGN { @@ -42,7 +42,7 @@ namespace Adaptive template Map3MR::Map3MR(typename PFP::MAP& map) : m_map(map), - shareVertexEmbeddings(false), + shareVertexEmbeddings(true), vertexVertexFunctor(NULL), edgeVertexFunctor(NULL), faceVertexFunctor(NULL), @@ -51,245 +51,6 @@ Map3MR::Map3MR(typename PFP::MAP& map) : } -/*! @name Topological helping functions - *************************************************************************/ -template -void Map3MR::swapEdges(Dart d, Dart e) -{ - if(!m_map.PFP::MAP::ParentMap::isBoundaryEdge(d) && !m_map.PFP::MAP::ParentMap::isBoundaryEdge(e)) - { - Dart d2 = m_map.phi2(d); - Dart e2 = m_map.phi2(e); - -// Map2::swapEdges(d,e); - - m_map.PFP::MAP::ParentMap::unsewFaces(d); - m_map.PFP::MAP::ParentMap::unsewFaces(e); - - m_map.PFP::MAP::ParentMap::sewFaces(d, e); - m_map.PFP::MAP::ParentMap::sewFaces(d2, e2); - - if(m_map.template isOrbitEmbedded()) - { - m_map.template copyDartEmbedding(d, m_map.phi2(m_map.phi_1(d))); - m_map.template copyDartEmbedding(e, m_map.phi2(m_map.phi_1(e))); - m_map.template copyDartEmbedding(d2, m_map.phi2(m_map.phi_1(d2))); - m_map.template copyDartEmbedding(e2, m_map.phi2(m_map.phi_1(e2))); - } - - if(m_map.template isOrbitEmbedded()) - { - - } - - if(m_map.template isOrbitEmbedded()) - m_map.template embedNewCell(d); - - -// propagateDartRelation(d, m_phi2) ; -// propagateDartRelation(d, m_phi3) ; -// //propagateDartRelation(d2, m_phi2) ; -// //propagateDartRelation(d2, m_phi3) ; -// propagateDartRelation(e, m_phi2) ; -// propagateDartRelation(e, m_phi3) ; -// //propagateDartRelation(e2, m_phi2) ; -// //propagateDartRelation(e2, m_phi3) ; - } -} - -template -void Map3MR::splitSurfaceInVolume(std::vector& vd, bool firstSideClosed, bool secondSideClosed) -{ - std::vector vd2 ; - vd2.reserve(vd.size()); - - // save the edge neighbors darts - for(std::vector::iterator it = vd.begin() ; it != vd.end() ; ++it) - { - vd2.push_back(phi2(*it)); - } - - assert(vd2.size() == vd.size()); - - Map2::splitSurface(vd, firstSideClosed, secondSideClosed); - - // follow the edge path a second time to embed the vertex, edge and volume orbits - for(unsigned int i = 0; i < vd.size(); ++i) - { - Dart dit = vd[i]; - Dart dit2 = vd2[i]; - -// propagateDartRelation(dit, m_phi1) ; -// propagateDartRelation(dit, m_phi_1) ; -// propagateDartRelation(dit, m_phi2) ; -// -// propagateDartRelation(dit2, m_phi1) ; -// propagateDartRelation(dit2, m_phi_1) ; -// propagateDartRelation(dit2, m_phi2) ; - - // embed the vertex embedded from the origin volume to the new darts - if(isOrbitEmbedded()) - { - copyDartEmbedding(phi2(dit), phi1(dit)); - copyDartEmbedding(phi2(dit2), phi1(dit2)); - } - } - -} - -template -void Map3MR::splitFaceInVolume(Dart d, Dart e) -{ - Dart dprev = phi_1(d) ; - Dart eprev = phi_1(e) ; - Map2::splitFace(d, e) ; - Dart dd = phi1(dprev) ; - Dart ee = phi1(eprev) ; - - propagateDartRelation(d, m_phi_1) ; - propagateDartRelation(e, m_phi_1) ; - propagateDartRelation(dd, m_phi1) ; - propagateDartRelation(dd, m_phi_1) ; - propagateDartRelation(dd, m_phi2) ; - propagateDartRelation(ee, m_phi1) ; - propagateDartRelation(ee, m_phi_1) ; - propagateDartRelation(ee, m_phi2) ; - propagateDartRelation(dprev, m_phi1) ; - propagateDartRelation(eprev, m_phi1) ; - - propagateDartEmbedding(dd) ; - propagateDartEmbedding(ee) ; -} - -template -Dart Map3MR::cutEdgeInVolume(Dart d) -{ - Dart dd = phi2(d) ; - Dart d1 = Map2::cutEdge(d) ; - Dart dd1 = phi1(dd) ; - Dart d11 = phi1(d1) ; - Dart dd11 = phi1(dd1) ; - - propagateDartRelation(d, m_phi1) ; - propagateDartRelation(d, m_phi2) ; - propagateDartRelation(dd, m_phi1) ; - propagateDartRelation(dd, m_phi2) ; - propagateDartRelation(d1, m_phi1) ; - propagateDartRelation(d1, m_phi_1) ; - propagateDartRelation(d1, m_phi2) ; - propagateDartRelation(dd1, m_phi1) ; - propagateDartRelation(dd1, m_phi_1) ; - propagateDartRelation(dd1, m_phi2) ; - propagateDartRelation(d11, m_phi_1) ; - propagateDartRelation(dd11, m_phi_1) ; - - return d1 ; -} - - -template -Dart Map3MR::cutEdge(Dart d) -{ - Dart dit = d; - do - { - Dart dd = phi2(dit) ; - Dart d1 = phi1(dit) ; - Dart dd1 = phi1(dd) ; - - Dart dd3 = phi3(dit); - Dart d13 = phi2(dd3) ; - Dart d13 = phi1(d13) ; - Dart dd13 = phi1(dd) ; - - - m_map.duplicateDart(d); - m_map.duplicateDart(dd); - m_map.duplicateDart(d1); - m_map.duplicateDart(dd1); - - m_map.duplicateDart(d); - m_map.duplicateDart(dd); - m_map.duplicateDart(d1); - m_map.duplicateDart(dd1); - - - dit = m_map.alpha2(dit); - }while(dit != d); - - Dart nd = EmbeddedMap3::cutEdge(d); - - return nd; -} - -template -void Map3MR::splitFace(Dart d, Dart e) -{ - Dart dprev = phi_1(d) ; - Dart eprev = phi_1(e) ; - EmbeddedMap3::splitFace(d,e); - Dart dd = phi1(dprev) ; - Dart ee = phi1(eprev) ; - - propagateDartRelation(d, m_phi_1) ; - propagateDartRelation(e, m_phi_1) ; - propagateDartRelation(dd, m_phi1) ; - propagateDartRelation(dd, m_phi_1) ; - propagateDartRelation(dd, m_phi2) ; - propagateDartRelation(ee, m_phi1) ; - propagateDartRelation(ee, m_phi_1) ; - propagateDartRelation(ee, m_phi2) ; - propagateDartRelation(dprev, m_phi1) ; - propagateDartRelation(eprev, m_phi1) ; - - propagateDartRelation(phi3(dprev), m_phi_1) ; - propagateDartRelation(phi3(eprev), m_phi_1) ; - propagateDartRelation(phi3(dd), m_phi1) ; - propagateDartRelation(phi3(dd), m_phi_1) ; - propagateDartRelation(phi3(dd), m_phi2) ; - propagateDartRelation(phi3(ee), m_phi1) ; - propagateDartRelation(phi3(ee), m_phi_1) ; - propagateDartRelation(phi3(ee), m_phi2) ; - propagateDartRelation(phi3(d), m_phi1) ; - propagateDartRelation(phi3(e), m_phi1) ; - - propagateDartEmbedding(dd) ; - propagateDartEmbedding(ee) ; - propagateDartEmbedding(phi3(dd)) ; - propagateDartEmbedding(phi3(ee)) ; -} - -template -void Map3MR::splitVolume(std::vector& vd) -{ - EmbeddedMap3::splitVolume(vd); - - for(std::vector::iterator it = vd.begin() ; it != vd.end() ; ++it) - { - Dart dit = *it; - propagateDartRelation(dit, m_phi2) ; - - Dart dit2 = phi2(dit); - propagateDartRelation(dit2, m_phi1) ; - propagateDartRelation(dit2, m_phi_1) ; - propagateDartRelation(dit2, m_phi2) ; - propagateDartRelation(dit2, m_phi3) ; - - Dart dit23 = phi3(dit2); - propagateDartRelation(dit23, m_phi1) ; - propagateDartRelation(dit23, m_phi_1) ; - propagateDartRelation(dit23, m_phi2) ; - propagateDartRelation(dit23, m_phi3) ; - - Dart dit232 = phi2(dit23); - propagateDartRelation(dit232, m_phi2) ; - - propagateDartEmbedding(dit2) ; - propagateDartEmbedding(dit23) ; - } - -} - /*! @name Cells informations *************************************************************************/ template @@ -297,15 +58,17 @@ unsigned int Map3MR::edgeLevel(Dart d) { assert(m_map.getDartLevel(d) <= m_map.getCurrentLevel() || !"edgeLevel : called with a dart inserted after current level") ; + // the level of an edge is the maximum of the + // insertion levels of its darts unsigned int r = 0; Dart dit = d; do { unsigned int ld = m_map.getDartLevel(dit) ; unsigned int ldd = m_map.getDartLevel(m_map.phi2(dit)) ; - unsigned int min = ld < ldd ? ldd : ld; + unsigned int max = ld < ldd ? ldd : ld; - r = r < min ? min : r ; + r = r < max ? max : r ; dit = m_map.alpha2(dit); } while(dit != d); @@ -503,6 +266,216 @@ bool Map3MR::volumeIsSubdivided(Dart d) return subd; } +/*! @name Topological helping functions + *************************************************************************/ +template +void Map3MR::swapEdges(Dart d, Dart e) +{ + if(!m_map.PFP::MAP::ParentMap::isBoundaryEdge(d) && !m_map.PFP::MAP::ParentMap::isBoundaryEdge(e)) + { + Dart d2 = m_map.phi2(d); + Dart e2 = m_map.phi2(e); + +// Map2::swapEdges(d,e); + + m_map.PFP::MAP::ParentMap::unsewFaces(d); + m_map.PFP::MAP::ParentMap::unsewFaces(e); + + m_map.PFP::MAP::ParentMap::sewFaces(d, e); + m_map.PFP::MAP::ParentMap::sewFaces(d2, e2); + + if(m_map.template isOrbitEmbedded()) + { + m_map.template copyDartEmbedding(d, m_map.phi2(m_map.phi_1(d))); + m_map.template copyDartEmbedding(e, m_map.phi2(m_map.phi_1(e))); + m_map.template copyDartEmbedding(d2, m_map.phi2(m_map.phi_1(d2))); + m_map.template copyDartEmbedding(e2, m_map.phi2(m_map.phi_1(e2))); + } + + if(m_map.template isOrbitEmbedded()) + { + + } + + if(m_map.template isOrbitEmbedded()) + m_map.template embedNewCell(d); + + +// propagateDartRelation(d, m_phi2) ; +// propagateDartRelation(d, m_phi3) ; +// //propagateDartRelation(d2, m_phi2) ; +// //propagateDartRelation(d2, m_phi3) ; +// propagateDartRelation(e, m_phi2) ; +// propagateDartRelation(e, m_phi3) ; +// //propagateDartRelation(e2, m_phi2) ; +// //propagateDartRelation(e2, m_phi3) ; + } +} + +template +void Map3MR::splitSurfaceInVolume(std::vector& vd, bool firstSideClosed, bool secondSideClosed) +{ + std::vector vd2 ; + vd2.reserve(vd.size()); + + // save the edge neighbors darts + for(std::vector::iterator it = vd.begin() ; it != vd.end() ; ++it) + { + vd2.push_back(phi2(*it)); + } + + assert(vd2.size() == vd.size()); + + Map2::splitSurface(vd, firstSideClosed, secondSideClosed); + + // follow the edge path a second time to embed the vertex, edge and volume orbits + for(unsigned int i = 0; i < vd.size(); ++i) + { + Dart dit = vd[i]; + Dart dit2 = vd2[i]; + +// propagateDartRelation(dit, m_phi1) ; +// propagateDartRelation(dit, m_phi_1) ; +// propagateDartRelation(dit, m_phi2) ; +// +// propagateDartRelation(dit2, m_phi1) ; +// propagateDartRelation(dit2, m_phi_1) ; +// propagateDartRelation(dit2, m_phi2) ; + + // embed the vertex embedded from the origin volume to the new darts + if(isOrbitEmbedded()) + { + copyDartEmbedding(phi2(dit), phi1(dit)); + copyDartEmbedding(phi2(dit2), phi1(dit2)); + } + } + +} + +template +void Map3MR::splitFaceInVolume(Dart d, Dart e) +{ + Dart dprev = phi_1(d) ; + Dart eprev = phi_1(e) ; + Map2::splitFace(d, e) ; + Dart dd = phi1(dprev) ; + Dart ee = phi1(eprev) ; + + propagateDartRelation(d, m_phi_1) ; + propagateDartRelation(e, m_phi_1) ; + propagateDartRelation(dd, m_phi1) ; + propagateDartRelation(dd, m_phi_1) ; + propagateDartRelation(dd, m_phi2) ; + propagateDartRelation(ee, m_phi1) ; + propagateDartRelation(ee, m_phi_1) ; + propagateDartRelation(ee, m_phi2) ; + propagateDartRelation(dprev, m_phi1) ; + propagateDartRelation(eprev, m_phi1) ; + + propagateDartEmbedding(dd) ; + propagateDartEmbedding(ee) ; +} + +template +Dart Map3MR::cutEdgeInVolume(Dart d) +{ + Dart dd = phi2(d) ; + Dart d1 = Map2::cutEdge(d) ; + Dart dd1 = phi1(dd) ; + Dart d11 = phi1(d1) ; + Dart dd11 = phi1(dd1) ; + + propagateDartRelation(d, m_phi1) ; + propagateDartRelation(d, m_phi2) ; + propagateDartRelation(dd, m_phi1) ; + propagateDartRelation(dd, m_phi2) ; + propagateDartRelation(d1, m_phi1) ; + propagateDartRelation(d1, m_phi_1) ; + propagateDartRelation(d1, m_phi2) ; + propagateDartRelation(dd1, m_phi1) ; + propagateDartRelation(dd1, m_phi_1) ; + propagateDartRelation(dd1, m_phi2) ; + propagateDartRelation(d11, m_phi_1) ; + propagateDartRelation(dd11, m_phi_1) ; + + return d1 ; +} + + +template +Dart Map3MR::cutEdge(Dart d) +{ + Dart dit = d; + do + { + Dart dd = m_map.phi2(dit) ; + Dart d1 = m_map.phi1(dit); + Dart dd1 = m_map.phi1(dd); + + m_map.duplicateDart(dit); + m_map.duplicateDart(dd); + m_map.duplicateDart(d1); + m_map.duplicateDart(dd1); + + dit = m_map.alpha2(dit); + }while(dit != d); + + Dart nd = m_map.cutEdge(d); + + return nd; +} + +template +void Map3MR::splitFace(Dart d, Dart e) +{ + Dart dprev = m_map.phi_1(d) ; + Dart eprev = m_map.phi_1(e) ; + + m_map.duplicateDart(d); + m_map.duplicateDart(e); + m_map.duplicateDart(dprev); + m_map.duplicateDart(eprev); + + m_map.duplicateDart(m_map.phi3(d)); + m_map.duplicateDart(m_map.phi3(e)); + m_map.duplicateDart(m_map.phi3(dprev)); + m_map.duplicateDart(m_map.phi3(eprev)); + + m_map.splitFace(d,e); +} + +template +void Map3MR::splitVolume(std::vector& vd) +{ + EmbeddedMap3::splitVolume(vd); + + for(std::vector::iterator it = vd.begin() ; it != vd.end() ; ++it) + { + Dart dit = *it; + propagateDartRelation(dit, m_phi2) ; + + Dart dit2 = phi2(dit); + propagateDartRelation(dit2, m_phi1) ; + propagateDartRelation(dit2, m_phi_1) ; + propagateDartRelation(dit2, m_phi2) ; + propagateDartRelation(dit2, m_phi3) ; + + Dart dit23 = phi3(dit2); + propagateDartRelation(dit23, m_phi1) ; + propagateDartRelation(dit23, m_phi_1) ; + propagateDartRelation(dit23, m_phi2) ; + propagateDartRelation(dit23, m_phi3) ; + + Dart dit232 = phi2(dit23); + propagateDartRelation(dit232, m_phi2) ; + + propagateDartEmbedding(dit2) ; + propagateDartEmbedding(dit23) ; + } + +} + + /*! @name Subdivision *************************************************************************/ template @@ -511,6 +484,8 @@ void Map3MR::subdivideEdge(Dart d) assert(m_map.getDartLevel(d) <= m_map.getCurrentLevel() || !"subdivideEdge : called with a dart inserted after current level") ; assert(!edgeIsSubdivided(d) || !"Trying to subdivide an already subdivided edge") ; + assert(m_map.getCurrentLevel() == edgeLevel(d) || !"Trying to subdivide an edge on a bad current level") ; + m_map.incCurrentLevel(); Dart nd = cutEdge(d); @@ -521,7 +496,7 @@ void Map3MR::subdivideEdge(Dart d) } template -void Map3MR::subdivideFace(Dart d, SubdivideType sType) +void Map3MR::subdivideFace(Dart d, bool triQuad) { assert(m_map.getDartLevel(d) <= m_map.getCurrentLevel() || !"subdivideFace : called with a dart inserted after current level") ; assert(!faceIsSubdivided(d) || !"Trying to subdivide an already subdivided face") ; @@ -545,7 +520,7 @@ void Map3MR::subdivideFace(Dart d, SubdivideType sType) m_map.setCurrentLevel(fLevel + 1) ; // go to the next level to perform face subdivision - if(degree == 3 && sType == S_TRI) // if subdividing a triangle + if(triQuad && degree == 3) // if subdividing a triangle { Dart dd = m_map.phi1(old) ; Dart e = m_map.phi1(dd) ; @@ -594,7 +569,7 @@ void Map3MR::subdivideFace(Dart d, SubdivideType sType) } template -void Map3MR::subdivideVolume(Dart d) +unsigned int Map3MR::subdivideVolume(Dart d, bool triQuad, bool OneLevelDifference) { assert(m_map.getDartLevel(d) <= m_map.getCurrentLevel() || !"subdivideVolume : called with a dart inserted after current level") ; assert(!volumeIsSubdivided(d) || !"Trying to subdivide an already subdivided face") ; @@ -616,7 +591,7 @@ void Map3MR::subdivideVolume(Dart d) { //if needed subdivide face if(!faceIsSubdivided(dit)) - subdivideFace(dit); + subdivideFace(dit,triQuad); } // @@ -626,12 +601,10 @@ void Map3MR::subdivideVolume(Dart d) std::vector > subdividedFaces; subdividedFaces.reserve(128); - Traversor3WV traWV(m_map, old); + Traversor3WV traWV(m_map, dit); for(Dart ditWV = traWV.begin(); ditWV != traWV.end(); ditWV = traWV.next()) { - m_map.setCurrentLevel(m_map.getMaxLevel()) ; - m_map.template embedOrbit(ditWV, EMBNULL); - (*vertexVertexFunctor)(ditWV) ; + m_map.setCurrentLevel(m_map.getMaxLevel()) ; //Utile ? Dart e = ditWV; std::vector v ; @@ -641,10 +614,10 @@ void Map3MR::subdivideVolume(Dart d) v.push_back(m_map.phi1(e)); v.push_back(m_map.phi1(m_map.phi1(e))); - if(!Map2::isBoundaryEdge(m_map.phi1(e))) + if(!m_map.PFP::MAP::ParentMap::isBoundaryEdge(m_map.phi1(e))) subdividedFaces.push_back(std::pair(m_map.phi1(e),m_map.phi2(m_map.phi1(e)))); - if(!Map2::isBoundaryEdge(m_map.phi1(m_map.phi1(e)))) + if(!m_map.PFP::MAP::ParentMap::isBoundaryEdge(m_map.phi1(m_map.phi1(e)))) subdividedFaces.push_back(std::pair(m_map.phi1(m_map.phi1(e)),m_map.phi2(m_map.phi1(m_map.phi1(e))))); e = m_map.phi2(m_map.phi_1(e)); @@ -655,32 +628,28 @@ void Map3MR::subdivideVolume(Dart d) Dart dd = m_map.phi2(m_map.phi1(ditWV)); Dart next = m_map.phi1(m_map.phi1(dd)) ; - //Map2::splitFace(dd, next) ; - splitFaceInVolume(dd, next); + m_map.PFP::MAP::ParentMap::splitFace(dd, next) ; Dart ne = m_map.phi2(m_map.phi_1(dd)); - //Map2::cutEdge(ne) ; - cutEdgeInVolume(ne); + m_map.PFP::MAP::ParentMap::cutEdge(ne) ; centralDart = m_map.phi1(ne); dd = m_map.phi1(m_map.phi1(next)) ; while(dd != ne) { Dart tmp = m_map.phi1(ne) ; - //Map2::splitFace(tmp, dd) ; - splitFaceInVolume(tmp, dd); + m_map.PFP::MAP::ParentMap::splitFace(tmp, dd) ; dd = m_map.phi1(m_map.phi1(dd)) ; } - m_map.setCurrentLevel(m_map.getMaxLevel() - 1) ; + m_map.setCurrentLevel(m_map.getMaxLevel() - 1) ; //Utile ? } - std::cout << "nb Vertices incident to volume : " << i << std::endl; - m_map.setCurrentLevel(m_map.getMaxLevel()) ; DartMarkerNoUnmark mf(*this); - std::cout << "subdFaces size = " << subdividedFaces.size() << std::endl; + m_map.setCurrentLevel(m_map.getMaxLevel()) ; + //4 couture des relations precedemment sauvegarde for (std::vector >::iterator it = subdividedFaces.begin(); it != subdividedFaces.end(); ++it) { Dart f1 = m_map.phi2((*it).first); @@ -688,25 +657,11 @@ void Map3MR::subdivideVolume(Dart d) //if(isBoundaryFace(f1) && isBoundaryFace(f2)) if(m_map.phi3(f1) == f1 && m_map.phi3(f2) == f2) - //if(!mf.isMarked(phi3(f1)) && !mf.isMarked(phi3(f2))) - { m_map.sewVolumes(f1, f2, false); - mf.markOrbit(f1); - - Dart temp = f1; - do - { - propagateDartRelation(temp, m_phi3) ; - propagateDartRelation(phi3(temp), m_phi3) ; - temp = m_map.phi1(temp); - }while(temp != f1); - } } - - m_map.template embedOrbit(centralDart, EMBNULL); - (*volumeVertexFunctor)(centralDart) ; m_map.template embedOrbit(centralDart, m_map.template getEmbedding(centralDart)); - propagateOrbitEmbedding(centralDart) ; + + m_map.setCurrentLevel(m_map.getMaxLevel() - 1) ; //A optimiser @@ -868,67 +823,6 @@ void Map3MR::subdivideVolumeTetOcta(Dart d) popLevel(); } -//void Map3MR_PrimalAdapt::subdivideVolumeTetOctaTemp(Dart d) -//{ -// assert(getDartLevel(d) <= getCurrentLevel() || !"subdivideVolumeTetOcta : called with a dart inserted after current level") ; -// assert(!volumeIsSubdivided(d) || !"Trying to subdivide an already subdivided face") ; -// -// unsigned int vLevel = volumeLevel(d); -// Dart old = volumeOldestDart(d); -// -// pushLevel() ; -// setCurrentLevel(vLevel) ; // go to the level of the face to subdivide its edges -// -// if(getCurrentLevel() == getMaxLevel()) -// addNewLevel() ; -// -// unsigned int j = 0; -// -// // -// // Subdivide Faces and Edges -// // -// Traversor3WF traF(*this, old); -// for(Dart dit = traF.begin(); dit != traF.end(); dit = traF.next()) -// { -// std::cout << "CurrentLevel = " << getCurrentLevel() << std::endl; -// std::cout << "face level = " << faceLevel(dit) << std::endl; -// -// //if needed subdivide face -// if(!faceIsSubdivided(dit)) -// { -// std::cout << "subdivide face = " << dit << std::endl; -// subdivideFace(dit, S_TRI); -// ++j; -// } -// } -// -// // -// // Create inside volumes -// // -// Dart centralDart = NIL; -// bool isNotTet = false; -// Traversor3WV traV(*this, old); -// setCurrentLevel(vLevel + 1) ; -// for(Dart dit = traV.begin(); dit != traV.end(); dit = traV.next()) -// { -// Dart f1 = phi1(dit); -// Dart e = dit; -// std::vector v ; -// -// do -// { -// v.push_back(phi1(e)); -// e = phi2(phi_1(e)); -// } -// while(e != dit); -// -// std::cout << "v size = " << v.size() << std::endl; -// -// } -// -// -// popLevel(); -//} } // namespace Adaptive diff --git a/include/Algo/Multiresolution/Map3MR/map3MR_PrimalRegular.h b/include/Algo/Multiresolution/Map3MR/map3MR_PrimalRegular.h index 2863cdbae13478e28d7a95b4ddcf4bc0a55c43c7..f721c02e31587b648515cfdb84d0ace49c2be6c9 100644 --- a/include/Algo/Multiresolution/Map3MR/map3MR_PrimalRegular.h +++ b/include/Algo/Multiresolution/Map3MR/map3MR_PrimalRegular.h @@ -70,6 +70,8 @@ protected: public: Map3MR(MAP& map); + ~Map3MR(); + private: /*! @name Topological helping functions * diff --git a/include/Algo/Multiresolution/Map3MR/map3MR_PrimalRegular.hpp b/include/Algo/Multiresolution/Map3MR/map3MR_PrimalRegular.hpp index 4b0f39e1d4f3c7b9925783cdd8142e93a5d6d1d2..e9c40572fcc4d29968f7b5abf89670e5718c0041 100644 --- a/include/Algo/Multiresolution/Map3MR/map3MR_PrimalRegular.hpp +++ b/include/Algo/Multiresolution/Map3MR/map3MR_PrimalRegular.hpp @@ -45,6 +45,19 @@ Map3MR::Map3MR(typename PFP::MAP& map) : } +template +Map3MR::~Map3MR() +{ + unsigned int level = m_map.getCurrentLevel(); + unsigned int maxL = m_map.getMaxLevel(); + + for(unsigned int i = maxL ; i > level ; --i) + m_map.removeLevelBack(); + + for(unsigned int i = 0 ; i < level ; ++i) + m_map.removeLevelFront(); +} + /************************************************************************ * Topological helping functions * @@ -138,6 +151,7 @@ void Map3MR::addNewLevelSqrt3(bool embedNewVertices) Algo::Modelisation::Tetrahedralization::flip1To4(m_map, dit); } +/* // // 2-3 swap of all old interior faces // @@ -150,7 +164,7 @@ void Map3MR::addNewLevelSqrt3(bool embedNewVertices) Algo::Modelisation::Tetrahedralization::swap2To3(m_map, dit); } } -/* + // // 1-3 flip of all boundary tetrahedra // @@ -364,37 +378,26 @@ void Map3MR::addNewLevelHexa(bool embedNewVertices) m_map.duplicateDarts(m_map.getMaxLevel()); m_map.setCurrentLevel(m_map.getMaxLevel()); -// if(!shareVertexEmbeddings) -// { -// //create the new level with the old one -// for(unsigned int i = m_mrattribs.begin(); i != m_mrattribs.end(); m_mrattribs.next(i)) -// { -// unsigned int index = (*m_mrDarts[m_mrCurrentLevel])[i] ; -// (*m_embeddings[VERTEX])[index] = EMBNULL ; // set vertex embedding to EMBNULL if no sharing -// } -// } - - //subdivision //1. cut edges TraversorE travE(m_map); for (Dart d = travE.begin(); d != travE.end(); d = travE.next()) { -// if(!shareVertexEmbeddings) -// { -// if(getEmbedding(d) == EMBNULL) -// embedNewCell(d) ; -// if(getEmbedding(phi1(d)) == EMBNULL) -// embedNewCell(phi1(d)) ; -// } + if(!shareVertexEmbeddings && embedNewVertices) + { + if(m_map.template getEmbedding(d) == EMBNULL) + m_map.template embedNewCell(d) ; + if(m_map.template getEmbedding(m_map.phi1(d)) == EMBNULL) + m_map.template embedNewCell(d) ; + } m_map.cutEdge(d) ; travE.skip(d) ; travE.skip(m_map.phi1(d)) ; -// When importing MR files : activated for DEBUG -// if(embedNewVertices) -// embedNewCell(phi1(d)) ; + if(embedNewVertices) + m_map.template embedNewCell(m_map.phi1(d)) ; } + std::cout << "current Level = " << m_map.getCurrentLevel() << std::endl; //2. split faces - quadrangule faces TraversorF travF(m_map) ; @@ -412,9 +415,8 @@ void Map3MR::addNewLevelHexa(bool embedNewVertices) m_map.cutEdge(ne) ; // cut the new edge to insert the central vertex travF.skip(dd) ; -// When importing MR files : activated for DEBUG -// if(embedNewVertices) -// embedNewCell(phi1(ne)) ; + if(embedNewVertices) + m_map.template embedNewCell(m_map.phi1(ne)) ; dd = m_map.phi1(m_map.phi1(next)) ; while(dd != ne) // turn around the face and insert new edges