appli/examples/example_ImageAlgorithmRegionGrow_MultipleThresholds.cxx

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