Commit 7cb42175 authored by Kenneth Vanhoey's avatar Kenneth Vanhoey
Browse files

clamping for cos_angle function

parent 52c137ae
...@@ -59,13 +59,24 @@ void sigmaBilateral(typename PFP::MAP& map, const VertexAttribute<typename PFP:: ...@@ -59,13 +59,24 @@ void sigmaBilateral(typename PFP::MAP& map, const VertexAttribute<typename PFP::
sigmaS = 2.5f * ( sumAngles / float(nbEdges) ) ; sigmaS = 2.5f * ( sumAngles / float(nbEdges) ) ;
} }
/**
* \brief Function applying a bilateral filter smoothing on the mesh.
* \param map the map of the mesh
* \param positionIn the current positions container of the mesh
* \param positionOut the smoothed positions after the function call
* \param normal the normals
*/
template <typename PFP> template <typename PFP>
void filterBilateral(typename PFP::MAP& map, const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& position, VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& position2, const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& normal) void filterBilateral(
typename PFP::MAP& map,
const VertexAttribute<typename PFP::VEC3,typename PFP::MAP>& positionIn,
VertexAttribute<typename PFP::VEC3,typename PFP::MAP>& positionOut,
const VertexAttribute<typename PFP::VEC3,typename PFP::MAP>& normal)
{ {
typedef typename PFP::VEC3 VEC3 ; typedef typename PFP::VEC3 VEC3 ;
float sigmaC, sigmaS ; float sigmaC, sigmaS ;
sigmaBilateral<PFP>(map, position, normal, sigmaC, sigmaS) ; sigmaBilateral<PFP>(map, positionIn, normal, sigmaC, sigmaS) ;
TraversorV<typename PFP::MAP> t(map) ; TraversorV<typename PFP::MAP> t(map) ;
for(Dart d = t.begin(); d != t.end(); d = t.next()) for(Dart d = t.begin(); d != t.end(); d = t.next())
...@@ -80,7 +91,7 @@ void filterBilateral(typename PFP::MAP& map, const VertexAttribute<typename PFP: ...@@ -80,7 +91,7 @@ void filterBilateral(typename PFP::MAP& map, const VertexAttribute<typename PFP:
Traversor2VE<typename PFP::MAP> te(map, d) ; Traversor2VE<typename PFP::MAP> te(map, d) ;
for(Dart it = te.begin(); it != te.end(); it = te.next()) for(Dart it = te.begin(); it != te.end(); it = te.next())
{ {
VEC3 vec = Algo::Surface::Geometry::vectorOutOfDart<PFP>(map, it, position) ; VEC3 vec = Algo::Surface::Geometry::vectorOutOfDart<PFP>(map, it, positionIn) ;
float h = normal_d * vec ; float h = normal_d * vec ;
float t = vec.norm() ; float t = vec.norm() ;
float wcs = exp( ( -1.0f * (t * t) / (2.0f * sigmaC * sigmaC) ) + ( -1.0f * (h * h) / (2.0f * sigmaS * sigmaS) ) ) ; float wcs = exp( ( -1.0f * (t * t) / (2.0f * sigmaC * sigmaC) ) + ( -1.0f * (h * h) / (2.0f * sigmaS * sigmaS) ) ) ;
...@@ -88,10 +99,10 @@ void filterBilateral(typename PFP::MAP& map, const VertexAttribute<typename PFP: ...@@ -88,10 +99,10 @@ void filterBilateral(typename PFP::MAP& map, const VertexAttribute<typename PFP:
normalizer += wcs ; normalizer += wcs ;
} }
position2[d] = position[d] + ((sum / normalizer) * normal_d) ; positionOut[d] = positionIn[d] + ((sum / normalizer) * normal_d) ;
} }
else else
position2[d] = position[d] ; positionOut[d] = positionIn[d] ;
} }
} }
......
...@@ -45,177 +45,177 @@ namespace Geometry ...@@ -45,177 +45,177 @@ namespace Geometry
template<typename PFP, typename V_ATT> template<typename PFP, typename V_ATT>
typename V_ATT::DATA_TYPE triangleNormal(typename PFP::MAP& map, Face f, const V_ATT& position) typename V_ATT::DATA_TYPE triangleNormal(typename PFP::MAP& map, Face f, const V_ATT& position)
{ {
typename V_ATT::DATA_TYPE N = Geom::triangleNormal( typename V_ATT::DATA_TYPE N = Geom::triangleNormal(
position[f.dart], position[f.dart],
position[map.phi1(f)], position[map.phi1(f)],
position[map.phi_1(f)] position[map.phi_1(f)]
) ; ) ;
N.normalize() ; N.normalize() ;
return N ; return N ;
} }
template<typename PFP, typename V_ATT> template<typename PFP, typename V_ATT>
typename V_ATT::DATA_TYPE newellNormal(typename PFP::MAP& map, Face f, const V_ATT& position) typename V_ATT::DATA_TYPE newellNormal(typename PFP::MAP& map, Face f, const V_ATT& position)
{ {
typedef typename V_ATT::DATA_TYPE VEC3; typedef typename V_ATT::DATA_TYPE VEC3;
VEC3 N(0); VEC3 N(0);
foreach_incident2<VERTEX>(map, f, [&] (Vertex v) foreach_incident2<VERTEX>(map, f, [&] (Vertex v)
{ {
const VEC3& P = position[v]; const VEC3& P = position[v];
const VEC3& Q = position[map.phi1(v)]; const VEC3& Q = position[map.phi1(v)];
N[0] += (P[1] - Q[1]) * (P[2] + Q[2]); N[0] += (P[1] - Q[1]) * (P[2] + Q[2]);
N[1] += (P[2] - Q[2]) * (P[0] + Q[0]); N[1] += (P[2] - Q[2]) * (P[0] + Q[0]);
N[2] += (P[0] - Q[0]) * (P[1] + Q[1]); N[2] += (P[0] - Q[0]) * (P[1] + Q[1]);
}); });
N.normalize(); N.normalize();
return N; return N;
} }
template<typename PFP, typename V_ATT> template<typename PFP, typename V_ATT>
typename V_ATT::DATA_TYPE faceNormal(typename PFP::MAP& map, Face f, const V_ATT& position) typename V_ATT::DATA_TYPE faceNormal(typename PFP::MAP& map, Face f, const V_ATT& position)
{ {
if(map.faceDegree(f) == 3) if(map.faceDegree(f) == 3)
return triangleNormal<PFP>(map, f, position) ; return triangleNormal<PFP>(map, f, position) ;
else else
return newellNormal<PFP>(map, f, position) ; return newellNormal<PFP>(map, f, position) ;
} }
template<typename PFP, typename V_ATT> template<typename PFP, typename V_ATT>
typename V_ATT::DATA_TYPE vertexNormal(typename PFP::MAP& map, Vertex v, const V_ATT& position) typename V_ATT::DATA_TYPE vertexNormal(typename PFP::MAP& map, Vertex v, const V_ATT& position)
{ {
typedef typename V_ATT::DATA_TYPE VEC3 ; typedef typename V_ATT::DATA_TYPE VEC3 ;
VEC3 N(0) ; VEC3 N(0) ;
foreach_incident2<FACE>(map, v, [&] (Face f) foreach_incident2<FACE>(map, v, [&] (Face f)
{ {
VEC3 n = faceNormal<PFP>(map, f, position) ; VEC3 n = faceNormal<PFP>(map, f, position) ;
if(!n.hasNan()) if(!n.hasNan())
{ {
VEC3 v1 = vectorOutOfDart<PFP>(map, f.dart, position) ; VEC3 v1 = vectorOutOfDart<PFP>(map, f.dart, position) ;
VEC3 v2 = vectorOutOfDart<PFP>(map, map.phi_1(f), position) ; VEC3 v2 = vectorOutOfDart<PFP>(map, map.phi_1(f), position) ;
n *= convexFaceArea<PFP>(map, f, position) / (v1.norm2() * v2.norm2()) ; n *= convexFaceArea<PFP>(map, f, position) / (v1.norm2() * v2.norm2()) ;
N += n ; N += n ;
} }
}); });
N.normalize() ; N.normalize() ;
return N ; return N ;
} }
template<typename PFP, typename V_ATT> template<typename PFP, typename V_ATT>
typename V_ATT::DATA_TYPE vertexBorderNormal(typename PFP::MAP& map, Vertex v, const V_ATT& position) typename V_ATT::DATA_TYPE vertexBorderNormal(typename PFP::MAP& map, Vertex v, const V_ATT& position)
{ {
assert(map.dimension() == 3); assert(map.dimension() == 3);
typedef typename V_ATT::DATA_TYPE VEC3 ; typedef typename V_ATT::DATA_TYPE VEC3 ;
VEC3 N(0) ; VEC3 N(0) ;
std::vector<Dart> faces; std::vector<Dart> faces;
faces.reserve(16); faces.reserve(16);
map.foreach_dart_of_vertex(v, [&] (Dart d) { faces.push_back(d); }); map.foreach_dart_of_vertex(v, [&] (Dart d) { faces.push_back(d); });
CellMarker<typename PFP::MAP, FACE> f(map); CellMarker<typename PFP::MAP, FACE> f(map);
for(std::vector<Dart>::iterator it = faces.begin() ; it != faces.end() ; ++it) for(std::vector<Dart>::iterator it = faces.begin() ; it != faces.end() ; ++it)
{ {
if(!f.isMarked(*it) && map.isBoundaryIncidentFace(*it)) if(!f.isMarked(*it) && map.isBoundaryIncidentFace(*it))
{ {
f.mark(*it); f.mark(*it);
VEC3 n = faceNormal<PFP>(map, *it, position); VEC3 n = faceNormal<PFP>(map, *it, position);
if(!n.hasNan()) if(!n.hasNan())
{ {
VEC3 v1 = vectorOutOfDart<PFP>(map, *it, position); VEC3 v1 = vectorOutOfDart<PFP>(map, *it, position);
VEC3 v2 = vectorOutOfDart<PFP>(map, map.phi_1(*it), position); VEC3 v2 = vectorOutOfDart<PFP>(map, map.phi_1(*it), position);
n *= convexFaceArea<PFP>(map, *it, position) / (v1.norm2() * v2.norm2()); n *= convexFaceArea<PFP>(map, *it, position) / (v1.norm2() * v2.norm2());
N += n ; N += n ;
} }
} }
} }
N.normalize() ; N.normalize() ;
return N ; return N ;
} }
template <typename PFP, typename V_ATT, typename F_ATT> template <typename PFP, typename V_ATT, typename F_ATT>
void computeNormalFaces(typename PFP::MAP& map, const V_ATT& position, F_ATT& face_normal, unsigned int thread) void computeNormalFaces(typename PFP::MAP& map, const V_ATT& position, F_ATT& face_normal, unsigned int thread)
{ {
if ((CGoGN::Parallel::NumberOfThreads > 1) && (thread == 0)) if ((CGoGN::Parallel::NumberOfThreads > 1) && (thread == 0))
{ {
Parallel::computeNormalFaces<PFP,V_ATT,F_ATT>(map, position, face_normal); Parallel::computeNormalFaces<PFP,V_ATT,F_ATT>(map, position, face_normal);
return; return;
} }
foreach_cell<FACE>(map, [&] (Face f) foreach_cell<FACE>(map, [&] (Face f)
{ {
face_normal[f] = faceNormal<PFP>(map, f, position) ; face_normal[f] = faceNormal<PFP>(map, f, position) ;
}, AUTO, thread); }, AUTO, thread);
} }
template <typename PFP, typename V_ATT> template <typename PFP, typename V_ATT>
void computeNormalVertices(typename PFP::MAP& map, const V_ATT& position, V_ATT& normal, unsigned int thread) void computeNormalVertices(typename PFP::MAP& map, const V_ATT& position, V_ATT& normal, unsigned int thread)
{ {
if ((CGoGN::Parallel::NumberOfThreads > 1) && (thread == 0)) if ((CGoGN::Parallel::NumberOfThreads > 1) && (thread == 0))
{ {
Parallel::computeNormalVertices<PFP,V_ATT>(map, position, normal); Parallel::computeNormalVertices<PFP,V_ATT>(map, position, normal);
return; return;
} }
foreach_cell<VERTEX>(map, [&] (Vertex v) foreach_cell<VERTEX>(map, [&] (Vertex v)
{ {
normal[v] = vertexNormal<PFP>(map, v, position) ; normal[v] = vertexNormal<PFP>(map, v, position) ;
}, FORCE_CELL_MARKING, thread); }, FORCE_CELL_MARKING, thread);
} }
template <typename PFP, typename V_ATT> template <typename PFP, typename V_ATT>
typename PFP::REAL computeAngleBetweenNormalsOnEdge(typename PFP::MAP& map, Edge e, const V_ATT& position) typename PFP::REAL computeAngleBetweenNormalsOnEdge(typename PFP::MAP& map, Edge e, const V_ATT& position)
{ {
typedef typename V_ATT::DATA_TYPE VEC3 ; typedef typename V_ATT::DATA_TYPE VEC3 ;
if(map.isBoundaryEdge(e)) if(map.isBoundaryEdge(e))
return 0 ; return 0 ;
Vertex v1(e.dart); Vertex v1(e.dart);
Vertex v2 = map.phi2(e) ; Vertex v2 = map.phi2(e) ;
const VEC3 n1 = faceNormal<PFP>(map, Face(v1), position) ; const VEC3 n1 = faceNormal<PFP>(map, Face(v1), position) ;
const VEC3 n2 = faceNormal<PFP>(map, Face(v2), position) ; const VEC3 n2 = faceNormal<PFP>(map, Face(v2), position) ;
VEC3 edge = position[v2] - position[v1] ; VEC3 edge = position[v2] - position[v1] ;
edge.normalize() ; edge.normalize() ;
typename PFP::REAL s = edge * (n1 ^ n2) ; typename PFP::REAL s = edge * (n1 ^ n2) ;
typename PFP::REAL c = n1 * n2 ; typename PFP::REAL c = n1 * n2 ;
typename PFP::REAL a(0) ; typename PFP::REAL a(0) ;
// the following trick is useful for avoiding NaNs (due to floating point errors) // the following trick is useful for avoiding NaNs (due to floating point errors)
if (c > 0.5) a = asin(s) ; if (c > 0.5) a = asin(s) ;
else else
{ {
if(c < -1) c = -1 ; if(c < -1) c = -1 ;
if (s >= 0) a = acos(c) ; if (s >= 0) a = acos(c) ;
else a = -acos(c) ; else a = -acos(c) ;
} }
// if (isnan(a)) // if (isnan(a))
if(a != a) if(a != a)
std::cerr<< "Warning : computeAngleBetweenNormalsOnEdge returns NaN on edge " << v1 << "-" << v2 << std::endl ; std::cerr<< "Warning : computeAngleBetweenNormalsOnEdge returns NaN on edge " << v1 << "-" << v2 << std::endl ;
return a ; return a ;
} }
template <typename PFP, typename V_ATT, typename E_ATT> template <typename PFP, typename V_ATT, typename E_ATT>
void computeAnglesBetweenNormalsOnEdges(typename PFP::MAP& map, const V_ATT& position, E_ATT& angles, unsigned int thread) void computeAnglesBetweenNormalsOnEdges(typename PFP::MAP& map, const V_ATT& position, E_ATT& angles, unsigned int thread)
{ {
if ((CGoGN::Parallel::NumberOfThreads > 1) && (thread == 0)) if ((CGoGN::Parallel::NumberOfThreads > 1) && (thread == 0))
{ {
Parallel::computeAnglesBetweenNormalsOnEdges<PFP,V_ATT,E_ATT>(map, position, angles); Parallel::computeAnglesBetweenNormalsOnEdges<PFP,V_ATT,E_ATT>(map, position, angles);
return; return;
} }
foreach_cell<EDGE>(map, [&] (Edge e) foreach_cell<EDGE>(map, [&] (Edge e)
{ {
angles[e] = computeAngleBetweenNormalsOnEdge<PFP>(map, e, position) ; angles[e] = computeAngleBetweenNormalsOnEdge<PFP>(map, e, position) ;
}, AUTO, thread); }, AUTO, thread);
} }
...@@ -225,28 +225,28 @@ namespace Parallel ...@@ -225,28 +225,28 @@ namespace Parallel
template <typename PFP, typename V_ATT> template <typename PFP, typename V_ATT>
void computeNormalVertices(typename PFP::MAP& map, const V_ATT& position, V_ATT& normal) void computeNormalVertices(typename PFP::MAP& map, const V_ATT& position, V_ATT& normal)
{ {
CGoGN::Parallel::foreach_cell<VERTEX>(map, [&] (Vertex v, unsigned int /*thr*/) CGoGN::Parallel::foreach_cell<VERTEX>(map, [&] (Vertex v, unsigned int /*thr*/)
{ {
normal[v] = vertexNormal<PFP>(map, v, position) ; normal[v] = vertexNormal<PFP>(map, v, position) ;
}, FORCE_CELL_MARKING); }, FORCE_CELL_MARKING);
} }
template <typename PFP, typename V_ATT, typename F_ATT> template <typename PFP, typename V_ATT, typename F_ATT>
void computeNormalFaces(typename PFP::MAP& map, const V_ATT& position, F_ATT& normal) void computeNormalFaces(typename PFP::MAP& map, const V_ATT& position, F_ATT& normal)
{ {
CGoGN::Parallel::foreach_cell<FACE>(map, [&] (Face f, unsigned int /*thr*/) CGoGN::Parallel::foreach_cell<FACE>(map, [&] (Face f, unsigned int /*thr*/)
{ {
normal[f] = faceNormal<PFP>(map, f, position) ; normal[f] = faceNormal<PFP>(map, f, position) ;
}); });
} }
template <typename PFP, typename V_ATT, typename E_ATT> template <typename PFP, typename V_ATT, typename E_ATT>
void computeAnglesBetweenNormalsOnEdges(typename PFP::MAP& map, const V_ATT& position, E_ATT& angles) void computeAnglesBetweenNormalsOnEdges(typename PFP::MAP& map, const V_ATT& position, E_ATT& angles)
{ {
CGoGN::Parallel::foreach_cell<EDGE>(map,[&](Edge e, unsigned int /*thr*/) CGoGN::Parallel::foreach_cell<EDGE>(map,[&](Edge e, unsigned int /*thr*/)
{ {
angles[e] = computeAngleBetweenNormalsOnEdge<PFP>(map, e, position) ; angles[e] = computeAngleBetweenNormalsOnEdge<PFP>(map, e, position) ;
}); });
} }
} // namespace Parallel } // namespace Parallel
......
...@@ -76,9 +76,11 @@ Vector<DIM,T> isobarycenter(const Vector<DIM,T>& v1, const Vector<DIM,T>& v2, co ...@@ -76,9 +76,11 @@ Vector<DIM,T> isobarycenter(const Vector<DIM,T>& v1, const Vector<DIM,T>& v2, co
template <typename VEC> template <typename VEC>
typename VEC::DATA_TYPE cos_angle(const VEC& a, const VEC& b) typename VEC::DATA_TYPE cos_angle(const VEC& a, const VEC& b)
{ {
typename VEC::DATA_TYPE na2 = a.norm2() ; typename VEC::DATA_TYPE na2 = a.norm2() ;
typename VEC::DATA_TYPE nb2 = b.norm2() ; typename VEC::DATA_TYPE nb2 = b.norm2() ;
return (a * b) / sqrt(na2 * nb2) ;
typename VEC::DATA_TYPE res = (a * b) / sqrt(na2 * nb2) ;
return res > 1.0 ? 1.0 : (res < -1.0 ? -1.0 : res) ;
} }
// angle formed by 2 vectors // angle formed by 2 vectors
......
Markdown is supported
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