appli/CTBronchi/MoriSegmentation.cxx

   1 // =========================================================================
   2 // @author Leonardo Florez Valencia
   3 // @email florez-l@javeriana.edu.co
   4 // =========================================================================
   5
   6 #include <chrono>
   7 #include <map>
   8 #include <itkImage.h>
   9 #include <itkImageFileReader.h>
  10 #include <itkImageFileWriter.h>
  11 #include <fpa/Filters/Image/Mori.h>
  12
  13 // -------------------------------------------------------------------------
  14 const unsigned int Dim = 3;
  15 typedef short         TPixel;
  16 typedef unsigned char TLabel;
  17 typedef itk::Image< TPixel, Dim > TImage;
  18 typedef itk::Image< TLabel, Dim > TLabels;
  19
  20 // -------------------------------------------------------------------------
  21 double MeasureTime( itk::ProcessObject* f )
  22 {
  23   std::chrono::time_point< std::chrono::high_resolution_clock > s, e;
  24   std::chrono::duration< double > t;
  25   s = std::chrono::high_resolution_clock::now( );
  26   f->Update( );
  27   e = std::chrono::high_resolution_clock::now( );
  28   t = e - s;
  29   return( t.count( ) );
  30 }
  31
  32 // -------------------------------------------------------------------------
  33 template< class _TImagePtr >
  34 void ReadImage( _TImagePtr& image, const std::string& fname )
  35 {
  36   typedef typename _TImagePtr::ObjectType _TImage;
  37   typedef itk::ImageFileReader< _TImage > _TReader;
  38   typename _TReader::Pointer reader = _TReader::New( );
  39   reader->SetFileName( fname );
  40   double t = MeasureTime( reader );
  41   std::cout << "Read " << fname << " in " << t << " s" << std::endl;
  42   image = reader->GetOutput( );
  43   image->DisconnectPipeline( );
  44 }
  45
  46 // -------------------------------------------------------------------------
  47 template< class _TImage >
  48 void WriteImage( const _TImage*image, const std::string& fname )
  49 {
  50   typedef itk::ImageFileWriter< _TImage > _TWriter;
  51   typename _TWriter::Pointer writer = _TWriter::New( );
  52   writer->SetFileName( fname );
  53   writer->SetInput( image );
  54   double t = MeasureTime( writer );
  55   std::cout << "Wrote " << fname << " in " << t << " s" << std::endl;
  56 }
  57
  58 // -------------------------------------------------------------------------
  59 bool ParseArgs(
  60   std::map< std::string, std::string >& args, int argc, char* argv[]
  61   )
  62 {
  63   std::set< std::string > mandatory;
  64   mandatory.insert( "in" );
  65   mandatory.insert( "out" );
  66   mandatory.insert( "seed" );
  67
  68   args[ "minimum_threshold" ] = "-850";
  69   args[ "signal_kernel_size" ] = "20";
  70   args[ "signal_threshold" ] = "100";
  71   args[ "signal_influence" ] = "0.5";
  72   args[ "lower_threshold" ] = "-1024";
  73   args[ "upper_threshold" ] = "0";
  74   args[ "delta_threshold" ] = "1";
  75
  76   int i = 1;
  77   while( i < argc )
  78   {
  79     std::string cmd = argv[ i ] + 1;
  80     if( cmd == "seed" )
  81     {
  82       std::stringstream seed;
  83       seed
  84         << argv[ i + 1 ] << ";"
  85         << argv[ i + 2 ] << ";"
  86         << argv[ i + 3 ];
  87       args[ cmd ] = seed.str( );
  88       i += 4;
  89     }
  90     else
  91     {
  92       args[ cmd ] = argv[ i + 1 ];
  93       i += 2;
  94
  95     } // fi
  96
  97   } // elihw
  98
  99   bool complete = true;
 100   for( std::string t: mandatory )
 101     complete &= ( args.find( t ) != args.end( ) );
 102
 103   if( !complete )
 104   {
 105     std::cerr
 106       << "Usage: " << argv[ 0 ] << std::endl
 107       << "\t-in filename" << std::endl
 108       << "\t-out filename" << std::endl
 109       << "\t-seed x y z" << std::endl
 110       << "\t[-out_signal filename]" << std::endl
 111       << "\t[-minimum_threshold value]" << std::endl
 112       << "\t[-signal_kernel_size value]" << std::endl
 113       << "\t[-signal_threshold value]" << std::endl
 114       << "\t[-signal_influence value]" << std::endl
 115       << "\t[-lower_threshold value]" << std::endl
 116       << "\t[-upper_threshold value]" << std::endl
 117       << "\t[-delta_threshold value]" << std::endl;
 118     return( false );
 119
 120   } // fi
 121   return( true );
 122 }
 123
 124 // -------------------------------------------------------------------------
 125 int main( int argc, char* argv[] )
 126 {
 127   std::map< std::string, std::string > args;
 128   try
 129   {
 130     if( ParseArgs( args, argc, argv ) )
 131     {
 132       // Parse seed
 133       TImage::PointType seed;
 134       char* str = new char[ args[ "seed" ].size( ) + 1 ];
 135       strcpy( str, args[ "seed" ].c_str( ) );
 136       seed[ 0 ] = std::atof( strtok( str, ";" ) );
 137       seed[ 1 ] = std::atof( strtok( NULL, ";" ) );
 138       seed[ 2 ] = std::atof( strtok( NULL, ";" ) );
 139
 140       // Read input image
 141       TImage::Pointer input_image;
 142       ReadImage( input_image, args[ "in" ] );
 143
 144       // Mori segmentation
 145       typedef fpa::Filters::Image::Mori< TImage, TLabels > TMori;
 146       TMori::Pointer mori = TMori::New( );
 147       mori->SetInput( input_image );
 148       mori->SetSeed( seed );
 149       mori->SetInsideValue( 1 );
 150       mori->SetOutsideValue( 0 );
 151       mori->SetMinimumThreshold(
 152         TPixel( std::atof( args[ "minimum_threshold" ].c_str( ) ) )
 153         );
 154       mori->SetSignalKernelSize(
 155         std::atoi( args[ "signal_kernel_size" ].c_str( ) )
 156         );
 157       mori->SetSignalThreshold(
 158         std::atof( args[ "signal_threshold" ].c_str( ) )
 159         );
 160       mori->SetSignalInfluence(
 161         std::atof( args[ "signal_influence" ].c_str( ) )
 162         );
 163       mori->SetThresholds(
 164         TPixel( std::atof( args[ "lower_threshold" ].c_str( ) ) ),
 165         TPixel( std::atof( args[ "upper_threshold" ].c_str( ) ) ),
 166         TPixel( std::atof( args[ "delta_threshold" ].c_str( ) ) )
 167         );
 168       double t = MeasureTime( mori );
 169       std::cout << "Mori executed in " << t << " s" << std::endl;
 170       WriteImage( mori->GetOutput( ), args[ "out" ] );
 171
 172       std::map< std::string, std::string >::const_iterator i =
 173         args.find( "out_signal" );
 174       if( i != args.end( ) )
 175       {
 176         std::stringstream signal;
 177         unsigned long nthr = mori->GetNumberOfEvaluatedThresholds( );
 178         signal << "# nThr = " << nthr << std::endl;
 179         signal << "# Opt  = " << mori->GetOptimumThreshold( ) << std::endl;
 180         for( unsigned long j = 0; j < nthr; ++j )
 181         {
 182           typename TMori::TPeak p;
 183           double x, y;
 184           mori->GetSignalValues( j, x, y, p );
 185           signal << x << " " << y << std::endl;
 186
 187         } // rof
 188         std::ofstream signals_str( i->second.c_str( ) );
 189         signals_str << signal.str( );
 190         signals_str.close( );
 191
 192       } // fi
 193     }
 194     else
 195       return( 1 );
 196   }
 197   catch( std::exception& err )
 198   {
 199     std::cerr
 200       << "===============================" << std::endl
 201       << "Error caught: " << std::endl
 202       << err.what( )
 203       << "===============================" << std::endl
 204       << std::endl;
 205     return( 1 );
 206
 207   } // yrt
 208   return( 0 );
 209 }
 210
 211 // eof - $RCSfile$