appli/CTBronchi/Process.cxx

   1 // =========================================================================
   2 // @author Leonardo Florez Valencia
   3 // @email florez-l@javeriana.edu.co
   4 // =========================================================================
   5
   6 #include <chrono>
   7 #include <fstream>
   8 #include <limits>
   9 #include <map>
  10 #include <set>
  11 #include <sstream>
  12 #include <string>
  13 #include <itkImage.h>
  14 #include <itkImageFileReader.h>
  15 #include <itkImageFileWriter.h>
  16 #include <fpa/Filters/Image/Mori.h>
  17 #include <CTBronchi/MoriLabelling.h>
  18
  19 // -------------------------------------------------------------------------
  20 const unsigned int Dim = 3;
  21 typedef short         TInputPixel;
  22 typedef unsigned char TLabelPixel;
  23 typedef itk::Image< TInputPixel, Dim > TInputImage;
  24 typedef itk::Image< TLabelPixel, Dim > TLabelImage;
  25 typedef std::map< std::string, std::string > TMap;
  26
  27 // -------------------------------------------------------------------------
  28 TMap Args;
  29 TInputImage::PointType global_seed;
  30 TLabelPixel inside_value = std::numeric_limits< TLabelPixel >::max( );
  31 TLabelPixel outside_value = TLabelPixel( 0 );
  32
  33 // -------------------------------------------------------------------------
  34 double MeasureTime( itk::ProcessObject* f )
  35 {
  36   std::chrono::time_point< std::chrono::high_resolution_clock > s, e;
  37   std::chrono::duration< double > t;
  38   s = std::chrono::high_resolution_clock::now( );
  39   f->Update( );
  40   e = std::chrono::high_resolution_clock::now( );
  41   t = e - s;
  42   return( t.count( ) );
  43 }
  44
  45 // -------------------------------------------------------------------------
  46 template< class _TImagePtr >
  47 void ReadImage( _TImagePtr& image, const std::string& fname )
  48 {
  49   typedef typename _TImagePtr::ObjectType _TImage;
  50   typedef itk::ImageFileReader< _TImage > _TReader;
  51   typename _TReader::Pointer reader = _TReader::New( );
  52   reader->SetFileName( fname );
  53   double t = MeasureTime( reader );
  54   std::cout << "Read " << fname << " in " << t << " s" << std::endl;
  55   image = reader->GetOutput( );
  56   image->DisconnectPipeline( );
  57 }
  58
  59 // -------------------------------------------------------------------------
  60 template< class _TImagePtr >
  61 void WriteImage( const _TImagePtr& image, const std::string& fname )
  62 {
  63   typedef typename _TImagePtr::ObjectType _TImage;
  64   typedef itk::ImageFileWriter< _TImage > _TWriter;
  65   typename _TWriter::Pointer writer = _TWriter::New( );
  66   writer->SetFileName( fname );
  67   writer->SetInput( image );
  68   double t = MeasureTime( writer );
  69   std::cout << "Wrote " << fname << " in " << t << " s" << std::endl;
  70 }
  71
  72 // -------------------------------------------------------------------------
  73 template< class _TInputPtr, class _TOutputPtr >
  74 void Mori( _TOutputPtr& output, const _TInputPtr& input )
  75 {
  76   typedef typename _TInputPtr::ObjectType  _TInput;
  77   typedef typename _TOutputPtr::ObjectType _TOutput;
  78   typedef fpa::Filters::Image::Mori< _TInput, _TOutput > _TMori;
  79
  80   typename _TMori::Pointer mori = _TMori::New( );
  81   mori->SetInput( input );
  82   mori->SetSeed( global_seed );
  83   mori->SetInsideValue( inside_value );
  84   mori->SetOutsideValue( outside_value );
  85   mori->SetMinimumThreshold(
  86     TInputPixel( std::atof( Args[ "mori_minimum_threshold" ].c_str( ) ) )
  87     );
  88   mori->SetSignalKernelSize(
  89     std::atoi( Args[ "mori_signal_kernel_size" ].c_str( ) )
  90     );
  91   mori->SetSignalThreshold(
  92     std::atof( Args[ "mori_signal_threshold" ].c_str( ) )
  93     );
  94   mori->SetSignalInfluence(
  95     std::atof( Args[ "mori_signal_influence" ].c_str( ) )
  96     );
  97   mori->SetThresholds(
  98     TInputPixel( std::atof( Args[ "mori_lower_threshold" ].c_str( ) ) ),
  99     TInputPixel( std::atof( Args[ "mori_upper_threshold" ].c_str( ) ) ),
 100     TInputPixel( std::atof( Args[ "mori_delta_threshold" ].c_str( ) ) )
 101     );
 102   double t = MeasureTime( mori );
 103   std::cout << "Mori executed in " << t << " s" << std::endl;
 104   output = mori->GetOutput( );
 105
 106   TMap::const_iterator i = Args.find( "out_mori" );
 107   if( i != Args.end( ) )
 108     WriteImage( output, i->second );
 109
 110   i = Args.find( "out_signal" );
 111   if( i != Args.end( ) )
 112   {
 113     std::stringstream signal;
 114     unsigned long nthr = mori->GetNumberOfEvaluatedThresholds( );
 115     signal << "# nThr = " << nthr << std::endl;
 116     signal << "# Opt  = " << mori->GetOptimumThreshold( ) << std::endl;
 117     for( unsigned long j = 0; j < nthr; ++j )
 118     {
 119       typename _TMori::TPeak p;
 120       double x, y;
 121       mori->GetSignalValues( j, x, y, p );
 122       signal << x << " " << y << std::endl;
 123
 124     } // rof
 125
 126     std::ofstream signals_str( i->second.c_str( ) );
 127     signals_str << signal.str( );
 128     signals_str.close( );
 129
 130   } // fi
 131   output->DisconnectPipeline( );
 132 }
 133
 134 // -------------------------------------------------------------------------
 135 template< class _TInputPtr, class _TOutputPtr >
 136 void Label( _TOutputPtr& output, const _TInputPtr& input, const _TOutputPtr& labels )
 137 {
 138   typedef typename _TInputPtr::ObjectType  _TInput;
 139   typedef typename _TOutputPtr::ObjectType _TOutput;
 140   typedef CTBronchi::MoriLabelling< _TInput, _TOutput > _TLabelling;
 141
 142   typename _TLabelling::Pointer labelling = _TLabelling::New( );
 143   labelling->SetInput( input );
 144   labelling->SetInputLabels( labels );
 145   // TODO: labelling->SetOutsideValue( ); // out label
 146   // TODO: labelling->SetInsideValue( ); // inside label
 147   // TODO: labelling->SetUpperThreshold( );
 148   double t = MeasureTime( labelling );
 149   std::cout << "Labelling executed in " << t << " s" << std::endl;
 150   output = labelling->GetOutput( );
 151   TMap::const_iterator i = Args.find( "out_labels" );
 152   if( i != Args.end( ) )
 153     WriteImage( output, i->second );
 154   output->DisconnectPipeline( );
 155 }
 156
 157 // -------------------------------------------------------------------------
 158 bool ParseArgs( int argc, char* argv[] )
 159 {
 160   std::set< std::string > mandatory;
 161   mandatory.insert( "in" );
 162   mandatory.insert( "out_segmentation" );
 163   mandatory.insert( "seed_x" );
 164   mandatory.insert( "seed_y" );
 165   mandatory.insert( "seed_z" );
 166
 167   Args[ "mori_minimum_threshold" ] = "-850";
 168   Args[ "mori_signal_kernel_size" ] = "20";
 169   Args[ "mori_signal_threshold" ] = "100";
 170   Args[ "mori_signal_influence" ] = "0.5";
 171   Args[ "mori_lower_threshold" ] = "-1024";
 172   Args[ "mori_upper_threshold" ] = "0";
 173   Args[ "mori_delta_threshold" ] = "1";
 174   for( int i = 1; i < argc; i += 2 )
 175     Args[ argv[ i ] + 1 ] = argv[ i + 1 ];
 176
 177   bool complete = true;
 178   for( std::string t: mandatory )
 179     complete &= ( Args.find( t ) != Args.end( ) );
 180
 181   if( !complete )
 182   {
 183     std::cerr
 184       << "Usage: " << argv[ 0 ] << std::endl
 185       << "\t-in filename" << std::endl
 186       << "\t-out_segmentation filename" << std::endl
 187       << "\t-seed_x value -seed_y value -seed_z value" << std::endl
 188       << "\t[-out_mori filename]" << std::endl
 189       << "\t[-out_signal filename]" << std::endl
 190       << "\t[-out_labels filename]" << std::endl
 191       << "\t[-mori_minimum_threshold value]" << std::endl
 192       << "\t[-mori_signal_kernel_size value]" << std::endl
 193       << "\t[-mori_signal_threshold value]" << std::endl
 194       << "\t[-mori_signal_influence value]" << std::endl
 195       << "\t[-mori_lower_threshold value]" << std::endl
 196       << "\t[-mori_upper_threshold value]" << std::endl
 197       << "\t[-mori_delta_threshold value]" << std::endl;
 198     return( false );
 199
 200   } // fi
 201   return( true );
 202 }
 203
 204 // -------------------------------------------------------------------------
 205 int main( int argc, char* argv[] )
 206 {
 207   try
 208   {
 209     if( ParseArgs( argc, argv ) )
 210     {
 211       // Parse seed
 212       global_seed[ 0 ] = std::atof( Args[ "seed_x" ].c_str( ) );
 213       global_seed[ 1 ] = std::atof( Args[ "seed_y" ].c_str( ) );
 214       global_seed[ 2 ] = std::atof( Args[ "seed_z" ].c_str( ) );
 215
 216       // Read input image
 217       TInputImage::Pointer input_image;
 218       ReadImage( input_image, Args[ "in" ] );
 219
 220       // Mori segmentation
 221       TLabelImage::Pointer mori;
 222       Mori( mori, input_image );
 223
 224       // Create labels
 225       TLabelImage::Pointer labels;
 226       Label( labels, input_image, mori );
 227
 228       return( 0 );
 229
 230     } // fi
 231     return( 1 );
 232   }
 233   catch( std::exception& err )
 234   {
 235     std::cerr
 236       << "===============================" << std::endl
 237       << "Error caught: " << std::endl
 238       << err.what( )
 239       << "===============================" << std::endl
 240       << std::endl;
 241     return( 1 );
 242
 243   } // yrt
 244   /* TODO
 245      if [ -z "$label_upper_threshold" ]; then label_upper_threshold="-600"; fi
 246      if [ -z "$label_inside" ]; then label_inside="1"; fi
 247      if [ -z "$label_outside" ]; then label_outside="2"; fi
 248      if [ -z "$random_walker_beta" ]; then random_walker_beta="20"; fi
 249   */
 250 }
 251
 252 // eof - $RCSfile$