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