reconstruct with dilatation (jef)

diff --git a/itk/clitkReconstructThroughDilationImageFilter.h b/itk/clitkReconstructThroughDilationImageFilter.h
new file mode 100644 (file)
index 0000000..23cc98b
--- /dev/null
+++ b/itk/clitkReconstructThroughDilationImageFilter.h
@@ -0,0 +1,124 @@
+#ifndef clitkReconstructThroughDilationImageFilter_h
+#define clitkReconstructThroughDilationImageFilter_h
+
+/* =================================================
+ * @file   clitkReconstructThroughDilationImageFilter.h
+ * @author 
+ * @date   
+ * 
+ * @brief 
+ * 
+ ===================================================*/
+
+
+// clitk include
+#include "clitkIO.h"
+#include "clitkCommon.h"
+#include "clitkConditionalBinaryDilateImageFilter.h"
+
+//itk include
+#include "itkImageToImageFilter.h"
+#include "itkBinaryBallStructuringElement.h"
+#include "itkConnectedComponentImageFilter.h"
+#include "itkStatisticsImageFilter.h"
+#include "itkCastImageFilter.h"
+#include "itkDifferenceImageFilter.h"
+#include "itkThresholdImageFilter.h"
+
+namespace clitk 
+{
+
+  template <class InputImageType, class OutputImageType>
+  class ITK_EXPORT ReconstructThroughDilationImageFilter :
+    public itk::ImageToImageFilter<InputImageType, OutputImageType>
+  {
+  public:
+    //----------------------------------------
+    // ITK
+    //----------------------------------------
+    typedef ReconstructThroughDilationImageFilter                                                 Self;
+    typedef itk::ImageToImageFilter<InputImageType, OutputImageType>  Superclass;
+    typedef itk::SmartPointer<Self>                                   Pointer;
+    typedef itk::SmartPointer<const Self>                             ConstPointer;
+   
+    // Method for creation through the object factory
+    itkNewMacro(Self);  
+
+    // Run-time type information (and related methods)
+    itkTypeMacro( ReconstructThroughDilationImageFilter, ImageToImageFilter );
+
+    /** Dimension of the domain space. */
+    itkStaticConstMacro(InputImageDimension, unsigned int, Superclass::InputImageDimension);
+    itkStaticConstMacro(OutputImageDimension, unsigned int, Superclass::OutputImageDimension);
+
+    //----------------------------------------
+    // Typedefs
+    //----------------------------------------
+    typedef typename OutputImageType::RegionType OutputImageRegionType;
+    typedef int InternalPixelType;
+    typedef typename InputImageType::PixelType InputPixelType;
+    typedef typename OutputImageType::PixelType OutputPixelType;
+    typedef typename InputImageType::SizeType SizeType;
+
+    //----------------------------------------
+    // Set & Get
+    //----------------------------------------    
+    itkBooleanMacro(Verbose);
+    itkSetMacro( Verbose, bool);
+    itkGetConstReferenceMacro( Verbose, bool);
+    void SetRadius ( const SizeType& s){ m_Radius=s; this->Modified();}
+    SizeType GetRadius(void){return m_Radius;}
+    itkSetMacro( ErosionPaddingValue, OutputPixelType);
+    itkGetConstMacro( ErosionPaddingValue, OutputPixelType)
+    itkSetMacro( MaximumNumberOfLabels, unsigned int);
+    itkGetConstMacro( MaximumNumberOfLabels, unsigned int);
+    itkSetMacro( BackgroundValue, InternalPixelType);
+    itkGetConstMacro( BackgroundValue, InternalPixelType);
+    itkSetMacro( ForegroundValue, InternalPixelType);
+    itkGetConstMacro( ForegroundValue, InternalPixelType);
+
+  protected:
+
+    //----------------------------------------  
+    // Constructor & Destructor
+    //----------------------------------------  
+    ReconstructThroughDilationImageFilter();
+    ~ReconstructThroughDilationImageFilter() {};
+
+    //----------------------------------------  
+    // Update
+    //----------------------------------------  
+    // Generate Data
+    void GenerateData(void);
+
+    //     // Threaded Generate Data
+    //     void BeforeThreadedGenerateData(void );
+    //     void ThreadedGenerateData(const OutputImageRegionType & outputRegionForThread, int threadId );
+    //     void AfterThreadedGenerateData(void );
+    //     // Override defaults
+    //     virtual void GenerateInputRequestedRegion();
+    //     virtual void GenerateOutputInformation (void);
+    //     virtual void EnlargeOutputRequestedRegion(DataObject *data);
+    //     void AllocateOutputs();
+    //----------------------------------------  
+    // Data members
+    //----------------------------------------
+    bool m_Verbose;
+    InternalPixelType m_BackgroundValue;
+    InternalPixelType m_ForegroundValue;
+    unsigned int m_MaximumNumberOfLabels;
+    OutputPixelType m_ErosionPaddingValue;
+    SizeType m_Radius;
+
+  };
+
+
+} // end namespace clitk
+
+#ifndef ITK_MANUAL_INSTANTIATION
+#include "clitkReconstructThroughDilationImageFilter.txx"
+#endif
+
+#endif // #define clitkReconstructThroughDilationImageFilter_h
+
+

diff --git a/itk/clitkReconstructThroughDilationImageFilter.txx b/itk/clitkReconstructThroughDilationImageFilter.txx
new file mode 100644 (file)
index 0000000..d2d3d31
--- /dev/null
+++ b/itk/clitkReconstructThroughDilationImageFilter.txx
@@ -0,0 +1,166 @@
+#ifndef clitkReconstructThroughDilationImageFilter_txx
+#define clitkReconstructThroughDilationImageFilter_txx
+
+/* =================================================
+ * @file   clitkReconstructThroughDilationImageFilter.txx
+ * @author 
+ * @date   
+ * 
+ * @brief 
+ * 
+ ===================================================*/
+
+
+namespace clitk
+{
+
+  //-------------------------------------------------------------------
+  // Update with the number of dimensions
+  //-------------------------------------------------------------------
+  template<class InputImageType, class OutputImageType>
+  ReconstructThroughDilationImageFilter<InputImageType, OutputImageType>::ReconstructThroughDilationImageFilter()
+  {
+   m_Verbose=false;
+   m_BackgroundValue=0;
+   m_ForegroundValue=1;
+   m_ErosionPaddingValue=static_cast<InputPixelType>(-1);
+   for (unsigned int i=0; i<InputImageDimension; i++)
+     m_Radius[i]=1;
+   m_MaximumNumberOfLabels=10;
+  }
+
+
+  //-------------------------------------------------------------------
+  // Update with the number of dimensions and the pixeltype
+  //-------------------------------------------------------------------
+  template<class InputImageType, class  OutputImageType> 
+  void 
+  ReconstructThroughDilationImageFilter<InputImageType, OutputImageType>::GenerateData()
+  {
+
+    //---------------------------------
+    // Typedefs 
+    //--------------------------------- 
+    
+    // Internal type
+    typedef itk::Image<InternalPixelType, InputImageDimension> InternalImageType;
+
+    // Filters used
+    typedef itk::CastImageFilter<InputImageType, InternalImageType> InputCastImageFilterType;   
+    typedef itk::ThresholdImageFilter<InternalImageType> InputThresholdImageFilterType;
+    typedef itk::StatisticsImageFilter<InternalImageType> StatisticsImageFilterType;
+    typedef itk::BinaryBallStructuringElement<InternalPixelType,InputImageDimension > KernelType;
+    typedef clitk::ConditionalBinaryDilateImageFilter<InternalImageType, InternalImageType , KernelType> ConditionalBinaryDilateImageFilterType;
+    typedef itk::DifferenceImageFilter<InternalImageType, InternalImageType> DifferenceImageFilterType;
+    typedef itk::CastImageFilter<InternalImageType, OutputImageType> OutputCastImageFilterType;
+    typedef clitk::SetBackgroundImageFilter<InternalImageType, InternalImageType, InternalImageType> SetBackgroundImageFilterType;
+
+    //---------------------------------
+    // Cast
+    //---------------------------------
+    typename InputCastImageFilterType::Pointer castImageFilter=InputCastImageFilterType::New();
+    castImageFilter->SetInput(this->GetInput());
+    castImageFilter->Update();
+
+    //---------------------------------
+    // Threshold
+    //---------------------------------
+    typename InputThresholdImageFilterType::Pointer thresholdImageFilter=InputThresholdImageFilterType::New();
+    thresholdImageFilter->SetInput(castImageFilter->GetOutput());
+    thresholdImageFilter->ThresholdAbove(m_MaximumNumberOfLabels);
+    thresholdImageFilter->SetOutsideValue(m_ForegroundValue);
+    if(m_Verbose) std::cout<<"Thresholding the input to "<<m_MaximumNumberOfLabels<<" labels ..."<<std::endl;
+    thresholdImageFilter->Update();
+
+    //---------------------------------
+    // Set -1 to padding value
+    //---------------------------------
+    typename SetBackgroundImageFilterType::Pointer setBackgroundFilter =SetBackgroundImageFilterType::New();
+    setBackgroundFilter->SetInput(thresholdImageFilter->GetOutput());
+    setBackgroundFilter->SetInput2(castImageFilter->GetOutput());
+    setBackgroundFilter->SetMaskValue(m_ErosionPaddingValue);
+    setBackgroundFilter->SetOutsideValue(-1);
+    if(m_Verbose) std::cout<<"Setting the eroded region from "<<m_ErosionPaddingValue<<" to -1..."<<std::endl;
+
+    
+    //---------------------------------
+    // Count the initial labels
+    //---------------------------------
+    typename StatisticsImageFilterType::Pointer inputStatisticsImageFilter=StatisticsImageFilterType::New();
+    inputStatisticsImageFilter->SetInput(setBackgroundFilter->GetOutput());
+    if(m_Verbose) std::cout<<"Counting the initial labels..."<<std::endl;
+    inputStatisticsImageFilter->Update();
+    unsigned int initialNumberOfLabels= inputStatisticsImageFilter->GetMaximum();
+    if(m_Verbose) std::cout<<"The input contained "<<initialNumberOfLabels<<" disctictive label(s)..."<<std::endl;
+    unsigned int numberOfConsideredLabels=std::min(initialNumberOfLabels, m_MaximumNumberOfLabels);
+    if(m_Verbose) std::cout<<"Performing dilation the first "<<numberOfConsideredLabels<<" disctictive labels..."<<std::endl;
+
+    //---------------------------------
+    // Dilate while change
+    //---------------------------------
+    typename itk::NumericTraits<InputPixelType>::AccumulateType difference=1;
+    typename InternalImageType::Pointer labelImage=inputStatisticsImageFilter->GetOutput();
+    typename InternalImageType::Pointer oldLabelImage=inputStatisticsImageFilter->GetOutput();
+
+    // element
+    KernelType structuringElement;
+    structuringElement.SetRadius(m_Radius);
+    structuringElement.CreateStructuringElement();
+
+    while( difference)
+      {
+       // Dilate all labels once
+       for ( int label=0; label<(int)numberOfConsideredLabels+1;label++)  
+         if ( m_BackgroundValue != label)
+           {
+             typename ConditionalBinaryDilateImageFilterType::Pointer dilateFilter=ConditionalBinaryDilateImageFilterType::New();
+             dilateFilter->SetBoundaryToForeground(false);
+             dilateFilter->SetKernel(structuringElement);
+             dilateFilter->SetBackgroundValue (-1);
+             dilateFilter->SetInput (labelImage);
+             dilateFilter->SetForegroundValue (label);
+             if(m_Verbose) std::cout<<"Dilating the label "<<label<<"..."<<std::endl;
+             dilateFilter->Update();
+             labelImage=dilateFilter->GetOutput();
+           }
+  
+       // Difference with previous labelImage
+       typename DifferenceImageFilterType::Pointer differenceFilter=DifferenceImageFilterType::New();
+       differenceFilter->SetValidInput(oldLabelImage);
+       differenceFilter->SetTestInput(labelImage);
+       differenceFilter->Update();
+       difference =differenceFilter->GetTotalDifference();
+       if(m_Verbose) std::cout<<"The change in this iteration was "<<difference<<"..."<<std::endl;     
+       oldLabelImage=labelImage;
+      }
+      
+    //---------------------------------
+    // Set -1 to padding value
+    //---------------------------------
+    typename SetBackgroundImageFilterType::Pointer setBackgroundFilter2 =SetBackgroundImageFilterType::New();
+    setBackgroundFilter2->SetInput(labelImage);
+    setBackgroundFilter2->SetInput2(labelImage);
+    setBackgroundFilter2->SetMaskValue(-1);
+    setBackgroundFilter2->SetOutsideValue(m_ErosionPaddingValue);
+    if(m_Verbose) std::cout<<"Setting the eroded region to "<<m_ErosionPaddingValue<<"..."<<std::endl;
+
+    //---------------------------------
+    // Cast
+    //---------------------------------
+    typename OutputCastImageFilterType::Pointer outputCastImageFilter=OutputCastImageFilterType::New();
+    outputCastImageFilter->SetInput(setBackgroundFilter2->GetOutput());
+    if(m_Verbose) std::cout<<"Casting the output..."<<std::endl;
+    outputCastImageFilter->Update();
+
+    //---------------------------------
+    // SetOutput
+    //---------------------------------
+    this->SetNthOutput(0, outputCastImageFilter->GetOutput());
+
+
+  }
+
+
+}//end clitk
+ 
+#endif //#define clitkReconstructThroughDilationImageFilter_txx