segmentation/clitkMotionMaskGenericFilter.txx

   1 /*=========================================================================
   2   Program:   vv                     http://www.creatis.insa-lyon.fr/rio/vv
   3
   4   Authors belong to:
   5   - University of LYON              http://www.universite-lyon.fr/
   6   - Léon Bérard cancer center       http://oncora1.lyon.fnclcc.fr
   7   - CREATIS CNRS laboratory         http://www.creatis.insa-lyon.fr
   8
   9   This software is distributed WITHOUT ANY WARRANTY; without even
  10   the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  11   PURPOSE.  See the copyright notices for more information.
  12
  13   It is distributed under dual licence
  14
  15   - BSD        See included LICENSE.txt file
  16   - CeCILL-B   http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
  17 ======================================================================-====*/
  18 #ifndef clitkMotionMaskGenericFilter_txx
  19 #define clitkMotionMaskGenericFilter_txx
  20
  21 /* =================================================
  22  * @file   clitkMotionMaskGenericFilter.txx
  23  * @author
  24  * @date
  25  *
  26  * @brief
  27  *
  28  ===================================================*/
  29
  30
  31 namespace clitk
  32 {
  33
  34 //-------------------------------------------------------------------
  35 // Update with the number of dimensions
  36 //-------------------------------------------------------------------
  37 template<unsigned int Dimension>
  38 void
  39 MotionMaskGenericFilter::UpdateWithDim(std::string PixelType)
  40 {
  41   if (m_Verbose) std::cout << "Image was detected to be "<<Dimension<<"D and "<< PixelType<<"..."<<std::endl;
  42
  43   if(PixelType == "short") {
  44     if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and signed short..." << std::endl;
  45     UpdateWithDimAndPixelType<Dimension, signed short>();
  46   }
  47   //    else if(PixelType == "unsigned_short"){
  48   //       if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and unsigned_short..." << std::endl;
  49   //       UpdateWithDimAndPixelType<Dimension, unsigned short>();
  50   //     }
  51
  52   //     else if (PixelType == "unsigned_char"){
  53   //       if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and unsigned_char..." << std::endl;
  54   //       UpdateWithDimAndPixelType<Dimension, unsigned char>();
  55   //     }
  56
  57   //     else if (PixelType == "char"){
  58   //       if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and signed_char..." << std::endl;
  59   //       UpdateWithDimAndPixelType<Dimension, signed char>();
  60   //     }
  61   else {
  62     if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and float..." << std::endl;
  63     UpdateWithDimAndPixelType<Dimension, float>();
  64   }
  65 }
  66
  67 //-------------------------------------------------------------------
  68 // Air
  69 //-------------------------------------------------------------------
  70 template <unsigned int Dimension, class  PixelType>
  71 typename itk::Image<MotionMaskGenericFilter::InternalPixelType, Dimension>::Pointer
  72 MotionMaskGenericFilter::GetAirImage(typename itk::Image<PixelType, Dimension>::Pointer input,
  73                                      typename itk::Image<MotionMaskGenericFilter::InternalPixelType, Dimension>::Pointer lungs)
  74 {
  75
  76   // ImageTypes
  77   typedef int InternalPixelType;
  78   typedef itk::Image<PixelType, Dimension> InputImageType;
  79   typedef itk::Image< InternalPixelType, Dimension> InternalImageType; //labeling doesn't work with unsigned char?
  80
  81   //----------------------------------------------------------------------------------------------------
  82   // Typedef for Processing
  83   //----------------------------------------------------------------------------------------------------
  84   typedef itk::ImageFileReader<InternalImageType> FeatureReaderType;
  85   typedef itk::BinaryThresholdImageFilter<InputImageType, InternalImageType> InputBinarizeFilterType;
  86   typedef itk::BinaryThresholdImageFilter<InternalImageType, InternalImageType> InversionFilterType;
  87   typedef itk::ThresholdImageFilter<InternalImageType> ThresholdFilterType;
  88   typedef itk::ConnectedComponentImageFilter<InternalImageType, InternalImageType> ConnectFilterType;
  89   typedef itk::RelabelComponentImageFilter<InternalImageType, InternalImageType> RelabelFilterType;
  90   typedef itk::MirrorPadImageFilter<InternalImageType, InternalImageType> MirrorPadImageFilterType;
  91
  92
  93   typename InternalImageType::Pointer air;
  94   if (m_ArgsInfo.featureAir_given) {
  95     typename FeatureReaderType::Pointer featureReader =FeatureReaderType::New();
  96     featureReader->SetFileName(m_ArgsInfo.featureAir_arg);
  97     if (m_Verbose) std::cout<<"Reading the air feature image..."<<std::endl;
  98     featureReader->Update();
  99     air=featureReader->GetOutput();
 100   } else {
 101     if (m_Verbose) std::cout<<"Extracting the air feature image..."<<std::endl;
 102     //---------------------------------
 103     // Binarize the image
 104     //---------------------------------
 105     typename InputBinarizeFilterType::Pointer binarizeFilter=InputBinarizeFilterType::New();
 106     binarizeFilter->SetInput(input);
 107     binarizeFilter->SetLowerThreshold(static_cast<PixelType>(m_ArgsInfo.lowerThresholdAir_arg));
 108     binarizeFilter->SetUpperThreshold(static_cast<PixelType>(m_ArgsInfo.upperThresholdAir_arg));
 109     if (m_Verbose) std::cout<<"Binarizing the image using thresholds "<<m_ArgsInfo.lowerThresholdAir_arg
 110                               <<", "<<m_ArgsInfo.upperThresholdAir_arg<<"..."<<std::endl;
 111     binarizeFilter->Update();
 112
 113     //---------------------------------
 114     // Remove lungs (in place)
 115     //---------------------------------
 116     typedef itk::ImageRegionIterator<InternalImageType> IteratorType;
 117     IteratorType itAir(binarizeFilter->GetOutput(), binarizeFilter->GetOutput()->GetLargestPossibleRegion());
 118     IteratorType itLungs(lungs, binarizeFilter->GetOutput()->GetLargestPossibleRegion()); //lungs padded, used air region
 119     itAir.GoToBegin();
 120     itLungs.GoToBegin();
 121     while(!itAir.IsAtEnd()) {
 122       if(!itLungs.Get())     // The lungs are set to 0 at that stage
 123         itAir.Set(0);
 124       ++itAir;
 125       ++itLungs;
 126     }
 127
 128     //---------------------------------
 129     // Label the connected components
 130     //---------------------------------
 131     typename ConnectFilterType::Pointer connectFilter=ConnectFilterType::New();
 132     connectFilter->SetInput(binarizeFilter->GetOutput());
 133     connectFilter->SetBackgroundValue(0);
 134     connectFilter->SetFullyConnected(false);
 135     if (m_Verbose) std::cout<<"Labeling the connected components..."<<std::endl;
 136
 137     //---------------------------------
 138     // Sort the labels according to size
 139     //---------------------------------
 140     typename RelabelFilterType::Pointer relabelFilter=RelabelFilterType::New();
 141     relabelFilter->SetInput(connectFilter->GetOutput());
 142     if (m_Verbose) std::cout<<"Sorting the labels..."<<std::endl;
 143
 144     //---------------------------------
 145     // Keep the label
 146     //---------------------------------
 147     typename ThresholdFilterType::Pointer thresholdFilter=ThresholdFilterType::New();
 148     thresholdFilter->SetInput(relabelFilter->GetOutput());
 149     thresholdFilter->SetUpper(1);
 150     if (m_Verbose) std::cout<<"Keeping the first label..."<<std::endl;
 151     thresholdFilter->Update();
 152     air=thresholdFilter->GetOutput();
 153   }
 154
 155   //---------------------------------
 156   // Invert
 157   //---------------------------------
 158   typename InversionFilterType::Pointer inversionFilter=InversionFilterType::New();
 159   inversionFilter->SetInput(air);
 160   inversionFilter->SetLowerThreshold(0);
 161   inversionFilter->SetUpperThreshold(0);
 162   inversionFilter->SetInsideValue(1);
 163   inversionFilter->Update();
 164
 165   //---------------------------------
 166   // Mirror
 167   //---------------------------------
 168   typename MirrorPadImageFilterType::Pointer mirrorPadImageFilter=MirrorPadImageFilterType::New();
 169   mirrorPadImageFilter->SetInput(inversionFilter->GetOutput());
 170   typename InternalImageType::SizeType padSize;
 171   padSize.Fill(0);
 172   padSize[2]=input->GetLargestPossibleRegion().GetSize()[2];
 173   mirrorPadImageFilter->SetPadLowerBound(padSize);
 174   if (m_Verbose) std::cout<<"Mirroring the entire image along the CC axis..."<<std::endl;
 175   mirrorPadImageFilter->Update();
 176   air=mirrorPadImageFilter->GetOutput();
 177   writeImage<InternalImageType>(air,"/home/srit/tmp/air.mhd");
 178
 179   //---------------------------------
 180   // Pad
 181   //---------------------------------
 182   if(m_ArgsInfo.pad_flag) {
 183     typedef itk::ImageRegionIteratorWithIndex<InternalImageType> IteratorType;
 184     IteratorType it(air, air->GetLargestPossibleRegion());
 185     typename InternalImageType::IndexType index;
 186     while(!it.IsAtEnd()) {
 187       index=it.GetIndex();
 188       for (unsigned int i=0; i<Dimension; i++)
 189         if(index[i]==air->GetLargestPossibleRegion().GetIndex()[i]
 190             || index[i]==(unsigned int)air->GetLargestPossibleRegion().GetIndex()[i]+ (unsigned int) air->GetLargestPossibleRegion().GetSize()[i]-1)
 191           it.Set(0);
 192       ++it;
 193     }
 194   }
 195
 196   return air;
 197 }
 198
 199
 200 //-------------------------------------------------------------------
 201 // Bones
 202 //-------------------------------------------------------------------
 203 template <unsigned int Dimension, class  PixelType>
 204 typename itk::Image<MotionMaskGenericFilter::InternalPixelType, Dimension>::Pointer
 205 MotionMaskGenericFilter::GetBonesImage(typename itk::Image<PixelType, Dimension>::Pointer input )
 206 {
 207
 208   // ImageTypes
 209   typedef int InternalPixelType;
 210   typedef itk::Image<PixelType, Dimension> InputImageType;
 211   typedef itk::Image< InternalPixelType, Dimension> InternalImageType; //labeling doesn't work with unsigned char?
 212
 213   //----------------------------------------------------------------------------------------------------
 214   // Typedef for Processing
 215   //----------------------------------------------------------------------------------------------------
 216   typedef itk::ImageFileReader<InternalImageType> FeatureReaderType;
 217   typedef itk::BinaryThresholdImageFilter<InputImageType, InternalImageType> InputBinarizeFilterType;
 218   typedef itk::BinaryThresholdImageFilter<InternalImageType, InternalImageType> InversionFilterType;
 219   typedef itk::ThresholdImageFilter<InternalImageType> ThresholdFilterType;
 220   typedef itk::ConnectedComponentImageFilter<InternalImageType, InternalImageType> ConnectFilterType;
 221   typedef itk::RelabelComponentImageFilter<InternalImageType, InternalImageType> RelabelFilterType;
 222   typedef itk::MirrorPadImageFilter<InternalImageType, InternalImageType> MirrorPadImageFilterType;
 223
 224
 225   typename InternalImageType::Pointer bones;
 226   if (m_ArgsInfo.featureBones_given) {
 227     typename FeatureReaderType::Pointer featureReader =FeatureReaderType::New();
 228     featureReader->SetFileName(m_ArgsInfo.featureBones_arg);
 229     if (m_Verbose) std::cout<<"Reading the bones feature image..."<<std::endl;
 230     featureReader->Update();
 231     bones=featureReader->GetOutput();
 232   } else {
 233     if (m_Verbose) std::cout<<"Extracting the bones feature image..."<<std::endl;
 234     //---------------------------------
 235     // Binarize the image
 236     //---------------------------------
 237     typename InputBinarizeFilterType::Pointer binarizeFilter=InputBinarizeFilterType::New();
 238     binarizeFilter->SetInput(input);
 239     binarizeFilter->SetLowerThreshold(static_cast<PixelType>(m_ArgsInfo.lowerThresholdBones_arg));
 240     binarizeFilter->SetUpperThreshold(static_cast<PixelType>(m_ArgsInfo.upperThresholdBones_arg));
 241     if (m_Verbose) std::cout<<"Binarizing the image using thresholds "<<m_ArgsInfo.lowerThresholdBones_arg
 242                               <<", "<<m_ArgsInfo.upperThresholdBones_arg<<"..."<<std::endl;
 243
 244     //---------------------------------
 245     // Label the connected components
 246     //---------------------------------
 247     typename ConnectFilterType::Pointer connectFilter=ConnectFilterType::New();
 248     connectFilter->SetInput(binarizeFilter->GetOutput());
 249     connectFilter->SetBackgroundValue(0);
 250     connectFilter->SetFullyConnected(false);
 251     if (m_Verbose) std::cout<<"Labeling the connected components..."<<std::endl;
 252
 253     //---------------------------------
 254     // Sort the labels according to size
 255     //---------------------------------
 256     typename RelabelFilterType::Pointer relabelFilter=RelabelFilterType::New();
 257     relabelFilter->SetInput(connectFilter->GetOutput());
 258     if (m_Verbose) std::cout<<"Sorting the labels..."<<std::endl;
 259
 260     //---------------------------------
 261     // Keep the label
 262     //---------------------------------
 263     typename ThresholdFilterType::Pointer thresholdFilter=ThresholdFilterType::New();
 264     thresholdFilter->SetInput(relabelFilter->GetOutput());
 265     thresholdFilter->SetUpper(1);
 266     if (m_Verbose) std::cout<<"Keeping the first label..."<<std::endl;
 267     thresholdFilter->Update();
 268     bones=thresholdFilter->GetOutput();
 269
 270   }
 271
 272   //---------------------------------
 273   // Invert
 274   //---------------------------------
 275   typename InversionFilterType::Pointer inversionFilter=InversionFilterType::New();
 276   inversionFilter->SetInput(bones);
 277   inversionFilter->SetLowerThreshold(0);
 278   inversionFilter->SetUpperThreshold(0);
 279   inversionFilter->SetInsideValue(1);
 280   inversionFilter->Update();
 281
 282   //---------------------------------
 283   // Mirror
 284   //---------------------------------
 285   typename MirrorPadImageFilterType::Pointer mirrorPadImageFilter=MirrorPadImageFilterType::New();
 286   mirrorPadImageFilter->SetInput(inversionFilter->GetOutput());
 287   typename InternalImageType::SizeType padSize;
 288   padSize.Fill(0);
 289   padSize[2]=input->GetLargestPossibleRegion().GetSize()[2];
 290   mirrorPadImageFilter->SetPadLowerBound(padSize);
 291   if (m_Verbose) std::cout<<"Mirroring the entire image along the CC axis..."<<std::endl;
 292   mirrorPadImageFilter->Update();
 293   bones=mirrorPadImageFilter->GetOutput();
 294   // writeImage<InternalImageType>(bones,"/home/jef/tmp/bones.mhd");
 295
 296   return bones;
 297 }
 298
 299
 300
 301
 302 //-------------------------------------------------------------------
 303 // Lungs
 304 //-------------------------------------------------------------------
 305 template <unsigned int Dimension, class  PixelType>
 306 typename itk::Image<MotionMaskGenericFilter::InternalPixelType, Dimension>::Pointer
 307 MotionMaskGenericFilter::GetLungsImage(typename itk::Image<PixelType, Dimension>::Pointer input )
 308 {
 309   // ImageTypes
 310   typedef int InternalPixelType;
 311   typedef itk::Image<PixelType, Dimension> InputImageType;
 312   typedef itk::Image< InternalPixelType, Dimension> InternalImageType; //labeling doesn't work with unsigned char?
 313
 314   //----------------------------------------------------------------------------------------------------
 315   // Typedef for Processing
 316   //----------------------------------------------------------------------------------------------------
 317   typedef itk::ImageFileReader<InternalImageType> FeatureReaderType;
 318   typedef itk::BinaryThresholdImageFilter<InputImageType, InternalImageType> InputBinarizeFilterType;
 319   typedef itk::BinaryThresholdImageFilter<InternalImageType, InternalImageType> InversionFilterType;
 320   typedef itk::ThresholdImageFilter<InternalImageType> ThresholdFilterType;
 321   typedef itk::ConnectedComponentImageFilter<InternalImageType, InternalImageType> ConnectFilterType;
 322   typedef itk::RelabelComponentImageFilter<InternalImageType, InternalImageType> RelabelFilterType;
 323   typedef itk::MirrorPadImageFilter<InternalImageType, InternalImageType> MirrorPadImageFilterType;
 324
 325   typename InternalImageType::Pointer lungs;
 326   if (m_ArgsInfo.featureLungs_given) {
 327     typename FeatureReaderType::Pointer featureReader =FeatureReaderType::New();
 328     featureReader->SetFileName(m_ArgsInfo.featureLungs_arg);
 329     if (m_Verbose) std::cout<<"Reading the lungs feature image..."<<std::endl;
 330     featureReader->Update();
 331     lungs=featureReader->GetOutput();
 332   } else {
 333     if (m_Verbose) std::cout<<"Extracting the lungs feature image..."<<std::endl;
 334     //---------------------------------
 335     // Binarize the image
 336     //---------------------------------
 337     typename InputBinarizeFilterType::Pointer binarizeFilter=InputBinarizeFilterType::New();
 338     binarizeFilter->SetInput(input);
 339     binarizeFilter->SetLowerThreshold(static_cast<PixelType>(m_ArgsInfo.lowerThresholdLungs_arg));
 340     binarizeFilter->SetUpperThreshold(static_cast<PixelType>(m_ArgsInfo.upperThresholdLungs_arg));
 341     if (m_Verbose) std::cout<<"Binarizing the image using a threshold "<<m_ArgsInfo.lowerThresholdLungs_arg
 342                               <<", "<<m_ArgsInfo.upperThresholdLungs_arg<<"..."<<std::endl;
 343
 344     //---------------------------------
 345     // Label the connected components
 346     //---------------------------------
 347     typename ConnectFilterType::Pointer connectFilter=ConnectFilterType::New();
 348     connectFilter->SetInput(binarizeFilter->GetOutput());
 349     connectFilter->SetBackgroundValue(0);
 350     connectFilter->SetFullyConnected(true);
 351     if (m_Verbose) std::cout<<"Labeling the connected components..."<<std::endl;
 352
 353     //---------------------------------
 354     // Sort the labels according to size
 355     //---------------------------------
 356     typename RelabelFilterType::Pointer relabelFilter=RelabelFilterType::New();
 357     relabelFilter->SetInput(connectFilter->GetOutput());
 358     if (m_Verbose) std::cout<<"Sorting the labels..."<<std::endl;
 359     // writeImage<InternalImageType> (relabelFilter->GetOutput(), "/home/jef/tmp/labels.mhd");
 360
 361     //---------------------------------
 362     // Keep the label
 363     //---------------------------------
 364     typename ThresholdFilterType::Pointer thresholdFilter=ThresholdFilterType::New();
 365     thresholdFilter->SetInput(relabelFilter->GetOutput());
 366     thresholdFilter->SetLower(1);
 367     thresholdFilter->SetUpper(2);
 368     if (m_Verbose) std::cout<<"Keeping the first two labels..."<<std::endl;
 369     thresholdFilter->Update();
 370     lungs=thresholdFilter->GetOutput();
 371
 372   }
 373
 374
 375   //---------------------------------
 376   // Invert
 377   //---------------------------------
 378   typename InversionFilterType::Pointer inversionFilter=InversionFilterType::New();
 379   inversionFilter->SetInput(lungs);
 380   inversionFilter->SetLowerThreshold(0);
 381   inversionFilter->SetUpperThreshold(0);
 382   inversionFilter->SetInsideValue(1);
 383   inversionFilter->Update();
 384
 385   //---------------------------------
 386   // Mirror
 387   //---------------------------------
 388   typename MirrorPadImageFilterType::Pointer mirrorPadImageFilter=MirrorPadImageFilterType::New();
 389   mirrorPadImageFilter->SetInput(inversionFilter->GetOutput());
 390   typename InternalImageType::SizeType padSize;
 391   padSize.Fill(0);
 392   padSize[2]=input->GetLargestPossibleRegion().GetSize()[2];
 393   mirrorPadImageFilter->SetPadLowerBound(padSize);
 394   if (m_Verbose) std::cout<<"Mirroring the entire image along the CC axis..."<<std::endl;
 395   mirrorPadImageFilter->Update();
 396   lungs=mirrorPadImageFilter->GetOutput();
 397   // writeImage<InternalImageType>(lungs,"/home/jef/tmp/lungs.mhd");
 398
 399   return lungs;
 400 }
 401
 402 //-------------------------------------------------------------------
 403 // Resample
 404 //-------------------------------------------------------------------
 405 template <unsigned int Dimension, class  PixelType>
 406 typename itk::Image<MotionMaskGenericFilter::InternalPixelType, Dimension>::Pointer
 407 MotionMaskGenericFilter::Resample( typename itk::Image<InternalPixelType,Dimension>::Pointer input )
 408 {
 409
 410   typedef int InternalPixelType;
 411   typedef itk::Image<PixelType, Dimension> InputImageType;
 412   typedef itk::Image< InternalPixelType, Dimension> InternalImageType; //labeling doesn't work with unsigned char?
 413
 414   //---------------------------------
 415   // Resample to new spacing
 416   //---------------------------------
 417   typename InternalImageType::SizeType size_low;
 418   typename InternalImageType::SpacingType spacing_low;
 419   for (unsigned int i=0; i<Dimension; i++)
 420     if (m_ArgsInfo.spacing_given==Dimension)
 421       for (unsigned int i=0; i<Dimension; i++) {
 422         spacing_low[i]=m_ArgsInfo.spacing_arg[i];
 423         size_low[i]=ceil(input->GetLargestPossibleRegion().GetSize()[i]*input->GetSpacing()[i]/spacing_low[i]);
 424       }
 425     else
 426       for (unsigned int i=0; i<Dimension; i++) {
 427         spacing_low[i]=m_ArgsInfo.spacing_arg[0];
 428         size_low[i]=ceil(input->GetLargestPossibleRegion().GetSize()[i]*input->GetSpacing()[i]/spacing_low[i]);
 429       }
 430
 431   typename InternalImageType::PointType origin;
 432   input->TransformIndexToPhysicalPoint(input->GetLargestPossibleRegion().GetIndex(), origin);
 433   typedef itk::ResampleImageFilter<InternalImageType, InternalImageType> ResampleImageFilterType;
 434   typename ResampleImageFilterType::Pointer resampler =ResampleImageFilterType::New();
 435   typedef itk::NearestNeighborInterpolateImageFunction<InternalImageType, double> InterpolatorType;
 436   typename InterpolatorType::Pointer interpolator=InterpolatorType::New();
 437   resampler->SetInterpolator(interpolator);
 438   resampler->SetOutputOrigin(origin);
 439   resampler->SetOutputSpacing(spacing_low);
 440   resampler->SetSize(size_low);
 441   resampler->SetInput(input);
 442   resampler->Update();
 443   typename InternalImageType::Pointer output =resampler->GetOutput();
 444   return output;
 445 }
 446
 447
 448 //-------------------------------------------------------------------
 449 // Initialize ellips
 450 //-------------------------------------------------------------------
 451 template <unsigned int Dimension, class  PixelType>
 452 typename itk::Image<MotionMaskGenericFilter::InternalPixelType, Dimension>::Pointer
 453 MotionMaskGenericFilter::InitializeEllips( typename itk::Vector<double,Dimension> center, typename itk::Image<InternalPixelType,Dimension>::Pointer bones_low )
 454 {
 455
 456   // ImageTypes
 457   typedef int InternalPixelType;
 458   typedef itk::Image<PixelType, Dimension> InputImageType;
 459   typedef itk::Image< InternalPixelType, Dimension> InternalImageType; //labeling doesn't work with unsigned char?
 460   typedef itk::Image< unsigned char , Dimension> OutputImageType;
 461   typedef itk::ImageRegionIteratorWithIndex<InternalImageType> IteratorType;
 462   typedef clitk::SetBackgroundImageFilter<InternalImageType,InternalImageType, InternalImageType> SetBackgroundFilterType;
 463   typedef itk::LabelStatisticsImageFilter<InternalImageType,InternalImageType> LabelStatisticsImageFilterType;
 464   typedef itk::CastImageFilter<InternalImageType,OutputImageType> CastImageFilterType;
 465
 466
 467   typename InternalImageType::Pointer ellips;
 468   if (m_ArgsInfo.ellips_given || m_ArgsInfo.grownEllips_given || m_ArgsInfo.filledRibs_given) {
 469     if(m_ArgsInfo.ellips_given) {
 470       typedef itk::ImageFileReader<InternalImageType> FeatureReaderType;
 471       typename FeatureReaderType::Pointer featureReader = FeatureReaderType::New();
 472       featureReader->SetFileName(m_ArgsInfo.ellips_arg);
 473       featureReader->Update();
 474       ellips=featureReader->GetOutput();
 475     }
 476   } else {
 477     if(m_Verbose) {
 478       std::cout<<std::endl;
 479       std::cout<<"=========================================="<<std::endl;
 480       std::cout<<"||       Initializing ellipsoide        ||"<<std::endl;
 481       std::cout<<"=========================================="<<std::endl;
 482     }
 483
 484     //---------------------------------
 485     // Offset from center
 486     //---------------------------------
 487     typename itk::Vector<double,Dimension> offset;
 488     if(m_ArgsInfo.offset_given== Dimension) {
 489       for(unsigned int i=0; i<Dimension; i++)
 490         offset[i]=m_ArgsInfo.offset_arg[i];
 491     } else {
 492       offset.Fill(0.);
 493       offset[1]=-50;
 494     }
 495     itk::Vector<double,Dimension> centerEllips=center+offset;
 496     if (m_Verbose) {
 497       std::cout<<"Placing the center of the initial ellipsoide at (mm) "<<std::endl;
 498       std::cout<<centerEllips[0];
 499       for (unsigned int i=1; i<Dimension; i++)
 500         std::cout<<" "<<centerEllips[i];
 501       std::cout<<std::endl;
 502     }
 503
 504     //---------------------------------
 505     // The ellips
 506     //---------------------------------
 507     ellips=InternalImageType::New();
 508     ellips->SetRegions (bones_low->GetLargestPossibleRegion());
 509     ellips->SetOrigin(bones_low->GetOrigin());
 510     ellips->SetSpacing(bones_low->GetSpacing());
 511     ellips->Allocate();
 512     ellips->FillBuffer(0);
 513
 514     // Axes
 515     typename itk::Vector<double, Dimension> axes;
 516     if (m_ArgsInfo.axes_given==Dimension)
 517       for(unsigned int i=0; i<Dimension; i++)
 518         axes[i]=m_ArgsInfo.axes_arg[i];
 519     else {
 520       axes[0]=100;
 521       axes[1]=30;
 522       axes[2]=150;
 523     }
 524
 525     // Draw an ellips
 526     IteratorType itEllips(ellips, ellips->GetLargestPossibleRegion());
 527     itEllips.GoToBegin();
 528     typename InternalImageType::PointType point;
 529     typename InternalImageType::IndexType index;
 530     double distance;
 531
 532     if (m_Verbose) std::cout<<"Drawing the initial ellipsoide..."<<std::endl;
 533     while (!itEllips.IsAtEnd()) {
 534       index=itEllips.GetIndex();
 535       ellips->TransformIndexToPhysicalPoint(index, point);
 536       distance=0.0;
 537       for(unsigned int i=0; i<Dimension; i++)
 538         distance+=powf( ( (centerEllips[i]-point[i])/axes[i] ), 2);
 539
 540       if (distance<1)
 541         itEllips.Set(1);
 542       ++itEllips;
 543     }
 544   }
 545
 546
 547   //---------------------------------
 548   // Write the ellips
 549   //---------------------------------
 550   if (m_ArgsInfo.writeEllips_given) {
 551     typename CastImageFilterType::Pointer caster=CastImageFilterType::New();
 552     caster->SetInput(ellips);
 553     writeImage<OutputImageType>(caster->GetOutput(), m_ArgsInfo.writeEllips_arg, m_Verbose);
 554   }
 555
 556   return ellips;
 557
 558 }
 559
 560
 561 //-------------------------------------------------------------------
 562 // Update with the number of dimensions and the pixeltype
 563 //-------------------------------------------------------------------
 564 template <unsigned int Dimension, class  PixelType>
 565 void
 566 MotionMaskGenericFilter::UpdateWithDimAndPixelType()
 567 {
 568
 569   // ImageTypes
 570   typedef int InternalPixelType;
 571   typedef itk::Image<PixelType, Dimension> InputImageType;
 572   typedef itk::Image< InternalPixelType, Dimension> InternalImageType; //labeling doesn't work with unsigned char?
 573   typedef itk::Image< float, Dimension> LevelSetImageType; //labeling doesn't work with unsigned char?
 574   typedef itk::Image<unsigned char, Dimension> OutputImageType;
 575
 576
 577   //----------------------------------------------------------------------------------------------------
 578   // Typedef for Processing
 579   //----------------------------------------------------------------------------------------------------
 580   typedef itk::ImageFileReader<InternalImageType> FeatureReaderType;
 581   typedef itk::BinaryThresholdImageFilter<InputImageType, InternalImageType> InputBinarizeFilterType;
 582   typedef itk::BinaryThresholdImageFilter< LevelSetImageType,InternalImageType>    LevelSetBinarizeFilterType;
 583   typedef itk::BinaryThresholdImageFilter<InternalImageType, InternalImageType> InversionFilterType;
 584   typedef itk::ThresholdImageFilter<InternalImageType> ThresholdFilterType;
 585   typedef itk::ConnectedComponentImageFilter<InternalImageType, InternalImageType> ConnectFilterType;
 586   typedef itk::RelabelComponentImageFilter<InternalImageType, InternalImageType> RelabelFilterType;
 587   typedef itk::MirrorPadImageFilter<InternalImageType, InternalImageType> MirrorPadImageFilterType;
 588   typedef itk::NearestNeighborInterpolateImageFunction<InternalImageType, double> InterpolatorType;
 589   typedef itk::ResampleImageFilter<InternalImageType, InternalImageType> ResampleImageFilterType;
 590   typedef itk::ImageRegionIteratorWithIndex<InternalImageType> IteratorType;
 591   typedef itk::GeodesicActiveContourLevelSetImageFilter<LevelSetImageType, InternalImageType> LevelSetImageFilterType;
 592   typedef itk::SignedDanielssonDistanceMapImageFilter<InternalImageType, LevelSetImageType> DistanceMapImageFilterType;
 593   typedef clitk::SetBackgroundImageFilter<InternalImageType,InternalImageType, InternalImageType> SetBackgroundFilterType;
 594   typedef itk::LabelStatisticsImageFilter<InternalImageType,InternalImageType> LabelStatisticsImageFilterType;
 595   typedef itk::CastImageFilter<InternalImageType,OutputImageType> CastImageFilterType;
 596
 597
 598   // Read the input
 599   typedef itk::ImageFileReader<InputImageType> InputReaderType;
 600   typename InputReaderType::Pointer reader = InputReaderType::New();
 601   reader->SetFileName( m_InputFileName);
 602   reader->Update();
 603   typename InputImageType::Pointer input= reader->GetOutput();
 604
 605   //     // globals for avi
 606   //     unsigned int number=0;
 607
 608
 609   if(m_Verbose) {
 610     std::cout<<std::endl;
 611     std::cout<<"=========================================="<<std::endl;
 612     std::cout<<"||         Making feature images        ||"<<std::endl;
 613     std::cout<<"=========================================="<<std::endl;
 614   }
 615
 616   //--------------------------------------------------------------------------------
 617   // Lungs
 618   //-------------------------------------------------------------------------------
 619   typename InternalImageType::Pointer lungs = this->GetLungsImage<Dimension, PixelType>(input);
 620
 621   //-------------------------------------------------------------------------------
 622   // Air
 623   //-------------------------------------------------------------------------------
 624   typename InternalImageType::Pointer air = this->GetAirImage<Dimension, PixelType>(input, lungs);
 625
 626   //-------------------------------------------------------------------------------
 627   // Bones
 628   //-------------------------------------------------------------------------------
 629   typename InternalImageType::Pointer bones = this->GetBonesImage<Dimension, PixelType>(input);
 630
 631   //----------------------------------------------------------------------------------------------------
 632   // Write feature images
 633   //----------------------------------------------------------------------------------------------------
 634   if(m_ArgsInfo.writeFeature_given) {
 635     typename SetBackgroundFilterType::Pointer setBackgroundFilter = SetBackgroundFilterType::New();
 636     setBackgroundFilter->SetInput(air);
 637     setBackgroundFilter->SetInput2(bones);
 638     setBackgroundFilter->SetMaskValue(0);
 639     setBackgroundFilter->SetOutsideValue(2);
 640     setBackgroundFilter->Update();
 641     typename InternalImageType::Pointer bones_air =setBackgroundFilter->GetOutput();
 642
 643     typename SetBackgroundFilterType::Pointer setBackgroundFilter2 = SetBackgroundFilterType::New();
 644     setBackgroundFilter2->SetInput(bones_air);
 645     setBackgroundFilter2->SetInput2(lungs);
 646     setBackgroundFilter2->SetMaskValue(0);
 647     setBackgroundFilter2->SetOutsideValue(3);
 648     setBackgroundFilter2->Update();
 649     typename InternalImageType::Pointer lungs_bones_air =setBackgroundFilter2->GetOutput();
 650
 651     typename CastImageFilterType::Pointer caster=CastImageFilterType::New();
 652     caster->SetInput(lungs_bones_air);
 653     writeImage<OutputImageType>(caster->GetOutput(), m_ArgsInfo.writeFeature_arg, m_Verbose);
 654   }
 655
 656   //----------------------------------------------------------------------------------------------------
 657   // Low dimensional versions
 658   //----------------------------------------------------------------------------------------------------
 659   typename InternalImageType::Pointer bones_low =Resample<Dimension , PixelType>(bones);
 660   typename InternalImageType::Pointer lungs_low =Resample<Dimension , PixelType>(lungs);
 661   typename InternalImageType::Pointer air_low =Resample<Dimension , PixelType>(air);
 662
 663   //---------------------------------
 664   // Pad
 665   //---------------------------------
 666   if(m_ArgsInfo.pad_flag) {
 667     typedef itk::ImageRegionIteratorWithIndex<InternalImageType> IteratorType;
 668     IteratorType it(air_low, air_low->GetLargestPossibleRegion());
 669     typename InternalImageType::IndexType index;
 670     while(!it.IsAtEnd()) {
 671       index=it.GetIndex();
 672       for (unsigned int i=0; i<Dimension; i++)
 673         if(index[i]==air_low->GetLargestPossibleRegion().GetIndex()[i]
 674             || index[i]==(unsigned int)air_low->GetLargestPossibleRegion().GetIndex()[i]+ (unsigned int) air_low->GetLargestPossibleRegion().GetSize()[i]-1)
 675           it.Set(0);
 676       ++it;
 677     }
 678   }
 679
 680   //---------------------------------
 681   // Center
 682   //---------------------------------
 683   typename itk::Vector<double,Dimension> center;
 684   typedef itk::ImageMomentsCalculator<InternalImageType> MomentsCalculatorType;
 685   typename    MomentsCalculatorType::Pointer momentsCalculator=MomentsCalculatorType::New();
 686   momentsCalculator->SetImage(air);
 687   momentsCalculator->Compute();
 688   center=momentsCalculator->GetCenterOfGravity();
 689   if (m_Verbose) {
 690     std::cout<<"The center of gravity of the patient body is located at (mm) "<<std::endl;
 691     std::cout<<center[0];
 692     for (unsigned int i=1; i<Dimension; i++)
 693       std::cout<<" "<<center[i];
 694     std::cout<<std::endl;
 695   }
 696
 697   //----------------------------------------------------------------------------------------------------
 698   // Ellipsoide initialization
 699   //----------------------------------------------------------------------------------------------------
 700   typename InternalImageType::Pointer m_Ellips= InitializeEllips<Dimension, PixelType>(center,bones_low);
 701
 702   //----------------------------------------------------------------------------------------------------
 703   // Grow ellips
 704   //----------------------------------------------------------------------------------------------------
 705   typename InternalImageType::Pointer grownEllips;
 706   if (m_ArgsInfo.grownEllips_given || m_ArgsInfo.filledRibs_given) {
 707     if (m_ArgsInfo.grownEllips_given) {
 708
 709       typename FeatureReaderType::Pointer featureReader = FeatureReaderType::New();
 710       featureReader->SetFileName(m_ArgsInfo.grownEllips_arg);
 711       featureReader->Update();
 712       grownEllips=featureReader->GetOutput();
 713     }
 714   } else {
 715
 716     if(m_Verbose) {
 717       std::cout<<std::endl;
 718       std::cout<<"=========================================="<<std::endl;
 719       std::cout<<"||          Growing ellipsoide           ||"<<std::endl;
 720       std::cout<<"=========================================="<<std::endl;
 721     }
 722
 723     //---------------------------------
 724     // Detect abdomen
 725     //---------------------------------
 726     typename InternalImageType::PointType dPoint;
 727     if (m_ArgsInfo.detectionPoint_given) {
 728       for (unsigned int i=0; i<Dimension; i++)
 729         dPoint[i]=m_ArgsInfo.detectionPoint_arg[i];
 730     } else {
 731       typename InternalImageType::RegionType searchRegion=air->GetLargestPossibleRegion();
 732       typename InternalImageType::SizeType searchRegionSize = searchRegion.GetSize();
 733       typename InternalImageType::IndexType searchRegionIndex = searchRegion.GetIndex();
 734       searchRegionIndex[2]+=searchRegionSize[2]/2;
 735       searchRegionSize[2]=1;
 736       searchRegion.SetSize(searchRegionSize);
 737       searchRegion.SetIndex(searchRegionIndex);
 738       IteratorType itAbdomen(air, searchRegion);
 739       itAbdomen.GoToBegin();
 740
 741       typename InternalImageType::PointType aPoint;
 742       typename InternalImageType::IndexType aIndex;
 743
 744       if (m_Verbose) std::cout<<"Detecting the abdomen..."<<std::endl;
 745       while (!itAbdomen.IsAtEnd()) {
 746         if(itAbdomen.Get()) break;
 747         ++itAbdomen;
 748       }
 749       aIndex=itAbdomen.GetIndex();
 750       air->TransformIndexToPhysicalPoint(aIndex,aPoint);
 751       if (m_Verbose) std::cout<<"Detected the abdomen at "<<aPoint<<".."<<std::endl;
 752
 753
 754       //---------------------------------
 755       // Detect abdomen in additional images?
 756       //---------------------------------
 757       if (m_ArgsInfo.detectionPairs_given) {
 758         for (unsigned int i=0; i<m_ArgsInfo.detectionPairs_given; i++) {
 759           typename InternalImageType::Pointer airAdd;
 760           //---------------------------------
 761           // Read the input
 762           //--------------------------------
 763           typedef itk::ImageFileReader<InputImageType> InputReaderType;
 764           typename InputReaderType::Pointer reader = InputReaderType::New();
 765           reader->SetFileName( m_ArgsInfo.detectionPairs_arg[i]);
 766           reader->Update();
 767           typename InputImageType::Pointer additional= reader->GetOutput();
 768           if (m_Verbose) std::cout<<"Extracting the air from additional image "<< i<<"..."<<std::endl;
 769
 770           //---------------------------------
 771           // Binarize the image
 772           //---------------------------------
 773           typename InputBinarizeFilterType::Pointer binarizeFilter=InputBinarizeFilterType::New();
 774           binarizeFilter->SetInput(additional);
 775           binarizeFilter->SetLowerThreshold(static_cast<PixelType>(m_ArgsInfo.lowerThresholdAir_arg));
 776           binarizeFilter->SetUpperThreshold(static_cast<PixelType>(m_ArgsInfo.upperThresholdAir_arg));
 777           if (m_Verbose) std::cout<<"Binarizing the image using thresholds "<<m_ArgsInfo.lowerThresholdAir_arg
 778                                     <<", "<<m_ArgsInfo.upperThresholdAir_arg<<"..."<<std::endl;
 779
 780           //---------------------------------
 781           // Label the connected components
 782           //---------------------------------
 783           typename ConnectFilterType::Pointer connectFilter=ConnectFilterType::New();
 784           connectFilter->SetInput(binarizeFilter->GetOutput());
 785           connectFilter->SetBackgroundValue(0);
 786           connectFilter->SetFullyConnected(false);
 787           if (m_Verbose) std::cout<<"Labeling the connected components..."<<std::endl;
 788
 789           //---------------------------------
 790           // Sort the labels according to size
 791           //---------------------------------
 792           typename RelabelFilterType::Pointer relabelFilter=RelabelFilterType::New();
 793           relabelFilter->SetInput(connectFilter->GetOutput());
 794           if (m_Verbose) std::cout<<"Sorting the labels..."<<std::endl;
 795
 796           //---------------------------------
 797           // Keep the label
 798           //---------------------------------
 799           typename ThresholdFilterType::Pointer thresholdFilter=ThresholdFilterType::New();
 800           thresholdFilter->SetInput(relabelFilter->GetOutput());
 801           thresholdFilter->SetUpper(1);
 802           if (m_Verbose) std::cout<<"Keeping the first label..."<<std::endl;
 803           thresholdFilter->Update();
 804           airAdd=thresholdFilter->GetOutput();
 805
 806
 807           //---------------------------------
 808           // Invert
 809           //---------------------------------
 810           typename InversionFilterType::Pointer inversionFilter=InversionFilterType::New();
 811           inversionFilter->SetInput(airAdd);
 812           inversionFilter->SetLowerThreshold(0);
 813           inversionFilter->SetUpperThreshold(0);
 814           inversionFilter->SetInsideValue(1);
 815           inversionFilter->Update();
 816
 817           //---------------------------------
 818           // Mirror
 819           //---------------------------------
 820           typename MirrorPadImageFilterType::Pointer mirrorPadImageFilter=MirrorPadImageFilterType::New();
 821           mirrorPadImageFilter->SetInput(inversionFilter->GetOutput());
 822           typename InternalImageType::SizeType padSize;
 823           padSize.Fill(0);
 824           padSize[2]=input->GetLargestPossibleRegion().GetSize()[2];
 825           mirrorPadImageFilter->SetPadLowerBound(padSize);
 826           if (m_Verbose) std::cout<<"Mirroring the entire image along the CC axis..."<<std::endl;
 827           mirrorPadImageFilter->Update();
 828           airAdd=mirrorPadImageFilter->GetOutput();
 829
 830           //---------------------------------
 831           // Detect abdomen
 832           //---------------------------------
 833           typename InternalImageType::RegionType searchRegion=airAdd->GetLargestPossibleRegion();
 834           typename InternalImageType::SizeType searchRegionSize = searchRegion.GetSize();
 835           typename InternalImageType::IndexType searchRegionIndex = searchRegion.GetIndex();
 836           searchRegionIndex[2]+=searchRegionSize[2]/2;
 837           searchRegionSize[2]=1;
 838           searchRegion.SetSize(searchRegionSize);
 839           searchRegion.SetIndex(searchRegionIndex);
 840           IteratorType itAbdomen(airAdd, searchRegion);
 841           itAbdomen.GoToBegin();
 842
 843           typename InternalImageType::PointType additionalPoint;
 844           typename InternalImageType::IndexType aIndex;
 845
 846           if (m_Verbose) std::cout<<"Detecting the abdomen..."<<std::endl;
 847           while (!itAbdomen.IsAtEnd()) {
 848             if(itAbdomen.Get()) break;
 849             ++itAbdomen;
 850           }
 851           aIndex=itAbdomen.GetIndex();
 852           airAdd->TransformIndexToPhysicalPoint(aIndex,additionalPoint);
 853           if (m_Verbose) std::cout<<"Detected the abdomen in the additional image at "<<additionalPoint<<".."<<std::endl;
 854
 855           if(additionalPoint[1]< aPoint[1]) {
 856             aPoint=additionalPoint;
 857             if (m_Verbose) std::cout<<"Modifying the detected abdomen to "<<aPoint<<".."<<std::endl;
 858
 859           }
 860         }
 861       }
 862
 863
 864       // Determine the detection point
 865       dPoint.Fill(0.0);
 866       dPoint+=center;
 867       dPoint[1]=aPoint[1];
 868       if(m_ArgsInfo.offsetDetect_given==Dimension)
 869         for(unsigned int i=0; i <Dimension; i++)
 870           dPoint[i]+=m_ArgsInfo.offsetDetect_arg[i];
 871       else
 872         dPoint[1]+=-10;
 873
 874     }
 875     if (m_Verbose) std::cout<<"Setting the detection point to "<<dPoint<<".."<<std::endl;
 876
 877
 878     //---------------------------------
 879     // Pad the rib image and ellips image
 880     //---------------------------------
 881     typename InternalImageType::Pointer padded_ellips;
 882     typename InternalImageType::Pointer padded_bones_low;
 883
 884     // If detection point not inside the image: pad
 885     typename InternalImageType::IndexType dIndex;
 886     if (!bones_low->TransformPhysicalPointToIndex(dPoint, dIndex)) {
 887       typename InternalImageType::SizeType padSize;
 888       padSize.Fill(0);
 889       padSize[1]=abs(dIndex[1])+1;
 890       if (m_Verbose) std::cout<<"Padding the images with padding size "<<padSize<<" to include the detection point..."<<std::endl;
 891
 892       typename MirrorPadImageFilterType::Pointer padBonesFilter=MirrorPadImageFilterType::New();
 893       padBonesFilter->SetInput(bones_low);
 894       padBonesFilter->SetPadLowerBound(padSize);
 895       padBonesFilter->Update();
 896       padded_bones_low=padBonesFilter->GetOutput();
 897
 898       typename MirrorPadImageFilterType::Pointer padEllipsFilter=MirrorPadImageFilterType::New();
 899       padEllipsFilter->SetInput(m_Ellips);
 900       padEllipsFilter->SetPadLowerBound(padSize);
 901       padEllipsFilter->Update();
 902       padded_ellips=padEllipsFilter->GetOutput();
 903     }
 904
 905     else {
 906       padded_bones_low=bones_low;
 907       padded_ellips=m_Ellips;
 908     }
 909
 910
 911     //---------------------------------
 912     // Calculate distance map
 913     //---------------------------------
 914     typename DistanceMapImageFilterType::Pointer distanceMapImageFilter = DistanceMapImageFilterType::New();
 915     distanceMapImageFilter->SetInput(padded_ellips);
 916     distanceMapImageFilter->SetInsideIsPositive(false);
 917     if (m_Verbose) std::cout<<"Calculating the distance map..."<<std::endl;
 918     distanceMapImageFilter->Update();
 919
 920     //---------------------------------
 921     // Grow while monitoring detection point
 922     //---------------------------------
 923     typename LevelSetImageFilterType::Pointer levelSetFilter=LevelSetImageFilterType::New();
 924     levelSetFilter->SetInput(  distanceMapImageFilter->GetOutput() );
 925     levelSetFilter->SetFeatureImage( padded_bones_low );
 926     levelSetFilter->SetPropagationScaling( 1 );
 927     levelSetFilter->SetCurvatureScaling( m_ArgsInfo.curve1_arg );
 928     levelSetFilter->SetAdvectionScaling( 0 );
 929     levelSetFilter->SetMaximumRMSError( m_ArgsInfo.maxRMS1_arg );
 930     levelSetFilter->SetNumberOfIterations( m_ArgsInfo.iter1_arg );
 931     levelSetFilter->SetUseImageSpacing(true);
 932
 933     //  //---------------------------------
 934     //  // Script for making movie
 935     //  //---------------------------------
 936     //          std::string script="source /home/jef/thesis/presentations/20091014_JDD/videos/motionmask/make_motion_mask_movie_4mm.sh ";
 937
 938
 939     if (m_Verbose) std::cout<<"Starting level set segmentation..."<<std::endl;
 940     unsigned int totalNumberOfIterations=0;
 941     while (true) {
 942       levelSetFilter->Update();
 943       totalNumberOfIterations+=levelSetFilter->GetElapsedIterations();
 944
 945       if ( levelSetFilter->GetOutput()->GetPixel(dIndex) < 0 ) {
 946         if (m_Verbose) std::cout<<"Detection point reached!"<<std::endl;
 947         break;
 948       } else {
 949         if (m_Verbose) std::cout<<"Detection point not reached after "<<totalNumberOfIterations<<" iterations..."<<std::endl;
 950         levelSetFilter->SetInput(levelSetFilter->GetOutput());
 951         if(m_ArgsInfo.monitor_given) writeImage<LevelSetImageType>(levelSetFilter->GetOutput(), m_ArgsInfo.monitor_arg, m_Verbose);
 952
 953         //              // script
 954         //              std::ostringstream number_str;
 955         //              number_str << number;
 956         //              std::string param =  number_str.str();
 957         //              system((script+param).c_str());
 958         //              number+=1;
 959
 960       }
 961       if ( (totalNumberOfIterations> 10000) || (levelSetFilter->GetRMSChange()< m_ArgsInfo.maxRMS1_arg) ) break;
 962     }
 963
 964     // Output values
 965     if (m_Verbose) std::cout<<"Level set segmentation stopped after "<<totalNumberOfIterations<<" iterations..."<<std::endl;
 966     std::cout << "Max. RMS error: " << levelSetFilter->GetMaximumRMSError() << std::endl;
 967     std::cout << "RMS change: " << levelSetFilter->GetRMSChange() << std::endl;
 968
 969     // Threshold
 970     typename LevelSetBinarizeFilterType::Pointer thresholder = LevelSetBinarizeFilterType::New();
 971     thresholder->SetUpperThreshold( 0.0 );
 972     thresholder->SetOutsideValue( 0 );
 973     thresholder->SetInsideValue( 1 );
 974     thresholder->SetInput( levelSetFilter->GetOutput() );
 975     if (m_Verbose) std::cout<<"Thresholding the output level set..."<<std::endl;
 976     thresholder->Update();
 977     grownEllips=thresholder->GetOutput();
 978   }
 979
 980   //---------------------------------
 981   // Write the grown ellips
 982   //---------------------------------
 983   if (m_ArgsInfo.writeGrownEllips_given) {
 984     typename CastImageFilterType::Pointer caster=CastImageFilterType::New();
 985     caster->SetInput(grownEllips);
 986     writeImage<OutputImageType>(caster->GetOutput(), m_ArgsInfo.writeGrownEllips_arg, m_Verbose);
 987   }
 988
 989
 990   //----------------------------------------------------------------------------------------------------
 991   // Grow inside ribs
 992   //----------------------------------------------------------------------------------------------------
 993   typename InternalImageType::Pointer filledRibs;
 994   if (m_ArgsInfo.filledRibs_given) {
 995     typename FeatureReaderType::Pointer featureReader = FeatureReaderType::New();
 996     featureReader->SetFileName(m_ArgsInfo.filledRibs_arg);
 997     if (m_Verbose) std::cout<<"Reading filled ribs mask..."<<std::endl;
 998     featureReader->Update();
 999     filledRibs=featureReader->GetOutput();
1000   } else {
1001     if(m_Verbose) {
1002       std::cout<<std::endl;
1003       std::cout<<"=========================================="<<std::endl;
1004       std::cout<<"||        Filling the ribs image         ||"<<std::endl;
1005       std::cout<<"=========================================="<<std::endl;
1006     }
1007     //---------------------------------
1008     // Make feature image air+bones
1009     //---------------------------------
1010     typename SetBackgroundFilterType::Pointer setBackgroundFilter = SetBackgroundFilterType::New();
1011     setBackgroundFilter->SetInput(air_low);
1012     setBackgroundFilter->SetInput2(bones_low);
1013     setBackgroundFilter->SetMaskValue(0);
1014     setBackgroundFilter->SetOutsideValue(0);
1015     setBackgroundFilter->Update();
1016     typename InternalImageType::Pointer bones_air_low =setBackgroundFilter->GetOutput();
1017
1018     //---------------------------------
1019     // Resampling previous solution
1020     //---------------------------------
1021     if (m_Verbose) std::cout<<"Resampling previous solution..."<<std::endl;
1022     typename ResampleImageFilterType::Pointer resampler =ResampleImageFilterType::New();
1023     typename InternalImageType::PointType origin;
1024     bones_low->TransformIndexToPhysicalPoint(bones_low->GetLargestPossibleRegion().GetIndex(), origin);
1025     resampler->SetOutputOrigin(origin);
1026     resampler->SetOutputSpacing(bones_low->GetSpacing());
1027     typename InternalImageType::SizeType size_low= bones_low->GetLargestPossibleRegion().GetSize();
1028     resampler->SetSize(size_low);
1029     typename InterpolatorType::Pointer interpolator=InterpolatorType::New();
1030     resampler->SetInterpolator(interpolator);
1031     resampler->SetInput(grownEllips);
1032     resampler->Update();
1033     typename InternalImageType::Pointer grownEllips =resampler->GetOutput();
1034
1035
1036     //---------------------------------
1037     // Calculate Distance Map
1038     //---------------------------------
1039     typename DistanceMapImageFilterType::Pointer distanceMapImageFilter = DistanceMapImageFilterType::New();
1040     distanceMapImageFilter->SetInput(grownEllips);
1041     distanceMapImageFilter->SetInsideIsPositive(false);
1042     if (m_Verbose) std::cout<<"Calculating distance map..."<<std::endl;
1043     distanceMapImageFilter->Update();
1044
1045     //---------------------------------
1046     // Grow while monitoring lung volume coverage
1047     //---------------------------------
1048     typename LevelSetImageFilterType::Pointer levelSetFilter=LevelSetImageFilterType::New();
1049     levelSetFilter->SetInput(  distanceMapImageFilter->GetOutput() );
1050     levelSetFilter->SetFeatureImage( bones_air_low );
1051     levelSetFilter->SetPropagationScaling( 1 );
1052     levelSetFilter->SetCurvatureScaling( m_ArgsInfo.curve2_arg );
1053     levelSetFilter->SetAdvectionScaling( 0 );
1054     levelSetFilter->SetMaximumRMSError( m_ArgsInfo.maxRMS2_arg );
1055     levelSetFilter->SetNumberOfIterations( m_ArgsInfo.iter2_arg );
1056     levelSetFilter->SetUseImageSpacing(true);
1057
1058     //---------------------------------
1059     // Invert lung mask
1060     //---------------------------------
1061     typename InversionFilterType::Pointer invertor= InversionFilterType::New();
1062     invertor->SetInput(lungs_low);
1063     invertor->SetLowerThreshold(0);
1064     invertor->SetUpperThreshold(0);
1065     invertor->SetInsideValue(1);
1066     invertor->Update();
1067     typename InternalImageType::Pointer lungs_low_inv=invertor->GetOutput();
1068
1069     //---------------------------------
1070     // Calculate lung volume
1071     //---------------------------------
1072     typename LabelStatisticsImageFilterType::Pointer statisticsFilter= LabelStatisticsImageFilterType::New();
1073     statisticsFilter->SetInput(lungs_low_inv);
1074     statisticsFilter->SetLabelInput(lungs_low);
1075     statisticsFilter->Update();
1076     unsigned int volume=statisticsFilter->GetSum(0);
1077
1078     // Prepare threshold
1079     typename LevelSetBinarizeFilterType::Pointer thresholder = LevelSetBinarizeFilterType::New();
1080     thresholder->SetUpperThreshold(     0.0 );
1081     thresholder->SetOutsideValue(  0  );
1082     thresholder->SetInsideValue(  1 );
1083
1084
1085     // Start monitoring
1086     unsigned int totalNumberOfIterations=0;
1087     unsigned int coverage=0;
1088
1089
1090     //  //---------------------------------
1091     //  // Script for making movie
1092     //  //---------------------------------
1093     //          std::string script="source /home/jef/thesis/presentations/20091014_JDD/videos/motionmask/make_motion_mask_movie_4mm.sh ";
1094
1095     while (true) {
1096       levelSetFilter->Update();
1097       totalNumberOfIterations+=levelSetFilter->GetElapsedIterations();
1098
1099       thresholder->SetInput( levelSetFilter->GetOutput() );
1100       thresholder->Update();
1101       statisticsFilter->SetInput(thresholder->GetOutput());
1102       statisticsFilter->Update();
1103       coverage=statisticsFilter->GetSum(0);
1104
1105       // Compare the volumes
1106       if ( coverage >= m_ArgsInfo.fillingLevel_arg * volume/100 ) {
1107         if (m_Verbose) std::cout<<"Lungs filled for "<< m_ArgsInfo.fillingLevel_arg<<"% !"<<std::endl;
1108         break;
1109       } else {
1110         if (m_Verbose) std::cout<<"After "<<totalNumberOfIterations<<" iterations, lungs are covered for "
1111                                   <<(double)coverage/volume*100.0<<"%..."<<std::endl;
1112         levelSetFilter->SetInput(levelSetFilter->GetOutput());
1113         if(m_ArgsInfo.monitor_given) writeImage<LevelSetImageType>(levelSetFilter->GetOutput(), m_ArgsInfo.monitor_arg, m_Verbose);
1114
1115         //              // script
1116         //              std::ostringstream number_str;
1117         //              number_str << number;
1118         //              std::string param =  number_str.str();
1119         //              system((script+param).c_str());
1120         //              number+=1;
1121
1122       }
1123       if ( (totalNumberOfIterations> 30000) || (levelSetFilter->GetRMSChange()< m_ArgsInfo.maxRMS2_arg) ) break;
1124     }
1125
1126     // Output values
1127     if (m_Verbose) std::cout<<"Level set segmentation stopped after "<<totalNumberOfIterations<<" iterations..."<<std::endl;
1128     std::cout << "Max. RMS error: " << levelSetFilter->GetMaximumRMSError() << std::endl;
1129     std::cout << "RMS change: " << levelSetFilter->GetRMSChange() << std::endl;
1130
1131     // Threshold
1132     thresholder->SetInput( levelSetFilter->GetOutput() );
1133     thresholder->Update();
1134     filledRibs=thresholder->GetOutput();
1135     // writeImage<InternalImageType>(filledRibs,"/home/jef/tmp/filled_ribs_before.mhd", m_Verbose);
1136
1137   }
1138
1139   //---------------------------------
1140   // Write the filled ribs
1141   //---------------------------------
1142   if (m_ArgsInfo.writeFilledRibs_given) {
1143     typename CastImageFilterType::Pointer caster=CastImageFilterType::New();
1144     caster->SetInput(filledRibs);
1145     writeImage<OutputImageType>(caster->GetOutput(), m_ArgsInfo.writeFilledRibs_arg, m_Verbose);
1146   }
1147
1148
1149   //----------------------------------------------------------------------------------------------------
1150   // Collapse to the lungs
1151   //----------------------------------------------------------------------------------------------------
1152   if(m_Verbose) {
1153     std::cout<<std::endl;
1154     std::cout<<"=========================================="<<std::endl;
1155     std::cout<<"||     Collapsing to the lung image     ||"<<std::endl;
1156     std::cout<<"=========================================="<<std::endl;
1157   }
1158
1159   //---------------------------------
1160   // Make feature image air+bones
1161   //---------------------------------
1162   if (m_Verbose) std::cout<<"Making feature images..."<<std::endl;
1163   typename SetBackgroundFilterType::Pointer setBackgroundFilter = SetBackgroundFilterType::New();
1164   setBackgroundFilter->SetInput(air);
1165   setBackgroundFilter->SetInput2(bones);
1166   setBackgroundFilter->SetMaskValue(0);
1167   setBackgroundFilter->SetOutsideValue(0);
1168   setBackgroundFilter->Update();
1169   typename InternalImageType::Pointer bones_air =setBackgroundFilter->GetOutput();
1170
1171   typename SetBackgroundFilterType::Pointer setBackgroundFilter2 = SetBackgroundFilterType::New();
1172   setBackgroundFilter2->SetInput(bones_air);
1173   setBackgroundFilter2->SetInput2(lungs);
1174   setBackgroundFilter2->SetMaskValue(0);
1175   setBackgroundFilter2->SetOutsideValue(0);
1176   setBackgroundFilter2->Update();
1177   typename InternalImageType::Pointer lungs_bones_air =setBackgroundFilter2->GetOutput();
1178
1179   //---------------------------------
1180   // Prepare previous solution
1181   //---------------------------------
1182   if (m_Verbose) std::cout<<"Resampling previous solution..."<<std::endl;
1183   typename ResampleImageFilterType::Pointer resampler =ResampleImageFilterType::New();
1184   typedef itk::NearestNeighborInterpolateImageFunction<InternalImageType, double> InterpolatorType;
1185   typename InterpolatorType::Pointer interpolator=InterpolatorType::New();
1186   resampler->SetOutputStartIndex(bones->GetLargestPossibleRegion().GetIndex());
1187   resampler->SetInput(filledRibs);
1188   resampler->SetInterpolator(interpolator);
1189   resampler->SetOutputParametersFromImage(bones);
1190   resampler->Update();
1191   filledRibs =resampler->GetOutput();
1192
1193   if (m_Verbose) std::cout<<"Setting lungs to 1..."<<std::endl;
1194   typename SetBackgroundFilterType::Pointer setBackgroundFilter3 = SetBackgroundFilterType::New();
1195   setBackgroundFilter3->SetInput(filledRibs);
1196   setBackgroundFilter3->SetInput2(lungs);
1197   setBackgroundFilter3->SetMaskValue(0);
1198   setBackgroundFilter3->SetOutsideValue(1);
1199   setBackgroundFilter3->Update();
1200   filledRibs=setBackgroundFilter3->GetOutput();
1201
1202   if (m_Verbose) std::cout<<"Removing overlap with bones..."<<std::endl;
1203   typename SetBackgroundFilterType::Pointer setBackgroundFilter4 = SetBackgroundFilterType::New();
1204   setBackgroundFilter4->SetInput(filledRibs);
1205   setBackgroundFilter4->SetInput2(bones);
1206   setBackgroundFilter4->SetMaskValue(0);
1207   setBackgroundFilter4->SetOutsideValue(0);
1208   setBackgroundFilter4->Update();
1209   filledRibs =setBackgroundFilter4->GetOutput();
1210   //   writeImage<InternalImageType>(filledRibs,"/home/jef/tmp/filledRibs_pp.mhd");
1211   //---------------------------------
1212   // Calculate Distance Map
1213   //---------------------------------
1214   //  typedef  itk::ApproximateSignedDistanceMapImageFilter <InternalImageType, LevelSetImageType> DistanceMapImageFilterType2;
1215   typename DistanceMapImageFilterType::Pointer distanceMapImageFilter = DistanceMapImageFilterType::New();
1216   distanceMapImageFilter->SetInput(filledRibs);
1217   distanceMapImageFilter->SetInsideIsPositive(false);
1218   // distanceMapImageFilter->SetInsideValue(0);
1219 //     distanceMapImageFilter->SetOutsideValue(1);
1220   if (m_Verbose) std::cout<<"Calculating distance map..."<<std::endl;
1221   distanceMapImageFilter->Update();
1222
1223   //---------------------------------
1224   // Collapse
1225   //---------------------------------
1226   typename LevelSetImageFilterType::Pointer levelSetFilter=LevelSetImageFilterType::New();
1227   levelSetFilter->SetInput(  distanceMapImageFilter->GetOutput() );
1228   levelSetFilter->SetFeatureImage( lungs_bones_air );
1229   levelSetFilter->SetPropagationScaling(m_ArgsInfo.prop3_arg);
1230   levelSetFilter->SetCurvatureScaling( m_ArgsInfo.curve3_arg );
1231   levelSetFilter->SetAdvectionScaling( 0 );
1232   levelSetFilter->SetMaximumRMSError( m_ArgsInfo.maxRMS3_arg );
1233   levelSetFilter->SetNumberOfIterations( m_ArgsInfo.iter3_arg );
1234   levelSetFilter->SetUseImageSpacing(true);
1235
1236   //     //script
1237   //     std::string script="source /home/jef/thesis/presentations/20091014_JDD/videos/motionmask/make_motion_mask_movie_4mm.sh ";
1238
1239   // Start monitoring
1240   if (m_Verbose) std::cout<<"Starting the levelset..."<<std::endl;
1241   unsigned int totalNumberOfIterations=0;
1242   while (true) {
1243     levelSetFilter->Update();
1244
1245     // monitor state
1246     totalNumberOfIterations+=levelSetFilter->GetElapsedIterations();
1247     levelSetFilter->SetInput(levelSetFilter->GetOutput());
1248     if(m_ArgsInfo.monitor_given) writeImage<LevelSetImageType>(levelSetFilter->GetOutput(), m_ArgsInfo.monitor_arg);
1249     std::cout << "After "<<totalNumberOfIterations<<" iteration the RMS change is " << levelSetFilter->GetRMSChange() <<"..."<< std::endl;
1250
1251     //  // script
1252     //  std::ostringstream number_str;
1253     //  number_str << number;
1254     //  std::string param =  number_str.str();
1255     //  system((script+param).c_str());
1256     //  number+=1;
1257
1258     // stopping condition
1259     if ( (totalNumberOfIterations>= (unsigned int) m_ArgsInfo.maxIter3_arg) || (levelSetFilter->GetRMSChange()< m_ArgsInfo.maxRMS3_arg) ) break;
1260     levelSetFilter->SetNumberOfIterations( std::min ((unsigned int) m_ArgsInfo.iter3_arg, (unsigned int) m_ArgsInfo.maxIter3_arg-totalNumberOfIterations ) );
1261   }
1262
1263   // Output values
1264   if (m_Verbose) {
1265     std::cout<<"Level set segmentation stopped after "<<totalNumberOfIterations<<" iterations..."<<std::endl;
1266     std::cout << "Max. RMS error: " << levelSetFilter->GetMaximumRMSError() << std::endl;
1267     std::cout << "RMS change: " << levelSetFilter->GetRMSChange() << std::endl;
1268   }
1269
1270   // Threshold
1271   typedef itk::BinaryThresholdImageFilter< LevelSetImageType,InternalImageType>    LevelSetBinarizeFilterType;
1272   typename LevelSetBinarizeFilterType::Pointer thresholder = LevelSetBinarizeFilterType::New();
1273   thresholder->SetUpperThreshold( 0.0 );
1274   thresholder->SetOutsideValue( 0 );
1275   thresholder->SetInsideValue( 1 );
1276   thresholder->SetInput( levelSetFilter->GetOutput() );
1277   thresholder->Update();
1278   typename InternalImageType::Pointer output = thresholder->GetOutput();
1279
1280
1281   //----------------------------------------------------------------------------------------------------
1282   // Prepare the output
1283   //----------------------------------------------------------------------------------------------------
1284
1285   //---------------------------------
1286   // Set Air to zero
1287   //---------------------------------
1288   if (m_Verbose) std::cout<<"Removing overlap mask with air..."<<std::endl;
1289   typename SetBackgroundFilterType::Pointer setBackgroundFilter5 = SetBackgroundFilterType::New();
1290   setBackgroundFilter5->SetInput(output);
1291   setBackgroundFilter5->SetInput2(air);
1292   setBackgroundFilter5->SetMaskValue(0);
1293   setBackgroundFilter5->SetOutsideValue(0);
1294   setBackgroundFilter5->Update();
1295   output=setBackgroundFilter5->GetOutput();
1296
1297   //---------------------------------
1298   // Open & close  to cleanup
1299   //---------------------------------
1300   if ( m_ArgsInfo.openClose_flag) {
1301
1302     //---------------------------------
1303     // Structuring element
1304     //---------------------------------
1305     typedef itk::BinaryBallStructuringElement<InternalPixelType,Dimension > KernelType;
1306     KernelType structuringElement;
1307     structuringElement.SetRadius(1);
1308     structuringElement.CreateStructuringElement();
1309
1310     //---------------------------------
1311     // Open
1312     //---------------------------------
1313     typedef itk::BinaryMorphologicalOpeningImageFilter<InternalImageType, InternalImageType , KernelType> OpenFilterType;
1314     typename OpenFilterType::Pointer openFilter = OpenFilterType::New();
1315     openFilter->SetInput(output);
1316     openFilter->SetBackgroundValue(0);
1317     openFilter->SetForegroundValue(1);
1318     openFilter->SetKernel(structuringElement);
1319     if(m_Verbose) std::cout<<"Opening the output image..."<<std::endl;
1320
1321     //---------------------------------
1322     // Close
1323     //---------------------------------
1324     typedef itk::BinaryMorphologicalClosingImageFilter<InternalImageType, InternalImageType , KernelType> CloseFilterType;
1325     typename CloseFilterType::Pointer closeFilter = CloseFilterType::New();
1326     closeFilter->SetInput(openFilter->GetOutput());
1327     closeFilter->SetSafeBorder(true);
1328     closeFilter->SetForegroundValue(1);
1329     closeFilter->SetKernel(structuringElement);
1330     if(m_Verbose) std::cout<<"Closing the output image..."<<std::endl;
1331     closeFilter->Update();
1332     output=closeFilter->GetOutput();
1333
1334   }
1335   //writeImage<InternalImageType>(output,"/home/jef/tmp/mm_double.mhd");
1336
1337   // Extract the upper part
1338   typedef  itk::CropImageFilter<InternalImageType, InternalImageType> CropImageFilterType;
1339   typename CropImageFilterType::Pointer cropFilter=CropImageFilterType::New();
1340   cropFilter->SetInput(output);
1341   typename InputImageType::SizeType lSize, uSize;
1342   uSize.Fill(0);
1343   lSize.Fill(0);
1344   lSize[2]=input->GetLargestPossibleRegion().GetSize()[2];
1345   cropFilter->SetLowerBoundaryCropSize(lSize);
1346   cropFilter->SetUpperBoundaryCropSize(uSize);
1347   cropFilter->Update();
1348
1349   // Output
1350   typename CastImageFilterType::Pointer caster=CastImageFilterType::New();
1351   caster->SetInput(cropFilter->GetOutput());
1352   writeImage<OutputImageType>(caster->GetOutput(), m_ArgsInfo.output_arg, m_Verbose);
1353
1354 }
1355
1356 }//end clitk
1357
1358 #endif //#define clitkMotionMaskGenericFilter_txx