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

Apps/Examples/Release/CMakeLists.txt

@@ -17,7 +17,6 @@ include_directories(
 	${CGoGN_ROOT_DIR}/ThirdParty/OpenCTM
 	${CGoGN_ROOT_DIR}/ThirdParty/Assimp/include
 	${CGoGN_ROOT_DIR}/ThirdParty/glm
-
 )
 
 # define libs path

Apps/Examples/Tests/CMakeLists.txt

@@ -9,6 +9,15 @@ SET(CMAKE_BUILD_TYPE Debug)
 # define includes path
 include_directories(
 	${CGoGN_ROOT_DIR}/include
+	/usr/include/libxml2/
+	${CGoGN_ROOT_DIR}/include
+	${CGoGN_ROOT_DIR}/ThirdParty/Numerical
+	${CGoGN_ROOT_DIR}/ThirdParty/Numerical/UFconfig
+	${CGoGN_ROOT_DIR}/ThirdParty/gzstream
+	${CGoGN_ROOT_DIR}/ThirdParty/Zinri
+	${CGoGN_ROOT_DIR}/ThirdParty/AntTweakBar/include
+	${CGoGN_ROOT_DIR}/ThirdParty/OpenCTM
+	${CGoGN_ROOT_DIR}/ThirdParty/Assimp/include
 )
 
 # define libs path
@@ -26,3 +35,7 @@ target_link_libraries( Geom_orientationD
 add_executable( Geom_inclusionD ./Geom_inclusion.cpp)
 target_link_libraries( Geom_inclusionD
 	${COMMON_LIBS} )
+
+add_executable( Algo_Geom_inclusionD ./Algo_Geom_inclusion.cpp)
+target_link_libraries( Algo_Geom_inclusionD
+	containerD topologyD utilsD algoD ${COMMON_LIBS} )

Apps/Examples/Tests/Geom_inclusion.cpp

@@ -39,6 +39,8 @@ int main()
 	VEC3 tc(0,5,0);
 	VEC3 pIn1(1,1,0);
 	VEC3 pIn2(2,0,0);
+	VEC3 pIn3 = (tb+tc)/2.0f;
+	VEC3 pIn4 = (tc+ta)/2.0f;
 	VEC3 pOut(6,0,0);
 
 	if(Geom::isPointInTriangle<VEC3>(pIn1,ta,tb,tc)!=Geom::FACE_INCLUSION)
@@ -51,6 +53,16 @@ int main()
 		std::cout << "ERROR : isPointInTriangle : edge inclusion" << std::endl;
 	}
 
+	if(Geom::isPointInTriangle<VEC3>(pIn3,ta,tb,tc)!=Geom::EDGE_INCLUSION)
+	{
+		std::cout << "ERROR : isPointInTriangle : edge inclusion" << std::endl;
+	}
+
+	if(Geom::isPointInTriangle<VEC3>(pIn4,ta,tb,tc)!=Geom::EDGE_INCLUSION)
+	{
+		std::cout << "ERROR : isPointInTriangle : edge inclusion" << std::endl;
+	}
+
 	if(Geom::isPointInTriangle<VEC3>(ta,ta,tb,tc)!=Geom::VERTEX_INCLUSION)
 	{
 		std::cout << "ERROR : isPointInTriangle : vertex inclusion" << std::endl;

include/Algo/Geometry/area.h

@@ -46,11 +46,11 @@ typename PFP::REAL totalArea(typename PFP::MAP& map, const typename PFP::TVEC3&
 template <typename PFP>
 void computeAreaFaces(typename PFP::MAP& map, const typename PFP::TVEC3& position, typename PFP::TREAL& face_area, const FunctorSelect& select = SelectorTrue()) ;
 
-}
+} // namespace Geometry
 
-}
+} // namespace Algo
 
-}
+} // namespace CGoGN
 
 #include "Algo/Geometry/area.hpp"
 

include/Algo/Geometry/basic.h

@@ -40,7 +40,7 @@ namespace Geometry
  * vectorOutOfDart return a dart from the position of vertex attribute of d to the position of vertex attribute of phi1(d)
  */
 template <typename PFP>
-typename PFP::VEC3 vectorOutOfDart(typename PFP::MAP& map, Dart d, const typename PFP::TVEC3& position)
+inline typename PFP::VEC3 vectorOutOfDart(typename PFP::MAP& map, Dart d, const typename PFP::TVEC3& position)
 {
 	typename PFP::VEC3 vec = position[map.phi1(d)] ;
 	vec -= position[d] ;
@@ -48,11 +48,17 @@ typename PFP::VEC3 vectorOutOfDart(typename PFP::MAP& map, Dart d, const typenam
 }
 
 template <typename PFP>
-float angle(typename PFP::MAP& map, Dart d1, Dart d2, const typename PFP::TVEC3& position)
+inline typename PFP::REAL edgeLength(typename PFP::MAP& map, Dart d, const typename PFP::TVEC3& position)
 {
-	typename PFP::VEC3 v1, v2 ;
-	vectorOutOfDart<PFP>(map, d1, position, v1) ;
-	vectorOutOfDart<PFP>(map, d2, position, v2) ;
+	typename PFP::VEC3 v = vectorOutOfDart<PFP>(map, d, position) ;
+	return v.norm() ;
+}
+
+template <typename PFP>
+inline float angle(typename PFP::MAP& map, Dart d1, Dart d2, const typename PFP::TVEC3& position)
+{
+	typename PFP::VEC3 v1 = vectorOutOfDart<PFP>(map, d1, position) ;
+	typename PFP::VEC3 v2 = vectorOutOfDart<PFP>(map, d2, position) ;
 	return Geom::angle(v1, v2) ;
 }
 

include/Algo/Geometry/centroid.hpp

@@ -119,9 +119,9 @@ EMB vertexNeighborhoodCentroidGen(typename PFP::MAP& map, Dart d, const EMBV& at
 	Dart it = d ;
 	do
 	{
-		center += attributs[map.phi2(it)];
+		center += attributs[map.phi1(it)];
 		++count ;
-		it = map.phi1(map.phi2(it)) ;
+		it = map.alpha1(it) ;
 	} while(it != d) ;
 	center /= count ;
 	return center ;

include/Algo/Geometry/feature.h

+/*******************************************************************************
+* CGoGN: Combinatorial and Geometric modeling with Generic N-dimensional Maps  *
+* version 0.1                                                                  *
+* Copyright (C) 2009, 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: https://iggservis.u-strasbg.fr/CGoGN/                              *
+* Contact information: cgogn@unistra.fr                                        *
+*                                                                              *
+*******************************************************************************/
+
+#ifndef __ALGO_GEOMETRY_FEATURE_H__
+#define __ALGO_GEOMETRY_FEATURE_H__
+
+namespace CGoGN
+{
+
+namespace Algo
+{
+
+namespace Geometry
+{
+
+template <typename PFP>
+void featureEdgeDetection(typename PFP::MAP& map, const typename PFP::TVEC3& position, DartMarker& feature) ;
+
+} // namespace Geometry
+
+} // namespace Algo
+
+} // namespace CGoGN
+
+#include "Algo/Geometry/feature.hpp"
+
+#endif

include/Algo/Geometry/feature.hpp

+/*******************************************************************************
+* CGoGN: Combinatorial and Geometric modeling with Generic N-dimensional Maps  *
+* version 0.1                                                                  *
+* Copyright (C) 2009, 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: https://iggservis.u-strasbg.fr/CGoGN/                              *
+* Contact information: cgogn@unistra.fr                                        *
+*                                                                              *
+*******************************************************************************/
+
+#include "Geometry/basic.h"
+#include "Algo/Geometry/normal.h"
+
+namespace CGoGN
+{
+
+namespace Algo
+{
+
+namespace Geometry
+{
+
+template <typename PFP>
+void featureEdgeDetection(typename PFP::MAP& map, typename PFP::TVEC3& position, DartMarker& feature)
+{
+	typedef typename PFP::VEC3 VEC3 ;
+	typedef typename PFP::REAL REAL ;
+
+	feature.unmarkAll() ;
+
+	AttributeHandler<VEC3> fNormal = map.template addAttribute<VEC3>(FACE_ORBIT, "fNormal") ;
+	Algo::Geometry::computeNormalFaces<PFP>(map, position, fNormal) ;
+
+	DartMarker m(map) ;
+	for(Dart d = map.begin(); d != map.end(); map.next(d))
+	{
+		if(!m.isMarked(d))
+		{
+			m.markOrbit(EDGE_ORBIT, d) ;
+			if(Geom::angle(fNormal[d], fNormal[map.phi2(d)]) > M_PI / REAL(6))
+				feature.markOrbit(EDGE_ORBIT, d) ;
+		}
+	}
+
+	map.template removeAttribute<VEC3>(fNormal) ;
+}
+
+} // namespace Geometry
+
+} // namespace Algo
+
+} // namespace CGoGN

include/Algo/Geometry/inclusion.h

@@ -73,7 +73,7 @@ bool isPointInConvexVolume(typename PFP::MAP& map, Dart d, const typename PFP::T
  * @param the point
  */
 template <typename PFP>
-bool isPointInConvexFace2D(typename PFP::MAP& map, Dart d, const typename PFP::TVEC3& positions, const typename PFP::VEC3& point, bool CCW );
+bool isPointInConvexFace2D(typename PFP::MAP& map, Dart d, const typename PFP::TVEC3& positions, const typename PFP::VEC3& point, bool CCW=true);
 
 /**
  * test if a point is inside a face

include/Algo/Geometry/localFrame.h

@@ -44,8 +44,8 @@ namespace Geometry
 template <typename PFP>
 void vertexLocalFrame(typename PFP::MAP& map, Dart d, const typename PFP::TVEC3& position, typename PFP::VEC3& X, typename PFP::VEC3& Y, typename PFP::VEC3& Z)
 {
-	Z = vertexNormal<PFP>(map,d,position) ;
-	X = vectorOutOfDart<PFP>(map,d,position) ;
+	Z = vertexNormal<PFP>(map, d, position) ;
+	X = vectorOutOfDart<PFP>(map, d, position) ;
 	Y = Z ^ X ;
 	Y.normalize() ;
 	X = Y ^ Z ;
@@ -56,7 +56,7 @@ template <typename PFP>
 typename PFP::MATRIX33 vertexLocalFrame(typename PFP::MAP& map, Dart d, const typename PFP::TVEC3& position)
 {
 	typename PFP::VEC3 X, Y, Z ;
-	vertexLocalFrame<PFP>(map,d,position,X,Y,Z) ;
+	vertexLocalFrame<PFP>(map, d, position, X, Y, Z) ;
 	typename PFP::MATRIX33 frame ;
 	frame(0,0) = X[0] ;	frame(0,1) = X[1] ;	frame(0,2) = X[2] ;
 	frame(1,0) = Y[0] ;	frame(1,1) = Y[1] ;	frame(1,2) = Y[2] ;
@@ -70,7 +70,7 @@ template <typename PFP>
 void vertexLocalFrame(typename PFP::MAP& map, Dart d, const typename PFP::TVEC3& position, typename PFP::VEC3& normal, typename PFP::VEC3& X, typename PFP::VEC3& Y, typename PFP::VEC3& Z)
 {
 	Z = normal ;
-	X = vectorOutOfDart<PFP>(map,d,position) ;
+	X = vectorOutOfDart<PFP>(map, d, position) ;
 	Y = Z ^ X ;
 	Y.normalize() ;
 	X = Y ^ Z ;
@@ -81,7 +81,7 @@ template <typename PFP>
 typename PFP::MATRIX33 vertexLocalFrame(typename PFP::MAP& map, Dart d, const typename PFP::TVEC3& position, typename PFP::VEC3& normal)
 {
 	typename PFP::VEC3 X, Y, Z ;
-	vertexLocalFrame<PFP>(map,d,position,normal,X,Y,Z) ;
+	vertexLocalFrame<PFP>(map, d, position, normal, X, Y, Z) ;
 	typename PFP::MATRIX33 frame ;
 	frame(0,0) = X[0] ;	frame(0,1) = X[1] ;	frame(0,2) = X[2] ;
 	frame(1,0) = Y[0] ;	frame(1,1) = Y[1] ;	frame(1,2) = Y[2] ;

include/Algo/Geometry/normal.hpp

@@ -40,9 +40,9 @@ typename PFP::VEC3 triangleNormal(typename PFP::MAP& map, Dart d, const typename
 {
 	typedef typename PFP::VEC3 VEC3 ;
 
-	VEC3 p1 = position[d];
-	VEC3 p2 = position[map.phi1(d)];
-	VEC3 p3 = position[map.phi_1(d)];
+	const VEC3& p1 = position[d];
+	const VEC3& p2 = position[map.phi1(d)];
+	const VEC3& p3 = position[map.phi_1(d)];
 
 	VEC3 N = Geom::triangleNormal(p1, p2, p3) ;
 	N.normalize() ;
@@ -64,8 +64,8 @@ typename PFP::VEC3 faceNormal(typename PFP::MAP& map, Dart d, const typename PFP
 		do
 		{
 			VEC3 n = triangleNormal<PFP>(map, it, position) ;
-			//if(!std::isnan(n[0]))
-			if (n[0] == n[0])
+			//if(!std::isnan(n[0]) && !std::isnan(n[1]) && !std::isnan(n[2]))
+			if (n[0] == n[0] && n[1] == n[1] && n[2] == n[2])
 				N += n ;
 			it = map.phi1(it) ;
 		} while (it != d) ;
@@ -84,10 +84,13 @@ typename PFP::VEC3 vertexNormal(typename PFP::MAP& map, Dart d, const typename P
 	do
 	{
 		VEC3 n = faceNormal<PFP>(map, it, position) ;
-		VEC3 v1 = vectorOutOfDart<PFP>(map, it, position) ;
-		VEC3 v2 = vectorOutOfDart<PFP>(map, map.phi_1(it), position) ;
-		n *= convexFaceArea<PFP>(map,it,position) / (v1.norm2() * v2.norm2()) ;
-		N += n ;
+		if(!n.hasNan())
+		{
+			VEC3 v1 = vectorOutOfDart<PFP>(map, it, position) ;
+			VEC3 v2 = vectorOutOfDart<PFP>(map, map.phi_1(it), position) ;
+			n *= convexFaceArea<PFP>(map, it, position) / (v1.norm2() * v2.norm2()) ;
+			N += n ;
+		}
 		it = map.phi1(map.phi2(it)) ;
 	} while (it != d) ;
 	N.normalize() ;

include/Algo/ImplicitHierarchicalMesh/ihm.hpp

@@ -347,11 +347,10 @@ T& AttributeHandler_IHM<T>::operator[](Dart d)
 	unsigned int nbSteps = m->m_curLevel - m->vertexInsertionLevel(d) ;
 	unsigned int index = m->getEmbedding(d, orbit) ;
 
-	if(index == EMBNULL && nbSteps == 0)
+	if(index == EMBNULL)
 	{
 		index = m->embedNewCell(orbit, d) ;
 		m->m_nextLevelCell[orbit]->operator[](index) = EMBNULL ;
-		return this->m_attrib->operator[](index);
 	}
 
 	AttributeContainer& cont = m->getAttributeContainer(orbit) ;

include/Algo/ImplicitHierarchicalMesh/subdivision.hpp

@@ -162,13 +162,39 @@ void coarsenFace(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& position)
 	assert(map.faceIsSubdividedOnce(d) || !"Trying to coarsen a non-subdivided face or a more than once subdivided face") ;
 
 	unsigned int cur = map.getCurrentLevel() ;
-	map.setCurrentLevel(cur + 1) ;
-	Dart cv = map.phi1(map.phi1(d)) ;
-	map.setCurrentLevel(map.getMaxLevel()) ;
-	map.deleteVertex(cv) ;
-	map.setCurrentLevel(cur) ;
+
+	unsigned int degree = 0 ;
 	Dart fit = d ;
 	do
+	{
+		++degree ;
+		fit = map.phi1(fit) ;
+	} while(fit != d) ;
+
+	if(degree == 3)
+	{
+		fit = d ;
+		do
+		{
+			map.setCurrentLevel(cur + 1) ;
+			Dart innerEdge = map.phi1(fit) ;
+			map.setCurrentLevel(map.getMaxLevel()) ;
+			map.mergeFaces(innerEdge) ;
+			map.setCurrentLevel(cur) ;
+			fit = map.phi1(fit) ;
+		} while(fit != d) ;
+	}
+	else
+	{
+		map.setCurrentLevel(cur + 1) ;
+		Dart centralV = map.phi1(map.phi1(d)) ;
+		map.setCurrentLevel(map.getMaxLevel()) ;
+		map.deleteVertex(centralV) ;
+		map.setCurrentLevel(cur) ;
+	}
+
+	fit = d ;
+	do
 	{
 		if(map.edgeCanBeCoarsened(fit))
 			coarsenEdge<PFP>(map, fit, position) ;

include/Algo/Modelisation/extrusion.h

@@ -194,9 +194,11 @@ template<typename PFP>
 Dart extrudeFace(typename PFP::MAP& the_map, typename PFP::TVEC3& positions,
 	Dart d,	float dist);
 
-}//end namespace
-}//end namespace
-}//end namespace
+} // namespace Modelisation
+
+} // namespace Algo
+
+} // namespace CGoGN
 
 #include "Algo/Modelisation/extrusion.hpp"
 

include/Algo/Modelisation/subdivision.h

@@ -54,6 +54,16 @@ void trianguleFaces(typename PFP::MAP& map, EMBV& attributs, const FunctorSelect
 template <typename PFP>
 void trianguleFaces(typename PFP::MAP& map, typename PFP::TVEC3& position, const FunctorSelect& selected = SelectorTrue()) ;
 
+/**
+ * Triangule all the faces of the mesh
+ * positions for new vertices are given as face attribute
+ */
+template <typename PFP>
+void trianguleFaces(
+	typename PFP::MAP& map,
+	typename PFP::TVEC3& position, typename PFP::TVEC3& positionF,
+	const FunctorSelect& selected = SelectorTrue()) ;
+
 /**
 * Quadrangule a face with central vertex
 * warning: edges are considered here as already cut !!

include/Algo/Modelisation/subdivision.hpp

@@ -84,6 +84,33 @@ void trianguleFaces(typename PFP::MAP& map, typename PFP::TVEC3& position, const
 	trianguleFaces<PFP, typename PFP::TVEC3, typename PFP::VEC3>(map, position, selected) ;
 }
 
+template <typename PFP>
+void trianguleFaces(
+		typename PFP::MAP& map,
+		typename PFP::TVEC3& position, typename PFP::TVEC3& positionF,
+		const FunctorSelect& selected)
+{
+	assert(position.getOrbit() == VERTEX_ORBIT) ;
+	assert(positionF.getOrbit() == FACE_ORBIT) ;
+
+	DartMarker m(map) ;
+	for (Dart d = map.begin(); d != map.end(); map.next(d))
+	{
+		if (selected(d) && !m.isMarked(d))
+		{
+			typename PFP::VEC3 p = positionF[d] ;
+			Dart cd = trianguleFace<PFP>(map, d) ;	// triangule the face
+			position[cd] = p ;						// affect the data to the central vertex
+			Dart fit = cd ;
+			do
+			{
+				m.markOrbit(FACE_ORBIT, fit);
+				fit = map.alpha1(fit);
+			} while(fit != cd);
+		}
+	}
+}
+
 template <typename PFP>
 Dart quadranguleFace(typename PFP::MAP& map, Dart d)
 {
@@ -428,6 +455,74 @@ void LoopSubdivision(typename PFP::MAP& map, EMBV& attributs, const FunctorSelec
 	}
 }
 
+template <typename PFP, typename EMBV, typename EMB>
+void TwoNPlusOneSubdivision(typename PFP::MAP& map, EMBV& attributs, const FunctorSelect& selected)
+{
+	AutoAttributeHandler< Dart > tablePred(map,EDGE_ORBIT);
+	CellMarker m0(map, EDGE_CELL);
+
+	std::vector<Dart> dOrig;
+	std::vector<EMB> cOrig;
+
+	//first pass cut edge
+	for (Dart d = map.begin(); d != map.end(); map.next(d))
+	{
+		if(selected(d))
+		{
+			if(!m0.isMarked(d)) {
+				dOrig.push_back(d);
+
+				Dart dd = d;
+				do {
+					if(!m0.isMarked(dd)) {
+						EMB e1 = attributs[dd];
+						EMB e2 = attributs[map.phi1(dd)];
+						map.cutEdge(dd);
+						attributs[map.phi1(dd)] = e1*2.0f/3.0f+e2/3.0f;
+						map.cutEdge(map.phi1(dd));
+						attributs[map.phi1(map.phi1(dd))] = e2*2.0f/3.0f+e1/3.0f;
+						m0.mark(dd);
+						m0.mark(map.phi1(dd));
+						m0.mark(map.template phi<11>(dd));
+					}
+					dd = map.template phi<111>(dd);
+				} while(dd!=d);
+
+			}
+		}
+	}
+
+	std::cout << "nb orig : " << dOrig.size() << std::endl;
+
+	DartMarkerNoUnmark mCorner(map);
+//	//second pass create corner face
+	for (std::vector<Dart>::iterator it = dOrig.begin(); it != dOrig.end(); ++it)
+	{
+		EMB c = Algo::Geometry::faceCentroid<PFP>(map,*it,attributs);
+		Dart dd = *it;
+		do
+		{
+			map.splitFace(map.phi1(dd),map.phi_1(dd));
+			map.cutEdge(map.phi1(dd));
+			mCorner.mark(map.phi2(map.phi1(dd)));
+			attributs[map.template phi<11>(dd)] = c*2.0/3.0f + attributs[dd]/3.0f;
+			dd = map.phi1(map.phi1(map.phi1(map.phi2(map.phi1(dd)))));
+		} while(!mCorner.isMarked(dd));
+	}
+
+	//third pass create center face
+	for (std::vector<Dart>::iterator it = dOrig.begin(); it != dOrig.end(); ++it)
+	{
+		Dart dd = map.phi2(map.phi1(*it));
+		do {
+			mCorner.unmark(dd);
+			Dart dNext = map.phi1(map.phi1(map.phi1(dd)));
+			map.splitFace(dd,dNext);
+			dd = dNext;
+		} while(mCorner.isMarked(dd));
+	}
+}
+
 template <typename PFP>
 void LoopSubdivision(typename PFP::MAP& map, typename PFP::TVEC3& position, const FunctorSelect& selected)
 {
@@ -478,11 +573,73 @@ void computeDual(typename PFP::MAP& map, const FunctorSelect& selected)
 	reverseOrientation<PFP>(map) ;
 }
 
+inline double sqrt3_K(unsigned int n)
+{
+	switch(n)
+	{
+		case 1: return 0.333333 ;
+		case 2: return 0.555556 ;
+		case 3: return 0.5 ;
+		case 4: return 0.444444 ;
+		case 5: return 0.410109 ;
+		case 6: return 0.388889 ;
+		case 7: return 0.375168 ;
+		case 8: return 0.365877 ;
+		case 9: return 0.359328 ;
+		case 10: return 0.354554 ;
+		case 11: return 0.350972 ;
+		case 12: return 0.348219 ;
+		default:
+			double t = cos((2.0*M_PI)/double(n)) ;
+			return (4.0 - t) / 9.0 ;
+	}
+}
+
 template <typename PFP>
 void Sqrt3Subdivision(typename PFP::MAP& map, typename PFP::TVEC3& position, const FunctorSelect& selected)
 {
+	typedef typename PFP::VEC3 VEC3 ;
+	typedef typename PFP::REAL REAL ;
+
+	AttributeHandler<VEC3> positionF = map.template getAttribute<VEC3>(FACE_ORBIT, "position") ;
+	if(!positionF.isValid())
+		positionF = map.template addAttribute<VEC3>(FACE_ORBIT, "position") ;
+	Algo::Geometry::computeCentroidFaces<PFP>(map, position, positionF) ;
+
 	computeDual<PFP>(map, selected);
-	trianguleFaces<PFP>(map, position, selected);
+
+	AttributeHandler<VEC3> tmp = position ;
+	position = positionF ;
+	positionF = tmp ;
+
+	CellMarker m(map, VERTEX_ORBIT) ;
+	m.markAll() ;
+
+	trianguleFaces<PFP>(map, position, positionF, selected);
+
+	for(Dart d = map.begin(); d != map.end(); map.next(d))
+	{
+		if(!m.isMarked(d))
+		{
+			m.mark(d) ;
+			VEC3 P = position[d] ;
+			VEC3 newP(0) ;
+			unsigned int val = 0 ;
+			Dart vit = d ;
+			do
+			{
+				newP += position[map.phi2(vit)] ;
+				++val ;
+				vit = map.alpha1(vit) ;
+			} while(vit != d) ;
+			REAL K = sqrt3_K(val) ;
+			newP *= REAL(3) ;
+			newP -= REAL(val) * P ;
+			newP *= K / REAL(2 * val) ;
+			newP += (REAL(1) - K) * P ;
+			position[d] = newP ;
+		}
+	}
 }
 
 template <typename PFP, typename EMBV, typename EMB>

include/Algo/MovingObjects/particle_cell_2D.hpp

@@ -174,8 +174,8 @@ void ParticleCell2D<PFP>::faceState(const VEC3& current)
 	std::cout << "faceState" <<  d << std::endl;
 	#endif
 
-// 	assert(std::isfinite(m_position[0]) && std::isfinite(m_position[1]) && std::isfinite(m_position[2]));
-// 	assert(std::isfinite(current[0]) && std::isfinite(current[1]) && std::isfinite(current[2]));
+ 	assert(std::isfinite(m_position[0]) && std::isfinite(m_position[1]) && std::isfinite(m_position[2]));
+ 	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));
 
 	Dart dd=d;
@@ -236,7 +236,7 @@ void ParticleCell2D<PFP>::faceState(const VEC3& current)
 		}
 
 		//in case of numerical incoherence
-		if(m.phi_1(d)==dd) {
+		if(m.phi_1(d)==dd && wsoe==Geom::RIGHT) {
 			d = m.phi_1(d);
 			do {
 				switch (getOrientationEdge(current,d)) {

include/Algo/Parallel/parallel_foreach.h

@@ -22,7 +22,6 @@
 *                                                                              *
 *******************************************************************************/
 
-
 #include "Topology/generic/functor.h"
 
 #ifndef __PARALLEL_FOREACH__
@@ -225,16 +224,12 @@ public:
 	virtual void operator()() =0;
 };
 
+} // namespace Parallel
 
+} // namespace Algo
 
-
-}
-}	// end namespace
-}
-
+} // namespace CGoGN
 
 #include "Algo/Parallel/parallel_foreach.hpp"
 
-
 #endif
-

include/Algo/Remeshing/pliant.h

+/*******************************************************************************
+* CGoGN: Combinatorial and Geometric modeling with Generic N-dimensional Maps  *
+* version 0.1                                                                  *
+* Copyright (C) 2009, 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: https://iggservis.u-strasbg.fr/CGoGN/                              *
+* Contact information: cgogn@unistra.fr                                        *
+*                                                                              *
+*******************************************************************************/
+
+#ifndef __ALGO_GEOMETRY_REMESHING_H__
+#define __ALGO_GEOMETRY_REMESHING_H__
+
+namespace CGoGN
+{
+
+namespace Algo
+{
+
+namespace Remeshing
+{
+
+template <typename PFP>
+void pliantRemeshing(typename PFP::MAP& map, typename PFP::TVEC3& position, typename PFP::TVEC3& normal) ;
+
+} // namespace Remeshing
+