Skip to content

GitLab

C
CGoGN
Commit 4e5e69b1 authored by Kenneth Vanhoey's avatar Kenneth Vanhoey
Browse files
Options

Merge cgogn:~cgogn/CGoGN

parents 4268ecf1 947fe6a0
Hide whitespace changes
Inline Side-by-side
Showing with 906 additions and 140 deletions
Apps/Examples/viewer.cpp
View file @ 4e5e69b1
......@@ -185,13 +185,9 @@ void Viewer::importMesh(std::string& filename)
m_render->initPrimitives<PFP>(myMap, allDarts, Algo::Render::GL2::LINES) ;
m_render->initPrimitives<PFP>(myMap, allDarts, Algo::Render::GL2::TRIANGLES) ;
Geom::BoundingBox<PFP::VEC3> bb = Algo::Geometry::computeBoundingBox<PFP>(myMap, position) ;
bb = Algo::Geometry::computeBoundingBox<PFP>(myMap, position) ;
gPosObj = bb.center() ;
float tailleX = bb.size(0) ;
float tailleY = bb.size(1) ;
float tailleZ = bb.size(2) ;
gWidthObj = std::max<float>(std::max<float>(tailleX, tailleY), tailleZ) ;
normalBaseSize = std::min<float>(std::min<float>(tailleX,tailleY),tailleZ) / 50.0f ;
normalBaseSize = bb.diagSize() / 100.0f ;
vertexBaseSize = normalBaseSize * 2.0f ;
if(!normal.isValid())
......@@ -202,7 +198,7 @@ void Viewer::importMesh(std::string& filename)
m_positionVBO->updateData(position) ;
m_normalVBO->updateData(normal) ;
setParamObject(gWidthObj, gPosObj.data()) ;
setParamObject(bb.maxSize(), gPosObj.data()) ;
updateGLMatrices() ;
}
......
Apps/Examples/viewer.h
View file @ 4e5e69b1
......@@ -77,8 +77,8 @@ public:
float shininess ;
Geom::BoundingBox<PFP::VEC3> bb ;
Geom::Vec3f gPosObj ;
float gWidthObj ;
float normalBaseSize ;
float normalScaleFactor ;
float vertexBaseSize ;
......
include/Algo/Geometry/curvature.h
View file @ 4e5e69b1
......@@ -43,30 +43,29 @@ namespace Geometry
typedef CPULinearSolverTraits< SparseMatrix<double>, FullVector<double> > CPUSolverTraits ;
template <typename PFP>
void computeCurvatureVertices(
void computeCurvatureVertices_QuadraticFitting(
typename PFP::MAP& map,
const typename PFP::TVEC3& position,
const typename PFP::TVEC3& normal,
typename PFP::TREAL& k1,
typename PFP::TREAL& k2,
typename PFP::TVEC3& K1,
typename PFP::TVEC3& K2,
typename PFP::TREAL& kmax,
typename PFP::TREAL& kmin,
typename PFP::TVEC3& Kmax,
typename PFP::TVEC3& Kmin,
const FunctorSelect& select = SelectorTrue()) ;
template <typename PFP>
void computeCurvatureVertex(
void computeCurvatureVertex_QuadraticFitting(
typename PFP::MAP& map,
Dart dart,
const typename PFP::TVEC3& position,
const typename PFP::TVEC3& normal,
typename PFP::TREAL& k1,
typename PFP::TREAL& k2,
typename PFP::TVEC3& K1,
typename PFP::TVEC3& K2) ;
typename PFP::TREAL& kmax,
typename PFP::TREAL& kmin,
typename PFP::TVEC3& Kmax,
typename PFP::TVEC3& Kmin) ;
template <typename PFP>
void vertexQuadricFitting(
void vertexQuadraticFitting(
typename PFP::MAP& map,
Dart dart,
typename PFP::MATRIX33& localFrame,
......@@ -75,10 +74,10 @@ void vertexQuadricFitting(
float& a, float& b, float& c, float& d, float& e) ;
template <typename PFP>
void quadricFittingAddVertexPos(typename PFP::VEC3& v, typename PFP::VEC3& p, typename PFP::MATRIX33& localFrame, LinearSolver<CPUSolverTraits>& solver) ;
void quadraticFittingAddVertexPos(typename PFP::VEC3& v, typename PFP::VEC3& p, typename PFP::MATRIX33& localFrame, LinearSolver<CPUSolverTraits>& solver) ;
template <typename PFP>
void quadricFittingAddVertexNormal(typename PFP::VEC3& v, typename PFP::VEC3& n, typename PFP::VEC3& p, typename PFP::MATRIX33& localFrame, LinearSolver<CPUSolverTraits>& solver) ;
void quadraticFittingAddVertexNormal(typename PFP::VEC3& v, typename PFP::VEC3& n, typename PFP::VEC3& p, typename PFP::MATRIX33& localFrame, LinearSolver<CPUSolverTraits>& solver) ;
/*
template <typename PFP>
void vertexCubicFitting(Dart dart, typename PFP::VEC3& normal, float& a, float& b, float& c, float& d, float& e, float& f, float& g, float& h, float& i) ;
......@@ -90,7 +89,35 @@ template <typename PFP>
void cubicFittingAddVertexNormal(typename PFP::VEC3& v, typename PFP::VEC3& n, typename PFP::VEC3& p, typename PFP::MATRIX33& localFrame) ;
*/
} // namespace Geoemtry
template <typename PFP>
void computeCurvatureVertices_NormalCycles(
typename PFP::MAP& map,
typename PFP::REAL radius,
const typename PFP::TVEC3& position,
const typename PFP::TVEC3& normal,
const typename PFP::TREAL& angles,
typename PFP::TREAL& kmax,
typename PFP::TREAL& kmin,
typename PFP::TVEC3& Kmax,
typename PFP::TVEC3& Kmin,
typename PFP::TVEC3& Knormal,
const FunctorSelect& select = SelectorTrue()) ;
template <typename PFP>
void computeCurvatureVertex_NormalCycles(
typename PFP::MAP& map,
Dart dart,
typename PFP::REAL radius,
const typename PFP::TVEC3& position,
const typename PFP::TVEC3& normal,
const typename PFP::TREAL& angles,
typename PFP::TREAL& kmax,
typename PFP::TREAL& kmin,
typename PFP::TVEC3& Kmax,
typename PFP::TVEC3& Kmin,
typename PFP::TVEC3& Knormal) ;
} // namespace Geometry
} // namespace Algo
......
include/Algo/Geometry/curvature.hpp
View file @ 4e5e69b1
......@@ -25,6 +25,7 @@
#include "Algo/Geometry/localFrame.h"
#include "Geometry/matrix.h"
#include "Topology/generic/cellmarker.h"
#include "Algo/Selection/collector.h"
extern "C"
{
......@@ -44,15 +45,15 @@ namespace Geometry
{
template <typename PFP>
void computeCurvatureVertices(
typename PFP::MAP& map,
const typename PFP::TVEC3& position,
const typename PFP::TVEC3& normal,
typename PFP::TREAL& k1,
typename PFP::TREAL& k2,
typename PFP::TVEC3& K1,
typename PFP::TVEC3& K2,
const FunctorSelect& select)
void computeCurvatureVertices_QuadraticFitting(
typename PFP::MAP& map,
const typename PFP::TVEC3& position,
const typename PFP::TVEC3& normal,
typename PFP::TREAL& kmax,
typename PFP::TREAL& kmin,
typename PFP::TVEC3& Kmax,
typename PFP::TVEC3& Kmin,
const FunctorSelect& select)
{
CellMarker marker(map, VERTEX);
for(Dart d = map.begin(); d != map.end(); map.next(d))
......@@ -60,61 +61,62 @@ void computeCurvatureVertices(
if(select(d) && !marker.isMarked(d))
{
marker.mark(d);
computeCurvatureVertex<PFP>(map, d, position, normal, k1, k2, K1, K2) ;
computeCurvatureVertex_QuadraticFitting<PFP>(map, d, position, normal, kmax, kmin, Kmax, Kmin) ;
}
}
}
template <typename PFP>
void computeCurvatureVertex(
typename PFP::MAP& map,
Dart dart,
const typename PFP::TVEC3& position,
const typename PFP::TVEC3& normal,
typename PFP::TREAL& k1,
typename PFP::TREAL& k2,
typename PFP::TVEC3& K1,
typename PFP::TVEC3& K2)
void computeCurvatureVertex_QuadraticFitting(
typename PFP::MAP& map,
Dart dart,
const typename PFP::TVEC3& position,
const typename PFP::TVEC3& normal,
typename PFP::TREAL& kmax,
typename PFP::TREAL& kmin,
typename PFP::TVEC3& Kmax,
typename PFP::TVEC3& Kmin)
{
typedef typename PFP::REAL REAL ;
typedef typename PFP::VEC3 VEC3 ;
typedef typename PFP::MATRIX33 MATRIX33 ;
VEC3 n = normal[dart] ;
MATRIX33 localFrame = Algo::Geometry::vertexLocalFrame<PFP>(map, dart, position, n) ;
MATRIX33 invLocalFrame ;
localFrame.invert(invLocalFrame) ;
REAL a, b, c, d, e;
//vertexCubicFitting(map,dart,localFrame,a,b,c,d,e,f,g,h,i) ;
vertexQuadricFitting<PFP>(map, dart, localFrame, position, normal, a, b, c, d, e) ;
vertexQuadraticFitting<PFP>(map, dart, localFrame, position, normal, a, b, c, d, e) ;
REAL k1_v, k2_v, K1x, K1y ;
REAL kmax_v, kmin_v, Kmax_x, Kmax_y ;
//int res = slaev2_(&e, &f, &g, &maxC, &minC, &dirX, &dirY) ;
/*int res = */slaev2_(&a, &b, &c, &k1_v, &k2_v, &K1x, &K1y) ;
/*int res = */slaev2_(&a, &b, &c, &kmax_v, &kmin_v, &Kmax_x, &Kmax_y) ;
VEC3 K1_v(K1x, K1y, 0.0f) ;
K1_v = invLocalFrame * K1_v ;
VEC3 K2_v = n ^ K1_v ;
VEC3 Kmax_v(Kmax_x, Kmax_y, 0.0f) ;
Kmax_v = invLocalFrame * Kmax_v ;
VEC3 Kmin_v = n ^ Kmax_v ;
if (k1_v < k2_v)
if (kmax_v < kmin_v)
{
k1[dart] = -k1_v ;
k2[dart] = -k2_v ;
K1[dart] = K1_v ;
K2[dart] = K2_v ;
kmax[dart] = -kmax_v ;
kmin[dart] = -kmin_v ;
Kmax[dart] = Kmax_v ;
Kmin[dart] = Kmin_v ;
}
else
{
k1[dart] = -k2_v ;
k2[dart] = -k1_v ;
K1[dart] = K2_v ;
K2[dart] = K1_v ;
kmax[dart] = -kmin_v ;
kmin[dart] = -kmax_v ;
Kmax[dart] = Kmin_v ;
Kmin[dart] = Kmax_v ;
}
}
template <typename PFP>
void vertexQuadricFitting(
void vertexQuadraticFitting(
typename PFP::MAP& map,
Dart dart,
typename PFP::MATRIX33& localFrame,
......@@ -122,11 +124,7 @@ void vertexQuadricFitting(
const typename PFP::TVEC3& normal,
float& a, float& b, float& c, float& d, float& e)
{
typedef typename PFP::REAL REAL ;
typedef typename PFP::VEC3 VEC3 ;
typedef typename PFP::MATRIX33 MATRIX33 ;
VEC3 p = position[dart] ;
typename PFP::VEC3 p = position[dart] ;
LinearSolver<CPUSolverTraits> solver(5) ;
solver.set_least_squares(true) ;
......@@ -135,10 +133,10 @@ void vertexQuadricFitting(
do
{
// 1-ring vertices
VEC3 v = position[map.phi2(it)] ;
quadricFittingAddVertexPos<PFP>(v, p, localFrame, solver) ;
VEC3 n = normal[map.phi2(it)] ;
quadricFittingAddVertexNormal<PFP>(v, n, p, localFrame, solver) ;
typename PFP::VEC3 v = position[map.phi2(it)] ;
quadraticFittingAddVertexPos<PFP>(v, p, localFrame, solver) ;
typename PFP::VEC3 n = normal[map.phi2(it)] ;
quadraticFittingAddVertexNormal<PFP>(v, n, p, localFrame, solver) ;
// 2-ring vertices
// Dart d2 = map.phi1(map.phi1(map.phi2(map.phi1(it)))) ;
// VEC3 v2 = position[d2] ;
......@@ -158,7 +156,7 @@ void vertexQuadricFitting(
}
template <typename PFP>
void quadricFittingAddVertexPos(typename PFP::VEC3& v, typename PFP::VEC3& p, typename PFP::MATRIX33& localFrame, LinearSolver<CPUSolverTraits>& solver)
void quadraticFittingAddVertexPos(typename PFP::VEC3& v, typename PFP::VEC3& p, typename PFP::MATRIX33& localFrame, LinearSolver<CPUSolverTraits>& solver)
{
typename PFP::VEC3 vec = v - p ;
vec = localFrame * vec ;
......@@ -175,33 +173,33 @@ void quadricFittingAddVertexPos(typename PFP::VEC3& v, typename PFP::VEC3& p, ty
}
template <typename PFP>
void quadricFittingAddVertexNormal(typename PFP::VEC3& v, typename PFP::VEC3& n, typename PFP::VEC3& p, typename PFP::MATRIX33& localFrame, LinearSolver<CPUSolverTraits>& solver)
void quadraticFittingAddVertexNormal(typename PFP::VEC3& v, typename PFP::VEC3& n, typename PFP::VEC3& p, typename PFP::MATRIX33& localFrame, LinearSolver<CPUSolverTraits>& solver)
{
typename PFP::VEC3 vec = v - p ;
vec = localFrame * vec ;
typename PFP::VEC3 norm = localFrame * n ;
solver.begin_row() ;
solver.add_coefficient(0, 2.0f*vec[0]) ;
solver.add_coefficient(0, 2.0f * vec[0]) ;
solver.add_coefficient(1, vec[1]) ;
solver.add_coefficient(2, 0) ;
solver.add_coefficient(3, 1.0f) ;
solver.add_coefficient(4, 0) ;
solver.set_right_hand_side(-1.0f*norm[0]/norm[2]) ;
solver.set_right_hand_side(-1.0f * norm[0] / norm[2]) ;
solver.end_row() ;
solver.begin_row() ;
solver.add_coefficient(0, 0) ;
solver.add_coefficient(1, vec[0]) ;
solver.add_coefficient(2, 2.0f*vec[1]) ;
solver.add_coefficient(2, 2.0f * vec[1]) ;
solver.add_coefficient(3, 0) ;
solver.add_coefficient(4, 1.0f) ;
solver.set_right_hand_side(-1.0f*norm[1]/norm[2]) ;
solver.set_right_hand_side(-1.0f * norm[1] / norm[2]) ;
solver.end_row() ;
}
/*
template <typename PFP>
void DifferentialProperties<PFP>::vertexCubicFitting(Dart dart, gmtl::Vec3f& normal, float& a, float& b, float& c, float& d, float& e, float& f, float& g, float& h, float& i)
void vertexCubicFitting(Dart dart, gmtl::Vec3f& normal, float& a, float& b, float& c, float& d, float& e, float& f, float& g, float& h, float& i)
{
gmtl::Matrix33f localFrame, invLocalFrame ;
Algo::Geometry::vertexLocalFrame<PFP>(m_map,dart,normal,localFrame) ;
......@@ -245,7 +243,7 @@ void DifferentialProperties<PFP>::vertexCubicFitting(Dart dart, gmtl::Vec3f& nor
}
template <typename PFP>
void DifferentialProperties<PFP>::cubicFittingAddVertexPos(gmtl::Vec3f& v, gmtl::Vec3f& p, gmtl::Matrix33f& localFrame)
void cubicFittingAddVertexPos(gmtl::Vec3f& v, gmtl::Vec3f& p, gmtl::Matrix33f& localFrame)
{
gmtl::Vec3f vec = v - p ;
vec = localFrame * vec ;
......@@ -266,7 +264,7 @@ void DifferentialProperties<PFP>::cubicFittingAddVertexPos(gmtl::Vec3f& v, gmtl:
}
template <typename PFP>
void DifferentialProperties<PFP>::cubicFittingAddVertexNormal(gmtl::Vec3f& v, gmtl::Vec3f& n, gmtl::Vec3f& p, gmtl::Matrix33f& localFrame)
void cubicFittingAddVertexNormal(gmtl::Vec3f& v, gmtl::Vec3f& n, gmtl::Vec3f& p, gmtl::Matrix33f& localFrame)
{
gmtl::Vec3f vec = v - p ;
vec = localFrame * vec ;
......@@ -299,6 +297,116 @@ void DifferentialProperties<PFP>::cubicFittingAddVertexNormal(gmtl::Vec3f& v, gm
solverC->end_row() ;
}
*/
template <typename PFP>
void computeCurvatureVertices_NormalCycles(
typename PFP::MAP& map,
typename PFP::REAL radius,
const typename PFP::TVEC3& position,
const typename PFP::TVEC3& normal,
const typename PFP::TREAL& angles,
typename PFP::TREAL& kmax,
typename PFP::TREAL& kmin,
typename PFP::TVEC3& Kmax,
typename PFP::TVEC3& Kmin,
typename PFP::TVEC3& Knormal,
const FunctorSelect& select = SelectorTrue())
{
CellMarker marker(map, VERTEX);
for(Dart d = map.begin(); d != map.end(); map.next(d))
{
if(select(d) && !marker.isMarked(d))
{
marker.mark(d);
computeCurvatureVertex_NormalCycles<PFP>(map, d, radius, position, normal, angles, kmax, kmin, Kmax, Kmin, Knormal) ;
}
}
}
template <typename PFP>
void computeCurvatureVertex_NormalCycles(
typename PFP::MAP& map,
Dart dart,
typename PFP::REAL radius,
const typename PFP::TVEC3& position,
const typename PFP::TVEC3& normal,
const typename PFP::TREAL& angles,
typename PFP::TREAL& kmax,
typename PFP::TREAL& kmin,
typename PFP::TVEC3& Kmax,
typename PFP::TVEC3& Kmin,
typename PFP::TVEC3& Knormal)
{
typedef typename PFP::VEC3 VEC3 ;
typedef typename PFP::REAL REAL ;
Algo::Selection::Collector_WithinSphere<PFP> neigh(map, position) ;
neigh.init(dart, radius) ;
neigh.collect() ;
neigh.computeArea() ;
VEC3 center = position[neigh.getCenter()] ;
typename PFP::MATRIX33 tensor(0) ;
// inside
const std::vector<Dart>& vd1 = neigh.getInsideEdges() ;
for (std::vector<Dart>::const_iterator it = vd1.begin(); it != vd1.end(); ++it)
{
const VEC3 e = position[map.phi2(*it)] - position[*it] ;
tensor += Geom::transposed_vectors_mult(e,e) * angles[*it] * (1 / e.norm()) ;
}
// border
const std::vector<Dart>& vd2 = neigh.getBorder() ;
for (std::vector<Dart>::const_iterator it = vd2.begin(); it != vd2.end(); ++it)
{
const VEC3 e = position[map.phi2(*it)] - position[*it] ;
const REAL alpha = neigh.intersect_SphereEdge(*it, map.phi2(*it)) ;
tensor += Geom::transposed_vectors_mult(e,e) * angles[*it] * (1 / e.norm()) * alpha ;
}
tensor /= neigh.getArea() ;
long int n = 3, lda = 3, info, lwork = 9 ;
char jobz='V', uplo = 'U' ;
float work[9] ;
float w[3] ;
float a[3*3] = {
tensor(0,0), 0.0f, 0.0f,
tensor(1,0), tensor(1,1), 0.0f,
tensor(2,0), tensor(2,1), tensor(2,2)
} ;
// Solve eigenproblem
ssyev_(&jobz, &uplo, (integer*)&n, a, (integer*)&lda, w, work, (integer*)&lwork, (integer*)&info) ;
// Check for convergence
if(info > 0)
std::cerr << "clapack ssyev_ failed to compute eigenvalues : exit with code " << info << std::endl ;
// sort eigen components : w[s[0]] has minimal absolute value ; kmin = w[s[1]] <= w[s[2]] = kmax
int s[3] = {0, 1, 2} ;
int tmp ;
if (abs(w[s[2]]) < abs(w[s[1]])) { tmp = s[1] ; s[1] = s[2] ; s[2] = tmp ; }
if (abs(w[s[1]]) < abs(w[s[0]])) { tmp = s[0] ; s[0] = s[1] ; s[1] = tmp ; }
if (w[s[2]] < w[s[1]]) { tmp = s[1] ; s[1] = s[2] ; s[2] = tmp ; }
kmin[dart] = w[s[1]] ;
kmax[dart] = w[s[2]] ;
VEC3& dirMin = Kmin[dart] ;
dirMin[0] = a[3*s[2]];
dirMin[1] = a[3*s[2]+1];
dirMin[2] = a[3*s[2]+2]; // warning : Kmin and Kmax are switched
VEC3& dirMax = Kmax[dart] ;
dirMax[0] = a[3*s[1]];
dirMax[1] = a[3*s[1]+1];
dirMax[2] = a[3*s[1]+2]; // warning : Kmin and Kmax are switched
VEC3& dirNormal = Knormal[dart] ;
dirNormal[0] = a[3*s[0]];
dirNormal[1] = a[3*s[0]+1];
dirNormal[2] = a[3*s[0]+2];
if (dirNormal * normal[dart] < 0)
dirNormal *= -1; // change orientation
}
} // namespace Geometry
} // namespace Algo
......
include/Algo/Geometry/normal.h
View file @ 4e5e69b1
......@@ -37,16 +37,16 @@ namespace Geometry
{
template <typename PFP>
typename PFP::VEC3 triangleNormal(typename PFP::MAP& map, Dart d, const typename PFP::TVEC3& position);
typename PFP::VEC3 triangleNormal(typename PFP::MAP& map, Dart d, const typename PFP::TVEC3& position) ;
template <typename PFP>
typename PFP::VEC3 faceNormal(typename PFP::MAP& map, Dart d, const typename PFP::TVEC3& position);
typename PFP::VEC3 faceNormal(typename PFP::MAP& map, Dart d, const typename PFP::TVEC3& position) ;
template <typename PFP>
typename PFP::VEC3 vertexNormal(typename PFP::MAP& map, Dart d, const typename PFP::TVEC3& position);
typename PFP::VEC3 vertexNormal(typename PFP::MAP& map, Dart d, const typename PFP::TVEC3& position) ;
template <typename PFP>
void computeNormalFaces(typename PFP::MAP& map, const typename PFP::TVEC3& position, typename PFP::TVEC3& face_normal, const FunctorSelect& select = SelectorTrue(), unsigned int thread=0) ;
void computeNormalFaces(typename PFP::MAP& map, const typename PFP::TVEC3& position, typename PFP::TVEC3& face_normal, const FunctorSelect& select = SelectorTrue(), unsigned int thread = 0) ;
/**
* compute normals of vertices
......@@ -57,7 +57,13 @@ void computeNormalFaces(typename PFP::MAP& map, const typename PFP::TVEC3& posit
* @ param th the thread number
*/
template <typename PFP>
void computeNormalVertices(typename PFP::MAP& map, const typename PFP::TVEC3& position, typename PFP::TVEC3& normal, const FunctorSelect& select = SelectorTrue(), unsigned int thread=0) ;
void computeNormalVertices(typename PFP::MAP& map, const typename PFP::TVEC3& position, typename PFP::TVEC3& normal, const FunctorSelect& select = SelectorTrue(), unsigned int thread = 0) ;
template <typename PFP>
typename PFP::REAL computeAngleBetweenNormalsOnEdge(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& position) ;
template <typename PFP>
void computeAnglesBetweenNormalsOnEdges(typename PFP::MAP& map, typename PFP::TVEC3& position, typename PFP::TREAL& angles, const FunctorSelect& select = SelectorTrue(), unsigned int thread = 0) ;
} // namespace Geometry
......
include/Algo/Geometry/normal.hpp
View file @ 4e5e69b1
......@@ -112,7 +112,7 @@ void computeNormalFaces(typename PFP::MAP& map, const typename PFP::TVEC3& posit
}
template <typename PFP>
class computeNormalVerticesFunctor: public FunctorMap<typename PFP::MAP>
class computeNormalVerticesFunctor : public FunctorMap<typename PFP::MAP>
{
protected:
typename PFP::MAP& m_map;
......@@ -120,8 +120,8 @@ protected:
typename PFP::TVEC3& m_normal;
public:
computeNormalVerticesFunctor(typename PFP::MAP& map, const typename PFP::TVEC3& position, typename PFP::TVEC3& normal):
m_map(map), m_position(position), m_normal(normal) {}
m_map(map), m_position(position), m_normal(normal)
{}
bool operator()(Dart d)
{
m_normal[d] = vertexNormal<PFP>(m_map, d, m_position) ;
......@@ -129,12 +129,10 @@ public:
}
};
template <typename PFP>
void computeNormalVertices(typename PFP::MAP& map, const typename PFP::TVEC3& position, typename PFP::TVEC3& normal, const FunctorSelect& select, unsigned int thread)
{
CellMarker marker(map, VERTEX,thread);
CellMarker marker(map, VERTEX, thread);
for(Dart d = map.begin(); d != map.end(); map.next(d))
{
if(select(d) && !marker.isMarked(d))
......@@ -145,6 +143,46 @@ void computeNormalVertices(typename PFP::MAP& map, const typename PFP::TVEC3& po
}
}
template <typename PFP>
typename PFP::REAL computeAngleBetweenNormalsOnEdge(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& position)
{
typedef typename PFP::VEC3 VEC3 ;
Dart dd = map.phi2(d) ;
const VEC3 n1 = faceNormal<PFP>(map, d, position) ;
const VEC3 n2 = faceNormal<PFP>(map, dd, position) ;
VEC3 e = position[dd] - position[d] ;
e.normalize() ;
typename PFP::REAL s = e * (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))
std::cerr<< "Warning : computeAngleBetweenNormalsOnEdge returns NaN on edge " << d << "-" << dd << std::endl ;
return a ;
}
template <typename PFP>
void computeAnglesBetweenNormalsOnEdges(typename PFP::MAP& map, typename PFP::TVEC3& position, typename PFP::TREAL& angles, const FunctorSelect& select, unsigned int thread)
{
CellMarker me(map, EDGE, thread) ;
for(Dart d = map.begin(); d != map.end(); map.next(d))
{
if(select(d) && !me.isMarked(d))
{
me.mark(d) ;
angles[d] = computeAngleBetweenNormalsOnEdge<PFP>(map, d, position) ;
}
}
}
} // namespace Geometry
} // namespace Algo
......
include/Algo/Selection/collector.h 0 → 100644
View file @ 4e5e69b1
/*******************************************************************************
* CGoGN: Combinatorial and Geometric modeling with Generic N-dimensional Maps *