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