]> Creatis software - FrontAlgorithms.git/blobdiff - appli/CTBronchi/MoriSegmentation.cxx
...
[FrontAlgorithms.git] / appli / CTBronchi / MoriSegmentation.cxx
index 39b9e0dcb3bcec9978554fa19e2c6bdec5a8f35d..90a1bf3333b608aa93392ab9210a82bde94310db 100644 (file)
-#include <chrono>
+// =========================================================================
+// @author Leonardo Florez Valencia
+// @email florez-l@javeriana.edu.co
+// =========================================================================
+
+#include <sstream>
+#include <string>
+#include <tclap/CmdLine.h>
 #include <itkImage.h>
-#include <itkImageFileReader.h>
-#include <itkImageFileWriter.h>
-#include <fpa/Image/Mori.h>
+#include <fpa/Filters/Image/Mori.h>
+#include "Functions.h"
 
 // -------------------------------------------------------------------------
 const unsigned int Dim = 3;
-typedef short          TPixel;
-typedef unsigned short TLabel;
-
-typedef itk::Image< TPixel, Dim > TInputImage;
-typedef itk::Image< TLabel, Dim > TLabelImage;
-typedef fpa::Image::Mori< TInputImage, TLabelImage > TFilter;
+typedef short                     TPixel;
+typedef unsigned char             TLabel;
+typedef itk::Image< TPixel, Dim > TImage;
+typedef itk::Image< TLabel, Dim > TLabels;
 
 // -------------------------------------------------------------------------
 int main( int argc, char* argv[] )
 {
-  // Get arguments
-  if( argc < 16 )
-  {
-    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 );
-
-  } // 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;
+  typedef TCLAP::ValueArg< std::string >   _TStringArg;
+  typedef TCLAP::ValueArg< double >        _TRealArg;
+  typedef TCLAP::ValueArg< unsigned long > _TUIntArg;
+  typedef TCLAP::ValueArg< TPixel >        _TPixelArg;
+
+  // Parse input line
+  _TStringArg in( "i", "input", "Input image", true, "", "file" );
+  _TStringArg out( "o", "output", "Output image", true, "", "file" );
+  _TStringArg signal( "s", "signal", "Output signal", false, "", "file" );
+  _TStringArg aSeed( "p", "seed", "Input seed", true, "", "\"x y z\"" );
+  _TPixelArg mThr( "t", "min_thr", "Minimum threshold", false, -850, "value (-850)" );
+  _TUIntArg sKernel( "k", "kernel", "Kernel size", false, 20, "value (20)" );
+  _TRealArg sThr( "r", "signal_thr", "Signal threshold", false, 100, "value (100)" );
+  _TRealArg sInfluence( "f", "signal_influence", "Signal influence", false, 0.5, "value (0.5)" );
+  _TPixelArg lower( "l", "lower", "Lower threshold", false, -1024, "value (-1024)" );
+  _TPixelArg upper( "u", "upper", "Upper threshold", false, 0, "value (0)" );
+  _TPixelArg delta( "d", "delta", "Delta threshold", false, 1, "value (1)" );
   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;
+    TCLAP::CmdLine cmd( "Mori segmentation", ' ', "1.0.0" );
+    cmd.add( delta );
+    cmd.add( upper );
+    cmd.add( lower );
+    cmd.add( sInfluence );
+    cmd.add( sThr );
+    cmd.add( sKernel );
+    cmd.add( mThr );
+    cmd.add( aSeed );
+    cmd.add( signal );
+    cmd.add( out );
+    cmd.add( in );
+    cmd.parse( argc, argv );
   }
-  catch( std::exception& err )
+  catch( TCLAP::ArgException& err )
   {
-    std::cerr << "Error caught: " << err.what( ) << std::endl;
+    std::cerr
+      << "===============================" << std::endl
+      << "Error caught: " << std::endl
+      << err.error( ) << " " << err.argId( )
+      << "===============================" << std::endl
+      << std::endl;
     return( 1 );
 
   } // yrt
 
-  // 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 );
+  // Get seed
+  std::istringstream sSeed( aSeed.getValue( ) );
+  TImage::PointType seed;
+  sSeed >> seed[ 0 ] >> seed[ 1 ] >> seed[ 2 ];
+
   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;
+    // Read image
+    TImage::Pointer input_image;
+    CTBronchi::ReadImage( input_image, in.getValue( ) );
+
+    // 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( mThr.getValue( ) );
+    mori->SetSignalKernelSize( sKernel.getValue( ) );
+    mori->SetSignalThreshold( sThr.getValue( ) );
+    mori->SetSignalInfluence( sInfluence.getValue( ) );
+    mori->SetThresholds(
+      lower.getValue( ), upper.getValue( ), delta.getValue( )
+      );
+    double t = CTBronchi::MeasureTime( mori );
+    std::cout << "Mori executed in " << t << " s" << std::endl;
+
+    // Write result
+    CTBronchi::WriteImage( mori->GetOutput( ), out.getValue( ) );
+
+    // Write signal
+    if( signal.getValue( ) != "" )
+    {
+      std::stringstream signalStr;
+      unsigned long nthr = mori->GetNumberOfEvaluatedThresholds( );
+      signalStr << "# nThr = " << nthr << std::endl;
+      signalStr << "# 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 );
+        signalStr << x << " " << y << std::endl;
+
+      } // rof
+      std::ofstream oStrSignals( signal.getValue( ).c_str( ) );
+      oStrSignals << signalStr.str( );
+      oStrSignals.close( );
+
+    } // fi
   }
   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 );
 }