#include <streambuf>
#include <tclap/CmdLine.h>
+#include <itkAndImageFilter.h>
+#include <itkBinaryThresholdImageFilter.h>
#include <itkMinimumMaximumImageCalculator.h>
#include <itkInvertIntensityImageFilter.h>
#include <itkHessianRecursiveGaussianImageFilter.h>
#include <itkHessian3DToVesselnessMeasureImageFilter.h>
+#include <fpa/Common/SliceBySliceRandomWalker.h>
+#include <fpa/Filters/Image/ExtractSkeleton.h>
#include <fpa/Filters/Image/Mori.h>
+#include <fpa/Filters/Image/RandomWalker.h>
+#include <fpa/Functors/Dijkstra/Image/Gaussian.h>
#include "MoriLabelling.h"
#include "Process.h"
m_InsideLabel( 1 ),
m_OutsideLabel( 2 )
{
- this->m_StrArg[ "input" ] = TArgument( "", true );
- this->m_StrArg[ "seed_file" ] = TArgument( "", false );
- this->m_StrArg[ "seed" ] = TArgument( "", false );
- this->m_StrArg[ "seed_type" ] = TArgument( "point", false );
-
- this->m_DblArg[ "beta" ] = TArgument( "2.5", false );
- this->m_DblArg[ "epsilon" ] = TArgument( "1e-5", false );
- this->m_DblArg[ "vesselness_sigma" ] = TArgument( "0.5", false );
- this->m_DblArg[ "vesselness_alpha1" ] = TArgument( "0.5", false );
- this->m_DblArg[ "vesselness_alpha2" ] = TArgument( "2", false );
- this->m_DblArg[ "mori_min_thr" ] = TArgument( "-850", false );
- this->m_DblArg[ "mori_kernel" ] = TArgument( "20", false );
- this->m_DblArg[ "mori_signal_thr" ] = TArgument( "100", false );
- this->m_DblArg[ "mori_signal_influence" ] = TArgument( "0.5", false );
- this->m_DblArg[ "mori_lower" ] = TArgument( "-1024", false );
- this->m_DblArg[ "mori_upper" ] = TArgument( "0", false );
- this->m_DblArg[ "mori_delta" ] = TArgument( "1", false );
- this->m_DblArg[ "slicerw_thr" ] = TArgument( "5", false );
- this->m_DblArg[ "fastrw_thr" ] = TArgument( "65", false );
- this->m_DblArg[ "labels_upper_thr" ] = TArgument( "-400", false );
+ this->m_StrArgs[ "input" ] = TStrArg( "", "I", true );
+ this->m_StrArgs[ "vesselness" ] = TStrArg( "", "V", false );
+ this->m_StrArgs[ "mori" ] = TStrArg( "", "M", false );
+ this->m_StrArgs[ "labels" ] = TStrArg( "", "L", false );
+ this->m_StrArgs[ "fastrw" ] = TStrArg( "", "F", false );
+ this->m_StrArgs[ "slicerw" ] = TStrArg( "", "R", false );
+ this->m_StrArgs[ "andrw" ] = TStrArg( "", "A", false );
+ this->m_StrArgs[ "fastrw_skeleton" ] = TStrArg( "", "G", false );
+ this->m_StrArgs[ "slicerw_skeleton" ] = TStrArg( "", "U", false );
+ this->m_StrArgs[ "andrw_skeleton" ] = TStrArg( "", "W", false );
+ this->m_StrArgs[ "fastrw_points" ] = TStrArg( "", "B", false );
+ this->m_StrArgs[ "slicerw_points" ] = TStrArg( "", "C", false );
+ this->m_StrArgs[ "andrw_points" ] = TStrArg( "", "D", false );
+ this->m_StrArgs[ "seed_file" ] = TStrArg( "", "P", false );
+ this->m_StrArgs[ "seed" ] = TStrArg( "", "S", false );
+ this->m_StrArgs[ "seed_type" ] = TStrArg( "point", "T", false );
+
+ this->m_DblArgs[ "beta" ] = TDblArg( 2.5, "b", false );
+ this->m_DblArgs[ "epsilon" ] = TDblArg( 1e-5, "e", false );
+ this->m_DblArgs[ "vesselness_sigma" ] = TDblArg( 0.5, "a", false );
+ this->m_DblArgs[ "vesselness_alpha1" ] = TDblArg( 0.5, "c", false );
+ this->m_DblArgs[ "vesselness_alpha2" ] = TDblArg( 2, "d", false );
+ this->m_DblArgs[ "mori_min_thr" ] = TDblArg( -850, "f", false );
+ this->m_DblArgs[ "mori_kernel" ] = TDblArg( 20, "g", false );
+ this->m_DblArgs[ "mori_signal_thr" ] = TDblArg( 100, "i", false );
+ this->m_DblArgs[ "mori_signal_influence" ] = TDblArg( 0.5, "j", false );
+ this->m_DblArgs[ "mori_lower" ] = TDblArg( -1024, "k", false );
+ this->m_DblArgs[ "mori_upper" ] = TDblArg( 0, "l", false );
+ this->m_DblArgs[ "mori_delta" ] = TDblArg( 1, "m", false );
+ this->m_DblArgs[ "slicerw_thr" ] = TDblArg( 5, "n", false );
+ this->m_DblArgs[ "fastrw_thr" ] = TDblArg( 65, "o", false );
+ this->m_DblArgs[ "labels_upper_thr" ] = TDblArg( -400, "p", false );
}
// -------------------------------------------------------------------------
void CTBronchi::Process::
ParseArguments( int argc, char* argv[] )
{
- typedef TCLAP::ValueArg< std::string > _TStrArg;
- typedef TCLAP::ValueArg< TScalar > _TDblArg;
- std::map< std::string, _TStrArg* > strArg;
- std::map< std::string, _TDblArg* > dblArg;
- unsigned char flag = 'a';
-
- TArguments::iterator sIt = this->m_StrArg.begin( );
- for( ; sIt != this->m_StrArg.end( ); ++sIt )
- {
- std::stringstream fStr;
- fStr << flag;
- flag++;
- if( flag == 'h' )
- flag++;
-
- _TStrArg* a = new _TStrArg(
- fStr.str( ), sIt->first, sIt->first,
+ typedef TCLAP::ValueArg< TString > _TStrArg;
+ typedef TCLAP::ValueArg< TScalar > _TDblArg;
+ std::map< TString, _TStrArg* > strArgs;
+ std::map< TString, _TDblArg* > dblArgs;
+
+ // Load arguments
+ TStrArgs::iterator sIt = this->m_StrArgs.begin( );
+ for( ; sIt != this->m_StrArgs.end( ); ++sIt )
+ strArgs[ sIt->first ] = new _TStrArg(
std::get< 1 >( sIt->second ),
+ sIt->first, sIt->first,
+ std::get< 2 >( sIt->second ),
std::get< 0 >( sIt->second ),
- "string"
+ TString( "string (" ) + std::get< 0 >( sIt->second ) + ")"
);
- strArg[ sIt->first ] = a;
-
- } // rof
- TArguments::iterator dIt = this->m_DblArg.begin( );
- for( ; dIt != this->m_DblArg.end( ); ++dIt )
+ TDblArgs::iterator dIt = this->m_DblArgs.begin( );
+ for( ; dIt != this->m_DblArgs.end( ); ++dIt )
{
- std::stringstream fStr;
- fStr << flag;
- flag++;
- if( flag == 'h' )
- flag++;
-
- std::istringstream vStr( std::get< 0 >( dIt->second ) );
- TScalar v;
- vStr >> v;
-
- _TDblArg* a = new _TDblArg(
- fStr.str( ), dIt->first, dIt->first,
+ std::stringstream v;
+ v << "value (" << std::get< 0 >( dIt->second ) << ")";
+ dblArgs[ dIt->first ] = new _TDblArg(
std::get< 1 >( dIt->second ),
- v, "value"
+ dIt->first, dIt->first,
+ std::get< 2 >( dIt->second ),
+ std::get< 0 >( dIt->second ),
+ v.str( )
);
- dblArg[ dIt->first ] = a;
} // rof
- std::map< std::string, _TStrArg* >::iterator saIt;
- std::map< std::string, _TDblArg* >::iterator daIt;
+ // Parse arguments
+ std::map< TString, _TStrArg* >::iterator saIt;
+ std::map< TString, _TDblArg* >::iterator daIt;
try
{
TCLAP::CmdLine cmd( "CTBronchi pipeline", ' ', "1.0.0" );
- for( saIt = strArg.begin( ); saIt != strArg.end( ); ++saIt )
+ for( saIt = strArgs.begin( ); saIt != strArgs.end( ); ++saIt )
cmd.add( *( saIt->second ) );
- for( daIt = dblArg.begin( ); daIt != dblArg.end( ); ++daIt )
+ for( daIt = dblArgs.begin( ); daIt != dblArgs.end( ); ++daIt )
cmd.add( *( daIt->second ) );
cmd.parse( argc, argv );
}
} // yrt
- for( saIt = strArg.begin( ); saIt != strArg.end( ); ++saIt )
+ // Get values
+ for( saIt = strArgs.begin( ); saIt != strArgs.end( ); ++saIt )
{
- std::get< 0 >( this->m_StrArg[ saIt->first ] ) = saIt->second->getValue( );
+ std::get< 0 >( this->m_StrArgs[ saIt->first ] ) =
+ saIt->second->getValue( );
delete saIt->second;
} // rof
- for( daIt = dblArg.begin( ); daIt != dblArg.end( ); ++daIt )
+ for( daIt = dblArgs.begin( ); daIt != dblArgs.end( ); ++daIt )
{
- std::stringstream vStr;
- vStr << daIt->second->getValue( );
- std::get< 0 >( this->m_DblArg[ daIt->first ] ) = vStr.str( );
+ std::get< 0 >( this->m_DblArgs[ daIt->first ] ) =
+ daIt->second->getValue( );
delete daIt->second;
} // rof
- // Compute base filename
- std::string iname = std::get< 0 >( this->m_StrArg[ "input" ] );
- this->m_BaseFileName = iname.substr( 0, iname.find_last_of( "." ) );
+ // Compute base filenames
+ TString bname = std::get< 0 >( this->m_StrArgs[ "input" ] );
+ bname = bname.substr( 0, bname.find_last_of( "." ) );
+ const unsigned int N = 6;
+ const unsigned int M = 6;
+ TString names[ N + M ] =
+ {
+ "vesselness",
+ "mori",
+ "labels",
+ "fastrw",
+ "slicerw",
+ "andrw",
+ "fastrw_skeleton",
+ "slicerw_skeleton",
+ "andrw_skeleton",
+ "fastrw_points",
+ "slicerw_points",
+ "andrw_points"
+ };
+ for( unsigned int i = 0; i < N; ++i )
+ if( std::get< 0 >( this->m_StrArgs[ names[ i ] ] ) == "" )
+ std::get< 0 >( this->m_StrArgs[ names[ i ] ] ) =
+ bname + "_" + names[ i ] + ".mha";
+ for( unsigned int i = 0; i < M; ++i )
+ if( std::get< 0 >( this->m_StrArgs[ names[ i + N ] ] ) == "" )
+ std::get< 0 >( this->m_StrArgs[ names[ i + N ] ] ) =
+ bname + "_" + names[ i + N ] + ".txt";
}
// -------------------------------------------------------------------------
void CTBronchi::Process::
Update( )
{
- this->_Input( );
- this->_Vesselness( );
- this->_Mori( );
- this->_MoriLabelling( );
+ this->_Input( this->m_Input );
+ this->_Vesselness( this->m_Input, this->m_Vesselness );
+ this->_Mori( this->m_Input, this->m_Mori, this->m_Seed );
+ this->_MoriLabelling(
+ this->m_Input, this->m_Mori, this->m_Vesselness, this->m_Labels
+ );
+ this->_FastRW( this->m_Input, this->m_Labels, this->m_FastRW );
+ this->_SliceRW(
+ this->m_Input, this->m_Mori, this->m_Vesselness, this->m_SliceRW
+ );
+ this->_AndImages( this->m_FastRW, this->m_SliceRW, this->m_AndRW );
+ this->_Skeleton(
+ this->m_SliceRW, this->m_SliceRWSkeleton,
+ std::get< 0 >( this->m_StrArgs[ "slicerw_skeleton" ] )
+ );
+ this->_Skeleton(
+ this->m_AndRW, this->m_AndRWSkeleton,
+ std::get< 0 >( this->m_StrArgs[ "andrw_skeleton" ] )
+ );
+ this->_Skeleton(
+ this->m_FastRW, this->m_FastRWSkeleton,
+ std::get< 0 >( this->m_StrArgs[ "fastrw_skeleton" ] )
+ );
}
// -------------------------------------------------------------------------
+template< class _TImage >
void CTBronchi::Process::
-_Input( )
+_Input( _TImage& input )
{
- std::string iname = std::get< 0 >( this->m_StrArg[ "input" ] );
- std::string sname = std::get< 0 >( this->m_StrArg[ "seed_file" ] );
- std::string seed = std::get< 0 >( this->m_StrArg[ "seed" ] );
- std::string stype = std::get< 0 >( this->m_StrArg[ "seed_type" ] );
+ TString sname = std::get< 0 >( this->m_StrArgs[ "seed_file" ] );
+ TString seed = std::get< 0 >( this->m_StrArgs[ "seed" ] );
+ TString stype = std::get< 0 >( this->m_StrArgs[ "seed_type" ] );
if( sname == "" && seed == "" )
throw std::runtime_error( "No seed given." );
if( sname != "" )
sSeed >> idx[ 0 ] >> idx[ 1 ] >> idx[ 2 ];
// Read image
- double t = this->m_Input.Load( iname );
+ TString iname = std::get< 0 >( this->m_StrArgs[ "input" ] );
+ double t = input.Load( iname );
if( t < 0 )
- std::runtime_error( std::string( "Could not read \"" ) + iname + "\"" );
+ std::runtime_error( TString( "Could not read \"" ) + iname + "\"" );
std::cout << "\"" << iname << "\" read in " << t << " s." << std::endl;
// Synch seed
if( stype == "point" )
- this->m_Input.Get( )->TransformPhysicalPointToIndex( pnt, idx );
+ input.Get( )->TransformPhysicalPointToIndex( pnt, idx );
else
- this->m_Input.Get( )->TransformIndexToPhysicalPoint( idx, pnt );
+ input.Get( )->TransformIndexToPhysicalPoint( idx, pnt );
this->m_Seed = TSeed( idx, pnt );
}
// -------------------------------------------------------------------------
+template< class _TInput, class _TVesselness >
void CTBronchi::Process::
-_Vesselness( )
+_Vesselness( _TInput& input, _TVesselness& vesselness )
{
- std::string vname = this->m_BaseFileName + "_vesselness.mha";
- double t = this->m_Vesselness.Load( vname );
+ TString iname = std::get< 0 >( this->m_StrArgs[ "vesselness" ] );
+ double t = vesselness.Load( iname );
if( t < 0 )
{
t = 0;
- // Arguments
- std::stringstream v;
- v << std::get< 0 >( this->m_DblArg[ "vesselness_sigma" ] ) << " "
- << std::get< 0 >( this->m_DblArg[ "vesselness_alpha1" ] ) << " "
- << std::get< 0 >( this->m_DblArg[ "vesselness_alpha2" ] );
- std::istringstream w( v.str( ) );
- TScalar s, a1, a2;
- w >> s >> a1 >> a2;
-
// Min-max
- typedef itk::MinimumMaximumImageCalculator< TPixelImage::TImage > _TMinMax;
- _TMinMax::Pointer minMax = _TMinMax::New( );
- minMax->SetImage( this->m_Input.Get( ) );
+ typedef itk::MinimumMaximumImageCalculator< typename _TInput::TImage > _TMinMax;
+ typename _TMinMax::Pointer minMax = _TMinMax::New( );
+ minMax->SetImage( input.Get( ) );
minMax->Compute( );
// Invert intensities
- typedef CTBronchi::Filter< itk::InvertIntensityImageFilter< TPixelImage::TImage > > _TInverter;
+ typedef CTBronchi::Filter< itk::InvertIntensityImageFilter< typename _TInput::TImage > > _TInverter;
_TInverter inverter;
- inverter.Get( )->SetInput( this->m_Input.Get( ) );
+ inverter.Get( )->SetInput( input.Get( ) );
inverter.Get( )->SetMaximum( minMax->GetMaximum( ) );
double t1 = inverter.Update( );
t += t1;
std::cout << "Inversion executed in " << t1 << " s" << std::endl;
// Compute hessian image
- typedef CTBronchi::Filter< itk::HessianRecursiveGaussianImageFilter< TPixelImage::TImage > > _THessian;
+ typedef CTBronchi::Filter< itk::HessianRecursiveGaussianImageFilter< typename _TInput::TImage > > _THessian;
_THessian hessian;
hessian.Get( )->SetInput( inverter.Get( )->GetOutput( ) );
- hessian.Get( )->SetSigma( s );
+ hessian.Get( )->SetSigma(
+ std::get< 0 >( this->m_DblArgs[ "vesselness_sigma" ] )
+ );
t1 = hessian.Update( );
t += t1;
std::cout << "Hessian executed in " << t1 << " s" << std::endl;
// Vesselness
- typedef CTBronchi::Filter< itk::Hessian3DToVesselnessMeasureImageFilter< TScalar > > _TVesselness;
- _TVesselness vesselness;
- vesselness.Get( )->SetInput( hessian.Get( )->GetOutput( ) );
- vesselness.Get( )->SetAlpha1( a1 );
- vesselness.Get( )->SetAlpha2( a2 );
- t1 = vesselness.Update( );
+ typedef CTBronchi::Filter< itk::Hessian3DToVesselnessMeasureImageFilter< typename _TVesselness::TImage::PixelType > > _TVesselnessFilter;
+ _TVesselnessFilter vFilter;
+ vFilter.Get( )->SetInput( hessian.Get( )->GetOutput( ) );
+ vFilter.Get( )->SetAlpha1(
+ std::get< 0 >( this->m_DblArgs[ "vesselness_alpha1" ] )
+ );
+ vFilter.Get( )->SetAlpha2(
+ std::get< 0 >( this->m_DblArgs[ "vesselness_alpha2" ] )
+ );
+ t1 = vFilter.Update( );
t += t1;
std::cout << "Hessian measure computed in " << t1 << " s." << std::endl;
- this->m_Vesselness.Set( vesselness.Get( )->GetOutput( ) );
- t1 = this->m_Vesselness.Save( vname );
+ vesselness.Set( vFilter.Get( )->GetOutput( ) );
+ t1 = vesselness.Save( iname );
t += t1;
- std::cout << "\"" << vname << "\" saved in " << t1 << " s." << std::endl;
+ std::cout << "\"" << iname << "\" saved in " << t1 << " s." << std::endl;
} // fi
std::cout << "Vesselness computed in " << t << " s." << std::endl;
}
// -------------------------------------------------------------------------
+template< class _TInput, class _TMori, class _TSeed >
void CTBronchi::Process::
-_Mori( )
+_Mori( _TInput& input, _TMori& mori, _TSeed& seed )
{
- std::string mname = this->m_BaseFileName + "_mori.mha";
- double t = this->m_Mori.Load( mname );
+ TString iname = std::get< 0 >( this->m_StrArgs[ "mori" ] );
+ double t = mori.Load( iname );
if( t < 0 )
{
t = 0;
- // Arguments
- std::stringstream v;
- v << std::get< 0 >( this->m_DblArg[ "mori_min_thr" ] ) << " "
- << std::get< 0 >( this->m_DblArg[ "mori_kernel" ] ) << " "
- << std::get< 0 >( this->m_DblArg[ "mori_signal_thr" ] ) << " "
- << std::get< 0 >( this->m_DblArg[ "mori_signal_influence" ] ) << " "
- << std::get< 0 >( this->m_DblArg[ "mori_lower" ] ) << " "
- << std::get< 0 >( this->m_DblArg[ "mori_upper" ] ) << " "
- << std::get< 0 >( this->m_DblArg[ "mori_delta" ] );
- std::istringstream w( v.str( ) );
- TScalar mThr, sKernel, sThr, sInfluence, lower, upper, delta;
- w >> mThr >> sKernel >> sThr >> sInfluence >> lower >> upper >> delta;
-
// Mori segmentation
- typedef CTBronchi::Filter< fpa::Filters::Image::Mori< TPixelImage::TImage, TLabelImage::TImage > > _TMori;
- _TMori mori;
- mori.Get( )->SetInput( this->m_Input.Get( ) );
- mori.Get( )->SetSeed( this->m_Seed.second );
- mori.Get( )->SetInsideValue( this->m_InsideLabel );
- mori.Get( )->SetOutsideValue( this->m_UndefinedLabel );
- mori.Get( )->SetMinimumThreshold( mThr );
- mori.Get( )->SetSignalKernelSize( sKernel );
- mori.Get( )->SetSignalThreshold( sThr );
- mori.Get( )->SetSignalInfluence( sInfluence );
- mori.Get( )->SetThresholds( lower, upper, delta );
- double t1 = mori.Update( );
+ typedef CTBronchi::Filter< fpa::Filters::Image::Mori< TPixelImage::TImage, TLabelImage::TImage > > _TMoriFilter;
+ _TMoriFilter mFilter;
+ mFilter.Get( )->SetInput( input.Get( ) );
+ mFilter.Get( )->SetSeed( seed.second );
+ mFilter.Get( )->SetInsideValue( this->m_InsideLabel );
+ mFilter.Get( )->SetOutsideValue( this->m_UndefinedLabel );
+ mFilter.Get( )->SetMinimumThreshold(
+ std::get< 0 >( this->m_DblArgs[ "mori_min_thr" ] )
+ );
+ mFilter.Get( )->SetSignalKernelSize(
+ std::get< 0 >( this->m_DblArgs[ "mori_kernel" ] )
+ );
+ mFilter.Get( )->SetSignalThreshold(
+ std::get< 0 >( this->m_DblArgs[ "mori_signal_thr" ] )
+ );
+ mFilter.Get( )->SetSignalInfluence(
+ std::get< 0 >( this->m_DblArgs[ "mori_signal_influence" ] )
+ );
+ mFilter.Get( )->SetThresholds(
+ std::get< 0 >( this->m_DblArgs[ "mori_lower" ] ),
+ std::get< 0 >( this->m_DblArgs[ "mori_upper" ] ),
+ std::get< 0 >( this->m_DblArgs[ "mori_delta" ] )
+ );
+ double t1 = mFilter.Update( );
t += t1;
std::cout << "Segmentation computed in " << t1 << " s." << std::endl;
- this->m_Mori.Set( mori.Get( )->GetOutput( ) );
- t1 = this->m_Mori.Save( mname );
+ mori.Set( mFilter.Get( )->GetOutput( ) );
+ t1 = mori.Save( iname );
t += t1;
- std::cout << "\"" << mname << "\" saved in " << t1 << " s." << std::endl;
+ std::cout << "\"" << iname << "\" saved in " << t1 << " s." << std::endl;
} // fi
std::cout << "Mori segmentation computed in " << t << " s." << std::endl;
}
// -------------------------------------------------------------------------
+template< class _TInput, class _TMori, class _TVesselness, class _TLabels >
void CTBronchi::Process::
-_MoriLabelling( )
+_MoriLabelling(
+ _TInput& input, _TMori& mori, _TVesselness& vesselness, _TLabels& labels
+ )
{
- std::string mname = this->m_BaseFileName + "_labels.mha";
- double t = this->m_Labels.Load( mname );
+ TString iname = std::get< 0 >( this->m_StrArgs[ "labels" ] );
+ double t = labels.Load( iname );
if( t < 0 )
{
t = 0;
- // Arguments
- std::stringstream v;
- v << std::get< 0 >( this->m_DblArg[ "fastrw_thr" ] ) << " "
- << std::get< 0 >( this->m_DblArg[ "labels_upper_thr" ] );
- std::istringstream w( v.str( ) );
- TScalar vThr, uThr;
- w >> vThr >> uThr;
-
// Mori labelling
- typedef CTBronchi::Filter< CTBronchi::MoriLabelling< TPixelImage::TImage, TLabelImage::TImage, TScalarImage::TImage > > _TLabelling;
+ typedef CTBronchi::Filter< CTBronchi::MoriLabelling< typename _TInput::TImage, typename _TLabels::TImage, typename _TVesselness::TImage > > _TLabelling;
_TLabelling labelling;
- labelling.Get( )->SetInput( this->m_Input.Get( ) );
- labelling.Get( )->SetInputLabels( this->m_Mori.Get( ) );
- labelling.Get( )->SetInputVesselness( this->m_Vesselness.Get( ) );
- labelling.Get( )->SetVesselnessThreshold( vThr );
- labelling.Get( )->SetUpperThreshold( uThr );
+ labelling.Get( )->SetInput( input.Get( ) );
+ labelling.Get( )->SetInputLabels( mori.Get( ) );
+ labelling.Get( )->SetInputVesselness( vesselness.Get( ) );
+ labelling.Get( )->SetVesselnessThreshold( std::get< 0 >( this->m_DblArgs[ "fastrw_thr" ] ) );
+ labelling.Get( )->SetUpperThreshold( std::get< 0 >( this->m_DblArgs[ "labels_upper_thr" ] ) );
labelling.Get( )->SetInsideValue( this->m_InsideLabel );
labelling.Get( )->SetOutsideValue( this->m_OutsideLabel );
labelling.Get( )->SetFillValue( this->m_OutsideLabel );
t += t1;
std::cout << "Labelling computed in " << t1 << " s." << std::endl;
- this->m_Labels.Set( labelling.Get( )->GetOutput( ) );
- t1 = this->m_Labels.Save( mname );
+ labels.Set( labelling.Get( )->GetOutput( ) );
+ t1 = labels.Save( iname );
t += t1;
- std::cout << "\"" << mname << "\" saved in " << t1 << " s." << std::endl;
+ std::cout << "\"" << iname << "\" saved in " << t1 << " s." << std::endl;
} // fi
std::cout << "Mori labelling computed in " << t << " s." << std::endl;
}
+// -------------------------------------------------------------------------
+template< class _TInput, class _TLabels, class _TFastRW >
+void CTBronchi::Process::
+_FastRW( _TInput& input, _TLabels& labels, _TFastRW& fastrw )
+{
+ TString iname = std::get< 0 >( this->m_StrArgs[ "fastrw" ] );
+ double t = fastrw.Load( iname );
+ if( t < 0 )
+ {
+ t = 0;
+
+ // Prepare weight functor
+ typedef fpa::Functors::Dijkstra::Image::Gaussian< typename _TInput::TImage, typename _TFastRW::TImage::PixelType > TWeight;
+ typename TWeight::Pointer weight = TWeight::New( );
+ weight->SetBeta( std::get< 0 >( this->m_DblArgs[ "beta" ] ) );
+ weight->SetEpsilon( std::get< 0 >( this->m_DblArgs[ "epsilon" ] ) );
+
+ // Random walk
+ typedef CTBronchi::Filter< fpa::Filters::Image::RandomWalker< typename _TInput::TImage, typename _TLabels::TImage, typename _TFastRW::TImage::PixelType > > TFilter;
+ TFilter filter;
+ filter.Get( )->SetInputImage( input.Get( ) );
+ filter.Get( )->SetInputLabels( labels.Get( ) );
+ filter.Get( )->SetWeightFunction( weight );
+ double t1 = filter.Update( );
+ t += t1;
+ std::cout << "Fast random walker executed in " << t1 << " s" << std::endl;
+
+ // Extract label
+ typedef CTBronchi::Filter< itk::BinaryThresholdImageFilter< typename _TLabels::TImage, typename _TLabels::TImage > > _TExtract;
+ _TExtract extract;
+ extract.Get( )->SetInput( filter.Get( )->GetOutputLabels( ) );
+ extract.Get( )->SetInsideValue( this->m_InsideLabel );
+ extract.Get( )->SetOutsideValue( this->m_UndefinedLabel );
+ extract.Get( )->SetLowerThreshold( this->m_InsideLabel );
+ extract.Get( )->SetUpperThreshold( this->m_InsideLabel );
+ t1 = extract.Update( );
+ t += t1;
+ std::cout << "Extract labels executed in " << t1 << " s" << std::endl;
+
+ fastrw.Set( extract.Get( )->GetOutput( ) );
+ t1 = fastrw.Save( iname );
+ t += t1;
+ std::cout << "\"" << iname << "\" saved in " << t1 << " s." << std::endl;
+
+ } // fi
+ std::cout << "Fast random walker computed in " << t << " s." << std::endl;
+}
+
+// -------------------------------------------------------------------------
+template< class _TInput, class _TLabels, class _TVesselness, class _TSliceRW >
+void CTBronchi::Process::
+_SliceRW(
+ _TInput& input, _TLabels& labels, _TVesselness& vesselness, _TSliceRW& slicerw
+ )
+{
+ TString iname = std::get< 0 >( this->m_StrArgs[ "slicerw" ] );
+ double t = slicerw.Load( iname );
+ if( t < 0 )
+ {
+ t = 0;
+
+ // Random walk
+ typedef CTBronchi::Filter< fpa::Common::SliceBySliceRandomWalker< typename _TInput::TImage, typename _TLabels::TImage, typename _TVesselness::TImage > > TFilter;
+ TFilter filter;
+ filter.Get( )->SetInput( input.Get( ) );
+ filter.Get( )->SetInputLabels( labels.Get( ) );
+ filter.Get( )->SetInputVesselness( vesselness.Get( ) );
+ filter.Get( )->SetBeta( std::get< 0 >( this->m_DblArgs[ "beta" ] ) );
+ filter.Get( )->SetVesselnessThreshold( std::get< 0 >( this->m_DblArgs[ "slicerw_thr" ] ) );
+ filter.Get( )->SetEpsilon( std::get< 0 >( this->m_DblArgs[ "epsilon" ] ) );
+ double t1 = filter.Update( );
+ t += t1;
+ std::cout << "Extract labels executed in " << t1 << " s" << std::endl;
+
+ slicerw.Set( filter.Get( )->GetOutput( ) );
+ t1 = slicerw.Save( iname );
+ t += t1;
+ std::cout << "\"" << iname << "\" saved in " << t1 << " s." << std::endl;
+
+ } // fi
+ std::cout << "Slice by slice random walker computed in " << t << " s." << std::endl;
+}
+
+// -------------------------------------------------------------------------
+template< class _TInput >
+void CTBronchi::Process::
+_AndImages( _TInput& a, _TInput& b, _TInput& c )
+{
+ TString iname = std::get< 0 >( this->m_StrArgs[ "andrw" ] );
+ double t = c.Load( iname );
+ if( t < 0 )
+ {
+ t = 0;
+
+ // Random walk
+ typedef CTBronchi::Filter< itk::AndImageFilter< typename _TInput::TImage > > TFilter;
+ TFilter filter;
+ filter.Get( )->SetInput( 0, a.Get( ) );
+ filter.Get( )->SetInput( 1, b.Get( ) );
+ double t1 = filter.Update( );
+ t += t1;
+ std::cout << "And segmentations executed in " << t1 << " s" << std::endl;
+
+ c.Set( filter.Get( )->GetOutput( ) );
+ t1 = c.Save( iname );
+ t += t1;
+ std::cout << "\"" << iname << "\" saved in " << t1 << " s." << std::endl;
+
+ } // fi
+ std::cout << "And segmentations computed in " << t << " s." << std::endl;
+}
+
+// -------------------------------------------------------------------------
+template< class _TInput, class _TSkeleton >
+void CTBronchi::Process::
+_Skeleton( _TInput& a, _TSkeleton& s, const TString& fname )
+{
+ double t = s.Load( fname );
+ if( t < 0 )
+ {
+ t = 0;
+
+ typedef CTBronchi::Filter< fpa::Filters::Image::ExtractSkeleton< typename _TInput::TImage > > TFilter;
+ TFilter filter;
+ filter.Get( )->SetInput( a.Get( ) );
+ filter.Get( )->SeedFromMaximumDistanceOff( );
+ filter.Get( )->SetSeed( this->m_Seed.first );
+ filter.Get( )->GetDistanceMap( )->InsideIsPositiveOn( );
+ filter.Get( )->GetDistanceMap( )->SquaredDistanceOff( );
+ filter.Get( )->GetDistanceMap( )->UseImageSpacingOn( );
+ double t1 = filter.Update( );
+ t += t1;
+ std::cout << "Skeleton executed in " << t1 << " s" << std::endl;
+
+ s.Set( filter.Get( )->GetOutput( ) );
+ t1 = s.Save( fname );
+ t += t1;
+ std::cout << "\"" << fname << "\" saved in " << t1 << " s." << std::endl;
+
+ } // fi
+ std::cout << "Skeleton computed in " << t << " s." << std::endl;
+}
+
// eof - $RCSfile$
#include <tuple>
#include <fpa_ctbronchi_export.h>
#include "Image.h"
+#include "Skeleton.h"
namespace CTBronchi
{
typedef short TPixel;
typedef unsigned char TLabel;
typedef float TScalar;
+ typedef std::string TString;
// Arguments
- typedef std::tuple< std::string, bool > TArgument;
- typedef std::map< std::string, TArgument > TArguments;
+ typedef std::tuple< TString, TString, bool > TStrArg;
+ typedef std::tuple< TScalar, TString, bool > TDblArg;
+ typedef std::map< TString, TStrArg > TStrArgs;
+ typedef std::map< TString, TDblArg > TDblArgs;
// Images
typedef CTBronchi::Image< TPixel, Dim > TPixelImage;
typedef CTBronchi::Image< TLabel, Dim > TLabelImage;
typedef CTBronchi::Image< TScalar, Dim > TScalarImage;
+ // Skeleton
+ typedef CTBronchi::Skeleton< Dim > TSkeleton;
+
// Seed
typedef TPixelImage::TImage::PointType TPoint;
typedef TPixelImage::TImage::IndexType TIndex;
void Update( );
protected:
- void _Input( );
- void _Vesselness( );
- void _Mori( );
- void _MoriLabelling( );
+ template< class _TImage >
+ void _Input( _TImage& input );
+
+ template< class _TInput, class _TVesselness >
+ void _Vesselness( _TInput& input, _TVesselness& vesselness );
+
+ template< class _TInput, class _TMori, class _TSeed >
+ void _Mori( _TInput& input, _TMori& mori, _TSeed& seed );
+
+ template< class _TInput, class _TMori, class _TVesselness, class _TLabels >
+ void _MoriLabelling(
+ _TInput& input, _TMori& mori, _TVesselness& vesselness, _TLabels& labels
+ );
+
+ template< class _TInput, class _TLabels, class _TFastRW >
+ void _FastRW( _TInput& input, _TLabels& labels, _TFastRW& fastrw );
+
+ template< class _TInput, class _TLabels, class _TVesselness, class _TSliceRW >
+ void _SliceRW(
+ _TInput& input, _TLabels& labels, _TVesselness& vesselness, _TSliceRW& slicerw
+ );
+
+ template< class _TInput >
+ void _AndImages( _TInput& a, _TInput& b, _TInput& c );
+
+ template< class _TInput, class _TSkeleton >
+ void _Skeleton( _TInput& a, _TSkeleton& s, const TString& fname );
protected:
- TArguments m_StrArg;
- TArguments m_DblArg;
+ TStrArgs m_StrArgs;
+ TDblArgs m_DblArgs;
- std::string m_BaseFileName;
- TSeed m_Seed;
+ TSeed m_Seed;
TLabel m_UndefinedLabel;
TLabel m_InsideLabel;
TLabel m_OutsideLabel;
TLabelImage m_Labels;
TLabelImage m_FastRW;
TLabelImage m_SliceRW;
+ TLabelImage m_AndRW;
+
+ TSkeleton m_FastRWSkeleton;
+ TSkeleton m_SliceRWSkeleton;
+ TSkeleton m_AndRWSkeleton;
};
} // ecapseman
-/* TODO
- #include <string>
- #include <tclap/CmdLine.h>
- #include <itkImage.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 itk::NumericTraits< TPixel >::RealType TScalar;
- typedef itk::Image< TPixel, Dim > TImage;
- typedef itk::Image< TScalar, Dim > TScalarImage;
-
- // -------------------------------------------------------------------------
- int main( int argc, char* argv[] )
- {
- 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
- {
- 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
- << "===============================" << std::endl
- << "Error caught: " << std::endl
- << err.error( ) << " " << err.argId( )
- << "===============================" << std::endl
- << std::endl;
- return( 1 );
-
- } // yrt
-
- try
- {
- // 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 )
- {
- std::cerr
- << "===============================" << std::endl
- << "Error caught: " << std::endl
- << err.what( )
- << "===============================" << std::endl
- << std::endl;
- return( 1 );
-
- } // yrt
- return( 0 );
- }
-*/
#endif // __CTBronchi__Process__h__
// eof - $RCSfile$
git://git.creatis.insa-lyon.fr / FrontAlgorithms.git / commitdiff