Merge branch 'master' of cgogn:~cgogn/CGoGN

include/Algo/Modelisation/polyhedron.hpp

@@ -431,7 +431,7 @@ Dart Polyhedron<PFP>::cylinder_topo(unsigned int n, unsigned int z, bool top_clo
 	if (top_closed)
 	{
 		Dart d =  m_map.phi_1(m_tableVertDarts[n*z]);
-		if(m_map.closeHole(d, true))
+		if(m_map.closeHole(d, false))
 		{
 			d = m_map.phi2(d);
 			if(m_map.faceDegree(d) > 3)
@@ -597,7 +597,7 @@ Dart Polyhedron<PFP>::cube_topo(unsigned int x, unsigned int y, unsigned int z)
 {
 	if (m_kind != NONE) return m_dart;
 
-	m_dart = cylinder_topo(2*(x+y),z, false,false);
+	m_dart = cylinder_topo(2*(x+y), z, false, false);
 	m_kind = CUBE;
 	m_nx = x;
 	m_ny = y;
@@ -697,17 +697,17 @@ Dart Polyhedron<PFP>::tore_topo(unsigned int m, unsigned int n)
 	if (m_kind != NONE) return m_dart;
 
 	m_dart = cylinder_topo(n, m, false, false);
-	m_nx=n;
-	m_ny=m;
+	m_nx = n;
+	m_ny = m;
 	m_kind = TORE;
 
-	// juste finish to sew
+	// just finish to sew
 	for(unsigned int i = 0; i < n; ++i)
 	{
 		Dart d = m_tableVertDarts[i];
 		Dart e = m_tableVertDarts[(m*n)+i];
 		e = m_map.phi_1(e);
-		m_map.sewFaces(d, e, true);
+		m_map.sewFaces(d, e);
 	}
 
 	// remove the last n vertex darts that are no more necessary (sewed with n first)
@@ -722,7 +722,8 @@ void Polyhedron<PFP>::embedGrid(float x, float y, float z)
 {
 	typedef typename PFP::VEC3 VEC3 ;
 
-	if (m_kind != GRID) {
+	if (m_kind != GRID)
+	{
 		CGoGNerr << "Warning try to embedGrid something that is not a grid"<<CGoGNendl;
 		return;
 	}

include/Algo/MovingObjects/particle_cell_2DandHalf.h

 #ifndef PARTCELL_H
 #define PARTCELL_H
 
-#include "particle_base.h"
+#include "Algo/MovingObjects/particle_base.h"
 
 #include "Algo/Geometry/inclusion.h"
+#include "Algo/Geometry/plane.h"
 #include "Geometry/intersection.h"
 #include "Geometry/orientation.h"
 #include <iostream>
@@ -20,7 +21,8 @@ namespace Algo
 namespace MovingObjects
 {
 
-enum {
+enum
+{
 	NO_CROSS,
 	CROSS_EDGE,
 	CROSS_OTHER
@@ -45,15 +47,19 @@ public :
 
 	unsigned int crossCell ;
 
+	float distance;
+
 	ParticleCell2DAndHalf(Map& map) : m(map)
 	{}
 
 	ParticleCell2DAndHalf(Map& map, Dart belonging_cell, VEC3 pos, const TAB_POS& tabPos) :
-		ParticleBase(pos), m(map), m_positions(tabPos), d(belonging_cell), lastCrossed(belonging_cell), state(FACE), crossCell(NO_CROSS)
+		ParticleBase(pos), m(map), m_positions(tabPos), d(belonging_cell), lastCrossed(belonging_cell), state(FACE), crossCell(NO_CROSS), distance(0)
 	{}
 
 	Dart getCell() { return d; }
 
+	float getDistance() { return distance; }
+
 	Geom::Orientation3D getOrientationEdge(const VEC3& point, Dart d);
 
 	void display();
@@ -77,28 +83,29 @@ public :
 
 	void move(const VEC3& newCurrent)
 	{
+		distance = 0 ;
 		crossCell = NO_CROSS ;
 		if(!Geom::arePointsEquals(newCurrent, m_position))
 		{
 			switch(state) {
-			case VERTEX : 	vertexState(newCurrent); break;
-			case EDGE : 	edgeState(newCurrent);   break;
-			case FACE : 	faceState(newCurrent);   break;
+				case VERTEX : 	vertexState(newCurrent); break;
+				case EDGE : 	edgeState(newCurrent);   break;
+				case FACE : 	faceState(newCurrent);   break;
 			}
 
-			display();
+//			display();
 		}
 		else
 			m_position = newCurrent;
 	}
 };
 
-#include "particle_cell_2DandHalf.hpp"
+} // namespace MovingObjects
 
-}
+} // namespace Algo
 
-}
+} // namespace CGoGN
 
-}
+#include "Algo/MovingObjects/particle_cell_2DandHalf.hpp"
 
 #endif

include/Algo/MovingObjects/particle_cell_2DandHalf.hpp

 //#define DEBUG
 
+#include "Geometry/frame.h"
+
+namespace CGoGN
+{
+
+namespace Algo
+{
+
+namespace MovingObjects
+{
+
 template <typename PFP>
 void ParticleCell2DAndHalf<PFP>::display()
 {
@@ -11,50 +22,50 @@ void ParticleCell2DAndHalf<PFP>::display()
 template <typename PFP>
 typename PFP::VEC3 ParticleCell2DAndHalf<PFP>::pointInFace(Dart d)
 {
-	const VEC3& p1(m_positions[d]);
-	Dart dd=m.phi1(d);
-	const VEC3& p2(m_positions[dd]);
-	dd=m.phi1(dd);
-	VEC3& p3(m_positions[dd]);
+	const VEC3& p1(m_positions[d]) ;
+	Dart dd = m.phi1(d) ;
+	const VEC3& p2(m_positions[dd]) ;
+	dd = m.phi1(dd) ;
+	VEC3& p3(m_positions[dd]) ;
 
-	VEC3 v1(p2-p1);
+	VEC3 v1(p2 - p1) ;
 
-	while((v1^VEC3(p3-p1)).norm2()==0.0f)
+	while ((v1 ^ VEC3(p3 - p1)).norm2() == 0.0f)
 	{
-		dd = m.phi1(dd);
-		p3 = m_positions[dd];
+		dd = m.phi1(dd) ;
+		p3 = m_positions[dd] ;
 	}
 
-	CGoGNout << "pointInFace " << (p1+p3)*0.5f << CGoGNendl;
+	CGoGNout << "pointInFace " << (p1 + p3) * 0.5f << CGoGNendl ;
 
-	return (p1+p3)*0.5f;
+	return (p1 + p3) * 0.5f ;
 }
 
 template <typename PFP>
 Geom::Orientation3D ParticleCell2DAndHalf<PFP>::getOrientationEdge(const VEC3& point, Dart d)
 {
-	const VEC3& endPoint = m_positions[m.phi2(d)];
+	const VEC3& endPoint = m_positions[m.phi1(d)];
 	const VEC3& vertexPoint = m_positions[d];
 
-	const VEC3& n1 = Algo::Geometry::faceNormal<PFP>(m,d,m_positions);
+	const VEC3& n1 = Algo::Geometry::faceNormal<PFP>(m, d, m_positions);
 
 	//orientation relative to the plane orthogonal to the face going through the edge
-	return Geom::testOrientation3D(point,vertexPoint, endPoint, vertexPoint+n1);
+	return Geom::testOrientation3D(point, vertexPoint, endPoint, vertexPoint+n1);
 }
 
 template <typename PFP>
 typename PFP::VEC3 ParticleCell2DAndHalf<PFP>::intersectLineEdge(const VEC3& pA, const VEC3& pB, Dart d)
 {
 	const VEC3& q1 = m_positions[d];
-	const VEC3& q2 = m_positions[m.phi2(d)];
+	const VEC3& q2 = m_positions[m.phi1(d)];
 	VEC3 Inter;
 
-	VEC3 n1 = Algo::Geometry::faceNormal<PFP>(m,d,m_positions);
-	VEC3 n = (q2-q1) ^ n1 ;
+	VEC3 n1 = Algo::Geometry::faceNormal<PFP>(m, d, m_positions);
+	VEC3 n = (q2 - q1) ^ n1 ;
 
-	Geom::intersectionLinePlane(pA,pB-pA,q1,n,Inter);
+	Geom::intersectionLinePlane(pA, pB - pA, q1, n, Inter) ;
 
-	Geom::Plane3D<float> pl = Algo::Geometry::facePlane<PFP>(m,d,m_positions);
+	Geom::Plane3D<float> pl = Algo::Geometry::facePlane<PFP>(m, d, m_positions);
 	pl.project(Inter);
 
 	return Inter;
@@ -65,7 +76,7 @@ Geom::Orientation3D ParticleCell2DAndHalf<PFP>::getOrientationFace(VEC3 point, V
 {
 	const VEC3& dPoint = m_positions[d];
 
-	VEC3 n1 = Algo::Geometry::faceNormal<PFP>(m,d,m_positions);
+	VEC3 n1 = Algo::Geometry::faceNormal<PFP>(m, d, m_positions);
 
 	return Geom::testOrientation3D(point, sourcePoint, dPoint+n1, dPoint);
 }
@@ -163,32 +174,31 @@ void ParticleCell2DAndHalf<PFP>::edgeState(VEC3 current, Geom::Orientation3D sid
 	switch(sideOfEdge)
 	{
 		case Geom::UNDER :
+		{
 			d = m.phi1(d);
 			faceState(current);
 			return;
+		}
+
 		case Geom::OVER:
 		{
 			//transform the displacement into the new entered face
-			VEC3 displ = current-m_position;
-			VEC3 edge = Algo::Geometry::vectorOutOfDart<PFP>(m,m.phi2(d),m_positions);
-
-			VEC3 n1 = Algo::Geometry::faceNormal<PFP>(m,d,m_positions);
-			VEC3 n2 = Algo::Geometry::faceNormal<PFP>(m,m.phi2(d),m_positions);
+			VEC3 displ = current - m_position;
 
-			float angle = acos(n1*n2);
+			VEC3 n1 = Algo::Geometry::faceNormal<PFP>(m, d, m_positions);
+			VEC3 n2 = Algo::Geometry::faceNormal<PFP>(m, m.phi2(d), m_positions);
+			VEC3 axis = n1 ^ n2 ;
 
-			Geom::Matrix<4,4,float> mRot;
-			mRot.identity();
+			float angle = Geom::angle(n1, n2) ;
 
-			Geom::rotate(edge[0],edge[1],edge[2],angle,mRot);
-
-			displ = Geom::transform(displ,mRot);
-			current = m_position+displ;
+			displ = Geom::rotate(axis, angle, displ) ;
+			current = m_position + displ;
 
 			d = m.phi1(m.phi2(d));
 			faceState(current);
 			return;
 		}
+
 		default :
 			state = EDGE;
 	}
@@ -221,14 +231,14 @@ void ParticleCell2DAndHalf<PFP>::faceState(VEC3 current)
  	assert(std::isfinite(current[0]) && std::isfinite(current[1]) && std::isfinite(current[2]));
 // 	assert(Algo::Geometry::isPointInConvexFace2D<PFP>(m,d,m_positions,m_position,true));
 
-	//project current within plane
+	//project current within face plane
 	VEC3 n1 = Algo::Geometry::faceNormal<PFP>(m,d,m_positions);
-	VEC3 n2 = current-m_position;
-	n1.normalize();
-	VEC3 n3 = n1^n2;
+	VEC3 n2 = current - m_position;
+//	n1.normalize();
+	VEC3 n3 = n1 ^ n2;
 	n3.normalize();
-	VEC3 n4 = n3^n1;
-	current = m_position+(n2*n4)*n4;
+	VEC3 n4 = n3 ^ n1;
+	current = m_position + (n2 * n4) * n4;
 
 	//track new position within map
 	Dart dd = d;
@@ -326,8 +336,9 @@ void ParticleCell2DAndHalf<PFP>::faceState(VEC3 current)
 	switch (getOrientationEdge(current, d))
 	{
 	case Geom::UNDER :
+		distance += (current - m_position).norm();
 		m_position = current;
-		state = FACE;;
+		state = FACE;
 		break;
 	default :
 		if(wsoe == Geom::ON)
@@ -338,14 +349,22 @@ void ParticleCell2DAndHalf<PFP>::faceState(VEC3 current)
 //
 //			vertexState(current);
 		}
-//		else
+		else
 		{
 // 			CGoGNout << "wsoe : " << wsoe << CGoGNendl;
 // 			CGoGNout << "current " << current << " " << m_position << CGoGNendl;
 // 			CGoGNout << "d " << d << "d1 " << m_positions[d] << " d2 " << m_positions[m.phi2(d)] << CGoGNendl;
-			m_position = intersectLineEdge(current, m_position, d);
+			VEC3 isect = intersectLineEdge(current, m_position, d);
+			distance += (isect - m_position).norm();
+			m_position = isect;
 // 			CGoGNout << " inter : " << m_position << CGoGNendl;
 			edgeState(current, Geom::OVER);
 		}
 	}
 }
+
+} // namespace MovingObjects
+
+} // namespace Algo
+
+} // namespace CGoGN

include/Geometry/frame.h

@@ -22,14 +22,16 @@
 *                                                                              *
 *******************************************************************************/
 
-#ifndef LOCALFRAME_H_
-#define LOCALFRAME_H_
+#ifndef _FRAME_H_
+#define _FRAME_H_
 
 #include <cmath>
 
-namespace CGoGN {
+namespace CGoGN
+{
 
-namespace Geom {
+namespace Geom
+{
 
 /**
  * Util for rotation of a 3D point (or vector) around a given line (going through the origin) and of a given angle
@@ -213,10 +215,10 @@ private : // private constants
 
 } ;
 
-} // Geom
+} // namespace Geom
 
-} // CGoGN
+} // namespace CGoGN
 
-#include "frame.hpp"
+#include "Geometry/frame.hpp"
 
-#endif /* LOCALFRAME_H_ */
+#endif /* _FRAME_H_ */

include/Geometry/intersection.h

@@ -43,7 +43,6 @@ enum Intersection
 	FACE_INTERSECTION = 3
 } ;
 
-
 /**
  * test the intersection between a line and a triangle
  * @param P a point on the line
@@ -54,9 +53,7 @@ enum Intersection
  * @return the intersection ( FACE_INTERSECTION = OK, EDGE_INTERSECTION = line inside of plane)
  */
 template <typename VEC3>
-Intersection intersectionLinePlane(const VEC3& P, const VEC3& Dir, const VEC3& PlaneP,  const VEC3& NormP, VEC3& Inter) ;
-
-
+Intersection intersectionLinePlane(const VEC3& P, const VEC3& Dir, const VEC3& PlaneP, const VEC3& NormP, VEC3& Inter) ;
 
 /**
  * test the intersection between a ray and a triangle (optimized version with triple product
@@ -84,8 +81,6 @@ Intersection intersectionRayTriangle(const VEC3& P, const VEC3& Dir, const VEC3&
 template <typename VEC3>
 Intersection intersectionRayTriangleOpt(const VEC3& P, const VEC3& Dir, const VEC3& Ta,  const VEC3& Tb, const VEC3& Tc, VEC3& Inter);
 
-
-
 /**
  * test the intersection between a ray and a triangle (optimized version with triple product
  * @param P a point on the line
@@ -98,7 +93,6 @@ Intersection intersectionRayTriangleOpt(const VEC3& P, const VEC3& Dir, const VE
 template <typename VEC3>
 Intersection intersectionRayTriangleOpt(const VEC3& P, const VEC3& Dir, const VEC3& Ta,  const VEC3& Tb, const VEC3& Tc);
 
-
 /**
  * test the intersection between a ray and a triangle, but test only face intersection and
  * @param P a point on the line

include/Geometry/intersection.hpp

@@ -28,32 +28,30 @@ namespace CGoGN
 namespace Geom
 {
 
-
 template <typename VEC3>
-Intersection intersectionLinePlane(const VEC3& P, const VEC3& Dir, const VEC3& PlaneP,  const VEC3& NormP, VEC3& Inter)
+Intersection intersectionLinePlane(const VEC3& P, const VEC3& Dir, const VEC3& PlaneP, const VEC3& NormP, VEC3& Inter)
 {
-	float b = NormP * Dir;
+	float b = NormP * Dir ;
 
-	#define PRECISION 1e-20
-	if(fabs(b) < PRECISION)		//ray parallel to triangle
+#define PRECISION 1e-20
+	if (fabs(b) < PRECISION)		//ray parallel to triangle
 	{
-		VEC3 v= PlaneP - P;
-		float c = NormP * v;
-		if(fabs(c) < PRECISION)
-			return EDGE_INTERSECTION;
+		VEC3 v = PlaneP - P ;
+		float c = NormP * v ;
+		if (fabs(c) < PRECISION )
+		return EDGE_INTERSECTION ;
 		return NO_INTERSECTION ;
 	}
-    #undef PRECISION
+#undef PRECISION
 
-	float a = NormP * (PlaneP - P);
+	float a = NormP * (PlaneP - P) ;
 
-	Inter = P + (a/b)*Dir;
-	return FACE_INTERSECTION;
+	Inter = P + (a / b) * Dir ;
+	return FACE_INTERSECTION ;
 }
 
-
 template <typename VEC3>
-Intersection intersectionRayTriangleOpt(const VEC3& P, const VEC3& Dir, const VEC3& Ta,  const VEC3& Tb, const VEC3& Tc, VEC3& Inter)
+Intersection intersectionRayTriangleOpt(const VEC3& P, const VEC3& Dir, const VEC3& Ta, const VEC3& Tb, const VEC3& Tc, VEC3& Inter)
 {
 	typedef typename VEC3::DATA_TYPE T ;
 
@@ -65,33 +63,32 @@ Intersection intersectionRayTriangleOpt(const VEC3& P, const VEC3& Dir, const VE
 	T y = tripleProduct(Dir, v, w) ;
 	T z = tripleProduct(Dir, w, u) ;
 
-	unsigned int np=0;
-	unsigned int nn=0;
-	unsigned int nz=0;
+	unsigned int np = 0 ;
+	unsigned int nn = 0 ;
+	unsigned int nz = 0 ;
 
-	if (x>T(0))
-		++np;
+	if (x > T(0))
+		++np ;
+	else if (x < T(0))
+		++nn ;
 	else
-		if (x<T(0))
-			++nn;
-		else ++nz;
+		++nz ;
 
-	if (y>T(0))
-		++np;
+	if (y > T(0))
+		++np ;
+	else if (y < T(0))
+		++nn ;
 	else
-		if (y<T(0))
-			++nn;
-		else ++nz;
+		++nz ;
 
-	if (z>T(0))
-		++np;
+	if (z > T(0))
+		++np ;
+	else if (z < T(0))
+		++nn ;
 	else
-		if (z<T(0))
-			++nn;
-		else ++nz;
+		++nz ;
 
-	if ((np !=0) && (nn!=0))
-		return NO_INTERSECTION;
+	if ((np != 0) && (nn != 0)) return NO_INTERSECTION ;
 
 	T sum = x + y + z ;
 	T alpha = y / sum ;
@@ -99,15 +96,11 @@ Intersection intersectionRayTriangleOpt(const VEC3& P, const VEC3& Dir, const VE
 	T gamma = T(1) - alpha - beta ;
 	Inter = Ta * alpha + Tb * beta + Tc * gamma ;
 
-
-	return Intersection(FACE_INTERSECTION-nz);
-
+	return Intersection(FACE_INTERSECTION - nz) ;
 }
 
-
-
 template <typename VEC3>
-Intersection intersectionRayTriangleOpt(const VEC3& P, const VEC3& Dir, const VEC3& Ta,  const VEC3& Tb, const VEC3& Tc)
+Intersection intersectionRayTriangleOpt(const VEC3& P, const VEC3& Dir, const VEC3& Ta, const VEC3& Tb, const VEC3& Tc)
 {
 	typedef typename VEC3::DATA_TYPE T ;
 
@@ -119,41 +112,38 @@ Intersection intersectionRayTriangleOpt(const VEC3& P, const VEC3& Dir, const VE
 	T y = tripleProduct(Dir, v, w) ;
 	T z = tripleProduct(Dir, w, u) ;
 
-	unsigned int np=0;
-	unsigned int nn=0;
-	unsigned int nz=0;
+	unsigned int np = 0 ;
+	unsigned int nn = 0 ;
+	unsigned int nz = 0 ;
 
-	if (x>T(0))
-		++np;
+	if (x > T(0))
+		++np ;
+	else if (x < T(0))
+		++nn ;
 	else
-		if (x<T(0))
-			++nn;
-		else ++nz;
+		++nz ;
 
-	if (y>T(0))
-		++np;
+	if (y > T(0))
+		++np ;
+	else if (y < T(0))
+		++nn ;
 	else
-		if (y<T(0))
-			++nn;
-		else ++nz;
+		++nz ;
 
-	if (z>T(0))
-		++np;
+	if (z > T(0))
+		++np ;
+	else if (z < T(0))
+		++nn ;
 	else
-		if (z<T(0))
-			++nn;
-		else ++nz;
+		++nz ;
 
-	if ((np !=0) && (nn!=0))
-		return NO_INTERSECTION;
-
-	return Intersection(FACE_INTERSECTION-nz);
+	if ((np != 0) && (nn != 0)) return NO_INTERSECTION ;
 
+	return Intersection(FACE_INTERSECTION - nz) ;
 }
 
-
 template <typename VEC3>
-Intersection intersectionRayTriangle(const VEC3& P, const VEC3& Dir, const VEC3& Ta,  const VEC3& Tb, const VEC3& Tc, VEC3& Inter)
+Intersection intersectionRayTriangle(const VEC3& P, const VEC3& Dir, const VEC3& Ta, const VEC3& Tb, const VEC3& Tc, VEC3& Inter)
 {
 	typedef typename VEC3::DATA_TYPE T ;
 
@@ -190,7 +180,7 @@ Intersection intersectionRayTriangle(const VEC3& P, const VEC3& Dir, const VEC3&
 
 
 template <typename VEC3>
-Intersection intersectionLineTriangle(const VEC3& P, const VEC3& Dir, const VEC3& Ta,  const VEC3& Tb, const VEC3& Tc, VEC3& Inter)
+Intersection intersectionLineTriangle(const VEC3& P, const VEC3& Dir, const VEC3& Ta, const VEC3& Tb, const VEC3& Tc, VEC3& Inter)
 {
 	typedef typename VEC3::DATA_TYPE T ;
 
@@ -408,48 +398,48 @@ Intersection intersectionSegmentTriangle(const VEC3& PA, const VEC3& PB, const V
 // }
 
 template <typename VEC3, typename PLANE3D>
-Intersection intersectionPlaneRay(const PLANE3D& pl,const VEC3& p1,const VEC3& dir, VEC3& Inter)
+Intersection intersectionPlaneRay(const PLANE3D& pl, const VEC3& p1, const VEC3& dir, VEC3& Inter)
 {
 	typename VEC3::DATA_TYPE denom = pl.normal()*dir;
 
-	if(denom==0)
+	if (denom == 0)
 	{
-		if(pl.distance(p1)==0)
+		if (pl.distance(p1) == 0)
 		{
-				Inter = p1;
-				return FACE_INTERSECTION;
+			Inter = p1 ;
+			return FACE_INTERSECTION ;
 		}
 		else
-			return NO_INTERSECTION;