Voronoi_Rn.cpp

Go to the documentation of this file.

// politopix allows to make computations on polytopes such as finding vertices, intersecting, Minkowski sums, ...
//     Copyright (C) 2011-2022 : Delos Vincent
//
//     This program is free software: you can redistribute it and/or modify
//     it under the terms of the GNU General Public License as published by
//     the Free Software Foundation, either version 3 of the License, or
//     (at your option) any later version.
//
//     This program 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 General Public License for more details.
//
//     You should have received a copy of the GNU General Public License
//     along with this program.  If not, see <http://www.gnu.org/licenses/>.
//
// I2M (UMR CNRS 5295 / University of Bordeaux)
 
#include "GeometricObjectIterator_Rn.h"
#include "PolyhedralAlgorithms_Rn.h"
#include "DoubleDescription_Rn.h"
#include "Voronoi_Rn.h"
#include "IO_Polytope.h"
#include <math.h>
#include <numeric>
#include <sstream>
#include <iostream>
#include <fstream>
#include <sstream>
 
 
// Constructor
Voronoi_Rn::Voronoi_Rn(const boost::shared_ptr<Polytope_Rn>& inputSpace, 
                       const std::vector<Point_Rn>& listOfSeeds):_inputSpace(inputSpace),_listOfSeeds(listOfSeeds) {
  _listOfVoronoiCells.resize(_listOfSeeds.size());
}
 
boost::shared_ptr<HalfSpace_Rn> Voronoi_Rn::computeMidPlane(std::vector<Point_Rn>::const_iterator iteSeed1, std::vector<Point_Rn>::const_iterator iteSeed2) {
  unsigned int dimension = iteSeed1->dimension();
  boost::shared_ptr<HalfSpace_Rn> HS;
  HS.reset(new HalfSpace_Rn(dimension));
  double sum_square = 0.;
  for (unsigned int coord_count=0; coord_count<dimension; coord_count++) {
    HS->setCoefficient(coord_count, iteSeed1->getCoordinate(coord_count)-iteSeed2->getCoordinate(coord_count));
    double sq = iteSeed1->getCoordinate(coord_count)*iteSeed1->getCoordinate(coord_count);
    sq -= iteSeed2->getCoordinate(coord_count)*iteSeed2->getCoordinate(coord_count);
    sq /= 2.;
    sum_square += sq;
  }
  HS->setConstant(-sum_square);
  return HS;
}
 
// CHECK POLYHEDRON
bool Voronoi_Rn::checkTopologyAndGeometry() const {
  return true;
}
 
// Dump the cell structure on the given output.
void Voronoi_Rn::dump(std::ostream& out) const {
  out << std::endl << "#VORONOI DIAGRAM" << std::endl;
  unsigned int nb=0;
  std::vector< boost::shared_ptr<Polytope_Rn> >::const_iterator itePC=_listOfVoronoiCells.begin();
  {for (; itePC!=_listOfVoronoiCells.end(); ++itePC) {
    out << "# Cell "<< nb << std::endl;
    (*itePC)->dump(out);
    out << std::endl;
#ifdef DEBUG
    //std::ostringstream stream_C;
    //stream_C << "cell";
    //stream_C << nb;
    //stream_C << ".ptop";
    //std::string cellfile = stream_C.str();
    //IO_Polytope::save(cellfile, *itePC);
#endif
   ++nb;
  }}
}
 
void Voronoi_Rn::gnuplot(std::ostream& out)  const {
  if (Rn::getDimension() != 2)
    throw std::domain_error("Voronoi_Rn::gnuplot(std::ostream& out) dimension is not 2 ");
  unsigned int nb = 1;
  // To code an unique color for each polygon.
  std::vector< boost::shared_ptr<Polytope_Rn> >::const_iterator iteVC=_listOfVoronoiCells.begin();
  {for (; iteVC!=_listOfVoronoiCells.end(); ++iteVC) {
    double frac = (double)nb / (double)(*iteVC)->numberOfGenerators();
    std::ostringstream stream_;
    stream_ << nb;
    std::string valString = stream_.str();
    Visualization::gnuplot2D(*iteVC, valString, frac, out);
    ++nb;
  }}
}
 
bool IncrementalVoronoi::compute() {
 
  if (_listOfSeeds.empty() == true)
    throw std::length_error("IncrementalVoronoi::compute() needs seeds as input!");
 
  if (_listOfSeeds.size() == 1) {
    // Make a deep copy of the current input space
    boost::shared_ptr<Polytope_Rn> is(new Polytope_Rn( *_inputSpace ));
    _listOfVoronoiCells[0] = is;
    return true;
  }
 
  _listOfVoronoiCells.clear();
  // The first cell is the whole input space.
  boost::shared_ptr<Polytope_Rn> is(new Polytope_Rn( *_inputSpace ));
  _listOfVoronoiCells.push_back(is);
  //unsigned int alreadyProcessedSeeds = 1;
  unsigned int counter = 0;
  {for (std::vector<Point_Rn>::const_iterator newVoronoiSeed = (_listOfSeeds.begin()+1); newVoronoiSeed!=_listOfSeeds.end(); ++newVoronoiSeed) {
 
#ifdef DEBUG
    std::cout << std::endl;
    std::cout << "Seed " << newVoronoiSeed-_listOfSeeds.begin() << ": ";
    (newVoronoiSeed)->save(std::cout);
    std::cout << std::endl;
#endif
    // NEW CELL //
    // Construction of the current cell as a first copy of the input space, then to be chopped.
    boost::shared_ptr<Polytope_Rn> newVoronoiCell(new Polytope_Rn(*_inputSpace));
    unsigned int newTruncationStep = newVoronoiCell->numberOfHalfSpaces();
    {for (std::vector<Point_Rn>::const_iterator processedSeed = _listOfSeeds.begin(); processedSeed!=newVoronoiSeed; ++processedSeed) {
      // Compute all the half-spaces of the current cell associated to newVoronoiSeed.
      boost::shared_ptr<HalfSpace_Rn> HS = computeMidPlane(newVoronoiSeed, processedSeed);
      newVoronoiCell->addHalfSpace(HS);
    }}
    // Now the truncation of the new cell
    constIteratorOfListOfGeometricObjects< boost::shared_ptr<HalfSpace_Rn> > lexmin_ite(newVoronoiCell->getListOfHalfSpaces());
    //StrongRedundancyProcessing< boost::shared_ptr<PolyhedralCone_Rn> > NRP(newVoronoiCell->neigbourhoodCondition());
    DoubleDescription<
      boost::shared_ptr<PolyhedralCone_Rn>,
      constIteratorOfListOfGeometricObjects< boost::shared_ptr<HalfSpace_Rn> > >
      //StrongRedundancyProcessing< boost::shared_ptr<PolyhedralCone_Rn> > >
        DD(newVoronoiCell, lexmin_ite, newTruncationStep);
 
    // OLD CELLS //
    // Chopping of the already created cells.
    std::vector<Point_Rn>::const_iterator processedSeed;
    std::vector< boost::shared_ptr<Polytope_Rn> >::iterator oldVoronoiCell;
    {for (oldVoronoiCell = _listOfVoronoiCells.begin(), processedSeed = _listOfSeeds.begin();
        oldVoronoiCell!=_listOfVoronoiCells.end(); ++oldVoronoiCell, ++processedSeed) {
#ifdef DEBUG
        //gnuplot(std::cout);
        //(newVoronoiSeed)->save(std::cout);
        //std::cout << "### Voronoi cell ###";
        //(*oldVoronoiCell)->dump(std::cout);
        //std::cout << std::endl;
#endif
        // Build the half-spaces the other way
        boost::shared_ptr<HalfSpace_Rn> HS = computeMidPlane(processedSeed, newVoronoiSeed);
        unsigned int oldTruncationStep = (*oldVoronoiCell)->numberOfHalfSpaces();
        (*oldVoronoiCell)->addHalfSpace(HS);
        // Now the truncation of the old cells
        constIteratorOfListOfGeometricObjects< boost::shared_ptr<HalfSpace_Rn> > lexmin_ite2((*oldVoronoiCell)->getListOfHalfSpaces());
        //StrongRedundancyProcessing< boost::shared_ptr<PolyhedralCone_Rn> > NRP2((*oldVoronoiCell)->neigbourhoodCondition());
        DoubleDescription<
          boost::shared_ptr<PolyhedralCone_Rn>,
          constIteratorOfListOfGeometricObjects< boost::shared_ptr<HalfSpace_Rn> > >
          //StrongRedundancyProcessing< boost::shared_ptr<PolyhedralCone_Rn> > >
            DD2(*oldVoronoiCell, lexmin_ite2, oldTruncationStep);
    }}
    _listOfVoronoiCells.push_back(newVoronoiCell);
    //std::cout << counter << std::endl;
    ++counter;
  }}
#ifdef DEBUG
  //gnuplot(std::cout);
  //(newVoronoiSeed)->save(std::cout);
  //std::cout << "### Voronoi cell ###";
  //(*oldVoronoiCell)->dump(std::cout);
  //std::cout << std::endl;
  //dump(std::cout);
#endif
 
  return true;
}
 
bool CellByCellVoronoi::compute() {
 
  if (_listOfSeeds.empty() == true)
    throw std::length_error("CellByCellVoronoi::compute() needs seeds as input!");
 
  if (_listOfSeeds.size() == 1) {
    // Make a deep copy of the current input space
    boost::shared_ptr<Polytope_Rn> is(new Polytope_Rn( *_inputSpace ));
    _listOfVoronoiCells[0] = is;
    return true;
  }
 
  //unsigned int alreadyProcessedSeeds = 1;
  unsigned int counter = 0;
  {for (std::vector<Point_Rn>::const_iterator newVoronoiSeed = _listOfSeeds.begin(); newVoronoiSeed!=_listOfSeeds.end(); ++newVoronoiSeed) {
 
#ifdef DEBUG
    std::cout << std::endl;
    std::cout << "Seed " << newVoronoiSeed-_listOfSeeds.begin() << ": ";
    (newVoronoiSeed)->save(std::cout);
    std::cout << std::endl;
#endif
    // NEW CELL //
    // Construction of the current cell as a first copy of the input space, then to be chopped.
    boost::shared_ptr<Polytope_Rn> newVoronoiCell(new Polytope_Rn(*_inputSpace));
    unsigned int newTruncationStep = newVoronoiCell->numberOfHalfSpaces();
    {for (std::vector<Point_Rn>::const_iterator currentSeed = _listOfSeeds.begin(); currentSeed!=_listOfSeeds.end(); ++currentSeed) {
      if (currentSeed != newVoronoiSeed) {
        // Compute all the half-spaces of the current cell associated to newVoronoiSeed.
        boost::shared_ptr<HalfSpace_Rn> HS = computeMidPlane(newVoronoiSeed, currentSeed);
        newVoronoiCell->addHalfSpace(HS);
      }
    }}
    // Now the truncation of the new cell
    constIteratorOfListOfGeometricObjects< boost::shared_ptr<HalfSpace_Rn> > lexmin_ite(newVoronoiCell->getListOfHalfSpaces());
    //StrongRedundancyProcessing< boost::shared_ptr<PolyhedralCone_Rn> > NRP(newVoronoiCell->neigbourhoodCondition());
    DoubleDescription<
      boost::shared_ptr<PolyhedralCone_Rn>,
      constIteratorOfListOfGeometricObjects< boost::shared_ptr<HalfSpace_Rn> > >
      //StrongRedundancyProcessing< boost::shared_ptr<PolyhedralCone_Rn> > >
        DD(newVoronoiCell, lexmin_ite, newTruncationStep);
 
    _listOfVoronoiCells[counter] = newVoronoiCell;
    //std::cout << counter << std::endl;
    ++counter;
  }}
#ifdef DEBUG
  //gnuplot(std::cout);
  //(newVoronoiSeed)->save(std::cout);
  //std::cout << "### Voronoi cell ###";
  //(*oldVoronoiCell)->dump(std::cout);
  //std::cout << std::endl;
  //dump(std::cout);
#endif
 
  return true;
}
 
//bool CellByCellVoronoiNgb::compute() {
 
  //if (_listOfSeeds.empty() == true)
    //throw std::length_error("ParallelVoronoi::compute() needs seeds as input!");
 
  //if (_listOfSeeds.size() == 1) {
    // Make a deep copy of the current input space
    //boost::shared_ptr<Polytope_Rn> is(new Polytope_Rn( *_inputSpace ));
    //_listOfVoronoiCells.push_back(is);
    //return true;
  //}
  //unsigned int RnDIM = Rn::getDimension();
  //std::vector< std::pair<unsigned int,unsigned int> > v(RnDIM);
  //_facetNeighbours.assign(_listOfSeeds.size(), v);
  //_facetNeighbours.reserve(_listOfSeeds.size());
 
  //unsigned int alreadyProcessedSeeds = 1;
  //unsigned int counter = 0;
  //{for (std::vector<Point_Rn>::const_iterator newVoronoiSeed = _listOfSeeds.begin(); newVoronoiSeed!=_listOfSeeds.end(); ++newVoronoiSeed) {
 
//#ifdef DEBUG
    //std::cout << std::endl;
    //std::cout << "Seed " << newVoronoiSeed-_listOfSeeds.begin() << ": ";
    //(newVoronoiSeed)->save(std::cout);
    //std::cout << std::endl;
//#endif
    // NEW CELL //
    // Build the array of all half-spaces.
    //std::vector< boost::shared_ptr<HalfSpace_Rn> > allHalfSpaces;
    // Construction of the current cell as a first copy of the input space, then to be chopped.
    //boost::shared_ptr<Polytope_Rn> newVoronoiCell(new Polytope_Rn(*_inputSpace));
    //unsigned int newTruncationStep = newVoronoiCell->numberOfHalfSpaces();
    //{for (std::vector<Point_Rn>::const_iterator currentSeed=_listOfSeeds.begin(); currentSeed!=_listOfSeeds.end(); ++currentSeed) {
      //if (currentSeed != newVoronoiSeed) {
        // Compute all the half-spaces of the current cell associated to newVoronoiSeed.
        //boost::shared_ptr<HalfSpace_Rn> HS = computeMidPlane(newVoronoiSeed, currentSeed);
        //newVoronoiCell->addHalfSpace(HS);
        //allHalfSpaces.push_back(HS);
      //}
      //else {
        // Build an artificial half-space that will always be respected: 1 + 0.x1 + ... + 0.xn >= 0.
        // This is not to change the numbering and should have no influence on the process.
        //boost::shared_ptr<HalfSpace_Rn> HS;
        //HS.reset(new HalfSpace_Rn(RnDIM));
        //HS->setConstant(1.);
        // The normal vector
        //for (unsigned int coord_count=0; coord_count<RnDIM; coord_count++)
          //HS->setCoefficient(coord_count, 0.);
        //allHalfSpaces.push_back(HS);
      //}
    //}}
    // Before we truncate the cell we want to swap half-spaces: put first the half-spaces that are known to give genuine facets.
    //if (!_facetNeighbours[counter].empty()) {
      // Half-spaces numbered from 0 to (newTruncationStep-1) belong to the input space only
      //std::vector< boost::shared_ptr<HalfSpace_Rn> >::iterator hs2translate = allHalfSpaces.begin()+newTruncationStep;
      //for (unsigned int i=0; i<_facetNeighbours[counter].size(); ++i) {
        // Reminder _facetNeighbours[counter][i].first refers to a half-space number and _facetNeighbours[counter][i].second to a neighbour polytope number
        //std::vector< boost::shared_ptr<HalfSpace_Rn> >::iterator hs_ngb = allHalfSpaces.begin() + _facetNeighbours[counter][i].first;
        //iter_swap(hs2translate, hs_ngb);
        //++hs2translate;
      //}
    //}
    // Now the truncation of the new cell
    //constIteratorOfListOfGeometricObjects< boost::shared_ptr<HalfSpace_Rn> > lexmin_ite(newVoronoiCell->getListOfHalfSpaces());
    //StrongRedundancyProcessing< boost::shared_ptr<PolyhedralCone_Rn> > NRP;
    //DoubleDescription<
      //boost::shared_ptr<PolyhedralCone_Rn>,
      //constIteratorOfListOfGeometricObjects< boost::shared_ptr<HalfSpace_Rn> >,
      //StrongRedundancyProcessing< boost::shared_ptr<PolyhedralCone_Rn> > >
        //DD(newVoronoiCell, lexmin_ite, NRP, newTruncationStep);
 
    //_listOfVoronoiCells.push_back(newVoronoiCell);
 
    //std::cout << counter << std::endl;
    //std::ostringstream stream_;
    //stream_ << "./ptop/cellVor";
    //stream_ << counter;
    //stream_ << ".ptop";
    //std::string valString = stream_.str();
    //std::cout << valString << std::endl;
    //IO_Polytope::save(valString, newVoronoiCell);
 
    // Compute all the facet neighbors of newVoronoiCell.
    // Don't forget that the first newTruncationStep half-spaces belong to the input space.
    // Half-spaces numbered from 0 to (newTruncationStep-1) belong to the input space only
    // and are not relevant from the neighborhood point of view.
    //updateFacetNeighbours(allHalfSpaces, newVoronoiCell, counter, newTruncationStep);
 
    //++counter;
 
  //}}
//#ifdef DEBUG
  //gnuplot(std::cout);
  //(newVoronoiSeed)->save(std::cout);
  //std::cout << "### Voronoi cell ###";
  //(*oldVoronoiCell)->dump(std::cout);
  //std::cout << std::endl;
  //dump(std::cout);
//#endif
 
  //return true;
//}
 
//void CellByCellVoronoiNgb::updateFacetNeighbours(
      //const std::vector< boost::shared_ptr<HalfSpace_Rn> >& allHalfSpaces,
      //boost::shared_ptr<Polytope_Rn> newVoronoiCell,
      //unsigned int counter,
      //unsigned int newTruncationStep) {
  //constIteratorOfListOfGeometricObjects< boost::shared_ptr<HalfSpace_Rn> > iteHS_Vorcell(newVoronoiCell->getListOfHalfSpaces());
  //unsigned int hs_vorcell_counter = 0;
  //unsigned int hs_all_counter = 0;
  //{for (iteHS_Vorcell.begin(); iteHS_Vorcell.end()!=true; iteHS_Vorcell.next()) {
    //const boost::shared_ptr<HalfSpace_Rn>& currentHalfSpace = iteHS_Vorcell.current();
    //std::vector< boost::shared_ptr<HalfSpace_Rn> >::const_iterator iteAll;
    //{for (iteAll=allHalfSpaces.begin(); iteAll!=allHalfSpaces.end(); ++iteAll) {
      //if (iteHS_Vorcell.current() == *iteAll)
        //if (hs_all_counter >= newTruncationStep) {
          //std::pair< unsigned int, unsigned int > curHalfspace_ngbPolytope(hs_vorcell_counter, hs_all_counter);
          // counter has already been incremented here so decrease it by one.
          //_facetNeighbours[counter].push_back(curHalfspace_ngbPolytope);//-newTruncationStep);
        //}
        //++hs_all_counter;
    //}}
    //++hs_vorcell_counter;
  //}}
//}
 
struct sortOnFirstCoord {
    bool operator() (std::pair<double, unsigned int> pt1, std::pair<double, unsigned int> pt2) { return (pt1.first < pt2.first);}
} thisSortOnFirstCoord;
 
void CellByCellVoronoiDist::processDistances(
    std::vector< std::pair< double, unsigned int > >& allDistances,
    std::vector<Point_Rn>::const_iterator newVoronoiSeed,
    boost::shared_ptr<Polytope_Rn> newVoronoiCell)
{
  // Determine all distances from currentSeed to the others
  allDistances.clear();
  std::pair< double, unsigned int > v;
  allDistances.assign(_listOfSeeds.size(), v);
 
  unsigned int distCounter = 0;
  {for (std::vector<Point_Rn>::const_iterator currentSeed = _listOfSeeds.begin(); currentSeed!=_listOfSeeds.end(); ++currentSeed) {
    if (currentSeed != newVoronoiSeed) {
      boost::numeric::ublas::vector<double> v1 = currentSeed->vect();
      v1 -= newVoronoiSeed->vect();
      allDistances[distCounter].first = norm_2(v1);          // The distance from newVoronoiSeed to currentSeed
    }
    else {
      allDistances[distCounter].first = 0.;          // The distance from newVoronoiSeed to itself
    }
    // We number the seeds because the array is going to be shuffled.
    allDistances[distCounter].second = distCounter;
    ++distCounter;
  }}
 
  // Sort allDistances using thisSortOnFirstCoord as comparator
  if (_percentageOfSortedHalfSpaces == 100.)
    std::sort(allDistances.begin(), allDistances.end(), thisSortOnFirstCoord);
  else {
    unsigned int nbOfSortedSeeds = _listOfSeeds.size() * (_percentageOfSortedHalfSpaces/100.);
    std::partial_sort(allDistances.begin(), allDistances.begin()+nbOfSortedSeeds, allDistances.end(), thisSortOnFirstCoord);
  }
 
  //std::cout << std::endl;
  //{for (unsigned int i=1; i<allDistances.size(); ++i) {
    //std::cout << allDistances[i][0] << " " << allDistances[i][1] << std::endl;
  //}}
  //std::cout << std::endl;
 
  // Put inside _allHalfSpaces all the half-space that could play a role in shaping newVoronoiSeed
  _allHalfSpaces.clear();
  constIteratorOfListOfGeometricObjects< boost::shared_ptr<HalfSpace_Rn> > iteHSnewVoronoiSeed(newVoronoiCell->getListOfHalfSpaces());
  for (iteHSnewVoronoiSeed.begin(); iteHSnewVoronoiSeed.end()!=true; iteHSnewVoronoiSeed.next())
    _allHalfSpaces.push_back(iteHSnewVoronoiSeed.current());
  {for (unsigned int i=1; i<allDistances.size(); ++i) {
    // i=0 <=> newVoronoiSeed
    std::vector<Point_Rn>::const_iterator currentSeed = _listOfSeeds.begin()+allDistances[i].second;
    // Compute all the half-spaces of the current cell associated to newVoronoiSeed.
    boost::shared_ptr<HalfSpace_Rn> HS = computeMidPlane(newVoronoiSeed, currentSeed);
    newVoronoiCell->addHalfSpace(HS);
    _allHalfSpaces.push_back(HS);
  }}
}
 
bool CellByCellVoronoiDist::compute() {
 
  if (_listOfSeeds.empty() == true)
    throw std::length_error("CellByCellVoronoiDist::compute() needs seeds as input!");
 
  if (_listOfSeeds.size() == 1) {
    // Make a deep copy of the current input space
    boost::shared_ptr<Polytope_Rn> is(new Polytope_Rn( *_inputSpace ));
    _listOfVoronoiCells[0] = is;
    return true;
  }
 
  //unsigned int alreadyProcessedSeeds = 1;
  unsigned int counter = 0;
  {for (std::vector<Point_Rn>::const_iterator newVoronoiSeed = _listOfSeeds.begin(); newVoronoiSeed!=_listOfSeeds.end(); ++newVoronoiSeed) {
 
#ifdef DEBUG
    std::cout << std::endl;
    std::cout << "Seed " << newVoronoiSeed-_listOfSeeds.begin() << ": ";
    (newVoronoiSeed)->save(std::cout);
    std::cout << std::endl;
#endif
    // NEW CELL //
    // Construction of the current cell as a first copy of the input space, then to be chopped.
    boost::shared_ptr<Polytope_Rn> newVoronoiCell(new Polytope_Rn(*_inputSpace));
    unsigned int newTruncationStep = newVoronoiCell->numberOfHalfSpaces();
 
    // Order the half-spaces inside newVoronoiCell according to sorting criteria based on distances.
    std::vector< std::pair< double, unsigned int > > allDistances;
    processDistances(allDistances, newVoronoiSeed, newVoronoiCell);
 
    // Now the truncation of the new cell
    constIteratorOfListOfGeometricObjects< boost::shared_ptr<HalfSpace_Rn> > lexmin_ite(newVoronoiCell->getListOfHalfSpaces());
    //StrongRedundancyProcessing< boost::shared_ptr<PolyhedralCone_Rn> > NRP(newVoronoiCell->neigbourhoodCondition());
    DoubleDescription<
      boost::shared_ptr<PolyhedralCone_Rn>,
      constIteratorOfListOfGeometricObjects< boost::shared_ptr<HalfSpace_Rn> > >
      //StrongRedundancyProcessing< boost::shared_ptr<PolyhedralCone_Rn> > >
        DD(newVoronoiCell, lexmin_ite, newTruncationStep);
 
    _listOfVoronoiCells[counter] = newVoronoiCell;
    //std::cout << counter << std::endl;
 
    ++counter;
  }}
#ifdef DEBUG
  //gnuplot(std::cout);
  //(newVoronoiSeed)->save(std::cout);
  //std::cout << "### Voronoi cell ###";
  //(*oldVoronoiCell)->dump(std::cout);
  //std::cout << std::endl;
  //dump(std::cout);
#endif
 
  return true;
}
 
bool CellByCellVoronoiDistNgb::compute() {
  std::vector< Point_Rn >::const_iterator iteBeg = _listOfSeeds.begin();
  std::vector< Point_Rn >::const_iterator iteEnd =  _listOfSeeds.end();
  return CellByCellVoronoiDistNgb::compute(iteBeg, iteEnd);
}
 
bool CellByCellVoronoiDistNgb::compute(
    std::vector< Point_Rn >::const_iterator iteBeg,
    std::vector< Point_Rn >::const_iterator iteEnd) {
  
  if (_listOfSeeds.empty() == true)
    throw std::length_error("CellByCellVoronoiDistNgb::compute() needs seeds as input!");
 
  if (_listOfSeeds.size() == 1) {
    // Make a deep copy of the current input space
    boost::shared_ptr<Polytope_Rn> is(new Polytope_Rn( *_inputSpace ));
    _listOfVoronoiCells[0] = is;
    return true;
  }
 
  allocateFacetNeighbours();
 
  //unsigned int alreadyProcessedSeeds = 1;
  unsigned int counter = iteBeg - _listOfSeeds.begin();
  {for (std::vector<Point_Rn>::const_iterator newVoronoiSeed = iteBeg; newVoronoiSeed!=iteEnd; ++newVoronoiSeed) {
 
#ifdef DEBUG
    std::cout << std::endl;
    std::cout << "Seed " << newVoronoiSeed-_listOfSeeds.begin() << ": ";
    (newVoronoiSeed)->save(std::cout);
    std::cout << std::endl;
#endif
    // NEW CELL //
    // Construction of the current cell as a first copy of the input space, then to be chopped.
    boost::shared_ptr<Polytope_Rn> newVoronoiCell(new Polytope_Rn(*_inputSpace));
    unsigned int newTruncationStep = newVoronoiCell->numberOfHalfSpaces();
 
    // Order the half-spaces inside newVoronoiCell according to sorting criteria based on distances.
    std::vector< std::pair< double, unsigned int > > allDistances;
    processDistances(allDistances, newVoronoiSeed, newVoronoiCell);
 
    // Now the truncation of the new cell
    constIteratorOfListOfGeometricObjects< boost::shared_ptr<HalfSpace_Rn> > lexmin_ite(newVoronoiCell->getListOfHalfSpaces());
    //StrongRedundancyProcessing< boost::shared_ptr<PolyhedralCone_Rn> > NRP(newVoronoiCell->neigbourhoodCondition());
    DoubleDescription<
      boost::shared_ptr<PolyhedralCone_Rn>,
      constIteratorOfListOfGeometricObjects< boost::shared_ptr<HalfSpace_Rn> > >
      //StrongRedundancyProcessing< boost::shared_ptr<PolyhedralCone_Rn> > >
        DD(newVoronoiCell, lexmin_ite, newTruncationStep);
 
    _listOfVoronoiCells[counter] = newVoronoiCell;
    //if (counter%100 == 0)
      //std::cout << counter << std::endl;
    // To stop the computations before the end.
    //if (counter == 500)
      //return true;
 
    updateFacetNeighbours(allDistances, newVoronoiCell, counter, newTruncationStep);
 
    ++counter;
  }}
#ifdef DEBUG
  //gnuplot(std::cout);
  //(newVoronoiSeed)->save(std::cout);
  //std::cout << "### Voronoi cell ###";
  //(*oldVoronoiCell)->dump(std::cout);
  //std::cout << std::endl;
  //dump(std::cout);
#endif
 
  return true;
}
 
void CellByCellVoronoiDistNgb::dumpFacetNeighbours(std::ostream &this_ostream) const {
  std::vector< std::vector< std::pair<boost::shared_ptr<HalfSpace_Rn>, unsigned int> > >::const_iterator iteSeed;
  // The list of seeds and their associated cells.
  unsigned int cellCounter = 0;
  for (iteSeed=_facetNeighbours.begin(); iteSeed!=_facetNeighbours.end(); ++iteSeed) {
    this_ostream << iteSeed->size() << " ";
    std::vector< std::pair<boost::shared_ptr<HalfSpace_Rn>, unsigned int> >::const_iterator iteCurHS_NgbNb;
    for (iteCurHS_NgbNb=iteSeed->begin(); iteCurHS_NgbNb!=iteSeed->end(); ++iteCurHS_NgbNb) {
      this_ostream << _listOfVoronoiCells[cellCounter]->getHalfSpaceNumber(iteCurHS_NgbNb->first) << " " << iteCurHS_NgbNb->second << " ";
    }
    this_ostream << std::endl;
    ++cellCounter;
  }
}
 
void CellByCellVoronoiDistNgb::updateFacetNeighbours(
      const std::vector< std::pair< double, unsigned int > >& allDistances,
      boost::shared_ptr<Polytope_Rn> newVoronoiCell,
      unsigned int counter, unsigned int newTruncationStep) {
  constIteratorOfListOfGeometricObjects< boost::shared_ptr<HalfSpace_Rn> > iteHS_Vorcell(newVoronoiCell->getListOfHalfSpaces());
  for (iteHS_Vorcell.begin(); iteHS_Vorcell.end()!=true; iteHS_Vorcell.next()) {
    //const boost::shared_ptr<HalfSpace_Rn>& currentHalfSpace = iteHS_Vorcell.current();
    std::vector< boost::shared_ptr<HalfSpace_Rn> >::const_iterator iteAll;
    unsigned int hs_all_counter = newTruncationStep;
    for (iteAll=_allHalfSpaces.begin()+newTruncationStep; iteAll!=_allHalfSpaces.end(); ++iteAll) {
      if (iteHS_Vorcell.current() == *iteAll) {
        // allHalfSpaces has been shuffled by the sorting process so we lost the information about which seeds are connected
        // to newVoronoiCell. However we can retrieve it in allDistances[hs_all_counter].second
        // Add -newTruncationStep+1 to hs_all_counter because there a shift between allDistances and _allHalfSpaces
        // _allHalfSpaces contains first all the input space half-spaces.
        // Example with a cube as the input space:
        // _allHalfSpaces: H0  ...  H7           H8   ... Hn
        // allDistances:        (0.,counter)  (d1,n1)
        unsigned int ngbSeedNumber = allDistances[hs_all_counter-newTruncationStep+1].second;
        std::pair< boost::shared_ptr<HalfSpace_Rn>, unsigned int > curHalfspace_ngbPolytope(iteHS_Vorcell.current(), ngbSeedNumber);
        _facetNeighbours[counter].push_back(curHalfspace_ngbPolytope);//-newTruncationStep);
        break;
      }
      ++hs_all_counter;
    }
  }
}
 
//========================================================================================================================
 
// void VoronoiWithProperties::loadProp(const std::string& filename, std::vector< int >& propList) {
//   std::ifstream file(filename.c_str(), std::ifstream::in);
//   if (!file) {
//     std::string s("VoronoiWithProperties::loadProp() Unable to open ");
//     s += filename;
//     s += "\n";
//     throw std::ios_base::failure(s);
//   }
 
//   propList.clear();
//   int prop;
//   unsigned int numberOfSeeds, counter=0;
//   std::string line;
//   std::getline(file, line);
//   std::istringstream iline(line);
//   iline >> numberOfSeeds;
//   while (counter != numberOfSeeds) {
//     std::getline(file, line);
//     std::istringstream iline2(line);
//     iline2 >> prop;
//     propList.push_back(prop);
//     ++counter;
//   }
//   file.close();
// }
 
// void VoronoiWithProperties::loadProperties(const std::string& filename) {
//   std::ifstream file(filename.c_str(), std::ifstream::in);
//   if (!file) {
//     std::string s("VoronoiWithProperties::loadProperties() Unable to open ");
//     s += filename;
//     s += "\n";
//     throw std::ios_base::failure(s);
//   }
 
//   _listOfSeedProperties.clear();
//   int prop;
//   unsigned int numberOfSeeds, counter=0;
//   std::string line;
//   std::getline(file, line);
//   std::istringstream iline(line);
//   iline >> numberOfSeeds;
//   while (counter != numberOfSeeds) {
//     std::getline(file, line);
//     std::istringstream iline2(line);
//     iline2 >> prop;
//     _listOfSeedProperties.push_back(prop);
//     ++counter;
//   }
//   file.close();
// }
 
// void VoronoiWithProperties::loadSeeds(const std::string& filename) {
//   boost::shared_ptr<Polytope_Rn> listOfPointsAsVDesc(new Polytope_Rn());
//   IO_Polytope::load(filename, listOfPointsAsVDesc);
//   constIteratorOfListOfGeometricObjects< boost::shared_ptr<Generator_Rn> > iteGN(listOfPointsAsVDesc->getListOfGenerators());
//   {for (iteGN.begin(); iteGN.end()!=true; iteGN.next()) {
//     _listOfSeeds.push_back(Point_Rn(
//           iteGN.current()->getCoordinate(0),
//           iteGN.current()->getCoordinate(1),
//           iteGN.current()->getCoordinate(2)) );
//   }}
// }
 
// void VoronoiWithProperties::loadCell(const std::string& filename) {
//   boost::shared_ptr<Polytope_Rn> thisCell(new Polytope_Rn());
//   IO_Polytope::load(filename, thisCell);
//   _listOfVoronoiCells.push_back(thisCell);
// }
 
// void VoronoiWithProperties::loadFacetNeighbours(const std::string& filename) {
//   std::ifstream file(filename.c_str(), std::ifstream::in);
//   if (!file) {
//     std::string s("VoronoiWithProperties::loadFacetNeighbours() Unable to open ");
//     s += filename;
//     s += "\n";
//     throw std::ios_base::failure(s);
//   }
//   if (_listOfVoronoiCells.empty() == true)
//     throw std::length_error("VoronoiWithProperties::loadFacetNeighbours() needs polytopes!");
 
//   allocateFacetNeighbours();
//   unsigned int numberOfSeeds, cellCounter = 0;
//   std::string line;
//   std::getline(file, line);
//   std::istringstream iline(line);
//   iline >> numberOfSeeds;
//   //std::vector< std::vector< std::pair<boost::shared_ptr<HalfSpace_Rn>, unsigned int> > >::iterator iteSeed;
//   while (cellCounter != numberOfSeeds) {
//     unsigned int nbnbNgb, nbFacet, nbNgb;
//     std::getline(file, line);
//     std::istringstream iline2(line);
//     iline2 >> nbnbNgb;
//     for (unsigned int coord_count=0; coord_count<nbnbNgb; coord_count++) {
//       iline2 >> nbFacet;
//       iline2 >> nbNgb;
//       addHalfSpaceAndNeighbour(cellCounter, _listOfVoronoiCells[cellCounter]->getHalfSpace(nbFacet), nbNgb);
//     }
//     ++cellCounter;
//   }
//   file.close();
// }
 
void VoronoiWithProperties::dumpNeighbours(std::ostream &this_ostream) const {
  std::vector< std::vector< std::pair<boost::shared_ptr<HalfSpace_Rn>, unsigned int> > >::const_iterator iteSeed;
  // The list of seeds and their associated cells.
  unsigned int cellCounter = 0;
  for (iteSeed=_facetNeighbours.begin(); iteSeed!=_facetNeighbours.end(); ++iteSeed) {
    this_ostream << iteSeed-_facetNeighbours.begin() << ": ";
    std::vector< std::pair<boost::shared_ptr<HalfSpace_Rn>, unsigned int> >::const_iterator iteCurHS_NgbNb;
    for (iteCurHS_NgbNb=iteSeed->begin(); iteCurHS_NgbNb!=iteSeed->end(); ++iteCurHS_NgbNb) {
      this_ostream << iteCurHS_NgbNb->second << " ";
    }
    this_ostream << std::endl;
    ++cellCounter;
  }
}
 
void VoronoiWithProperties::dumpFacetNeighbours(std::ostream &this_ostream) const {
  std::vector< std::vector< std::pair<boost::shared_ptr<HalfSpace_Rn>, unsigned int> > >::const_iterator iteSeed;
  // The list of seeds and their associated cells.
  unsigned int cellCounter = 0;
  for (iteSeed=_facetNeighbours.begin(); iteSeed!=_facetNeighbours.end(); ++iteSeed) {
    this_ostream << iteSeed->size() << " ";
    std::vector< std::pair<boost::shared_ptr<HalfSpace_Rn>, unsigned int> >::const_iterator iteCurHS_NgbNb;
    for (iteCurHS_NgbNb=iteSeed->begin(); iteCurHS_NgbNb!=iteSeed->end(); ++iteCurHS_NgbNb) {
      this_ostream << _listOfVoronoiCells[cellCounter]->getHalfSpaceNumber(iteCurHS_NgbNb->first) << " " << iteCurHS_NgbNb->second << " ";
    }
    this_ostream << std::endl;
    ++cellCounter;
  }
}
 
bool VoronoiWithProperties::fuseCellsWithSameProperty(std::vector< std::vector< unsigned int > >& newListOfNonConvexPolytopes) {
  if (_listOfSeedProperties.empty() == true)
    throw std::length_error("VoronoiWithProperties::fuseCellsWithSameProperty() needs seeds as input!");
  unsigned int currentSeedNb = 0;
  std::vector< boost::shared_ptr<Polytope_Rn> >::iterator iteCell;
  std::vector< std::vector< std::pair<boost::shared_ptr<HalfSpace_Rn>, unsigned int> > >::iterator iteSeed;
  // During the first step we go through the list of seeds and their associated cells to delete
  // useless half-spaces separating two cells sharing the same property.
  for (iteSeed=_facetNeighbours.begin(), iteCell =_listOfVoronoiCells.begin();
      iteSeed!=_facetNeighbours.end() && iteCell!=_listOfVoronoiCells.end(); ++iteSeed, ++iteCell) {
    std::vector< std::pair<boost::shared_ptr<HalfSpace_Rn>, unsigned int> >::iterator iteCurHS_NgbNb;
    for (iteCurHS_NgbNb=iteSeed->begin(); iteCurHS_NgbNb!=iteSeed->end(); ++iteCurHS_NgbNb) {
      // All the neighbours of iteSeed (iteCurHS_NgbNb->second) and the separating hyperplane in iteCell (iteCurHS_NgbNb->first)
      unsigned int ngbNb = iteCurHS_NgbNb->second;
      // Check whether both cells share the same property
      if (_listOfSeedProperties[currentSeedNb] == _listOfSeedProperties[ngbNb]) {
        // Some cells may have been already emptied
        if ((*iteCell)->numberOfHalfSpaces() > 0) {
          // In that case we need to remove in the current cell the frontier half-space shared with the neighbour cell.
          boost::shared_ptr<HalfSpace_Rn> HS = iteCurHS_NgbNb->first;
          // Remove HS from the current polytope
          (*iteCell)->removeHalfSpaceFromListAndGenerators(HS);
        }
      }
    }
    ++currentSeedNb;
  }
 
  // At this step some cells have been reduced to the null polytope, others have lost facets and only the isolated ones
  //  remain unchanged. Now use a stack to pick up all the neighbours cells having the same property.
  currentSeedNb = 0;
  iteSeed = _facetNeighbours.begin();
  iteCell = _listOfVoronoiCells.begin();
  std::set< unsigned int > allProcessedSeeds;
  newListOfNonConvexPolytopes.clear();
  while (iteSeed!=_facetNeighbours.end() && iteCell!=_listOfVoronoiCells.end()) {
    // Do nothing if the seed and its cell have already been taken into account.
    if (allProcessedSeeds.find(currentSeedNb) == allProcessedSeeds.end()) {
      //allProcessedSeeds.insert(currentSeedNb);// Not needed as allProcessedSeeds.insert(seedNgbNb) will do the job
      std::vector< unsigned int > v;
      newListOfNonConvexPolytopes.push_back(v);
      // Build the stack of neighbours having the same property
      std::stack< unsigned int > allNeighboursWithSameProperty;
      allNeighboursWithSameProperty.push(currentSeedNb);
      unsigned int currentProperty = _listOfSeedProperties[currentSeedNb];
      // Mark it as processed now.
      if (allProcessedSeeds.insert(currentSeedNb).second == true) {
        // Not inserted yet
        if (_listOfVoronoiCells[currentSeedNb]->numberOfHalfSpaces() != 0) {
          // Insert here in the main DS if the associated polytope is not null
          newListOfNonConvexPolytopes[newListOfNonConvexPolytopes.size()-1].push_back(currentSeedNb);
        }
      }
      while (!allNeighboursWithSameProperty.empty()) {
        unsigned int topSeed = allNeighboursWithSameProperty.top();
        allNeighboursWithSameProperty.pop();
        std::vector< std::pair<boost::shared_ptr<HalfSpace_Rn>, unsigned int> >::iterator iteCurHS_NgbNb;
        {for (iteCurHS_NgbNb=_facetNeighbours[topSeed].begin(); iteCurHS_NgbNb!=_facetNeighbours[topSeed].end(); ++iteCurHS_NgbNb) {
          unsigned int ngbSeed = iteCurHS_NgbNb->second;
          if (allProcessedSeeds.find(ngbSeed) == allProcessedSeeds.end()) {
            if (_listOfSeedProperties[ngbSeed] == (int)currentProperty) {
              allNeighboursWithSameProperty.push(ngbSeed);
              allProcessedSeeds.insert(ngbSeed);
              if (_listOfVoronoiCells[ngbSeed]->numberOfHalfSpaces() != 0) {
                // Insert here in the main DS if the associated polytope is not null
                newListOfNonConvexPolytopes[newListOfNonConvexPolytopes.size()-1].push_back(ngbSeed);
              }
            }
          }
        }}
      } // while (!allNeighboursWithSameProperty.empty()) {
    } // if (allProcessedSeeds.find(currentSeedNb) == allProcessedSeeds.end()) {
    ++currentSeedNb;
    ++iteSeed;
    ++iteCell;
  }
 
  return true;
}
 
bool VoronoiWithProperties::fuseCellsWithSameProperty(std::vector< std::vector< boost::shared_ptr<Polytope_Rn> > >& newListOfNonConvexPolytopes) {
  if (_listOfSeedProperties.empty() == true)
    throw std::length_error("VoronoiWithProperties::fuseCellsWithSameProperty() needs seeds as input!");
  std::vector< std::vector< unsigned int > > allNonConvexPolytopeIds;
  // Compute the ids of the polytopes to be fused.
  fuseCellsWithSameProperty(allNonConvexPolytopeIds);
  newListOfNonConvexPolytopes.resize(allNonConvexPolytopeIds.size());
  std::vector< std::vector< unsigned int > >::const_iterator iteNCPIds = allNonConvexPolytopeIds.begin();
  unsigned int currentSeedNb = 0;
  for (; iteNCPIds != allNonConvexPolytopeIds.end(); ++iteNCPIds) {
    std::vector< unsigned int >::const_iterator iteCurrentNCPId = iteNCPIds->begin();
    for (; iteCurrentNCPId != iteNCPIds->end(); ++iteCurrentNCPId) {
      newListOfNonConvexPolytopes[currentSeedNb].push_back(_listOfVoronoiCells[*iteCurrentNCPId]);
    }
    ++currentSeedNb;
  }
 
  return true;
}
 
void VoronoiWithProperties::propagate(std::vector< bool >& listOfAlreadyProcessedSeeds, unsigned int nbCell, unsigned int nbProperty) {
  std::set< boost::shared_ptr<HalfSpace_Rn> > alreadyProcessedFacets;
  // First of all, explore the marked facets of the Voronoi cell to see if their neighbours share the same property.
  {
    std::vector< std::pair< boost::shared_ptr<HalfSpace_Rn>, unsigned int > >::iterator iteNgb = _facetNeighbours[nbCell].begin();
    for (; iteNgb != _facetNeighbours[nbCell].end(); ++iteNgb) {
      // Check all neighbours to see if they share the same property.
      unsigned int newProperty = _listOfSeedProperties[iteNgb->second];
      //std::cout << "iteNgb->second = " << iteNgb->second << std::endl;
      //std::cout << "newProperty = " << newProperty << std::endl;
      alreadyProcessedFacets.insert(iteNgb->first);
      if (nbProperty != newProperty) {
        // This facet has to be stored into the non-convex polytope with its vertices, so collect them first.
        std::vector< boost::shared_ptr< Generator_Rn > > listOfVtcesPerFacet;
        constIteratorOfListOfGeometricObjects< boost::shared_ptr<Generator_Rn> > iteGN(_listOfVoronoiCells[nbCell]->getListOfGenerators());
        for (iteGN.begin(); iteGN.end() != true; iteGN.next()) {
          for (unsigned int i = 0; i < iteGN.current()->numberOfFacets(); i++) {
            if (iteGN.current()->getFacet(i) == iteNgb->first)
              listOfVtcesPerFacet.push_back( iteGN.current() );
          }
        } // for (iteGN.begin(); iteGN.end()!=true; iteGN.next()) {
        _listOfNonConvexPolytopes.back()->addFacetAndItsVertices(iteNgb->first, listOfVtcesPerFacet);
      }
      else if (listOfAlreadyProcessedSeeds[iteNgb->second] == false) {
        // Run recursively from the neighbor cell.
        listOfAlreadyProcessedSeeds[iteNgb->second] = true;
        propagate(listOfAlreadyProcessedSeeds, iteNgb->second, nbProperty);
      }
    }
  }
  // The remaining facets are connected to the outside of the input space.
  {
    for (unsigned int i = 0; i < _listOfVoronoiCells[nbCell]->numberOfHalfSpaces(); ++i) {
      if (alreadyProcessedFacets.find(_listOfVoronoiCells[nbCell]->getHalfSpace(i)) == alreadyProcessedFacets.end()) {
        // This facet is a border facet of the Voronoi diagram.
        std::vector< boost::shared_ptr< Generator_Rn > > listOfVtcesPerFacet;
        constIteratorOfListOfGeometricObjects< boost::shared_ptr<Generator_Rn> > iteGN(_listOfVoronoiCells[nbCell]->getListOfGenerators());
        for (iteGN.begin(); iteGN.end() != true; iteGN.next()) {
          for (unsigned int j = 0; j < iteGN.current()->numberOfFacets(); ++j) {
            if (iteGN.current()->getFacet(j) == _listOfVoronoiCells[nbCell]->getHalfSpace(i))
              listOfVtcesPerFacet.push_back( iteGN.current() );
          }
        } // for (iteGN.begin(); iteGN.end()!=true; iteGN.next()) {
        _listOfNonConvexPolytopes.back()->addFacetAndItsVertices(_listOfVoronoiCells[nbCell]->getHalfSpace(i), listOfVtcesPerFacet);
      }
    }
  }
}
 
bool VoronoiWithProperties::fuseCellsWithSameProperty() {
  if (_listOfSeedProperties.empty() == true)
    throw std::length_error("VoronoiWithProperties::fuseCellsWithSameProperty() needs seeds as input!");
  //std::copy(_listOfSeedProperties.begin(), _listOfSeedProperties.end(), std::ostream_iterator<unsigned int>(std::cout, " ") );
  std::set< unsigned int > listOfAlreadyProcessedLayers;
  std::vector< bool > listOfAlreadyProcessedSeeds(_listOfSeedProperties.size(), false);
  {
    std::vector< std::vector< std::pair< boost::shared_ptr<HalfSpace_Rn>, unsigned int > > >::const_iterator iteCell = _facetNeighbours.begin();
    for (; iteCell != _facetNeighbours.end(); ++iteCell) {
      // Get the current cell property
      unsigned int nbCell = iteCell - _facetNeighbours.begin();
      unsigned int nbProperty = _listOfSeedProperties[nbCell];
      if (listOfAlreadyProcessedSeeds[nbCell] == false) {
        std::cout << listOfAlreadyProcessedSeeds[nbCell] << std::endl;
        listOfAlreadyProcessedSeeds[nbCell] = true;
        boost::shared_ptr< NonConvexPolytope > NCP(new NonConvexPolytope());
        _listOfNonConvexPolytopes.push_back(NCP);
        _listOfNonConvexPolytopeProperties.push_back(nbProperty);
        propagate(listOfAlreadyProcessedSeeds, nbCell, nbProperty);
      }
    }
  }
 
  return true;
}