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