diff --git a/include/Algo/Geometry/voronoiDiagrams.h b/include/Algo/Geometry/voronoiDiagrams.h index e1f22066ea30d9118e30c4f9783975db2b357a59..f359ff8f9b9104f53a8a49b4a8aa9207abad2b14 100644 --- a/include/Algo/Geometry/voronoiDiagrams.h +++ b/include/Algo/Geometry/voronoiDiagrams.h @@ -67,6 +67,8 @@ private : typedef typename PFP::REAL REAL; typedef typename PFP::VEC3 VEC3; + double globalEnergy; + VertexAttribute& distances; // distances from the seed VertexAttribute& pathOrigins; // previous vertex on the shortest path from origin VertexAttribute& areaElts; // area element attached to each vertex @@ -80,13 +82,14 @@ public : VertexAttribute& a); ~CentroidalVoronoiDiagram (); - void cumulateDistancesOnPaths(); + void cumulateEnergyOnPaths(); unsigned int moveSeeds(); // returns the number of seeds that did move + REAL getGlobalEnergy() {return globalEnergy;} protected : void clear(); void collectVertexFromFront(Dart e); - REAL cumulateDistancesFromRoot(Dart e); + REAL cumulateEnergyFromRoot(Dart e); unsigned int moveSeed(unsigned int numSeed); }; diff --git a/include/Algo/Geometry/voronoiDiagrams.hpp b/include/Algo/Geometry/voronoiDiagrams.hpp index aad48229e9c190e3d15f7dfcd7955e10980d67ea..e68108a96a08050d9939fc6e7ba69f87acb1840b 100644 --- a/include/Algo/Geometry/voronoiDiagrams.hpp +++ b/include/Algo/Geometry/voronoiDiagrams.hpp @@ -172,39 +172,43 @@ void CentroidalVoronoiDiagram::collectVertexFromFront(Dart e){ } template -void CentroidalVoronoiDiagram::cumulateDistancesOnPaths(){ +void CentroidalVoronoiDiagram::cumulateEnergyOnPaths(){ + globalEnergy = 0.0; for (unsigned int i = 0; i < this->seeds.size(); i++) { - cumulateDistancesFromRoot(this->seeds[i]); + cumulateEnergyFromRoot(this->seeds[i]); + globalEnergy += distances[this->seeds[i]]; } } template unsigned int CentroidalVoronoiDiagram::moveSeeds(){ unsigned int m = 0; + globalEnergy = 0.0; for (unsigned int i = 0; i < this->seeds.size(); i++) { m += moveSeed(i); + globalEnergy += distances[this->seeds[i]]; } return m; } template -typename PFP::REAL CentroidalVoronoiDiagram::cumulateDistancesFromRoot(Dart e){ - REAL area = areaElts[e]; - distances[e] *= areaElts[e]; +typename PFP::REAL CentroidalVoronoiDiagram::cumulateEnergyFromRoot(Dart e){ + REAL distArea = areaElts[e] * distances[e]; + distances[e] = areaElts[e] * distances[e] * distances[e]; Traversor2VVaE tv (this->map, e); for (Dart f = tv.begin(); f != tv.end(); f=tv.next()) { if ( pathOrigins[f] == this->map.phi2(f)) { - area += cumulateDistancesFromRoot(f); + distArea += cumulateEnergyFromRoot(f); distances[e] += distances[f]; } } - return area; + return distArea; } template @@ -215,10 +219,9 @@ unsigned int CentroidalVoronoiDiagram::moveSeed(unsigned int numSeed){ std::vector v; v.reserve(8); - std::vector a; - a.reserve(8); + std::vector da; + da.reserve(8); - float regionArea = areaElts[e]; distances[e] = 0.0; Traversor2VVaE tv (this->map, e); @@ -226,46 +229,51 @@ unsigned int CentroidalVoronoiDiagram::moveSeed(unsigned int numSeed){ { if ( pathOrigins[f] == this->map.phi2(f)) { - float area = cumulateDistancesFromRoot(f); + float distArea = cumulateEnergyFromRoot(f); + da.push_back(distArea); distances[e] += distances[f]; v.push_back(f); - a.push_back(area); - regionArea += area; } } /* second attempt */ + // TODO : check if the computation of grad and proj is still valid for other edgeCost than geodesic distances VEC3 grad (0.0); const VertexAttribute& pos = this->map.template getAttribute("position"); + // compute the gradient for (unsigned int j = 0; j(this->map,this->map.phi2(f),pos); VEC3 edgeV = pos[f] - pos[this->map.phi2(f)]; edgeV.normalize(); - grad += distances[f] * edgeV; + grad += da[j] * edgeV; } + grad /= 2.0; float maxProj = 0.0; +// float memoForTest = 0.0; for (unsigned int j = 0; j(this->map,this->map.phi2(f),pos); VEC3 edgeV = pos[f] - pos[this->map.phi2(f)]; // edgeV.normalize(); float proj = edgeV * grad; - proj -= areaElts[e] * this->edgeCost[f] * this->edgeCost[f]; +// proj -= areaElts[e] * this->edgeCost[f] * this->edgeCost[f]; if (proj > maxProj) { - if (numSeed==1 && j==0) - CGoGNout << (edgeV * grad) / (areaElts[e] * this->edgeCost[f] * this->edgeCost[f]) * this->seeds.size() << CGoGNendl; +// 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]; } } + +// if (numSeed==1 && seedMoved==1) +// CGoGNout << memoForTest * this->seeds.size() << CGoGNendl; + /* end second attempt */ /* first attempt by approximating the energy change