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 < 16 )
  21   {
  22     std::cerr
  23       << "Usage: " << argv[ 0 ] << std::endl
  24       << "   input_image output_image output_signal" << std::endl
  25       << "   init_threshold end_threshold delta" << std::endl
  26       << "   minimum_threshold" << std::endl
  27       << "   inside_value outside_value" << std::endl
  28       << "   signal_kernel_size signal_threshold signal_influence"
  29       << std::endl
  30       << "   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   TInputImage::PointType pseed;
  48   for( int i = 0; i < Dim; ++i )
  49     pseed[ i ] = std::atof( argv[ 13 + i ] );
  50
  51   // Read image
  52   itk::ImageFileReader< TInputImage >::Pointer input_image_reader =
  53     itk::ImageFileReader< TInputImage >::New( );
  54   input_image_reader->SetFileName( input_image_filename );
  55
  56   // Prepare filter
  57   TFilter::Pointer filter = TFilter::New( );
  58   filter->SetInput( input_image_reader->GetOutput( ) );
  59   filter->SetSeed( pseed );
  60   filter->SetThresholds( init_threshold, end_threshold, delta );
  61   filter->SetMinimumThreshold( minimum_threshold );
  62   filter->SetInsideValue( inside_value );
  63   filter->SetOutsideValue( outside_value );
  64   filter->SetSignalKernelSize( signal_kernel_size );
  65   filter->SetSignalThreshold( signal_threshold );
  66   filter->SetSignalInfluence( signal_influence );
  67
  68   // Show some information
  69   std::cout << "----------------------------------------------" << std::endl;
  70   std::cout << "Image: " << input_image_filename << std::endl;
  71
  72   // Execute pipeline
  73   std::chrono::time_point< std::chrono::high_resolution_clock > ts, te;
  74   std::chrono::duration< double > tr;
  75   try
  76   {
  77     ts = std::chrono::high_resolution_clock::now( );
  78     input_image_reader->Update( );
  79     te = std::chrono::high_resolution_clock::now( );
  80     tr = te - ts;
  81     std::cout << "Read time: " << tr.count( ) << " s" << std::endl;
  82
  83     ts = std::chrono::high_resolution_clock::now( );
  84     filter->Update( );
  85     te = std::chrono::high_resolution_clock::now( );
  86     tr = te - ts;
  87     std::cout
  88       << "Mori time: " << tr.count( ) << " s" << std::endl
  89       << "Optimum threshold: " << filter->GetOptimumThreshold( ) << std::endl
  90       << "Number of thresholds: "
  91       << filter->GetNumberOfEvaluatedThresholds( ) << std::endl;
  92   }
  93   catch( std::exception& err )
  94   {
  95     std::cerr << "Error caught: " << err.what( ) << std::endl;
  96     return( 1 );
  97
  98   } // yrt
  99
 100   // Save output image
 101   itk::ImageFileWriter< TLabelImage >::Pointer output_image_writer =
 102     itk::ImageFileWriter< TLabelImage >::New( );
 103   output_image_writer->SetInput( filter->GetThresholdedOutput( ) );
 104   output_image_writer->SetFileName( output_image_filename );
 105   try
 106   {
 107     ts = std::chrono::high_resolution_clock::now( );
 108     output_image_writer->Update( );
 109     te = std::chrono::high_resolution_clock::now( );
 110     tr = te - ts;
 111     std::cout << "Write time: " << tr.count( ) << " s" << std::endl;
 112   }
 113   catch( std::exception& err )
 114   {
 115     std::cerr << "Error caught: " << err.what( ) << std::endl;
 116     return( 1 );
 117
 118   } // yrt
 119
 120   // Save output signal
 121   std::ofstream osignal( output_signal_filename.c_str( ) );
 122   for(
 123     unsigned long i = 0; i < filter->GetNumberOfEvaluatedThresholds( ); ++i
 124     )
 125   {
 126     double x, y;
 127     TFilter::TPeak p;
 128     filter->GetSignalValues( i, x, y, p );
 129     osignal << x << " " << y << std::endl;
 130
 131   } // rof
 132   osignal.close( );
 133   std::cout << "----------------------------------------------" << std::endl;
 134
 135   return( 0 );
 136 }
 137
 138 // eof - $RCSfile$