correction of isPointInTetrahedron

include/Algo/Geometry/orientation.h

@@ -34,6 +34,9 @@ namespace CGoGN
 namespace Algo
 {
 
+namespace Surface
+{
+
 namespace Geometry
 {
 
@@ -53,6 +56,8 @@ bool isTetrahedronWellOriented(typename PFP::MAP& map, Dart d, const VertexAttri
 
 } // namespace Geometry
 
+} // namespace Surface
+
 } // namespace Algo
 
 } // namespace CGoGN

include/Algo/Geometry/plane.h

@@ -35,6 +35,9 @@ namespace CGoGN
 namespace Algo
 {
 
+namespace Surface
+{
+
 namespace Geometry
 {
 
@@ -61,6 +64,8 @@ Geom::Plane3D<typename PFP::REAL> vertexTangentPlane(typename PFP::MAP& map, Dar
 
 } // namespace Geometry
 
+} // namespace Surface
+
 } // namespace Algo
 
 } // namespace CGoGN

include/Geometry/inclusion.hpp

@@ -22,6 +22,9 @@
  *                                                                              *
  *******************************************************************************/
 
+#include <Eigen/Dense>
+#include "Utils/convertType.h"
+
 namespace CGoGN
 {
 
@@ -128,23 +131,18 @@ bool isPointInTetrahedron(VEC3 points[4], VEC3& point, bool CCW)
 	VEC3 AB = points[1] - points[0] ;
 	VEC3 AC = points[2] - points[0] ;
 	VEC3 AD = points[3] - points[0] ;
-	VEC3 BD = points[3] - points[1] ;
-	VEC3 BC = points[2] - points[1] ;
-
-	VEC3 nABC = AB ^ AC ;
-	VEC3 nACD = AC ^ AD ;
-	VEC3 nADB = AD ^ AB ;
-	VEC3 nBDC = BD ^ BC ;
-
-	T d0 = nABC * (point - points[0]) ;
-	T d1 = nACD * (point - points[0]) ;
-	T d2 = nADB * (point - points[0]) ;
-	T d3 = nBDC * (point - points[1]) ;
-
-	if(CCW)
-		return (d0 < 0 && d1 < 0 && d2 < 0 && d3 < 0) ;
-	else
-		return (d0 > 0 && d1 > 0 && d2 > 0 && d3 > 0) ;
+
+	Eigen::Matrix3f A;
+	A << AB[0], AB[1], AB[2],
+			AC[0], AC[1], AC[2],
+			AD[0], AD[1], AD[2];
+
+	A = A.inverse();
+
+	VEC3 v1(point-points[0]);
+	Eigen::Vector3f& v = Utils::convertRef<Eigen::Vector3f>(v1);
+	Eigen::Vector3f beta = A* v;
+	return (beta[0] >= 0.0f && beta[1] >= 0.0f && beta[2] >= 0.0f && beta[0]+beta[1]+beta[2]<=1.0f);
 }
 
 template <typename VEC3>

include/Topology/generic/dartmarker.h

@@ -312,8 +312,7 @@ public:
 	void unmarkAll()
 	{
 		assert(m_map.getMarkerSet<DART>(m_thread).testMark(m_mark));
-		for (std::vector<unsigned int>::iterator it = m_markedDarts.begin();
-		    it != m_markedDarts.end(); ++it)
+		for (std::vector<unsigned int>::iterator it = m_markedDarts.begin(); it != m_markedDarts.end(); ++it)
 			m_markVector->operator[](*it).unsetMark(m_mark) ;
 	}
 } ;