From 60dcdb0f627009029db67dfb6470f1fae7c9b51c Mon Sep 17 00:00:00 2001 From: Axel Garcia Date: Thu, 23 Oct 2025 17:14:22 +0200 Subject: [PATCH 1/3] STYLE: Rename SmallHoleFiller to pctHoleFillingImageFilter --- include/{SmallHoleFiller.h => pctHoleFillingImageFilter.h} | 0 include/{SmallHoleFiller.hxx => pctHoleFillingImageFilter.hxx} | 0 .../{pctSmallHoleFiller.wrap => pctHoleFillingImageFilter.wrap} | 0 3 files changed, 0 insertions(+), 0 deletions(-) rename include/{SmallHoleFiller.h => pctHoleFillingImageFilter.h} (100%) rename include/{SmallHoleFiller.hxx => pctHoleFillingImageFilter.hxx} (100%) rename wrapping/{pctSmallHoleFiller.wrap => pctHoleFillingImageFilter.wrap} (100%) diff --git a/include/SmallHoleFiller.h b/include/pctHoleFillingImageFilter.h similarity index 100% rename from include/SmallHoleFiller.h rename to include/pctHoleFillingImageFilter.h diff --git a/include/SmallHoleFiller.hxx b/include/pctHoleFillingImageFilter.hxx similarity index 100% rename from include/SmallHoleFiller.hxx rename to include/pctHoleFillingImageFilter.hxx diff --git a/wrapping/pctSmallHoleFiller.wrap b/wrapping/pctHoleFillingImageFilter.wrap similarity index 100% rename from wrapping/pctSmallHoleFiller.wrap rename to wrapping/pctHoleFillingImageFilter.wrap From fc45319b513b73f556132b4e4b0d317e49882f7c Mon Sep 17 00:00:00 2001 From: Axel Garcia Date: Thu, 23 Oct 2025 17:19:20 +0200 Subject: [PATCH 2/3] ENH: Adapt pctHoleFillingImageFilter to inherits from itkImageToImageFilter --- .../pctbackprojectionbinning.cxx | 11 +- applications/pctbinning/pctbinning.cxx | 11 +- include/pctHoleFillingImageFilter.h | 83 ++++--- include/pctHoleFillingImageFilter.hxx | 215 ++++++------------ wrapping/pctHoleFillingImageFilter.wrap | 2 +- 5 files changed, 124 insertions(+), 198 deletions(-) diff --git a/applications/pctbackprojectionbinning/pctbackprojectionbinning.cxx b/applications/pctbackprojectionbinning/pctbackprojectionbinning.cxx index 855db23..707a6b1 100644 --- a/applications/pctbackprojectionbinning/pctbackprojectionbinning.cxx +++ b/applications/pctbackprojectionbinning/pctbackprojectionbinning.cxx @@ -6,7 +6,7 @@ #include #include "pctProtonPairsToBackProjection.h" -#include "SmallHoleFiller.h" +#include "pctHoleFillingImageFilter.h" #include #include @@ -124,17 +124,16 @@ main(int argc, char * argv[]) TRY_AND_EXIT_ON_ITK_EXCEPTION(projection->Update()); - SmallHoleFiller filler; + auto holeFilter = pct::HoleFillingImageFilter::New(); if (args_info.fill_flag) { - filler.SetImage(projection->GetOutput()); - filler.SetHolePixel(0.); - filler.Fill(); + holeFilter->SetInput(projection->GetOutput()); + TRY_AND_EXIT_ON_ITK_EXCEPTION(holeFilter->Update()); } auto cii = itk::ChangeInformationImageFilter::New(); if (args_info.fill_flag) - cii->SetInput(filler.GetOutput()); + cii->SetInput(holeFilter->GetOutput()); else cii->SetInput(projection->GetOutput()); cii->ChangeOriginOn(); diff --git a/applications/pctbinning/pctbinning.cxx b/applications/pctbinning/pctbinning.cxx index 5f51925..5921b6f 100644 --- a/applications/pctbinning/pctbinning.cxx +++ b/applications/pctbinning/pctbinning.cxx @@ -5,7 +5,7 @@ #include #include "pctProtonPairsToDistanceDrivenProjection.h" -#include "SmallHoleFiller.h" +#include "pctHoleFillingImageFilter.h" #include #include @@ -87,17 +87,16 @@ main(int argc, char * argv[]) TRY_AND_EXIT_ON_ITK_EXCEPTION(projection->Update()); - SmallHoleFiller filler; + auto holeFilter = pct::HoleFillingImageFilter::New(); if (args_info.fill_flag) { - filler.SetImage(projection->GetOutput()); - filler.SetHolePixel(0.); - filler.Fill(); + holeFilter->SetInput(projection->GetOutput()); + TRY_AND_EXIT_ON_ITK_EXCEPTION(holeFilter->Update()); } auto cii = itk::ChangeInformationImageFilter::New(); if (args_info.fill_flag) - cii->SetInput(filler.GetOutput()); + cii->SetInput(holeFilter->GetOutput()); else cii->SetInput(projection->GetOutput()); cii->ChangeOriginOn(); diff --git a/include/pctHoleFillingImageFilter.h b/include/pctHoleFillingImageFilter.h index 6a4e35f..5f4cf4a 100644 --- a/include/pctHoleFillingImageFilter.h +++ b/include/pctHoleFillingImageFilter.h @@ -1,54 +1,63 @@ -#ifndef __SmallHoleFiller_h -#define __SmallHoleFiller_h +#ifndef __pctHoleFillingImageFilter_h +#define __pctHoleFillingImageFilter_h -template -class SmallHoleFiller +#include +#include +#include + +namespace pct +{ +/** + * + * A simple iterative hole filling filter. A pixel equal to the HolePixel is + * considered a hole. On each iteration, every hole pixel is replaced by the + * average of its non-hole neighbors in a radius=1 neighborhood (center + * excluded). Iterations repeat until no hole pixels remain or MaxIterations + * is reached or an iteration makes no progress. + * + * This filter reproduces the behavior of the project's SmallHoleFiller + * helper. For the original reference see: + * https://www.insight-journal.org/browse/publication/835/ + */ + +template +class HoleFillingImageFilter : public itk::ImageToImageFilter { public: - // Constructor - SmallHoleFiller(); + using Self = HoleFillingImageFilter; + using Superclass = itk::ImageToImageFilter; + using Pointer = itk::SmartPointer; + using ConstPointer = itk::SmartPointer; - // Inputs - void - SetImage(const TImage * image); - void - SetHolePixel(typename TImage::PixelType pixel); + itkNewMacro(Self); + itkTypeMacro(HoleFillingImageFilter, ImageToImageFilter); - // Outputs - TImage * - GetOutput(); + using InputImageType = TInputImage; + using OutputImageType = TOutputImage; + using PixelType = typename OutputImageType::PixelType; - // This is the main loop. It simply calls Iterate() until complete. - void - Fill(); + itkSetMacro(HolePixelValue, PixelType); + itkGetConstMacro(HolePixelValue, PixelType); + itkSetMacro(MaximumNumberOfIterations, unsigned int); + itkGetConstMacro(MaximumNumberOfIterations, unsigned int); - // This is the core functionality. - void - Iterate(); +protected: + HoleFillingImageFilter() = default; + ~HoleFillingImageFilter() override = default; - // This function returns true if any of the Output pixels match the HolePixel. This indicates there is more work to be - // done. - bool - HasEmptyPixels(); + void + GenerateData() override; private: - // The input image. - typename TImage::ConstPointer 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; + ITK_DISALLOW_COPY_AND_MOVE(HoleFillingImageFilter); + PixelType m_HolePixelValue = itk::NumericTraits::Zero; + unsigned int m_MaximumNumberOfIterations = 20; }; -// This function copies the data from 'input' to 'output' -template -void -DeepCopy(const TImage * input, TImage * output); +} // namespace pct #ifndef ITK_MANUAL_INSTANTIATION -# include "SmallHoleFiller.hxx" +# include "pctHoleFillingImageFilter.hxx" #endif #endif diff --git a/include/pctHoleFillingImageFilter.hxx b/include/pctHoleFillingImageFilter.hxx index 6fdfc99..a685d60 100644 --- a/include/pctHoleFillingImageFilter.hxx +++ b/include/pctHoleFillingImageFilter.hxx @@ -1,172 +1,91 @@ -#ifndef __SmallHoleFiller_hxx -#define __SmallHoleFiller_hxx +#include "pctHoleFillingImageFilter.h" -#include "itkImageFileWriter.h" // For intermediate debugging output -#include "itkImageRegionConstIterator.h" -#include "itkImageRegionIterator.h" -#include "itkNeighborhoodIterator.h" -#include "itkNumericTraits.h" -#include "rtkMacro.h" +#include +#include +#include +#include +#include +#include -template -SmallHoleFiller::SmallHoleFiller() +namespace pct { - this->HolePixel = itk::NumericTraits::Zero; - this->Image = NULL; - this->Output = NULL; -} -template +template void -SmallHoleFiller::SetImage(const TImage * image) +HoleFillingImageFilter::GenerateData() { - this->Image = image; - this->Output = TImage::New(); -} - -template -void -SmallHoleFiller::SetHolePixel(typename TImage::PixelType pixel) -{ - this->HolePixel = pixel; -} - -template -TImage * -SmallHoleFiller::GetOutput() -{ - return this->Output.GetPointer(); -} - -template -void -SmallHoleFiller::Fill() -{ - if (!this->Image) - { - std::cerr << "You must first call SetImage!" << std::endl; - return; - } + using InputImageType = TInputImage; + using OutputImageType = TOutputImage; + using PixelType = typename OutputImageType::PixelType; + + // Allocate output as a copy of input + auto initialDuplicator = itk::ImageDuplicator::New(); + initialDuplicator->SetInputImage(this->GetInput()); + initialDuplicator->Update(); + this->GetOutput()->Graft(initialDuplicator->GetOutput()); + + unsigned int iteration = 0; + typename OutputImageType::SizeType radius; + radius.Fill(1); - // Initialize by setting the output image to the input image. - DeepCopy(this->Image, this->Output); - unsigned int numberOfIterations = 0; - while (HasEmptyPixels()) + while (iteration < this->m_MaximumNumberOfIterations) { - 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); + // Fast presence check: if there are no hole pixels left, exit early. + bool hasHole = false; + itk::ImageRegionConstIterator checkIt(this->GetOutput(), + this->GetOutput()->GetLargestPossibleRegion()); + for (checkIt.GoToBegin(); !checkIt.IsAtEnd(); ++checkIt) + { + if (checkIt.Get() == this->m_HolePixelValue) + { + hasHole = true; + break; + } + } + if (!hasHole) + break; - // 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); + auto duplicator = itk::ImageDuplicator::New(); + duplicator->SetInputImage(this->GetOutput()); + duplicator->Update(); + typename OutputImageType::Pointer snapshot = duplicator->GetOutput(); - itk::NeighborhoodIterator intermediateNeighborhoodIterator( - radius, temporaryImage, temporaryImage->GetLargestPossibleRegion()); + itk::NeighborhoodIterator neighbor(radius, snapshot, snapshot->GetLargestPossibleRegion()); + itk::ImageRegionIterator writer(this->GetOutput(), this->GetOutput()->GetLargestPossibleRegion()); - // This iterator is used to set the output pixels. - itk::ImageRegionIterator outputIterator(this->Output, this->Output->GetLargestPossibleRegion()); + const unsigned int neighborSize = neighbor.Size(); + const unsigned int centerIndex = neighborSize / 2; - while (!intermediateNeighborhoodIterator.IsAtEnd()) - { - if (intermediateNeighborhoodIterator.GetCenterPixel() == this->HolePixel) // We want to fill this pixel + for (neighbor.GoToBegin(), writer.GoToBegin(); !neighbor.IsAtEnd(); ++neighbor, ++writer) { - typename TImage::PixelType pixelSum = itk::NumericTraits::Zero; - // Loop over the 8-neighorhood - unsigned int validPixels = 0; - for (unsigned int i = 0; i < 27; i++) + PixelType centerPixel = neighbor.GetCenterPixel(); + if (centerPixel == this->m_HolePixelValue) { - 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) + PixelType sum = itk::NumericTraits::Zero; + unsigned int valid = 0; + + for (unsigned int i = 0; i < neighborSize; ++i) { - if (currentPixel != this->HolePixel) + if (i == centerIndex) + continue; + bool inBounds = false; + PixelType p = neighbor.GetPixel(i, inBounds); + if (inBounds && p != this->m_HolePixelValue) { - validPixels++; - pixelSum += currentPixel; + sum += p; + ++valid; } } - } // 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; + if (valid > 0) + { + PixelType outValue = sum / valid; + writer.Set(outValue); + } } - } // 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(const TImage * input, TImage * 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; + ++iteration; } } -#endif +} // namespace pct diff --git a/wrapping/pctHoleFillingImageFilter.wrap b/wrapping/pctHoleFillingImageFilter.wrap index e6a4e40..cccc899 100644 --- a/wrapping/pctHoleFillingImageFilter.wrap +++ b/wrapping/pctHoleFillingImageFilter.wrap @@ -1,4 +1,4 @@ -itk_wrap_class("SmallHoleFiller") +itk_wrap_class("pct::HoleFillingImageFilter" POINTER) foreach(t ${WRAP_ITK_REAL}) foreach(d 3 4) itk_wrap_template("I${ITKM_${t}}${d}" "itk::Image<${ITKT_${t}}, ${d}>") From fda079266c3f4b01c3b9976e69fe577cb6461b97 Mon Sep 17 00:00:00 2001 From: Axel Garcia Date: Wed, 22 Oct 2025 13:47:10 +0200 Subject: [PATCH 3/3] ENH: Add test for pct::HoleFillingImageFilter --- test/CMakeLists.txt | 8 +++- test/pctHoleFillingImageFilterTest.cxx | 51 ++++++++++++++++++++++++++ 2 files changed, 58 insertions(+), 1 deletion(-) create mode 100644 test/pctHoleFillingImageFilterTest.cxx diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index f9ffca5..18ee405 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -3,7 +3,11 @@ itk_module_test() #----------------------------------------------------------------------------- # CXX tests -set(PCTTests pctProtonPairsToDistanceDrivenProjectionTest.cxx) +set( + PCTTests + pctProtonPairsToDistanceDrivenProjectionTest.cxx + pctHoleFillingImageFilterTest.cxx +) createtestdriver(PCT "${PCT-Test_LIBRARIES}" "${PCTTests}") @@ -11,6 +15,8 @@ itk_add_test(NAME pctProtonPairsToDistanceDrivenProjectionTest COMMAND PCTTestDriver pctProtonPairsToDistanceDrivenProjectionTest ) +itk_add_test(NAME pctHoleFillingImageFilterTest COMMAND PCTTestDriver pctHoleFillingImageFilterTest) + #----------------------------------------------------------------------------- # Python tests if(ITK_WRAP_PYTHON) diff --git a/test/pctHoleFillingImageFilterTest.cxx b/test/pctHoleFillingImageFilterTest.cxx new file mode 100644 index 0000000..25d706f --- /dev/null +++ b/test/pctHoleFillingImageFilterTest.cxx @@ -0,0 +1,51 @@ +#include "pctHoleFillingImageFilter.h" + +#include "itkImage.h" +#include "itkImageRegionIterator.h" +#include "itkTestingMacros.h" + +int +pctHoleFillingImageFilterTest(int argc, char * argv[]) +{ + using ImageType = itk::Image; + + auto img = ImageType::New(); + ImageType::RegionType region; + region.SetIndex(itk::MakeIndex(0, 0, 0)); + region.SetSize(itk::MakeSize(8, 8, 8)); + + img->SetRegions(region); + img->Allocate(); + img->FillBuffer(42); + + // Insert some holes (use 0 as hole value) + itk::ImageRegionIterator it(img, region); + for (it.GoToBegin(); !it.IsAtEnd(); ++it) + { + ImageType::IndexType idx = it.GetIndex(); + if ((idx[0] + idx[1] + idx[2]) % 10 == 0) + { + it.Set(0); + } + } + + using FilterType = pct::HoleFillingImageFilter; + FilterType::Pointer filter = FilterType::New(); + filter->SetInput(img); + filter->SetHolePixelValue(0); + filter->SetMaximumNumberOfIterations(10); + filter->Update(); + + ImageType::Pointer out = filter->GetOutput(); + + // Check no zeros remain and values are 42 + itk::ImageRegionIterator outIt(out, region); + for (outIt.GoToBegin(); !outIt.IsAtEnd(); ++outIt) + { + ITK_TEST_EXPECT_EQUAL(outIt.Get(), 42); + } + + std::cout << "\n\nTest PASSED! " << std::endl; + + return EXIT_SUCCESS; +}