Added filler for small holes taken from http://www.insight-journal.or…

…g/browse/publication/835/

SmallHoleFiller.h

@@ -0,0 +1,47 @@
+#ifndef SmallHoleFiller_H
+#define SmallHoleFiller_H
+
+template<typename TImage>
+class SmallHoleFiller
+{
+public:
+  // Constructor
+  SmallHoleFiller();
+
+  // Inputs
+  void SetImage(typename TImage::Pointer image);
+  void SetHolePixel(typename TImage::PixelType pixel);
+
+  // Outputs
+  typename TImage::Pointer GetOutput();
+
+  // This is the main loop. It simply calls Iterate() until complete.
+  void Fill();
+
+  // This is the core functionality.
+  void Iterate();
+
+  // This function returns true if any of the Output pixels match the HolePixel. This indicates there is more work to be done.
+  bool HasEmptyPixels();
+
+private:
+  // The input image.
+  typename TImage::Pointer Image;
+
+  // The intermediate and eventually output image.
+  typename TImage::Pointer Output;
+
+  // The value of a pixel to be considered a hole (to be filled).
+  typename TImage::PixelType HolePixel;
+
+
+};
+
+// This function copies the data from 'input' to 'output'
+template<typename TImage>
+void DeepCopy(typename TImage::Pointer input, typename TImage::Pointer output);
+
+
+#include "SmallHoleFiller.hxx"
+
+#endif

SmallHoleFiller.hxx

@@ -0,0 +1,157 @@
+#include "itkImageFileWriter.h" // For intermediate debugging output
+#include "itkImageRegionConstIterator.h"
+#include "itkImageRegionIterator.h"
+#include "itkNeighborhoodIterator.h"
+#include "itkNumericTraits.h"
+#include "rtkMacro.h"
+
+template <typename TImage>
+SmallHoleFiller<TImage>::SmallHoleFiller()
+{
+  this->HolePixel = itk::NumericTraits< typename TImage::PixelType >::Zero;
+  this->Image = NULL;
+  this->Output = NULL;
+}
+
+template <typename TImage>
+void SmallHoleFiller<TImage>::SetImage(typename TImage::Pointer image)
+{
+  this->Image = image;
+  this->Output = TImage::New();
+}
+
+template <typename TImage>
+void SmallHoleFiller<TImage>::SetHolePixel(typename TImage::PixelType pixel)
+{
+  this->HolePixel = pixel;
+}
+
+template <typename TImage>
+typename TImage::Pointer SmallHoleFiller<TImage>::GetOutput()
+{
+  return this->Output;
+}
+
+template <typename TImage>
+void SmallHoleFiller<TImage>::Fill()
+{
+  if(!this->Image)
+    {
+    std::cerr << "You must first call SetImage!" << std::endl;
+    return;
+    }
+
+  // Initialize by setting the output image to the input image.
+  DeepCopy<TImage>(this->Image, this->Output);
+  unsigned int numberOfIterations = 0;
+  while(HasEmptyPixels())
+    {
+    std::cout << "Iteration " << numberOfIterations << "..." << std::endl;
+    Iterate();
+    numberOfIterations++;
+
+//    typedef  itk::ImageFileWriter< TImage  > WriterType;
+//    typename WriterType::Pointer writer = WriterType::New();
+//    std::stringstream ss;
+//    ss << "/tmp/intermediate_" << numberOfIterations << ".mha";
+//    writer->SetFileName(ss.str());
+//    writer->SetInput(this->Output);
+//    writer->Update();
+    }
+
+  std::cout << "Filling completed in " << numberOfIterations << " iterations." << std::endl;
+}
+
+template <typename TImage>
+void SmallHoleFiller<TImage>::Iterate()
+{
+  // Make a copy of the current intermediate output. We use this to determine which pixels were holes at the beginning of this iteration.
+  typename TImage::Pointer temporaryImage = TImage::New();
+  DeepCopy<TImage>(this->Output, temporaryImage);
+
+  // Traverse the image. When a pixel is encountered that is a hole, fill it with the average of its non-hole neighbors.
+  // Do not mark pixels filled on this iteration as known. This will result in a less accurate filling, as it favors the colors
+  // of pixels that occur earlier in the raster scan order.
+  typename TImage::SizeType radius;
+  radius.Fill(1);
+
+  itk::NeighborhoodIterator<TImage> intermediateNeighborhoodIterator(radius, temporaryImage, temporaryImage->GetLargestPossibleRegion());
+
+  // This iterator is used to set the output pixels.
+  itk::ImageRegionIterator<TImage> outputIterator(this->Output, this->Output->GetLargestPossibleRegion());
+
+  while(!intermediateNeighborhoodIterator.IsAtEnd())
+    {
+    if(intermediateNeighborhoodIterator.GetCenterPixel() == this->HolePixel) // We want to fill this pixel
+      {
+      typename TImage::PixelType pixelSum = itk::NumericTraits< typename TImage::PixelType >::Zero;
+      // Loop over the 8-neighorhood
+      unsigned int validPixels = 0;
+      for(unsigned int i = 0; i < 27; i++)
+	{
+	if(i==13) // this is the center (current) pixel, so skip it
+	  {
+	  continue;
+	  }
+	bool isInBounds = false;
+	typename TImage::PixelType currentPixel = intermediateNeighborhoodIterator.GetPixel(i, isInBounds);
+	if(isInBounds)
+	  {
+	  if(currentPixel != this->HolePixel)
+	    {
+	    validPixels++;
+	    pixelSum += currentPixel;
+	    }
+	  }
+	} // end 8-connected neighbor for
+
+      if(validPixels > 0)
+	{
+	//typename TImage::PixelType pixelAverage = static_cast<typename TImage::PixelType>(pixelSum / validPixels);
+	typename TImage::PixelType pixelAverage = static_cast<typename TImage::PixelType>(pixelSum * (1.0/ validPixels)); // We multiply by the reciprocal because operator/ is not defined for all types.
+	outputIterator.Set(pixelAverage);
+	//std::cout << "Set " << outputIterator.GetIndex() << " to " << pixelAverage << std::endl;
+	}
+      } // end "fill this pixel?" conditional
+
+    ++intermediateNeighborhoodIterator;
+    ++outputIterator;
+    } // end main traversal loop
+
+}
+
+template <typename TImage>
+bool SmallHoleFiller<TImage>::HasEmptyPixels()
+{
+  itk::ImageRegionConstIterator<TImage> imageIterator(this->Output, this->Output->GetLargestPossibleRegion());
+
+  while(!imageIterator.IsAtEnd())
+    {
+    if(imageIterator.Get() == this->HolePixel)
+      {
+      //std::cout << "Pixel " << imageIterator.GetIndex() << " is empty." << std::endl;
+      return true;
+      }
+    ++imageIterator;
+    }
+  return false;
+}
+
+///////// This is not a member function /////////////
+/** Copy the input to the output*/
+template<typename TImage>
+void DeepCopy(typename TImage::Pointer input, typename TImage::Pointer output)
+{
+  output->SetRegions(input->GetLargestPossibleRegion());
+  output->Allocate();
+
+  itk::ImageRegionConstIterator<TImage> inputIterator(input, input->GetLargestPossibleRegion());
+  itk::ImageRegionIterator<TImage> outputIterator(output, output->GetLargestPossibleRegion());
+
+  while(!inputIterator.IsAtEnd())
+    {
+    outputIterator.Set(inputIterator.Get());
+    ++inputIterator;
+    ++outputIterator;
+    }
+}

pctbinning.cxx

@@ -5,10 +5,12 @@
 #include <rtkConstantImageSource.h>
 
 #include "pctProtonPairsToDistanceDrivenProjection.h"
+#include "SmallHoleFiller.h"
 
 #include <itkImageFileWriter.h>
 #include <itkRegularExpressionSeriesFileNames.h>
 #include <itkTimeProbe.h>
+#include <itkChangeInformationImageFilter.h>
 
 int main(int argc, char * argv[])
 {
@@ -70,11 +72,26 @@ int main(int argc, char * argv[])
 
   TRY_AND_EXIT_ON_ITK_EXCEPTION( projection->Update() );
 
+  SmallHoleFiller< OutputImageType > filler;
+  filler.SetImage( projection->GetOutput() );
+  filler.SetHolePixel(0.);
+  filler.Fill();
+
+  typedef itk::ChangeInformationImageFilter< OutputImageType > CIIType;
+  CIIType::Pointer cii = CIIType::New();
+  cii->SetInput(filler.GetOutput());
+  cii->ChangeOriginOn();
+  cii->ChangeDirectionOn();
+  cii->ChangeSpacingOn();
+  cii->SetOutputDirection( projection->GetOutput()->GetDirection() );
+  cii->SetOutputOrigin(    projection->GetOutput()->GetOrigin() );
+  cii->SetOutputSpacing(   projection->GetOutput()->GetSpacing() );
+
   // Write
   typedef itk::ImageFileWriter<  OutputImageType > WriterType;
   WriterType::Pointer writer = WriterType::New();
   writer->SetFileName( args_info.output_arg );
-  writer->SetInput( projection->GetOutput() );
+  writer->SetInput( cii->GetOutput() );
   TRY_AND_EXIT_ON_ITK_EXCEPTION( writer->Update() );
 
   if(args_info.count_given)