]> Creatis software - FrontAlgorithms.git/blobdiff - appli/CTBronchi/MoriSegmentation.cxx
...
[FrontAlgorithms.git] / appli / CTBronchi / MoriSegmentation.cxx
index 39b9e0dcb3bcec9978554fa19e2c6bdec5a8f35d..dcbae23addd5d556e37814759cf05673691874e9 100644 (file)
+// =========================================================================
+// @author Leonardo Florez Valencia
+// @email florez-l@javeriana.edu.co
+// =========================================================================
+
 #include <chrono>
+#include <map>
 #include <itkImage.h>
 #include <itkImageFileReader.h>
 #include <itkImageFileWriter.h>
-#include <fpa/Image/Mori.h>
+#include <fpa/Filters/Image/Mori.h>
 
 // -------------------------------------------------------------------------
 const unsigned int Dim = 3;
-typedef short          TPixel;
-typedef unsigned short TLabel;
+typedef short         TPixel;
+typedef unsigned char TLabel;
+typedef itk::Image< TPixel, Dim > TImage;
+typedef itk::Image< TLabel, Dim > TLabels;
 
-typedef itk::Image< TPixel, Dim > TInputImage;
-typedef itk::Image< TLabel, Dim > TLabelImage;
-typedef fpa::Image::Mori< TInputImage, TLabelImage > TFilter;
+// -------------------------------------------------------------------------
+double MeasureTime( itk::ProcessObject* f )
+{
+  std::chrono::time_point< std::chrono::high_resolution_clock > s, e;
+  std::chrono::duration< double > t;
+  s = std::chrono::high_resolution_clock::now( );
+  f->Update( );
+  e = std::chrono::high_resolution_clock::now( );
+  t = e - s;
+  return( t.count( ) );
+}
 
 // -------------------------------------------------------------------------
-int main( int argc, char* argv[] )
+template< class _TImagePtr >
+void ReadImage( _TImagePtr& image, const std::string& fname )
+{
+  typedef typename _TImagePtr::ObjectType _TImage;
+  typedef itk::ImageFileReader< _TImage > _TReader;
+  typename _TReader::Pointer reader = _TReader::New( );
+  reader->SetFileName( fname );
+  double t = MeasureTime( reader );
+  std::cout << "Read " << fname << " in " << t << " s" << std::endl;
+  image = reader->GetOutput( );
+  image->DisconnectPipeline( );
+}
+
+// -------------------------------------------------------------------------
+template< class _TImage >
+void WriteImage( const _TImage*image, const std::string& fname )
 {
-  // Get arguments
-  if( argc < 16 )
+  typedef itk::ImageFileWriter< _TImage > _TWriter;
+  typename _TWriter::Pointer writer = _TWriter::New( );
+  writer->SetFileName( fname );
+  writer->SetInput( image );
+  double t = MeasureTime( writer );
+  std::cout << "Wrote " << fname << " in " << t << " s" << std::endl;
+}
+
+// -------------------------------------------------------------------------
+bool ParseArgs(
+  std::map< std::string, std::string >& args, int argc, char* argv[]
+  )
+{
+  std::set< std::string > mandatory;
+  mandatory.insert( "in" );
+  mandatory.insert( "out" );
+  mandatory.insert( "seed" );
+
+  args[ "minimum_threshold" ] = "-850";
+  args[ "signal_kernel_size" ] = "20";
+  args[ "signal_threshold" ] = "100";
+  args[ "signal_influence" ] = "0.5";
+  args[ "lower_threshold" ] = "-1024";
+  args[ "upper_threshold" ] = "0";
+  args[ "delta_threshold" ] = "1";
+
+  int i = 1;
+  while( i < argc )
+  {
+    std::string cmd = argv[ i ] + 1;
+    if( cmd == "seed" )
+    {
+      std::stringstream seed;
+      seed
+        << argv[ i + 1 ] << ";"
+        << argv[ i + 2 ] << ";"
+        << argv[ i + 3 ];
+      args[ cmd ] = seed.str( );
+      i += 4;
+    }
+    else
+    {
+      args[ cmd ] = argv[ i + 1 ];
+      i += 2;
+
+    } // fi
+
+  } // elihw
+
+  bool complete = true;
+  for( std::string t: mandatory )
+    complete &= ( args.find( t ) != args.end( ) );
+
+  if( !complete )
   {
     std::cerr
       << "Usage: " << argv[ 0 ] << std::endl
-      << "   input_image output_image output_signal" << std::endl
-      << "   init_threshold end_threshold delta" << std::endl
-      << "   minimum_threshold" << std::endl
-      << "   inside_value outside_value" << std::endl
-      << "   signal_kernel_size signal_threshold signal_influence"
-      << std::endl
-      << "   seed_x seed_y seed_z"
-      << std::endl;
-    return( 1 );
+      << "\t-in filename" << std::endl
+      << "\t-out filename" << std::endl
+      << "\t-seed x y z" << std::endl
+      << "\t[-out_signal filename]" << std::endl
+      << "\t[-minimum_threshold value]" << std::endl
+      << "\t[-signal_kernel_size value]" << std::endl
+      << "\t[-signal_threshold value]" << std::endl
+      << "\t[-signal_influence value]" << std::endl
+      << "\t[-lower_threshold value]" << std::endl
+      << "\t[-upper_threshold value]" << std::endl
+      << "\t[-delta_threshold value]" << std::endl;
+    return( false );
 
   } // fi
-  std::string input_image_filename = argv[ 1 ];
-  std::string output_image_filename = argv[ 2 ];
-  std::string output_signal_filename = argv[ 3 ];
-  TPixel init_threshold = std::atoi( argv[ 4 ] );
-  TPixel end_threshold = std::atoi( argv[ 5 ] );
-  TPixel delta = std::atoi( argv[ 6 ] );
-  TPixel minimum_threshold = std::atoi( argv[ 7 ] );
-  TLabel inside_value = std::atoi( argv[ 8 ] );
-  TLabel outside_value = std::atoi( argv[ 9 ] );
-  unsigned long signal_kernel_size = std::atoi( argv[ 10 ] );
-  double signal_threshold = std::atof( argv[ 11 ] );
-  double signal_influence = std::atof( argv[ 12 ] );
-  TInputImage::PointType pseed;
-  for( int i = 0; i < Dim; ++i )
-    pseed[ i ] = std::atof( argv[ 13 + i ] );
-
-  // Read image
-  itk::ImageFileReader< TInputImage >::Pointer input_image_reader =
-    itk::ImageFileReader< TInputImage >::New( );
-  input_image_reader->SetFileName( input_image_filename );
-
-  // Prepare filter
-  TFilter::Pointer filter = TFilter::New( );
-  filter->SetInput( input_image_reader->GetOutput( ) );
-  filter->SetSeed( pseed );
-  filter->SetThresholds( init_threshold, end_threshold, delta );
-  filter->SetMinimumThreshold( minimum_threshold );
-  filter->SetInsideValue( inside_value );
-  filter->SetOutsideValue( outside_value );
-  filter->SetSignalKernelSize( signal_kernel_size );
-  filter->SetSignalThreshold( signal_threshold );
-  filter->SetSignalInfluence( signal_influence );
-
-  // Show some information
-  std::cout << "----------------------------------------------" << std::endl;
-  std::cout << "Image: " << input_image_filename << std::endl;
-
-  // Execute pipeline
-  std::chrono::time_point< std::chrono::high_resolution_clock > ts, te;
-  std::chrono::duration< double > tr;
-  try
-  {
-    ts = std::chrono::high_resolution_clock::now( );
-    input_image_reader->Update( );
-    te = std::chrono::high_resolution_clock::now( );
-    tr = te - ts;
-    std::cout << "Read time: " << tr.count( ) << " s" << std::endl;
-
-    ts = std::chrono::high_resolution_clock::now( );
-    filter->Update( );
-    te = std::chrono::high_resolution_clock::now( );
-    tr = te - ts;
-    std::cout
-      << "Mori time: " << tr.count( ) << " s" << std::endl
-      << "Optimum threshold: " << filter->GetOptimumThreshold( ) << std::endl
-      << "Number of thresholds: "
-      << filter->GetNumberOfEvaluatedThresholds( ) << std::endl;
-  }
-  catch( std::exception& err )
-  {
-    std::cerr << "Error caught: " << err.what( ) << std::endl;
-    return( 1 );
-
-  } // yrt
+  return( true );
+}
 
-  // Save output image
-  itk::ImageFileWriter< TLabelImage >::Pointer output_image_writer =
-    itk::ImageFileWriter< TLabelImage >::New( );
-  output_image_writer->SetInput( filter->GetThresholdedOutput( ) );
-  output_image_writer->SetFileName( output_image_filename );
+// -------------------------------------------------------------------------
+int main( int argc, char* argv[] )
+{
+  std::map< std::string, std::string > args;
   try
   {
-    ts = std::chrono::high_resolution_clock::now( );
-    output_image_writer->Update( );
-    te = std::chrono::high_resolution_clock::now( );
-    tr = te - ts;
-    std::cout << "Write time: " << tr.count( ) << " s" << std::endl;
+    if( ParseArgs( args, argc, argv ) )
+    {
+      // Parse seed
+      TImage::PointType seed;
+      char* str = new char[ args[ "seed" ].size( ) + 1 ];
+      strcpy( str, args[ "seed" ].c_str( ) );
+      seed[ 0 ] = std::atof( strtok( str, ";" ) );
+      seed[ 1 ] = std::atof( strtok( NULL, ";" ) );
+      seed[ 2 ] = std::atof( strtok( NULL, ";" ) );
+
+      // Read input image
+      TImage::Pointer input_image;
+      ReadImage( input_image, args[ "in" ] );
+
+      // Mori segmentation
+      typedef fpa::Filters::Image::Mori< TImage, TLabels > TMori;
+      TMori::Pointer mori = TMori::New( );
+      mori->SetInput( input_image );
+      mori->SetSeed( seed );
+      mori->SetInsideValue( 1 );
+      mori->SetOutsideValue( 0 );
+      mori->SetMinimumThreshold(
+        TPixel( std::atof( args[ "minimum_threshold" ].c_str( ) ) )
+        );
+      mori->SetSignalKernelSize(
+        std::atoi( args[ "signal_kernel_size" ].c_str( ) )
+        );
+      mori->SetSignalThreshold(
+        std::atof( args[ "signal_threshold" ].c_str( ) )
+        );
+      mori->SetSignalInfluence(
+        std::atof( args[ "signal_influence" ].c_str( ) )
+        );
+      mori->SetThresholds(
+        TPixel( std::atof( args[ "lower_threshold" ].c_str( ) ) ),
+        TPixel( std::atof( args[ "upper_threshold" ].c_str( ) ) ),
+        TPixel( std::atof( args[ "delta_threshold" ].c_str( ) ) )
+        );
+      double t = MeasureTime( mori );
+      std::cout << "Mori executed in " << t << " s" << std::endl;
+      WriteImage( mori->GetOutput( ), args[ "out" ] );
+
+      std::map< std::string, std::string >::const_iterator i =
+        args.find( "out_signal" );
+      if( i != args.end( ) )
+      {
+        std::stringstream signal;
+        unsigned long nthr = mori->GetNumberOfEvaluatedThresholds( );
+        signal << "# nThr = " << nthr << std::endl;
+        signal << "# Opt  = " << mori->GetOptimumThreshold( ) << std::endl;
+        for( unsigned long j = 0; j < nthr; ++j )
+        {
+          typename TMori::TPeak p;
+          double x, y;
+          mori->GetSignalValues( j, x, y, p );
+          signal << x << " " << y << std::endl;
+
+        } // rof
+        std::ofstream signals_str( i->second.c_str( ) );
+        signals_str << signal.str( );
+        signals_str.close( );
+
+      } // fi
+    }
+    else
+      return( 1 );
   }
   catch( std::exception& err )
   {
-    std::cerr << "Error caught: " << err.what( ) << std::endl;
+    std::cerr
+      << "===============================" << std::endl
+      << "Error caught: " << std::endl
+      << err.what( )
+      << "===============================" << std::endl
+      << std::endl;
     return( 1 );
 
   } // yrt
-
-  // Save output signal
-  std::ofstream osignal( output_signal_filename.c_str( ) );
-  for(
-    unsigned long i = 0; i < filter->GetNumberOfEvaluatedThresholds( ); ++i
-    )
-  {
-    double x, y;
-    TFilter::TPeak p;
-    filter->GetSignalValues( i, x, y, p );
-    osignal << x << " " << y << std::endl;
-
-  } // rof
-  osignal.close( );
-  std::cout << "----------------------------------------------" << std::endl;
-
   return( 0 );
 }