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