itk/clitkReconstructThroughDilationImageFilter.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 clitkReconstructThroughDilationImageFilter_txx
  19 #define clitkReconstructThroughDilationImageFilter_txx
  20
  21 /* =================================================
  22  * @file   clitkReconstructThroughDilationImageFilter.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<class InputImageType, class OutputImageType>
  38   ReconstructThroughDilationImageFilter<InputImageType, OutputImageType>::ReconstructThroughDilationImageFilter()
  39   {
  40    m_Verbose=false;
  41    m_BackgroundValue=0;
  42    m_ForegroundValue=1;
  43    m_ErosionPaddingValue=static_cast<InputPixelType>(-1);
  44    for (unsigned int i=0; i<InputImageDimension; i++)
  45      m_Radius[i]=1;
  46    m_MaximumNumberOfLabels=10;
  47   }
  48
  49
  50   //-------------------------------------------------------------------
  51   // Update with the number of dimensions and the pixeltype
  52   //-------------------------------------------------------------------
  53   template<class InputImageType, class  OutputImageType>
  54   void
  55   ReconstructThroughDilationImageFilter<InputImageType, OutputImageType>::GenerateData()
  56   {
  57
  58     //---------------------------------
  59     // Typedefs
  60     //---------------------------------
  61
  62     // Internal type
  63     typedef itk::Image<InternalPixelType, InputImageDimension> InternalImageType;
  64
  65     // Filters used
  66     typedef itk::CastImageFilter<InputImageType, InternalImageType> InputCastImageFilterType;
  67     typedef itk::ThresholdImageFilter<InternalImageType> InputThresholdImageFilterType;
  68     typedef itk::StatisticsImageFilter<InternalImageType> StatisticsImageFilterType;
  69     typedef itk::BinaryBallStructuringElement<InternalPixelType,InputImageDimension > KernelType;
  70     typedef clitk::ConditionalBinaryDilateImageFilter<InternalImageType, InternalImageType , KernelType> ConditionalBinaryDilateImageFilterType;
  71     typedef itk::Testing::ComparisonImageFilter<InternalImageType, InternalImageType> DifferenceImageFilterType;
  72     typedef itk::CastImageFilter<InternalImageType, OutputImageType> OutputCastImageFilterType;
  73     typedef clitk::SetBackgroundImageFilter<InternalImageType, InternalImageType, InternalImageType> SetBackgroundImageFilterType;
  74
  75     //---------------------------------
  76     // Cast
  77     //---------------------------------
  78     typename InputCastImageFilterType::Pointer castImageFilter=InputCastImageFilterType::New();
  79     castImageFilter->SetInput(this->GetInput());
  80     castImageFilter->Update();
  81
  82     //---------------------------------
  83     // Threshold
  84     //---------------------------------
  85     typename InputThresholdImageFilterType::Pointer thresholdImageFilter=InputThresholdImageFilterType::New();
  86     thresholdImageFilter->SetInput(castImageFilter->GetOutput());
  87     thresholdImageFilter->ThresholdAbove(m_MaximumNumberOfLabels);
  88     thresholdImageFilter->SetOutsideValue(m_ForegroundValue);
  89     if(m_Verbose) std::cout<<"Thresholding the input to "<<m_MaximumNumberOfLabels<<" labels ..."<<std::endl;
  90     thresholdImageFilter->Update();
  91
  92     //---------------------------------
  93     // Set -1 to padding value
  94     //---------------------------------
  95     typename SetBackgroundImageFilterType::Pointer setBackgroundFilter =SetBackgroundImageFilterType::New();
  96     setBackgroundFilter->SetInput(thresholdImageFilter->GetOutput());
  97     setBackgroundFilter->SetInput2(castImageFilter->GetOutput());
  98     setBackgroundFilter->SetMaskValue(m_ErosionPaddingValue);
  99     setBackgroundFilter->SetOutsideValue(-1);
 100     if(m_Verbose) std::cout<<"Setting the eroded region from "<<m_ErosionPaddingValue<<" to -1..."<<std::endl;
 101
 102
 103     //---------------------------------
 104     // Count the initial labels
 105     //---------------------------------
 106     typename StatisticsImageFilterType::Pointer inputStatisticsImageFilter=StatisticsImageFilterType::New();
 107     inputStatisticsImageFilter->SetInput(setBackgroundFilter->GetOutput());
 108     if(m_Verbose) std::cout<<"Counting the initial labels..."<<std::endl;
 109     inputStatisticsImageFilter->Update();
 110     unsigned int initialNumberOfLabels= inputStatisticsImageFilter->GetMaximum();
 111     if(m_Verbose) std::cout<<"The input contained "<<initialNumberOfLabels<<" disctictive label(s)..."<<std::endl;
 112     unsigned int numberOfConsideredLabels=std::min(initialNumberOfLabels, m_MaximumNumberOfLabels);
 113     if(m_Verbose) std::cout<<"Performing dilation the first "<<numberOfConsideredLabels<<" disctictive labels..."<<std::endl;
 114
 115     //---------------------------------
 116     // Dilate while change
 117     //---------------------------------
 118     typename itk::NumericTraits<InputPixelType>::AccumulateType difference=1;
 119     typename InternalImageType::Pointer labelImage=inputStatisticsImageFilter->GetOutput();
 120     typename InternalImageType::Pointer oldLabelImage=inputStatisticsImageFilter->GetOutput();
 121
 122     // element
 123     KernelType structuringElement;
 124     structuringElement.SetRadius(m_Radius);
 125     structuringElement.CreateStructuringElement();
 126
 127     while( difference)
 128       {
 129         // Dilate all labels once
 130         for ( int label=0; label<(int)numberOfConsideredLabels+1;label++)
 131           if ( m_BackgroundValue != label)
 132             {
 133               typename ConditionalBinaryDilateImageFilterType::Pointer dilateFilter=ConditionalBinaryDilateImageFilterType::New();
 134               dilateFilter->SetBoundaryToForeground(false);
 135               dilateFilter->SetKernel(structuringElement);
 136               dilateFilter->SetBackgroundValue (-1);
 137               dilateFilter->SetInput (labelImage);
 138               dilateFilter->SetForegroundValue (label);
 139               if(m_Verbose) std::cout<<"Dilating the label "<<label<<"..."<<std::endl;
 140               dilateFilter->Update();
 141               labelImage=dilateFilter->GetOutput();
 142             }
 143
 144         // Difference with previous labelImage
 145         typename DifferenceImageFilterType::Pointer differenceFilter=DifferenceImageFilterType::New();
 146         differenceFilter->SetValidInput(oldLabelImage);
 147         differenceFilter->SetTestInput(labelImage);
 148         differenceFilter->Update();
 149         difference =differenceFilter->GetTotalDifference();
 150         if(m_Verbose) std::cout<<"The change in this iteration was "<<difference<<"..."<<std::endl;
 151         oldLabelImage=labelImage;
 152       }
 153
 154     //---------------------------------
 155     // Set -1 to padding value
 156     //---------------------------------
 157     typename SetBackgroundImageFilterType::Pointer setBackgroundFilter2 =SetBackgroundImageFilterType::New();
 158     setBackgroundFilter2->SetInput(labelImage);
 159     setBackgroundFilter2->SetInput2(labelImage);
 160     setBackgroundFilter2->SetMaskValue(-1);
 161     setBackgroundFilter2->SetOutsideValue(m_ErosionPaddingValue);
 162     if(m_Verbose) std::cout<<"Setting the eroded region to "<<m_ErosionPaddingValue<<"..."<<std::endl;
 163
 164     //---------------------------------
 165     // Cast
 166     //---------------------------------
 167     typename OutputCastImageFilterType::Pointer outputCastImageFilter=OutputCastImageFilterType::New();
 168     outputCastImageFilter->SetInput(setBackgroundFilter2->GetOutput());
 169     if(m_Verbose) std::cout<<"Casting the output..."<<std::endl;
 170     outputCastImageFilter->Update();
 171
 172     //---------------------------------
 173     // SetOutput
 174     //---------------------------------
 175     this->SetNthOutput(0, outputCastImageFilter->GetOutput());
 176
 177
 178   }
 179
 180
 181 }//end clitk
 182
 183 #endif //#define clitkReconstructThroughDilationImageFilter_txx