appli/CTBronchi/MoriSegmentation.cxx

   1 #include <chrono>
   2 #include <itkImage.h>
   3 #include <itkImageFileReader.h>
   4 #include <itkImageFileWriter.h>
   5 #include <fpa/Image/Mori.h>
   6
   7 // -------------------------------------------------------------------------
   8 const unsigned int Dim = 3;
   9 typedef short          TPixel;
  10 typedef unsigned short TLabel;
  11
  12 typedef itk::Image< TPixel, Dim > TInputImage;
  13 typedef itk::Image< TLabel, Dim > TLabelImage;
  14 typedef fpa::Image::Mori< TInputImage, TLabelImage > TFilter;
  15
  16 // -------------------------------------------------------------------------
  17 int main( int argc, char* argv[] )
  18 {
  19   // Get arguments
  20   if( argc < 17 )
  21   {
  22     std::cerr
  23       << "Usage: " << argv[ 0 ] << std::endl
  24       << "   input_image output_image output_signal" << std::endl
  25       << "   init_threshold(-1024) end_threshold(0) delta(1)" << std::endl
  26       << "   minimum_threshold(-850)" << std::endl
  27       << "   inside_value(255) outside_value(0)" << std::endl
  28       << "   signal_kernel_size(20) signal_threshold(500) signal_influence(0.5)"
  29       << std::endl
  30       << "   [index/point] seed_x seed_y seed_z"
  31       << std::endl;
  32     return( 1 );
  33
  34   } // fi
  35   std::string input_image_filename = argv[ 1 ];
  36   std::string output_image_filename = argv[ 2 ];
  37   std::string output_signal_filename = argv[ 3 ];
  38   TPixel init_threshold = std::atoi( argv[ 4 ] );
  39   TPixel end_threshold = std::atoi( argv[ 5 ] );
  40   TPixel delta = std::atoi( argv[ 6 ] );
  41   TPixel minimum_threshold = std::atoi( argv[ 7 ] );
  42   TLabel inside_value = std::atoi( argv[ 8 ] );
  43   TLabel outside_value = std::atoi( argv[ 9 ] );
  44   unsigned long signal_kernel_size = std::atoi( argv[ 10 ] );
  45   double signal_threshold = std::atof( argv[ 11 ] );
  46   double signal_influence = std::atof( argv[ 12 ] );
  47   std::string seed_type = argv[ 13 ];
  48   TInputImage::IndexType iseed;
  49   TInputImage::PointType pseed;
  50   for( unsigned int i = 0; i < Dim; ++i )
  51   {
  52     if( seed_type == "index" )
  53       iseed[ i ] = std::atoi( argv[ 14 + i ] );
  54     else
  55       pseed[ i ] = std::atof( argv[ 14 + i ] );
  56
  57   } // rof
  58
  59   // Read image
  60   itk::ImageFileReader< TInputImage >::Pointer input_image_reader =
  61     itk::ImageFileReader< TInputImage >::New( );
  62   input_image_reader->SetFileName( input_image_filename );
  63
  64   // Prepare filter
  65   TFilter::Pointer filter = TFilter::New( );
  66   filter->SetInput( input_image_reader->GetOutput( ) );
  67   if( seed_type == "index" )
  68     filter->SetSeed( iseed );
  69   else
  70     filter->SetSeed( pseed );
  71   filter->SetThresholds( init_threshold, end_threshold, delta );
  72   filter->SetMinimumThreshold( minimum_threshold );
  73   filter->SetInsideValue( inside_value );
  74   filter->SetOutsideValue( outside_value );
  75   filter->SetSignalKernelSize( signal_kernel_size );
  76   filter->SetSignalThreshold( signal_threshold );
  77   filter->SetSignalInfluence( signal_influence );
  78
  79   // Show some information
  80   std::cout << "----------------------------------------------" << std::endl;
  81   std::cout << "Image: " << input_image_filename << std::endl;
  82
  83   // Execute pipeline
  84   std::chrono::time_point< std::chrono::high_resolution_clock > ts, te;
  85   std::chrono::duration< double > tr;
  86   try
  87   {
  88     ts = std::chrono::high_resolution_clock::now( );
  89     input_image_reader->Update( );
  90     te = std::chrono::high_resolution_clock::now( );
  91     tr = te - ts;
  92     std::cout << "Read time: " << tr.count( ) << " s" << std::endl;
  93
  94     ts = std::chrono::high_resolution_clock::now( );
  95     filter->Update( );
  96     te = std::chrono::high_resolution_clock::now( );
  97     tr = te - ts;
  98     std::cout
  99       << "Mori time: " << tr.count( ) << " s" << std::endl
 100       << "Optimum threshold: " << filter->GetOptimumThreshold( ) << std::endl
 101       << "Number of thresholds: "
 102       << filter->GetNumberOfEvaluatedThresholds( ) << std::endl;
 103   }
 104   catch( std::exception& err )
 105   {
 106     std::cerr << "Error caught: " << err.what( ) << std::endl;
 107     return( 1 );
 108
 109   } // yrt
 110
 111   // Save output image
 112   itk::ImageFileWriter< TLabelImage >::Pointer output_image_writer =
 113     itk::ImageFileWriter< TLabelImage >::New( );
 114   output_image_writer->SetInput( filter->GetThresholdedOutput( ) );
 115   output_image_writer->SetFileName( output_image_filename );
 116   try
 117   {
 118     ts = std::chrono::high_resolution_clock::now( );
 119     output_image_writer->Update( );
 120     te = std::chrono::high_resolution_clock::now( );
 121     tr = te - ts;
 122     std::cout << "Write time: " << tr.count( ) << " s" << std::endl;
 123   }
 124   catch( std::exception& err )
 125   {
 126     std::cerr << "Error caught: " << err.what( ) << std::endl;
 127     return( 1 );
 128
 129   } // yrt
 130
 131   // Save output signal
 132   std::ofstream osignal( output_signal_filename.c_str( ) );
 133   for(
 134     unsigned long i = 0; i < filter->GetNumberOfEvaluatedThresholds( ); ++i
 135     )
 136   {
 137     double x, y;
 138     TFilter::TPeak p;
 139     filter->GetSignalValues( i, x, y, p );
 140     osignal << x << " " << y << std::endl;
 141
 142   } // rof
 143   osignal.close( );
 144   std::cout << "----------------------------------------------" << std::endl;
 145
 146   return( 0 );
 147 }
 148
 149 // eof - $RCSfile$