3dKSSurfaceExtraction.cpp

#include <iostream>
#include "DGtal/base/Common.h"

#include <QtGui/qapplication.h>
#include "DGtal/base/Common.h"
#include "DGtal/io/readers/VolReader.h"
#include "DGtal/io/viewers/Viewer3D.h"
#include "DGtal/io/DrawWithDisplay3DModifier.h"
#include "DGtal/io/Color.h"

#include "DGtal/images/ImageSelector.h"
#include "DGtal/helpers/StdDefs.h"
#include "ConfigExamples.h"
#include "DGtal/io/Color.h"
#include "DGtal/io/colormaps/GradientColorMap.h"
#include "DGtal/topology/KhalimskySpaceND.h"
#include "DGtal/topology/helpers/Surfaces.h"



using namespace std;
using namespace DGtal;
using namespace Z3i;


int main( int argc, char** argv )
{
  
  Point p1( 0, 0,  0 );
  Point p2( 20, 20, 20 );
  Point c( 10, 10, 10 );
  Domain domain( p1, p2);

  // Generate the digital set from randam seeds and distance threshold.
  DigitalSet diamond_set( domain );
  //srand ( time(NULL) );
  uint nbSeeds = 35;
  vector<Point> vCenters;
  vector<uint> vRad;
  for(uint i=0;i<nbSeeds; i++){
    vCenters.push_back(Point(rand()%p2[0], rand()%p2[1], 
          rand()%p2[2]));
    vRad.push_back(rand()%7);
  }
  for ( Domain::ConstIterator it = domain.begin(); it != domain.end(); ++it ){
    for(unsigned int i=0;i<nbSeeds; i++){
      if ( (*it - vCenters.at(i)  ).norm1() <= vRad.at(i) && domain.isInside(*it) &&
     domain.isInside(*it+Point(1,1,1)) && domain.isInside(*it-Point(1,1,1)) ){ 
  diamond_set.insertNew( *it );
  break;
      }
    }
  }
  
  
  //A KhalimskySpace is constructed from the domain boundary points.
  KSpace K;
  K.init(p1, p2, true);
  
  SurfelAdjacency<3> SAdj( true );
  vector<vector<SCell> > vectConnectedSCell;
  
  
  //Here since the last argument is set to true, the resulting
  //SignedKhalimskySpaceND are signed in order to indicate the direction
  //of exterior. You can also get the SignefKhalimskySpaceND with default
  //sign:

  SetPredicate<DigitalSet> shape_set_predicate( diamond_set );
  Surfaces<KSpace>::extractAllConnectedSCell(vectConnectedSCell,K, SAdj, shape_set_predicate, false);
  
  
  QApplication application(argc,argv);
  Viewer3D viewer;
  viewer.show(); 


  // Each connected compoments are simply displayed with a specific color.
  GradientColorMap<long> gradient(0, (const long)vectConnectedSCell.size());
  gradient.addColor(Color::Red);
  gradient.addColor(Color::Yellow);
  gradient.addColor(Color::Green);
  gradient.addColor(Color::Cyan);
  gradient.addColor(Color::Blue);
  gradient.addColor(Color::Magenta);
  gradient.addColor(Color::Red);  
 
  
  for(uint i=0; i< vectConnectedSCell.size();i++){
    DGtal::Color col= gradient(i);
    viewer << CustomColors3D(Color(250, 0,0), Color(col.red(), 
                                                    col.green(),
                                                    col.blue()));
    
    for(uint j=0; j< vectConnectedSCell.at(i).size();j++){
      viewer << vectConnectedSCell.at(i).at(j);
    }    
  }

  
  viewer << CustomColors3D(Color(250, 0,0),Color(250, 200,200, 200));
  viewer << diamond_set;
  //viewer << ClippingPlane(0,1,0.0,-2);
  viewer << Viewer3D::updateDisplay;
  return application.exec();
}
//                                                                           //