add vertexVoronoiArea + isTriangleObtuse

include/Algo/Geometry/area.h

@@ -46,6 +46,9 @@ typename PFP::REAL totalArea(typename PFP::MAP& map, const typename PFP::TVEC3&
 template <typename PFP>
 typename PFP::REAL vertexOneRingArea(typename PFP::MAP& map, Dart d, const typename PFP::TVEC3& position) ;
 
+template <typename PFP>
+typename PFP::REAL vertexBarycentricArea(typename PFP::MAP& map, Dart d, const typename PFP::TVEC3& position) ;
+
 template <typename PFP>
 typename PFP::REAL vertexVoronoiArea(typename PFP::MAP& map, Dart d, const typename PFP::TVEC3& position) ;
 
@@ -55,6 +58,9 @@ void computeAreaFaces(typename PFP::MAP& map, const typename PFP::TVEC3& positio
 template <typename PFP>
 void computeOneRingAreaVertices(typename PFP::MAP& map, const typename PFP::TVEC3& position, typename PFP::TREAL& vertex_area, const FunctorSelect& select = SelectorTrue()) ;
 
+template <typename PFP>
+void computeBarycentricAreaVertices(typename PFP::MAP& map, const typename PFP::TVEC3& position, typename PFP::TREAL& vertex_area, const FunctorSelect& select = SelectorTrue()) ;
+
 template <typename PFP>
 void computeVoronoiAreaVertices(typename PFP::MAP& map, const typename PFP::TVEC3& position, typename PFP::TREAL& vertex_area, const FunctorSelect& select = SelectorTrue()) ;
 

include/Algo/Geometry/area.hpp

@@ -97,7 +97,7 @@ typename PFP::REAL vertexOneRingArea(typename PFP::MAP& map, Dart d, const typen
 }
 
 template <typename PFP>
-typename PFP::REAL vertexVoronoiArea(typename PFP::MAP& map, Dart d, const typename PFP::TVEC3& position)
+typename PFP::REAL vertexBarycentricArea(typename PFP::MAP& map, Dart d, const typename PFP::TVEC3& position)
 {
 	typename PFP::REAL area(0) ;
 	Dart it = d ;
@@ -109,6 +109,32 @@ typename PFP::REAL vertexVoronoiArea(typename PFP::MAP& map, Dart d, const typen
 	return area ;
 }
 
+template <typename PFP>
+typename PFP::REAL vertexVoronoiArea(typename PFP::MAP& map, Dart d, const typename PFP::TVEC3& position)
+{
+	typename PFP::REAL area(0) ;
+	Dart it = d ;
+	do
+	{
+		if(!isTriangleObtuse<PFP>(map, it, position))
+		{
+			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 ;
+		}
+		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)
+				area += tArea / 2 ;
+			else
+				area += tArea / 4 ;
+		}
+		it = map.alpha1(it) ;
+	} while(it != d) ;
+	return area ;
+}
+
 template <typename PFP>
 void computeAreaFaces(typename PFP::MAP& map, const typename PFP::TVEC3& position, typename PFP::TREAL& face_area, const FunctorSelect& select)
 {
@@ -137,6 +163,20 @@ void computeOneRingAreaVertices(typename PFP::MAP& map, const typename PFP::TVEC
 	}
 }
 
+template <typename PFP>
+void computeBarycentricAreaVertices(typename PFP::MAP& map, const typename PFP::TVEC3& position, typename PFP::TREAL& vertex_area, const FunctorSelect& select)
+{
+	CellMarker marker(map, VERTEX);
+	for(Dart d = map.begin(); d != map.end(); map.next(d))
+	{
+		if(select(d) && !marker.isMarked(d))
+		{
+			marker.mark(d);
+			vertex_area[d] = vertexBarycentricArea<PFP>(map, d, position) ;
+		}
+	}
+}
+
 template <typename PFP>
 void computeVoronoiAreaVertices(typename PFP::MAP& map, const typename PFP::TVEC3& position, typename PFP::TREAL& vertex_area, const FunctorSelect& select)
 {

include/Algo/Geometry/basic.h

@@ -62,6 +62,12 @@ inline float angle(typename PFP::MAP& map, Dart d1, Dart d2, const typename PFP:
 	return Geom::angle(v1, v2) ;
 }
 
+template <typename PFP>
+bool isTriangleObtuse(typename PFP::MAP& map, Dart d, const typename PFP::TVEC3& position)
+{
+	return Geom::isTriangleObtuse(position[d], position[map.phi1(d)], position[map.phi_1(d)]) ;
+}
+
 } // namespace Geometry
 
 } // namespace Algo

include/Geometry/basic.h

@@ -95,12 +95,26 @@ typename VEC3::DATA_TYPE triangleArea(const VEC3& p1, const VEC3& p2, const VEC3
 	return 0.5 * ((p2 - p1) ^ (p3 - p1)).norm() ;
 }
 
+// normal of the plane spanned by 3 points in 3D
 template <typename VEC3>
 VEC3 triangleNormal(const VEC3& p1, const VEC3& p2, const VEC3& p3)
 {
 	return (p2 - p1) ^ (p3 - p1) ;
 }
 
+// return true if the triangle formed by 3 points in 3D is obtuse, false otherwise
+template <typename VEC3>
+bool isTriangleObtuse(const VEC3& p1, const VEC3& p2, const VEC3& p3)
+{
+	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) ;
+	if(a2 > M_PI / 2 || a1 + a2 < M_PI / 2)
+		return true ;
+	return false ;
+}
+
 // signed volume of the tetrahedron formed by 4 points in 3D
 template <typename VEC3>
 typename VEC3::DATA_TYPE tetraSignedVolume(const VEC3& p1, const VEC3& p2, const VEC3& p3, const VEC3& p4)

include/Geometry/transfo.h

@@ -22,15 +22,15 @@
  *                                                                              *
  *******************************************************************************/
 
-#ifndef __TRANSFO__
-#define __TRANSFO__
+#ifndef __TRANSFO__H__
+#define __TRANSFO__H__
 
 #include "Geometry/matrix.h"
 #include <cmath>
 
-
 namespace CGoGN
 {
+
 namespace Geom
 {
 
@@ -44,7 +44,6 @@ namespace Geom
 template <typename T>
 void scale(T sx, T sy, T sz, Matrix<4,4,T>& mat);
 
-
 /**
  * Apply a translation to matrix
  * @param tx scale in x axis
@@ -55,16 +54,14 @@ void scale(T sx, T sy, T sz, Matrix<4,4,T>& mat);
 template <typename T>
 void translate(T tx, T ty, T tz, Matrix<4,4,T>& mat);
 
-
 /**
- * Apply a rotation around Z axis  to matrix
+ * Apply a rotation around Z axis to matrix
  * @param angle angle of rotation in radian
  * @param mat current matrix
  */
 template <typename T>
 void rotateZ(T angle, Matrix<4,4,T>& mat);
 
-
 /**
  * Apply a rotation around Y axis to matrix
  * @param angle angle of rotation in radian
@@ -73,7 +70,6 @@ void rotateZ(T angle, Matrix<4,4,T>& mat);
 template <typename T>
 void rotateY(T angle, Matrix<4,4,T>& mat);
 
-
 /**
  * Apply a rotation around X axis to matrix
  * @param angle angle of rotation in radian
@@ -93,7 +89,6 @@ void rotateX(T angle, Matrix<4,4,T>& mat);
 template <typename T>
 void rotate(T axis_x, T axis_y, T axis_z, T angle, Matrix<4,4,T>& mat);
 
-
 /**
  * Apply a transformation (stored in matrix) to a 3D point
  * @param P the point to transfo
@@ -102,9 +97,10 @@ void rotate(T axis_x, T axis_y, T axis_z, T angle, Matrix<4,4,T>& mat);
 template <typename T>
 Vector<3,T> transform(const Vector<3,T>& P,const Matrix<4,4,T>& mat);
 
-}
-}
+} // namespace Geom
+
+} // namespace CGoGN
 
-#include "transfo.hpp"
+#include "Geometry/transfo.hpp"
 
 #endif