Commit bd0de59b authored by Pierre Kraemer's avatar Pierre Kraemer

merge no_boundary + GMap without fixed point

parent 7054966d
......@@ -57,16 +57,7 @@ struct PFP: public PFP_STANDARD
typedef Map2TP MAP;
};
/// definition de la carte en global, plus facile
PFP::MAP myMap;
// handler d'attribut de position par sommet
AttributeHandler<PFP::VEC3> position;
// handler d'attribut de normale par sommet
AttributeHandler<PFP::VEC3> normal;
/// encore 1 typedef pour simplifier l'ecriture du code
/// typedef pour simplifier l'ecriture du code
typedef PFP::VEC3 Point3D;
// Variables pour la gestion des plongements
......
......@@ -51,30 +51,35 @@ void filterAverageAttribute_OneRing(
TraversorV<typename PFP::MAP> t(map, select) ;
for(Dart d = t.begin(); d != t.end(); d = t.next())
{
if (neigh & INSIDE)
col.collectAll(d) ;
else
col.collectBorder(d) ;
fa.reset() ;
if (neigh & INSIDE)
if(!map.isBoundaryVertex(d))
{
switch (attIn.getOrbit())
if (neigh & INSIDE)
col.collectAll(d) ;
else
col.collectBorder(d) ;
fa.reset() ;
if (neigh & INSIDE)
{
case VERTEX :
col.applyOnInsideVertices(fa) ;
break;
case EDGE :
col.applyOnInsideEdges(fa) ;
break;
case FACE :
col.applyOnInsideFaces(fa) ;
break;
switch (attIn.getOrbit())
{
case VERTEX :
col.applyOnInsideVertices(fa) ;
break;
case EDGE :
col.applyOnInsideEdges(fa) ;
break;
case FACE :
col.applyOnInsideFaces(fa) ;
break;
}
}
if (neigh & BORDER)
col.applyOnBorder(fa) ;
attOut[d] = fa.getAverage() ;
}
if (neigh & BORDER)
col.applyOnBorder(fa) ;
attOut[d] = fa.getAverage() ;
else
attOut[d] = attIn[d] ;
}
}
......@@ -95,23 +100,28 @@ void filterAverageVertexAttribute_WithinSphere(
TraversorV<typename PFP::MAP> t(map, select) ;
for(Dart d = t.begin(); d != t.end(); d = t.next())
{
if (neigh & INSIDE)
col.collectAll(d) ;
else
col.collectBorder(d) ;
attOut[d] = T(0);
if (neigh & INSIDE){
faInside.reset() ;
col.applyOnInsideVertices(faInside) ;
attOut[d] += faInside.getSum();
}
if (neigh & BORDER){
faBorder.reset(position[d], radius);
col.applyOnBorder(faBorder) ;
attOut[d] += faBorder.getSum();
if(!map.isBoundaryVertex(d))
{
if (neigh & INSIDE)
col.collectAll(d) ;
else
col.collectBorder(d) ;
attOut[d] = T(0);
if (neigh & INSIDE){
faInside.reset() ;
col.applyOnInsideVertices(faInside) ;
attOut[d] += faInside.getSum();
}
if (neigh & BORDER){
faBorder.reset(position[d], radius);
col.applyOnBorder(faBorder) ;
attOut[d] += faBorder.getSum();
}
attOut[d] /= faInside.getCount() + faBorder.getCount() ;
}
attOut[d] /= faInside.getCount() + faBorder.getCount() ;
else
attOut[d] = attIn[d] ;
}
}
......@@ -128,7 +138,7 @@ void filterAverageEdgeAttribute_WithinSphere(
FunctorAverage<T> fa(attIn) ;
Algo::Selection::Collector_WithinSphere<PFP> col(map, position, radius) ;
TraversorV<typename PFP::MAP> t(map, select) ;
TraversorE<typename PFP::MAP> t(map, select) ;
for(Dart d = t.begin(); d != t.end(); d = t.next())
{
if (neigh & INSIDE)
......@@ -156,7 +166,7 @@ void filterAverageFaceAttribute_WithinSphere(
FunctorAverage<T> fa(attIn) ;
Algo::Selection::Collector_WithinSphere<PFP> col(map, position, radius) ;
TraversorV<typename PFP::MAP> t(map, select) ;
TraversorF<typename PFP::MAP> t(map, select) ;
for(Dart d = t.begin(); d != t.end(); d = t.next())
{
if (neigh & INSIDE)
......
......@@ -344,8 +344,6 @@ void filterVNBA(typename PFP::MAP& map, float sigmaN2, float SUSANthreshold, con
long nbAdapt = 0 ;
long nbSusan = 0 ;
std::cout << "compute new vertices normals.." << std::endl ;
TraversorV<typename PFP::MAP> tv(map, select) ;
for(Dart d = tv.begin(); d != tv.end(); d = tv.next())
{
......@@ -434,8 +432,6 @@ void filterVNBA(typename PFP::MAP& map, float sigmaN2, float SUSANthreshold, con
}
}
std::cout << "update face normals.." << std::endl ;
// Compute face normals from vertex normals
TraversorF<typename PFP::MAP> tf(map, select) ;
for(Dart d = tf.begin(); d != tf.end(); d = tf.next())
......@@ -455,8 +451,6 @@ void filterVNBA(typename PFP::MAP& map, float sigmaN2, float SUSANthreshold, con
faceNewNormal[d] = newNormal ;
}
std::cout << "update vertices positions.." << std::endl ;
// Compute new vertices position
computeNewPositionsFromFaceNormals<PFP>(
map, position, position2, faceArea, faceCentroid, faceNormal, faceNewNormal, select
......
......@@ -67,24 +67,25 @@ void filterBilateral(typename PFP::MAP& map, const typename PFP::TVEC3& position
TraversorV<typename PFP::MAP> t(map, select) ;
for(Dart d = t.begin(); d != t.end(); d = t.next())
{
// get position & normal of vertex
const VEC3& pos_d = position[d] ;
const VEC3& normal_d = normal[d] ;
// traversal of incident edges
float sum = 0.0f, normalizer = 0.0f ;
Traversor2VE<typename PFP::MAP> te(map, d) ;
for(Dart it = te.begin(); it != te.end(); it = te.next())
if(!map.isBoundaryVertex(d))
{
VEC3 vec = Algo::Geometry::vectorOutOfDart<PFP>(map, it, position) ;
float h = normal_d * vec ;
float t = vec.norm() ;
float wcs = exp( ( -1.0f * (t * t) / (2.0f * sigmaC * sigmaC) ) + ( -1.0f * (h * h) / (2.0f * sigmaS * sigmaS) ) ) ;
sum += wcs * h ;
normalizer += wcs ;
}
// traversal of incident edges
float sum = 0.0f, normalizer = 0.0f ;
Traversor2VE<typename PFP::MAP> te(map, d) ;
for(Dart it = te.begin(); it != te.end(); it = te.next())
{
VEC3 vec = Algo::Geometry::vectorOutOfDart<PFP>(map, it, position) ;
float h = normal_d * vec ;
float t = vec.norm() ;
float wcs = exp( ( -1.0f * (t * t) / (2.0f * sigmaC * sigmaC) ) + ( -1.0f * (h * h) / (2.0f * sigmaS * sigmaS) ) ) ;
sum += wcs * h ;
normalizer += wcs ;
}
position2[d] = pos_d + ((sum / normalizer) * normal_d) ;
position2[d] = position[d] + ((sum / normalizer) * normal[d]) ;
}
else
position2[d] = position[d] ;
}
}
......@@ -102,39 +103,44 @@ void filterSUSAN(typename PFP::MAP& map, float SUSANthreshold, const typename PF
TraversorV<typename PFP::MAP> t(map, select) ;
for(Dart d = t.begin(); d != t.end(); d = t.next())
{
// get position & normal of vertex
const VEC3& pos_d = position[d] ;
const VEC3& normal_d = normal[d] ;
// traversal of incident edges
float sum = 0.0f, normalizer = 0.0f ;
bool SUSANregion = false ;
Traversor2VE<typename PFP::MAP> te(map, d) ;
for(Dart it = te.begin(); it != te.end(); it = te.next())
if(!map.isBoundaryVertex(d))
{
const VEC3& neighborNormal = normal[map.phi1(it)] ;
float angle = Geom::angle(normal_d, neighborNormal) ;
if( angle <= SUSANthreshold )
// get position & normal of vertex
const VEC3& pos_d = position[d] ;
const VEC3& normal_d = normal[d] ;
// traversal of incident edges
float sum = 0.0f, normalizer = 0.0f ;
bool SUSANregion = false ;
Traversor2VE<typename PFP::MAP> te(map, d) ;
for(Dart it = te.begin(); it != te.end(); it = te.next())
{
VEC3 vec = Algo::Geometry::vectorOutOfDart<PFP>(map, it, position) ;
float h = normal_d * vec ;
float t = vec.norm() ;
float wcs = exp( ( -1.0f * (t * t) / (2.0f * sigmaC * sigmaC) ) + ( -1.0f * (h * h) / (2.0f * sigmaS * sigmaS) ) );
sum += wcs * h ;
normalizer += wcs ;
const VEC3& neighborNormal = normal[map.phi1(it)] ;
float angle = Geom::angle(normal_d, neighborNormal) ;
if( angle <= SUSANthreshold )
{
VEC3 vec = Algo::Geometry::vectorOutOfDart<PFP>(map, it, position) ;
float h = normal_d * vec ;
float t = vec.norm() ;
float wcs = exp( ( -1.0f * (t * t) / (2.0f * sigmaC * sigmaC) ) + ( -1.0f * (h * h) / (2.0f * sigmaS * sigmaS) ) );
sum += wcs * h ;
normalizer += wcs ;
}
else
SUSANregion = true ;
}
else
SUSANregion = true ;
}
if(SUSANregion)
nbSusan++ ;
nbTot++ ;
if(SUSANregion)
nbSusan++ ;
nbTot++ ;
if (normalizer != 0.0f)
position2[d] = pos_d + ((sum / normalizer) * normal_d) ;
if (normalizer != 0.0f)
position2[d] = pos_d + ((sum / normalizer) * normal_d) ;
else
position2[d] = pos_d ;
}
else
position2[d] = pos_d ;
position2[d] = position[d] ;
}
// CGoGNout <<" susan rate = "<< float(nbSusan)/float(nbTot)<<CGoGNendl;
......
......@@ -53,13 +53,18 @@ void filterTaubin(typename PFP::MAP& map, typename PFP::TVEC3& position, typenam
{
mv.mark(d);
c.collectBorder(d) ;
fa1.reset() ;
c.applyOnBorder(fa1) ;
VEC3 p = position[d] ;
VEC3 displ = fa1.getAverage() - p ;
displ *= lambda ;
position2[d] = p + displ ;
if(!map.isBoundaryVertex(d))
{
c.collectBorder(d) ;
fa1.reset() ;
c.applyOnBorder(fa1) ;
VEC3 p = position[d] ;
VEC3 displ = fa1.getAverage() - p ;
displ *= lambda ;
position2[d] = p + displ ;
}
else
position2[d] = position[d] ;
}
}
......@@ -71,13 +76,18 @@ void filterTaubin(typename PFP::MAP& map, typename PFP::TVEC3& position, typenam
{
mv.unmark(d);
c.collectBorder(d) ;
fa2.reset() ;
c.applyOnBorder(fa2) ;
VEC3 p = position2[d] ;
VEC3 displ = fa2.getAverage() - p ;
displ *= mu ;
position[d] = p + displ ;
if(!map.isBoundaryVertex(d))
{
c.collectBorder(d) ;
fa2.reset() ;
c.applyOnBorder(fa2) ;
VEC3 p = position2[d] ;
VEC3 displ = fa2.getAverage() - p ;
displ *= mu ;
position[d] = p + displ ;
}
else
position[d] = position2[d] ;
}
}
}
......@@ -93,22 +103,28 @@ void filterTaubin_modified(typename PFP::MAP& map, typename PFP::TVEC3& position
const float lambda = 0.6307 ;
const float mu = -0.6732 ;
CellMarkerNoUnmark mv(map, VERTEX) ;
FunctorAverageOnSphereBorder<PFP, typename PFP::VEC3> fa1(map, position, position) ;
Algo::Selection::Collector_WithinSphere<PFP> c1(map, position, radius) ;
CellMarkerNoUnmark mv(map, VERTEX) ;
for(Dart d = map.begin(); d != map.end(); map.next(d))
{
if(select(d) && !mv.isMarked(d))
{
mv.mark(d);
c1.collectBorder(d) ;
VEC3 center = position[d] ;
fa1.reset(center, radius) ;
c1.applyOnBorder(fa1) ;
VEC3 displ = fa1.getAverage() - center ;
displ *= lambda ;
position2[d] = center + displ ;
if(!map.isBoundaryVertex(d))
{
c1.collectBorder(d) ;
VEC3 center = position[d] ;
fa1.reset(center, radius) ;
c1.applyOnBorder(fa1) ;
VEC3 displ = fa1.getAverage() - center ;
displ *= lambda ;
position2[d] = center + displ ;
}
else
position2[d] = position[d] ;
}
}
......@@ -121,13 +137,18 @@ void filterTaubin_modified(typename PFP::MAP& map, typename PFP::TVEC3& position
{
mv.unmark(d);
c2.collectBorder(d) ;
VEC3 center = position2[d] ;
fa2.reset(center, radius) ;
c2.applyOnBorder(fa2) ;
VEC3 displ = fa2.getAverage() - center ;
displ *= mu ;
position[d] = center + displ ;
if(!map.isBoundaryVertex(d))
{
c2.collectBorder(d) ;
VEC3 center = position2[d] ;
fa2.reset(center, radius) ;
c2.applyOnBorder(fa2) ;
VEC3 displ = fa2.getAverage() - center ;
displ *= mu ;
position[d] = center + displ ;
}
else
position[d] = position2[d] ;
}
}
}
......
......@@ -81,11 +81,11 @@ public:
Dart phi2(Dart d) ;
Dart alpha0(Dart d);
Dart alpha0(Dart d) ;
Dart alpha1(Dart d);
Dart alpha1(Dart d) ;
Dart alpha_1(Dart d);
Dart alpha_1(Dart d) ;
virtual Dart begin() ;
......@@ -97,10 +97,10 @@ public:
virtual bool foreach_dart_of_edge(Dart d, FunctorType& f, unsigned int thread = 0) ;
bool foreach_dart_of_oriented_face(Dart d, FunctorType& f, unsigned int thread = 0);
virtual bool foreach_dart_of_oriented_face(Dart d, FunctorType& f, unsigned int thread = 0) ;
virtual bool foreach_dart_of_face(Dart d, FunctorType& f, unsigned int thread = 0) ;
bool foreach_dart_of_oriented_volume(Dart d, FunctorType& f, unsigned int thread = 0);
virtual bool foreach_dart_of_oriented_volume(Dart d, FunctorType& f, unsigned int thread = 0) ;
virtual bool foreach_dart_of_volume(Dart d, FunctorType& f, unsigned int thread = 0) ;
virtual bool foreach_dart_of_cc(Dart d, FunctorType& f, unsigned int thread = 0) ;
......
......@@ -142,7 +142,7 @@ inline Dart ImplicitHierarchicalMap::alpha_1(Dart d)
inline Dart ImplicitHierarchicalMap::begin()
{
Dart d = Map2::begin() ;
while(m_dartLevel[d] > m_curLevel)
while(d != Map2::end() && m_dartLevel[d] > m_curLevel)
Map2::next(d) ;
return d ;
}
......
......@@ -184,21 +184,14 @@ public:
bool operator()(Dart d)
{
solver->begin_row() ;
Dart it = d ;
REAL aii = 0 ;
do
Traversor2VE<typename PFP::MAP> t(this->m_map, d) ;
for(Dart it = t.begin(); it != t.end(); it = t.next())
{
REAL aij = 1 ;
aii += aij ;
solver->add_coefficient(indexTable[this->m_map.phi1(it)], aij) ;
Dart dboundary = this->m_map.phi_1(it) ;
if(this->m_map.phi2(dboundary) == dboundary)
{
aii += aij ;
solver->add_coefficient(indexTable[dboundary], aij) ;
}
it = this->m_map.alpha1(it) ;
} while(it != d) ;
}
solver->add_coefficient(indexTable[d], -aii) ;
solver->normalize_row() ;
solver->set_right_hand_side(0) ;
......@@ -234,21 +227,14 @@ public:
bool operator()(Dart d)
{
solver->begin_row() ;
Dart it = d ;
REAL aii = 0 ;
do
Traversor2VE<typename PFP::MAP> t(this->m_map, d) ;
for(Dart it = t.begin(); it != t.end(); it = t.next())
{
REAL aij = 1 ;
aii += aij ;
solver->add_coefficient(indexTable[this->m_map.phi1(it)], aij) ;
Dart dboundary = this->m_map.phi_1(it) ;
if(this->m_map.phi2(dboundary) == dboundary)
{
aii += aij ;
solver->add_coefficient(indexTable[dboundary], aij) ;
}
it = this->m_map.alpha1(it) ;
} while(it != d) ;
}
solver->add_coefficient(indexTable[d], -aii) ;
solver->normalize_row() ;
solver->set_right_hand_side(attrTable[d]) ;
......@@ -286,21 +272,14 @@ public:
bool operator()(Dart d)
{
solver->begin_row() ;
Dart it = d ;
REAL aii = 0 ;
do
Traversor2VE<typename PFP::MAP> t(this->m_map, d) ;
for(Dart it = t.begin(); it != t.end(); it = t.next())
{
REAL aij = 1 ;
aii += aij ;
solver->add_coefficient(indexTable[this->m_map.phi1(it)], aij) ;
Dart dboundary = this->m_map.phi_1(it) ;
if(this->m_map.phi2(dboundary) == dboundary)
{
aii += aij ;
solver->add_coefficient(indexTable[dboundary], aij) ;
}
it = this->m_map.alpha1(it) ;
} while(it != d) ;
}
solver->add_coefficient(indexTable[d], -aii) ;
solver->normalize_row() ;
solver->set_right_hand_side((attrTable[d])[coord]) ;
......@@ -338,21 +317,16 @@ public:
bool operator()(Dart d)
{
solver->begin_row() ;
Dart it = d ;
REAL vArea = vertexArea[d] ;
REAL aii = 0 ;
Dart it = d ;
// Traversor2VE<typename PFP::MAP> t(this->m_map, d) ;
// for(Dart it = t.begin(); it != t.end(); it = t.next())
do
{
REAL aij = edgeWeight[it] / vArea ;
aii += aij ;
solver->add_coefficient(indexTable[this->m_map.phi1(it)], aij) ;
Dart dboundary = this->m_map.phi_1(it) ;
if(this->m_map.phi2(dboundary) == dboundary)
{
aij = edgeWeight[dboundary] / vArea ;
aii += aij ;
solver->add_coefficient(indexTable[dboundary], aij) ;
}
it = this->m_map.alpha1(it) ;
} while(it != d) ;
solver->add_coefficient(indexTable[d], -aii) ;
......@@ -394,23 +368,15 @@ public:
bool operator()(Dart d)
{
solver->begin_row() ;
Dart it = d ;
REAL vArea = vertexArea[d] ;
REAL aii = 0 ;
do
Traversor2VE<typename PFP::MAP> t(this->m_map, d) ;
for(Dart it = t.begin(); it != t.end(); it = t.next())
{
REAL aij = edgeWeight[it] / vArea ;
aii += aij ;
solver->add_coefficient(indexTable[this->m_map.phi1(it)], aij) ;
Dart dboundary = this->m_map.phi_1(it) ;
if(this->m_map.phi2(dboundary) == dboundary)
{
aij = edgeWeight[dboundary] / vArea ;
aii += aij ;
solver->add_coefficient(indexTable[dboundary], aij) ;
}
it = this->m_map.alpha1(it) ;
} while(it != d) ;
}
solver->add_coefficient(indexTable[d], -aii) ;
solver->normalize_row() ;
solver->set_right_hand_side(attrTable[d]) ;
......@@ -452,23 +418,15 @@ public:
bool operator()(Dart d)
{
solver->begin_row() ;
Dart it = d ;
REAL vArea = vertexArea[d] ;
REAL aii = 0 ;
do
Traversor2VE<typename PFP::MAP> t(this->m_map, d) ;
for(Dart it = t.begin(); it != t.end(); it = t.next())
{
REAL aij = edgeWeight[it] / vArea ;
aii += aij ;
solver->add_coefficient(indexTable[this->m_map.phi1(it)], aij) ;
Dart dboundary = this->m_map.phi_1(it) ;
if(this->m_map.phi2(dboundary) == dboundary)
{
aij = edgeWeight[dboundary] / vArea ;
aii += aij ;
solver->add_coefficient(indexTable[dboundary], aij) ;
}
it = this->m_map.alpha1(it) ;
} while(it != d) ;
}
solver->add_coefficient(indexTable[d], -aii) ;
solver->normalize_row() ;
solver->set_right_hand_side((attrTable[d])[coord]) ;
......
......@@ -299,18 +299,8 @@ Dart extrudeFace(typename PFP::MAP& the_map, typename PFP::TVEC3& positions, Dar
cc = the_map.alpha1(cc);
}while (cc != c);
//merge central faces by removing edges
bool notFinished=true;
do
{
Dart d1 = the_map.alpha1(cc);
if (d1 == cc) // last edge is pending edge inside of face
notFinished = false;
the_map.deleteFace(cc);
cc = d1;
} while (notFinished);
the_map.closeHole(the_map.phi1(the_map.phi1(d)));
// delete the central vertex
the_map.deleteVertex(c) ;
// embedding of new vertices
Dart dd = d;
......