SOURCE ${_e}.cxx
INSTALL
RECURRENT
- LINKS fpa
+ LINKS fpa cpPlugins::tclap
)
endforeach(_e)
endif(fpa_BUILD_CTBronchi)
#define __CTBronchi__Functions__h__
#include <chrono>
-#include <map>
-#include <sstream>
#include <itkImageFileReader.h>
#include <itkImageFileWriter.h>
typename _TReader::Pointer reader = _TReader::New( );
reader->SetFileName( fname );
double t = MeasureTime( reader );
- std::cout << "Read " << fname << " in " << t << " s" << std::endl;
+ std::cout << "Read \"" << fname << "\" in " << t << " s" << std::endl;
image = reader->GetOutput( );
image->DisconnectPipeline( );
}
// -----------------------------------------------------------------------
- template< class _TImagePtr >
- void WriteImage( const _TImagePtr& image, const std::string& fname )
+ template< class _TImage >
+ void WriteImage( const _TImage* image, const std::string& fname )
{
- typedef typename _TImagePtr::ObjectType _TImage;
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[ "mori_minimum_threshold" ] = "-850";
- args[ "mori_signal_kernel_size" ] = "20";
- args[ "mori_signal_threshold" ] = "100";
- args[ "mori_signal_influence" ] = "0.5";
- args[ "mori_lower_threshold" ] = "-1024";
- args[ "mori_upper_threshold" ] = "0";
- args[ "mori_delta_threshold" ] = "1";
- args[ "labelling_upper_threshold" ] = "-600";
- args[ "random_walker_alpha" ] = "1";
- args[ "random_walker_beta" ] = "100";
-
- 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
- << "\t-in filename" << std::endl
- << "\t-out filename" << std::endl
- << "\t-seed x y z" << std::endl
- << "\t[-out_mori filename]" << std::endl
- << "\t[-out_signal filename]" << std::endl
- << "\t[-out_labels filename]" << std::endl
- << "\t[-out_diff filename]" << std::endl
- << "\t[-mori_minimum_threshold value]" << std::endl
- << "\t[-mori_signal_kernel_size value]" << std::endl
- << "\t[-mori_signal_threshold value]" << std::endl
- << "\t[-mori_signal_influence value]" << std::endl
- << "\t[-mori_lower_threshold value]" << std::endl
- << "\t[-mori_upper_threshold value]" << std::endl
- << "\t[-mori_delta_threshold value]" << std::endl
- << "\t[-labelling_upper_threshold value]" << std::endl
- << "\t[-random_walker_alpha value]" << std::endl
- << "\t[-random_walker_beta value]" << std::endl;
- return( false );
-
- } // fi
- return( true );
+ std::cout << "Wrote \"" << fname << "\" in " << t << " s" << std::endl;
}
} // ecapseman
// @email florez-l@javeriana.edu.co
// =========================================================================
-#include <chrono>
-#include <map>
-
+#include <sstream>
+#include <string>
+#include <tclap/CmdLine.h>
#include <itkImage.h>
-#include <itkImageFileReader.h>
-#include <itkImageFileWriter.h>
-#include <itkRegionOfInterestImageFilter.h>
-#include <itkMinimumMaximumImageCalculator.h>
-
-#include <fpa/Filters/BaseMarksInterface.h>
-#include <fpa/Filters/Image/SeedsFromLabelsInterface.h>
-#include <fpa/Filters/Image/DefaultTraits.h>
-#include <fpa/Filters/Image/RegionGrow.h>
-#include <fpa/Functors/RegionGrow/BinaryThreshold.h>
+#include "MoriLabelling.h"
+#include "Functions.h"
// -------------------------------------------------------------------------
const unsigned int Dim = 3;
typedef itk::Image< TLabel, Dim > TLabels;
typedef itk::Image< TScalar, Dim > TScalarImage;
-/**
- */
-class MoriLabellingTraits
- : public fpa::Filters::Image::DefaultTraits< TImage, TLabels, TLabel >
-{
-public:
- typedef fpa::Filters::Image::DefaultTraits< TImage, TLabels, TLabel > Superclass;
- typedef typename Superclass::TInternalTraits TInternalTraits;
- typedef typename Superclass::TMarksImage TMarksImage;
- typedef typename Superclass::TFilterInterface TFilterInterface;
-
- typedef fpa::Filters::BaseMarksInterface< TInternalTraits > TMarksInterface;
- typedef fpa::Filters::Image::SeedsFromLabelsInterface< TInternalTraits > TSeedsInterface;
-};
-
-/**
- */
-class MoriLabelling
- : public fpa::Filters::Image::RegionGrow< TImage, TLabels, TLabel, MoriLabellingTraits >
-{
-public:
- typedef MoriLabellingTraits TTraits;
- fpaTraitsMacro( TTraits );
-
- typedef fpa::Filters::Image::RegionGrow< TImage, TLabels, TMark, TTraits > Superclass;
- typedef MoriLabelling Self;
- typedef itk::SmartPointer< Self > Pointer;
- typedef itk::SmartPointer< const Self > ConstPointer;
- typedef fpa::Functors::RegionGrow::BinaryThreshold< TInputValue > TFunctor;
-
-public:
- itkNewMacro( Self );
- itkTypeMacro( MoriLabelling, fpa::Filters::Image::RegionGrow );
-
- itkGetConstMacro( LastThreshold, TInputValue );
- itkSetMacro( LastThreshold, TInputValue );
-
- itkGetConstMacro( MinVertex, TVertex );
- itkGetConstMacro( MaxVertex, TVertex );
-
- ivqITKInputMacro( InputLabels, TLabels );
- ivqITKInputMacro( InputVesselness, TScalarImage );
-
-public:
- TInputValue GetUpperThreshold( ) const
- {
- return( this->m_Functor->GetUpperThreshold( ) );
- }
- void SetUpperThreshold( TInputValue t )
- {
- this->m_Functor->SetUpperThreshold( t );
- }
-
-protected:
- MoriLabelling( )
- : Superclass( ),
- m_LastThreshold( TInputValue( 0 ) ),
- m_VesselnessThr( TScalar( 0.05 ) )
- {
- ivqITKInputConfigureMacro( InputLabels, TLabels );
- ivqITKInputConfigureMacro( InputVesselness, TScalarImage );
- this->m_Functor = TFunctor::New( );
- this->SetPredicate( this->m_Functor );
- }
-
- virtual ~MoriLabelling( )
- {
- }
-
- virtual const itk::DataObject* _GetReferenceInput( ) const override
- {
- return( this->GetInputLabels( ) );
- }
- virtual void _BeforeGenerateData( ) override
- {
- this->Superclass::_BeforeGenerateData( );
- this->m_FirstVertex = true;
-
- typedef itk::MinimumMaximumImageCalculator< TScalarImage > _TMinMax;
- _TMinMax::Pointer minMax = _TMinMax::New( );
- minMax->SetImage( this->GetInputVesselness( ) );
- minMax->Compute( );
- this->m_MaxVesselness = ( 1.0 - this->m_VesselnessThr ) * minMax->GetMaximum( );
- }
-
- virtual void _PostComputeOutputValue( TNode& n ) override
- {
- this->Superclass::_PostComputeOutputValue( n );
- if( n.Value == this->GetInsideValue( ) )
- {
- const TImage* input = this->GetInput( );
- const TLabels* labels = this->GetInputLabels( );
- const TScalarImage* vesselness = this->GetInputVesselness( );
- double x = input->GetPixel( n.Vertex );
- /* TODO
- double a = std::fabs( x - double( this->m_Functor->GetUpperThreshold( ) ) );
- double b = std::fabs( x - double( this->m_LastThreshold ) );
- */
- if( labels->GetPixel( n.Vertex ) == 0 )
- {
- if( vesselness->GetPixel( n.Vertex ) > this->m_MaxVesselness )
- n.Value = this->GetInsideValue( );
- else
- n.Value = 0;
-
- } // fi
-
- if( !( this->m_FirstVertex ) )
- {
- for( unsigned int d = 0; d < TImage::ImageDimension; ++d )
- {
- if( n.Vertex[ d ] < this->m_MinVertex[ d ] )
- this->m_MinVertex[ d ] = n.Vertex[ d ];
- if( this->m_MaxVertex[ d ] < n.Vertex[ d ] )
- this->m_MaxVertex[ d ] = n.Vertex[ d ];
-
- } // rof
- }
- else
- {
- this->m_MinVertex = n.Vertex;
- this->m_MaxVertex = n.Vertex;
- this->m_FirstVertex = false;
-
- } // fi
-
- } // fi
- }
-
-private:
- // Purposely not implemented.
- MoriLabelling( const Self& other );
- Self& operator=( const Self& other );
-
-protected:
- TFunctor::Pointer m_Functor;
- TInputValue m_LastThreshold;
- bool m_FirstVertex;
- TVertex m_MinVertex;
- TVertex m_MaxVertex;
-
- TScalar m_MaxVesselness;
- TScalar m_VesselnessThr;
-};
-
-// -------------------------------------------------------------------------
-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 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 )
-{
- 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;
-}
-
-// -------------------------------------------------------------------------
-template< class _TImagePtr, class _TRegion >
-void ROI(
- _TImagePtr& output,
- const typename _TImagePtr::ObjectType* input,
- const _TRegion& roi
- )
-{
- typedef typename _TImagePtr::ObjectType _TImage;
- typedef itk::RegionOfInterestImageFilter< _TImage, _TImage > _TFilter;
-
- typename _TFilter::Pointer filter = _TFilter::New( );
- filter->SetInput( input );
- filter->SetRegionOfInterest( roi );
- double t = MeasureTime( filter );
- std::cout << "ROI computed in " << t << " s" << std::endl;
- output = filter->GetOutput( );
- output->DisconnectPipeline( );
-}
-
// -------------------------------------------------------------------------
-bool ParseArgs(
- std::map< std::string, std::string >& args, int argc, char* argv[]
- )
+int main( int argc, char* argv[] )
{
- std::set< std::string > mandatory;
- mandatory.insert( "in" );
- mandatory.insert( "out" );
- mandatory.insert( "labels" );
- mandatory.insert( "vesselness" );
-
- args[ "upper_threshold" ] = "-650";
-
- int i = 1;
- while( i < argc )
+ typedef TCLAP::ValueArg< std::string > _TStringArg;
+ typedef TCLAP::ValueArg< double > _TRealArg;
+ typedef TCLAP::ValueArg< TPixel > _TPixelArg;
+
+ // Parse input line
+ _TStringArg in( "i", "input", "Input image", true, "", "file" );
+ _TStringArg labels( "l", "labels", "Input labels", true, "", "file" );
+ _TStringArg vesselness( "v", "vesselness", "Input vesselness", true, "", "file" );
+ _TStringArg out( "o", "output", "Output image", true, "", "file" );
+ _TRealArg vThr( "a", "vesselness_thr", "Vesselness threshold", false, 0.05, "value (0.05)" );
+ _TPixelArg uThr( "u", "upper_thr", "Upper threshold", false, -650, "value (-650)" );
+ try
{
- 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 )
+ TCLAP::CmdLine cmd( "Labelling", ' ', "1.0.0" );
+ cmd.add( uThr );
+ cmd.add( vThr );
+ cmd.add( out );
+ cmd.add( vesselness );
+ cmd.add( labels );
+ cmd.add( in );
+ cmd.parse( argc, argv );
+ }
+ catch( TCLAP::ArgException& err )
{
std::cerr
- << "Usage: " << argv[ 0 ] << std::endl
- << "\t-in filename" << std::endl
- << "\t-out filename" << std::endl
- << "\t-labels filename" << std::endl
- << "\t-vesselness filename" << std::endl
- << "\t[-roi] filename" << std::endl
- << "\t[-upper_threshold value]" << std::endl;
- return( false );
-
- } // fi
- return( true );
-}
+ << "===============================" << std::endl
+ << "Error caught: " << std::endl
+ << err.error( ) << " " << err.argId( ) << std::endl
+ << "===============================" << std::endl
+ << std::endl;
+ return( 1 );
+
+ } // yrt
-// -------------------------------------------------------------------------
-int main( int argc, char* argv[] )
-{
- std::map< std::string, std::string > args;
try
{
- 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" ] );
-
- // Read vesselness image
- TScalarImage::Pointer input_vesselness;
- ReadImage( input_vesselness, args[ "vesselness" ] );
-
- // Mori labelling
- MoriLabelling::Pointer labelling = MoriLabelling::New( );
- labelling->SetInput( input_image );
- labelling->SetInputLabels( input_labels );
- labelling->SetInputVesselness( input_vesselness );
- labelling->SetOutsideValue( 2 );
- labelling->SetFillValue( 2 );
- labelling->SetInsideValue( 1 );
- labelling->SetUpperThreshold(
- TPixel( std::atof( args[ "upper_threshold" ].c_str( ) ) )
- );
- /* TODO
- labelling->SetLastThreshold( last_thr );
- */
- double t = MeasureTime( labelling );
- std::cout << "Labelling executed in " << t << " s" << std::endl;
-
- std::map< std::string, std::string >::const_iterator i =
- args.find( "roi" );
- if( i != args.end( ) )
- {
- TImage::IndexType minV = labelling->GetMinVertex( );
- TImage::IndexType maxV = labelling->GetMaxVertex( );
- TImage::SizeType roiS;
- for( unsigned d = 0; d < TImage::ImageDimension; ++d )
- roiS[ d ] = maxV[ d ] - minV[ d ] + 1;
- TImage::RegionType roi;
- roi.SetIndex( minV );
- roi.SetSize( roiS );
-
- // ROI input image
- TImage::Pointer input_image_roi;
- ROI( input_image_roi, input_image.GetPointer( ), roi );
- WriteImage( input_image_roi.GetPointer( ), i->second );
-
- // ROI output image
- TLabels::Pointer output_roi;
- ROI( output_roi, labelling->GetOutput( ), roi );
- WriteImage( output_roi.GetPointer( ), args[ "out" ] );
- }
- else
- WriteImage( labelling->GetOutput( ), args[ "out" ] );
- }
- else
- return( 1 );
+ // Read input image
+ TImage::Pointer input_image;
+ CTBronchi::ReadImage( input_image, in.getValue( ) );
+
+ // Read input labels
+ TLabels::Pointer input_labels;
+ CTBronchi::ReadImage( input_labels, labels.getValue( ) );
+
+ // Read input vesselness
+ TScalarImage::Pointer input_vesselness;
+ CTBronchi::ReadImage( input_vesselness, vesselness.getValue( ) );
+
+ // Create labels
+ typedef CTBronchi::MoriLabelling< TImage, TLabels, TScalarImage > _TLabelling;
+ _TLabelling::Pointer labelling = _TLabelling::New( );
+ labelling->SetInput( input_image );
+ labelling->SetInputLabels( input_labels );
+ labelling->SetInputVesselness( input_vesselness );
+ labelling->SetVesselnessThreshold( vThr.getValue( ) );
+ labelling->SetUpperThreshold( uThr.getValue( ) );
+ labelling->SetInsideValue( 1 );
+ labelling->SetOutsideValue( 2 );
+ labelling->SetFillValue( 2 );
+ double t = CTBronchi::MeasureTime( labelling );
+ std::cout << "Labelling executed in " << t << " s" << std::endl;
+
+ // Write result
+ CTBronchi::WriteImage( labelling->GetOutput( ), out.getValue( ) );
}
catch( std::exception& err )
{
std::cerr
<< "===============================" << std::endl
<< "Error caught: " << std::endl
- << err.what( )
+ << err.what( ) << std::endl
<< "===============================" << std::endl
<< std::endl;
return( 1 );
return( 0 );
}
+/* TODO
+ #include <chrono>
+ #include <map>
+
+ #include <itkImage.h>
+ #include <itkImageFileReader.h>
+ #include <itkImageFileWriter.h>
+ #include <itkRegionOfInterestImageFilter.h>
+ #include <itkMinimumMaximumImageCalculator.h>
+
+ #include <fpa/Filters/BaseMarksInterface.h>
+ #include <fpa/Filters/Image/SeedsFromLabelsInterface.h>
+ #include <fpa/Filters/Image/DefaultTraits.h>
+ #include <fpa/Filters/Image/RegionGrow.h>
+ #include <fpa/Functors/RegionGrow/BinaryThreshold.h>
+
+ // -------------------------------------------------------------------------
+ const unsigned int Dim = 3;
+ typedef short TPixel;
+ typedef unsigned char TLabel;
+ typedef itk::NumericTraits< TPixel >::RealType TScalar;
+ typedef itk::Image< TPixel, Dim > TImage;
+ typedef itk::Image< TLabel, Dim > TLabels;
+ typedef itk::Image< TScalar, Dim > TScalarImage;
+
+ class MoriLabellingTraits
+ : public fpa::Filters::Image::DefaultTraits< TImage, TLabels, TLabel >
+ {
+ public:
+ typedef fpa::Filters::Image::DefaultTraits< TImage, TLabels, TLabel > Superclass;
+ typedef typename Superclass::TInternalTraits TInternalTraits;
+ typedef typename Superclass::TMarksImage TMarksImage;
+ typedef typename Superclass::TFilterInterface TFilterInterface;
+
+ typedef fpa::Filters::BaseMarksInterface< TInternalTraits > TMarksInterface;
+ typedef fpa::Filters::Image::SeedsFromLabelsInterface< TInternalTraits > TSeedsInterface;
+ };
+
+ class MoriLabelling
+ : public fpa::Filters::Image::RegionGrow< TImage, TLabels, TLabel, MoriLabellingTraits >
+ {
+ public:
+ typedef MoriLabellingTraits TTraits;
+ fpaTraitsMacro( TTraits );
+
+ typedef fpa::Filters::Image::RegionGrow< TImage, TLabels, TMark, TTraits > Superclass;
+ typedef MoriLabelling Self;
+ typedef itk::SmartPointer< Self > Pointer;
+ typedef itk::SmartPointer< const Self > ConstPointer;
+ typedef fpa::Functors::RegionGrow::BinaryThreshold< TInputValue > TFunctor;
+
+ public:
+ itkNewMacro( Self );
+ itkTypeMacro( MoriLabelling, fpa::Filters::Image::RegionGrow );
+
+ itkGetConstMacro( LastThreshold, TInputValue );
+ itkSetMacro( LastThreshold, TInputValue );
+
+ itkGetConstMacro( MinVertex, TVertex );
+ itkGetConstMacro( MaxVertex, TVertex );
+
+ ivqITKInputMacro( InputLabels, TLabels );
+ ivqITKInputMacro( InputVesselness, TScalarImage );
+
+ public:
+ TInputValue GetUpperThreshold( ) const
+ {
+ return( this->m_Functor->GetUpperThreshold( ) );
+ }
+ void SetUpperThreshold( TInputValue t )
+ {
+ this->m_Functor->SetUpperThreshold( t );
+ }
+
+ protected:
+ MoriLabelling( )
+ : Superclass( ),
+ m_LastThreshold( TInputValue( 0 ) ),
+ m_VesselnessThr( TScalar( 0.05 ) )
+ {
+ ivqITKInputConfigureMacro( InputLabels, TLabels );
+ ivqITKInputConfigureMacro( InputVesselness, TScalarImage );
+ this->m_Functor = TFunctor::New( );
+ this->SetPredicate( this->m_Functor );
+ }
+
+ virtual ~MoriLabelling( )
+ {
+ }
+
+ virtual const itk::DataObject* _GetReferenceInput( ) const override
+ {
+ return( this->GetInputLabels( ) );
+ }
+ virtual void _BeforeGenerateData( ) override
+ {
+ this->Superclass::_BeforeGenerateData( );
+ this->m_FirstVertex = true;
+
+ typedef itk::MinimumMaximumImageCalculator< TScalarImage > _TMinMax;
+ _TMinMax::Pointer minMax = _TMinMax::New( );
+ minMax->SetImage( this->GetInputVesselness( ) );
+ minMax->Compute( );
+ this->m_MaxVesselness = ( 1.0 - this->m_VesselnessThr ) * minMax->GetMaximum( );
+ }
+
+ virtual void _PostComputeOutputValue( TNode& n ) override
+ {
+ this->Superclass::_PostComputeOutputValue( n );
+ if( n.Value == this->GetInsideValue( ) )
+ {
+ const TImage* input = this->GetInput( );
+ const TLabels* labels = this->GetInputLabels( );
+ const TScalarImage* vesselness = this->GetInputVesselness( );
+ double x = input->GetPixel( n.Vertex );
+ if( labels->GetPixel( n.Vertex ) == 0 )
+ {
+ if( vesselness->GetPixel( n.Vertex ) > this->m_MaxVesselness )
+ n.Value = this->GetInsideValue( );
+ else
+ n.Value = 0;
+
+ } // fi
+
+ if( !( this->m_FirstVertex ) )
+ {
+ for( unsigned int d = 0; d < TImage::ImageDimension; ++d )
+ {
+ if( n.Vertex[ d ] < this->m_MinVertex[ d ] )
+ this->m_MinVertex[ d ] = n.Vertex[ d ];
+ if( this->m_MaxVertex[ d ] < n.Vertex[ d ] )
+ this->m_MaxVertex[ d ] = n.Vertex[ d ];
+
+ } // rof
+ }
+ else
+ {
+ this->m_MinVertex = n.Vertex;
+ this->m_MaxVertex = n.Vertex;
+ this->m_FirstVertex = false;
+
+ } // fi
+
+ } // fi
+ }
+
+ private:
+ // Purposely not implemented.
+ MoriLabelling( const Self& other );
+ Self& operator=( const Self& other );
+
+ protected:
+ TFunctor::Pointer m_Functor;
+ TInputValue m_LastThreshold;
+ bool m_FirstVertex;
+ TVertex m_MinVertex;
+ TVertex m_MaxVertex;
+
+ TScalar m_MaxVesselness;
+ TScalar m_VesselnessThr;
+ };
+
+ // -------------------------------------------------------------------------
+ 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 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 )
+ {
+ 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;
+ }
+
+ // -------------------------------------------------------------------------
+ template< class _TImagePtr, class _TRegion >
+ void ROI(
+ _TImagePtr& output,
+ const typename _TImagePtr::ObjectType* input,
+ const _TRegion& roi
+ )
+ {
+ typedef typename _TImagePtr::ObjectType _TImage;
+ typedef itk::RegionOfInterestImageFilter< _TImage, _TImage > _TFilter;
+
+ typename _TFilter::Pointer filter = _TFilter::New( );
+ filter->SetInput( input );
+ filter->SetRegionOfInterest( roi );
+ double t = MeasureTime( filter );
+ std::cout << "ROI computed in " << t << " s" << std::endl;
+ output = filter->GetOutput( );
+ output->DisconnectPipeline( );
+ }
+
+ // -------------------------------------------------------------------------
+ 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( "labels" );
+ mandatory.insert( "vesselness" );
+
+ args[ "upper_threshold" ] = "-650";
+
+ 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
+ << "\t-in filename" << std::endl
+ << "\t-out filename" << std::endl
+ << "\t-labels filename" << std::endl
+ << "\t-vesselness filename" << std::endl
+ << "\t[-roi] filename" << std::endl
+ << "\t[-upper_threshold value]" << std::endl;
+ return( false );
+
+ } // fi
+ return( true );
+ }
+
+ // -------------------------------------------------------------------------
+ int main( int argc, char* argv[] )
+ {
+ std::map< std::string, std::string > args;
+ try
+ {
+ 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" ] );
+
+ // Read vesselness image
+ TScalarImage::Pointer input_vesselness;
+ ReadImage( input_vesselness, args[ "vesselness" ] );
+
+ // Mori labelling
+ MoriLabelling::Pointer labelling = MoriLabelling::New( );
+ labelling->SetInput( input_image );
+ labelling->SetInputLabels( input_labels );
+ labelling->SetInputVesselness( input_vesselness );
+ labelling->SetOutsideValue( 2 );
+ labelling->SetFillValue( 2 );
+ labelling->SetInsideValue( 1 );
+ labelling->SetUpperThreshold(
+ TPixel( std::atof( args[ "upper_threshold" ].c_str( ) ) )
+ );
+ double t = MeasureTime( labelling );
+ std::cout << "Labelling executed in " << t << " s" << std::endl;
+
+ std::map< std::string, std::string >::const_iterator i =
+ args.find( "roi" );
+ if( i != args.end( ) )
+ {
+ TImage::IndexType minV = labelling->GetMinVertex( );
+ TImage::IndexType maxV = labelling->GetMaxVertex( );
+ TImage::SizeType roiS;
+ for( unsigned d = 0; d < TImage::ImageDimension; ++d )
+ roiS[ d ] = maxV[ d ] - minV[ d ] + 1;
+ TImage::RegionType roi;
+ roi.SetIndex( minV );
+ roi.SetSize( roiS );
+
+ // ROI input image
+ TImage::Pointer input_image_roi;
+ ROI( input_image_roi, input_image.GetPointer( ), roi );
+ WriteImage( input_image_roi.GetPointer( ), i->second );
+
+ // ROI output image
+ TLabels::Pointer output_roi;
+ ROI( output_roi, labelling->GetOutput( ), roi );
+ WriteImage( output_roi.GetPointer( ), args[ "out" ] );
+ }
+ else
+ WriteImage( labelling->GetOutput( ), args[ "out" ] );
+ }
+ else
+ return( 1 );
+ }
+ catch( std::exception& err )
+ {
+ std::cerr
+ << "===============================" << std::endl
+ << "Error caught: " << std::endl
+ << err.what( )
+ << "===============================" << std::endl
+ << std::endl;
+ return( 1 );
+
+ } // yrt
+ return( 0 );
+ }
+*/
+
// eof - $RCSfile$
--- /dev/null
+// =========================================================================
+// @author Leonardo Florez Valencia
+// @email florez-l@javeriana.edu.co
+// =========================================================================
+#ifndef __CTBronchi__MoriLabelling__h__
+#define __CTBronchi__MoriLabelling__h__
+
+#include <fpa/Filters/BaseMarksInterface.h>
+#include <fpa/Filters/Image/SeedsFromLabelsInterface.h>
+#include <fpa/Filters/Image/DefaultTraits.h>
+#include <fpa/Filters/Image/RegionGrow.h>
+#include <fpa/Functors/RegionGrow/BinaryThreshold.h>
+
+namespace CTBronchi
+{
+ /**
+ */
+ template< class _TImage, class _TLabels >
+ class MoriLabellingTraits
+ : public fpa::Filters::Image::DefaultTraits< _TImage, _TLabels, typename _TLabels::PixelType >
+ {
+ public:
+ typedef fpa::Filters::Image::DefaultTraits< _TImage, _TLabels, typename _TLabels::PixelType > Superclass;
+ typedef typename Superclass::TInternalTraits TInternalTraits;
+ typedef typename Superclass::TMarksImage TMarksImage;
+ typedef typename Superclass::TFilterInterface TFilterInterface;
+
+ typedef fpa::Filters::BaseMarksInterface< TInternalTraits > TMarksInterface;
+ typedef fpa::Filters::Image::SeedsFromLabelsInterface< TInternalTraits > TSeedsInterface;
+ };
+
+ /**
+ */
+ template< class _TImage, class _TLabels, class _TScalarImage >
+ class MoriLabelling
+ : public fpa::Filters::Image::RegionGrow< _TImage, _TLabels, typename _TLabels::PixelType, CTBronchi::MoriLabellingTraits< _TImage, _TLabels > >
+ {
+ public:
+ typedef CTBronchi::MoriLabellingTraits< _TImage, _TLabels > TTraits;
+ fpaTraitsMacro( typename TTraits );
+
+ typedef fpa::Filters::Image::RegionGrow< _TImage, _TLabels, TMark, TTraits > Superclass;
+ typedef MoriLabelling Self;
+ typedef itk::SmartPointer< Self > Pointer;
+ typedef itk::SmartPointer< const Self > ConstPointer;
+
+ typedef fpa::Functors::RegionGrow::BinaryThreshold< TInputValue > TFunctor;
+
+ public:
+ itkNewMacro( Self );
+ itkTypeMacro( MoriLabelling, fpa::Filters::Image::RegionGrow );
+
+ itkGetConstMacro( VesselnessThreshold, double );
+ itkSetMacro( VesselnessThreshold, double );
+
+ itkGetConstMacro( UpperThreshold, TInputValue );
+ itkSetMacro( UpperThreshold, TInputValue );
+
+ ivqITKInputMacro( InputLabels, _TLabels );
+ ivqITKInputMacro( InputVesselness, _TScalarImage );
+
+ protected:
+ MoriLabelling( )
+ : Superclass( ),
+ m_VesselnessThreshold( 0.05 ),
+ m_UpperThreshold( -650 )
+ {
+ ivqITKInputConfigureMacro( InputLabels, _TLabels );
+ ivqITKInputConfigureMacro( InputVesselness, _TScalarImage );
+ this->m_Functor = TFunctor::New( );
+ this->SetPredicate( this->m_Functor );
+ }
+
+ virtual ~MoriLabelling( )
+ {
+ }
+
+ private:
+ // Purposely not implemented.
+ MoriLabelling( const Self& other );
+ Self& operator=( const Self& other );
+
+ protected:
+ typename TFunctor::Pointer m_Functor;
+ double m_VesselnessThreshold;
+ TInputValue m_UpperThreshold;
+ };
+
+} // ecapseman
+
+#endif // __CTBronchi__Functions__h__
+
+// eof - $RCSfile$
// @email florez-l@javeriana.edu.co
// =========================================================================
-#include <chrono>
-#include <map>
+#include <sstream>
+#include <string>
+#include <tclap/CmdLine.h>
#include <itkImage.h>
-#include <itkImageFileReader.h>
-#include <itkImageFileWriter.h>
#include <fpa/Filters/Image/Mori.h>
+#include "Functions.h"
// -------------------------------------------------------------------------
const unsigned int Dim = 3;
-typedef short TPixel;
-typedef unsigned char TLabel;
+typedef short TPixel;
+typedef unsigned char TLabel;
typedef itk::Image< TPixel, Dim > TImage;
typedef itk::Image< TLabel, Dim > TLabels;
// -------------------------------------------------------------------------
-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 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 )
-{
- 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[]
- )
+int main( 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 )
+ 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
{
- 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 )
+ 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( TCLAP::ArgException& err )
{
std::cerr
- << "Usage: " << argv[ 0 ] << std::endl
- << "\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 );
+ << "===============================" << std::endl
+ << "Error caught: " << std::endl
+ << err.error( ) << " " << err.argId( )
+ << "===============================" << std::endl
+ << std::endl;
+ return( 1 );
- } // fi
- return( true );
-}
+ } // yrt
+
+ // Get seed
+ std::istringstream sSeed( aSeed.getValue( ) );
+ TImage::PointType seed;
+ sSeed >> seed[ 0 ] >> seed[ 1 ] >> seed[ 2 ];
-// -------------------------------------------------------------------------
-int main( int argc, char* argv[] )
-{
- std::map< std::string, std::string > args;
try
{
- if( ParseArgs( args, argc, argv ) )
+ // 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( ) != "" )
{
- // 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 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 )
{
- 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;
+ typename _TMori::TPeak p;
+ double x, y;
+ mori->GetSignalValues( j, x, y, p );
+ signalStr << x << " " << y << std::endl;
- } // rof
- std::ofstream signals_str( i->second.c_str( ) );
- signals_str << signal.str( );
- signals_str.close( );
+ } // rof
+ std::ofstream oStrSignals( signal.getValue( ).c_str( ) );
+ oStrSignals << signalStr.str( );
+ oStrSignals.close( );
- } // fi
- }
- else
- return( 1 );
+ } // fi
}
catch( std::exception& err )
{
// @email florez-l@javeriana.edu.co
// =========================================================================
-#include <chrono>
-#include <map>
+#include <string>
+#include <tclap/CmdLine.h>
#include <itkImage.h>
-#include <itkImageFileReader.h>
-#include <itkImageFileWriter.h>
#include <itkMinimumMaximumImageCalculator.h>
#include <itkInvertIntensityImageFilter.h>
#include <itkHessianRecursiveGaussianImageFilter.h>
#include <itkHessian3DToVesselnessMeasureImageFilter.h>
+#include "Functions.h"
// -------------------------------------------------------------------------
const unsigned int Dim = 3;
-typedef short TPixel;
+typedef short TPixel;
typedef itk::NumericTraits< TPixel >::RealType TScalar;
-typedef itk::Image< TPixel, Dim > TImage;
-typedef itk::Image< TScalar, Dim > TScalarImage;
+typedef itk::Image< TPixel, Dim > TImage;
+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 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 )
-{
- 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[]
- )
+int main( int argc, char* argv[] )
{
- std::set< std::string > mandatory;
- mandatory.insert( "in" );
- mandatory.insert( "out" );
-
- args[ "sigma" ] = "0.5";
- args[ "alpha1" ] = "0.5";
- args[ "alpha2" ] = "2";
-
- int i = 1;
- while( i < argc )
+ typedef TCLAP::ValueArg< std::string > _TStringArg;
+ typedef TCLAP::ValueArg< TScalar > _TRealArg;
+
+ // Parse input line
+ _TStringArg in( "i", "input", "Input image", true, "", "file" );
+ _TStringArg out( "o", "output", "Output image", true, "", "file" );
+ _TRealArg s( "s", "sigma", "Sigma", false, 0.5, "value (0.5)" );
+ _TRealArg a1( "a", "alpha1", "Alpha1", false, 0.5, "value (0.5)" );
+ _TRealArg a2( "b", "alpha2", "Alpha2", false, 2, "value (2)" );
+ try
{
- 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 )
+ TCLAP::CmdLine cmd( "Vesselness computation", ' ', "1.0.0" );
+ cmd.add( a2 );
+ cmd.add( a1 );
+ cmd.add( s );
+ cmd.add( out );
+ cmd.add( in );
+ cmd.parse( argc, argv );
+ }
+ catch( TCLAP::ArgException& err )
{
std::cerr
- << "Usage: " << argv[ 0 ] << std::endl
- << "\t-in filename" << std::endl
- << "\t-out filename" << std::endl
- << "\t[-sigma val]" << std::endl
- << "\t[-alpha1 val]" << std::endl
- << "\t[-alpha2 val]" << std::endl;
- return( false );
+ << "===============================" << std::endl
+ << "Error caught: " << std::endl
+ << err.error( ) << " " << err.argId( )
+ << "===============================" << std::endl
+ << std::endl;
+ return( 1 );
- } // fi
- return( true );
-}
+ } // yrt
-// -------------------------------------------------------------------------
-int main( int argc, char* argv[] )
-{
- std::map< std::string, std::string > args;
try
{
- if( ParseArgs( args, argc, argv ) )
- {
- // Read input image
- TImage::Pointer input_image;
- ReadImage( input_image, args[ "in" ] );
-
- // Vesselness
- typedef itk::MinimumMaximumImageCalculator< TImage > TMinMax;
- typedef itk::InvertIntensityImageFilter< TImage > TInverter;
- typedef itk::HessianRecursiveGaussianImageFilter< TImage > THessian;
- typedef itk::Hessian3DToVesselnessMeasureImageFilter< TScalar > TVesselness;
-
- TMinMax::Pointer minMax = TMinMax::New( );
- minMax->SetImage( input_image );
- minMax->Compute( );
-
- TInverter::Pointer inverter = TInverter::New( );
- inverter->SetInput( input_image );
- inverter->SetMaximum( minMax->GetMaximum( ) );
- double t = MeasureTime( inverter );
- std::cout << "Inversion executed in " << t << " s" << std::endl;
-
- THessian::Pointer hessian = THessian::New( );
- hessian->SetInput( inverter->GetOutput( ) );
- hessian->SetSigma(
- std::atof( args[ "sigma" ].c_str( ) )
- );
- t = MeasureTime( hessian );
- std::cout << "Hessian executed in " << t << " s" << std::endl;
-
- TVesselness::Pointer vesselness = TVesselness::New( );
- vesselness->SetInput( hessian->GetOutput( ) );
- vesselness->SetAlpha1(
- std::atof( args[ "alpha1" ].c_str( ) )
- );
- vesselness->SetAlpha2(
- std::atof( args[ "alpha2" ].c_str( ) )
- );
- t = MeasureTime( vesselness );
- std::cout << "Vesselness executed in " << t << " s" << std::endl;
-
- WriteImage( vesselness->GetOutput( ), args[ "out" ] );
- }
- else
- return( 1 );
+ // Read image
+ TImage::Pointer input_image;
+ CTBronchi::ReadImage( input_image, in.getValue( ) );
+
+ // Min-max
+ typedef itk::MinimumMaximumImageCalculator< TImage > _TMinMax;
+ _TMinMax::Pointer minMax = _TMinMax::New( );
+ minMax->SetImage( input_image );
+ minMax->Compute( );
+
+ // Invert intensities
+ typedef itk::InvertIntensityImageFilter< TImage > _TInverter;
+ _TInverter::Pointer inverter = _TInverter::New( );
+ inverter->SetInput( input_image );
+ inverter->SetMaximum( minMax->GetMaximum( ) );
+ double t = CTBronchi::MeasureTime( inverter );
+ std::cout << "Inversion executed in " << t << " s" << std::endl;
+
+ // Compute hessian image
+ typedef itk::HessianRecursiveGaussianImageFilter< TImage > _THessian;
+ _THessian::Pointer hessian = _THessian::New( );
+ hessian->SetInput( inverter->GetOutput( ) );
+ hessian->SetSigma( s.getValue( ) );
+ t = CTBronchi::MeasureTime( hessian );
+ std::cout << "Hessian executed in " << t << " s" << std::endl;
+
+ // Vesselness
+ typedef
+ itk::Hessian3DToVesselnessMeasureImageFilter< TScalar >
+ _TVesselness;
+ _TVesselness::Pointer vesselness = _TVesselness::New( );
+ vesselness->SetInput( hessian->GetOutput( ) );
+ vesselness->SetAlpha1( a1.getValue( ) );
+ vesselness->SetAlpha2( a2.getValue( ) );
+ t = CTBronchi::MeasureTime( vesselness );
+ std::cout << "Vesselness executed in " << t << " s" << std::endl;
+
+ // Write result
+ CTBronchi::WriteImage( vesselness->GetOutput( ), out.getValue( ) );
}
catch( std::exception& err )
{
git://git.creatis.insa-lyon.fr / FrontAlgorithms.git / commitdiff