appli/CTBronchi/Process.cxx

   1 // =========================================================================
   2 // @author Leonardo Florez Valencia (florez-l@javeriana.edu.co)
   3 // =========================================================================
   4
   5 #include <fstream>
   6 #include <streambuf>
   7 #include <tclap/CmdLine.h>
   8
   9 #include <itkMinimumMaximumImageCalculator.h>
  10 #include <itkInvertIntensityImageFilter.h>
  11 #include <itkHessianRecursiveGaussianImageFilter.h>
  12 #include <itkHessian3DToVesselnessMeasureImageFilter.h>
  13
  14 #include <fpa/Filters/Image/Mori.h>
  15
  16 #include "MoriLabelling.h"
  17 #include "Process.h"
  18
  19 // -------------------------------------------------------------------------
  20 CTBronchi::Process::
  21 Process( )
  22   : m_UndefinedLabel( 0 ),
  23     m_InsideLabel( 1 ),
  24     m_OutsideLabel( 2 )
  25 {
  26   this->m_StrArg[ "input" ] = TArgument( "", true );
  27   this->m_StrArg[ "seed_file" ] = TArgument( "", false );
  28   this->m_StrArg[ "seed" ] = TArgument( "", false );
  29   this->m_StrArg[ "seed_type" ] = TArgument( "point", false );
  30
  31   this->m_DblArg[ "beta" ] = TArgument( "2.5", false );
  32   this->m_DblArg[ "epsilon" ] = TArgument( "1e-5", false );
  33   this->m_DblArg[ "vesselness_sigma" ] = TArgument( "0.5", false );
  34   this->m_DblArg[ "vesselness_alpha1" ] = TArgument( "0.5", false );
  35   this->m_DblArg[ "vesselness_alpha2" ] = TArgument(  "2", false );
  36   this->m_DblArg[ "mori_min_thr" ] = TArgument( "-850", false );
  37   this->m_DblArg[ "mori_kernel" ] = TArgument( "20", false );
  38   this->m_DblArg[ "mori_signal_thr" ] = TArgument( "100", false );
  39   this->m_DblArg[ "mori_signal_influence" ] = TArgument( "0.5", false );
  40   this->m_DblArg[ "mori_lower" ] = TArgument( "-1024", false );
  41   this->m_DblArg[ "mori_upper" ] = TArgument( "0", false );
  42   this->m_DblArg[ "mori_delta" ] = TArgument( "1", false );
  43   this->m_DblArg[ "slicerw_thr" ] = TArgument( "5", false );
  44   this->m_DblArg[ "fastrw_thr" ] = TArgument( "65", false );
  45   this->m_DblArg[ "labels_upper_thr" ] = TArgument( "-400", false );
  46 }
  47
  48 // -------------------------------------------------------------------------
  49 CTBronchi::Process::
  50 ~Process( )
  51 {
  52 }
  53
  54 // -------------------------------------------------------------------------
  55 void CTBronchi::Process::
  56 ParseArguments( int argc, char* argv[] )
  57 {
  58   typedef TCLAP::ValueArg< std::string > _TStrArg;
  59   typedef TCLAP::ValueArg< TScalar >     _TDblArg;
  60   std::map< std::string, _TStrArg* > strArg;
  61   std::map< std::string, _TDblArg* > dblArg;
  62   unsigned char flag = 'a';
  63
  64   TArguments::iterator sIt = this->m_StrArg.begin( );
  65   for( ; sIt != this->m_StrArg.end( ); ++sIt )
  66   {
  67     std::stringstream fStr;
  68     fStr << flag;
  69     flag++;
  70     if( flag == 'h' )
  71       flag++;
  72
  73     _TStrArg* a = new _TStrArg(
  74       fStr.str( ), sIt->first, sIt->first,
  75       std::get< 1 >( sIt->second ),
  76       std::get< 0 >( sIt->second ),
  77       "string"
  78       );
  79     strArg[ sIt->first ] = a;
  80
  81   } // rof
  82
  83   TArguments::iterator dIt = this->m_DblArg.begin( );
  84   for( ; dIt != this->m_DblArg.end( ); ++dIt )
  85   {
  86     std::stringstream fStr;
  87     fStr << flag;
  88     flag++;
  89     if( flag == 'h' )
  90       flag++;
  91
  92     std::istringstream vStr( std::get< 0 >( dIt->second ) );
  93     TScalar v;
  94     vStr >> v;
  95
  96     _TDblArg* a = new _TDblArg(
  97       fStr.str( ), dIt->first, dIt->first,
  98       std::get< 1 >( dIt->second ),
  99       v, "value"
 100       );
 101     dblArg[ dIt->first ] = a;
 102
 103   } // rof
 104
 105   std::map< std::string, _TStrArg* >::iterator saIt;
 106   std::map< std::string, _TDblArg* >::iterator daIt;
 107   try
 108   {
 109     TCLAP::CmdLine cmd( "CTBronchi pipeline", ' ', "1.0.0" );
 110     for( saIt = strArg.begin( ); saIt != strArg.end( ); ++saIt )
 111       cmd.add( *( saIt->second ) );
 112     for( daIt = dblArg.begin( ); daIt != dblArg.end( ); ++daIt )
 113       cmd.add( *( daIt->second ) );
 114     cmd.parse( argc, argv );
 115   }
 116   catch( TCLAP::ArgException& err )
 117   {
 118     std::stringstream msg;
 119     msg << err.error( ) << " " << err.argId( );
 120     throw std::runtime_error( msg.str( ) );
 121
 122   } // yrt
 123
 124   for( saIt = strArg.begin( ); saIt != strArg.end( ); ++saIt )
 125   {
 126     std::get< 0 >( this->m_StrArg[ saIt->first ] ) = saIt->second->getValue( );
 127     delete saIt->second;
 128
 129   } // rof
 130   for( daIt = dblArg.begin( ); daIt != dblArg.end( ); ++daIt )
 131   {
 132     std::stringstream vStr;
 133     vStr << daIt->second->getValue( );
 134     std::get< 0 >( this->m_DblArg[ daIt->first ] ) = vStr.str( );
 135     delete daIt->second;
 136
 137   } // rof
 138
 139   // Compute base filename
 140   std::string iname = std::get< 0 >( this->m_StrArg[ "input" ] );
 141   this->m_BaseFileName = iname.substr( 0, iname.find_last_of( "." ) );
 142 }
 143
 144 // -------------------------------------------------------------------------
 145 void CTBronchi::Process::
 146 Update( )
 147 {
 148   this->_Input( );
 149   this->_Vesselness( );
 150   this->_Mori( );
 151   this->_MoriLabelling( );
 152 }
 153
 154 // -------------------------------------------------------------------------
 155 void CTBronchi::Process::
 156 _Input( )
 157 {
 158   std::string iname = std::get< 0 >( this->m_StrArg[ "input" ] );
 159   std::string sname = std::get< 0 >( this->m_StrArg[ "seed_file" ] );
 160   std::string seed = std::get< 0 >( this->m_StrArg[ "seed" ] );
 161   std::string stype = std::get< 0 >( this->m_StrArg[ "seed_type" ] );
 162   if( sname == "" && seed == "" )
 163     throw std::runtime_error( "No seed given." );
 164   if( sname != "" )
 165   {
 166     std::ifstream str( sname.c_str( ) );
 167     str.seekg( 0, std::ios::end );
 168     seed.reserve( str.tellg( ) );
 169     str.seekg( 0, std::ios::beg );
 170     seed.assign(
 171       std::istreambuf_iterator< char >( str ),
 172       std::istreambuf_iterator< char >( )
 173       );
 174     str.close( );
 175
 176   } // fi
 177
 178   // Get seed
 179   std::istringstream sSeed( seed );
 180   TPoint pnt;
 181   TIndex idx;
 182   if( stype == "point" )
 183     sSeed >> pnt[ 0 ] >> pnt[ 1 ] >> pnt[ 2 ];
 184   else
 185     sSeed >> idx[ 0 ] >> idx[ 1 ] >> idx[ 2 ];
 186
 187   // Read image
 188   double t = this->m_Input.Load( iname );
 189   if( t < 0 )
 190     std::runtime_error( std::string( "Could not read \"" ) + iname + "\"" );
 191   std::cout << "\"" << iname << "\" read in " << t << " s." << std::endl;
 192
 193   // Synch seed
 194   if( stype == "point" )
 195     this->m_Input.Get( )->TransformPhysicalPointToIndex( pnt, idx );
 196   else
 197     this->m_Input.Get( )->TransformIndexToPhysicalPoint( idx, pnt );
 198   this->m_Seed = TSeed( idx, pnt );
 199 }
 200
 201 // -------------------------------------------------------------------------
 202 void CTBronchi::Process::
 203 _Vesselness( )
 204 {
 205   std::string vname = this->m_BaseFileName + "_vesselness.mha";
 206   double t = this->m_Vesselness.Load( vname );
 207   if( t < 0 )
 208   {
 209     t = 0;
 210
 211     // Arguments
 212     std::stringstream v;
 213     v << std::get< 0 >( this->m_DblArg[ "vesselness_sigma" ] ) << " "
 214       << std::get< 0 >( this->m_DblArg[ "vesselness_alpha1" ] ) << " "
 215       << std::get< 0 >( this->m_DblArg[ "vesselness_alpha2" ] );
 216     std::istringstream w( v.str( ) );
 217     TScalar s, a1, a2;
 218     w >> s >> a1 >> a2;
 219
 220     // Min-max
 221     typedef itk::MinimumMaximumImageCalculator< TPixelImage::TImage > _TMinMax;
 222     _TMinMax::Pointer minMax = _TMinMax::New( );
 223     minMax->SetImage( this->m_Input.Get( ) );
 224     minMax->Compute( );
 225
 226     // Invert intensities
 227     typedef CTBronchi::Filter< itk::InvertIntensityImageFilter< TPixelImage::TImage > > _TInverter;
 228     _TInverter inverter;
 229     inverter.Get( )->SetInput( this->m_Input.Get( ) );
 230     inverter.Get( )->SetMaximum( minMax->GetMaximum( ) );
 231     double t1 = inverter.Update( );
 232     t += t1;
 233     std::cout << "Inversion executed in " << t1 << " s" << std::endl;
 234
 235     // Compute hessian image
 236     typedef CTBronchi::Filter< itk::HessianRecursiveGaussianImageFilter< TPixelImage::TImage > > _THessian;
 237     _THessian hessian;
 238     hessian.Get( )->SetInput( inverter.Get( )->GetOutput( ) );
 239     hessian.Get( )->SetSigma( s );
 240     t1 = hessian.Update( );
 241     t += t1;
 242     std::cout << "Hessian executed in " << t1 << " s" << std::endl;
 243
 244     // Vesselness
 245     typedef CTBronchi::Filter< itk::Hessian3DToVesselnessMeasureImageFilter< TScalar > > _TVesselness;
 246     _TVesselness vesselness;
 247     vesselness.Get( )->SetInput( hessian.Get( )->GetOutput( ) );
 248     vesselness.Get( )->SetAlpha1( a1 );
 249     vesselness.Get( )->SetAlpha2( a2 );
 250     t1 = vesselness.Update( );
 251     t += t1;
 252     std::cout << "Hessian measure computed in " << t1 << " s." << std::endl;
 253
 254     this->m_Vesselness.Set( vesselness.Get( )->GetOutput( ) );
 255     t1 = this->m_Vesselness.Save( vname );
 256     t += t1;
 257     std::cout << "\"" << vname << "\" saved in " << t1 << " s." << std::endl;
 258
 259   } // fi
 260   std::cout << "Vesselness computed in " << t << " s." << std::endl;
 261 }
 262
 263 // -------------------------------------------------------------------------
 264 void CTBronchi::Process::
 265 _Mori( )
 266 {
 267   std::string mname = this->m_BaseFileName + "_mori.mha";
 268   double t = this->m_Mori.Load( mname );
 269   if( t < 0 )
 270   {
 271     t = 0;
 272
 273     // Arguments
 274     std::stringstream v;
 275     v << std::get< 0 >( this->m_DblArg[ "mori_min_thr" ] ) << " "
 276       << std::get< 0 >( this->m_DblArg[ "mori_kernel" ] ) << " "
 277       << std::get< 0 >( this->m_DblArg[ "mori_signal_thr" ] ) << " "
 278       << std::get< 0 >( this->m_DblArg[ "mori_signal_influence" ] ) << " "
 279       << std::get< 0 >( this->m_DblArg[ "mori_lower" ] ) << " "
 280       << std::get< 0 >( this->m_DblArg[ "mori_upper" ] ) << " "
 281       << std::get< 0 >( this->m_DblArg[ "mori_delta" ] );
 282     std::istringstream w( v.str( ) );
 283     TScalar mThr, sKernel, sThr, sInfluence, lower, upper, delta;
 284     w >> mThr >> sKernel >> sThr >> sInfluence >> lower >> upper >> delta;
 285
 286     // Mori segmentation
 287     typedef CTBronchi::Filter< fpa::Filters::Image::Mori< TPixelImage::TImage, TLabelImage::TImage > > _TMori;
 288     _TMori mori;
 289     mori.Get( )->SetInput( this->m_Input.Get( ) );
 290     mori.Get( )->SetSeed( this->m_Seed.second );
 291     mori.Get( )->SetInsideValue( this->m_InsideLabel );
 292     mori.Get( )->SetOutsideValue( this->m_UndefinedLabel );
 293     mori.Get( )->SetMinimumThreshold( mThr );
 294     mori.Get( )->SetSignalKernelSize( sKernel );
 295     mori.Get( )->SetSignalThreshold( sThr );
 296     mori.Get( )->SetSignalInfluence( sInfluence );
 297     mori.Get( )->SetThresholds( lower, upper, delta );
 298     double t1 = mori.Update( );
 299     t += t1;
 300     std::cout << "Segmentation computed in " << t1 << " s." << std::endl;
 301
 302     this->m_Mori.Set( mori.Get( )->GetOutput( ) );
 303     t1 = this->m_Mori.Save( mname );
 304     t += t1;
 305     std::cout << "\"" << mname << "\" saved in " << t1 << " s." << std::endl;
 306
 307   } // fi
 308   std::cout << "Mori segmentation computed in " << t << " s." << std::endl;
 309 }
 310
 311 // -------------------------------------------------------------------------
 312 void CTBronchi::Process::
 313 _MoriLabelling( )
 314 {
 315   std::string mname = this->m_BaseFileName + "_labels.mha";
 316   double t = this->m_Labels.Load( mname );
 317   if( t < 0 )
 318   {
 319     t = 0;
 320
 321     // Arguments
 322     std::stringstream v;
 323     v << std::get< 0 >( this->m_DblArg[ "fastrw_thr" ] ) << " "
 324       << std::get< 0 >( this->m_DblArg[ "labels_upper_thr" ] );
 325     std::istringstream w( v.str( ) );
 326     TScalar vThr, uThr;
 327     w >> vThr >> uThr;
 328
 329     // Mori labelling
 330     typedef CTBronchi::Filter< CTBronchi::MoriLabelling< TPixelImage::TImage, TLabelImage::TImage, TScalarImage::TImage > > _TLabelling;
 331     _TLabelling labelling;
 332     labelling.Get( )->SetInput( this->m_Input.Get( ) );
 333     labelling.Get( )->SetInputLabels( this->m_Mori.Get( ) );
 334     labelling.Get( )->SetInputVesselness( this->m_Vesselness.Get( ) );
 335     labelling.Get( )->SetVesselnessThreshold( vThr );
 336     labelling.Get( )->SetUpperThreshold( uThr );
 337     labelling.Get( )->SetInsideValue( this->m_InsideLabel );
 338     labelling.Get( )->SetOutsideValue( this->m_OutsideLabel );
 339     labelling.Get( )->SetFillValue( this->m_OutsideLabel );
 340     double t1 = labelling.Update( );
 341     t += t1;
 342     std::cout << "Labelling computed in " << t1 << " s." << std::endl;
 343
 344     this->m_Labels.Set( labelling.Get( )->GetOutput( ) );
 345     t1 = this->m_Labels.Save( mname );
 346     t += t1;
 347     std::cout << "\"" << mname << "\" saved in " << t1 << " s." << std::endl;
 348
 349   } // fi
 350   std::cout << "Mori labelling computed in " << t << " s." << std::endl;
 351 }
 352
 353 // eof - $RCSfile$