From 554d37d2e44e150788d6acd40b48cfcfc4da0f34 Mon Sep 17 00:00:00 2001 From: Niels Dekker Date: Fri, 16 Feb 2024 10:49:36 +0100 Subject: [PATCH 1/2] STYLE: Remove `UseMask` property from ComputeImageExtremaFilter Assumed that ComputeImageExtremaFilter should use its mask, if and only if the mask is not null. --- .../itkAdvancedImageToImageMetric.hxx | 13 ++----------- .../GTesting/itkComputeImageExtremaFilterGTest.cxx | 4 ---- Common/itkComputeImageExtremaFilter.h | 3 --- Common/itkComputeImageExtremaFilter.hxx | 11 ++++------- .../itkAdvancedMeanSquaresImageToImageMetric.hxx | 14 ++------------ 5 files changed, 8 insertions(+), 37 deletions(-) diff --git a/Common/CostFunctions/itkAdvancedImageToImageMetric.hxx b/Common/CostFunctions/itkAdvancedImageToImageMetric.hxx index eb9190276..6c0736732 100644 --- a/Common/CostFunctions/itkAdvancedImageToImageMetric.hxx +++ b/Common/CostFunctions/itkAdvancedImageToImageMetric.hxx @@ -192,12 +192,7 @@ AdvancedImageToImageMetric::InitializeLimiters() const auto computeFixedImageExtrema = ComputeImageExtremaFilter::New(); computeFixedImageExtrema->SetInput(this->GetFixedImage()); - if (const auto * const fMask = this->GetFixedImageMask()) - { - computeFixedImageExtrema->SetUseMask(true); - computeFixedImageExtrema->SetImageSpatialMask(fMask); - } - + computeFixedImageExtrema->SetImageSpatialMask(this->GetFixedImageMask()); computeFixedImageExtrema->Update(); this->m_FixedImageTrueMax = computeFixedImageExtrema->GetMaximum(); @@ -228,11 +223,7 @@ AdvancedImageToImageMetric::InitializeLimiters() const auto computeMovingImageExtrema = ComputeImageExtremaFilter::New(); computeMovingImageExtrema->SetInput(this->GetMovingImage()); - if (const auto * const mask = this->GetMovingImageMask()) - { - computeMovingImageExtrema->SetUseMask(true); - computeMovingImageExtrema->SetImageSpatialMask(mask); - } + computeMovingImageExtrema->SetImageSpatialMask(this->GetMovingImageMask()); computeMovingImageExtrema->Update(); this->m_MovingImageTrueMax = computeMovingImageExtrema->GetMaximum(); diff --git a/Common/GTesting/itkComputeImageExtremaFilterGTest.cxx b/Common/GTesting/itkComputeImageExtremaFilterGTest.cxx index 30758b085..ae5a481de 100644 --- a/Common/GTesting/itkComputeImageExtremaFilterGTest.cxx +++ b/Common/GTesting/itkComputeImageExtremaFilterGTest.cxx @@ -196,7 +196,6 @@ Expect_one_non_zero_pixel_value_masked_in(const typename TImage::SizeType & imag const auto filter = CheckNew(); filter->SetInput(image); filter->SetImageSpatialMask(maskSpatialObject); - filter->SetUseMask(true); filter->Update(); EXPECT_EQ(filter->GetMinimum(), pixelValue); @@ -247,7 +246,6 @@ Expect_one_positive_pixel_value_all_pixels_masked_in(const typename TImage::Size const auto filter = CheckNew(); filter->SetInput(image); filter->SetImageSpatialMask(maskSpatialObject); - filter->SetUseMask(true); filter->Update(); EXPECT_EQ(filter->GetMinimum(), 0); @@ -356,8 +354,6 @@ GTEST_TEST(ComputeImageExtremaFilter, MaskSpacingAndDirectionAffectResults) elastix::DefaultConstruct filter{}; filter.SetInput(image); filter.SetImageSpatialMask(maskSpatialObject); - filter.SetUseMask(true); - filter.Update(); return filter.GetMaximum(); }; diff --git a/Common/itkComputeImageExtremaFilter.h b/Common/itkComputeImageExtremaFilter.h index 7ca068cf1..1d7c554ef 100644 --- a/Common/itkComputeImageExtremaFilter.h +++ b/Common/itkComputeImageExtremaFilter.h @@ -76,8 +76,6 @@ class ITK_TEMPLATE_EXPORT ComputeImageExtremaFilter : public StatisticsImageFilt /** Type to use for computations. */ using typename Superclass::RealType; - itkSetMacro(UseMask, bool); - using ImageSpatialMaskType = ImageMaskSpatialObject; using ImageSpatialMaskPointer = typename ImageSpatialMaskType::Pointer; using ImageSpatialMaskConstPointer = typename ImageSpatialMaskType::ConstPointer; @@ -105,7 +103,6 @@ class ITK_TEMPLATE_EXPORT ComputeImageExtremaFilter : public StatisticsImageFilt private: ImageSpatialMaskConstPointer m_ImageSpatialMask{}; - bool m_UseMask{ false }; bool m_SameGeometry{ false }; CompensatedSummation m_ThreadSum{ 1 }; diff --git a/Common/itkComputeImageExtremaFilter.hxx b/Common/itkComputeImageExtremaFilter.hxx index 04fdee20f..ec9352ca4 100644 --- a/Common/itkComputeImageExtremaFilter.hxx +++ b/Common/itkComputeImageExtremaFilter.hxx @@ -29,7 +29,7 @@ template void ComputeImageExtremaFilter::BeforeStreamedGenerateData() { - if (!this->m_UseMask) + if (m_ImageSpatialMask == nullptr) { Superclass::BeforeStreamedGenerateData(); } @@ -57,7 +57,7 @@ template void ComputeImageExtremaFilter::AfterStreamedGenerateData() { - if (!this->m_UseMask) + if (m_ImageSpatialMask == nullptr) { Superclass::AfterStreamedGenerateData(); } @@ -89,16 +89,13 @@ template void ComputeImageExtremaFilter::ThreadedStreamedGenerateData(const RegionType & regionForThread) { - if (!this->m_UseMask) + if (m_ImageSpatialMask == nullptr) { Superclass::ThreadedStreamedGenerateData(regionForThread); } else { - if (this->GetImageSpatialMask()) - { - this->ThreadedGenerateDataImageSpatialMask(regionForThread); - } + this->ThreadedGenerateDataImageSpatialMask(regionForThread); } } // end ThreadedGenerateData() diff --git a/Components/Metrics/AdvancedMeanSquares/itkAdvancedMeanSquaresImageToImageMetric.hxx b/Components/Metrics/AdvancedMeanSquares/itkAdvancedMeanSquaresImageToImageMetric.hxx index 2085bfdc7..78f463894 100644 --- a/Components/Metrics/AdvancedMeanSquares/itkAdvancedMeanSquaresImageToImageMetric.hxx +++ b/Components/Metrics/AdvancedMeanSquares/itkAdvancedMeanSquaresImageToImageMetric.hxx @@ -56,12 +56,7 @@ AdvancedMeanSquaresImageToImageMetric::Initialize() /** Try to guess a normalization factor. */ const auto computeFixedImageExtrema = ComputeImageExtremaFilter::New(); computeFixedImageExtrema->SetInput(this->GetFixedImage()); - if (const auto * const mask = this->GetFixedImageMask()) - { - computeFixedImageExtrema->SetUseMask(true); - computeFixedImageExtrema->SetImageSpatialMask(mask); - } - + computeFixedImageExtrema->SetImageSpatialMask(this->GetFixedImageMask()); computeFixedImageExtrema->Update(); this->m_FixedImageTrueMax = computeFixedImageExtrema->GetMaximum(); @@ -76,12 +71,7 @@ AdvancedMeanSquaresImageToImageMetric::Initialize() const auto computeMovingImageExtrema = ComputeImageExtremaFilter::New(); computeMovingImageExtrema->SetInput(this->GetMovingImage()); - if (const auto * const mask = this->GetMovingImageMask()) - { - computeMovingImageExtrema->SetUseMask(true); - computeMovingImageExtrema->SetImageSpatialMask(mask); - } - + computeMovingImageExtrema->SetImageSpatialMask(this->GetMovingImageMask()); computeMovingImageExtrema->Update(); this->m_MovingImageTrueMax = computeMovingImageExtrema->GetMaximum(); From 476969aad7ce57f0e09f4b0fee264e292630e78e Mon Sep 17 00:00:00 2001 From: Niels Dekker Date: Fri, 16 Feb 2024 11:09:36 +0100 Subject: [PATCH 2/2] PERF: Let ComputeImageExtremaFilter only compute min and max No longer the mean, sigma, variance, sum, or sum of squares. Removed m_ThreadSum, m_SumOfSquares, and m_Count. Removed inheritance from StatisticsImageFilter (which computed the mean, sigma, variance, sum, and sum of squares for images without a mask). ComputeImageExtremaFilter is only being used in elastix by AdvancedImageToImageMetric and AdvancedMeanSquaresImageToImageMetric, and both of them only need to have the minimum and the maximum from ComputeImageExtremaFilter. A performance improvement of ~8 percent was observed, from more than 0.13 sec. (before this commit) to less than 0.12 sec., on a 8192x8192 image with mask. Without a mask, the computation appears even more than six times as fast as before, from more than 0.25 sec. (before this commit) to less than 0.04 sec. (after this commit), on a 16384x16384 image. Using Visual Studio 2022 (Release). --- Common/itkComputeImageExtremaFilter.h | 72 +++++----- Common/itkComputeImageExtremaFilter.hxx | 174 +++++++++--------------- 2 files changed, 93 insertions(+), 153 deletions(-) diff --git a/Common/itkComputeImageExtremaFilter.h b/Common/itkComputeImageExtremaFilter.h index 1d7c554ef..f84d7ca58 100644 --- a/Common/itkComputeImageExtremaFilter.h +++ b/Common/itkComputeImageExtremaFilter.h @@ -18,40 +18,26 @@ #ifndef itkComputeImageExtremaFilter_h #define itkComputeImageExtremaFilter_h -#include "itkStatisticsImageFilter.h" +#include "itkImageSink.h" #include "itkImageMaskSpatialObject.h" namespace itk { /** \class ComputeImageExtremaFilter - * \brief Compute min. max, variance and mean of an Image. - * - * StatisticsImageFilter computes the minimum, maximum, sum, mean, variance - * sigma of an image. The filter needs all of its input image. It - * behaves as a filter with an input and output. Thus it can be inserted - * in a pipline with other filters and the statistics will only be - * recomputed if a downstream filter changes. - * - * The filter passes its input through unmodified. The filter is - * threaded. It computes statistics in each thread then combines them in - * its AfterThreadedGenerate method. + * \brief Compute minimum and maximum pixel value of an Image. * * \ingroup MathematicalStatisticsImageFilters * \ingroup ITKImageStatistics - * - * \wiki - * \wikiexample{Statistics/StatisticsImageFilter,Compute min\, max\, variance and mean of an Image.} - * \endwiki */ template -class ITK_TEMPLATE_EXPORT ComputeImageExtremaFilter : public StatisticsImageFilter +class ITK_TEMPLATE_EXPORT ComputeImageExtremaFilter : public ImageSink { public: ITK_DISALLOW_COPY_AND_MOVE(ComputeImageExtremaFilter); /** Standard Self typedef */ using Self = ComputeImageExtremaFilter; - using Superclass = StatisticsImageFilter; + using Superclass = ImageSink; using Pointer = SmartPointer; using ConstPointer = SmartPointer; @@ -64,52 +50,58 @@ class ITK_TEMPLATE_EXPORT ComputeImageExtremaFilter : public StatisticsImageFilt /** Image related typedefs. */ using InputImagePointer = typename TInputImage::Pointer; - using typename Superclass::RegionType; - using typename Superclass::SizeType; - using typename Superclass::IndexType; - using typename Superclass::PixelType; - using PointType = typename TInputImage::PointType; + using Superclass::InputImageDimension; + using typename Superclass::InputImageRegionType; + using PixelType = typename Superclass::InputImagePixelType; /** Image related typedefs. */ itkStaticConstMacro(ImageDimension, unsigned int, TInputImage::ImageDimension); - /** Type to use for computations. */ - using typename Superclass::RealType; - using ImageSpatialMaskType = ImageMaskSpatialObject; using ImageSpatialMaskPointer = typename ImageSpatialMaskType::Pointer; using ImageSpatialMaskConstPointer = typename ImageSpatialMaskType::ConstPointer; itkSetConstObjectMacro(ImageSpatialMask, ImageSpatialMaskType); itkGetConstObjectMacro(ImageSpatialMask, ImageSpatialMaskType); + PixelType + GetMinimum() const + { + return m_ThreadMin; + } + + PixelType + GetMaximum() const + { + return m_ThreadMax; + } + protected: ComputeImageExtremaFilter() = default; ~ComputeImageExtremaFilter() override = default; - /** Initialize some accumulators before the threads run. */ + /** Initialize minimum and maximum before the threads run. */ void BeforeStreamedGenerateData() override; - /** Do final mean and variance computation from data accumulated in threads. - */ - void - AfterStreamedGenerateData() override; - /** Multi-thread version GenerateData. */ void - ThreadedStreamedGenerateData(const RegionType &) override; - virtual void - ThreadedGenerateDataImageSpatialMask(const RegionType &); + ThreadedStreamedGenerateData(const InputImageRegionType &) override; private: + struct MinMaxResult + { + PixelType Min; + PixelType Max; + }; + + static MinMaxResult + RetrieveMinMax(const TInputImage &, const InputImageRegionType &, const ImageSpatialMaskType *, bool); + ImageSpatialMaskConstPointer m_ImageSpatialMask{}; bool m_SameGeometry{ false }; - CompensatedSummation m_ThreadSum{ 1 }; - CompensatedSummation m_SumOfSquares{ 1 }; - SizeValueType m_Count{ 1 }; - PixelType m_ThreadMin{ 1 }; - PixelType m_ThreadMax{ 1 }; + PixelType m_ThreadMin{ 1 }; + PixelType m_ThreadMax{ 1 }; std::mutex m_Mutex{}; }; // end of class diff --git a/Common/itkComputeImageExtremaFilter.hxx b/Common/itkComputeImageExtremaFilter.hxx index ec9352ca4..3199277da 100644 --- a/Common/itkComputeImageExtremaFilter.hxx +++ b/Common/itkComputeImageExtremaFilter.hxx @@ -17,9 +17,12 @@ *=========================================================================*/ #ifndef itkComputeImageExtremaFilter_hxx #define itkComputeImageExtremaFilter_hxx + #include "itkComputeImageExtremaFilter.h" #include +#include + #include "elxMaskHasSameImageDomain.h" namespace itk @@ -29,141 +32,86 @@ template void ComputeImageExtremaFilter::BeforeStreamedGenerateData() { - if (m_ImageSpatialMask == nullptr) - { - Superclass::BeforeStreamedGenerateData(); - } - else - { - // Resize the thread temporaries - m_Count = SizeValueType{}; - m_SumOfSquares = RealType{}; - m_ThreadSum = RealType{}; - m_ThreadMin = NumericTraits::max(); - m_ThreadMax = NumericTraits::NonpositiveMin(); - - if (this->GetImageSpatialMask()) - { - this->m_SameGeometry = elastix::MaskHasSameImageDomain(*m_ImageSpatialMask, *(this->GetInput())); - } - else - { - this->m_SameGeometry = false; - } - } -} - -template -void -ComputeImageExtremaFilter::AfterStreamedGenerateData() -{ - if (m_ImageSpatialMask == nullptr) - { - Superclass::AfterStreamedGenerateData(); - } - else - { - const SizeValueType count = m_Count; - const RealType sumOfSquares(m_SumOfSquares); - const PixelType minimum = m_ThreadMin; - const PixelType maximum = m_ThreadMax; - const RealType sum(m_ThreadSum); - - const RealType mean = sum / static_cast(count); - const RealType variance = - (sumOfSquares - (sum * sum / static_cast(count))) / (static_cast(count) - 1); - const RealType sigma = std::sqrt(variance); - - // Set the outputs - this->SetMinimum(minimum); - this->SetMaximum(maximum); - this->SetMean(mean); - this->SetSigma(sigma); - this->SetVariance(variance); - this->SetSum(sum); - this->SetSumOfSquares(sumOfSquares); - } + m_ThreadMin = NumericTraits::max(); + m_ThreadMax = NumericTraits::NonpositiveMin(); + m_SameGeometry = + (m_ImageSpatialMask != nullptr) && elastix::MaskHasSameImageDomain(*m_ImageSpatialMask, *(this->GetInput())); } -template -void -ComputeImageExtremaFilter::ThreadedStreamedGenerateData(const RegionType & regionForThread) -{ - if (m_ImageSpatialMask == nullptr) - { - Superclass::ThreadedStreamedGenerateData(regionForThread); - } - else - { - this->ThreadedGenerateDataImageSpatialMask(regionForThread); - } -} // end ThreadedGenerateData() template -void -ComputeImageExtremaFilter::ThreadedGenerateDataImageSpatialMask(const RegionType & regionForThread) +auto +ComputeImageExtremaFilter::RetrieveMinMax(const TInputImage & inputImage, + const InputImageRegionType & regionForThread, + const ImageSpatialMaskType * const imageSpatialMask, + const bool sameGeometry) -> MinMaxResult { - if (regionForThread.GetSize(0) == 0) - { - return; - } - RealType sum{}; - RealType sumOfSquares{}; - SizeValueType count{}; - PixelType min = NumericTraits::max(); - PixelType max = NumericTraits::NonpositiveMin(); + PixelType min = NumericTraits::max(); + PixelType max = NumericTraits::NonpositiveMin(); - const auto & inputImage = *(this->GetInput()); - - if (this->m_SameGeometry) + if (imageSpatialMask) { - const auto & maskImage = *(this->m_ImageSpatialMask->GetImage()); + if (sameGeometry) + { + const auto & maskImage = *(imageSpatialMask->GetImage()); - for (ImageRegionConstIterator it(&inputImage, regionForThread); !it.IsAtEnd(); ++it) + for (ImageRegionConstIterator it(&inputImage, regionForThread); !it.IsAtEnd(); ++it) + { + if (maskImage.GetPixel(it.GetIndex()) != PixelType{}) + { + const PixelType value = it.Get(); + min = std::min(min, value); + max = std::max(max, value); + } + } + } + else { - if (maskImage.GetPixel(it.GetIndex()) != PixelType{}) + for (ImageRegionConstIterator it(&inputImage, regionForThread); !it.IsAtEnd(); ++it) { - const PixelType value = it.Get(); - const auto realValue = static_cast(value); - - min = std::min(min, value); - max = std::max(max, value); - - sum += realValue; - sumOfSquares += (realValue * realValue); - ++count; + typename ImageSpatialMaskType::PointType point; + inputImage.TransformIndexToPhysicalPoint(it.GetIndex(), point); + if (imageSpatialMask->IsInsideInWorldSpace(point)) + { + const PixelType value = it.Get(); + min = std::min(min, value); + max = std::max(max, value); + } } - } // end for + } } else { - for (ImageRegionConstIterator it(&inputImage, regionForThread); !it.IsAtEnd(); ++it) + for (ImageScanlineConstIterator it(&inputImage, regionForThread); !it.IsAtEnd(); it.NextLine()) { - PointType point; - inputImage.TransformIndexToPhysicalPoint(it.GetIndex(), point); - if (this->m_ImageSpatialMask->IsInsideInWorldSpace(point)) + while (!it.IsAtEndOfLine()) { const PixelType value = it.Get(); - const auto realValue = static_cast(value); - min = std::min(min, value); max = std::max(max, value); - - sum += realValue; - sumOfSquares += (realValue * realValue); - ++count; + ++it; } - } // end for - } // end if + } + } + return { min, max }; +} - const std::lock_guard lockGuard(m_Mutex); - m_ThreadSum += sum; - m_SumOfSquares += sumOfSquares; - m_Count += count; - m_ThreadMin = std::min(min, m_ThreadMin); - m_ThreadMax = std::max(max, m_ThreadMax); -} // end ThreadedGenerateDataImageSpatialMask() +template +void +ComputeImageExtremaFilter::ThreadedStreamedGenerateData(const InputImageRegionType & regionForThread) +{ + if (regionForThread.GetSize(0) > 0) + { + const MinMaxResult minMaxResult = + RetrieveMinMax(*(this->GetInput()), regionForThread, m_ImageSpatialMask, m_SameGeometry); + + // Lock after calling RetrieveMinMax. + const std::lock_guard lockGuard(m_Mutex); + m_ThreadMin = std::min(minMaxResult.Min, m_ThreadMin); + m_ThreadMax = std::max(minMaxResult.Max, m_ThreadMax); + } +} } // end namespace itk #endif