segmentation/clitkMotionMaskGenericFilter.txx

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