Merge branch 'master' of CGoGN:CGoGN

include/Algo/Modelisation/planeCutting.h

+/*******************************************************************************
+* CGoGN: Combinatorial and Geometric modeling with Generic N-dimensional Maps  *
+* version 0.1                                                                  *
+* Copyright (C) 2009-2012, IGG Team, LSIIT, University of Strasbourg           *
+*                                                                              *
+* This library is free software; you can redistribute it and/or modify it      *
+* under the terms of the GNU Lesser General Public License as published by the *
+* Free Software Foundation; either version 2.1 of the License, or (at your     *
+* option) any later version.                                                   *
+*                                                                              *
+* This library is distributed in the hope that it will be useful, but WITHOUT  *
+* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or        *
+* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License  *
+* for more details.                                                            *
+*                                                                              *
+* You should have received a copy of the GNU Lesser General Public License     *
+* 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/                                           *
+* Contact information: cgogn@unistra.fr                                        *
+*                                                                              *
+*******************************************************************************/
+
+#ifndef __PLANE_CUTTING_H__
+#define __PLANE_CUTTING_H__
+
+#include <math.h>
+#include <vector>
+#include "Geometry/plane_3d.h"
+#include "Topology/generic/cellmarker.h"
+#include "Topology/generic/traversorCell.h"
+#include "Algo/Modelisation/triangulation.h"
+
+namespace CGoGN
+{
+
+namespace Algo
+{
+
+namespace Surface
+{
+
+namespace Modelisation
+{
+
+/**
+*/
+template <typename PFP>
+void planeCut(typename PFP::MAP& map, VertexAttribute<typename PFP::VEC3>& position, const Geom::Plane3D<typename PFP::REAL>& plane,
+			  CellMarker<FACE>& cmf_over, bool keepTriangles=false, bool with_unsew = true);
+
+
+} // namespace Modelisation
+
+} // namespace Surface
+
+} // namespace Algo
+
+} // namespace CGoGN
+
+#include "Algo/Modelisation/planeCutting.hpp"
+
+#endif

include/Algo/Modelisation/planeCutting.hpp

+namespace CGoGN
+{
+
+namespace Algo
+{
+
+namespace Surface
+{
+
+namespace Modelisation
+{
+
+/*
+template <typename PFP>
+void planeCut(typename PFP::MAP& map, VertexAttribute<typename PFP::VEC3>& position, const Geom::Plane3D<typename PFP::REAL>& plane,
+			  CellMarker<FACE>& cmf_over, bool keepTriangles=false, bool with_unsew = true)
+{
+	typedef typename PFP::REAL REAL;
+
+	//marker for vertices on the plane
+	CellMarker<VERTEX> cmv(map);
+
+	// marker for vertices over the plane
+	CellMarker<VERTEX> cmv_over(map);
+
+
+	TraversorV<typename PFP::MAP> traVert(map);
+	for (Dart d=traVert.begin(); d!=traVert.end();d=traVert.next())
+	{
+		Geom::Orientation3D or1 = plane.orient(position[d]);
+		if (or1 == Geom::ON)
+			cmv.mark(d);
+		if (or1 == Geom::OVER)
+			cmv_over.mark(d);
+	}
+
+
+
+	TraversorE<typename PFP::MAP> traEdg(map);
+	for (Dart d=traEdg.begin(); d!=traEdg.end();d=traEdg.next())
+	{
+		Dart dd = map.phi1(d);
+
+
+		if (!cmv.isMarked(d) && !cmv.isMarked(dd) && (cmv_over.isMarked(d) != cmv_over.isMarked(dd)))
+		{
+			Dart x = map.cutEdge(d);
+			REAL dist1 = plane.distance(position[d]);
+			REAL dist2 = plane.distance(position[dd]);
+
+			if (dist1<0.0)
+				dist1 = -dist1;
+			if (dist2<0.0)			// no abs() to avoid type problem with REAL template
+				dist2 = -dist2;
+
+			position[x] = (position[d]*dist2 + position[dd]*dist1)/(dist1+dist2);
+
+			traEdg.skip(x);
+			traEdg.skip(map.phi_1(x));
+			cmv.mark(x);
+		}
+	}
+
+	Algo::Surface::Modelisation::EarTriangulation<PFP>* triangulator=NULL;
+	if (keepTriangles)	// triangule faces if needed
+	{
+		triangulator = new Algo::Surface::Modelisation::EarTriangulation<PFP>(map);
+	}
+
+	TraversorF<typename PFP::MAP> traFac(map);
+	for (Dart d=traFac.begin(); d!=traFac.end();d=traFac.next())
+	{
+		// turn in the face to search if there are 2 vertices marked as on the plane
+		Traversor2FV<typename PFP::MAP> traV(map,d);
+		Dart e=traV.begin();
+		while ((e!=traV.end())&&(!cmv.isMarked(e)))
+			e=traV.next();
+
+		Dart V1=NIL;
+		if (e!=traV.end())
+			V1 = e;
+
+		e=traV.next();
+		while ((e!=traV.end())&&(!cmv.isMarked(e)))
+			e=traV.next();
+
+		Dart V2=NIL;
+		if (e!=traV.end())
+			V2 = e;
+
+		// is there 2 vertices in the plane (but not consecutive)
+		if ((V1!=NIL) && (V2!=NIL) && (V2!=map.phi1(V1)) && (V1!=map.phi1(V2)))
+		{
+			map.splitFace(V1,V2);
+			if (with_unsew)
+				map.unsewFaces(map.phi_1(V1)); // ne marche pas !
+
+			// ensure to not scan this two new faces
+			traFac.skip(V1);
+			traFac.skip(V2);
+			// mark face of V1 if it is over
+			if (cmv_over.isMarked(map.phi1(V1)))
+				cmf_over.mark(V1);
+			// mark face of V2 if it is over
+			if (cmv_over.isMarked(map.phi1(V2)))
+				cmf_over.mark(V2);
+
+			if (keepTriangles)	// triangule faces if needed
+			{
+				triangulator->trianguleFace(V1);
+				triangulator->trianguleFace(V2);
+			}
+		}
+		else
+		{
+			// find the first vertex not on the plane
+			Dart e = d;
+			while (cmv.isMarked(e))
+				e = map.phi1(e);
+			// face is all on same side than vertex
+			if (cmv_over.isMarked(e))
+				cmf_over.mark(e);
+		}
+	}
+
+	if (triangulator != NULL)
+		delete triangulator;
+}
+*/
+
+template <typename PFP>
+void planeCut(typename PFP::MAP& map, VertexAttribute<typename PFP::VEC3>& position, const Geom::Plane3D<typename PFP::REAL>& plane,
+			  CellMarker<FACE>& cmf_over, bool keepTriangles=false, bool with_unsew = true)
+{
+	typedef typename PFP::REAL REAL;
+
+	//marker for vertices on the plane
+	CellMarker<VERTEX> cmv(map);
+	// marker for vertices over the plane
+	CellMarker<VERTEX> cmv_over(map);
+
+	TraversorE<typename PFP::MAP> traEdg(map);
+	for (Dart d=traEdg.begin(); d!=traEdg.end();d=traEdg.next())
+	{
+		Dart dd = map.phi1(d);
+
+		REAL dist1;
+		REAL dist2;
+
+		Geom::Orientation3D or1 = plane.orient(position[d],dist1);
+		Geom::Orientation3D or2 = plane.orient(position[dd],dist2);
+
+		if (or1 == Geom::ON)
+			cmv.mark(d);
+
+		if (or2 == Geom::ON)
+			cmv.mark(dd);
+
+		if ((or1!=Geom::ON) && (or2!=Geom::ON) && (or1 != or2))
+		{
+			Dart x = map.cutEdge(d);
+
+			if (dist1<0.0)
+				dist1 = -dist1;
+			if (dist2<0.0)			// no abs() to avoid type problem with REAL template
+				dist2 = -dist2;
+
+			position[x] = (position[d]*dist2 + position[dd]*dist1)/(dist1+dist2);
+
+			traEdg.skip(x);
+			traEdg.skip(map.phi_1(x));
+			cmv.mark(x);
+
+			if (or1 == Geom::OVER)
+				cmv_over.mark(d);
+			else
+				cmv_over.mark(dd);
+		}
+		else
+		{
+			if (or1 == Geom::OVER)
+			{
+				cmv_over.mark(d);
+				cmv_over.mark(dd);
+			}
+		}
+	}
+
+	Algo::Surface::Modelisation::EarTriangulation<PFP>* triangulator=NULL;
+	if (keepTriangles)	// triangule faces if needed
+	{
+		triangulator = new Algo::Surface::Modelisation::EarTriangulation<PFP>(map);
+	}
+
+	TraversorF<typename PFP::MAP> traFac(map);
+	for (Dart d=traFac.begin(); d!=traFac.end();d=traFac.next())
+	{
+		// turn in the face to search if there are 2 vertices marked as on the plane
+		Traversor2FV<typename PFP::MAP> traV(map,d);
+		Dart e=traV.begin();
+		while ((e!=traV.end())&&(!cmv.isMarked(e)))
+			e=traV.next();
+
+		Dart V1=NIL;
+		if (e!=traV.end())
+			V1 = e;
+
+		e=traV.next();
+		while ((e!=traV.end())&&(!cmv.isMarked(e)))
+			e=traV.next();
+
+		Dart V2=NIL;
+		if (e!=traV.end())
+			V2 = e;
+
+		// is there 2 vertices in the plane (but not consecutive)
+		if ((V1!=NIL) && (V2!=NIL) && (V2!=map.phi1(V1)) && (V1!=map.phi1(V2)))
+		{
+			map.splitFace(V1,V2);
+			if (with_unsew)
+				map.unsewFaces(map.phi_1(V1)); // ne marche pas !
+
+			// ensure to not scan this two new faces
+			traFac.skip(V1);
+			traFac.skip(V2);
+			// mark face of V1 if it is over
+			if (cmv_over.isMarked(map.phi1(V1)))
+				cmf_over.mark(V1);
+			// mark face of V2 if it is over
+			if (cmv_over.isMarked(map.phi1(V2)))
+				cmf_over.mark(V2);
+
+			if (keepTriangles)	// triangule faces if needed
+			{
+				triangulator->trianguleFace(V1);
+				triangulator->trianguleFace(V2);
+			}
+		}
+		else
+		{
+			// find the first vertex not on the plane
+			Dart e = d;
+			while (cmv.isMarked(e))
+				e = map.phi1(e);
+			// face is all on same side than vertex
+			if (cmv_over.isMarked(e))
+				cmf_over.mark(e);
+		}
+	}
+
+	if (triangulator != NULL)
+		delete triangulator;
+}
+
+
+} // namespace Modelisation
+
+} // namespace Surface
+
+} // namespace Algo
+
+} // namespace CGoGN

include/Geometry/distances.h

@@ -89,6 +89,11 @@ typename VEC3::DATA_TYPE squaredDistanceLine2Point(const VEC3& A, const VEC3& B,
 
 /**
 * compute squared distance from line to line
+* @param A point of first line
+* @param AB vector of first line
+* @param AB2 AB*AB (for optimization if call several times with AB
+* @param P first point of second line
+* @param Q second point of second line
 * @return the squared distance
 */
 template <typename VEC3>
@@ -96,6 +101,11 @@ typename VEC3::DATA_TYPE squaredDistanceLine2Line(const VEC3& A, const VEC3& AB,
 
 /**
 * compute squared distance from line to segment
+* @param A point of line
+* @param AB vector of line
+* @param AB2 AB*AB (for optimization if call several times with AB
+* @param P first point of segment
+* @param Q second point of segment
 * @return the squared distance
 */
 template <typename VEC3>
@@ -103,6 +113,9 @@ typename VEC3::DATA_TYPE squaredDistanceLine2Seg(const VEC3& A, const VEC3& AB,
 
 /**
 * compute squared distance from segment to point
+* @param A point of segment
+* @param AB vector of segment
+* @param AB2 AB*AB (for optimization if call several times with AB
 * @return the squared distance
 */
 template <typename VEC3>

include/Geometry/intersection.h

@@ -166,6 +166,10 @@ Intersection intersection2DSegmentSegment(const VEC3& PA, const VEC3& PB, const
 template <typename VEC3>
 Intersection intersectionSegmentPlan(const VEC3& PA, const VEC3& PB, const VEC3& PlaneP, const VEC3& NormP); //, VEC3& Inter) ;
 
+template <typename VEC3>
+Intersection intersectionSegmentPlan(const VEC3& PA, const VEC3& PB, const VEC3& PlaneP, const VEC3& NormP, VEC3& Inter) ;
+
+
 template <typename VEC3>
 Intersection intersectionTrianglePlan(const VEC3& Ta, const VEC3& Tb, const VEC3& Tc, const VEC3& PlaneP, const VEC3& NormP);//, VEC3& Inter) ;
 
@@ -177,6 +181,20 @@ template <typename VEC3>
 Intersection intersectionTriangleHalfPlan(const VEC3& Ta, const VEC3& Tb, const VEC3& Tc,
 		const VEC3& P, const VEC3& DirP, const VEC3& OrientP); //, VEC3& Inter) ;
 
+
+/**
+* compute intersection between line and segment
+* @param A point of line
+* @param AB vector of line
+* @param AB2 AB*AB (for optimization if call several times with AB
+* @param P first point of segment
+* @param Q second point of segment
+* @param inter the computed intersection
+* @return intersection of not
+*/
+template <typename VEC3>
+bool interLineSeg(const VEC3& A, const VEC3& AB, typename VEC3::DATA_TYPE AB2,
+				  const VEC3& P, const VEC3& Q, VEC3& inter);
 } // namespace Geom
 
 } // namespace CGoGN

include/Geometry/intersection.hpp

@@ -22,6 +22,8 @@
  *                                                                              *
  *******************************************************************************/
 
+#include "Geometry/distances.h"
+
 namespace CGoGN
 {
 
@@ -487,19 +489,61 @@ Intersection intersection2DSegmentSegment(const VEC3& PA, const VEC3& PB, const
 template <typename VEC3>
 Intersection intersectionSegmentPlan(const VEC3& PA, const VEC3& PB, const VEC3& PlaneP, const VEC3& NormP)//, VEC3& Inter)
 {
-	typename VEC3::DATA_TYPE panp = NormP * PA;
-	typename VEC3::DATA_TYPE pbnp = NormP * PB;
+//	typename VEC3::DATA_TYPE panp = NormP * PA;
+//	typename VEC3::DATA_TYPE pbnp = NormP * PB;
+
+//	if(panp == 0 || pbnp == 0)
+//		return VERTEX_INTERSECTION;
+//	else if((panp < 0 && pbnp > 0) || (panp > 0 && pbnp < 0))
+//		return EDGE_INTERSECTION;
+//	else
+//		return NO_INTERSECTION;
+
+#define EPSILON 1e-12
 
-	if(panp == 0 || pbnp == 0)
+	typename VEC3::DATA_TYPE panp = NormP * (PA-PlaneP);
+	typename VEC3::DATA_TYPE pbnp = NormP * (PB-PlaneP);
+
+	if(abs(panp) < EPSILON || abs(pbnp) < EPSILON)
 		return VERTEX_INTERSECTION;
-	else if((panp < 0 && pbnp > 0) || (panp > 0 && pbnp < 0))
-		return EDGE_INTERSECTION;
+//	else if((panp < 0 && pbnp > 0) || (panp > 0 && pbnp < 0))
+	else if (panp*pbnp < 0)
+			return EDGE_INTERSECTION;
 	else
 		return NO_INTERSECTION;
+#undef EPSILON
+}
+
+template <typename VEC3>
+Intersection intersectionSegmentPlan(const VEC3& PA, const VEC3& PB, const VEC3& PlaneP, const VEC3& NormP, VEC3& Inter)
+{
+#define EPSILON 1e-12
+
+	typename VEC3::DATA_TYPE panp = NormP * (PA-PlaneP);
+	typename VEC3::DATA_TYPE pbnp = NormP * (PB-PlaneP);
 
+	if(abs(panp) < EPSILON)
+	{
+		Inter = PA;
+		return VERTEX_INTERSECTION;
+	}
+	else if(abs(pbnp) < EPSILON)
+	{
+		Inter = PB;
+		return VERTEX_INTERSECTION;
+	}
+	else if (panp*pbnp < 0)
+	{
+		Inter = (abs(panp)*PB + abs(pbnp)*PA)/(abs(panp)+abs(pbnp)) ;
+		return EDGE_INTERSECTION;
+	}
+	else
+		return NO_INTERSECTION;
+#undef EPSILON
 }
 
 
+
 template <typename VEC3>
 Intersection intersectionTrianglePlan(const VEC3& Ta, const VEC3& Tb, const VEC3& Tc, const VEC3& PlaneP, const VEC3& NormP) //, VEC3& Inter) ;
 {
@@ -566,6 +610,35 @@ Intersection intersectionTriangleHalfPlan(const VEC3& Ta, const VEC3& Tb, const
 	}
 }
 
+template <typename VEC3>
+bool interLineSeg(const VEC3& A, const VEC3& AB, typename VEC3::DATA_TYPE AB2,
+				  const VEC3& P, const VEC3& Q, VEC3& inter)
+{
+#define EPSILON (1.0e-5)
+	typedef typename VEC3::DATA_TYPE T ;
+
+	T dist = Geom::distancePoint2TrianglePlane(AB-A,A,P,Q);
+
+//	std::cout << "dist "<<  dist << std::endl;
+
+	if (dist>EPSILON)
+		return false;
+
+	VEC3 AP = P - A ;
+	VEC3 PQ = Q - P ;
+	T X = AB * PQ ;
+	T beta = ( AB2 * (AP*PQ) - X * (AP*AB) ) / ( X*X - AB2 * PQ.norm2() ) ;
+
+//	std::cout << "beta "<<  beta << std::endl;
+
+	if ((beta<0.0) || (beta>1.0))
+		return false;
+
+	inter = beta*Q +(1.0-beta)*P;
+	return true;
+#undef EPSILON
+}
+
 }
 
 }

include/Geometry/plane_3d.h

@@ -92,6 +92,11 @@ class Plane3D
 		// return on/over/under according to the side of the plane where point p is
 		Orientation3D orient(const Vector<3,T>& p) const;
 
+		// return on/over/under according to the side of the plane where point p is
+		// return the distance to the point in dist reference
+		Orientation3D orient(const Vector<3,T>& p, T& dist) const;
+
+
 		/**********************************************/
 		/*             STREAM OPERATORS               */
 		/**********************************************/

include/Geometry/plane_3d.hpp

@@ -139,6 +139,23 @@ Orientation3D Plane3D<T>::orient(const Vector<3,T>& p) const
 #undef PRECISION
 }
 
+
+template <typename T>
+Orientation3D Plane3D<T>::orient(const Vector<3,T>& p, T& dist) const
+{
+#define PRECISION 1e-6
+	dist = distance(p) ;
+
+	if(dist < -PRECISION)
+		return UNDER ;
+	if(dist > PRECISION)
+		return OVER ;
+	return ON ;
+#undef PRECISION
+}
+
+
+
 /**********************************************/
 /*             STREAM OPERATORS               */
 /**********************************************/

include/Utils/GLSLShader.h

@@ -105,7 +105,7 @@ protected:
 	/**
 	 * handle of program
 	 */
-	CGoGNGLhandleARB m_program_object;
+	CGoGNGLint m_program_object;
 
 	CGoGNGLint m_uniMat_Proj;
 	CGoGNGLint m_uniMat_Model;
@@ -178,7 +178,7 @@ protected:
 	 * @param obj what log do you want ?
 	 * @return the log
 	 */
-	char* getInfoLog( GLhandleARB obj );
+	char* getInfoLog( GLuint obj );
 
 public:
 	/**
@@ -304,7 +304,7 @@ public:
 	/**
 	 * get handler of program for external use og gl functions
 	 */
-	GLhandleARB program_handler() { return *m_program_object;}
+	GLuint program_handler() { return *m_program_object;}
 
 	/**
 	 * check shader validity width official GLSL syntax

include/Utils/gl_def.h

@@ -61,7 +61,7 @@ public:
 
 typedef FalsePtr<GLint> CGoGNGLint;
 typedef FalsePtr<GLuint> CGoGNGLuint;
-typedef FalsePtr<GLhandleARB> CGoGNGLhandleARB;
+typedef FalsePtr<GLuint> CGoGNGLhandleARB;
 typedef FalsePtr<GLenum> CGoGNGLenum;
 typedef FalsePtr<GLenum*> CGoGNGLenumTable;
 

src/Topology/map/embeddedMap2.cpp

@@ -440,6 +440,9 @@ void EmbeddedMap2::unsewFaces(Dart d, bool withBoundary)
 
 	if (isOrbitEmbedded<VERTEX>())
 	{
+		copyDartEmbedding<VERTEX>(phi2(e), d) ;
+		copyDartEmbedding<VERTEX>(phi2(d), e) ;
+
 		Dart ee = phi1(e) ;
 		if(!sameVertex(d, ee))
 		{

src/Utils/Qt/qtQGLV.cpp

@@ -305,13 +305,13 @@ void SimpleQGLV::keyPressEvent(QKeyEvent *e)
     else
         QWidget::keyPressEvent(e);
 
-	m_qglWidget->keyPressEvent(e); // ?
+//	m_qglWidget->keyPressEvent(e); // ?
 }
 
 void SimpleQGLV::keyReleaseEvent(QKeyEvent *e)
 {
 	QWidget::keyReleaseEvent(e);
-	m_qglWidget->keyReleaseEvent(e);
+//	m_qglWidget->keyReleaseEvent(e);
 }
 
 void SimpleQGLV::glMousePosition(int& x, int& y)