Merge cgogn:~untereiner/CGoGN

Conflicts:
	include/Algo/Render/GL2/explodeVolumeRender.hpp

README.TXT

 Dépendences Linux:
 installer les paquets suivants:
-cmake libXi-dev libXmu-dev libglew-dev libxml2-dev libboost-all-dev zlib1g-dev qt4-designer qt4-dev-tools uuid-dev libgsl0-dev
+cmake libXi-dev libXmu-dev libglew-dev libxml2-dev libboost-all-dev zlib1g-dev qt4-designer qt4-dev-tools uuid-dev libgsl0-dev libsuitesparse-dev
 
 Pour compiler CGoGN:
 - aller dans ThirdParty, cd build, taper "cmake .", puis make ( avec -j x si vous avez x core sur votre machine)
@@ -41,7 +41,7 @@ https://iggservis.u-strasbg.fr/Data/data.zip
 
 Linux dependencies:
 install the following packages:
-ccmake libXi-dev libXmu-dev libglew-dev libxml2-dev libboost-all-dev zlib1g-dev qt4-designer qt4-dev-tools uuid-dev
+ccmake libXi-dev libXmu-dev libglew-dev libxml2-dev libboost-all-dev zlib1g-dev qt4-designer qt4-dev-tools uuid-dev libgsl0-dev libsuitesparse-dev
 
 To compile CGoGN:
 - Go ThirdParty, cd build, type "cmake .." and then make (with -j x if you have x core on your machine)

include/Algo/Decimation/halfEdgeSelector.h

@@ -517,8 +517,8 @@ private:
 
 
 public:
-	HalfEdgeSelector_ColorExperimental(MAP& m, VertexAttribute<typename PFP::VEC3>& pos, std::vector<ApproximatorGen<PFP>*>& approx, const FunctorSelect& select = allDarts) :
-		EdgeSelector<PFP>(m, pos, approx, select),
+	HalfEdgeSelector_ColorExperimental(MAP& m, VertexAttribute<typename PFP::VEC3>& pos, std::vector<ApproximatorGen<PFP>*>& approx) :
+		EdgeSelector<PFP>(m, pos, approx),
 		m_approxindex_pos(-1),
 		m_attrindex_pos(-1),
 		m_approxindex_color(-1),
@@ -612,8 +612,8 @@ private:
 	REAL computeSquaredLightfieldDifference(const Dart& d1, const Dart& d2) ;
 
 public:
-	HalfEdgeSelector_LFexperimental(MAP& m, VertexAttribute<typename PFP::VEC3>& pos, std::vector<ApproximatorGen<PFP>*>& approx, const FunctorSelect& select = allDarts) :
-		EdgeSelector<PFP>(m, pos, approx, select),
+	HalfEdgeSelector_LFexperimental(MAP& m, VertexAttribute<typename PFP::VEC3>& pos, std::vector<ApproximatorGen<PFP>*>& approx) :
+		EdgeSelector<PFP>(m, pos, approx),
 		m_approxindex_pos(-1),
 		m_attrindex_pos(-1),
 		m_approxindex_FT(-1),

include/Algo/Filtering/tools.h

@@ -134,7 +134,7 @@ void computeNoise(typename PFP::MAP& map, long amount, const VertexAttribute<typ
 
 } //namespace Filtering
 
-}
+} //namespace Surface
 
 } //namespace Algo
 

include/Algo/Geometry/centroid.h

@@ -303,7 +303,6 @@ void computeCentroidELWVolumes(typename PFP::MAP& map,
  * @param map the map
  * @param position vertex attribute of position
  * @param vertex_centroid vertex attribute to store the centroids
- * @param select the selector
  */
 template <typename PFP>
 void computeNeighborhoodCentroidVertices(typename PFP::MAP& map,

include/Algo/Geometry/laplacian.h

@@ -83,6 +83,26 @@ void computeCotanWeightEdges(
 
 } // namespace Surface
 
+namespace Volume
+{
+
+namespace Geometry
+{
+
+template <typename PFP, typename ATTR_TYPE>
+ATTR_TYPE computeLaplacianTopoVertex( typename PFP::MAP& map, Dart d, const VertexAttribute<ATTR_TYPE>& attr) ;
+
+template <typename PFP, typename ATTR_TYPE>
+void computeLaplacianTopoVertices(typename PFP::MAP& map, const VertexAttribute<ATTR_TYPE>& attr,
+									VertexAttribute<ATTR_TYPE>& laplacian) ;
+
+
+
+} // namespace Geometry
+
+} // namespace Volume
+
+
 } // namespace Algo
 
 } // namespace CGoGN

include/Algo/Geometry/laplacian.hpp

@@ -151,6 +151,46 @@ void computeCotanWeightEdges(
 
 } // namespace Surface
 
+
+namespace Volume
+{
+
+namespace Geometry
+{
+
+template <typename PFP, typename ATTR_TYPE>
+ATTR_TYPE computeLaplacianTopoVertex(typename PFP::MAP& map, Dart d, const VertexAttribute<ATTR_TYPE>& attr)
+{
+	ATTR_TYPE l(0) ;
+	ATTR_TYPE value = attr[d] ;
+	unsigned int wSum = 0 ;
+
+	Traversor3VE<typename PFP::MAP> t(map, d) ;
+	for(Dart it = t.begin(); it != t.end(); it = t.next())
+	{
+		l += attr[map.phi1(it)] - value ;
+		++wSum ;
+	}
+
+	l /= wSum ;
+	return l ;
+}
+
+template <typename PFP, typename ATTR_TYPE>
+void computeLaplacianTopoVertices(typename PFP::MAP& map, const VertexAttribute<ATTR_TYPE>& attr,
+									VertexAttribute<ATTR_TYPE>& laplacian)
+{
+	TraversorV<typename PFP::MAP> t(map) ;
+	for(Dart d = t.begin(); d != t.end(); d = t.next())
+		laplacian[d] = computeLaplacianTopoVertex<PFP, ATTR_TYPE>(map, d, attr) ;
+}
+
+
+} // namespace Geometry
+
+} // namespace Volume
+
+
 } // namespace Algo
 
 } // namespace CGoGN

include/Algo/Modelisation/subdivision3.hpp

@@ -26,6 +26,9 @@
 #include "Algo/Modelisation/subdivision.h"
 #include "Algo/Modelisation/extrusion.h"
 #include "Geometry/intersection.h"
+#include "NL/nl.h"
+#include "Algo/LinearSolving/basic.h"
+#include "Algo/Geometry/laplacian.h"
 
 namespace CGoGN
 {
@@ -631,12 +634,35 @@ void catmullClarkVol(typename PFP::MAP& map, EMBV& attributs)
 //	}
 }
 
+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 sqrt3Vol(typename PFP::MAP& map, VertexAttribute<typename PFP::VEC3>& position)
 {
 	DartMarkerStore m(map);
 
+	DartMarkerStore newBoundaryV(map);
+
 	//
 	// 1-4 flip of all tetrahedra
 	//
@@ -697,6 +723,8 @@ void sqrt3Vol(typename PFP::MAP& map, VertexAttribute<typename PFP::VEC3>& posit
 
 			Dart dres = Volume::Modelisation::Tetrahedralization::flip1To3<PFP>(map, dit);
 			position[dres] = faceCenter;
+
+			newBoundaryV.markOrbit<VERTEX>(dres);
 		}
 	}
 
@@ -714,6 +742,147 @@ void sqrt3Vol(typename PFP::MAP& map, VertexAttribute<typename PFP::VEC3>& posit
 
 		}
 	}
+
+	TraversorV<typename PFP::MAP> tVg(map,selected);
+	for(Dart dit = tVg.begin() ; dit != tVg.end() ; dit = tVg.next())
+	{
+		if(map.isBoundaryVertex(dit) && !newBoundaryV.isMarked(dit))
+		{
+			typename PFP::VEC3 P = position[dit] ;
+			typename PFP::VEC3 newP(0) ;
+			unsigned int val = 0 ;
+			Dart vit = dit ;
+			do
+			{
+				newP += position[map.phi_1(map.phi2(map.phi1(vit)))] ;
+				++val ;
+				vit = map.phi2(map.phi_1(vit)) ;
+			} while(vit != dit) ;
+			typename PFP::REAL K = sqrt3_K(val) ;
+			newP *= typename PFP::REAL(3) ;
+			newP -= typename PFP::REAL(val) * P ;
+			newP *= K / typename PFP::REAL(2 * val) ;
+			newP += (typename PFP::REAL(1) - K) * P ;
+			position[dit] = newP ;
+		}
+	}
+
+	//AutoVertexAttribute laplacian qui est une copie de position
+
+//	TraversorV<typename PFP::MAP> tVg2(map,selected);
+//	for(Dart dit = tVg2.begin() ; dit != tVg2.end() ; dit = tVg2.next())
+//	{
+//		if(!map.isBoundaryVertex(dit))
+//		{
+//			//modifie laplacian
+//		}
+//	}
+
+	//echange lapaclian et position
+
+
+	VertexAutoAttribute<typename PFP::VEC3> diffCoord(map, "diffCoord");
+	Algo::Volume::Geometry::computeLaplacianTopoVertices<PFP>(map, position, diffCoord) ;
+
+	VertexAutoAttribute<unsigned int> vIndex(map, "vIndex");
+
+	unsigned int nb_vertices = map.template computeIndexCells<VERTEX>(vIndex);
+
+
+	CellMarker<VERTEX> lockingMarker(map);
+
+	TraversorV<typename PFP::MAP> tv(map);
+	for(Dart dit = tv.begin() ; dit != tv.end() ; dit = tv.next())
+	{
+		if(!lockingMarker.isMarked(dit) && map.isBoundaryVertex(dit))
+			lockingMarker.mark(dit);
+	}
+
+
+	NLContext nlContext = nlNewContext();
+	nlSolverParameteri(NL_NB_VARIABLES, nb_vertices);
+	nlSolverParameteri(NL_LEAST_SQUARES, NL_TRUE);
+	nlSolverParameteri(NL_SOLVER, NL_CHOLMOD_EXT);
+
+
+	nlMakeCurrent(nlContext);
+	if(nlGetCurrentState() == NL_STATE_INITIAL)
+		nlBegin(NL_SYSTEM) ;
+
+	for(int coord = 0; coord < 3; ++coord)
+	{
+		LinearSolving::setupVariables<PFP>(map, vIndex, lockingMarker, position, coord) ;
+		nlBegin(NL_MATRIX) ;
+		LinearSolving::addRowsRHS_Laplacian_Topo<PFP>(map, vIndex, diffCoord, coord) ;
+//		LinearSolving::addRowsRHS_Laplacian_Cotan<PFP>(*map, perMap->vIndex, perMap->edgeWeight, perMap->vertexArea, perMap->diffCoord, coord) ;
+		nlEnd(NL_MATRIX) ;
+		nlEnd(NL_SYSTEM) ;
+		nlSolve() ;
+		LinearSolving::getResult<PFP>(map, vIndex, position, coord) ;
+		nlReset(NL_TRUE) ;
+	}
+
+	nlDeleteContext(nlContext);
+
+
+//
+//	float weight = 1.0;
+
+//	LinearSolving::LinearSolver<PFP::REAL> ls(nbV);
+//	ls.set_least_squares(true);
+
+//	for(unsigned int coord = 0 ; coord < 3 ; ++coord)
+//	{
+//		std::cout << "coord " << coord << std::flush;
+
+//		TraversorV<PFP::MAP> tv(map);
+//		for(Dart dit = tv.begin() ; dit != tv.end() ; dit = tv.next())
+//		{
+//			ls.variable(indexV[dit]).set_value(position[dit][coord]);
+
+//			if(map.isBoundaryVertex(dit))
+//				ls.variable(indexV[dit]).lock();
+//		}
+//		std::cout << "... variables set... " << std::flush;
+
+//		ls.begin_system();
+
+//		TraversorV<PFP::MAP> tv2(map);
+//		for(Dart dit = tv2.begin() ; dit != tv2.end() ; dit = tv2.next())
+//		{
+//			if(!map.isBoundaryVertex(dit))
+//			{
+//				ls.begin_row();
+
+//				float sum = 0;
+//				Traversor3VVaE<PFP::MAP> tvvae(map, dit);
+//				for(Dart ditvvae = tvvae.begin() ; ditvvae != tvvae.end() ; ditvvae = tvvae.next())
+//				{
+//					ls.add_coefficient(indexV[ditvvae],weight);
+//					sum += weight;
+//				}
+
+//				ls.add_coefficient(indexV[dit],-sum);
+//				ls.normalize_row();
+
+//				ls.end_row();
+//			}
+//		}
+
+//		ls.end_system();
+//		std::cout << "... system built... " << std::flush;
+//		ls.solve();
+//		std::cout << "... system solved... " << std::flush;
+
+//		TraversorV<PFP::MAP> tv3(map);
+//		for(Dart dit = tv3.begin() ; dit != tv3.end() ; dit = tv3.next())
+//		{
+//			position[dit][coord] = ls.variable(indexV[dit]).value();
+//		}
+
+//		ls.reset(false);
+//		std::cout << "... done" << std::endl;
+//	}
 }
 
 template <typename PFP>
@@ -733,11 +902,11 @@ void computeDual(typename PFP::MAP& map, VertexAttribute<typename PFP::VEC3>& po
 }
 
 
-} //namespace Modelisation
+} // namespace Modelisation
 
-} // volume
+} // namespace volume
 
-} //namespace Algo
+} // namespace Algo
 
-} //namespace CGoGN
+} // namespace CGoGN
 

include/Algo/Multiresolution/Map2MR/Filters/bertram.h

@@ -48,55 +48,30 @@ namespace Primal
 namespace Filters
 {
 
-//
-// Synthesis
-//
+/*********************************************************************************
+ *                           ANALYSIS FILTERS
+ *********************************************************************************/
 
 //w-lift(a)
 template <typename PFP>
-class Ber02OddSynthesisFilter : public Algo::MR::Filter
+class Ber02OddAnalysisFilter : public Algo::MR::Filter
 {
 protected:
 	typename PFP::MAP& m_map ;
 	VertexAttribute<typename PFP::VEC3>& m_position ;
 	typename PFP::VEC3::DATA_TYPE m_a;
 
-
 public:
-	Ber02OddSynthesisFilter(typename PFP::MAP& m, VertexAttribute<typename PFP::VEC3>& p, typename PFP::VEC3::DATA_TYPE a) : m_map(m), m_position(p), m_a(a)
+	Ber02OddAnalysisFilter(typename PFP::MAP& m, VertexAttribute<typename PFP::VEC3>& p, typename PFP::VEC3::DATA_TYPE a) : m_map(m), m_position(p), m_a(a)
 	{}
 
 	void operator() ()
 	{
-		TraversorF<typename PFP::MAP> travF(m_map) ;
-		for (Dart d = travF.begin(); d != travF.end(); d = travF.next())
-		{
-			typename PFP::VEC3 vf(0.0);
-			typename PFP::VEC3 ef(0.0);
-
-			unsigned int count = 0;
-			Traversor2FE<typename PFP::MAP> travFE(m_map, d);
-			for (Dart dit = travFE.begin(); dit != travFE.end(); dit = travFE.next())
-			{
-				vf += m_position[dit];
-				m_map.incCurrentLevel();
-				ef += m_position[m_map.phi1(dit)];
-				m_map.decCurrentLevel();
-				++count;
-			}
-			ef /= count;
-			ef *= 4.0 * m_a;
 
-			vf /= count;
-			vf *= 4.0 * m_a * m_a;
-
-			m_map.incCurrentLevel() ;
-			Dart midF = m_map.phi1(m_map.phi1(d));
-			m_position[midF] += vf + ef ;
-			m_map.decCurrentLevel() ;
-
-		}
+	}
 
+	void operator() (bool filtering)
+	{
 		TraversorE<typename PFP::MAP> travE(m_map) ;
 		for (Dart d = travE.begin(); d != travE.end(); d = travE.next())
 		{
@@ -105,13 +80,10 @@ public:
 
 			m_map.incCurrentLevel() ;
 			Dart midV = m_map.phi1(d) ;
-			m_position[midV] += ve;
+			m_position[midV] -= ve;
 			m_map.decCurrentLevel() ;
 		}
-	}
 
-	void operator() (bool filtering)
-	{
 		TraversorF<typename PFP::MAP> travF(m_map) ;
 		for (Dart d = travF.begin(); d != travF.end(); d = travF.next())
 		{
@@ -136,46 +108,15 @@ public:
 
 			m_map.incCurrentLevel() ;
 			Dart midF = m_map.phi1(m_map.phi1(d));
-			if(filtering)
-			{
-				std::cout << "sans +=" << std::endl;
-				m_position[midF] = vf + ef ;
-			}
-			else
-			{
-				std::cout << "avec +=" << std::endl;
-				m_position[midF] += vf + ef ;
-			}
-			m_map.decCurrentLevel() ;
-
-		}
-
-		TraversorE<typename PFP::MAP> travE(m_map) ;
-		for (Dart d = travE.begin(); d != travE.end(); d = travE.next())
-		{
-			typename PFP::VEC3 ve = (m_position[d] + m_position[m_map.phi1(d)]) * typename PFP::REAL(0.5);
-			ve *= 2.0 * m_a;
-
-			m_map.incCurrentLevel() ;
-			Dart midV = m_map.phi1(d) ;
-			if(filtering)
-			{
-				std::cout << "sans +=" << std::endl;
-				m_position[midV] = ve;
-			}
-			else
-			{
-				std::cout << "avec +=" << std::endl;
-				m_position[midV] += ve;
-			}
+			m_position[midF] -= vf + ef ;
 			m_map.decCurrentLevel() ;
 		}
 	}
-} ;
+};
 
 // s-lift(a)
 template <typename PFP>
-class Ber02EvenSynthesisFilter : public Algo::MR::Filter
+class Ber02EvenAnalysisFilter : public Algo::MR::Filter
 {
 protected:
 	typename PFP::MAP& m_map ;
@@ -183,50 +124,16 @@ protected:
 	typename PFP::VEC3::DATA_TYPE m_a;
 
 public:
-	Ber02EvenSynthesisFilter(typename PFP::MAP& m, VertexAttribute<typename PFP::VEC3>& p, typename PFP::VEC3::DATA_TYPE a) : m_map(m), m_position(p), m_a(a)
+	Ber02EvenAnalysisFilter(typename PFP::MAP& m, VertexAttribute<typename PFP::VEC3>& p, typename PFP::VEC3::DATA_TYPE a) : m_map(m), m_position(p), m_a(a)
 	{}
 
 	void operator() ()
 	{
-		TraversorV<typename PFP::MAP> travV(m_map);
-		for(Dart d = travV.begin() ; d != travV.end() ; d = travV.next())
-		{
-			typename PFP::VEC3 ev(0.0);
-			typename PFP::VEC3 fv(0.0);
-			if(m_map.isBoundaryVertex(d))
-			{
-				Dart db = m_map.findBoundaryEdgeOfVertex(d);
-				m_map.incCurrentLevel() ;
-				ev += (m_position[m_map.phi1(db)] + m_position[m_map.phi_1(db)]) * typename PFP::REAL(0.5);
-				m_map.decCurrentLevel() ;
-				ev *= 2 * m_a;
-
-				m_position[db] += ev;
-			}
-			else
-			{
-				unsigned int count = 0;
-				Traversor2VF<typename PFP::MAP> travVF(m_map,d);
-				for(Dart dit = travVF.begin(); dit != travVF.end() ; dit = travVF.next())
-				{
-					m_map.incCurrentLevel() ;
-
-					Dart midEdgeV = m_map.phi1(dit);
-					ev += m_position[midEdgeV];
-					fv += m_position[m_map.phi1(midEdgeV)];
-
-					m_map.decCurrentLevel() ;
-					++count;
-				}
-				fv /= count;
-				fv *= 4 * m_a * m_a;
 
-				ev /= count;
-				ev *= 4 * m_a;
-				m_position[d] += fv + ev;
-			}
-		}
+	}
 
+	void operator() (bool filtering)
+	{
 		TraversorE<typename PFP::MAP> travE(m_map);
 		for(Dart d = travE.begin() ; d != travE.end() ; d = travE.next())
 		{
@@ -234,7 +141,7 @@ public:
 			{
 				unsigned int count = 0;
 
-				typename PFP::VEC3 fe(0.0);
+				typename PFP::VEC3 fe(0);
 				Traversor2EF<typename PFP::MAP> travEF(m_map, d);
 				for(Dart dit = travEF.begin() ; dit != travEF.end() ; dit = travEF.next())
 				{
@@ -250,14 +157,11 @@ public:
 
 				m_map.incCurrentLevel() ;
 				Dart midF = m_map.phi1(d);
-				m_position[midF] += fe;
+				m_position[midF] -= fe;
 				m_map.decCurrentLevel() ;
 			}
 		}
-	}
 
-	void operator() (bool filtering)
-	{
 		TraversorV<typename PFP::MAP> travV(m_map);
 		for(Dart d = travV.begin() ; d != travV.end() ; d = travV.next())
 		{
@@ -267,24 +171,16 @@ public:
 			{
 				Dart db = m_map.findBoundaryEdgeOfVertex(d);
 				m_map.incCurrentLevel() ;
-				ev += (m_position[m_map.phi1(db)] + m_position[m_map.phi_1(db)]) * typename PFP::REAL(0.5);
+				ev = (m_position[m_map.phi1(db)] + m_position[m_map.phi_1(db)]);
 				m_map.decCurrentLevel() ;
-				ev *= 2 * m_a;
+				ev *= m_a;
 
-				if(filtering)
-				{
-					std::cout << "sans +=" << std::endl;
-					m_position[db] = ev;
-				}
-				else
-				{
-					std::cout << "avec +=" << std::endl;
-					m_position[db] += ev;
-				}
+				m_position[d] -= ev;
 			}
 			else
 			{
 				unsigned int count = 0;
+
 				Traversor2VF<typename PFP::MAP> travVF(m_map,d);
 				for(Dart dit = travVF.begin(); dit != travVF.end() ; dit = travVF.next())
 				{
@@ -302,61 +198,16 @@ public:
 
 				ev /= count;
 				ev *= 4 * m_a;
-				if(filtering)
-				{
-					std::cout << "sans +=" << std::endl;
-					m_position[d] = fv + ev;
-				}
-				else
-				{
-					std::cout << "avec +=" << std::endl;
-					m_position[d] += fv + ev;