diff --git a/SmallHoleFiller.h b/SmallHoleFiller.h new file mode 100644 index 0000000..2b92c7f --- /dev/null +++ b/SmallHoleFiller.h @@ -0,0 +1,47 @@ +#ifndef SmallHoleFiller_H +#define SmallHoleFiller_H + +template +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 +void DeepCopy(typename TImage::Pointer input, typename TImage::Pointer output); + + +#include "SmallHoleFiller.hxx" + +#endif diff --git a/SmallHoleFiller.hxx b/SmallHoleFiller.hxx new file mode 100644 index 0000000..562568c --- /dev/null +++ b/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 +SmallHoleFiller::SmallHoleFiller() +{ + this->HolePixel = itk::NumericTraits< typename TImage::PixelType >::Zero; + this->Image = NULL; + this->Output = NULL; +} + +template +void SmallHoleFiller::SetImage(typename TImage::Pointer image) +{ + this->Image = image; + this->Output = TImage::New(); +} + +template +void SmallHoleFiller::SetHolePixel(typename TImage::PixelType pixel) +{ + this->HolePixel = pixel; +} + +template +typename TImage::Pointer SmallHoleFiller::GetOutput() +{ + return this->Output; +} + +template +void SmallHoleFiller::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(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 +void SmallHoleFiller::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(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 intermediateNeighborhoodIterator(radius, temporaryImage, temporaryImage->GetLargestPossibleRegion()); + + // This iterator is used to set the output pixels. + itk::ImageRegionIterator 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(pixelSum / validPixels); + typename TImage::PixelType pixelAverage = static_cast(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 +bool SmallHoleFiller::HasEmptyPixels() +{ + itk::ImageRegionConstIterator 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 +void DeepCopy(typename TImage::Pointer input, typename TImage::Pointer output) +{ + output->SetRegions(input->GetLargestPossibleRegion()); + output->Allocate(); + + itk::ImageRegionConstIterator inputIterator(input, input->GetLargestPossibleRegion()); + itk::ImageRegionIterator outputIterator(output, output->GetLargestPossibleRegion()); + + while(!inputIterator.IsAtEnd()) + { + outputIterator.Set(inputIterator.Get()); + ++inputIterator; + ++outputIterator; + } +} diff --git a/pctbinning.cxx b/pctbinning.cxx index 1adf723..8447fd6 100644 --- a/pctbinning.cxx +++ b/pctbinning.cxx @@ -5,10 +5,12 @@ #include #include "pctProtonPairsToDistanceDrivenProjection.h" +#include "SmallHoleFiller.h" #include #include #include +#include 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)