appli/CTBronchi/Process.cxx

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