From 82984859f1e6e1f9a3d53940c3b540e83c263856 Mon Sep 17 00:00:00 2001 From: Niels Dekker Date: Tue, 9 May 2023 17:12:03 +0200 Subject: [PATCH] PERF: Make EvaluateParzenValues calls faster, using raw buffer of values Added a `ParzenWindowHistogramImageToImageMetric::EvaluateParzenValues` overload that has a raw pointer to the buffer of values as parameter. Stored the local Parzen values of both fixed and moving image together in one array, instead of having separate arrays for fixed and moving images, in both ParzenWindowMutualInformationImageToImageMetric::UpdateDerivativeLowMemory and ParzenWindowHistogramImageToImageMetric::UpdateJointPDFAndDerivatives. A performance improvement of ~5% was observed, executing an ElastixRegistrationMethod on VS2019 Release, using the default registration parameter maps. --- ...kParzenWindowHistogramImageToImageMetric.h | 7 +++ ...arzenWindowHistogramImageToImageMetric.hxx | 43 +++++++++++++------ ...dowMutualInformationImageToImageMetric.hxx | 25 +++++++---- 3 files changed, 55 insertions(+), 20 deletions(-) diff --git a/Common/CostFunctions/itkParzenWindowHistogramImageToImageMetric.h b/Common/CostFunctions/itkParzenWindowHistogramImageToImageMetric.h index 544b76549..8628fbb9d 100644 --- a/Common/CostFunctions/itkParzenWindowHistogramImageToImageMetric.h +++ b/Common/CostFunctions/itkParzenWindowHistogramImageToImageMetric.h @@ -329,6 +329,13 @@ class ITK_TEMPLATE_EXPORT ParzenWindowHistogramImageToImageMetric const KernelFunctionType * kernel, ParzenValueContainerType & parzenValues) const; + /** Compute the Parzen values. Overload that places the values in the buffer, pointed to by the last argument. */ + static void + EvaluateParzenValues(double parzenWindowTerm, + OffsetValueType parzenWindowIndex, + const KernelFunctionType & kernel, + PDFValueType * parzenValues); + /** Update the joint PDF with a pixel pair; on demand also updates the * pdf derivatives (if the Jacobian pointers are nonzero). */ diff --git a/Common/CostFunctions/itkParzenWindowHistogramImageToImageMetric.hxx b/Common/CostFunctions/itkParzenWindowHistogramImageToImageMetric.hxx index da4ddb44c..0e50e78a6 100644 --- a/Common/CostFunctions/itkParzenWindowHistogramImageToImageMetric.hxx +++ b/Common/CostFunctions/itkParzenWindowHistogramImageToImageMetric.hxx @@ -518,7 +518,19 @@ ParzenWindowHistogramImageToImageMetric::EvaluateParz const KernelFunctionType * kernel, ParzenValueContainerType & parzenValues) const { - kernel->Evaluate(static_cast(parzenWindowIndex) - parzenWindowTerm, parzenValues.data_block()); + EvaluateParzenValues(parzenWindowTerm, parzenWindowIndex, *kernel, parzenValues.data_block()); +} // end EvaluateParzenValues() + + +template +void +ParzenWindowHistogramImageToImageMetric::EvaluateParzenValues( + double parzenWindowTerm, + OffsetValueType parzenWindowIndex, + const KernelFunctionType & kernel, + PDFValueType * const parzenValues) +{ + kernel.Evaluate(static_cast(parzenWindowIndex) - parzenWindowTerm, parzenValues); } // end EvaluateParzenValues() @@ -550,12 +562,19 @@ ParzenWindowHistogramImageToImageMetric::UpdateJointP static_cast(std::floor(movingImageParzenWindowTerm + this->m_MovingParzenTermToIndexOffset)); /** The Parzen values. */ - ParzenValueContainerType fixedParzenValues(this->m_JointPDFWindow.GetSize()[1]); - ParzenValueContainerType movingParzenValues(this->m_JointPDFWindow.GetSize()[0]); - this->EvaluateParzenValues( - fixedImageParzenWindowTerm, fixedImageParzenWindowIndex, this->m_FixedKernel, fixedParzenValues); - this->EvaluateParzenValues( - movingImageParzenWindowTerm, movingImageParzenWindowIndex, this->m_MovingKernel, movingParzenValues); + const auto numberOfFixedParzenValues = m_JointPDFWindow.GetSize()[1]; + const auto numberOfMovingParzenValues = m_JointPDFWindow.GetSize()[0]; + + // Create a buffer of Parzen values for both the fixed and the moving image. + const auto parzenValues = std::make_unique(numberOfFixedParzenValues + numberOfMovingParzenValues); + + PDFValueType * const fixedParzenValues = parzenValues.get(); + PDFValueType * const movingParzenValues = parzenValues.get() + numberOfFixedParzenValues; + + Self::EvaluateParzenValues( + fixedImageParzenWindowTerm, fixedImageParzenWindowIndex, *m_FixedKernel, fixedParzenValues); + Self::EvaluateParzenValues( + movingImageParzenWindowTerm, movingImageParzenWindowIndex, *m_MovingKernel, movingParzenValues); /** Position the JointPDFWindow. */ JointPDFIndexType pdfWindowIndex; @@ -572,10 +591,10 @@ ParzenWindowHistogramImageToImageMetric::UpdateJointP if (!imageJacobian) { /** Loop over the Parzen window region and increment the values. */ - for (unsigned int f = 0; f < fixedParzenValues.GetSize(); ++f) + for (unsigned int f = 0; f < numberOfFixedParzenValues; ++f) { const double fv = fixedParzenValues[f]; - for (unsigned int m = 0; m < movingParzenValues.GetSize(); ++m) + for (unsigned int m = 0; m < numberOfMovingParzenValues; ++m) { it.Value() += static_cast(fv * movingParzenValues[m]); ++it; @@ -586,7 +605,7 @@ ParzenWindowHistogramImageToImageMetric::UpdateJointP else { /** Compute the derivatives of the moving Parzen window. */ - ParzenValueContainerType derivativeMovingParzenValues(this->m_JointPDFWindow.GetSize()[0]); + ParzenValueContainerType derivativeMovingParzenValues(numberOfMovingParzenValues); this->EvaluateParzenValues(movingImageParzenWindowTerm, movingImageParzenWindowIndex, this->m_DerivativeMovingKernel, @@ -597,11 +616,11 @@ ParzenWindowHistogramImageToImageMetric::UpdateJointP /** Loop over the Parzen window region and increment the values * Also update the pdf derivatives. */ - for (unsigned int f = 0; f < fixedParzenValues.GetSize(); ++f) + for (unsigned int f = 0; f < numberOfFixedParzenValues; ++f) { const double fv = fixedParzenValues[f]; const double fv_et = fv / et; - for (unsigned int m = 0; m < movingParzenValues.GetSize(); ++m) + for (unsigned int m = 0; m < numberOfMovingParzenValues; ++m) { it.Value() += static_cast(fv * movingParzenValues[m]); this->UpdateJointPDFDerivatives(it.GetIndex(), fv_et * derivativeMovingParzenValues[m], *imageJacobian, *nzji); diff --git a/Components/Metrics/AdvancedMattesMutualInformation/itkParzenWindowMutualInformationImageToImageMetric.hxx b/Components/Metrics/AdvancedMattesMutualInformation/itkParzenWindowMutualInformationImageToImageMetric.hxx index 52e456e45..005c4ec8d 100644 --- a/Components/Metrics/AdvancedMattesMutualInformation/itkParzenWindowMutualInformationImageToImageMetric.hxx +++ b/Components/Metrics/AdvancedMattesMutualInformation/itkParzenWindowMutualInformationImageToImageMetric.hxx @@ -765,25 +765,34 @@ ParzenWindowMutualInformationImageToImageMetric::Upda const int movingParzenWindowIndex = static_cast(std::floor(movingImageParzenWindowTerm + this->m_MovingParzenTermToIndexOffset)); + const auto numberOfFixedParzenValues = Superclass::m_JointPDFWindow.GetSize()[1]; + const auto numberOfDerivedMovingParzenValues = Superclass::m_JointPDFWindow.GetSize()[0]; + + // Create a buffer of Parzen values for both the fixed and the moving image. + const auto parzenValues = + std::make_unique(numberOfFixedParzenValues + numberOfDerivedMovingParzenValues); + /** Compute the fixed Parzen values. */ - ParzenValueContainerType fixedParzenValues(this->m_JointPDFWindow.GetSize()[1]); - this->EvaluateParzenValues( - fixedImageParzenWindowTerm, fixedParzenWindowIndex, this->m_FixedKernel, fixedParzenValues); + PDFValueType * const fixedParzenValues = parzenValues.get(); + Superclass::EvaluateParzenValues( + fixedImageParzenWindowTerm, fixedParzenWindowIndex, *Superclass::m_FixedKernel, fixedParzenValues); /** Compute the derivatives of the moving Parzen window. */ - ParzenValueContainerType derivativeMovingParzenValues(this->m_JointPDFWindow.GetSize()[0]); - this->EvaluateParzenValues( - movingImageParzenWindowTerm, movingParzenWindowIndex, this->m_DerivativeMovingKernel, derivativeMovingParzenValues); + PDFValueType * const derivativeMovingParzenValues = parzenValues.get() + numberOfFixedParzenValues; + Superclass::EvaluateParzenValues(movingImageParzenWindowTerm, + movingParzenWindowIndex, + *Superclass::m_DerivativeMovingKernel, + derivativeMovingParzenValues); /** Get the moving image bin size. */ const double et = static_cast(this->m_MovingImageBinSize); /** Loop over the Parzen window region and increment sum. */ PDFValueType sum = 0.0; - for (unsigned int f = 0; f < fixedParzenValues.GetSize(); ++f) + for (unsigned int f = 0; f < numberOfFixedParzenValues; ++f) { const double fv_et = fixedParzenValues[f] / et; - for (unsigned int m = 0; m < derivativeMovingParzenValues.GetSize(); ++m) + for (unsigned int m = 0; m < numberOfDerivedMovingParzenValues; ++m) { sum += this->m_PRatioArray[f + fixedParzenWindowIndex][m + movingParzenWindowIndex] * fv_et * derivativeMovingParzenValues[m];