Merge cgogn:~kraemer/CGoGN

Apps/Examples/clipping.h

@@ -28,11 +28,11 @@
 #include <iostream>
 
 #include "Utils/Qt/qtSimple.h"
-#include "Utils/Qt/qtInputs.h"
+#include "Utils/qtInputs.h"
 
 #include "ui_clipping.h"
 // inclure qtui.h juste après le ui_xxx.h
-#include "Utils/Qt/qtui.h"
+#include "Utils/qtui.h"
 
 #include "Topology/generic/parameters.h"
 #include "Topology/map/embeddedMap3.h"

Apps/Examples/mcmesh.h

@@ -17,7 +17,7 @@
 * along with this library; if not, write to the Free Software Foundation,      *
 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301 USA.           *
 *                                                                              *
-* Web site: https://iggservis.u-strasbg.fr/CGoGN/                              *
+* Web site: http://cgogn.unistra.fr                                            *
 * Contact information: cgogn@unistra.fr                                        *
 *                                                                              *
 *******************************************************************************/
@@ -26,7 +26,7 @@
 
 #include "Utils/Qt/qtSimple.h"
 #include "ui_mcmesh.h"
-#include "Utils/Qt/qtui.h"
+#include "Utils/qtui.h"
 
 #include "Topology/generic/parameters.h"
 #include "Topology/map/map2.h"

Apps/Examples/viewer.h

@@ -17,7 +17,7 @@
 * along with this library; if not, write to the Free Software Foundation,      *
 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301 USA.           *
 *                                                                              *
-* Web site: http://cgogn.unistra.fr/                                  *
+* Web site: http://cgogn.unistra.fr/                                           *
 * Contact information: cgogn@unistra.fr                                        *
 *                                                                              *
 *******************************************************************************/
@@ -26,7 +26,7 @@
 
 #include "Utils/Qt/qtSimple.h"
 #include "ui_viewer.h"
-#include "Utils/Qt/qtui.h"
+#include "Utils/qtui.h"
 
 #include "Topology/generic/parameters.h"
 #include "Topology/map/embeddedMap2.h"

Apps/Examples/volumeExplorer.h

@@ -46,7 +46,7 @@
 #include "Utils/frameManipulator.h"
 
 #include "ui_volumeExplorer.h"
-#include "Utils/Qt/qtui.h"
+#include "Utils/qtui.h"
 
 
 using namespace CGoGN ;

Apps/Tuto/show_traversors.h

@@ -56,7 +56,7 @@
 
 #include "ui_show_traversors.h"
 // inclure qtui.h juste après le ui_xxx.h
-#include "Utils/Qt/qtui.h"
+#include "Utils/qtui.h"
 
 
 using namespace CGoGN ;

Apps/Tuto/tuto5.h

@@ -57,7 +57,7 @@
 
 #include "ui_tuto5.h"
 // inclure qtui.h juste après le ui_xxx.h
-#include "Utils/Qt/qtui.h"
+#include "Utils/qtui.h"
 
 
 using namespace CGoGN ;

Apps/Tuto/tuto_mt.cpp

@@ -42,7 +42,7 @@
 #include "Algo/Parallel/parallel_foreach.h"
 
 // for file input
- #include "Utils/Qt/qtInputs.h"
+#include "Utils/qtInputs.h"
 
 using namespace CGoGN ;
 

Apps/Tuto/tuto_oper2.h

@@ -40,7 +40,7 @@
 
 
 #include "ui_tuto_oper2.h"
-#include "Utils/Qt/qtui.h"
+#include "Utils/qtui.h"
 #include "Utils/Qt/qtSimple.h"
 #include "Utils/cgognStream.h"
 

Apps/Tuto/tuto_oper3.h

@@ -40,7 +40,7 @@
 
 
 #include "ui_tuto_oper3.h"
-#include "Utils/Qt/qtui.h"
+#include "Utils/qtui.h"
 #include "Utils/Qt/qtSimple.h"
 #include "Utils/cgognStream.h"
 

Apps/Tuto/tuto_orbits.h

@@ -57,7 +57,7 @@
 
 #include "ui_tuto_orbits.h"
 // inclure qtui.h juste après le ui_xxx.h
-#include "Utils/Qt/qtui.h"
+#include "Utils/qtui.h"
 
 
 using namespace CGoGN ;

include/Algo/Geometry/feature.h

@@ -25,6 +25,8 @@
 #ifndef __ALGO_GEOMETRY_FEATURE_H__
 #define __ALGO_GEOMETRY_FEATURE_H__
 
+
+
 namespace CGoGN
 {
 
@@ -34,44 +36,144 @@ namespace Algo
 namespace Geometry
 {
 
+enum
+{
+	EMPTY,
+	SEGMENT,
+	BARY
+};
+
+//#define EMPTY 0
+//#define SEGMENT 1
+//#define BARY 2
+
+typedef struct { Dart d ; float w ; } e0point ;
+typedef struct { e0point p1 ; e0point p2 ; unsigned char type ; } e0segment ;
+typedef NoMathIONameAttribute<e0segment> ridgeSegment ;
+
 template <typename PFP>
-void featureEdgeDetection(typename PFP::MAP& map, const typename PFP::TVEC3& position, CellMarker& featureEdge) ;
+void featureEdgeDetection(
+	typename PFP::MAP& map,
+	const typename PFP::TVEC3& position,
+	CellMarker& featureEdge) ;
 
 template <typename PFP>
-void computeArea(typename PFP::MAP& map, const typename PFP::TVEC3& position, typename PFP::TREAL& area, const FunctorSelect& select = allDarts, unsigned int thread = 0) ;
+void computeFaceGradient(
+	typename PFP::MAP& map,
+	const typename PFP::TVEC3& position,
+	typename PFP::TVEC3& face_gradient,
+	const typename PFP::TVEC3& face_normal,
+	const typename PFP::TREAL& kmax,
+	const typename PFP::TREAL& area,
+	const FunctorSelect& select = allDarts,
+	unsigned int thread = 0) ;
+
+//template <typename PFP>
+//void computeVertexGradient(
+//	typename PFP::MAP& map,
+//	const typename PFP::TVEC3& position,
+//	typename PFP::TVEC3& gradient,
+//	typename PFP::TVEC3& face_gradient,
+//	const typename PFP::TREAL& area,
+//	const FunctorSelect& select = allDarts,
+//	unsigned int thread = 0) ;
 
 template <typename PFP>
-void computeFaceGradient(typename PFP::MAP& map, const typename PFP::TVEC3& position, typename PFP::TVEC3& face_gradient, const typename PFP::TVEC3& face_normal, const typename PFP::TREAL& kmax, const typename PFP::TREAL& area, const FunctorSelect& select = allDarts, unsigned int thread = 0) ;
+void computeTriangleType(
+	typename PFP::MAP& map,
+	const typename PFP::TVEC3& Kmax,
+	CellMarker& regularMarker,
+	const FunctorSelect& select = allDarts,
+	unsigned int thread = 0) ;
 
 template <typename PFP>
-void computeGradient(typename PFP::MAP& map, const typename PFP::TVEC3& position, typename PFP::TVEC3& gradient, typename PFP::TVEC3& face_gradient, const typename PFP::TREAL& area, const FunctorSelect& select = allDarts, unsigned int thread = 0) ;
+void computeRidgeLines(
+	typename PFP::MAP& map,
+	const typename PFP::TVEC3& position,
+	CellMarker& regularMarker,
+	const typename PFP::TVEC3& faceGradient,
+	const typename PFP::TREAL& area,
+	const typename PFP::TVEC3& K,
+	AttributeHandler<ridgeSegment>& ridge_segments,
+	const FunctorSelect& select = allDarts,
+	unsigned int thread = 0) ;
 
 template <typename PFP>
-void computeTriangleType(typename PFP::MAP& map, const typename PFP::TVEC3& position, CellMarker& regularMarker, const FunctorSelect& select = allDarts, unsigned int thread = 0) ;
+void initRidgeSegments(
+	typename PFP::MAP& map,
+	const typename PFP::TVEC3& position,
+	AttributeHandler<ridgeSegment>& ridge_segments,
+	const FunctorSelect& select = allDarts,
+	unsigned int thread = 0) ;
 
 template <typename PFP>
-void computeCurvatureSign(typename PFP::MAP& map, const typename PFP::TVEC3& position, typename PFP::TVEC3& triangle_type, typename PFP::TVEC3& k, const FunctorSelect& select = allDarts, unsigned int thread = 0) ;
+void computeSingularTriangle(
+	typename PFP::MAP& map,
+	const typename PFP::TVEC3& position,
+	CellMarker& regularMarker,
+	typename PFP::TVEC3& k,
+	AttributeHandler<ridgeSegment>& ridge_segments,
+	const FunctorSelect& select = allDarts,
+	unsigned int thread = 0) ;
 
 template <typename PFP>
-std::vector<typename PFP::VEC3> occludingContoursDetection(typename PFP::MAP& map, const typename PFP::VEC3& cameraPosition, const typename PFP::TVEC3& position, const typename PFP::TVEC3& normal) ;
+std::vector<typename PFP::VEC3> occludingContoursDetection(
+	typename PFP::MAP& map,
+	const typename PFP::VEC3& cameraPosition,
+	const typename PFP::TVEC3& position,
+	const typename PFP::TVEC3& normal) ;
 
 template <typename PFP>
-typename PFP::TREAL faceArea(typename PFP::MAP& map, Dart d, const typename PFP::TVEC3& position) ;
+typename PFP::VEC3 faceGradient(
+	typename PFP::MAP& map,
+	Dart d,
+	const typename PFP::TVEC3& position,
+	const typename PFP::TVEC3& face_normal,
+	const typename PFP::TREAL& kmax,
+	const typename PFP::TREAL& area) ;
 
 template <typename PFP>
-typename PFP::VEC3 faceGradient(typename PFP::MAP& map, Dart d, const typename PFP::TVEC3& position, const typename PFP::TVEC3& face_normal, const typename PFP::TREAL& kmax, const typename PFP::TREAL& area) ;
+bool isTriangleRegular(typename PFP::MAP& map, Dart d, const typename PFP::TVEC3& Kmax) ;
+
+//template <typename PFP>
+//bool isInSameOctant(const typename PFP::VEC3& pos1, const typename PFP::VEC3& pos2) ;
 
 template <typename PFP>
-bool triangleType(typename PFP::MAP& map, Dart d, const typename PFP::TVEC3& position) ;
+void ridgeLines(
+	typename PFP::MAP& map,
+	Dart d,
+	const typename PFP::TVEC3& position,
+	const typename PFP::TVEC3& K,
+	const typename PFP::TVEC3& faceGradient,
+	const typename PFP::TREAL& area,
+	AttributeHandler<ridgeSegment>& ridge_segments) ;
 
 template <typename PFP>
-bool isInSameOctant(const typename PFP::VEC3& pos1, const typename PFP::VEC3& pos2) ;
+typename PFP::REAL extremality(
+	typename PFP::MAP& map,
+	Dart d,
+	const typename PFP::VEC3& K,
+	const typename PFP::TVEC3& faceGradient,
+	const typename PFP::TREAL& area) ;
+
+//template <typename PFP>
+//typename PFP::TVEC3 vertexGradient(
+//	typename PFP::MAP& map,
+//	Dart d,
+//	const typename PFP::TVEC3& position,
+//	const typename PFP::TVEC3& face_gradient,
+//	const typename PFP::TREAL& area) ;
 
 template <typename PFP>
-void curvatureSign(typename PFP::MAP& map, Dart d, const typename PFP::TVEC3& position, typename PFP::TVEC3& triangle_type, typename PFP::TVEC3& k) ;
+void singularTriangle(
+	typename PFP::MAP& map,
+	Dart d,
+	const typename PFP::TVEC3& position,
+	CellMarker& regularMarker,
+	AttributeHandler<ridgeSegment>& ridge_segments) ;
 
 template <typename PFP>
-typename PFP::TVEC3 vertexGradient(typename PFP::MAP& map, Dart d, const typename PFP::TVEC3& position, const typename PFP::TVEC3& face_gradient, const typename PFP::TREAL& area) ;
+bool isEdgeInTriangle(typename PFP::MAP& map, Dart edge, Dart triangle) ;
 
 } // namespace Geometry
 

include/Algo/Geometry/feature.hpp

@@ -36,7 +36,10 @@ namespace Geometry
 {
 
 template <typename PFP>
-void featureEdgeDetection(typename PFP::MAP& map, typename PFP::TVEC3& position, CellMarker& featureEdge)
+void featureEdgeDetection(
+	typename PFP::MAP& map,
+	typename PFP::TVEC3& position,
+	CellMarker& featureEdge)
 {
 	typedef typename PFP::VEC3 VEC3 ;
 	typedef typename PFP::REAL REAL ;
@@ -59,54 +62,106 @@ void featureEdgeDetection(typename PFP::MAP& map, typename PFP::TVEC3& position,
 }
 
 template <typename PFP>
-void computeArea(typename PFP::MAP& map, const typename PFP::TVEC3& position, typename PFP::TREAL& area, const FunctorSelect& select = allDarts, unsigned int thread = 0)
+void computeFaceGradient(
+	typename PFP::MAP& map,
+	const typename PFP::TVEC3& position,
+	typename PFP::TVEC3& face_gradient,
+	const typename PFP::TVEC3& face_normal,
+	const typename PFP::TREAL& kmax,
+	const typename PFP::TREAL& area,
+	const FunctorSelect& select,
+	unsigned int thread)
 {
 	TraversorF<typename PFP::MAP> trav(map, select, thread);
 	for (Dart d = trav.begin(); d != trav.end(); d = trav.next())
-		area[d] = faceArea<PFP>(map, d, position)[1] ;
+		face_gradient[d] = faceGradient<PFP>(map, d, position, face_normal, kmax, area) ;
 }
 
+//template <typename PFP>
+//void computeVertexGradient(
+//	typename PFP::MAP& map,
+//	const typename PFP::TVEC3& position,
+//	typename PFP::TVEC3& gradient,
+//	typename PFP::TVEC3& face_gradient,
+//	const typename PFP::TREAL& area,
+//	const FunctorSelect& select,
+//	unsigned int thread)
+//{
+//	TraversorV<typename PFP::MAP> trav(map, select, thread);
+//	for (Dart d = trav.begin(); d != trav.end(); d = trav.next())
+//		gradient[d] = vertexGradient<PFP>(map, d, position, face_gradient, area) ;
+//}
+
 template <typename PFP>
-void computeFaceGradient(typename PFP::MAP& map, const typename PFP::TVEC3& position, typename PFP::TVEC3& face_gradient, const typename PFP::TVEC3& face_normal, const typename PFP::TREAL& kmax, const typename PFP::TREAL& area, const FunctorSelect& select = allDarts, unsigned int thread = 0)
+void computeTriangleType(
+	typename PFP::MAP& map,
+	const typename PFP::TVEC3& Kmax,
+	CellMarker& regularMarker,
+	const FunctorSelect& select,
+	unsigned int thread)
 {
 	TraversorF<typename PFP::MAP> trav(map, select, thread);
 	for (Dart d = trav.begin(); d != trav.end(); d = trav.next())
-		face_gradient[d] = faceGradient<PFP>(map, d, position, face_normal, kmax, area) ;
+		if(isTriangleRegular<PFP>(map, d, Kmax))
+			regularMarker.mark(d) ;
 }
 
 template <typename PFP>
-void computeGradient(typename PFP::MAP& map, const typename PFP::TVEC3& position, typename PFP::TVEC3& gradient, typename PFP::TVEC3& face_gradient, const typename PFP::TREAL& area, const FunctorSelect& select = allDarts, unsigned int thread = 0)
+void computeRidgeLines(
+	typename PFP::MAP& map,
+	const typename PFP::TVEC3& position,
+	CellMarker& regularMarker,
+	const typename PFP::TVEC3& faceGradient,
+	const typename PFP::TREAL& area,
+	const typename PFP::TVEC3& K,
+	AttributeHandler<ridgeSegment>& ridge_segments,
+	const FunctorSelect& select,
+	unsigned int thread)
 {
-	TraversorV<typename PFP::MAP> trav(map, select, thread);
+	TraversorF<typename PFP::MAP> trav(map, select, thread);
 	for (Dart d = trav.begin(); d != trav.end(); d = trav.next())
-		gradient[d] = vertexGradient<PFP>(map, d, position, face_gradient, area) ;
+	{
+		if (regularMarker.isMarked(d))
+			ridgeLines<PFP>(map, d, position, K, faceGradient, area, ridge_segments) ;
+	}
 }
 
 template <typename PFP>
-void computeTriangleType(typename PFP::MAP& map, const typename PFP::TVEC3& position, CellMarker& regularMarker, const FunctorSelect& select = allDarts, unsigned int thread = 0)
+void initRidgeSegments(
+	typename PFP::MAP& map,
+	const typename PFP::TVEC3& position,
+	AttributeHandler<ridgeSegment>& ridge_segments,
+	const FunctorSelect& select,
+	unsigned int thread)
 {
 	TraversorF<typename PFP::MAP> trav(map, select, thread);
 	for (Dart d = trav.begin(); d != trav.end(); d = trav.next())
-		if(triangleType<PFP>(map, d, position))
-			regularMarker.mark(d) ;
+		ridge_segments[d].type = EMPTY ;
 }
 
 template <typename PFP>
-void computeCurvatureSign(typename PFP::MAP& map, const typename PFP::TVEC3& position, typename PFP::TVEC3& triangle_type, typename PFP::TVEC3& k, const FunctorSelect& select = allDarts, unsigned int thread = 0)
+void computeSingularTriangle(
+	typename PFP::MAP& map,
+	const typename PFP::TVEC3& position,
+	CellMarker& regularMarker,
+	AttributeHandler<ridgeSegment>& ridge_segments,
+	const FunctorSelect& select = allDarts,
+	unsigned int thread = 0)
 {
-	CellMarker mv(map, VERTEX, thread);
-	TraversorV<typename PFP::MAP> trav(map, select, thread);
+	TraversorF<typename PFP::MAP> trav(map, select, thread);
 	for (Dart d = trav.begin(); d != trav.end(); d = trav.next())
 	{
-		if (!mv.isMarked(d))
-		{
-			curvatureSign<PFP>(map, d, position, triangle_type, k) ;
-		}
+		if (! regularMarker.isMarked(d))
+			singularTriangle<PFP>(map, d, position, regularMarker, ridge_segments) ;
 	}
 }
 
 template <typename PFP>
-std::vector<typename PFP::VEC3> occludingContoursDetection(typename PFP::MAP& map, const typename PFP::VEC3& cameraPosition, const typename PFP::TVEC3& position, const typename PFP::TVEC3& normal)
+std::vector<typename PFP::VEC3> occludingContoursDetection(
+	typename PFP::MAP& map,
+	const typename PFP::VEC3& cameraPosition,
+	const typename PFP::TVEC3& position,
+	const typename PFP::TVEC3& normal)
 {
 	typedef typename PFP::VEC3 VEC3 ;
 	typedef typename PFP::REAL REAL ;
@@ -160,185 +215,329 @@ std::vector<typename PFP::VEC3> occludingContoursDetection(typename PFP::MAP& ma
 }
 
 template <typename PFP>
-float faceArea(typename PFP::MAP& map, Dart d, const typename PFP::TVEC3& position)
+typename PFP::VEC3 faceGradient(
+	typename PFP::MAP& map,
+	Dart d,
+	const typename PFP::TVEC3& position,
+	const typename PFP::TVEC3& face_normal,
+	const typename PFP::TREAL& kmax,
+	const typename PFP::TREAL& area)
 {
 	typedef typename PFP::REAL REAL ;
 	typedef typename PFP::VEC3 VEC3 ;
 
-	float A [1] ;
-
-	VEC3 pos1(0) ;
-	VEC3 pos2(0) ;
-	VEC3 pos3(0) ;
-
 	Traversor2FV<typename PFP::MAP> t(map, d) ;
+
 	Dart it = t.begin() ;
+	VEC3 pos1 = position[it] ;
+	REAL k1 = kmax[it] ;
+
+	it = t.next() ;
+	VEC3 pos2 = position[it] ;
+	REAL k2 = kmax[it] ;
+
+	it = t.next() ;
+	VEC3 pos3 = position[it] ;
+	REAL k3 = kmax[it] ;
 
-	pos1 += position[it] ;
-	it = t.next ;
-	pos2 += position[it] ;
-	it = t.next ;
-	pos3 += position[it] ;
+	VEC3 n = face_normal[d] ;
+	REAL a = area[d] ;
+	VEC3 G = k1 * ( ( n ^ ( pos3 - pos2 ) ) / ( 2 * a ) ) +
+			 k2 * ( ( n ^ ( pos1 - pos3 ) ) / ( 2 * a ) ) +
+			 k3 * ( ( n ^ ( pos2 - pos1 ) ) / ( 2 * a ) ) ;
 
-	A = 1/2 * abs(
-				( (pos2.x - pos1.x) * (pos3.y - pos1.y) )
-				-
-				( (pos3.x - pos1.x) * (pos2.y - pos1.y) )
-				) ;
+	G.normalize() ;
+	return G ;
+}
 
-	return A ;
+template <typename PFP>
+bool mutuallyPositive(typename PFP::VEC3& v1, typename PFP::VEC3& v2, typename PFP::VEC3& v3)
+{
+	typename PFP::REAL v1v2 = v1 * v2 ;
+	typename PFP::REAL v1v3 = v1 * v3 ;
+	typename PFP::REAL v2v3 = v2 * v3 ;
+	if(v1v2 > 0 && v1v3 > 0 && v2v3 > 0)
+		return true ;
+	return false ;
 }
 
 template <typename PFP>
-typename PFP::VEC3 faceGradient(typename PFP::MAP& map, Dart d, const typename PFP::TVEC3& position, const typename PFP::TVEC3& face_normal, const typename PFP::TREAL& kmax, const typename PFP::TREAL& area)
+bool isTriangleRegular(typename PFP::MAP& map, Dart d, const typename PFP::TVEC3& Kmax)
 {
 	typedef typename PFP::REAL REAL ;
 	typedef typename PFP::VEC3 VEC3 ;
 
-	VEC3 G(0) ;
-
-	VEC3 pos1(0) ;
-	VEC3 pos2(0) ;
-	VEC3 pos3(0) ;
-
-	VEC3 n = face_normal[d] ;
-
-	REAL k1 ;
-	REAL k2 ;
-	REAL k3 ;
+	Dart v1 = d ;
+	Dart v2 = map.phi1(v1) ;
+	Dart v3 = map.phi1(v2) ;
 
-	REAL a = area[d] ;
+	VEC3 K1 = Kmax[v1] ;
+	VEC3 K2 = Kmax[v2] ;
+	VEC3 K3 = Kmax[v3] ;
 
-	Traversor2FV<typename PFP::MAP> t(map, d) ;
-	Dart it = t.begin() ;
-
-	pos1 += position[it] ;
-	k1 = kmax[it] ;
-	it = t.next() ;
-	pos2 += position[it] ;
-	k2 = kmax[it] ;
-	it = t.next() ;
-	pos3 += position[it] ;
-	k3 = kmax[it] ;
+//	VEC3 K1n = typename VEC3::DATA_TYPE(-1) * K1 ;
+	VEC3 K2n = typename VEC3::DATA_TYPE(-1) * K2 ;
+	VEC3 K3n = typename VEC3::DATA_TYPE(-1) * K3 ;
 
-	G += k1 * ( ( n ^ ( pos3 - pos2 ) ) / ( 2 * a ) ) +
-			k2 * ( ( n ^ ( pos1 - pos3 ) ) / ( 2 * a ) ) +
-			k3 * ( ( n ^ ( pos2 - pos1 ) ) / ( 2 * a ) ) ;
+	if(mutuallyPositive<PFP>(K1, K2, K3))
+		return true ;
+	if(mutuallyPositive<PFP>(K1, K2, K3n))
+		return true ;
+	if(mutuallyPositive<PFP>(K1, K2n, K3))
+		return true ;
+	if(mutuallyPositive<PFP>(K1, K2n, K3n))
+		return true ;
 
-	G.normalize() ;
-	return G ;
+//	if(mutuallyPositive<PFP>(K1n, K2, K3))
+//		return true ;
+//	if(mutuallyPositive<PFP>(K1n, K2, K3n))
+//		return true ;
+//	if(mutuallyPositive<PFP>(K1n, K2n, K3))
+//		return true ;
+//	if(mutuallyPositive<PFP>(K1n, K2n, K3n))
+//		return true ;
+
+	return false ;
+
+//	if ( isInSameOctant<PFP>(K1, K2) && isInSameOctant<PFP>(K1, K3) )
+//	{
+//		//regular triangle
+//		return true ;
+//	}
+//	else
+//	{
+//		//singular triangle
+//		return false ;
+//	}
 }
 
+//template <typename PFP>
+//bool isInSameOctant(const typename PFP::VEC3& pos1, const typename PFP::VEC3& pos2)
+//{
+//	return ( ( pos1[0] > 0 && pos1[1] > 0 && pos1[2] > 0 && //one of the eight octant
+//			( ( pos2[0] > 0 && pos2[1] > 0 && pos2[2] > 0 ) or //one choice of sign
+//			( pos2[0] < 0 && pos2[1] < 0 && pos2[2] < 0 ) ) ) or //the other choice
+//
+//			( pos1[0] > 0 && pos1[1] > 0 && pos1[2] < 0 &&
+//			( ( pos2[0] > 0 && pos2[1] > 0 && pos2[2] < 0 ) or
+//			( pos2[0] < 0 && pos2[1] < 0 && pos2[2] > 0 ) ) ) or
+//
+//			( pos1[0] > 0 && pos1[1] < 0 && pos1[2] > 0 &&
+//			( ( pos2[0] > 0 && pos2[1] < 0 && pos2[2] > 0 ) or
+//			( pos2[0] < 0 && pos2[1] > 0 && pos2[2] < 0 ) ) ) or
+//
+//			( pos1[0] > 0 && pos1[1] < 0 && pos1[2] < 0 &&
+//			( ( pos2[0] > 0 && pos2[1] < 0 && pos2[2] < 0 ) or
+//			( pos2[0] < 0 && pos2[1] > 0 && pos2[2] > 0 ) ) ) or
+//
+//			( pos1[0] < 0 && pos1[1] > 0 && pos1[2] > 0 &&
+//			( ( pos2[0] < 0 && pos2[1] > 0 && pos2[2] > 0 ) or