-#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 );
}