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

include/Algo/Geometry/voronoiDiagrams.h

@@ -3,7 +3,9 @@
 
 #include <vector>
 #include <map>
+#include <set>
 
+//#include "Topology/map/map2.h"
 #include "Topology/generic/traversor2.h"
 
 namespace CGoGN
@@ -45,12 +47,13 @@ public :
 	~VoronoiDiagram ();
 
 	const std::vector<Dart>& getSeeds (){return seeds;}
-	virtual void setSeeds (const std::vector<Dart>&);
-	virtual void setRandomSeeds (unsigned int nbseeds);
+	virtual void setSeeds_fromVector (const std::vector<Dart>&);
+	virtual void setSeeds_random (unsigned int nbseeds);
 	const std::vector<Dart>& getBorder (){return border;}
 	void setCost (const EdgeAttribute<REAL>& c);
 
-	void computeDiagram ();
+	Dart computeDiagram ();
+	virtual void computeDiagram_incremental (unsigned int nbseeds);
 	void computeDistancesWithinRegion (Dart seed);
 
 protected :
@@ -84,8 +87,9 @@ public :
 			VertexAttribute<REAL>& a);
 	~CentroidalVoronoiDiagram ();
 
-	void setSeeds (const std::vector<Dart>&);
-	void setRandomSeeds (unsigned int nbseeds);
+	void setSeeds_fromVector (const std::vector<Dart>&);
+	void setSeeds_random (unsigned int nbseeds);
+	void computeDiagram_incremental (unsigned int nbseeds);
 	void cumulateEnergy();
 	void cumulateEnergyAndGradients();
 	unsigned int moveSeedsOneEdgeNoCheck(); // returns the number of seeds that did move

include/Algo/Geometry/voronoiDiagrams.hpp

@@ -33,29 +33,48 @@ void VoronoiDiagram<PFP>::clear ()
 }
 
 template <typename PFP>
-void VoronoiDiagram<PFP>::setSeeds (const std::vector<Dart>& s)
+void VoronoiDiagram<PFP>::setSeeds_fromVector (const std::vector<Dart>& s)
 {
 	seeds.clear();
 	seeds = s;
 }
 
 template <typename PFP>
-void VoronoiDiagram<PFP>::setRandomSeeds (unsigned int nseeds)
+void VoronoiDiagram<PFP>::setSeeds_random (unsigned int nseeds)
 {
 	seeds.clear();
-	vmReached.unmarkAll();
-
 	srand ( time(NULL) );
-	unsigned int n = nseeds;
-	while (n > 0)
-	{ // TODO : correct this random init which assumes contiguous Dart table
-		Dart d = rand() % map.getNbDarts() ;
-		if (! vmReached.isMarked(d))
+	const unsigned int nbv = map.getNbCells(VERTEX);
+
+	std::set<unsigned int> myVertices ;;
+	while (myVertices.size() < nseeds)
+	{
+		myVertices.insert(rand() % nbv);
+	}
+
+	std::set<unsigned int>::iterator it = myVertices.begin();
+	unsigned int n = 0;
+	TraversorV<typename PFP::MAP> tv (map);
+	Dart dit = tv.begin();
+
+	while (it != myVertices.end())
+	{
+		while(n<*it)
 		{
-			vmReached.mark(d);
-			seeds.push_back(d);
-			n--;
+			dit = tv.next();
+			++n;
 		}
+		seeds.push_back(dit);
+		it++;
+	}
+
+	// random permutation = un-sort the seeds
+	for (unsigned int i=0; i<nseeds; i++)
+	{
+		unsigned int j = i + rand() % (nseeds - i);
+		Dart d = seeds[i];
+		seeds[i] = seeds[j];
+		seeds[j] = d;
 	}
 }
 
@@ -109,13 +128,14 @@ void VoronoiDiagram<PFP>::updateVertexInFront(Dart f, float d){
 }
 
 template <typename PFP>
-void VoronoiDiagram<PFP>::computeDiagram ()
+Dart VoronoiDiagram<PFP>::computeDiagram ()
 {
 	initFrontWithSeeds();
 
+	Dart e;
 	while ( !front.empty() )
 	{
-		Dart e = front.begin()->second;
+		e = front.begin()->second;
 		float d = front.begin()->first;
 
 		collectVertexFromFront(e);
@@ -136,6 +156,35 @@ void VoronoiDiagram<PFP>::computeDiagram ()
 			}
 		}
 	}
+	return e;
+}
+
+template <typename PFP>
+void VoronoiDiagram<PFP>::computeDiagram_incremental (unsigned int nseeds)
+{
+	seeds.clear();
+
+	// first seed
+	srand ( time(NULL) );
+	unsigned int s = rand() % map.getNbCells(VERTEX);
+	unsigned int n = 0;
+	TraversorV<typename PFP::MAP> tv (map);
+	Dart dit = tv.begin();
+	while(n<s)
+	{
+		dit = tv.next();
+		++n;
+	}
+	seeds.push_back(dit);
+
+	// add other seeds one by one
+	Dart e = computeDiagram();
+
+	for(unsigned int i = 1; i< nseeds ; i++)
+	{
+		seeds.push_back(e);
+		e = computeDiagram();
+	}
 }
 
 template <typename PFP>
@@ -205,16 +254,23 @@ void CentroidalVoronoiDiagram<PFP>::collectVertexFromFront(Dart e){
 
 
 template <typename PFP>
-void CentroidalVoronoiDiagram<PFP>::setSeeds (const std::vector<Dart>& s)
+void CentroidalVoronoiDiagram<PFP>::setSeeds_fromVector (const std::vector<Dart>& s)
+{
+	VoronoiDiagram<PFP>::setSeeds_fromVector (s);
+	energyGrad.resize(this->seeds.size());
+}
+
+template <typename PFP>
+void CentroidalVoronoiDiagram<PFP>::setSeeds_random (unsigned int nseeds)
 {
-	VoronoiDiagram<PFP>::setSeeds (s);
+	VoronoiDiagram<PFP>::setSeeds_random (nseeds);
 	energyGrad.resize(this->seeds.size());
 }
 
 template <typename PFP>
-void CentroidalVoronoiDiagram<PFP>::setRandomSeeds (unsigned int nseeds)
+void CentroidalVoronoiDiagram<PFP>::computeDiagram_incremental (unsigned int nseeds)
 {
-	VoronoiDiagram<PFP>::setRandomSeeds (nseeds);
+	VoronoiDiagram<PFP>::computeDiagram_incremental (nseeds);
 	energyGrad.resize(this->seeds.size());
 }
 
@@ -259,7 +315,6 @@ unsigned int CentroidalVoronoiDiagram<PFP>::moveSeedsOneEdgeNoCheck(){
 
 template <typename PFP>
 unsigned int CentroidalVoronoiDiagram<PFP>::moveSeedsOneEdgeCheck(){
-	// TODO : probable bug (memoire ?) car les iterations ralentissent inexplicablement
 	unsigned int m = 0;
 	for (unsigned int i = 0; i < this->seeds.size(); i++)
 	{
@@ -334,6 +389,7 @@ typename PFP::REAL CentroidalVoronoiDiagram<PFP>::cumulateEnergyFromRoot(Dart e)
 
 template <typename PFP>
 void CentroidalVoronoiDiagram<PFP>::cumulateEnergyAndGradientFromSeed(unsigned int numSeed){
+	// precondition : energyGrad.size() > numSeed
 	Dart e = this->seeds[numSeed];
 
 	std::vector<Dart> v;

include/Algo/Selection/collector.h

@@ -139,6 +139,7 @@ public:
  * collect all primitives of the connected component containing "centerDart"
  * within the sphere of radius "radius" and center "position[centerDart]"
  * (hopefully) it defines a 2-manifold (if inserting border-vertices along the border-edges)
+ * NB : is equivalent to Collector_Vertices with CollectorCriterion_VertexWithinSphere
  */
 template <typename PFP>
 class Collector_WithinSphere : public Collector<PFP>
@@ -174,6 +175,7 @@ public:
  * collect all primitives of the connected component containing "centerDart"
  * the angle between the included vertices normal vectors and the central normal vector
  * stays under a given threshold
+ * NB : is equivalent to Collector_Vertices with CollectorCriterion_VertexNormalAngle
  */
 template <typename PFP>
 class Collector_NormalAngle : public Collector<PFP>
@@ -198,6 +200,105 @@ public:
 	void collectBorder(Dart d) ;
 };
 
+/*********************************************************
+ * Collector Criterions
+ *********************************************************/
+class CollectorCriterion
+{
+public :
+	CollectorCriterion() {};
+	virtual void init(Dart center) = 0;
+	virtual bool isInside(Dart d) = 0;
+
+};
+
+template <typename PFP>
+class CollectorCriterion_VertexNormalAngle : public CollectorCriterion
+{ // tests if the angle between vertex normals is below some threshold
+private :
+	typedef typename PFP::VEC3 VEC3;
+	typedef typename PFP::REAL REAL;
+
+	const VertexAttribute<VEC3> & vertexNormals;
+	REAL threshold;
+	VEC3 centerNormal;
+public :
+	CollectorCriterion_VertexNormalAngle(const VertexAttribute<VEC3> & n, REAL th) :
+		vertexNormals(n), threshold(th), centerNormal(0) {}
+
+	void init (Dart center) {centerNormal = vertexNormals[center];}
+	bool isInside (Dart d) {
+		return ( Geom::angle(centerNormal, vertexNormals[d]) < threshold);
+	}
+};
+
+template <typename PFP>
+class CollectorCriterion_TriangleNormalAngle : public CollectorCriterion
+{ // tests if the angle between vertex normals is below some threshold
+private :
+	typedef typename PFP::VEC3 VEC3;
+	typedef typename PFP::REAL REAL;
+
+	const FaceAttribute<VEC3> & faceNormals;
+	REAL threshold;
+	VEC3 centerNormal;
+public :
+	CollectorCriterion_TriangleNormalAngle(const FaceAttribute<VEC3> & n, REAL th) :
+		faceNormals(n), threshold(th), centerNormal(0) {}
+
+	void init (Dart center) {centerNormal = faceNormals[center];}
+	bool isInside (Dart d) {
+		return ( Geom::angle(centerNormal, faceNormals[d]) < threshold);
+	}
+};
+
+template <typename PFP>
+class CollectorCriterion_VertexWithinSphere : public CollectorCriterion
+{ // tests if the distance between vertices is below some threshold
+private :
+	typedef typename PFP::VEC3 VEC3;
+	typedef typename PFP::REAL REAL;
+
+	const VertexAttribute<VEC3> & vertexPositions;
+	REAL threshold;
+	VEC3 centerPosition;
+public :
+	CollectorCriterion_VertexWithinSphere(const VertexAttribute<VEC3> & p, REAL th) :
+		vertexPositions(p), threshold(th), centerPosition(0) {}
+
+	void init (Dart center) {centerPosition = vertexPositions[center];}
+	bool isInside (Dart d) {
+		return (vertexPositions[d] - centerPosition).norm() < threshold ;
+	}
+};
+
+
+/*********************************************************
+ * Collector Vertices
+ *********************************************************/
+
+/*
+ * collect all vertices of the connected component containing "centerDart"
+ * within a distance to centerDart defined by the CollectorCriterion
+ * (hopefully) it defines a 2-manifold (if inserting border-vertices along the border-edges)
+ */
+template <typename PFP>
+class Collector_Vertices : public Collector<PFP>
+{
+protected:
+	CollectorCriterion & crit;
+
+public:
+	Collector_Vertices(typename PFP::MAP& m, CollectorCriterion& c, unsigned int thread=0) :
+		Collector<PFP>(m, thread),
+		crit(c)
+	{}
+
+	void collectAll(Dart d);
+	void collectBorder(Dart d);
+};
+
+
 /*********************************************************
  * Collector Normal Angle (Triangles)
  *********************************************************/
@@ -206,6 +307,7 @@ public:
  * collect all primitives of the connected component containing "centerDart"
  * the angle between the included triangles normal vectors and the central normal vector
  * stays under a given threshold
+ * NB : is equivalent to Collector_Triangles with CollectorCriterion_TriangleNormalAngle
  */
 template <typename PFP>
 class Collector_NormalAngle_Triangles : public Collector<PFP>
@@ -230,6 +332,78 @@ public:
 	void collectBorder(Dart d) ;
 };
 
+/*********************************************************
+ * Collector Triangles
+ *********************************************************/
+
+/*
+ * collect all triangles of the connected component containing "centerDart"
+ * within a distance to centerDart defined by the CollectorCriterion
+ */
+template <typename PFP>
+class Collector_Triangles : public Collector<PFP>
+{
+protected:
+	CollectorCriterion & crit;
+
+public:
+	Collector_Triangles(typename PFP::MAP& m, CollectorCriterion& c, unsigned int thread=0) :
+		Collector<PFP>(m,thread), crit(c)
+	{}
+
+	void collectAll(Dart d) ;
+	void collectBorder(Dart d) ;
+};
+
+/*********************************************************
+ * Collector Dijkstra_Vertices
+ *********************************************************/
+
+/*
+ * collect all primitives of the connected component containing "centerDart"
+ * within a distance < maxDist (the shortest path follows edges)
+ * the edge length is specified in edge_cost attribute
+ */
+
+template <typename PFP>
+class Collector_Dijkstra_Vertices : public Collector<PFP>
+{
+protected:
+	const EdgeAttribute<typename PFP::REAL>& edge_cost;
+	typename PFP::REAL maxDist;
+
+	typedef struct
+	{
+		typename std::multimap<float,Dart>::iterator it ;
+		bool valid ;
+		static std::string CGoGNnameOfType() { return "DijkstraVertexInfo" ; }
+	} DijkstraVertexInfo ;
+	typedef NoMathIOAttribute<DijkstraVertexInfo> VertexInfo ;
+
+	VertexAttribute<VertexInfo> vertexInfo ;
+
+	std::multimap<float,Dart> front ;
+
+public:
+	Collector_Dijkstra_Vertices(typename PFP::MAP& m, const EdgeAttribute<typename PFP::REAL>& c, typename PFP::REAL d = 0, unsigned int thread=0) :
+		Collector<PFP>(m,thread),
+		edge_cost(c),
+		maxDist(d)
+	{
+		vertexInfo = m.template addAttribute<VertexInfo, VERTEX>("vertexInfo");
+	}
+	~Collector_Dijkstra_Vertices(){
+		this->map.removeAttribute(vertexInfo);
+	}
+	inline void init (Dart d) {Collector<PFP>::init(d); front.clear();}
+	inline void setMaxDistance(typename PFP::REAL d) { maxDist = d; }
+	inline typename PFP::REAL getMaxDist() const { return maxDist; }
+
+	void collectAll(Dart d);
+	void collectBorder(Dart d);
+};
+
+
 /*********************************************************
  * Collector Dijkstra
  *********************************************************/
@@ -281,6 +455,7 @@ private :
 };
 
 
+
 } // namespace Selection
 
 } // namespace Algo

include/Algo/Selection/collector.hpp

@@ -402,6 +402,121 @@ void Collector_NormalAngle<PFP>::collectBorder(Dart d)
 	this->insideVertices.clear();
 }
 
+/*********************************************************
+ * Collector Vertices
+ *********************************************************/
+
+template <typename PFP>
+void Collector_Vertices<PFP>::collectAll(Dart d)
+{
+	typedef typename PFP::VEC3 VEC3;
+	typedef typename PFP::REAL REAL;
+
+	crit.init(d);
+	this->init(d);
+	this->insideEdges.reserve(32);
+	this->insideFaces.reserve(32);
+	this->border.reserve(32);
+
+	CellMarkerStore<VERTEX> vm(this->map, this->m_thread);	// mark the collected inside-vertices
+	CellMarkerStore<EDGE> em(this->map, this->m_thread);	// mark the collected inside-edges + border-edges
+	CellMarkerStore<FACE> fm(this->map, this->m_thread);	// mark the collected inside-faces + border-faces
+
+	this->insideVertices.push_back(this->centerDart);
+	vm.mark(this->centerDart);
+
+	unsigned int i = 0;
+	while (i < this->insideVertices.size())
+	{
+		Dart end = this->insideVertices[i];
+		Dart e = end;
+		do
+		{
+			if (! em.isMarked(e) || ! fm.isMarked(e)) // are both tests useful ?
+			{
+				const Dart f = this->map.phi1(e);
+				const Dart g = this->map.phi1(f);
+
+				if (! crit.isInside(f))
+				{
+					this->border.push_back(e); // add to border
+					em.mark(e);
+					fm.mark(e); // is it useful ?
+				}
+				else
+				{
+					if (! vm.isMarked(f))
+					{
+						this->insideVertices.push_back(f);
+						vm.mark(f);
+					}
+					if (! em.isMarked(e))
+					{
+						this->insideEdges.push_back(e);
+						em.mark(e);
+					}
+					if (! fm.isMarked(e) && crit.isInside(g))
+					{
+						this->insideFaces.push_back(e);
+						fm.mark(e);
+					}
+				}
+			}
+			e = this->map.phi2_1(e);
+		} while (e != end);
+		++i;
+	}
+}
+
+template <typename PFP>
+void Collector_Vertices<PFP>::collectBorder(Dart d)
+{
+	typedef typename PFP::VEC3 VEC3;
+	typedef typename PFP::REAL REAL;
+
+	crit.init(d);
+	this->init(d);
+	this->border.reserve(128);
+	this->insideVertices.reserve(128);
+
+	CellMarkerStore<VERTEX> vm(this->map, this->m_thread);	// mark the collected inside-vertices
+	CellMarkerStore<EDGE> em(this->map, this->m_thread);	// mark the collected inside-edges + border-edges
+
+	this->insideVertices.push_back(this->centerDart);
+	vm.mark(this->centerDart);
+
+	unsigned int i = 0;
+	while (i < this->insideVertices.size())
+	{
+		Dart end = this->insideVertices[i];
+		Dart e = end;
+		do
+		{
+			if ( ! em.isMarked(e) )
+			{
+				const Dart f = this->map.phi1(e);
+
+				if (! crit.isInside(f))
+				{
+					this->border.push_back(e); // add to border
+				}
+				else
+				{
+					if (! vm.isMarked(f))
+					{
+						this->insideVertices.push_back(f);
+						vm.mark(f);
+					}
+				}
+				em.mark(e);
+			}
+			e = this->map.phi2_1(e);
+		} while (e != end);
+		++i;
+	}
+	this->insideVertices.clear();
+}
+
 
 /*********************************************************
  * Collector Normal Angle (Triangles)
@@ -553,6 +668,307 @@ void Collector_NormalAngle_Triangles<PFP>::collectBorder(Dart d)
 	this->insideFaces.clear();
 }
 
+/*********************************************************
+ * Collector Triangles
+ *********************************************************/
+
+template <typename PFP>
+void Collector_Triangles<PFP>::collectAll(Dart d)
+{
+	typedef typename PFP::VEC3 VEC3;
+	typedef typename PFP::REAL REAL;
+
+	crit.init(d);
+	this->init(d);
+	this->insideVertices.reserve(32);
+	this->insideEdges.reserve(32);
+	this->insideFaces.reserve(32);
+	this->border.reserve(32);
+
+	CellMarkerStore<FACE> fm(this->map, this->m_thread);	// mark the collected inside-faces + front-faces
+	CellMarkerStore<FACE> fminside(this->map, this->m_thread);	// mark the collected inside-faces
+
+	std::queue<Dart> front;
+	front.push(this->centerDart);
+	fm.mark(this->centerDart);
+
+	while ( !front.empty() ) // collect inside faces
+	{
+		Dart f = front.front();
+		front.pop();
+
+		if (crit.isInside(f))
+		{ // collect this face and add adjacent faces to the front
+			this->insideFaces.push_back(f);
+			fminside.mark(f);
+			Traversor2FFaE<typename PFP::MAP> t (this->map, f) ;
+			for (Dart it = t.begin(); it != t.end(); it=t.next())
+			{
+				if (!fm.isMarked(it))
+				{
+					front.push(it);
+					fm.mark(it);
+				}
+			}
+		}
+	}
+
+	CellMarkerStore<VERTEX> vm(this->map, this->m_thread);	// mark inside-vertices and border-vertices
+	CellMarkerStore<EDGE> em(this->map, this->m_thread);	// mark inside-edges and border-edges
+	std::vector<Dart>::iterator f_it;
+	for (f_it = this->insideFaces.begin(); f_it != this->insideFaces.end(); f_it++)
+	{ // collect insideVertices, insideEdges, and border
+		Traversor2FE<typename PFP::MAP> te (this->map, *f_it) ;
+		for (Dart it = te.begin(); it != te.end(); it=te.next())
+		{ // collect insideEdges and border
+			if (!em.isMarked(it))
+			{
+				em.mark(it);
+				if (this->map.isBoundaryEdge(it))
+					this->border.push_back(it);
+				else if ( fminside.isMarked(it) && fminside.isMarked(this->map.phi2(it)) )
+					this->insideEdges.push_back(it);
+				else
+					this->border.push_back(it);
+			}
+		}
+
+		Traversor2FV<typename PFP::MAP> tv (this->map, *f_it) ;
+		for (Dart it = tv.begin(); it != tv.end(); it=tv.next())
+		{ // collect insideVertices
+			if (!vm.isMarked(it))
+			{
+				vm.mark(it);
+				this->insideVertices.push_back(it);
+			}