Commit 8750e5fd authored by Sylvain Thery's avatar Sylvain Thery
Browse files

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

parents 6afbce7f 4b7fbbef
......@@ -17,22 +17,21 @@ namespace Geometry
template <typename PFP>
class VoronoiDiagram {
private :
protected :
typedef typename PFP::REAL REAL;
typedef typename PFP::VEC3 VEC3;
typedef struct
{
typename std::multimap<float,Dart>::iterator it ;
bool valid ;
unsigned int region;
// unsigned int region;
Dart pathOrigin;
static std::string CGoGNnameOfType() { return "VoronoiVertexInfo" ; }
} VoronoiVertexInfo ;
typedef NoMathIOAttribute<VoronoiVertexInfo> VertexInfo ;
typename PFP::MAP& map;
const EdgeAttribute<REAL>& edgeCost; // weights on the graph edges
// const VertexAttribute<VEC3>& position; // positions
VertexAttribute<unsigned int>& regions; // region labels
std::vector<Dart> border;
std::vector<Dart> seeds;
......@@ -46,18 +45,67 @@ public :
~VoronoiDiagram ();
const std::vector<Dart>& getSeeds (){return seeds;}
void setSeeds (const std::vector<Dart>&);
void setRandomSeeds (unsigned int nbseeds);
virtual void setSeeds (const std::vector<Dart>&);
virtual void setRandomSeeds (unsigned int nbseeds);
const std::vector<Dart>& getBorder (){return border;}
void setCost (const EdgeAttribute<REAL>& c);
void computeDiagram ();
void computeDistancesWithinRegion (Dart seed);
private :
void clear ();
protected :
virtual void clear ();
void initFrontWithSeeds();
virtual void collectVertexFromFront(Dart e);
void addVertexToFront(Dart f, float d);
void updateVertexInFront(Dart f, float d);
};
template <typename PFP>
class CentroidalVoronoiDiagram : public VoronoiDiagram<PFP> {
private :
typedef typename PFP::REAL REAL;
typedef typename PFP::VEC3 VEC3;
double globalEnergy;
std::vector<VEC3> energyGrad; // gradient of the region energy at seed
VertexAttribute<REAL>& distances; // distances from the seed
VertexAttribute<Dart>& pathOrigins; // previous vertex on the shortest path from origin
VertexAttribute<REAL>& areaElts; // area element attached to each vertex
public :
CentroidalVoronoiDiagram (typename PFP::MAP& m,
const EdgeAttribute<REAL>& c,
VertexAttribute<unsigned int>& r,
VertexAttribute<REAL>& d,
VertexAttribute<Dart>& o,
VertexAttribute<REAL>& a);
~CentroidalVoronoiDiagram ();
void setSeeds (const std::vector<Dart>&);
void setRandomSeeds (unsigned int nbseeds);
void cumulateEnergy();
void cumulateEnergyAndGradients();
unsigned int moveSeedsOneEdgeNoCheck(); // returns the number of seeds that did move
// move each seed along one edge according to the energy gradient
unsigned int moveSeedsOneEdgeCheck(); // returns the number of seeds that did move
// move each seed along one edge according to the energy gradient + check that the energy decreases
unsigned int moveSeedsToMedioid(); // returns the number of seeds that did move
// move each seed to the medioid of its region
REAL getGlobalEnergy() {return globalEnergy;}
protected :
void clear();
void collectVertexFromFront(Dart e);
REAL cumulateEnergyFromRoot(Dart e);
void cumulateEnergyAndGradientFromSeed(unsigned int numSeed);
Dart selectBestNeighborFromSeed(unsigned int numSeed);
// unsigned int moveSeed(unsigned int numSeed);
};
}// end namespace Geometry
}// end namespace Algo
}// end namespace CGoGN
......
......@@ -7,6 +7,10 @@ namespace Algo
namespace Geometry
{
/***********************************************************
* class VoronoiDiagram
***********************************************************/
template <typename PFP>
VoronoiDiagram<PFP>::VoronoiDiagram (typename PFP::MAP& m, const EdgeAttribute<REAL>& p, VertexAttribute<unsigned int>& r) : map(m), edgeCost (p), regions (r), vmReached(m)
{
......@@ -22,10 +26,8 @@ VoronoiDiagram<PFP>::~VoronoiDiagram ()
template <typename PFP>
void VoronoiDiagram<PFP>::clear ()
{
border.clear();
regions.setAllValues(0);
border.clear();
seeds.clear();
front.clear();
vmReached.unmarkAll();
}
......@@ -33,14 +35,15 @@ void VoronoiDiagram<PFP>::clear ()
template <typename PFP>
void VoronoiDiagram<PFP>::setSeeds (const std::vector<Dart>& s)
{
clear();
seeds.clear();
seeds = s;
}
template <typename PFP>
void VoronoiDiagram<PFP>::setRandomSeeds (unsigned int nseeds)
{
clear();
seeds.clear();
vmReached.unmarkAll();
srand ( time(NULL) );
unsigned int n = nseeds;
......@@ -59,14 +62,16 @@ void VoronoiDiagram<PFP>::setRandomSeeds (unsigned int nseeds)
template <typename PFP>
void VoronoiDiagram<PFP>::initFrontWithSeeds ()
{
vmReached.unmarkAll();
// vmReached.unmarkAll();
clear();
for (unsigned int i = 0; i < seeds.size(); i++)
{
Dart d = seeds[i];
vmReached.mark(d);
vertexInfo[d].it = front.insert(std::pair<float,Dart>(0.0, d));
vertexInfo[d].valid = true;
vertexInfo[d].region = i;
regions[d] = i;
vertexInfo[d].pathOrigin = d;
}
}
......@@ -75,6 +80,33 @@ void VoronoiDiagram<PFP>::setCost (const EdgeAttribute<typename PFP::REAL>& c){
edgeCost = c;
}
template <typename PFP>
void VoronoiDiagram<PFP>::collectVertexFromFront(Dart e){
regions[e] = regions[vertexInfo[e].pathOrigin];
front.erase(vertexInfo[e].it);
vertexInfo[e].valid=false;
}
template <typename PFP>
void VoronoiDiagram<PFP>::addVertexToFront(Dart f, float d){
VertexInfo& vi (vertexInfo[f]);
vi.it = front.insert(std::pair<float,Dart>(d + edgeCost[f], f));
vi.valid=true;
vi.pathOrigin = map.phi2(f);
vmReached.mark(f);
}
template <typename PFP>
void VoronoiDiagram<PFP>::updateVertexInFront(Dart f, float d){
VertexInfo& vi (vertexInfo[f]);
float dist = d + edgeCost[f];
if (dist < vi.it->first)
{
front.erase(vi.it);
vi.it = front.insert(std::pair<float,Dart>(dist, f));
vi.pathOrigin = map.phi2(f);
}
}
template <typename PFP>
void VoronoiDiagram<PFP>::computeDiagram ()
......@@ -85,45 +117,392 @@ void VoronoiDiagram<PFP>::computeDiagram ()
{
Dart e = front.begin()->second;
float d = front.begin()->first;
front.erase(vertexInfo[e].it);
vertexInfo[e].valid=false;
regions[e] = vertexInfo[e].region;
collectVertexFromFront(e);
Traversor2VVaE<typename PFP::MAP> tv (map, e);
for (Dart f = tv.begin(); f != tv.end(); f=tv.next())
{
VertexInfo& vi (vertexInfo[f]);
if (vmReached.isMarked(f))
{ // f has been reached
if (vertexInfo[f].valid) // f is in the front : update
updateVertexInFront(f,d);
else // f is not in the front any more (already collected) : detect a border edge
if ( regions[f] != regions[e] ) border.push_back(f);
}
else
{ // f has not been reached : add it to the front
addVertexToFront(f,d);
}
}
}
}
template <typename PFP>
void VoronoiDiagram<PFP>::computeDistancesWithinRegion (Dart seed)
{ // TODO : this algo has not yet been tested
// init
front.clear();
vmReached.unmarkAll();
vmReached.mark(seed);
vertexInfo[seed].it = front.insert(std::pair<float,Dart>(0.0, seed));
vertexInfo[seed].valid = true;
vertexInfo[seed].pathOrigin = seed;
//compute
while ( !front.empty() )
{
if (vi.valid) // probably useless (because of distance test) but faster
Dart e = front.begin()->second;
float d = front.begin()->first;
collectVertexFromFront(e);
Traversor2VVaE<typename PFP::MAP> tv (map, e);
for (Dart f = tv.begin(); f != tv.end(); f=tv.next())
{
// float dist = d + Algo::Geometry::edgeLength<PFP>(map,f,position);
float dist = d + edgeCost[f];
if (dist < vi.it->first)
if (vmReached.isMarked(f))
{ // f has been reached
if (vertexInfo[f].valid) updateVertexInFront(f,d); // f is in the front : update
}
else
{ // f has not been reached : add it to the front
if ( regions[f] == regions[e] ) addVertexToFront(f,d);
}
}
}
}
/***********************************************************
* class CentroidalVoronoiDiagram
***********************************************************/
template <typename PFP>
CentroidalVoronoiDiagram<PFP>::CentroidalVoronoiDiagram (typename PFP::MAP& m, const EdgeAttribute<REAL>& c, VertexAttribute<unsigned int>& r, VertexAttribute<REAL>& d, VertexAttribute<Dart>& o, VertexAttribute<REAL>& a) :
VoronoiDiagram<PFP>(m,c,r), distances(d), pathOrigins(o), areaElts(a)
{
}
template <typename PFP>
CentroidalVoronoiDiagram<PFP>::~CentroidalVoronoiDiagram ()
{
}
template <typename PFP>
void CentroidalVoronoiDiagram<PFP>::clear ()
{
VoronoiDiagram<PFP>::clear();
distances.setAllValues(0.0);
}
template <typename PFP>
void CentroidalVoronoiDiagram<PFP>::collectVertexFromFront(Dart e){
distances[e] = this->vertexInfo[e].it->first;
pathOrigins[e] = this->vertexInfo[e].pathOrigin;
VoronoiDiagram<PFP>::collectVertexFromFront(e);
}
template <typename PFP>
void CentroidalVoronoiDiagram<PFP>::setSeeds (const std::vector<Dart>& s)
{
VoronoiDiagram<PFP>::setSeeds (s);
energyGrad.resize(this->seeds.size());
}
template <typename PFP>
void CentroidalVoronoiDiagram<PFP>::setRandomSeeds (unsigned int nseeds)
{
VoronoiDiagram<PFP>::setRandomSeeds (nseeds);
energyGrad.resize(this->seeds.size());
}
template <typename PFP>
void CentroidalVoronoiDiagram<PFP>::cumulateEnergy(){
globalEnergy = 0.0;
for (unsigned int i = 0; i < this->seeds.size(); i++)
{
front.erase(vi.it);
vi.it = front.insert(std::pair<float,Dart>(dist, f));
vi.region = regions[e];
cumulateEnergyFromRoot(this->seeds[i]);
globalEnergy += distances[this->seeds[i]];
}
}
template <typename PFP>
void CentroidalVoronoiDiagram<PFP>::cumulateEnergyAndGradients(){
globalEnergy = 0.0;
for (unsigned int i = 0; i < this->seeds.size(); i++)
{
cumulateEnergyAndGradientFromSeed(i);
globalEnergy += distances[this->seeds[i]];
}
}
template <typename PFP>
unsigned int CentroidalVoronoiDiagram<PFP>::moveSeedsOneEdgeNoCheck(){
unsigned int m = 0;
for (unsigned int i = 0; i < this->seeds.size(); i++)
{
Dart oldSeed = this->seeds[i];
Dart newSeed = selectBestNeighborFromSeed(i);
// move the seed
if (newSeed != oldSeed)
{
this->seeds[i] = newSeed;
m++;
}
else
}
return m;
}
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++)
{
if ( regions[f] != regions[e] )
border.push_back(f);
Dart oldSeed = this->seeds[i];
Dart newSeed = selectBestNeighborFromSeed(i);
// move the seed
if (newSeed != oldSeed)
{
REAL regionEnergy = distances[oldSeed];
this->computeDistancesWithinRegion(newSeed);
cumulateEnergyFromRoot(newSeed);
if (distances[newSeed] < regionEnergy)
{
this->seeds[i] = newSeed;
m++;
}
}
}
return m;
}
/*
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++)
{
Dart oldSeed = this->seeds[i];
Dart newSeed = selectBestNeighborFromSeed(i);
// move the seed
if (newSeed != oldSeed)
{
REAL regionEnergy = distances[oldSeed];
this->seeds[i] = newSeed;
this->computeDistancesWithinRegion(newSeed);
cumulateEnergyAndGradientFromSeed(i);
if (distances[newSeed] < regionEnergy)
m++;
else
this->seeds[i] = newSeed;
}
}
return m;
}
*/
template <typename PFP>
unsigned int CentroidalVoronoiDiagram<PFP>::moveSeedsToMedioid(){
unsigned int m = 0;
/* for (unsigned int i = 0; i < this->seeds.size(); i++)
{
// vi.it = front.insert(std::pair<float,Dart>(d + Algo::Geometry::edgeLength<PFP>(map,f,position), f));
vi.it = front.insert(std::pair<float,Dart>(d + edgeCost[f], f));
vi.valid=true;
vi.region = regions[e];
vmReached.mark(f);
Dart oldSeed = this->seeds[i];
Dart newSeed = selectBestNeighborFromSeed(i);
unsigned int seedMoved = 0;
REAL regionEnergy = distances[oldSeed];
// move the seed
while (newSeed != oldSeed)
{
this->cumulateEnergyAndGradientFromSeed(numSeed);
cumulateEnergyFromRoot(newSeed);
if (distances[newSeed] < regionEnergy)
{
this->seeds[i] = newSeed;
m++;
}
}
} */
return m;
}
template <typename PFP>
typename PFP::REAL CentroidalVoronoiDiagram<PFP>::cumulateEnergyFromRoot(Dart e){
REAL distArea = areaElts[e] * distances[e];
distances[e] = areaElts[e] * distances[e] * distances[e];
Traversor2VVaE<typename PFP::MAP> tv (this->map, e);
for (Dart f = tv.begin(); f != tv.end(); f=tv.next())
{
if ( pathOrigins[f] == this->map.phi2(f))
{
distArea += cumulateEnergyFromRoot(f);
distances[e] += distances[f];
}
}
return distArea;
}
template <typename PFP>
void CentroidalVoronoiDiagram<PFP>::cumulateEnergyAndGradientFromSeed(unsigned int numSeed){
Dart e = this->seeds[numSeed];
std::vector<Dart> v;
v.reserve(8);
std::vector<float> da;
da.reserve(8);
distances[e] = 0.0;
Traversor2VVaE<typename PFP::MAP> tv (this->map, e);
for (Dart f = tv.begin(); f != tv.end(); f=tv.next())
{
if ( pathOrigins[f] == this->map.phi2(f))
{
float distArea = cumulateEnergyFromRoot(f);
da.push_back(distArea);
distances[e] += distances[f];
v.push_back(f);
}
}
// compute the gradient
// TODO : check if the computation of grad and proj is still valid for other edgeCost than geodesic distances
VEC3 grad (0.0);
const VertexAttribute<VEC3>& pos = this->map.template getAttribute<VEC3,VERTEX>("position");
for (unsigned int j = 0; j<v.size(); ++j)
{
Dart f = v[j];
VEC3 edgeV = pos[f] - pos[this->map.phi2(f)];
edgeV.normalize();
grad += da[j] * edgeV;
}
grad /= 2.0;
energyGrad[numSeed] = grad;
}
template <typename PFP>
Dart CentroidalVoronoiDiagram<PFP>::selectBestNeighborFromSeed(unsigned int numSeed)
{
Dart e = this->seeds[numSeed];
Dart newSeed = e;
const VertexAttribute<VEC3>& pos = this->map.template getAttribute<VEC3,VERTEX>("position");
// TODO : check if the computation of grad and proj is still valid for other edgeCost than geodesic distances
float maxProj = 0.0;
Traversor2VVaE<typename PFP::MAP> tv (this->map, e);
for (Dart f = tv.begin(); f != tv.end(); f=tv.next())
{
if ( pathOrigins[f] == this->map.phi2(f))
{
VEC3 edgeV = pos[f] - pos[this->map.phi2(f)];
// edgeV.normalize();
float proj = edgeV * energyGrad[numSeed];
if (proj > maxProj)
{
maxProj = proj;
newSeed = f;
}
}
}
return newSeed;
}
/*
template <typename PFP>
unsigned int CentroidalVoronoiDiagram<PFP>::moveSeed(unsigned int numSeed){
// collect energy and compute the gradient
// cumulateEnergyAndGradientFromSeed(numSeed);
// select the new seed
Dart newSeed = selectBestNeighborFromSeed(numSeed);
// move the seed
Dart oldSeed = this->seeds[numSeed];
unsigned int seedMoved = 0;
if (newSeed != oldSeed)
{
this->seeds[numSeed] = newSeed;
seedMoved = 1;
}
return seedMoved;
}
*/
/*
template <typename PFP>
unsigned int CentroidalVoronoiDiagram<PFP>::moveSeed(unsigned int numSeed){
Dart e = this->seeds[numSeed];
unsigned int seedMoved = 0;
std::vector<Dart> v;
v.reserve(8);
std::vector<float> da;
da.reserve(8);
distances[e] = 0.0;
Traversor2VVaE<typename PFP::MAP> tv (this->map, e);
for (Dart f = tv.begin(); f != tv.end(); f=tv.next())
{
if ( pathOrigins[f] == this->map.phi2(f))
{
float distArea = cumulateEnergyFromRoot(f);
da.push_back(distArea);
distances[e] += distances[f];
v.push_back(f);
}
}
// TODO : check if the computation of grad and proj is still valid for other edgeCost than geodesic distances
VEC3 grad (0.0);
const VertexAttribute<VEC3>& pos = this->map.template getAttribute<VEC3,VERTEX>("position");
// compute the gradient
for (unsigned int j = 0; j<v.size(); ++j)
{
Dart f = v[j];
VEC3 edgeV = pos[f] - pos[this->map.phi2(f)];
edgeV.normalize();
grad += da[j] * edgeV;
}
grad /= 2.0;
float maxProj = 0.0;
// float memoForTest = 0.0;
for (unsigned int j = 0; j<v.size(); ++j)
{
Dart f = v[j];
VEC3 edgeV = pos[f] - pos[this->map.phi2(f)];
// edgeV.normalize();
float proj = edgeV * grad;
// proj -= areaElts[e] * this->edgeCost[f] * this->edgeCost[f];
if (proj > maxProj)
{
// if (numSeed==1) memoForTest = (edgeV * grad) / (areaElts[e] * this->edgeCost[f] * this->edgeCost[f]);
// CGoGNout << (edgeV * grad) / (areaElts[e] * this->edgeCost[f] * this->edgeCost[f]) * this->seeds.size() << CGoGNendl;
// CGoGNout << edgeV * grad << "\t - \t" << areaElts[e] * this->edgeCost[f] * this->edgeCost[f] << CGoGNendl;
maxProj = proj;
seedMoved = 1;
this->seeds[numSeed] = v[j];
}
}
return seedMoved;
}
*/
}// end namespace Geometry
}// end namespace Algo
}// end namespace CGoGN
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment