...

[FrontAlgorithms.git] / appli / CTBronchi / MoriSegmentation.cxx

diff --git a/appli/CTBronchi/MoriSegmentation.cxx b/appli/CTBronchi/MoriSegmentation.cxx
new file mode 100644 (file)
index 0000000..0b26649
--- /dev/null
+++ b/appli/CTBronchi/MoriSegmentation.cxx
@@ -0,0 +1,149 @@
+#include <chrono>
+#include <itkImage.h>
+#include <itkImageFileReader.h>
+#include <itkImageFileWriter.h>
+#include <fpa/Image/Mori.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;
+
+// -------------------------------------------------------------------------
+int main( int argc, char* argv[] )
+{
+  // Get arguments
+  if( argc < 17 )
+  {
+    std::cerr
+      << "Usage: " << argv[ 0 ] << std::endl
+      << "   input_image output_image output_signal" << std::endl
+      << "   init_threshold(-1024) end_threshold(0) delta(1)" << std::endl
+      << "   minimum_threshold(-850)" << std::endl
+      << "   inside_value(255) outside_value(0)" << std::endl
+      << "   signal_kernel_size(20) signal_threshold(500) signal_influence(0.5)"
+      << std::endl
+      << "   [index/point] 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 ] );
+  std::string seed_type = argv[ 13 ];
+  TInputImage::IndexType iseed;
+  TInputImage::PointType pseed;
+  for( unsigned int i = 0; i < Dim; ++i )
+  {
+    if( seed_type == "index" )
+      iseed[ i ] = std::atoi( argv[ 14 + i ] );
+    else
+      pseed[ i ] = std::atof( argv[ 14 + i ] );
+
+  } // rof
+
+  // 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( ) );
+  if( seed_type == "index" )
+    filter->SetSeed( iseed );
+  else
+    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
+
+  // 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 );
+  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;
+  }
+  catch( std::exception& err )
+  {
+    std::cerr << "Error caught: " << err.what( ) << 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 );
+}
+
+// eof - $RCSfile$