set(_pfx fpa_CTBronchi_)
set(
_examples
- Process
+ MoriSegmentation
+ MoriLabelling
)
foreach(_e ${_examples})
add_executable(${_pfx}${_e} ${_e}.cxx)
+// =========================================================================
+// @author Leonardo Florez Valencia
+// @email florez-l@javeriana.edu.co
+// =========================================================================
+
#include <chrono>
+#include <map>
#include <itkImage.h>
#include <itkImageFileReader.h>
#include <itkImageFileWriter.h>
-#include <CTBronchi/MoriLabelling.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;
+// -------------------------------------------------------------------------
+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 );
}
#ifndef ITK_MANUAL_INSTANTIATION
# include <CTBronchi/MoriLabelling.hxx>
#endif // ITK_MANUAL_INSTANTIATION
-
#endif // __CTBronchi__MoriLabelling__h__
// eof - $RCSfile$
// @author Leonardo Florez Valencia
// @email florez-l@javeriana.edu.co
// =========================================================================
-
#ifndef __CTBronchi__MoriLabelling__hxx__
#define __CTBronchi__MoriLabelling__hxx__
}
} // ecapseman
-
#endif // __CTBronchi__MoriLabelling__hxx__
// eof - $RCSfile$
+// =========================================================================
+// @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 );
}
#ifndef __fpa__Common__OriginalRandomWalker__hxx__
#define __fpa__Common__OriginalRandomWalker__hxx__
-#include <cmath>
#include <itkImageRegionConstIteratorWithIndex.h>
#include <itkImageRegionIteratorWithIndex.h>
#include <Eigen/Sparse>
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(
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 );
{
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 );
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;
// 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;
_tests
Gaussian
Original
+ CompareRandomWalkers
)
include_directories(${PROJECT_SOURCE_DIR}/lib ${PROJECT_BINARY_DIR}/lib)
--- /dev/null
+// =========================================================================
+// @author Leonardo Florez Valencia
+// @email florez-l@javeriana.edu.co
+// =========================================================================
+
+#include <chrono>
+#include <iostream>
+#include <string>
+#include <sstream>
+#include <itkImage.h>
+#include <itkImageFileReader.h>
+#include <itkImageFileWriter.h>
+
+#include <fpa/Common/OriginalRandomWalker.h>
+#include <fpa/Filters/Image/RandomWalker.h>
+#include <fpa/Functors/Dijkstra/Image/Gaussian.h>
+
+// -------------------------------------------------------------------------
+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$
{
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 );
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;
// 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;
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 );
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;
git://git.creatis.insa-lyon.fr / FrontAlgorithms.git / commitdiff