From: Leonardo Flórez-Valencia Date: Thu, 21 Sep 2017 16:53:32 +0000 (+0200) Subject: ... X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=commitdiff_plain;h=d2f5533aef455b7d862dbfbe72c42c5497eaf98b;p=FrontAlgorithms.git ... --- diff --git a/appli/CTBronchi/CMakeLists.txt b/appli/CTBronchi/CMakeLists.txt index 1851e55..fccfb97 100644 --- a/appli/CTBronchi/CMakeLists.txt +++ b/appli/CTBronchi/CMakeLists.txt @@ -14,7 +14,8 @@ if(fpa_BUILD_CTBronchi) set(_pfx fpa_CTBronchi_) set( _examples - Process + MoriSegmentation + MoriLabelling ) foreach(_e ${_examples}) add_executable(${_pfx}${_e} ${_e}.cxx) diff --git a/appli/CTBronchi/MoriLabelling.cxx b/appli/CTBronchi/MoriLabelling.cxx index 0915cd5..0ad694d 100644 --- a/appli/CTBronchi/MoriLabelling.cxx +++ b/appli/CTBronchi/MoriLabelling.cxx @@ -1,117 +1,181 @@ +// ========================================================================= +// @author Leonardo Florez Valencia +// @email florez-l@javeriana.edu.co +// ========================================================================= + #include +#include #include #include #include -#include +#include // ------------------------------------------------------------------------- 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; +// ------------------------------------------------------------------------- +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( ) ); +} -typedef CTBronchi::MoriLabelling< TInputImage, TLabelImage > TFilter; +// ------------------------------------------------------------------------- +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( ); +} // ------------------------------------------------------------------------- -int main( int argc, char* argv[] ) +template< class _TImage > +void WriteImage( const _TImage*image, const std::string& fname ) +{ + 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[] + ) { - // Get arguments - if( argc < 8 ) + std::set< std::string > mandatory; + mandatory.insert( "in" ); + mandatory.insert( "out" ); + mandatory.insert( "labels" ); + + args[ "upper_threshold" ] = "-600"; + + int i = 1; + while( i < argc ) + { + std::string cmd = argv[ i ] + 1; + args[ cmd ] = argv[ i + 1 ]; + i += 2; + + } // 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 label_image output_image" << std::endl - << " upper_threshold(-400)" << std::endl - << " inside_value(255)" << std::endl - << " inside_label(1) outside_label(2)" - << std::endl; - return( 1 ); + << "\t-in filename" << std::endl + << "\t-out filename" << std::endl + << "\t-labels filename" << std::endl + << "\t[-upper_threshold value]" << std::endl; + return( false ); } // fi - std::string input_image_filename = argv[ 1 ]; - std::string label_image_filename = argv[ 2 ]; - std::string output_image_filename = argv[ 3 ]; - TLabel upper_threshold = std::atoi( argv[ 4 ] ); - TLabel inside_value = std::atoi( argv[ 5 ] ); - TLabel inside_label = std::atoi( argv[ 6 ] ); - TLabel outside_label = std::atoi( argv[ 7 ] ); - - // Read images - itk::ImageFileReader< TInputImage >::Pointer input_image_reader = - itk::ImageFileReader< TInputImage >::New( ); - input_image_reader->SetFileName( input_image_filename ); - - itk::ImageFileReader< TLabelImage >::Pointer label_image_reader = - itk::ImageFileReader< TLabelImage >::New( ); - label_image_reader->SetFileName( label_image_filename ); - - // Prepare filter - TFilter::Pointer filter = TFilter::New( ); - filter->SetInputLabelImage( label_image_reader->GetOutput( ) ); - filter->SetInputRawImage( input_image_reader->GetOutput( ) ); - filter->SetUpperThreshold( upper_threshold ); - filter->SetInsideValue( inside_value ); - filter->SetInsideLabel( inside_label ); - filter->SetOutsideLabel( outside_label ); - - // 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 << "Raw read time: " << tr.count( ) << " s" << std::endl; - - ts = std::chrono::high_resolution_clock::now( ); - label_image_reader->Update( ); - te = std::chrono::high_resolution_clock::now( ); - tr = te - ts; - std::cout << "Label 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 - << "Labelling time: " << tr.count( ) << " s" << 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->GetOutput( ) ); - 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 ) ) + { + // Read input image + TImage::Pointer input_image; + ReadImage( input_image, args[ "in" ] ); + + // Read labels image + TLabels::Pointer input_labels; + ReadImage( input_labels, args[ "labels" ] ); + + // Mori segmentation + /* TODO + 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 - std::cout << "----------------------------------------------" << std::endl; - return( 0 ); } diff --git a/appli/CTBronchi/MoriLabelling.h b/appli/CTBronchi/MoriLabelling.h index b599310..d07968b 100644 --- a/appli/CTBronchi/MoriLabelling.h +++ b/appli/CTBronchi/MoriLabelling.h @@ -86,7 +86,6 @@ namespace CTBronchi #ifndef ITK_MANUAL_INSTANTIATION # include #endif // ITK_MANUAL_INSTANTIATION - #endif // __CTBronchi__MoriLabelling__h__ // eof - $RCSfile$ diff --git a/appli/CTBronchi/MoriLabelling.hxx b/appli/CTBronchi/MoriLabelling.hxx index b6f982a..f786948 100644 --- a/appli/CTBronchi/MoriLabelling.hxx +++ b/appli/CTBronchi/MoriLabelling.hxx @@ -2,7 +2,6 @@ // @author Leonardo Florez Valencia // @email florez-l@javeriana.edu.co // ========================================================================= - #ifndef __CTBronchi__MoriLabelling__hxx__ #define __CTBronchi__MoriLabelling__hxx__ @@ -107,7 +106,6 @@ namespace CTBronchi } } // ecapseman - #endif // __CTBronchi__MoriLabelling__hxx__ // eof - $RCSfile$ diff --git a/appli/CTBronchi/MoriSegmentation.cxx b/appli/CTBronchi/MoriSegmentation.cxx index 39b9e0d..dcbae23 100644 --- a/appli/CTBronchi/MoriSegmentation.cxx +++ b/appli/CTBronchi/MoriSegmentation.cxx @@ -1,137 +1,210 @@ +// ========================================================================= +// @author Leonardo Florez Valencia +// @email florez-l@javeriana.edu.co +// ========================================================================= + #include +#include #include #include #include -#include +#include // ------------------------------------------------------------------------- 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 ); } diff --git a/lib/fpa/Common/OriginalRandomWalker.hxx b/lib/fpa/Common/OriginalRandomWalker.hxx index 4b60c1b..46dd935 100644 --- a/lib/fpa/Common/OriginalRandomWalker.hxx +++ b/lib/fpa/Common/OriginalRandomWalker.hxx @@ -5,7 +5,6 @@ #ifndef __fpa__Common__OriginalRandomWalker__hxx__ #define __fpa__Common__OriginalRandomWalker__hxx__ -#include #include #include #include @@ -247,7 +246,7 @@ _Laplacian( _TTriplets& A, _TTriplets& R, const _TTriplets& B ) thrStr.B = reinterpret_cast< const void* >( &B ); // Configure threader - const TImage* in = this->GetInputLabels( ); + const TImage* in = this->GetInput( ); const itk::ImageRegionSplitterBase* split = this->GetImageRegionSplitter( ); const unsigned int nThreads = split->GetNumberOfSplits( diff --git a/lib/fpa/Functors/Dijkstra/Image/Gaussian.h b/lib/fpa/Functors/Dijkstra/Image/Gaussian.h index 2effb51..b819f01 100644 --- a/lib/fpa/Functors/Dijkstra/Image/Gaussian.h +++ b/lib/fpa/Functors/Dijkstra/Image/Gaussian.h @@ -63,7 +63,7 @@ namespace fpa d -= TValue( image->GetPixel( p ) ); d /= this->m_Beta; d *= d; - if( this->m_TreatAsWeight ) d = std::exp( d ); + if( this->m_TreatAsWeight ) d = std::exp( d ) - TValue( 1 ); else d = std::exp( -d ); if( d < this->m_Epsilon ) return( this->m_Epsilon ); diff --git a/tests/image/Dijkstra/Gaussian.cxx b/tests/image/Dijkstra/Gaussian.cxx index 3139e83..ef6d827 100644 --- a/tests/image/Dijkstra/Gaussian.cxx +++ b/tests/image/Dijkstra/Gaussian.cxx @@ -24,7 +24,8 @@ int main( int argc, char* argv[] ) { std::cerr << "Usage: " << argv[ 0 ] - << " input_image output_image output_marks output_mst alpha beta [seeds]" + << " input_image output_image output_marks output_mst" + << " beta epsilon [seeds]" << std::endl; return( 1 ); @@ -33,8 +34,8 @@ int main( int argc, char* argv[] ) std::string output_image_filename = argv[ 2 ]; std::string output_marks_filename = argv[ 3 ]; std::string output_mst_filename = argv[ 4 ]; - double alpha = std::atof( argv[ 5 ] ); - double beta = std::atof( argv[ 6 ] ); + double beta = std::atof( argv[ 5 ] ); + double eps = std::atof( argv[ 6 ] ); // Read image typedef itk::ImageFileReader< TInputImage > TInputImageReader; @@ -44,8 +45,8 @@ int main( int argc, char* argv[] ) // Prepare weight functor typedef fpa::Functors::Dijkstra::Image::Gaussian< TInputImage, TOutputPixel > TWeight; TWeight::Pointer weight = TWeight::New( ); - weight->SetAlpha( alpha ); weight->SetBeta( beta ); + weight->SetEpsilon( eps ); // Prepare filter typedef fpa::Filters::Image::Dijkstra< TInputImage, TOutputImage > TFilter; diff --git a/tests/image/RandomWalker/CMakeLists.txt b/tests/image/RandomWalker/CMakeLists.txt index 859b3a6..e0aef14 100644 --- a/tests/image/RandomWalker/CMakeLists.txt +++ b/tests/image/RandomWalker/CMakeLists.txt @@ -8,6 +8,7 @@ set( _tests Gaussian Original + CompareRandomWalkers ) include_directories(${PROJECT_SOURCE_DIR}/lib ${PROJECT_BINARY_DIR}/lib) diff --git a/tests/image/RandomWalker/CompareRandomWalkers.cxx b/tests/image/RandomWalker/CompareRandomWalkers.cxx new file mode 100644 index 0000000..67afb48 --- /dev/null +++ b/tests/image/RandomWalker/CompareRandomWalkers.cxx @@ -0,0 +1,192 @@ +// ========================================================================= +// @author Leonardo Florez Valencia +// @email florez-l@javeriana.edu.co +// ========================================================================= + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +// ------------------------------------------------------------------------- +const unsigned int Dim = 3; +typedef short TPixel; +typedef unsigned char TLabel; +typedef double TScalar; + +typedef itk::Image< TPixel, Dim > TImage; +typedef itk::Image< TLabel, Dim > TLabels; +typedef itk::Image< TScalar, Dim > TScalarImage; + +// ------------------------------------------------------------------------- +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( ) ); +} + +// ------------------------------------------------------------------------- +template< class _TImagePtr > +void ReadImage( _TImagePtr& image, const std::string& fname ) +{ + typedef itk::ImageFileReader< typename _TImagePtr::ObjectType > _TReader; + typename _TReader::Pointer r = _TReader::New( ); + r->SetFileName( fname ); + double t = MeasureTime( r ); + image = r->GetOutput( ); + image->DisconnectPipeline( ); + std::cout << "Image \"" << fname << "\" read in " << t << "s" << std::endl; +} + +// ------------------------------------------------------------------------- +template< class _TImagePtr > +void WriteImage( const _TImagePtr& image, const std::string& fname ) +{ + typedef itk::ImageFileWriter< typename _TImagePtr::ObjectType > _TWriter; + typename _TWriter::Pointer w = _TWriter::New( ); + w->SetInput( image ); + w->SetFileName( fname ); + double t = MeasureTime( w ); + std::cout << "Image \"" << fname << "\" written in " << t << "s" << std::endl; +} + +// ------------------------------------------------------------------------- +template< class _TImagePtr, class _TLabelsPtr, class _TScalarImagePtr > +void Original( + _TLabelsPtr& output_labels, + _TScalarImagePtr& output_probabilities, + const _TImagePtr& input, + const _TLabelsPtr& labels, + const typename _TScalarImagePtr::ObjectType::PixelType& beta, + const typename _TScalarImagePtr::ObjectType::PixelType& epsilon + ) +{ + typedef typename _TImagePtr::ObjectType _TImage; + typedef typename _TLabelsPtr::ObjectType _TLabels; + typedef typename _TScalarImagePtr::ObjectType _TScalarImage; + typedef typename _TScalarImage::PixelType _TScalar; + typedef + fpa::Functors::Dijkstra::Image::Gaussian< _TImage, _TScalar > + _TFunction; + typedef + fpa::Common::OriginalRandomWalker< _TImage, _TLabels, _TScalar > + _TFilter; + + // Random walker function + typename _TFunction::Pointer edge = _TFunction::New( ); + edge->SetBeta( beta ); + edge->SetEpsilon( epsilon ); + edge->TreatAsWeightOff( ); + + // Random walker + typename _TFilter::Pointer filter = _TFilter::New( ); + filter->SetInput( input ); + filter->SetInputLabels( labels ); + filter->SetEdgeFunction( edge ); + double t = MeasureTime( filter ); + output_labels = filter->GetOutput( ); + output_probabilities = filter->GetOutputProbabilities( ); + output_labels->DisconnectPipeline( ); + output_probabilities->DisconnectPipeline( ); + std::cout << "Original RW executed in " << t << "s" << std::endl; +} + +// ------------------------------------------------------------------------- +template< class _TImagePtr, class _TLabelsPtr, class _TScalarImagePtr > +void FrontPropagation( + _TLabelsPtr& output_labels, + _TScalarImagePtr& output_probabilities, + const _TImagePtr& input, + const _TLabelsPtr& labels, + const typename _TScalarImagePtr::ObjectType::PixelType& beta, + const typename _TScalarImagePtr::ObjectType::PixelType& epsilon + ) +{ + typedef typename _TImagePtr::ObjectType _TImage; + typedef typename _TLabelsPtr::ObjectType _TLabels; + typedef typename _TScalarImagePtr::ObjectType _TScalarImage; + typedef typename _TScalarImage::PixelType _TScalar; + typedef + fpa::Functors::Dijkstra::Image::Gaussian< _TImage, _TScalar > + _TFunction; + typedef + fpa::Filters::Image::RandomWalker< _TImage, _TLabels, _TScalar > + _TFilter; + + // Random walker function + typename _TFunction::Pointer edge = _TFunction::New( ); + edge->SetBeta( beta ); + edge->SetEpsilon( epsilon ); + edge->TreatAsWeightOn( ); + + // Random walker + typename _TFilter::Pointer filter = _TFilter::New( ); + filter->SetInputImage( input ); + filter->SetInputLabels( labels ); + filter->SetWeightFunction( edge ); + double t = MeasureTime( filter ); + output_labels = filter->GetOutputLabels( ); + output_probabilities = filter->GetOutputCosts( ); + output_labels->DisconnectPipeline( ); + output_probabilities->DisconnectPipeline( ); + std::cout << "Front propagation RW executed in " << t << "s" << std::endl; +} + +// ------------------------------------------------------------------------- +int main( int argc, char* argv[] ) +{ + // Get input arguments + if( argc < 5 ) + { + std::cerr + << "Usage: " << argv[ 0 ] + << " input_image input_labels beta epsilon" + << std::endl; + return( 1 ); + + } // fi + std::string input_image_filename = argv[ 1 ]; + std::string input_labels_filename = argv[ 2 ]; + std::istringstream values( std::string( argv[ 3 ] ) + " " + argv[ 4 ] ); + TScalar beta, epsilon; + values >> beta >> epsilon; + + // Read inputs + TImage::Pointer input_image; + TLabels::Pointer input_labels; + ReadImage( input_image, input_image_filename ); + ReadImage( input_labels, input_labels_filename ); + + // Execute algorithms + TLabels::Pointer original_labels, fp_labels; + TScalarImage::Pointer original_probabilities, fp_probabilities; + Original( + original_labels, original_probabilities, + input_image, input_labels, + beta, epsilon + ); + FrontPropagation( + fp_labels, fp_probabilities, + input_image, input_labels, + beta, epsilon + ); + + WriteImage( original_labels, "orig.png" ); + WriteImage( fp_labels, "fp.png" ); + + return( 0 ); +} + +// eof - $RCSfile$ diff --git a/tests/image/RandomWalker/Gaussian.cxx b/tests/image/RandomWalker/Gaussian.cxx index 82d44ae..5c2c798 100644 --- a/tests/image/RandomWalker/Gaussian.cxx +++ b/tests/image/RandomWalker/Gaussian.cxx @@ -25,7 +25,7 @@ int main( int argc, char* argv[] ) { std::cerr << "Usage: " << argv[ 0 ] - << " input_image input_label output_image output_costs alpha beta" + << " input_image input_label output_image output_costs beta epsilon" << std::endl; return( 1 ); @@ -34,8 +34,8 @@ int main( int argc, char* argv[] ) std::string input_label_filename = argv[ 2 ]; std::string output_image_filename = argv[ 3 ]; std::string output_costs_filename = argv[ 4 ]; - double alpha = std::atof( argv[ 5 ] ); - double beta = std::atof( argv[ 6 ] ); + double beta = std::atof( argv[ 5 ] ); + double eps = std::atof( argv[ 6 ] ); // Read image typedef itk::ImageFileReader< TImage > TImageReader; @@ -50,8 +50,8 @@ int main( int argc, char* argv[] ) // Prepare weight functor typedef fpa::Functors::Dijkstra::Image::Gaussian< TImage, TCost > TWeight; TWeight::Pointer weight = TWeight::New( ); - weight->SetAlpha( alpha ); weight->SetBeta( beta ); + weight->SetEpsilon( eps ); // Prepare filter typedef fpa::Filters::Image::RandomWalker< TImage, TLabelImage, TCost > TFilter; diff --git a/tests/image/RandomWalker/Original.cxx b/tests/image/RandomWalker/Original.cxx index 0ea2bd0..4bc9fbe 100644 --- a/tests/image/RandomWalker/Original.cxx +++ b/tests/image/RandomWalker/Original.cxx @@ -22,12 +22,12 @@ typedef itk::Image< TLabel, Dim > TLabels; int main( int argc, char* argv[] ) { // Get input arguments - if( argc < 8 ) + if( argc < 7 ) { std::cerr << "Usage: " << argv[ 0 ] << " input_image input_labels output_labels output_probabilities" - << " beta eps normalize" + << " beta eps" << std::endl; return( 1 ); @@ -35,14 +35,12 @@ int main( int argc, char* argv[] ) std::string input_image_filename, input_labels_filename; std::string output_labels_filename, output_probabilities_filename; double beta, eps; - bool normalize; input_image_filename = argv[ 1 ]; input_labels_filename = argv[ 2 ]; output_labels_filename = argv[ 3 ]; output_probabilities_filename = argv[ 4 ]; beta = std::atof( argv[ 5 ] ); eps = std::atof( argv[ 6 ] ); - normalize = ( argv[ 7 ][ 0 ] == '1' ); // Read image typedef itk::ImageFileReader< TImage > TImageReader;