Skip to content

Commit

Permalink
STYLE: Consolidate duplicated code
Browse files Browse the repository at this point in the history 
Place common code in a private utility member
function that can be called in both places.
  • Loading branch information
@hjmjohnson
hjmjohnson committed Oct 20, 2022
1 parent 0d98b79 commit 6b03920
Show file tree
Hide file tree
Showing 2 changed files with 59 additions and 80 deletions.
4 changes: 4 additions & 0 deletions Modules/Registration/Common/include/itkMattesMutualInformationImageToImageMetric.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -278,6 +278,10 @@ class ITK_TEMPLATE_EXPORT MattesMutualInformationImageToImageMetric
using CubicBSplineFunctionType = BSplineKernelFunction<3, PDFValueType>;
using CubicBSplineDerivativeFunctionType = BSplineDerivativeKernelFunction<3, PDFValueType>;

/** Extract common processing for both GetValueAndDerivative and GetValue functions */
void
CommonGetValueProcessing() const;

/** Precompute fixed image parzen window indices. */
void
ComputeFixedImageParzenWindowIndices(FixedImageSampleContainer & samples);
Expand Down
135 changes: 55 additions & 80 deletions Modules/Registration/Common/include/itkMattesMutualInformationImageToImageMetric.hxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -511,48 +511,7 @@ MattesMutualInformationImageToImageMetric<TFixedImage, TMovingImage>::GetValue(c
itkExceptionMacro("Joint PDF summed to zero\n" << this->m_MMIMetricPerThreadVariables[0].JointPDF);
}

std::fill(this->m_MovingImageMarginalPDF.begin(), this->m_MovingImageMarginalPDF.end(), 0.0F);

// NOTE: Since the m_ThreaderFixedImageMarginalPDF is the sum of mass
// in the fixed image dimension, accumulating these values gives the
// same answer as computing the the sum of individual values over
// the the entire histogram. IMPORTANT NOTICE: THIS MAKES AN
// ASSUMPTION OF CONSERVATION OF MASS OF THE BSPLINE SMOOTHING. The
// sum of all the values should equal the number of samples being
// used, since each sample contributes only one sample somewhere in
// the histogram
PDFValueType totalMassOfPDF = 0.0;
for (unsigned int i = 0; i < this->m_NumberOfHistogramBins; ++i)
{
totalMassOfPDF += this->m_MMIMetricPerThreadVariables[0].FixedImageMarginalPDF[i];
}
const PDFValueType normalizationFactor = 1.0 / this->m_MMIMetricPerThreadVariables[0].JointPDFSum;
JointPDFValueType * pdfPtr = this->m_MMIMetricPerThreadVariables[0].JointPDF->GetBufferPointer();
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++);
}
}

if (this->m_NumberOfPixelsCounted < this->m_NumberOfFixedImageSamples / 16)
{
itkExceptionMacro("Too many samples map outside moving image buffer: "
<< this->m_NumberOfPixelsCounted << " / " << this->m_NumberOfFixedImageSamples << std::endl);
}

// Normalize the fixed image marginal PDF
if (totalMassOfPDF == 0.0)
{
itkExceptionMacro("Fixed image marginal PDF summed to zero");
}
for (unsigned int bin = 0; bin < this->m_NumberOfHistogramBins; ++bin)
{
this->m_MMIMetricPerThreadVariables[0].FixedImageMarginalPDF[bin] /= totalMassOfPDF;
}
this->CommonGetValueProcessing();

/**
* Compute the metric by double sumdation over histogram.
Expand Down Expand Up @@ -776,52 +735,20 @@ MattesMutualInformationImageToImageMetric<TFixedImage, TMovingImage>::GetValueAn
// MUST BE CALLED TO INITIATE PROCESSING ON SAMPLES
this->GetValueAndDerivativeMultiThreadedInitiate();

// MUST BE CALLED TO INITIATE PROCESSING
// CALL IF DOING THREADED POST PROCESSING
this->GetValueAndDerivativeMultiThreadedPostProcessInitiate();
// Consolidate to the first element in the vector
for (ThreadIdType workUnitID = 1; workUnitID < this->m_NumberOfWorkUnits; ++workUnitID)
{
this->m_MMIMetricPerThreadVariables[0].JointPDFSum += this->m_MMIMetricPerThreadVariables[workUnitID].JointPDFSum;
}
if (this->m_MMIMetricPerThreadVariables[0].JointPDFSum < itk::NumericTraits<PDFValueType>::epsilon())
{
itkExceptionMacro("Joint PDF summed to zero");
}

std::fill(this->m_MovingImageMarginalPDF.begin(), this->m_MovingImageMarginalPDF.end(), 0.0F);

PDFValueType totalMassOfPDF = 0.0;
for (unsigned int i = 0; i < this->m_NumberOfHistogramBins; ++i)
{
totalMassOfPDF += this->m_MMIMetricPerThreadVariables[0].FixedImageMarginalPDF[i];
}

const PDFValueType normalizationFactor = 1.0 / this->m_MMIMetricPerThreadVariables[0].JointPDFSum;
JointPDFValueType * pdfPtr = this->m_MMIMetricPerThreadVariables[0].JointPDF->GetBufferPointer();
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++);
}
}

if (this->m_NumberOfPixelsCounted < this->m_NumberOfFixedImageSamples / 16)
{
itkExceptionMacro("Too many samples map outside moving image buffer: "
<< this->m_NumberOfPixelsCounted << " / " << this->m_NumberOfFixedImageSamples << std::endl);
itkExceptionMacro("Joint PDF summed to zero\n" << this->m_MMIMetricPerThreadVariables[0].JointPDF);
}

// Normalize the fixed image marginal PDF
if (totalMassOfPDF == 0.0)
{
itkExceptionMacro("Fixed image marginal PDF summed to zero");
}
for (unsigned int bin = 0; bin < this->m_NumberOfHistogramBins; ++bin)
{
this->m_MMIMetricPerThreadVariables[0].FixedImageMarginalPDF[bin] /= totalMassOfPDF;
}
this->CommonGetValueProcessing();

/**
* Compute the metric by double summation over histogram.
Expand All @@ -831,8 +758,7 @@ MattesMutualInformationImageToImageMetric<TFixedImage, TMovingImage>::GetValueAn
JointPDFValueType * jointPDFPtr = this->m_MMIMetricPerThreadVariables[0].JointPDF->GetBufferPointer();

// Initialize sum to zero
PDFValueType sum = 0.0;

PDFValueType sum = 0.0;
const PDFValueType nFactor = 1.0 / (this->m_MovingImageBinSize * this->m_NumberOfPixelsCounted);
for (unsigned int fixedIndex = 0; fixedIndex < this->m_NumberOfHistogramBins; ++fixedIndex)
{
Expand Down Expand Up @@ -903,6 +829,55 @@ MattesMutualInformationImageToImageMetric<TFixedImage, TMovingImage>::GetValueAn
value = static_cast<MeasureType>(-1.0 * sum);
}

template <typename TFixedImage, typename TMovingImage>
void
MattesMutualInformationImageToImageMetric<TFixedImage, TMovingImage>::CommonGetValueProcessing() const
{
std::fill(this->m_MovingImageMarginalPDF.begin(), this->m_MovingImageMarginalPDF.end(), 0.0F);

// NOTE: Since the m_ThreaderFixedImageMarginalPDF is the sum of mass
// in the fixed image dimension, accumulating these values gives the
// same answer as computing the the sum of individual values over
// the the entire histogram. IMPORTANT NOTICE: THIS MAKES AN
// ASSUMPTION OF CONSERVATION OF MASS OF THE BSPLINE SMOOTHING. The
// sum of all the values should equal the number of samples being
// used, since each sample contributes only one sample somewhere in
// the histogram
PDFValueType totalMassOfPDF = 0.0;
for (unsigned int i = 0; i < this->m_NumberOfHistogramBins; ++i)
{
totalMassOfPDF += this->m_MMIMetricPerThreadVariables[0].FixedImageMarginalPDF[i];
}

const PDFValueType normalizationFactor = 1.0 / this->m_MMIMetricPerThreadVariables[0].JointPDFSum;
JointPDFValueType * pdfPtr = this->m_MMIMetricPerThreadVariables[0].JointPDF->GetBufferPointer();
for (unsigned int i = 0; i < this->m_NumberOfHistogramBins; ++i)
{
PDFValueType * movingMarginalPtr = &(this->m_MovingImageMarginalPDF[0]);
for (unsigned int j = 0; j < this->m_NumberOfHistogramBins; ++j)
{
*(pdfPtr) *= normalizationFactor;
*(movingMarginalPtr++) += *(pdfPtr++);
}
}

if (this->m_NumberOfPixelsCounted < this->m_NumberOfFixedImageSamples / 16)
{
itkExceptionMacro("Too many samples map outside moving image buffer: "
<< this->m_NumberOfPixelsCounted << " / " << this->m_NumberOfFixedImageSamples << std::endl);
}

// Normalize the fixed image marginal PDF
if (totalMassOfPDF == 0.0)
{
itkExceptionMacro("Fixed image marginal PDF summed to zero");
}
for (unsigned int bin = 0; bin < this->m_NumberOfHistogramBins; ++bin)
{
this->m_MMIMetricPerThreadVariables[0].FixedImageMarginalPDF[bin] /= totalMassOfPDF;
}
}

/**
* Get the match measure derivative
*/
Expand Down

0 comments on commit 6b03920

Please sign in to comment.