ajout collector + debut integration normal cycles

include/Algo/Geometry/curvature.h

@@ -116,7 +116,7 @@ void computeCurvatureVertex_NormalCycles(
 	typename PFP::TVEC3& Kmin,
 	typename PFP::TVEC3& Knormal) ;
 
-} // namespace Geoemtry
+} // namespace Geometry
 
 } // namespace Algo
 

include/Algo/Geometry/curvature.hpp

@@ -337,10 +337,74 @@ void computeCurvatureVertex_NormalCycles(
 	typename PFP::TVEC3& Kmin,
 	typename PFP::TVEC3& Knormal)
 {
-	Collector_WithinSphere<PFP> neigh(map, position);
+	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

include/Algo/Geometry/normal.h

@@ -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

@@ -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,50 @@ 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

@@ -35,6 +35,9 @@
 namespace CGoGN
 {
 
+namespace Algo
+{
+
 namespace Selection
 {
 
@@ -154,6 +157,8 @@ public:
 
 } // namespace Selection
 
+} // namespace Algo
+
 } // namespace CGoGN
 
 #include "Algo/Selection/collector.hpp"

include/Algo/Selection/collector.hpp

@@ -25,6 +25,9 @@
 namespace CGoGN
 {
 
+namespace Algo
+{
+
 namespace Selection
 {
 
@@ -201,4 +204,6 @@ void Collector_WithinSphere<PFP>::computeArea()
 
 } // namespace Selection
 
+} // namespace Algo
+
 } // namespace CGoGN