appli/examples/example_ImageAlgorithm_Skeletonization.cxx

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