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 <itkAndImageFilter.h>
  10 #include <itkBinaryThresholdImageFilter.h>
  11 #include <itkMinimumMaximumImageCalculator.h>
  12 #include <itkInvertIntensityImageFilter.h>
  13 #include <itkHessianRecursiveGaussianImageFilter.h>
  14 #include <itkHessian3DToVesselnessMeasureImageFilter.h>
  15
  16 #include <fpa/Common/SliceBySliceRandomWalker.h>
  17 #include <fpa/Filters/Image/ExtractSkeleton.h>
  18 #include <fpa/Filters/Image/Mori.h>
  19 #include <fpa/Filters/Image/RandomWalker.h>
  20 #include <fpa/Functors/Dijkstra/Image/Gaussian.h>
  21
  22 #include "MoriLabelling.h"
  23 #include "Process.h"
  24
  25 // -------------------------------------------------------------------------
  26 CTBronchi::Process::
  27 Process( )
  28   : m_UndefinedLabel( 0 ),
  29     m_InsideLabel( 1 ),
  30     m_OutsideLabel( 2 )
  31 {
  32   this->m_StrArgs[ "input" ] = TStrArg( "", "I", true );
  33   this->m_StrArgs[ "vesselness" ] = TStrArg( "", "V", false );
  34   this->m_StrArgs[ "mori" ] = TStrArg( "", "M", false );
  35   this->m_StrArgs[ "labels" ] = TStrArg( "", "L", false );
  36   this->m_StrArgs[ "fastrw" ] = TStrArg( "", "F", false );
  37   this->m_StrArgs[ "slicerw" ] = TStrArg( "", "R", false );
  38   this->m_StrArgs[ "andrw" ] = TStrArg( "", "A", false );
  39   this->m_StrArgs[ "fastrw_skeleton" ] = TStrArg( "", "G", false );
  40   this->m_StrArgs[ "slicerw_skeleton" ] = TStrArg( "", "U", false );
  41   this->m_StrArgs[ "andrw_skeleton" ] = TStrArg( "", "W", false );
  42   this->m_StrArgs[ "fastrw_points" ] = TStrArg( "", "B", false );
  43   this->m_StrArgs[ "slicerw_points" ] = TStrArg( "", "C", false );
  44   this->m_StrArgs[ "andrw_points" ] = TStrArg( "", "D", false );
  45   this->m_StrArgs[ "seed_file" ] = TStrArg( "", "P", false );
  46   this->m_StrArgs[ "seed" ] = TStrArg( "", "S", false );
  47   this->m_StrArgs[ "seed_type" ] = TStrArg( "point", "T", false );
  48
  49   this->m_DblArgs[ "beta" ] = TDblArg( 2.5, "b", false );
  50   this->m_DblArgs[ "epsilon" ] = TDblArg( 1e-5, "e", false );
  51   this->m_DblArgs[ "vesselness_sigma" ] = TDblArg( 0.5, "a", false );
  52   this->m_DblArgs[ "vesselness_alpha1" ] = TDblArg( 0.5, "c", false );
  53   this->m_DblArgs[ "vesselness_alpha2" ] = TDblArg(  2, "d", false );
  54   this->m_DblArgs[ "mori_min_thr" ] = TDblArg( -850, "f", false );
  55   this->m_DblArgs[ "mori_kernel" ] = TDblArg( 20, "g", false );
  56   this->m_DblArgs[ "mori_signal_thr" ] = TDblArg( 100, "i", false );
  57   this->m_DblArgs[ "mori_signal_influence" ] = TDblArg( 0.5, "j", false );
  58   this->m_DblArgs[ "mori_lower" ] = TDblArg( -1024, "k", false );
  59   this->m_DblArgs[ "mori_upper" ] = TDblArg( 0, "l", false );
  60   this->m_DblArgs[ "mori_delta" ] = TDblArg( 1, "m", false );
  61   this->m_DblArgs[ "slicerw_thr" ] = TDblArg( 5, "n", false );
  62   this->m_DblArgs[ "fastrw_thr" ] = TDblArg( 5, "o", false );
  63   this->m_DblArgs[ "labels_upper_thr" ] = TDblArg( -400, "p", false );
  64 }
  65
  66 // -------------------------------------------------------------------------
  67 CTBronchi::Process::
  68 ~Process( )
  69 {
  70 }
  71
  72 // -------------------------------------------------------------------------
  73 void CTBronchi::Process::
  74 ParseArguments( int argc, char* argv[] )
  75 {
  76   typedef TCLAP::ValueArg< TString > _TStrArg;
  77   typedef TCLAP::ValueArg< TScalar > _TDblArg;
  78   std::map< TString, _TStrArg* > strArgs;
  79   std::map< TString, _TDblArg* > dblArgs;
  80
  81   // Load arguments
  82   TStrArgs::iterator sIt = this->m_StrArgs.begin( );
  83   for( ; sIt != this->m_StrArgs.end( ); ++sIt )
  84     strArgs[ sIt->first ] = new _TStrArg(
  85       std::get< 1 >( sIt->second ),
  86       sIt->first, sIt->first,
  87       std::get< 2 >( sIt->second ),
  88       std::get< 0 >( sIt->second ),
  89       TString( "string (" ) + std::get< 0 >( sIt->second ) + ")"
  90       );
  91
  92   TDblArgs::iterator dIt = this->m_DblArgs.begin( );
  93   for( ; dIt != this->m_DblArgs.end( ); ++dIt )
  94   {
  95     std::stringstream v;
  96     v << "value (" << std::get< 0 >( dIt->second ) << ")";
  97     dblArgs[ dIt->first ] = new _TDblArg(
  98       std::get< 1 >( dIt->second ),
  99       dIt->first, dIt->first,
 100       std::get< 2 >( dIt->second ),
 101       std::get< 0 >( dIt->second ),
 102       v.str( )
 103       );
 104
 105   } // rof
 106
 107   // Parse arguments
 108   std::map< TString, _TStrArg* >::iterator saIt;
 109   std::map< TString, _TDblArg* >::iterator daIt;
 110   try
 111   {
 112     TCLAP::CmdLine cmd( "CTBronchi pipeline", ' ', "1.0.0" );
 113     for( saIt = strArgs.begin( ); saIt != strArgs.end( ); ++saIt )
 114       cmd.add( *( saIt->second ) );
 115     for( daIt = dblArgs.begin( ); daIt != dblArgs.end( ); ++daIt )
 116       cmd.add( *( daIt->second ) );
 117     cmd.parse( argc, argv );
 118   }
 119   catch( TCLAP::ArgException& err )
 120   {
 121     std::stringstream msg;
 122     msg << err.error( ) << " " << err.argId( );
 123     throw std::runtime_error( msg.str( ) );
 124
 125   } // yrt
 126
 127   // Get values
 128   for( saIt = strArgs.begin( ); saIt != strArgs.end( ); ++saIt )
 129   {
 130     std::get< 0 >( this->m_StrArgs[ saIt->first ] ) =
 131       saIt->second->getValue( );
 132     delete saIt->second;
 133
 134   } // rof
 135   for( daIt = dblArgs.begin( ); daIt != dblArgs.end( ); ++daIt )
 136   {
 137     std::get< 0 >( this->m_DblArgs[ daIt->first ] ) =
 138       daIt->second->getValue( );
 139     delete daIt->second;
 140
 141   } // rof
 142
 143   // Compute base filenames
 144   TString bname = std::get< 0 >( this->m_StrArgs[ "input" ] );
 145   bname = bname.substr( 0, bname.find_last_of( "." ) );
 146   const unsigned int N = 6;
 147   const unsigned int M = 6;
 148   TString names[ N + M ] =
 149     {
 150       "vesselness",
 151       "mori",
 152       "labels",
 153       "fastrw",
 154       "slicerw",
 155       "andrw",
 156       "fastrw_skeleton",
 157       "slicerw_skeleton",
 158       "andrw_skeleton",
 159       "fastrw_points",
 160       "slicerw_points",
 161       "andrw_points"
 162     };
 163   for( unsigned int i = 0; i < N; ++i )
 164     if( std::get< 0 >( this->m_StrArgs[ names[ i ] ] ) == "" )
 165       std::get< 0 >( this->m_StrArgs[ names[ i ] ] ) =
 166         bname + "_" + names[ i ] + ".mha";
 167   for( unsigned int i = 0; i < M; ++i )
 168     if( std::get< 0 >( this->m_StrArgs[ names[ i + N ] ] ) == "" )
 169       std::get< 0 >( this->m_StrArgs[ names[ i + N ] ] ) =
 170         bname + "_" + names[ i + N ] + ".txt";
 171 }
 172
 173 // -------------------------------------------------------------------------
 174 void CTBronchi::Process::
 175 Update( )
 176 {
 177   this->_Input( this->m_Input );
 178   this->_Vesselness( this->m_Input, this->m_Vesselness );
 179   this->_Mori( this->m_Input, this->m_Mori, this->m_Seed );
 180   this->_MoriLabelling(
 181     this->m_Input, this->m_Mori, this->m_Vesselness, this->m_Labels
 182     );
 183   this->_FastRW( this->m_Input, this->m_Labels, this->m_FastRW );
 184   this->_SliceRW(
 185     this->m_Input, this->m_Mori, this->m_Vesselness, this->m_SliceRW
 186     );
 187   this->_AndImages( this->m_FastRW, this->m_SliceRW, this->m_AndRW );
 188   this->_Skeleton(
 189     this->m_SliceRW, this->m_SliceRWSkeleton,
 190     std::get< 0 >( this->m_StrArgs[ "slicerw_skeleton" ] ),
 191     std::get< 0 >( this->m_StrArgs[ "slicerw_points" ] )
 192     );
 193   this->_Skeleton(
 194     this->m_AndRW, this->m_AndRWSkeleton,
 195     std::get< 0 >( this->m_StrArgs[ "andrw_skeleton" ] ),
 196     std::get< 0 >( this->m_StrArgs[ "andrw_points" ] )
 197     );
 198   this->_Skeleton(
 199     this->m_FastRW, this->m_FastRWSkeleton,
 200     std::get< 0 >( this->m_StrArgs[ "fastrw_skeleton" ] ),
 201     std::get< 0 >( this->m_StrArgs[ "fastrw_points" ] )
 202     );
 203 }
 204
 205 // -------------------------------------------------------------------------
 206 template< class _TImage >
 207 void CTBronchi::Process::
 208 _Input( _TImage& input )
 209 {
 210   TString sname = std::get< 0 >( this->m_StrArgs[ "seed_file" ] );
 211   TString seed = std::get< 0 >( this->m_StrArgs[ "seed" ] );
 212   TString stype = std::get< 0 >( this->m_StrArgs[ "seed_type" ] );
 213   if( sname == "" && seed == "" )
 214     throw std::runtime_error( "No seed given." );
 215   if( sname != "" )
 216   {
 217     std::ifstream str( sname.c_str( ) );
 218     str.seekg( 0, std::ios::end );
 219     seed.reserve( str.tellg( ) );
 220     str.seekg( 0, std::ios::beg );
 221     seed.assign(
 222       std::istreambuf_iterator< char >( str ),
 223       std::istreambuf_iterator< char >( )
 224       );
 225     str.close( );
 226
 227   } // fi
 228
 229   // Get seed
 230   std::istringstream sSeed( seed );
 231   TPoint pnt;
 232   TIndex idx;
 233   if( stype == "point" )
 234     sSeed >> pnt[ 0 ] >> pnt[ 1 ] >> pnt[ 2 ];
 235   else
 236     sSeed >> idx[ 0 ] >> idx[ 1 ] >> idx[ 2 ];
 237
 238   // Read image
 239   TString iname = std::get< 0 >( this->m_StrArgs[ "input" ] );
 240   double t = input.Load( iname );
 241   if( t < 0 )
 242     std::runtime_error( TString( "Could not read \"" ) + iname + "\"" );
 243   std::cout << "\"" << iname << "\" read in " << t << " s." << std::endl;
 244
 245   // Synch seed
 246   if( stype == "point" )
 247     input.Get( )->TransformPhysicalPointToIndex( pnt, idx );
 248   else
 249     input.Get( )->TransformIndexToPhysicalPoint( idx, pnt );
 250   this->m_Seed = TSeed( idx, pnt );
 251 }
 252
 253 // -------------------------------------------------------------------------
 254 template< class _TInput, class _TVesselness >
 255 void CTBronchi::Process::
 256 _Vesselness( _TInput& input, _TVesselness& vesselness )
 257 {
 258   TString iname = std::get< 0 >( this->m_StrArgs[ "vesselness" ] );
 259   double t = vesselness.Load( iname );
 260   if( t < 0 )
 261   {
 262     t = 0;
 263
 264     // Min-max
 265     typedef itk::MinimumMaximumImageCalculator< typename _TInput::TImage > _TMinMax;
 266     typename _TMinMax::Pointer minMax = _TMinMax::New( );
 267     minMax->SetImage( input.Get( ) );
 268     minMax->Compute( );
 269
 270     // Invert intensities
 271     typedef CTBronchi::Filter< itk::InvertIntensityImageFilter< typename _TInput::TImage > > _TInverter;
 272     _TInverter inverter;
 273     inverter.Get( )->SetInput( input.Get( ) );
 274     inverter.Get( )->SetMaximum( minMax->GetMaximum( ) );
 275     double t1 = inverter.Update( );
 276     t += t1;
 277     std::cout << "Inversion executed in " << t1 << " s" << std::endl;
 278
 279     // Compute hessian image
 280     typedef CTBronchi::Filter< itk::HessianRecursiveGaussianImageFilter< typename _TInput::TImage > > _THessian;
 281     _THessian hessian;
 282     hessian.Get( )->SetInput( inverter.Get( )->GetOutput( ) );
 283     hessian.Get( )->SetSigma(
 284       std::get< 0 >( this->m_DblArgs[ "vesselness_sigma" ] )
 285       );
 286     t1 = hessian.Update( );
 287     t += t1;
 288     std::cout << "Hessian executed in " << t1 << " s" << std::endl;
 289
 290     // Vesselness
 291     typedef CTBronchi::Filter< itk::Hessian3DToVesselnessMeasureImageFilter< typename _TVesselness::TPixel > > _TVesselnessFilter;
 292     _TVesselnessFilter vFilter;
 293     vFilter.Get( )->SetInput( hessian.Get( )->GetOutput( ) );
 294     vFilter.Get( )->SetAlpha1(
 295       std::get< 0 >( this->m_DblArgs[ "vesselness_alpha1" ] )
 296       );
 297     vFilter.Get( )->SetAlpha2(
 298       std::get< 0 >( this->m_DblArgs[ "vesselness_alpha2" ] )
 299       );
 300     t1 = vFilter.Update( );
 301     t += t1;
 302     std::cout << "Hessian measure computed in " << t1 << " s." << std::endl;
 303
 304     vesselness.Set( vFilter.Get( )->GetOutput( ) );
 305     t1 = vesselness.Save( iname );
 306     t += t1;
 307     std::cout << "\"" << iname << "\" saved in " << t1 << " s." << std::endl;
 308
 309   } // fi
 310   std::cout << "Vesselness computed in " << t << " s." << std::endl;
 311 }
 312
 313 // -------------------------------------------------------------------------
 314 template< class _TInput, class _TMori, class _TSeed >
 315 void CTBronchi::Process::
 316 _Mori( _TInput& input, _TMori& mori, _TSeed& seed )
 317 {
 318   TString iname = std::get< 0 >( this->m_StrArgs[ "mori" ] );
 319   double t = mori.Load( iname );
 320   if( t < 0 )
 321   {
 322     t = 0;
 323
 324     // Mori segmentation
 325     typedef CTBronchi::Filter< fpa::Filters::Image::Mori< TPixelImage::TImage, TLabelImage::TImage > > _TMoriFilter;
 326     _TMoriFilter mFilter;
 327     mFilter.Get( )->SetInput( input.Get( ) );
 328     mFilter.Get( )->SetSeed( seed.second );
 329     mFilter.Get( )->SetInsideValue( this->m_InsideLabel );
 330     mFilter.Get( )->SetOutsideValue( this->m_UndefinedLabel );
 331     mFilter.Get( )->SetMinimumThreshold(
 332       std::get< 0 >( this->m_DblArgs[ "mori_min_thr" ] )
 333       );
 334     mFilter.Get( )->SetSignalKernelSize(
 335       std::get< 0 >( this->m_DblArgs[ "mori_kernel" ] )
 336       );
 337     mFilter.Get( )->SetSignalThreshold(
 338       std::get< 0 >( this->m_DblArgs[ "mori_signal_thr" ] )
 339       );
 340     mFilter.Get( )->SetSignalInfluence(
 341       std::get< 0 >( this->m_DblArgs[ "mori_signal_influence" ] )
 342       );
 343     mFilter.Get( )->SetThresholds(
 344       std::get< 0 >( this->m_DblArgs[ "mori_lower" ] ),
 345       std::get< 0 >( this->m_DblArgs[ "mori_upper" ] ),
 346       std::get< 0 >( this->m_DblArgs[ "mori_delta" ] )
 347       );
 348     double t1 = mFilter.Update( );
 349     t += t1;
 350     std::cout << "Segmentation computed in " << t1 << " s." << std::endl;
 351
 352     mori.Set( mFilter.Get( )->GetOutput( ) );
 353     t1 = mori.Save( iname );
 354     t += t1;
 355     std::cout << "\"" << iname << "\" saved in " << t1 << " s." << std::endl;
 356
 357   } // fi
 358   std::cout << "Mori segmentation computed in " << t << " s." << std::endl;
 359 }
 360
 361 // -------------------------------------------------------------------------
 362 template< class _TInput, class _TMori, class _TVesselness, class _TLabels >
 363 void CTBronchi::Process::
 364 _MoriLabelling(
 365   _TInput& input, _TMori& mori, _TVesselness& vesselness, _TLabels& labels
 366   )
 367 {
 368   TString iname = std::get< 0 >( this->m_StrArgs[ "labels" ] );
 369   double t = labels.Load( iname );
 370   if( t < 0 )
 371   {
 372     t = 0;
 373
 374     // Mori labelling
 375     typedef CTBronchi::Filter< CTBronchi::MoriLabelling< typename _TInput::TImage, typename _TLabels::TImage, typename _TVesselness::TImage > > _TLabelling;
 376     _TLabelling labelling;
 377     labelling.Get( )->SetInput( input.Get( ) );
 378     labelling.Get( )->SetInputLabels( mori.Get( ) );
 379     labelling.Get( )->SetInputVesselness( vesselness.Get( ) );
 380     labelling.Get( )->SetVesselnessThreshold( std::get< 0 >( this->m_DblArgs[ "fastrw_thr" ] ) );
 381     labelling.Get( )->SetUpperThreshold( std::get< 0 >( this->m_DblArgs[ "labels_upper_thr" ] ) );
 382     labelling.Get( )->SetInsideValue( this->m_InsideLabel );
 383     labelling.Get( )->SetOutsideValue( this->m_OutsideLabel );
 384     labelling.Get( )->SetFillValue( this->m_OutsideLabel );
 385     double t1 = labelling.Update( );
 386     t += t1;
 387     std::cout << "Labelling computed in " << t1 << " s." << std::endl;
 388
 389     labels.Set( labelling.Get( )->GetOutput( ) );
 390     t1 = labels.Save( iname );
 391     t += t1;
 392     std::cout << "\"" << iname << "\" saved in " << t1 << " s." << std::endl;
 393
 394   } // fi
 395   std::cout << "Mori labelling computed in " << t << " s." << std::endl;
 396 }
 397
 398 // -------------------------------------------------------------------------
 399 template< class _TInput, class _TLabels, class _TFastRW >
 400 void CTBronchi::Process::
 401 _FastRW( _TInput& input, _TLabels& labels, _TFastRW& fastrw )
 402 {
 403   TString iname = std::get< 0 >( this->m_StrArgs[ "fastrw" ] );
 404   double t = fastrw.Load( iname );
 405   if( t < 0 )
 406   {
 407     t = 0;
 408
 409     // Prepare weight functor
 410     typedef fpa::Functors::Dijkstra::Image::Gaussian< typename _TInput::TImage, typename _TFastRW::TPixel > TWeight;
 411     typename TWeight::Pointer weight = TWeight::New( );
 412     weight->SetBeta( std::get< 0 >( this->m_DblArgs[ "beta" ] ) );
 413     weight->SetEpsilon( std::get< 0 >( this->m_DblArgs[ "epsilon" ] ) );
 414
 415     // Random walk
 416     typedef CTBronchi::Filter< fpa::Filters::Image::RandomWalker< typename _TInput::TImage, typename _TLabels::TImage, typename _TFastRW::TPixel > > TFilter;
 417     TFilter filter;
 418     filter.Get( )->SetInputImage( input.Get( ) );
 419     filter.Get( )->SetInputLabels( labels.Get( ) );
 420     filter.Get( )->SetWeightFunction( weight );
 421     double t1 = filter.Update( );
 422     t += t1;
 423     std::cout << "Fast random walker executed in " << t1 << " s" << std::endl;
 424
 425     // Extract label
 426     typedef CTBronchi::Filter< itk::BinaryThresholdImageFilter< typename _TLabels::TImage, typename _TLabels::TImage > > _TExtract;
 427     _TExtract extract;
 428     extract.Get( )->SetInput( filter.Get( )->GetOutputLabels( ) );
 429     extract.Get( )->SetInsideValue( this->m_InsideLabel );
 430     extract.Get( )->SetOutsideValue( this->m_UndefinedLabel );
 431     extract.Get( )->SetLowerThreshold( this->m_InsideLabel );
 432     extract.Get( )->SetUpperThreshold( this->m_InsideLabel );
 433     t1 = extract.Update( );
 434     t += t1;
 435     std::cout << "Extract labels executed in " << t1 << " s" << std::endl;
 436
 437     fastrw.Set( extract.Get( )->GetOutput( ) );
 438     t1 = fastrw.Save( iname );
 439     t += t1;
 440     std::cout << "\"" << iname << "\" saved in " << t1 << " s." << std::endl;
 441
 442   } // fi
 443   std::cout << "Fast random walker computed in " << t << " s." << std::endl;
 444 }
 445
 446 // -------------------------------------------------------------------------
 447 template< class _TInput, class _TLabels, class _TVesselness, class _TSliceRW >
 448 void CTBronchi::Process::
 449 _SliceRW(
 450   _TInput& input, _TLabels& labels, _TVesselness& vesselness, _TSliceRW& slicerw
 451   )
 452 {
 453   TString iname = std::get< 0 >( this->m_StrArgs[ "slicerw" ] );
 454   double t = slicerw.Load( iname );
 455   if( t < 0 )
 456   {
 457     t = 0;
 458
 459     // Random walk
 460     typedef CTBronchi::Filter< fpa::Common::SliceBySliceRandomWalker< typename _TInput::TImage, typename _TLabels::TImage, typename _TVesselness::TImage > > TFilter;
 461     TFilter filter;
 462     filter.Get( )->SetInput( input.Get( ) );
 463     filter.Get( )->SetInputLabels( labels.Get( ) );
 464     filter.Get( )->SetInputVesselness( vesselness.Get( ) );
 465     filter.Get( )->SetBeta( std::get< 0 >( this->m_DblArgs[ "beta" ] ) );
 466     filter.Get( )->SetVesselnessThreshold( std::get< 0 >( this->m_DblArgs[ "slicerw_thr" ] ) );
 467     filter.Get( )->SetEpsilon( std::get< 0 >( this->m_DblArgs[ "epsilon" ] ) );
 468     double t1 = filter.Update( );
 469     t += t1;
 470     std::cout << "Extract labels executed in " << t1 << " s" << std::endl;
 471
 472     slicerw.Set( filter.Get( )->GetOutput( ) );
 473     t1 = slicerw.Save( iname );
 474     t += t1;
 475     std::cout << "\"" << iname << "\" saved in " << t1 << " s." << std::endl;
 476
 477   } // fi
 478   std::cout << "Slice by slice random walker computed in " << t << " s." << std::endl;
 479 }
 480
 481 // -------------------------------------------------------------------------
 482 template< class _TInput >
 483 void CTBronchi::Process::
 484 _AndImages( _TInput& a, _TInput& b, _TInput& c )
 485 {
 486   TString iname = std::get< 0 >( this->m_StrArgs[ "andrw" ] );
 487   double t = c.Load( iname );
 488   if( t < 0 )
 489   {
 490     t = 0;
 491
 492     // Random walk
 493     typedef CTBronchi::Filter< itk::AndImageFilter< typename _TInput::TImage > > TFilter;
 494     TFilter filter;
 495     filter.Get( )->SetInput( 0, a.Get( ) );
 496     filter.Get( )->SetInput( 1, b.Get( ) );
 497     double t1 = filter.Update( );
 498     t += t1;
 499     std::cout << "And segmentations executed in " << t1 << " s" << std::endl;
 500
 501     c.Set( filter.Get( )->GetOutput( ) );
 502     t1 = c.Save( iname );
 503     t += t1;
 504     std::cout << "\"" << iname << "\" saved in " << t1 << " s." << std::endl;
 505
 506   } // fi
 507   std::cout << "And segmentations computed in " << t << " s." << std::endl;
 508 }
 509
 510 // -------------------------------------------------------------------------
 511 template< class _TInput, class _TSkeleton >
 512 void CTBronchi::Process::
 513 _Skeleton(
 514   _TInput& a, _TSkeleton& s,
 515   const TString& fname,
 516   const TString& pname
 517   )
 518 {
 519   double t = s.Load( fname );
 520   if( t < 0 )
 521   {
 522     t = 0;
 523
 524     typedef CTBronchi::Filter< fpa::Filters::Image::ExtractSkeleton< typename _TInput::TImage > > TFilter;
 525     TFilter filter;
 526     filter.Get( )->SetInput( a.Get( ) );
 527     filter.Get( )->SeedFromMaximumDistanceOff( );
 528     filter.Get( )->SetSeed( this->m_Seed.first );
 529     filter.Get( )->GetDistanceMap( )->InsideIsPositiveOn( );
 530     filter.Get( )->GetDistanceMap( )->SquaredDistanceOff( );
 531     filter.Get( )->GetDistanceMap( )->UseImageSpacingOn( );
 532     double t1 = filter.Update( );
 533     t += t1;
 534     std::cout << "Skeleton executed in " << t1 << " s" << std::endl;
 535
 536     s.Set( filter.Get( )->GetOutput( ) );
 537     t1 = s.Save( fname );
 538     t += t1;
 539     std::cout << "\"" << fname << "\" saved in " << t1 << " s." << std::endl;
 540
 541     // End points
 542     std::stringstream ePointsStr;
 543     std::vector< typename _TInput::TImage::IndexType > ePoints =
 544       filter.Get( )->GetEndPoints( );
 545     for( typename _TInput::TImage::IndexType i: ePoints )
 546     {
 547       for( unsigned int d = 0; d < _TInput::VDim; ++d )
 548         ePointsStr << i[ d ] << " ";
 549       ePointsStr << std::endl;
 550
 551     } // rof
 552
 553     std::ofstream ePointsF( pname.c_str( ), std::ofstream::binary );
 554     if( !ePointsF )
 555       throw std::runtime_error(
 556         TString( "Unable to write skeleton to \"" ) + pname + "\""
 557         );
 558     ePointsF.write( ePointsStr.str( ).c_str( ), ePointsStr.str( ).size( ) );
 559
 560   } // fi
 561   std::cout << "Skeleton computed in " << t << " s." << std::endl;
 562 }
 563
 564 // eof - $RCSfile$