+// =========================================================================
+// @author Leonardo Florez Valencia (florez-l@javeriana.edu.co)
+// =========================================================================
+
+#include <fstream>
+#include <streambuf>
+#include <tclap/CmdLine.h>
+
+#include <itkMinimumMaximumImageCalculator.h>
+#include <itkInvertIntensityImageFilter.h>
+#include <itkHessianRecursiveGaussianImageFilter.h>
+#include <itkHessian3DToVesselnessMeasureImageFilter.h>
+
+#include <fpa/Filters/Image/Mori.h>
+
+#include "MoriLabelling.h"
+#include "Process.h"
+
+// -------------------------------------------------------------------------
+CTBronchi::Process::
+Process( )
+ : m_UndefinedLabel( 0 ),
+ 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 );
+}
+
+// -------------------------------------------------------------------------
+CTBronchi::Process::
+~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,
+ std::get< 1 >( sIt->second ),
+ std::get< 0 >( sIt->second ),
+ "string"
+ );
+ strArg[ sIt->first ] = a;
+
+ } // rof
+
+ TArguments::iterator dIt = this->m_DblArg.begin( );
+ for( ; dIt != this->m_DblArg.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::get< 1 >( dIt->second ),
+ v, "value"
+ );
+ dblArg[ dIt->first ] = a;
+
+ } // rof
+
+ std::map< std::string, _TStrArg* >::iterator saIt;
+ std::map< std::string, _TDblArg* >::iterator daIt;
+ try
+ {
+ TCLAP::CmdLine cmd( "CTBronchi pipeline", ' ', "1.0.0" );
+ for( saIt = strArg.begin( ); saIt != strArg.end( ); ++saIt )
+ cmd.add( *( saIt->second ) );
+ for( daIt = dblArg.begin( ); daIt != dblArg.end( ); ++daIt )
+ cmd.add( *( daIt->second ) );
+ cmd.parse( argc, argv );
+ }
+ catch( TCLAP::ArgException& err )
+ {
+ std::stringstream msg;
+ msg << err.error( ) << " " << err.argId( );
+ throw std::runtime_error( msg.str( ) );
+
+ } // yrt
+
+ for( saIt = strArg.begin( ); saIt != strArg.end( ); ++saIt )
+ {
+ std::get< 0 >( this->m_StrArg[ saIt->first ] ) = saIt->second->getValue( );
+ delete saIt->second;
+
+ } // rof
+ for( daIt = dblArg.begin( ); daIt != dblArg.end( ); ++daIt )
+ {
+ std::stringstream vStr;
+ vStr << daIt->second->getValue( );
+ std::get< 0 >( this->m_DblArg[ daIt->first ] ) = vStr.str( );
+ 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( "." ) );
+}
+
+// -------------------------------------------------------------------------
+void CTBronchi::Process::
+Update( )
+{
+ this->_Input( );
+ this->_Vesselness( );
+ this->_Mori( );
+ this->_MoriLabelling( );
+}
+
+// -------------------------------------------------------------------------
+void CTBronchi::Process::
+_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" ] );
+ if( sname == "" && seed == "" )
+ throw std::runtime_error( "No seed given." );
+ if( sname != "" )
+ {
+ std::ifstream str( sname.c_str( ) );
+ str.seekg( 0, std::ios::end );
+ seed.reserve( str.tellg( ) );
+ str.seekg( 0, std::ios::beg );
+ seed.assign(
+ std::istreambuf_iterator< char >( str ),
+ std::istreambuf_iterator< char >( )
+ );
+ str.close( );
+
+ } // fi
+
+ // Get seed
+ std::istringstream sSeed( seed );
+ TPoint pnt;
+ TIndex idx;
+ if( stype == "point" )
+ sSeed >> pnt[ 0 ] >> pnt[ 1 ] >> pnt[ 2 ];
+ else
+ sSeed >> idx[ 0 ] >> idx[ 1 ] >> idx[ 2 ];
+
+ // Read image
+ double t = this->m_Input.Load( iname );
+ if( t < 0 )
+ std::runtime_error( std::string( "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 );
+ else
+ this->m_Input.Get( )->TransformIndexToPhysicalPoint( idx, pnt );
+ this->m_Seed = TSeed( idx, pnt );
+}
+
+// -------------------------------------------------------------------------
+void CTBronchi::Process::
+_Vesselness( )
+{
+ std::string vname = this->m_BaseFileName + "_vesselness.mha";
+ double t = this->m_Vesselness.Load( vname );
+ 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( ) );
+ minMax->Compute( );
+
+ // Invert intensities
+ typedef CTBronchi::Filter< itk::InvertIntensityImageFilter< TPixelImage::TImage > > _TInverter;
+ _TInverter inverter;
+ inverter.Get( )->SetInput( this->m_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;
+ _THessian hessian;
+ hessian.Get( )->SetInput( inverter.Get( )->GetOutput( ) );
+ hessian.Get( )->SetSigma( s );
+ 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( );
+ 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 );
+ t += t1;
+ std::cout << "\"" << vname << "\" saved in " << t1 << " s." << std::endl;
+
+ } // fi
+ std::cout << "Vesselness computed in " << t << " s." << std::endl;
+}
+
+// -------------------------------------------------------------------------
+void CTBronchi::Process::
+_Mori( )
+{
+ std::string mname = this->m_BaseFileName + "_mori.mha";
+ double t = this->m_Mori.Load( mname );
+ 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( );
+ t += t1;
+ std::cout << "Segmentation computed in " << t1 << " s." << std::endl;
+
+ this->m_Mori.Set( mori.Get( )->GetOutput( ) );
+ t1 = this->m_Mori.Save( mname );
+ t += t1;
+ std::cout << "\"" << mname << "\" saved in " << t1 << " s." << std::endl;
+
+ } // fi
+ std::cout << "Mori segmentation computed in " << t << " s." << std::endl;
+}
+
+// -------------------------------------------------------------------------
+void CTBronchi::Process::
+_MoriLabelling( )
+{
+ std::string mname = this->m_BaseFileName + "_labels.mha";
+ double t = this->m_Labels.Load( mname );
+ 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;
+ _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( )->SetInsideValue( this->m_InsideLabel );
+ labelling.Get( )->SetOutsideValue( this->m_OutsideLabel );
+ labelling.Get( )->SetFillValue( this->m_OutsideLabel );
+ double t1 = labelling.Update( );
+ t += t1;
+ std::cout << "Labelling computed in " << t1 << " s." << std::endl;
+
+ this->m_Labels.Set( labelling.Get( )->GetOutput( ) );
+ t1 = this->m_Labels.Save( mname );
+ t += t1;
+ std::cout << "\"" << mname << "\" saved in " << t1 << " s." << std::endl;
+
+ } // fi
+ std::cout << "Mori labelling computed in " << t << " s." << std::endl;
+}
+
+// eof - $RCSfile$
git://git.creatis.insa-lyon.fr / FrontAlgorithms.git / blobdiff