From 3d7baeb20fa756d8f57654b7ff75d3be03ee1c1f Mon Sep 17 00:00:00 2001 From: =?utf8?q?Leonardo=20Fl=C3=B3rez-Valencia?= Date: Wed, 6 Dec 2017 15:14:38 -0500 Subject: [PATCH] ... --- appli/CTBronchi/Process.cxx | 509 +++++++++++++++------- appli/CTBronchi/Process.h | 161 ++----- appli/CTBronchi/Skeleton.h | 45 ++ appli/CTBronchi/Skeleton.hxx | 130 ++++++ lib/fpa/Filters/Image/ExtractSkeleton.hxx | 3 - 5 files changed, 566 insertions(+), 282 deletions(-) create mode 100644 appli/CTBronchi/Skeleton.h create mode 100644 appli/CTBronchi/Skeleton.hxx diff --git a/appli/CTBronchi/Process.cxx b/appli/CTBronchi/Process.cxx index 9f29ffc..d1498ea 100644 --- a/appli/CTBronchi/Process.cxx +++ b/appli/CTBronchi/Process.cxx @@ -6,12 +6,18 @@ #include #include +#include +#include #include #include #include #include +#include +#include #include +#include +#include #include "MoriLabelling.h" #include "Process.h" @@ -23,26 +29,38 @@ Process( ) 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 ); } // ------------------------------------------------------------------------- @@ -55,61 +73,46 @@ CTBronchi::Process:: 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 ); } @@ -121,44 +124,89 @@ ParseArguments( int argc, char* 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 != "" ) @@ -185,155 +233,149 @@ _Input( ) 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 ); @@ -341,13 +383,156 @@ _MoriLabelling( ) 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$ diff --git a/appli/CTBronchi/Process.h b/appli/CTBronchi/Process.h index bae6f43..d9b939b 100644 --- a/appli/CTBronchi/Process.h +++ b/appli/CTBronchi/Process.h @@ -8,6 +8,7 @@ #include #include #include "Image.h" +#include "Skeleton.h" namespace CTBronchi { @@ -19,16 +20,22 @@ 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; @@ -42,17 +49,39 @@ namespace CTBronchi 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; @@ -63,117 +92,15 @@ namespace CTBronchi TLabelImage m_Labels; TLabelImage m_FastRW; TLabelImage m_SliceRW; + TLabelImage m_AndRW; + + TSkeleton m_FastRWSkeleton; + TSkeleton m_SliceRWSkeleton; + TSkeleton m_AndRWSkeleton; }; } // ecapseman -/* TODO - #include - #include - #include - #include - #include - #include - #include - #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$ diff --git a/appli/CTBronchi/Skeleton.h b/appli/CTBronchi/Skeleton.h new file mode 100644 index 0000000..50893d5 --- /dev/null +++ b/appli/CTBronchi/Skeleton.h @@ -0,0 +1,45 @@ +// ========================================================================= +// @author Leonardo Florez Valencia (florez-l@javeriana.edu.co) +// ========================================================================= +#ifndef __CTBronchi__Skeleton__h__ +#define __CTBronchi__Skeleton__h__ + +#include + +namespace CTBronchi +{ + /** + */ + template< unsigned int _VDim > + class Skeleton + { + public: + static const unsigned int VDim = _VDim; + typedef ivq::ITK::ImageSkeleton< VDim > TSkeleton; + + public: + Skeleton( ); + virtual ~Skeleton( ); + + bool IsNotNull( ) const; + bool IsNull( ) const; + + TSkeleton* Get( ); + const TSkeleton* Get( ) const; + void Set( TSkeleton* sk ); + void Set( const typename TSkeleton::Pointer& sk ); + + double Load( const std::string& fname ); + double Save( const std::string& fname ); + + protected: + typename TSkeleton::Pointer m_Skeleton; + }; + +} // ecapseman + +#include "Skeleton.hxx" + +#endif // __CTBronchi__Skeleton__h__ + +// eof - $RCSfile$ diff --git a/appli/CTBronchi/Skeleton.hxx b/appli/CTBronchi/Skeleton.hxx new file mode 100644 index 0000000..2b293e5 --- /dev/null +++ b/appli/CTBronchi/Skeleton.hxx @@ -0,0 +1,130 @@ +// ========================================================================= +// @author Leonardo Florez Valencia (florez-l@javeriana.edu.co) +// ========================================================================= +#ifndef __CTBronchi__Skeleton__hxx__ +#define __CTBronchi__Skeleton__hxx__ + +#include "Filter.h" +#include +#include + +// ------------------------------------------------------------------------- +template< unsigned int _VDim > +CTBronchi::Skeleton< _VDim >:: +Skeleton( ) +{ +} + +// ------------------------------------------------------------------------- +template< unsigned int _VDim > +CTBronchi::Skeleton< _VDim >:: +~Skeleton( ) +{ +} + +// ------------------------------------------------------------------------- +template< unsigned int _VDim > +bool CTBronchi::Skeleton< _VDim >:: +IsNotNull( ) const +{ + return( this->m_Skeleton.IsNotNull( ) ); +} + +// ------------------------------------------------------------------------- +template< unsigned int _VDim > +bool CTBronchi::Skeleton< _VDim >:: +IsNull( ) const +{ + return( this->m_Skeleton.IsNull( ) ); +} + +// ------------------------------------------------------------------------- +template< unsigned int _VDim > +typename CTBronchi::Skeleton< _VDim >:: +TSkeleton* CTBronchi::Skeleton< _VDim >:: +Get( ) +{ + return( this->m_Skeleton.GetPointer( ) ); +} + +// ------------------------------------------------------------------------- +template< unsigned int _VDim > +const typename CTBronchi::Skeleton< _VDim >:: +TSkeleton* CTBronchi::Skeleton< _VDim >:: +Get( ) const +{ + return( this->m_Skeleton.GetPointer( ) ); +} + +// ------------------------------------------------------------------------- +template< unsigned int _VDim > +void CTBronchi::Skeleton< _VDim >:: +Set( TSkeleton* sk ) +{ + this->m_Skeleton = sk; + /* TODO + if( this->m_Skeleton.IsNotNull( ) ) + this->m_Skeleton->DisconnectPipeline( ); + */ +} + +// ------------------------------------------------------------------------- +template< unsigned int _VDim > +void CTBronchi::Skeleton< _VDim >:: +Set( const typename TSkeleton::Pointer& sk ) +{ + this->m_Skeleton = sk.GetPointer( ); + /* TODO + if( this->m_Skeleton.IsNotNull( ) ) + this->m_Skeleton->DisconnectPipeline( ); + */ +} + +// ------------------------------------------------------------------------- +template< unsigned int _VDim > +double CTBronchi::Skeleton< _VDim >:: +Load( const std::string& fname ) +{ + typedef CTBronchi::Filter< ivq::ITK::ImageSkeletonReader< TSkeleton > > _TReader; + _TReader r; + r.Get( )->SetFileName( fname ); + try + { + double t = r.Update( ); + this->Set( r.Get( )->GetOutput( ) ); + return( t ); + } + catch( ... ) + { + this->Set( NULL ); + return( -1 ); + + } // yrt +} + +// ------------------------------------------------------------------------- +template< unsigned int _VDim > +double CTBronchi::Skeleton< _VDim >:: +Save( const std::string& fname ) +{ + typedef CTBronchi::Filter< ivq::ITK::ImageSkeletonWriter< TSkeleton > > _TWriter; + + double t = -1; + if( this->IsNotNull( ) ) + { + _TWriter w; + w.Get( )->SetFileName( fname ); + w.Get( )->SetInput( this->m_Skeleton ); + try + { + t = w.Update( ); + } + catch( ... ) { } + + } // fi + return( t ); +} + +#endif // __CTBronchi__Skeleton__hxx__ + +// eof - $RCSfile$ diff --git a/lib/fpa/Filters/Image/ExtractSkeleton.hxx b/lib/fpa/Filters/Image/ExtractSkeleton.hxx index 9d2dec0..cffdb0e 100644 --- a/lib/fpa/Filters/Image/ExtractSkeleton.hxx +++ b/lib/fpa/Filters/Image/ExtractSkeleton.hxx @@ -159,9 +159,6 @@ GenerateData( ) } // fi - std::cout << int( this->GetInput( )->GetPixel( this->m_Seed ) ) << std::endl; - std::cout << this->m_DistanceMap->GetOutput( )->GetPixel( this->m_Seed ) << std::endl; - // Compute MST typename _TDijkstra::Pointer dijkstra = _TDijkstra::New( ); dijkstra->SetInput( this->m_DistanceMap->GetOutput( ) ); -- 2.47.1