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 #if ITK_VERSION_MAJOR >= 4
  72     typedef itk::Testing::ComparisonImageFilter<InternalImageType, InternalImageType> DifferenceImageFilterType;
  73 #else
  74     typedef itk::DifferenceImageFilter<InternalImageType, InternalImageType> DifferenceImageFilterType;
  75 #endif
  76     typedef itk::CastImageFilter<InternalImageType, OutputImageType> OutputCastImageFilterType;
  77     typedef clitk::SetBackgroundImageFilter<InternalImageType, InternalImageType, InternalImageType> SetBackgroundImageFilterType;
  78
  79     //---------------------------------
  80     // Cast
  81     //---------------------------------
  82     typename InputCastImageFilterType::Pointer castImageFilter=InputCastImageFilterType::New();
  83     castImageFilter->SetInput(this->GetInput());
  84     castImageFilter->Update();
  85
  86     //---------------------------------
  87     // Threshold
  88     //---------------------------------
  89     typename InputThresholdImageFilterType::Pointer thresholdImageFilter=InputThresholdImageFilterType::New();
  90     thresholdImageFilter->SetInput(castImageFilter->GetOutput());
  91     thresholdImageFilter->ThresholdAbove(m_MaximumNumberOfLabels);
  92     thresholdImageFilter->SetOutsideValue(m_ForegroundValue);
  93     if(m_Verbose) std::cout<<"Thresholding the input to "<<m_MaximumNumberOfLabels<<" labels ..."<<std::endl;
  94     thresholdImageFilter->Update();
  95
  96     //---------------------------------
  97     // Set -1 to padding value
  98     //---------------------------------
  99     typename SetBackgroundImageFilterType::Pointer setBackgroundFilter =SetBackgroundImageFilterType::New();
 100     setBackgroundFilter->SetInput(thresholdImageFilter->GetOutput());
 101     setBackgroundFilter->SetInput2(castImageFilter->GetOutput());
 102     setBackgroundFilter->SetMaskValue(m_ErosionPaddingValue);
 103     setBackgroundFilter->SetOutsideValue(-1);
 104     if(m_Verbose) std::cout<<"Setting the eroded region from "<<m_ErosionPaddingValue<<" to -1..."<<std::endl;
 105
 106
 107     //---------------------------------
 108     // Count the initial labels
 109     //---------------------------------
 110     typename StatisticsImageFilterType::Pointer inputStatisticsImageFilter=StatisticsImageFilterType::New();
 111     inputStatisticsImageFilter->SetInput(setBackgroundFilter->GetOutput());
 112     if(m_Verbose) std::cout<<"Counting the initial labels..."<<std::endl;
 113     inputStatisticsImageFilter->Update();
 114     unsigned int initialNumberOfLabels= inputStatisticsImageFilter->GetMaximum();
 115     if(m_Verbose) std::cout<<"The input contained "<<initialNumberOfLabels<<" disctictive label(s)..."<<std::endl;
 116     unsigned int numberOfConsideredLabels=std::min(initialNumberOfLabels, m_MaximumNumberOfLabels);
 117     if(m_Verbose) std::cout<<"Performing dilation the first "<<numberOfConsideredLabels<<" disctictive labels..."<<std::endl;
 118
 119     //---------------------------------
 120     // Dilate while change
 121     //---------------------------------
 122     typename itk::NumericTraits<InputPixelType>::AccumulateType difference=1;
 123     typename InternalImageType::Pointer labelImage=inputStatisticsImageFilter->GetOutput();
 124     typename InternalImageType::Pointer oldLabelImage=inputStatisticsImageFilter->GetOutput();
 125
 126     // element
 127     KernelType structuringElement;
 128     structuringElement.SetRadius(m_Radius);
 129     structuringElement.CreateStructuringElement();
 130
 131     while( difference)
 132       {
 133         // Dilate all labels once
 134         for ( int label=0; label<(int)numberOfConsideredLabels+1;label++)
 135           if ( m_BackgroundValue != label)
 136             {
 137               typename ConditionalBinaryDilateImageFilterType::Pointer dilateFilter=ConditionalBinaryDilateImageFilterType::New();
 138               dilateFilter->SetBoundaryToForeground(false);
 139               dilateFilter->SetKernel(structuringElement);
 140               dilateFilter->SetBackgroundValue (-1);
 141               dilateFilter->SetInput (labelImage);
 142               dilateFilter->SetForegroundValue (label);
 143               if(m_Verbose) std::cout<<"Dilating the label "<<label<<"..."<<std::endl;
 144               dilateFilter->Update();
 145               labelImage=dilateFilter->GetOutput();
 146             }
 147
 148         // Difference with previous labelImage
 149         typename DifferenceImageFilterType::Pointer differenceFilter=DifferenceImageFilterType::New();
 150         differenceFilter->SetValidInput(oldLabelImage);
 151         differenceFilter->SetTestInput(labelImage);
 152         differenceFilter->Update();
 153         difference =differenceFilter->GetTotalDifference();
 154         if(m_Verbose) std::cout<<"The change in this iteration was "<<difference<<"..."<<std::endl;
 155         oldLabelImage=labelImage;
 156       }
 157
 158     //---------------------------------
 159     // Set -1 to padding value
 160     //---------------------------------
 161     typename SetBackgroundImageFilterType::Pointer setBackgroundFilter2 =SetBackgroundImageFilterType::New();
 162     setBackgroundFilter2->SetInput(labelImage);
 163     setBackgroundFilter2->SetInput2(labelImage);
 164     setBackgroundFilter2->SetMaskValue(-1);
 165     setBackgroundFilter2->SetOutsideValue(m_ErosionPaddingValue);
 166     if(m_Verbose) std::cout<<"Setting the eroded region to "<<m_ErosionPaddingValue<<"..."<<std::endl;
 167
 168     //---------------------------------
 169     // Cast
 170     //---------------------------------
 171     typename OutputCastImageFilterType::Pointer outputCastImageFilter=OutputCastImageFilterType::New();
 172     outputCastImageFilter->SetInput(setBackgroundFilter2->GetOutput());
 173     if(m_Verbose) std::cout<<"Casting the output..."<<std::endl;
 174     outputCastImageFilter->Update();
 175
 176     //---------------------------------
 177     // SetOutput
 178     //---------------------------------
 179     this->SetNthOutput(0, outputCastImageFilter->GetOutput());
 180
 181
 182   }
 183
 184
 185 }//end clitk
 186
 187 #endif //#define clitkReconstructThroughDilationImageFilter_txx