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