Skip to content
GitLab
Commit d204e178 authored by Pierre Kraemer's avatar Pierre Kraemer
Browse files

amelioration gestion des surfaces a bord

parent 1a289716
Hide whitespace changes
Inline Side-by-side
include/Algo/Geometry/area.hpp
View file @ d204e178
......@@ -116,16 +116,19 @@ typename PFP::REAL vertexVoronoiArea(typename PFP::MAP& map, Dart d, const typen
Dart it = d ;
do
{
if(!isTriangleObtuse<PFP>(map, it, position))
const typename PFP::VEC3& p1 = position[it] ;
const typename PFP::VEC3& p2 = position[map.phi1(it)] ;
const typename PFP::VEC3& p3 = position[map.phi_1(it)] ;
if(!Geom::isTriangleObtuse(p1, p2, p3))
{
typename PFP::REAL a = angle<PFP>(map, map.phi1(it), map.phi2(it), position) ;
typename PFP::REAL b = angle<PFP>(map, map.phi_1(it), map.phi2(map.phi1(it)), position) ;
area += (vectorOutOfDart<PFP>(map, it, position).norm2() / tan(a) + vectorOutOfDart<PFP>(map, map.phi_1(it), position).norm2() / tan(b)) / 8 ;
typename PFP::REAL a = Geom::angle(p3 - p2, p1 - p2) ;
typename PFP::REAL b = Geom::angle(p1 - p3, p2 - p3) ;
area += ( (p2 - p1).norm2() / tan(b) + (p3 - p1).norm2() / tan(a) ) / 8 ;
}
else
{
typename PFP::REAL tArea = convexFaceArea<PFP>(map, it, position) ;
if(angle<PFP>(map, it, map.phi2(map.phi_1(it)), position) > M_PI / 2)
typename PFP::REAL tArea = Geom::triangleArea(p1, p2, p3) ;
if(Geom::angle(p2 - p1, p3 - p1) > M_PI / 2)
area += tArea / 2 ;
else
area += tArea / 4 ;
......
include/Algo/Geometry/laplacian.hpp
View file @ d204e178
......@@ -40,15 +40,22 @@ ATTR_TYPE computeLaplacianTopoVertex(
const AttributeHandler<ATTR_TYPE>& attr)
{
ATTR_TYPE l(0) ;
unsigned int val = 0 ;
ATTR_TYPE value = attr[d] ;
unsigned int wSum = 0 ;
Dart it = d ;
do
{
l += attr[map.phi1(it)] - attr[it] ;
val++ ;
l += attr[map.phi1(it)] - value ;
++wSum ;
Dart dboundary = map.phi_1(it) ;
if(map.phi2(dboundary) == dboundary)
{
l += attr[dboundary] - value ;
++wSum ;
}
it = map.alpha1(it) ;
} while(it != d) ;
l /= val ;
l /= wSum ;
return l ;
}
......@@ -64,15 +71,23 @@ ATTR_TYPE computeLaplacianCotanVertex(
ATTR_TYPE l(0) ;
Dart it = d ;
REAL vArea = vertexArea[d] ;
REAL val = 0 ;
ATTR_TYPE value = attr[d] ;
REAL wSum = 0 ;
do
{
REAL w = edgeWeight[it] / vArea ;
l += (attr[map.phi1(it)] - attr[it]) * w ;
val += w ;
l += (attr[map.phi1(it)] - value) * w ;
wSum += w ;
Dart dboundary = map.phi_1(it) ;
if(map.phi2(dboundary) == dboundary)
{
w = edgeWeight[dboundary] / vArea ;
l += (attr[dboundary] - value) * w ;
wSum += w ;
}
it = map.alpha1(it) ;
} while(it != d) ;
l /= val ;
l /= wSum ;
return l ;
}
......@@ -120,11 +135,17 @@ typename PFP::REAL computeCotanWeightEdge(
Dart d,
const typename PFP::TVEC3& position)
{
typename PFP::REAL cot_alpha = 1 / tan(angle<PFP>(map, map.phi_1(d), map.phi2(map.phi1(d)), position)) ;
const typename PFP::VEC3& p1 = position[d] ;
const typename PFP::VEC3& p2 = position[map.phi1(d)] ;
const typename PFP::VEC3& p3 = position[map.phi_1(d)] ;
typename PFP::REAL cot_alpha = 1 / tan(Geom::angle(p1 - p3, p2 - p3)) ;
typename PFP::REAL cot_beta = 0 ;
Dart dd = map.phi2(d) ;
if(dd != d)
cot_beta = 1 / tan(angle<PFP>(map, map.phi_1(dd), map.phi2(map.phi1(dd)), position)) ;
{
const typename PFP::VEC3& p4 = position[map.phi_1(map.phi2(d))] ;
cot_beta = 1 / tan(Geom::angle(p2 - p4, p1 - p4)) ;
}
return 0.5 * ( cot_alpha + cot_beta ) ;
}
......
include/Algo/Geometry/normal.hpp
View file @ d204e178
......@@ -168,6 +168,9 @@ typename PFP::REAL computeAngleBetweenNormalsOnEdge(typename PFP::MAP& map, Dart
typedef typename PFP::VEC3 VEC3 ;
Dart dd = map.phi2(d) ;
if(dd == d)
return 0 ;
const VEC3 n1 = faceNormal<PFP>(map, d, position) ;
const VEC3 n2 = faceNormal<PFP>(map, dd, position) ;
VEC3 e = position[dd] - position[d] ;
......
include/Algo/LinearSolving/matrixSetup.h
View file @ d204e178
......@@ -131,6 +131,12 @@ public:
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) ;
......@@ -175,6 +181,12 @@ public:
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) ;
......@@ -221,6 +233,12 @@ public:
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) ;
......@@ -268,6 +286,13 @@ public:
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) ;
......@@ -317,6 +342,13 @@ public:
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) ;
......@@ -368,6 +400,13 @@ public:
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) ;
......
include/Geometry/basic.h
View file @ d204e178
......@@ -106,10 +106,10 @@ VEC3 triangleNormal(const VEC3& p1, const VEC3& p2, const VEC3& p3)
template <typename VEC3>
bool isTriangleObtuse(const VEC3& p1, const VEC3& p2, const VEC3& p3)
{
typename VEC3::DATA_TYPE a1 = angle(p2-p1, p3-p1) ;
typename VEC3::DATA_TYPE a1 = angle(p2 - p1, p3 - p1) ;
if(a1 > M_PI / 2)
return true ;
typename VEC3::DATA_TYPE a2 = angle(p3-p2, p1-p2) ;
typename VEC3::DATA_TYPE a2 = angle(p3 - p2, p1 - p2) ;
if(a2 > M_PI / 2 || a1 + a2 < M_PI / 2)
return true ;
return false ;
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment