From 7cb4217512a146f19a34cf78d673809baf61c2b8 Mon Sep 17 00:00:00 2001 From: Kenneth Vanhoey Date: Thu, 10 Jul 2014 14:44:39 +0200 Subject: [PATCH] clamping for cos_angle function --- include/Algo/Filtering/bilateral.h | 21 ++- include/Algo/Geometry/normal.hpp | 286 ++++++++++++++--------------- include/Geometry/basic.h | 8 +- 3 files changed, 164 insertions(+), 151 deletions(-) diff --git a/include/Algo/Filtering/bilateral.h b/include/Algo/Filtering/bilateral.h index 9492eac4..35209cc2 100644 --- a/include/Algo/Filtering/bilateral.h +++ b/include/Algo/Filtering/bilateral.h @@ -59,13 +59,24 @@ void sigmaBilateral(typename PFP::MAP& map, const VertexAttribute -void filterBilateral(typename PFP::MAP& map, const VertexAttribute& position, VertexAttribute& position2, const VertexAttribute& normal) +void filterBilateral( + typename PFP::MAP& map, + const VertexAttribute& positionIn, + VertexAttribute& positionOut, + const VertexAttribute& normal) { typedef typename PFP::VEC3 VEC3 ; float sigmaC, sigmaS ; - sigmaBilateral(map, position, normal, sigmaC, sigmaS) ; + sigmaBilateral(map, positionIn, normal, sigmaC, sigmaS) ; TraversorV t(map) ; for(Dart d = t.begin(); d != t.end(); d = t.next()) @@ -80,7 +91,7 @@ void filterBilateral(typename PFP::MAP& map, const VertexAttribute te(map, d) ; for(Dart it = te.begin(); it != te.end(); it = te.next()) { - VEC3 vec = Algo::Surface::Geometry::vectorOutOfDart(map, it, position) ; + VEC3 vec = Algo::Surface::Geometry::vectorOutOfDart(map, it, positionIn) ; 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) ) ) ; @@ -88,10 +99,10 @@ void filterBilateral(typename PFP::MAP& map, const VertexAttribute typename V_ATT::DATA_TYPE triangleNormal(typename PFP::MAP& map, Face f, const V_ATT& position) { - typename V_ATT::DATA_TYPE N = Geom::triangleNormal( - position[f.dart], - position[map.phi1(f)], - position[map.phi_1(f)] - ) ; - N.normalize() ; - return N ; + typename V_ATT::DATA_TYPE N = Geom::triangleNormal( + position[f.dart], + position[map.phi1(f)], + position[map.phi_1(f)] + ) ; + N.normalize() ; + return N ; } template typename V_ATT::DATA_TYPE newellNormal(typename PFP::MAP& map, Face f, const V_ATT& position) { - typedef typename V_ATT::DATA_TYPE VEC3; - VEC3 N(0); - - foreach_incident2(map, f, [&] (Vertex v) - { - const VEC3& P = position[v]; - const VEC3& Q = position[map.phi1(v)]; - N[0] += (P[1] - Q[1]) * (P[2] + Q[2]); - N[1] += (P[2] - Q[2]) * (P[0] + Q[0]); - N[2] += (P[0] - Q[0]) * (P[1] + Q[1]); - }); - - N.normalize(); - return N; + typedef typename V_ATT::DATA_TYPE VEC3; + VEC3 N(0); + + foreach_incident2(map, f, [&] (Vertex v) + { + const VEC3& P = position[v]; + const VEC3& Q = position[map.phi1(v)]; + N[0] += (P[1] - Q[1]) * (P[2] + Q[2]); + N[1] += (P[2] - Q[2]) * (P[0] + Q[0]); + N[2] += (P[0] - Q[0]) * (P[1] + Q[1]); + }); + + N.normalize(); + return N; } template typename V_ATT::DATA_TYPE faceNormal(typename PFP::MAP& map, Face f, const V_ATT& position) { - if(map.faceDegree(f) == 3) - return triangleNormal(map, f, position) ; - else - return newellNormal(map, f, position) ; + if(map.faceDegree(f) == 3) + return triangleNormal(map, f, position) ; + else + return newellNormal(map, f, position) ; } template typename V_ATT::DATA_TYPE vertexNormal(typename PFP::MAP& map, Vertex v, const V_ATT& position) { - typedef typename V_ATT::DATA_TYPE VEC3 ; - - VEC3 N(0) ; - - foreach_incident2(map, v, [&] (Face f) - { - VEC3 n = faceNormal(map, f, position) ; - if(!n.hasNan()) - { - VEC3 v1 = vectorOutOfDart(map, f.dart, position) ; - VEC3 v2 = vectorOutOfDart(map, map.phi_1(f), position) ; - n *= convexFaceArea(map, f, position) / (v1.norm2() * v2.norm2()) ; - N += n ; - } - }); - - N.normalize() ; - return N ; + typedef typename V_ATT::DATA_TYPE VEC3 ; + + VEC3 N(0) ; + + foreach_incident2(map, v, [&] (Face f) + { + VEC3 n = faceNormal(map, f, position) ; + if(!n.hasNan()) + { + VEC3 v1 = vectorOutOfDart(map, f.dart, position) ; + VEC3 v2 = vectorOutOfDart(map, map.phi_1(f), position) ; + n *= convexFaceArea(map, f, position) / (v1.norm2() * v2.norm2()) ; + N += n ; + } + }); + + N.normalize() ; + return N ; } template typename V_ATT::DATA_TYPE vertexBorderNormal(typename PFP::MAP& map, Vertex v, const V_ATT& position) { - assert(map.dimension() == 3); - - typedef typename V_ATT::DATA_TYPE VEC3 ; - - VEC3 N(0) ; - - std::vector faces; - faces.reserve(16); - map.foreach_dart_of_vertex(v, [&] (Dart d) { faces.push_back(d); }); - - CellMarker f(map); - - for(std::vector::iterator it = faces.begin() ; it != faces.end() ; ++it) - { - if(!f.isMarked(*it) && map.isBoundaryIncidentFace(*it)) - { - f.mark(*it); - VEC3 n = faceNormal(map, *it, position); - if(!n.hasNan()) - { - VEC3 v1 = vectorOutOfDart(map, *it, position); - VEC3 v2 = vectorOutOfDart(map, map.phi_1(*it), position); - n *= convexFaceArea(map, *it, position) / (v1.norm2() * v2.norm2()); - N += n ; - } - } - } - - N.normalize() ; - return N ; + assert(map.dimension() == 3); + + typedef typename V_ATT::DATA_TYPE VEC3 ; + + VEC3 N(0) ; + + std::vector faces; + faces.reserve(16); + map.foreach_dart_of_vertex(v, [&] (Dart d) { faces.push_back(d); }); + + CellMarker f(map); + + for(std::vector::iterator it = faces.begin() ; it != faces.end() ; ++it) + { + if(!f.isMarked(*it) && map.isBoundaryIncidentFace(*it)) + { + f.mark(*it); + VEC3 n = faceNormal(map, *it, position); + if(!n.hasNan()) + { + VEC3 v1 = vectorOutOfDart(map, *it, position); + VEC3 v2 = vectorOutOfDart(map, map.phi_1(*it), position); + n *= convexFaceArea(map, *it, position) / (v1.norm2() * v2.norm2()); + N += n ; + } + } + } + + N.normalize() ; + return N ; } template void computeNormalFaces(typename PFP::MAP& map, const V_ATT& position, F_ATT& face_normal, unsigned int thread) { - if ((CGoGN::Parallel::NumberOfThreads > 1) && (thread == 0)) - { - Parallel::computeNormalFaces(map, position, face_normal); - return; - } - - foreach_cell(map, [&] (Face f) - { - face_normal[f] = faceNormal(map, f, position) ; - }, AUTO, thread); + if ((CGoGN::Parallel::NumberOfThreads > 1) && (thread == 0)) + { + Parallel::computeNormalFaces(map, position, face_normal); + return; + } + + foreach_cell(map, [&] (Face f) + { + face_normal[f] = faceNormal(map, f, position) ; + }, AUTO, thread); } template void computeNormalVertices(typename PFP::MAP& map, const V_ATT& position, V_ATT& normal, unsigned int thread) { - if ((CGoGN::Parallel::NumberOfThreads > 1) && (thread == 0)) - { - Parallel::computeNormalVertices(map, position, normal); - return; - } - - foreach_cell(map, [&] (Vertex v) - { - normal[v] = vertexNormal(map, v, position) ; - }, FORCE_CELL_MARKING, thread); + if ((CGoGN::Parallel::NumberOfThreads > 1) && (thread == 0)) + { + Parallel::computeNormalVertices(map, position, normal); + return; + } + + foreach_cell(map, [&] (Vertex v) + { + normal[v] = vertexNormal(map, v, position) ; + }, FORCE_CELL_MARKING, thread); } template typename PFP::REAL computeAngleBetweenNormalsOnEdge(typename PFP::MAP& map, Edge e, const V_ATT& position) { - typedef typename V_ATT::DATA_TYPE VEC3 ; - - if(map.isBoundaryEdge(e)) - return 0 ; - - Vertex v1(e.dart); - Vertex v2 = map.phi2(e) ; - const VEC3 n1 = faceNormal(map, Face(v1), position) ; - const VEC3 n2 = faceNormal(map, Face(v2), position) ; - VEC3 edge = position[v2] - position[v1] ; - edge.normalize() ; - typename PFP::REAL s = edge * (n1 ^ n2) ; - typename PFP::REAL c = n1 * n2 ; - typename PFP::REAL a(0) ; - - // the following trick is useful for avoiding NaNs (due to floating point errors) - if (c > 0.5) a = asin(s) ; - else - { - if(c < -1) c = -1 ; - if (s >= 0) a = acos(c) ; - else a = -acos(c) ; - } - // if (isnan(a)) - if(a != a) - std::cerr<< "Warning : computeAngleBetweenNormalsOnEdge returns NaN on edge " << v1 << "-" << v2 << std::endl ; - - return a ; + typedef typename V_ATT::DATA_TYPE VEC3 ; + + if(map.isBoundaryEdge(e)) + return 0 ; + + Vertex v1(e.dart); + Vertex v2 = map.phi2(e) ; + const VEC3 n1 = faceNormal(map, Face(v1), position) ; + const VEC3 n2 = faceNormal(map, Face(v2), position) ; + VEC3 edge = position[v2] - position[v1] ; + edge.normalize() ; + typename PFP::REAL s = edge * (n1 ^ n2) ; + typename PFP::REAL c = n1 * n2 ; + typename PFP::REAL a(0) ; + + // the following trick is useful for avoiding NaNs (due to floating point errors) + if (c > 0.5) a = asin(s) ; + else + { + if(c < -1) c = -1 ; + if (s >= 0) a = acos(c) ; + else a = -acos(c) ; + } + // if (isnan(a)) + if(a != a) + std::cerr<< "Warning : computeAngleBetweenNormalsOnEdge returns NaN on edge " << v1 << "-" << v2 << std::endl ; + + return a ; } template void computeAnglesBetweenNormalsOnEdges(typename PFP::MAP& map, const V_ATT& position, E_ATT& angles, unsigned int thread) { - if ((CGoGN::Parallel::NumberOfThreads > 1) && (thread == 0)) - { - Parallel::computeAnglesBetweenNormalsOnEdges(map, position, angles); - return; - } - - foreach_cell(map, [&] (Edge e) - { - angles[e] = computeAngleBetweenNormalsOnEdge(map, e, position) ; - }, AUTO, thread); + if ((CGoGN::Parallel::NumberOfThreads > 1) && (thread == 0)) + { + Parallel::computeAnglesBetweenNormalsOnEdges(map, position, angles); + return; + } + + foreach_cell(map, [&] (Edge e) + { + angles[e] = computeAngleBetweenNormalsOnEdge(map, e, position) ; + }, AUTO, thread); } @@ -225,28 +225,28 @@ namespace Parallel template void computeNormalVertices(typename PFP::MAP& map, const V_ATT& position, V_ATT& normal) { - CGoGN::Parallel::foreach_cell(map, [&] (Vertex v, unsigned int /*thr*/) - { - normal[v] = vertexNormal(map, v, position) ; - }, FORCE_CELL_MARKING); + CGoGN::Parallel::foreach_cell(map, [&] (Vertex v, unsigned int /*thr*/) + { + normal[v] = vertexNormal(map, v, position) ; + }, FORCE_CELL_MARKING); } template void computeNormalFaces(typename PFP::MAP& map, const V_ATT& position, F_ATT& normal) { - CGoGN::Parallel::foreach_cell(map, [&] (Face f, unsigned int /*thr*/) - { - normal[f] = faceNormal(map, f, position) ; - }); + CGoGN::Parallel::foreach_cell(map, [&] (Face f, unsigned int /*thr*/) + { + normal[f] = faceNormal(map, f, position) ; + }); } template void computeAnglesBetweenNormalsOnEdges(typename PFP::MAP& map, const V_ATT& position, E_ATT& angles) { - CGoGN::Parallel::foreach_cell(map,[&](Edge e, unsigned int /*thr*/) - { - angles[e] = computeAngleBetweenNormalsOnEdge(map, e, position) ; - }); + CGoGN::Parallel::foreach_cell(map,[&](Edge e, unsigned int /*thr*/) + { + angles[e] = computeAngleBetweenNormalsOnEdge(map, e, position) ; + }); } } // namespace Parallel diff --git a/include/Geometry/basic.h b/include/Geometry/basic.h index e830d1f9..70ddd30c 100644 --- a/include/Geometry/basic.h +++ b/include/Geometry/basic.h @@ -76,9 +76,11 @@ Vector isobarycenter(const Vector& v1, const Vector& v2, co template typename VEC::DATA_TYPE cos_angle(const VEC& a, const VEC& b) { - typename VEC::DATA_TYPE na2 = a.norm2() ; - typename VEC::DATA_TYPE nb2 = b.norm2() ; - return (a * b) / sqrt(na2 * nb2) ; + typename VEC::DATA_TYPE na2 = a.norm2() ; + typename VEC::DATA_TYPE nb2 = b.norm2() ; + + 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 -- 2.24.1