distancetransform3D.cpp

#include <iostream>

#include <QtGui/qapplication.h>

#include "DGtal/base/Common.h"
#include "DGtal/kernel/SpaceND.h"
#include "DGtal/kernel/domains/HyperRectDomain.h"
#include "DGtal/images/ImageSelector.h"
#include "DGtal/io/readers/VolReader.h"
#include "DGtal/io/viewers/Viewer3D.h"

#include "DGtal/io/Color.h"
#include "DGtal/io/colormaps/GradientColorMap.h"
#include "DGtal/io/DrawWithDisplay3DModifier.h"

#include "DGtal/kernel/SpaceND.h"
#include "DGtal/kernel/domains/HyperRectDomain.h"
#include "DGtal/images/ImageSelector.h"

#include "DGtal/geometry/volumes/distance/SeparableMetricHelper.h"
#include "DGtal/geometry/volumes/distance/DistanceTransformation.h"
#include "DGtal/images/imagesSetsUtils/SimpleThresholdForegroundPredicate.h"

#include "ConfigExamples.h"



using namespace std;
using namespace DGtal;



template<typename Image>
void randomSeeds(Image &image, const unsigned int nb, const int value)
{
  typename Image::Point p, low = image.domain().lowerBound();
  typename Image::Vector ext;
  srand ( time(NULL) );

  ext = image.extent();

  for (unsigned int k = 0 ; k < nb; k++)
  {
    for (unsigned int dim = 0; dim < Image::dimension; dim++)
      p[dim] = rand() % (ext[dim]) +  low[dim];
    
    image.setValue(p, value);
  }
}

int main( int argc, char** argv )
{
  
  std::string inputFilename = examplesPath + "samples/Al.100.vol";
  
 //------------

 typedef SpaceND<3> Space4Type;
 typedef HyperRectDomain<Space4Type> TDomain;

 typedef TDomain::Point Point;

  
 QApplication application(argc,argv);
 Viewer3D viewer;
 viewer.setWindowTitle("simpleViewer");
 viewer.show();


 
 //Default image selector = STLVector
 typedef ImageSelector<TDomain, unsigned char>::Type Image;
 Image image = VolReader<Image>::importVol( inputFilename );
 TDomain domain = image.domain();


 Image imageSeeds ( domain);
 for ( Image::Iterator it = imageSeeds.begin(), itend = imageSeeds.end();it != itend; ++it)
   (*it)=1;
 Z3i::Point p0(10,10,10);
 //imageSeeds.setValue(p0, 0 );
 randomSeeds(imageSeeds, 70, 0);


 //Distance transformation computation
 typedef ImageSelector<TDomain, long int>::Type ImageLong;

 typedef SimpleThresholdForegroundPredicate<Image> Predicate;
 Predicate aPredicate(imageSeeds,0);

 typedef  DistanceTransformation<Space4Type,Predicate, 2> DTL2;
 typedef  DistanceTransformation<Space4Type,Predicate, 0> DTLInf;
 typedef  DistanceTransformation<Space4Type,Predicate, 1> DTL1;
 
 DTL2 dtL2(domain, aPredicate);
 DTLInf dtLinf(domain, aPredicate);
 DTL1 dtL1(domain, aPredicate);

 
 DTL1::OutputImage resultL1 = dtL1.compute (  );

 unsigned int min = 0;
 unsigned int max = 0;
 for(DTL1::OutputImage::ConstIterator it = resultL1.begin(), itend=resultL1.end();
     it!=itend;
     ++it)
   {
     if(  (*it) < min )   
       min=(*it);
     if( (*it) > max ) 
       max=(*it);
   }
     
     
  GradientColorMap<long> gradient( 0,30);
  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);  
 

  viewer << SetMode3D( (*(domain.begin())).className(), "Paving" );
  
  for(TDomain::ConstIterator it = domain.begin(), itend=domain.end();
     it!=itend;
     ++it){
   
   unsigned int valDist= resultL1( (*it) );     
   Color c= gradient(valDist);
   
   if(resultL1(*it)<=30 ){
     viewer << CustomColors3D(Color((float)(c.red()), 
            (float)(c.green()),
            (float)(c.blue(),205)), 
            Color((float)(c.red()), 
            (float)(c.green()),
            (float)(c.blue()),205));
     viewer << *it ;
   }     
 }
 
  //viewer << ClippingPlane(1,0,0,-60);
 viewer<< Viewer3D::updateDisplay;
 
 return application.exec();
}
//                                                                           //