appli/examples/example_ImageAlgorithm_Skeletonization.cxx

   1 #include <ctime>
   2 #include <iostream>
   3 #include <string>
   4
   5 #include <itkDanielssonDistanceMapImageFilter.h>
   6 #include <itkImage.h>
   7 #include <itkImageFileReader.h>
   8 #include <itkImageToVTKImageFilter.h>
   9
  10 #include <fpa/Base/Functors/InvertCostFunction.h>
  11 #include <fpa/Image/Dijkstra.h>
  12 #include <fpa/Image/RegionGrowWithMultipleThresholds.h>
  13 #include <fpa/VTK/ImageMPR.h>
  14 #include <fpa/VTK/Image3DObserver.h>
  15
  16 /*
  17   #include <limits>
  18   #include <set>
  19
  20   #include <itkImageFileWriter.h>
  21
  22
  23   #include <vtkImageActor.h>
  24   #include <vtkImageMarchingCubes.h>
  25   #include <vtkProperty.h>
  26   #include <vtkRenderer.h>
  27   #include <vtkRenderWindow.h>
  28   #include <vtkRenderWindowInteractor.h>
  29   #include <vtkSmartPointer.h>
  30   #include <vtkSphereSource.h>
  31
  32 */
  33
  34 // -------------------------------------------------------------------------
  35 const unsigned int Dim = 3;
  36 typedef short                                TPixel;
  37 typedef double                               TScalar;
  38 typedef itk::Image< TPixel, Dim >            TImage;
  39 typedef itk::Image< TScalar, Dim >           TDistanceMap;
  40
  41 /*
  42   typedef itk::ImageToVTKImageFilter< TImage > TVTKImage;
  43   typedef itk::ImageFileReader< TImage >   TImageReader;
  44   typedef itk::ImageFileWriter< TImage >   TImageWriter;
  45   typedef
  46   fpa::Image::RegionGrowWithMultipleThresholds< TImage >
  47   TSegmentor;
  48   typedef
  49
  50   TObserver;
  51 */
  52
  53 // -------------------------------------------------------------------------
  54 int main( int argc, char* argv[] )
  55 {
  56   if( argc < 6 )
  57   {
  58     std::cerr
  59       << "Usage: " << argv[ 0 ]
  60       << " input_image thr_0 thr_1 step output_image"
  61       << " visual_debug"
  62       << std::endl;
  63     return( 1 );
  64
  65   } // fi
  66   std::string input_image_fn = argv[ 1 ];
  67   TPixel thr_0 = TPixel( std::atof( argv[ 2 ] ) );
  68   TPixel thr_1 = TPixel( std::atof( argv[ 3 ] ) );
  69   unsigned int step = std::atoi( argv[ 4 ] );
  70   std::string output_image_fn = argv[ 5 ];
  71   bool visual_debug = false;
  72   if( argc > 6 )
  73     visual_debug = ( std::atoi( argv[ 6 ] ) == 1 );
  74
  75   // Read image
  76   itk::ImageFileReader< TImage >::Pointer input_image_reader =
  77     itk::ImageFileReader< TImage >::New( );
  78   input_image_reader->SetFileName( input_image_fn );
  79   try
  80   {
  81     input_image_reader->Update( );
  82   }
  83   catch( itk::ExceptionObject& err )
  84   {
  85     std::cerr << "Error caught: " << err << std::endl;
  86     return( 1 );
  87
  88   } // yrt
  89   TImage::ConstPointer input_image = input_image_reader->GetOutput( );
  90
  91   // Show image
  92   itk::ImageToVTKImageFilter< TImage >::Pointer vtk_image =
  93     itk::ImageToVTKImageFilter< TImage >::New( );
  94   vtk_image->SetInput( input_image );
  95   vtk_image->Update( );
  96
  97   fpa::VTK::ImageMPR view;
  98   view.SetBackground( 0.3, 0.2, 0.8 );
  99   view.SetSize( 800, 800 );
 100   view.SetImage( vtk_image->GetOutput( ) );
 101
 102   // Wait for a seed to be given
 103   while( view.GetNumberOfSeeds( ) == 0 )
 104     view.Start( );
 105
 106   // Get seed from user
 107   double seed[ 3 ];
 108   view.GetSeed( 0, seed );
 109   TImage::PointType seed_pnt;
 110   seed_pnt[ 0 ] = seed[ 0 ];
 111   seed_pnt[ 1 ] = seed[ 1 ];
 112   seed_pnt[ 2 ] = seed[ 2 ];
 113   TImage::IndexType seed_idx;
 114   input_image->TransformPhysicalPointToIndex( seed_pnt, seed_idx );
 115
 116   // Segment input image
 117   fpa::Image::RegionGrowWithMultipleThresholds< TImage >::Pointer segmentor =
 118     fpa::Image::RegionGrowWithMultipleThresholds< TImage >::New( );
 119   segmentor->AddThresholds( thr_0, thr_1, step );
 120   segmentor->AddSeed( seed_idx, 0 );
 121   segmentor->SetInput( input_image );
 122   segmentor->SetNeighborhoodOrder( 1 );
 123   segmentor->SetDifferenceThreshold( double( 3 ) );
 124   segmentor->SetInsideValue( TPixel( 0 ) );  // WARNING: This is to optimize
 125   segmentor->SetOutsideValue( TPixel( 1 ) ); // distance map calculation
 126   if( visual_debug )
 127   {
 128     // Configure observer
 129     fpa::VTK::Image3DObserver< fpa::Image::RegionGrowWithMultipleThresholds< TImage >, vtkRenderWindow >::Pointer obs =
 130       fpa::VTK::Image3DObserver< fpa::Image::RegionGrowWithMultipleThresholds< TImage >, vtkRenderWindow >::New( );
 131     obs->SetImage( input_image, view.GetWindow( ) );
 132     segmentor->AddObserver( itk::AnyEvent( ), obs );
 133     segmentor->ThrowEventsOn( );
 134   }
 135   else
 136     segmentor->ThrowEventsOff( );
 137   std::clock_t start = std::clock( );
 138   segmentor->Update( );
 139   std::clock_t end = std::clock( );
 140   double seconds = double( end - start ) / double( CLOCKS_PER_SEC );
 141   std::cout << "Segmentation time = " << seconds << std::endl;
 142
 143   // Extract distance map
 144   itk::DanielssonDistanceMapImageFilter< TImage, TDistanceMap >::Pointer
 145     distanceMap =
 146     itk::DanielssonDistanceMapImageFilter< TImage, TDistanceMap >::New( );
 147   distanceMap->SetInput( segmentor->GetOutput( ) );
 148   distanceMap->InputIsBinaryOn( );
 149   start = std::clock( );
 150   distanceMap->Update( );
 151   end = std::clock( );
 152   seconds = double( end - start ) / double( CLOCKS_PER_SEC );
 153   std::cout << "Distance map time = " << seconds << std::endl;
 154
 155   // Extract paths
 156   fpa::Image::Dijkstra< TDistanceMap, TScalar >::Pointer paths =
 157     fpa::Image::Dijkstra< TDistanceMap, TScalar >::New( );
 158   paths->AddSeed( segmentor->GetSeed( 0 ), TScalar( 0 ) );
 159   paths->SetInput( distanceMap->GetOutput( ) );
 160   paths->SetNeighborhoodOrder( 1 );
 161
 162   // Associate an inversion cost function
 163   fpa::Base::Functors::InvertCostFunction< TScalar >::Pointer
 164     cost_function =
 165     fpa::Base::Functors::InvertCostFunction< TScalar >::New( );
 166   paths->SetCostConversion( cost_function );
 167
 168   if( visual_debug )
 169   {
 170     // Configure observer
 171     fpa::VTK::Image3DObserver< fpa::Image::Dijkstra< TDistanceMap, TScalar >, vtkRenderWindow >::Pointer obs =
 172       fpa::VTK::Image3DObserver< fpa::Image::Dijkstra< TDistanceMap, TScalar >, vtkRenderWindow >::New( );
 173     obs->SetImage( distanceMap->GetOutput( ), view.GetWindow( ) );
 174     paths->AddObserver( itk::AnyEvent( ), obs );
 175     paths->ThrowEventsOn( );
 176   }
 177   else
 178     paths->ThrowEventsOff( );
 179   start = std::clock( );
 180   paths->Update( );
 181   end = std::clock( );
 182   seconds = double( end - start ) / double( CLOCKS_PER_SEC );
 183   std::cout << "Paths extraction time = " << seconds << std::endl;
 184
 185   // Show result
 186   /*
 187     TVTKImage::Pointer output_image_vtk = TVTKImage::New( );
 188     output_image_vtk->SetInput( segmentor->GetOutput( ) );
 189     output_image_vtk->Update( );
 190
 191     vtkSmartPointer< vtkImageMarchingCubes > mc =
 192     vtkSmartPointer< vtkImageMarchingCubes >::New( );
 193     mc->SetInputData( output_image_vtk->GetOutput( ) );
 194     mc->SetValue(
 195     0,
 196     double( segmentor->GetInsideValue( ) ) * double( 0.95 )
 197     );
 198     mc->Update( );
 199
 200     // Let some interaction and close program
 201     view.AddPolyData( mc->GetOutput( ), 0.1, 0.6, 0.8, 0.5 );
 202     view.Start( );
 203
 204     // Write resulting image
 205     TImageWriter::Pointer output_image_writer = TImageWriter::New( );
 206     output_image_writer->SetInput( segmentor->GetOutput( ) );
 207     output_image_writer->SetFileName( output_image_fn );
 208     try
 209     {
 210     output_image_writer->Update( );
 211     }
 212     catch( itk::ExceptionObject& err )
 213     {
 214     std::cerr << "Error caught: " << err << std::endl;
 215     return( 1 );
 216
 217     } // yrt
 218   */
 219   return( 0 );
 220 }
 221
 222 // eof - $RCSfile$