From 814c79fa9354ae6145ecb034a2fc5950d869aa45 Mon Sep 17 00:00:00 2001 From: Hans Johnson Date: Thu, 9 Apr 2020 10:21:44 -0500 Subject: [PATCH 01/10] DOC: Fix typo in comment User -> Uses --- .../Registration/Metricsv4/include/itkImageToImageMetricv4.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Modules/Registration/Metricsv4/include/itkImageToImageMetricv4.h b/Modules/Registration/Metricsv4/include/itkImageToImageMetricv4.h index 27cf52599e5..0c96c228774 100644 --- a/Modules/Registration/Metricsv4/include/itkImageToImageMetricv4.h +++ b/Modules/Registration/Metricsv4/include/itkImageToImageMetricv4.h @@ -734,7 +734,7 @@ class ITK_TEMPLATE_EXPORT ImageToImageMetricv4 FixedImageGradientCalculatorPointer m_FixedImageGradientCalculator; MovingImageGradientCalculatorPointer m_MovingImageGradientCalculator; - /** Derivative results holder. User a raw pointer so we can point it + /** Derivative results holder. Uses a raw pointer so we can point it * to a user-provided object. This is used in internal methods so * the user-provided variable does not have to be passed around. It also enables * safely sharing a derivative object between metrics during multi-variate From f7a8395b008148d28bca87a11b139008084bfbb8 Mon Sep 17 00:00:00 2001 From: Hans Johnson Date: Thu, 9 Apr 2020 10:23:25 -0500 Subject: [PATCH 02/10] BUG: Pointer to temporary variable needs resetting to null. When the pointer to the local variable persists outside of the local function. When the local variable goes out of scope the mutable class variable should be set to nullptr. --- .../Metricsv4/include/itkImageToImageMetricv4.hxx | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Modules/Registration/Metricsv4/include/itkImageToImageMetricv4.hxx b/Modules/Registration/Metricsv4/include/itkImageToImageMetricv4.hxx index 2fcf966f16d..0d91677ca1d 100644 --- a/Modules/Registration/Metricsv4/include/itkImageToImageMetricv4.hxx +++ b/Modules/Registration/Metricsv4/include/itkImageToImageMetricv4.hxx @@ -206,15 +206,15 @@ typename ImageToImageMetricv4m_ComputeDerivative = false; - DerivativeType derivative; - this->m_DerivativeResult = &derivative; + DerivativeType local_derivative; + this->m_DerivativeResult = &local_derivative; this->InitializeForIteration(); // Do the threaded processing using the appropriate // GetValueAndDerivativeThreader. Results get written to // member vars. this->GetValueAndDerivativeExecute(); - + this->m_DerivativeResult = nullptr; // Pointer to temporary local_derivative is invalid after function return return this->m_Value; } From f2bd92a89ed994f316534a562b56e140659c7edd Mon Sep 17 00:00:00 2001 From: Hans Johnson Date: Thu, 9 Apr 2020 11:02:21 -0500 Subject: [PATCH 03/10] STYLE: Reorder commands for locallity of reference Pre-compute commmon logFixedImagePDFValue outside of inner loop. --- ...sMutualInformationImageToImageMetricv4.hxx | 24 +++++++++---------- 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/Modules/Registration/Metricsv4/include/itkMattesMutualInformationImageToImageMetricv4.hxx b/Modules/Registration/Metricsv4/include/itkMattesMutualInformationImageToImageMetricv4.hxx index d6999b83669..bbf3402d829 100644 --- a/Modules/Registration/Metricsv4/include/itkMattesMutualInformationImageToImageMetricv4.hxx +++ b/Modules/Registration/Metricsv4/include/itkMattesMutualInformationImageToImageMetricv4.hxx @@ -352,23 +352,23 @@ MattesMutualInformationImageToImageMetricv4m_ThreaderFixedImageMarginalPDF[0][bin] /= totalMassOfPDF; } - /** - * Compute the metric by double summation over histogram. - */ + static constexpr PDFValueType closeToZero = std::numeric_limits::epsilon(); + const PDFValueType nFactor = 1.0 / (this->m_MovingImageBinSize * this->GetNumberOfValidPoints()); // Setup pointer to point to the first bin - JointPDFValueType * jointPDFPtr = this->m_ThreaderJointPDF[0]->GetBufferPointer(); + const JointPDFValueType * jointPDFPtr = this->m_ThreaderJointPDF[0]->GetBufferPointer(); - // Initialize sum to zero + auto const temp_num_histogram_bins = this->m_NumberOfHistogramBins; + /** + * Compute the metric by double summation over histogram. + */ PDFValueType sum = 0.0; - - const PDFValueType nFactor = 1.0 / (this->m_MovingImageBinSize * this->GetNumberOfValidPoints()); - - static constexpr PDFValueType closeToZero = std::numeric_limits::epsilon(); - for (unsigned int fixedIndex = 0; fixedIndex < this->m_NumberOfHistogramBins; ++fixedIndex) + for (unsigned int fixedIndex = 0; fixedIndex < temp_num_histogram_bins; ++fixedIndex) { const PDFValueType fixedImagePDFValue = this->m_ThreaderFixedImageMarginalPDF[0][fixedIndex]; - for (unsigned int movingIndex = 0; movingIndex < this->m_NumberOfHistogramBins; ++movingIndex, jointPDFPtr++) + const PDFValueType logFixedImagePDFValue = std::log(fixedImagePDFValue); + for (unsigned int movingIndex = 0; movingIndex < temp_num_histogram_bins; + ++movingIndex, ++jointPDFPtr /* GOTO NEXT BIN */) { const PDFValueType movingImagePDFValue = this->m_MovingImageMarginalPDF[movingIndex]; const PDFValueType jointPDFValue = *(jointPDFPtr); @@ -382,7 +382,7 @@ MattesMutualInformationImageToImageMetricv4 closeToZero) { - sum += jointPDFValue * (pRatio - std::log(fixedImagePDFValue)); + sum += jointPDFValue * (pRatio - logFixedImagePDFValue); } if (this->GetComputeDerivative()) From 5e2bfa5e107334e4f01786e778949fa446a22727 Mon Sep 17 00:00:00 2001 From: Hans Johnson Date: Thu, 9 Apr 2020 14:38:41 -0500 Subject: [PATCH 04/10] ENH: Prefer multiplication to division Multiplication is more efficient than division, so use multiply by reciprical to make 1 division. --- .../include/itkMattesMutualInformationImageToImageMetricv4.hxx | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/Modules/Registration/Metricsv4/include/itkMattesMutualInformationImageToImageMetricv4.hxx b/Modules/Registration/Metricsv4/include/itkMattesMutualInformationImageToImageMetricv4.hxx index bbf3402d829..13931dbbcd3 100644 --- a/Modules/Registration/Metricsv4/include/itkMattesMutualInformationImageToImageMetricv4.hxx +++ b/Modules/Registration/Metricsv4/include/itkMattesMutualInformationImageToImageMetricv4.hxx @@ -347,9 +347,10 @@ MattesMutualInformationImageToImageMetricv4m_NumberOfHistogramBins; ++bin) { - this->m_ThreaderFixedImageMarginalPDF[0][bin] /= totalMassOfPDF; + this->m_ThreaderFixedImageMarginalPDF[0][bin] *= inv_totalMassOfPDF; } static constexpr PDFValueType closeToZero = std::numeric_limits::epsilon(); From 4383ea07eba937b96fc21e8e3c0d22c9dba80e92 Mon Sep 17 00:00:00 2001 From: Hans Johnson Date: Fri, 10 Apr 2020 08:52:44 -0500 Subject: [PATCH 05/10] ENH: Variable names updated to better indicate their pupopse Replace 'if(! condition) {continue;}' with 'if( condition ) { large block of code }' to remove early loop termination. --- ...sMutualInformationImageToImageMetricv4.hxx | 87 +++++++++---------- 1 file changed, 43 insertions(+), 44 deletions(-) diff --git a/Modules/Registration/Metricsv4/include/itkMattesMutualInformationImageToImageMetricv4.hxx b/Modules/Registration/Metricsv4/include/itkMattesMutualInformationImageToImageMetricv4.hxx index 13931dbbcd3..d5a80355638 100644 --- a/Modules/Registration/Metricsv4/include/itkMattesMutualInformationImageToImageMetricv4.hxx +++ b/Modules/Registration/Metricsv4/include/itkMattesMutualInformationImageToImageMetricv4.hxx @@ -366,54 +366,53 @@ MattesMutualInformationImageToImageMetricv4m_ThreaderFixedImageMarginalPDF[0][fixedIndex]; - const PDFValueType logFixedImagePDFValue = std::log(fixedImagePDFValue); - for (unsigned int movingIndex = 0; movingIndex < temp_num_histogram_bins; - ++movingIndex, ++jointPDFPtr /* GOTO NEXT BIN */) + const PDFValueType fixedImageMarginalPDFValue = this->m_ThreaderFixedImageMarginalPDF[0][fixedIndex]; { - const PDFValueType movingImagePDFValue = this->m_MovingImageMarginalPDF[movingIndex]; - const PDFValueType jointPDFValue = *(jointPDFPtr); - - // check for non-zero bin contribution - if (!(jointPDFValue > closeToZero && movingImagePDFValue > closeToZero)) - { - continue; - } - const PDFValueType pRatio = std::log(jointPDFValue / movingImagePDFValue); - - if (fixedImagePDFValue > closeToZero) + for (unsigned int movingIndex = 0; movingIndex < temp_num_histogram_bins; + ++movingIndex, ++jointPDFPtr /* GOTO NEXT BIN */) { - sum += jointPDFValue * (pRatio - logFixedImagePDFValue); - } - - if (this->GetComputeDerivative()) - { - if (!this->HasLocalSupport()) + const PDFValueType movingImageMarginalPDF = this->m_MovingImageMarginalPDF[movingIndex]; + const PDFValueType jointPDFValue = *(jointPDFPtr); + // check for non-zero bin contribution + if (jointPDFValue > closeToZero && movingImageMarginalPDF > closeToZero) { - // Collect global derivative contributions - - // move joint pdf derivative pointer to the right position - JointPDFValueType const * derivPtr = this->m_JointPDFDerivatives->GetBufferPointer() + - (fixedIndex * this->m_JointPDFDerivatives->GetOffsetTable()[2]) + - (movingIndex * this->m_JointPDFDerivatives->GetOffsetTable()[1]); - for (unsigned int parameter = 0, lastParameter = this->GetNumberOfLocalParameters(); - parameter < lastParameter; - ++parameter, derivPtr++) + const PDFValueType pRatio = std::log(jointPDFValue / movingImageMarginalPDF); + if (fixedImageMarginalPDFValue > closeToZero) { - // Ref: eqn 23 of Thevenaz & Unser paper [3] - (*(this->m_DerivativeResult))[parameter] += (*derivPtr) * pRatio; - } // end for-loop over parameters - } - else - { - // Collect the pRatio per pdf indices. - // Will be applied subsequently to local-support derivative - const OffsetValueType index = movingIndex + (fixedIndex * this->m_NumberOfHistogramBins); - this->m_PRatioArray[index] = pRatio * nFactor; - } - } - } // end for-loop over moving index - } // end for-loop over fixed index + const PDFValueType logfixedImageMarginalPDFValue = std::log(fixedImageMarginalPDFValue); + sum += jointPDFValue * (pRatio - logfixedImageMarginalPDFValue); + } + + if (this->GetComputeDerivative()) + { + if (!this->HasLocalSupport()) + { + // Collect global derivative contributions + + // move joint pdf derivative pointer to the right position + JointPDFValueType const * derivPtr = this->m_JointPDFDerivatives->GetBufferPointer() + + (fixedIndex * this->m_JointPDFDerivatives->GetOffsetTable()[2]) + + (movingIndex * this->m_JointPDFDerivatives->GetOffsetTable()[1]); + for (unsigned int parameter = 0, lastParameter = this->GetNumberOfLocalParameters(); + parameter < lastParameter; + ++parameter, derivPtr++) + { + // Ref: eqn 23 of Thevenaz & Unser paper [3] + (*(this->m_DerivativeResult))[parameter] += (*derivPtr) * pRatio; + } // end for-loop over parameters + } + else + { + // Collect the pRatio per pdf indices. + // Will be applied subsequently to local-support derivative + const OffsetValueType index = movingIndex + (fixedIndex * this->m_NumberOfHistogramBins); + this->m_PRatioArray[index] = pRatio * nFactor; + } + } + } // end if( jointPDFValue > closeToZero && movingImageMarginalPDF > closeToZero ) + } // end for-loop over moving index + } // end conditional for fixedmarginalPDF > close to zero + } // end for-loop over fixed index // Apply the pRatio and sum the per-window derivative // contributions, in the local-support case. From 1db27b6da6c5be4efcd6d96b165c5d1a2ff806a8 Mon Sep 17 00:00:00 2001 From: Hans Johnson Date: Fri, 10 Apr 2020 11:16:52 -0500 Subject: [PATCH 06/10] STYLE: Make alias variables to assist with clear documentation Use alias variable names that more clearly document what the intent of the workspace ivars that are re-used memory this stage of the processing. The Per-thread arrays of array accumulators have been reduced into thread 0 at this point, and will be normalized into PDF values by dividing by the total mass of the joint accumulator PDF. Also move exception checking for pre-existing states to the top of the function call ( there is no need to do computations before checking for pre-state exception cases). --- ...sMutualInformationImageToImageMetricv4.hxx | 71 +++++++++++++------ 1 file changed, 48 insertions(+), 23 deletions(-) diff --git a/Modules/Registration/Metricsv4/include/itkMattesMutualInformationImageToImageMetricv4.hxx b/Modules/Registration/Metricsv4/include/itkMattesMutualInformationImageToImageMetricv4.hxx index d5a80355638..20dcd52fd18 100644 --- a/Modules/Registration/Metricsv4/include/itkMattesMutualInformationImageToImageMetricv4.hxx +++ b/Modules/Registration/Metricsv4/include/itkMattesMutualInformationImageToImageMetricv4.hxx @@ -307,22 +307,40 @@ MattesMutualInformationImageToImageMetricv4::ComputeResults() const { + if (this->GetNumberOfValidPoints() == 0) + { + itkExceptionMacro("All samples map outside moving image buffer. " + "The images do not sufficiently " + "overlap. They need to be initialized to have more overlap before this " + "metric will work. For instance, you can align the image centers by translation." + << std::endl); + } if (this->m_JointPDFSum < itk::NumericTraits::epsilon()) { itkExceptionMacro("Joint PDF summed to zero"); } + const PDFValueType normalizationFactor = 1.0 / this->m_JointPDFSum; - std::fill(this->m_MovingImageMarginalPDF.begin(), this->m_MovingImageMarginalPDF.end(), 0.0F); + // Create aliases to make variable name intent more clear + // At this point the multiple thread partial values have been merged into + // the zero'th element of the m_ThreaderJointPDF and m_ThreaderFixedImageMarginalPDF. + const auto & l_JointPDF = this->m_ThreaderJointPDF[0]; + std::vector & l_FixedImageMarginalPDF = this->m_ThreaderFixedImageMarginalPDF[0]; - // Collect some results - PDFValueType totalMassOfPDF = 0.0; - for (unsigned int i = 0; i < this->m_NumberOfHistogramBins; ++i) - { - totalMassOfPDF += this->m_ThreaderFixedImageMarginalPDF[0][i]; - } + std::fill(this->m_MovingImageMarginalPDF.begin(), this->m_MovingImageMarginalPDF.end(), 0.0F); - const PDFValueType normalizationFactor = 1.0 / this->m_JointPDFSum; - JointPDFValueType * pdfPtr = this->m_ThreaderJointPDF[0]->GetBufferPointer(); + /* MovingMarginalPDF JointPDF + * ------------------- + * [ 3 ] | 1 | 2 | 0 | + * --- --- --- + * [ 5 ] | 4 | 1 | 0 | + * --- --- --- + * [ 0 ] | 0 | 0 | 0 | + * ------------------- + * + * [ 5 ] [ 3 ] [ 0 ] <-- FixedMarginalPDF + */ + JointPDFValueType * pdfPtr = l_JointPDF->GetBufferPointer(); for (unsigned int i = 0; i < this->m_NumberOfHistogramBins; ++i) { PDFValueType * movingMarginalPtr = &(m_MovingImageMarginalPDF[0]); @@ -332,32 +350,36 @@ MattesMutualInformationImageToImageMetricv4GetNumberOfValidPoints() == 0) + // Collect some results + PDFValueType totalMassOfPDF = 0.0; + for (unsigned int i = 0; i < this->m_NumberOfHistogramBins; ++i) { - itkExceptionMacro("All samples map outside moving image buffer. " - "The images do not sufficiently " - "overlap. They need to be initialized to have more overlap before this " - "metric will work. For instance, you can align the image centers by translation." - << std::endl); + totalMassOfPDF += l_FixedImageMarginalPDF[i]; + } + if (totalMassOfPDF != this->m_JointPDFSum) + { + std::cout << "ERROR " << totalMassOfPDF << " != " << this->m_JointPDFSum << std::endl; } // Normalize the fixed image marginal PDF - if (totalMassOfPDF == 0.0) + if (totalMassOfPDF > 0.0) { - itkExceptionMacro("Fixed image marginal PDF summed to zero"); + const PDFValueType inv_totalMassOfPDF = 1.0 / totalMassOfPDF; + for (unsigned int bin = 0; bin < this->m_NumberOfHistogramBins; ++bin) + { + l_FixedImageMarginalPDF[bin] *= inv_totalMassOfPDF; + } } - const PDFValueType inv_totalMassOfPDF = 1.0 / totalMassOfPDF; - for (unsigned int bin = 0; bin < this->m_NumberOfHistogramBins; ++bin) + else { - this->m_ThreaderFixedImageMarginalPDF[0][bin] *= inv_totalMassOfPDF; + itkExceptionMacro("Fixed image marginal PDF summed to zero"); } static constexpr PDFValueType closeToZero = std::numeric_limits::epsilon(); const PDFValueType nFactor = 1.0 / (this->m_MovingImageBinSize * this->GetNumberOfValidPoints()); // Setup pointer to point to the first bin - const JointPDFValueType * jointPDFPtr = this->m_ThreaderJointPDF[0]->GetBufferPointer(); + const JointPDFValueType * jointPDFPtr = l_JointPDF->GetBufferPointer(); auto const temp_num_histogram_bins = this->m_NumberOfHistogramBins; /** @@ -366,7 +388,10 @@ MattesMutualInformationImageToImageMetricv4m_ThreaderFixedImageMarginalPDF[0][fixedIndex]; + const PDFValueType fixedImageMarginalPDFValue = l_FixedImageMarginalPDF[fixedIndex]; + + // const bool isNotNearZerofixedImageMarginalPDFValue = (fixedImageMarginalPDFValue > closeToZero); + // if (isNotNearZerofixedImageMarginalPDFValue) { for (unsigned int movingIndex = 0; movingIndex < temp_num_histogram_bins; ++movingIndex, ++jointPDFPtr /* GOTO NEXT BIN */) From 397dbc019e4ff0634f27a3f15aac7600bb803363 Mon Sep 17 00:00:00 2001 From: Hans Johnson Date: Fri, 10 Apr 2020 13:34:20 -0500 Subject: [PATCH 07/10] STYLE: totalMassOfPDF == this->m_JointPDFSum by definition Simplify computations avoiding duplicate work to compute same value for totalMassOfPDF that is already available from this->m_JointPDFSum. This also allows for simplifying loops and the use of std::for_each transforms Avoid unnecessary indexing by looping over vectors with range-based loops. --- ...sMutualInformationImageToImageMetricv4.hxx | 55 +++++++------------ 1 file changed, 21 insertions(+), 34 deletions(-) diff --git a/Modules/Registration/Metricsv4/include/itkMattesMutualInformationImageToImageMetricv4.hxx b/Modules/Registration/Metricsv4/include/itkMattesMutualInformationImageToImageMetricv4.hxx index 20dcd52fd18..bc274dd60e8 100644 --- a/Modules/Registration/Metricsv4/include/itkMattesMutualInformationImageToImageMetricv4.hxx +++ b/Modules/Registration/Metricsv4/include/itkMattesMutualInformationImageToImageMetricv4.hxx @@ -327,10 +327,8 @@ MattesMutualInformationImageToImageMetricv4m_ThreaderJointPDF[0]; std::vector & l_FixedImageMarginalPDF = this->m_ThreaderFixedImageMarginalPDF[0]; - std::fill(this->m_MovingImageMarginalPDF.begin(), this->m_MovingImageMarginalPDF.end(), 0.0F); - - /* MovingMarginalPDF JointPDF - * ------------------- + /* FixedMarginalPDF JointPDF + * (j) ------------------- * [ 3 ] | 1 | 2 | 0 | * --- --- --- * [ 5 ] | 4 | 1 | 0 | @@ -338,42 +336,31 @@ MattesMutualInformationImageToImageMetricv4GetBufferPointer(); - for (unsigned int i = 0; i < this->m_NumberOfHistogramBins; ++i) - { - PDFValueType * movingMarginalPtr = &(m_MovingImageMarginalPDF[0]); - for (unsigned int j = 0; j < this->m_NumberOfHistogramBins; j++) - { - *(pdfPtr) *= normalizationFactor; - *(movingMarginalPtr++) += *(pdfPtr++); - } - } - // Collect some results - PDFValueType totalMassOfPDF = 0.0; - for (unsigned int i = 0; i < this->m_NumberOfHistogramBins; ++i) - { - totalMassOfPDF += l_FixedImageMarginalPDF[i]; - } - if (totalMassOfPDF != this->m_JointPDFSum) - { - std::cout << "ERROR " << totalMassOfPDF << " != " << this->m_JointPDFSum << std::endl; - } - // Normalize the fixed image marginal PDF - if (totalMassOfPDF > 0.0) + auto nomalize_labmda = [&normalizationFactor](const JointPDFValueType & c) -> JointPDFValueType { + return c * normalizationFactor; + }; + const size_t number_of_JointPDF_bins = this->m_NumberOfHistogramBins * this->m_NumberOfHistogramBins; + JointPDFValueType * pdfPtr = l_JointPDF->GetBufferPointer(); + std::transform(pdfPtr, pdfPtr + number_of_JointPDF_bins, pdfPtr, nomalize_labmda); + std::transform( + l_FixedImageMarginalPDF.cbegin(), l_FixedImageMarginalPDF.cend(), l_FixedImageMarginalPDF.begin(), nomalize_labmda); { - const PDFValueType inv_totalMassOfPDF = 1.0 / totalMassOfPDF; - for (unsigned int bin = 0; bin < this->m_NumberOfHistogramBins; ++bin) + size_t curr_col = 0; + for (auto & currMarginalBin : this->m_MovingImageMarginalPDF) { - l_FixedImageMarginalPDF[bin] *= inv_totalMassOfPDF; + currMarginalBin = 0; + const auto start = pdfPtr + curr_col++; + const auto stop = pdfPtr + number_of_JointPDF_bins; + for (auto colptr = start; colptr < stop; colptr += this->m_NumberOfHistogramBins) + { + currMarginalBin += *colptr; + } + currMarginalBin *= normalizationFactor; } } - else - { - itkExceptionMacro("Fixed image marginal PDF summed to zero"); - } static constexpr PDFValueType closeToZero = std::numeric_limits::epsilon(); const PDFValueType nFactor = 1.0 / (this->m_MovingImageBinSize * this->GetNumberOfValidPoints()); From 6790207e6a9749f9bd2f5074613820faba611254 Mon Sep 17 00:00:00 2001 From: Hans Johnson Date: Fri, 10 Apr 2020 14:02:29 -0500 Subject: [PATCH 08/10] PERF: Only compute log(fixedImageMarginalPDFValue once outside loop A dummy value is used in the case where the log computation is not used and would otherwise be invalid so that the log() can be precomputed outside the looping structure. --- ...sMutualInformationImageToImageMetricv4.hxx | 19 +++++++++++-------- 1 file changed, 11 insertions(+), 8 deletions(-) diff --git a/Modules/Registration/Metricsv4/include/itkMattesMutualInformationImageToImageMetricv4.hxx b/Modules/Registration/Metricsv4/include/itkMattesMutualInformationImageToImageMetricv4.hxx index bc274dd60e8..ec23b0fb534 100644 --- a/Modules/Registration/Metricsv4/include/itkMattesMutualInformationImageToImageMetricv4.hxx +++ b/Modules/Registration/Metricsv4/include/itkMattesMutualInformationImageToImageMetricv4.hxx @@ -377,21 +377,24 @@ MattesMutualInformationImageToImageMetricv4 closeToZero); - // if (isNotNearZerofixedImageMarginalPDFValue) + const bool isNotNearZerofixedImageMarginalPDFValue = (fixedImageMarginalPDFValue > closeToZero); + // NOTE: if isNotNearZerofixedImageMarginalPDFValue is false, logfixedImageMarginalPDFValue is never used + // The common case is that it is used, so only perform std::log(fixedImageMarginalPDFValue one time + const PDFValueType logfixedImageMarginalPDFValue = (isNotNearZerofixedImageMarginalPDFValue) + ? std::log(fixedImageMarginalPDFValue) + : std::numeric_limits::max(); { for (unsigned int movingIndex = 0; movingIndex < temp_num_histogram_bins; ++movingIndex, ++jointPDFPtr /* GOTO NEXT BIN */) { - const PDFValueType movingImageMarginalPDF = this->m_MovingImageMarginalPDF[movingIndex]; - const PDFValueType jointPDFValue = *(jointPDFPtr); - // check for non-zero bin contribution - if (jointPDFValue > closeToZero && movingImageMarginalPDF > closeToZero) + const PDFValueType & movingImageMarginalPDF = this->m_MovingImageMarginalPDF[movingIndex]; + const PDFValueType & jointPDFValue = *(jointPDFPtr); + // check for non-zero bin contribution, if movingImageMarginalPDF <= closeToZero, then so is joinPDFValue + if (movingImageMarginalPDF > closeToZero && jointPDFValue > closeToZero) { const PDFValueType pRatio = std::log(jointPDFValue / movingImageMarginalPDF); - if (fixedImageMarginalPDFValue > closeToZero) + if (isNotNearZerofixedImageMarginalPDFValue) { - const PDFValueType logfixedImageMarginalPDFValue = std::log(fixedImageMarginalPDFValue); sum += jointPDFValue * (pRatio - logfixedImageMarginalPDFValue); } From c7502c55ca862dbe9a1f204b7e6f02953b8835a8 Mon Sep 17 00:00:00 2001 From: Hans Johnson Date: Fri, 10 Apr 2020 17:38:06 -0500 Subject: [PATCH 09/10] PERF: Avoid processing entire rows based on mathematical certainties When either the fixedMarginalPDF or movingMarginalPDF (which are sums of rows or columns of the 2D joint histogram) then the intersection point must also be zero. Both the contribution to the final values of such a point is no change, so one can avoid much testing and subsequent inner conditional processing in such cases. --- ...sMutualInformationImageToImageMetricv4.hxx | 34 ++++++++----------- 1 file changed, 15 insertions(+), 19 deletions(-) diff --git a/Modules/Registration/Metricsv4/include/itkMattesMutualInformationImageToImageMetricv4.hxx b/Modules/Registration/Metricsv4/include/itkMattesMutualInformationImageToImageMetricv4.hxx index ec23b0fb534..0cd542ed3b9 100644 --- a/Modules/Registration/Metricsv4/include/itkMattesMutualInformationImageToImageMetricv4.hxx +++ b/Modules/Registration/Metricsv4/include/itkMattesMutualInformationImageToImageMetricv4.hxx @@ -365,9 +365,6 @@ MattesMutualInformationImageToImageMetricv4::epsilon(); const PDFValueType nFactor = 1.0 / (this->m_MovingImageBinSize * this->GetNumberOfValidPoints()); - // Setup pointer to point to the first bin - const JointPDFValueType * jointPDFPtr = l_JointPDF->GetBufferPointer(); - auto const temp_num_histogram_bins = this->m_NumberOfHistogramBins; /** * Compute the metric by double summation over histogram. @@ -377,37 +374,36 @@ MattesMutualInformationImageToImageMetricv4 closeToZero); - // NOTE: if isNotNearZerofixedImageMarginalPDFValue is false, logfixedImageMarginalPDFValue is never used - // The common case is that it is used, so only perform std::log(fixedImageMarginalPDFValue one time - const PDFValueType logfixedImageMarginalPDFValue = (isNotNearZerofixedImageMarginalPDFValue) - ? std::log(fixedImageMarginalPDFValue) - : std::numeric_limits::max(); + // If fixedImageMarginalPDFValue is <= closeToZero, then the entire row of values is <= closeToZero, so + // by definition, each of the individual jointPDFVlaue's must also be close to zero in the line below! + if (fixedImageMarginalPDFValue > closeToZero) { - for (unsigned int movingIndex = 0; movingIndex < temp_num_histogram_bins; - ++movingIndex, ++jointPDFPtr /* GOTO NEXT BIN */) + // NOTE: If fixedImageMarginalPDFValue <=> closeToZero, logfixedImageMarginalPDFValue is never used + // The common case is that it is used, so perform std::log(fixedImageMarginalPDFValue) one time + const PDFValueType logfixedImageMarginalPDFValue = std::log(fixedImageMarginalPDFValue); + // Setup pointer to point to the first bin of this row in the jointPDF + const JointPDFValueType * jointPDFRowPtr = l_JointPDF->GetBufferPointer() + fixedIndex * temp_num_histogram_bins; + for (unsigned int movingIndex = 0; movingIndex < temp_num_histogram_bins; ++movingIndex) { const PDFValueType & movingImageMarginalPDF = this->m_MovingImageMarginalPDF[movingIndex]; - const PDFValueType & jointPDFValue = *(jointPDFPtr); + const PDFValueType jointPDFValue = *(jointPDFRowPtr++); /* Get Value and goto next bin in row */ + // check for non-zero bin contribution, if movingImageMarginalPDF <= closeToZero, then so is joinPDFValue - if (movingImageMarginalPDF > closeToZero && jointPDFValue > closeToZero) + if (movingImageMarginalPDF > closeToZero && + jointPDFValue > closeToZero) //<-- This check is always false if not isNotNearZerofixedImageMarginalPDFValue { const PDFValueType pRatio = std::log(jointPDFValue / movingImageMarginalPDF); - if (isNotNearZerofixedImageMarginalPDFValue) - { - sum += jointPDFValue * (pRatio - logfixedImageMarginalPDFValue); - } + sum += jointPDFValue * (pRatio - logfixedImageMarginalPDFValue); if (this->GetComputeDerivative()) { if (!this->HasLocalSupport()) { // Collect global derivative contributions - - // move joint pdf derivative pointer to the right position JointPDFValueType const * derivPtr = this->m_JointPDFDerivatives->GetBufferPointer() + (fixedIndex * this->m_JointPDFDerivatives->GetOffsetTable()[2]) + (movingIndex * this->m_JointPDFDerivatives->GetOffsetTable()[1]); + // move joint pdf derivative pointer to the right position for (unsigned int parameter = 0, lastParameter = this->GetNumberOfLocalParameters(); parameter < lastParameter; ++parameter, derivPtr++) From aa534653dd586dd5d7b29bea27d4ee72e17bd627 Mon Sep 17 00:00:00 2001 From: Hans Johnson Date: Fri, 10 Apr 2020 21:53:00 -0500 Subject: [PATCH 10/10] ENH: Provide new baseline for ImageRegisration4Test This test started failing due to very small numerical accumulation changes. Other packages (i.e. ANTs & BRAINSTools) registration tools that are more complete versions of this example all pass their more extensive testing. --- .../Data/Baseline/Registration/ImageRegistration4Test.png.md5 | 1 - .../Baseline/Registration/ImageRegistration4Test.png.sha512 | 2 +- 2 files changed, 1 insertion(+), 2 deletions(-) delete mode 100644 Testing/Data/Baseline/Registration/ImageRegistration4Test.png.md5 diff --git a/Testing/Data/Baseline/Registration/ImageRegistration4Test.png.md5 b/Testing/Data/Baseline/Registration/ImageRegistration4Test.png.md5 deleted file mode 100644 index b0a71bccc8d..00000000000 --- a/Testing/Data/Baseline/Registration/ImageRegistration4Test.png.md5 +++ /dev/null @@ -1 +0,0 @@ -1fa40e8d13e22f06fefb518a9115acbd diff --git a/Testing/Data/Baseline/Registration/ImageRegistration4Test.png.sha512 b/Testing/Data/Baseline/Registration/ImageRegistration4Test.png.sha512 index 2c4e50f5f6d..9e2656fbbb8 100644 --- a/Testing/Data/Baseline/Registration/ImageRegistration4Test.png.sha512 +++ b/Testing/Data/Baseline/Registration/ImageRegistration4Test.png.sha512 @@ -1 +1 @@ -227d438599c355b42bbfc913ed46ea05d52b0605dc756c9560e849e11235b82c77d3651ed280d5a12c1e08e57f157ef9c852e31d94fc7331096219695f2a9f3b +f9fe1b3d00a4ff711683d69672f36577bbcba4438f104aa99518907db4fac53d55dbbdac157ade1c1c90f1cc0a8614dad40a7aa140b7deef374a61310c0cfcd5