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