diff --git a/Documentation/docs/migration_guides/index.md b/Documentation/docs/migration_guides/index.md index b8af3118227..7af7696e24d 100644 --- a/Documentation/docs/migration_guides/index.md +++ b/Documentation/docs/migration_guides/index.md @@ -8,4 +8,5 @@ These migration guides explain how to update major versions of ITK, whic may con itk_6_migration_guide itk_5_migration_guide +joint_histogram_mutual_information_metric_correction ``` diff --git a/Documentation/docs/migration_guides/itk_6_migration_guide.md b/Documentation/docs/migration_guides/itk_6_migration_guide.md index 3660a323336..ccdcebf121c 100644 --- a/Documentation/docs/migration_guides/itk_6_migration_guide.md +++ b/Documentation/docs/migration_guides/itk_6_migration_guide.md @@ -954,3 +954,15 @@ unchanged; consumers compile without modification. re-read. - The vendored `gdcm::SerieHelper` is untouched (GDCM still uses it internally); it may be removed once upstream GDCM drops it. + +## `JointHistogramMutualInformationImageToImageMetricv4` corrected + +The JointHistogram mutual-information metric's marginal PDFs were swapped +(since introduction) and its analytic derivative was not an MI gradient; both +are fixed, and the derivative is now unit-consistent with a bounded per-sample +weight. The public API is unchanged and code compiles unmodified, but +registration values and convergence with this metric change (toward the true +optimum), and hand-tuned fixed learning rates must be re-tuned. Details, the +downstream-consumer impact survey (SimpleITK, BRAINSFit, ANTs), and the +re-baselining checklist are in the companion guide: +[JointHistogramMutualInformation metric correction](./joint_histogram_mutual_information_metric_correction.md). diff --git a/Documentation/docs/migration_guides/joint_histogram_mutual_information_metric_correction.md b/Documentation/docs/migration_guides/joint_histogram_mutual_information_metric_correction.md new file mode 100644 index 00000000000..545bed82701 --- /dev/null +++ b/Documentation/docs/migration_guides/joint_histogram_mutual_information_metric_correction.md @@ -0,0 +1,181 @@ +JointHistogramMutualInformation Metric Correction — Downstream Impact +==================================================================== + +This report documents a behavior change in +`itk::JointHistogramMutualInformationImageToImageMetricv4` introduced by +[PR #6569](https://github.com/InsightSoftwareConsortium/ITK/pull/6569), and +the downstream consumers whose numerical output — and stored regression +baselines — change as a result. + +It is a companion to the [ITK v6 Migration Guide](./itk_6_migration_guide.md): +these are *correctness* fixes, but registration results computed with the +affected metric change (they move toward the true optimum), and any test that +pins the old numbers must be re-baselined. One protected API is removed: +`JointHistogramMutualInformationGetValueAndDerivativeThreader::ComputeFixedImageMarginalPDFDerivative` +(and its per-thread `FixedImageMarginalPDFInterpolator`) had no remaining +callers after the marginal-pairing fix; subclasses that called it should use +`ComputeMovingImageMarginalPDFDerivative` against the correct marginal. + + +Summary +------- + +Three correctness defects are corrected (the related +`ScaleLogarithmicTransform` Jacobian fix was split out to +[PR #6572](https://github.com/InsightSoftwareConsortium/ITK/pull/6572)): + +1. **`JointHistogramMutualInformation` marginal PDFs were swapped.** + The joint PDF is indexed `(fixed bin, moving bin)`. Summing over the fixed + axis yields the *moving* marginal, but that result was stored into + `m_FixedImageMarginalPDF` (and vice-versa). `GetValue()` then divided + `p(a,b)` by `p_moving(a) * p_fixed(b)`. The value is therefore wrong for any + image pair whose two intensity marginals differ; it is (nearly) invisible + when the two marginals are the same, e.g. registering an image to a shifted + or rotated copy of itself. + +2. **The metric derivative was not an MI gradient.** + The per-sample weight is now the standard Parzen-window (Viola & Wells) + form `(dJPDF/J - dMmPDF/Pm) / log(2)`, replacing a weight that mixed a + log-ratio-weighted joint term with a wrong-signed marginal term. Because the + derivative reads the *moving* marginal, fix (2) depends on fix (1): the two + must land together. On the ITK registration test the corrected derivative + converges to within <1% of the true shift versus ~35% for the old formula, + and its direction agrees with the independently-correct Mattes metric + (cosine 0.95–0.99). + +3. **The derivative was not unit-consistent and could spike on sparse bins.** + `ComputeJointPDFPoint()` maps intensities to a `[0,1]` bin axis, so the + per-sample gradient needs the chain-rule factor + `d(bin)/d(intensity) = 1/(TrueMax - TrueMin)`. With it, `GetDerivative()` + matches a central finite difference of `GetValue()` and no longer scales + with the moving image's intensity range (previously an 8-bit vs 16-bit + encoding of the same data produced a 256x different derivative). Each + log-slope term is additionally clamped to `|ln(eps)|/binwidth` — the + steepest slope the `eps`-guarded PDF can represent across one bin — which + never engages on well-posed data and only bounds the worst-case per-sample + weight. + +Fixes (1) and (2) are coupled and change the JointHistogram MI metric's value +**and** gradient; fix (3) further rescales the gradient. + + +Derivative magnitude and learning rates +--------------------------------------- + +The corrected, unit-consistent derivative differs in magnitude from the old +formula — typically by orders of magnitude, depending on the moving image's +intensity range. + +- Pipelines using a **scales estimator** (e.g. + `RegistrationParameterScalesFromPhysicalShift`, the + `ImageRegistrationMethodv4` default, or SimpleITK's + `SetOptimizerScalesFromPhysicalShift()` + + `EstimateLearningRate`) absorb the rescale automatically — ITK's own + registration tests converge to results within existing tolerances on all + platforms except one float configuration. +- Pipelines with a **hand-tuned fixed learning rate** and no scales estimator + must re-tune it: the old rate was calibrated against a derivative that was + neither a gradient nor unit-consistent. Prefer switching to a scales + estimator over re-tuning a constant. +- The derivative remains the gradient of the per-sample average estimator; it + can differ from a finite difference of the plug-in histogram sum reported by + `GetValue()` by the Parzen-smoothing mismatch (measured 0.68x–1.03x across + histogram configurations). Direction is exact. + +Who is affected +--------------- + +Only workflows that **explicitly select the JointHistogram MI metric** change. +The Mattes MI metric (`itk::MattesMutualInformationImageToImageMetricv4`), +correlation, mean-squares, and demons metrics are untouched. No public +signature changes, so all downstream code compiles unchanged. + +### ITK + +| Test | Nature | Status | +|------|--------|--------| +| `itkSimpleImageRegistrationTestFloat` / `Double` | baseline-image compare | **Baselines regenerated in PR #6569** (primary double, `.1` float, `.2` x86_64 float). | +| `itkJointHistogramMutualInformationImageToImageMetricv4Test` | value/derivative consistency + finite-difference gradient check (hardened in PR #6569; failures now fail the test) | Passes. | +| `itkImageToImageMetricv4RegistrationTest`, `…Test2`, `itkSimpleImageRegistrationTest2/3`, `…WithMaskAndSampling`, `itkQuasiNewtonOptimizerv4RegistrationTest`, `itkMultiStartImageToImageMetricv4RegistrationTest` | convergence / no-exception | Pass (convergence improves). | + +### SimpleITK + +`SetMetricAsJointHistogramMutualInformation()` maps directly to the affected +class. + +- `Testing/Unit/sitkImageRegistrationMethodTests.cxx` pins a metric value: + `EXPECT_NEAR(-0.52624100016564002, MetricEvaluate, 2e-6)` for + `SetMetricAsJointHistogramMutualInformation(20, 1.5)`. The pair is asymmetric + (its JHMI value differs from its Mattes value), so the marginal fix moves the + value past the `2e-6` tolerance. **The expected constant must be updated.** +- The `ImageRegistrationMethod2`, `ImageRegistrationMethodBSpline3`, and + `ImageRegistrationMethodDisplacement1` examples (C++/Python/Java/C#/R) are the + canonical JHMI examples; their final metric/transform shifts. +- Requires a follow-up SimpleITK change after this lands in ITK. + +### BRAINSTools (BRAINSFit `--costMetric MIH`) + +Three baseline-comparing tests in `BRAINSFit/TestSuite/CMakeLists.txt`: + +- `BRAINSFitTest_MIHAffineRotationMasks` +- `BRAINSFitTest_MIHScaleSkewVersorRotationMasks` +- `BRAINSFitTest_MIHMetricBrainToItself` + +These are mono-modal (brain to a rotated brain / to itself), so the two +marginals are nearly identical and fix (2) barely moves them; fix (3) may still +shift the converged transform. Run them against a before/after ITK to see +whether they breach tolerance. The BRAINSFit default `MMI` (Mattes) and the +`NMIH` metric (which uses the v3 +`itk::NormalizedMutualInformationHistogramImageToImageMetric`, a different +class) are **not** affected. + +### ANTs (`antsRegistration --metric MI[...]`) + +In `itkantsRegistrationHelper.hxx` the `MI` metric enum builds the affected +JointHistogram class (the separate `Mattes` metric does not). However: + +- ANTs' own CTest regression suite exercises only the **legacy `ANTS` + binary** (`ANTS_MI_1/2_test.cmake`, `COMMAND ANTS … -m MI[…] -r Gauss`), + which uses ANTs' *internal* MI implementation, not the ITK class — **not + affected**. No `antsRegistration`/`antsAI`/`MeasureImageSimilarity` + exact-baseline test was found. +- The user-facing impact is in the pipeline scripts that call + `antsRegistration MI` — `antsRegistrationSyNQuick.sh`, + `antsBrainExtraction.sh`, `antsCorticalThickness.sh`, + `antsMultivariateTemplateConstruction2.sh`, `antsJointLabelFusion2.sh` — whose + outputs change (converging closer to the true alignment). These are not + pinned baselines, so no CI turns red, but results drift for users. + +### Not affected + +- **elastix / ITKElastix** — no use of the ITK JointHistogram metric (their MI + is `AdvancedMattesMutualInformation`). +- Non-registration toolkits (teem, VTK, DCMTK, …) — no metric consumers. +- `ScaleLogarithmicTransform` — its Jacobian fix moved to + [PR #6572](https://github.com/InsightSoftwareConsortium/ITK/pull/6572); no + forest consumer optimizes it (BRAINSTools references it only for + transform-file I/O). + + +Re-baselining checklist +----------------------- + +| Project | Item | Action | +|---------|------|--------| +| ITK | `itkSimpleImageRegistrationTest{Float,Double}` baselines | Done in PR #6569 (`.nii` content links regenerated, blobs in ITKTestingData). | +| SimpleITK | `sitkImageRegistrationMethodTests.cxx` JHMI constant | Recompute `MetricEvaluate` and update the expected value. | +| SimpleITK | `ImageRegistrationMethod2` / `BSpline3` / `Displacement1` example outputs | Refresh any pinned example results. | +| BRAINSTools | `BRAINSFitTest_MIH{AffineRotationMasks,ScaleSkewVersorRotationMasks,MetricBrainToItself}` | Run before/after; regenerate `*.result.nii.gz` if tolerance is breached. | +| ANTs | pipeline-script outputs (no pinned baselines) | Optional smoke comparison to quantify the (beneficial) drift. | + + +Verification notes +------------------ + +The correction moves results **toward** the true optimum, not away — the +corrected derivative agrees in direction with ITK's independently-correct +Mattes metric, and on the tested pairs the metric value's argmin is essentially +unchanged (the marginal swap offsets the value by a term that varies only +weakly with alignment). Downstream +baselines that shift are encoding the *previous, incorrect* behavior and should +be regenerated rather than treated as regressions. diff --git a/Modules/Registration/Metricsv4/include/itkJointHistogramMutualInformationGetValueAndDerivativeThreader.h b/Modules/Registration/Metricsv4/include/itkJointHistogramMutualInformationGetValueAndDerivativeThreader.h index f5459f02f99..ef58b3ba218 100644 --- a/Modules/Registration/Metricsv4/include/itkJointHistogramMutualInformationGetValueAndDerivativeThreader.h +++ b/Modules/Registration/Metricsv4/include/itkJointHistogramMutualInformationGetValueAndDerivativeThreader.h @@ -103,9 +103,6 @@ class ITK_TEMPLATE_EXPORT JointHistogramMutualInformationGetValueAndDerivativeTh DerivativeType & localDerivativeReturn, const ThreadIdType threadId) const override; - inline InternalComputationValueType - ComputeFixedImageMarginalPDFDerivative(const MarginalPDFPointType & margPDFpoint, const ThreadIdType threadId) const; - inline InternalComputationValueType ComputeMovingImageMarginalPDFDerivative(const MarginalPDFPointType & margPDFpoint, const ThreadIdType threadId) const; @@ -116,7 +113,6 @@ class ITK_TEMPLATE_EXPORT JointHistogramMutualInformationGetValueAndDerivativeTh struct JointHistogramMIPerThreadStruct { JointPDFInterpolatorPointer JointPDFInterpolator; - MarginalPDFInterpolatorPointer FixedImageMarginalPDFInterpolator; MarginalPDFInterpolatorPointer MovingImageMarginalPDFInterpolator; }; itkPadStruct(ITK_CACHE_LINE_ALIGNMENT, JointHistogramMIPerThreadStruct, PaddedJointHistogramMIPerThreadStruct); diff --git a/Modules/Registration/Metricsv4/include/itkJointHistogramMutualInformationGetValueAndDerivativeThreader.hxx b/Modules/Registration/Metricsv4/include/itkJointHistogramMutualInformationGetValueAndDerivativeThreader.hxx index b8d3633a24d..5f27739ee4e 100644 --- a/Modules/Registration/Metricsv4/include/itkJointHistogramMutualInformationGetValueAndDerivativeThreader.hxx +++ b/Modules/Registration/Metricsv4/include/itkJointHistogramMutualInformationGetValueAndDerivativeThreader.hxx @@ -20,6 +20,9 @@ #include "itkMakeUniqueForOverwrite.h" +#include // For clamp. +#include // For log. + namespace itk { @@ -60,13 +63,6 @@ JointHistogramMutualInformationGetValueAndDerivativeThreaderm_JointHistogramMIPerThreadVariables[i].JointPDFInterpolator->SetInputImage( this->m_JointAssociate->m_JointPDF); - if (this->m_JointHistogramMIPerThreadVariables[i].FixedImageMarginalPDFInterpolator.IsNull()) - { - this->m_JointHistogramMIPerThreadVariables[i].FixedImageMarginalPDFInterpolator = - MarginalPDFInterpolatorType::New(); - } - this->m_JointHistogramMIPerThreadVariables[i].FixedImageMarginalPDFInterpolator->SetInputImage( - this->m_JointAssociate->m_FixedImageMarginalPDF); if (this->m_JointHistogramMIPerThreadVariables[i].MovingImageMarginalPDFInterpolator.IsNull()) { this->m_JointHistogramMIPerThreadVariables[i].MovingImageMarginalPDFInterpolator = @@ -144,11 +140,14 @@ JointHistogramMutualInformationGetValueAndDerivativeThreader< constexpr InternalComputationValueType eps{ 1.e-16 }; if (jointPDFValue > eps && movingImagePDFValue > eps) { - const InternalComputationValueType pRatio = std::log(jointPDFValue) - std::log(movingImagePDFValue); - const InternalComputationValueType & term1 = dJPDF * pRatio; - const InternalComputationValueType & term2 = - this->m_JointAssociate->m_Log2 * dMmPDF * jointPDFValue / movingImagePDFValue; - scalingfactor = (term2 - term1); + // Steepest log-slope the eps-guarded PDF can represent across one bin. + const InternalComputationValueType maxSlope = -std::log(eps) / this->m_JointAssociate->m_JointPDFSpacing[ind]; + const InternalComputationValueType term1 = std::clamp(dJPDF / jointPDFValue, -maxSlope, maxSlope); + const InternalComputationValueType term2 = std::clamp(dMmPDF / movingImagePDFValue, -maxSlope, maxSlope); + // Chain rule for ComputeJointPDFPoint's intensity -> [0,1] bin-axis normalization. + const InternalComputationValueType movingIntensityRange = + this->m_JointAssociate->m_MovingImageTrueMax - this->m_JointAssociate->m_MovingImageTrueMin; + scalingfactor = (term1 - term2) / (this->m_JointAssociate->m_Log2 * movingIntensityRange); } // end if-block to check non-zero bin contribution else { @@ -177,51 +176,6 @@ JointHistogramMutualInformationGetValueAndDerivativeThreader< return true; } -template -typename JointHistogramMutualInformationGetValueAndDerivativeThreader< - TDomainPartitioner, - TImageToImageMetric, - TJointHistogramMetric>::InternalComputationValueType -JointHistogramMutualInformationGetValueAndDerivativeThreader< - TDomainPartitioner, - TImageToImageMetric, - TJointHistogramMetric>::ComputeFixedImageMarginalPDFDerivative(const MarginalPDFPointType & margPDFpoint, - const ThreadIdType threadId) const -{ - InternalComputationValueType offset = 0.5 * this->m_JointPDFSpacing[0]; - InternalComputationValueType eps = this->m_JointPDFSpacing[0]; - MarginalPDFPointType leftpoint = margPDFpoint; - leftpoint[0] -= offset; - MarginalPDFPointType rightpoint = margPDFpoint; - rightpoint[0] += offset; - if (leftpoint[0] < eps) - { - leftpoint[0] = eps; - } - if (rightpoint[0] < eps) - { - rightpoint[0] = eps; - } - if (leftpoint[0] > 1.0) - { - leftpoint[0] = 1.0; - } - if (rightpoint[0] > 1.0) - { - rightpoint[0] = 1.0; - } - InternalComputationValueType delta = rightpoint[0] - leftpoint[0]; - if (delta > InternalComputationValueType{}) - { - InternalComputationValueType deriv = - this->m_ThreaderFixedImageMarginalPDFInterpolator[threadId]->Evaluate(rightpoint) - - this->m_ThreaderFixedImageMarginalPDFInterpolator[threadId]->Evaluate(leftpoint); - return deriv / delta; - } - - return InternalComputationValueType{}; -} - template typename JointHistogramMutualInformationGetValueAndDerivativeThreader< TDomainPartitioner, diff --git a/Modules/Registration/Metricsv4/include/itkJointHistogramMutualInformationImageToImageMetricv4.hxx b/Modules/Registration/Metricsv4/include/itkJointHistogramMutualInformationImageToImageMetricv4.hxx index d71a504fd1a..96fe04c28b4 100644 --- a/Modules/Registration/Metricsv4/include/itkJointHistogramMutualInformationImageToImageMetricv4.hxx +++ b/Modules/Registration/Metricsv4/include/itkJointHistogramMutualInformationImageToImageMetricv4.hxx @@ -251,11 +251,11 @@ JointHistogramMutualInformationImageToImageMetricv4m_JointPDF = (dg->GetOutput()); } - // Compute moving image marginal PDF by summing over fixed image bins. + // Summing joint-PDF axis 0 (fixed bins) yields the moving marginal. ImageLinearIteratorWithIndex linearIter(m_JointPDF, m_JointPDF->GetBufferedRegion()); linearIter.SetDirection(0); linearIter.GoToBegin(); - unsigned int fixedIndex = 0; + unsigned int movingIndex = 0; CompensatedSummation sum; while (!linearIter.IsAtEnd()) { @@ -266,15 +266,16 @@ JointHistogramMutualInformationImageToImageMetricv4SetPixel(mind, static_cast(sum.GetSum())); + mind[0] = movingIndex; + m_MovingImageMarginalPDF->SetPixel(mind, static_cast(sum.GetSum())); linearIter.NextLine(); - ++fixedIndex; + ++movingIndex; } + // Compute fixed image marginal PDF by summing over moving image bins. linearIter.SetDirection(1); linearIter.GoToBegin(); - unsigned int movingIndex = 0; + unsigned int fixedIndex = 0; while (!linearIter.IsAtEnd()) { sum.ResetToZero(); @@ -284,10 +285,10 @@ JointHistogramMutualInformationImageToImageMetricv4SetPixel(mind, static_cast(sum.GetSum())); + mind[0] = fixedIndex; + m_FixedImageMarginalPDF->SetPixel(mind, static_cast(sum.GetSum())); linearIter.NextLine(); - ++movingIndex; + ++fixedIndex; } } diff --git a/Modules/Registration/Metricsv4/test/itkJointHistogramMutualInformationImageToImageMetricv4Test.cxx b/Modules/Registration/Metricsv4/test/itkJointHistogramMutualInformationImageToImageMetricv4Test.cxx index ef301bcc55d..59aa55bb3e9 100644 --- a/Modules/Registration/Metricsv4/test/itkJointHistogramMutualInformationImageToImageMetricv4Test.cxx +++ b/Modules/Registration/Metricsv4/test/itkJointHistogramMutualInformationImageToImageMetricv4Test.cxx @@ -21,17 +21,15 @@ #include "itkTranslationTransform.h" #include "itkTestingMacros.h" -/* Simple test to verify that class builds and runs. - * Results are not verified. See ImageToImageMetricv4Test - * for verification of basic metric functionality. - * - * TODO Numerical verification. - */ +/* Checks value/derivative path consistency, agreement of the analytic + * derivative with a central finite difference of the value, and the + * non-overlap contract. See ImageToImageMetricv4Test for basic metric + * functionality. */ int itkJointHistogramMutualInformationImageToImageMetricv4Test(int, char *[]) { - constexpr unsigned int imageSize{ 10 }; + constexpr unsigned int imageSize{ 20 }; constexpr unsigned int imageDimensionality{ 3 }; using ImageType = itk::Image; @@ -47,24 +45,22 @@ itkJointHistogramMutualInformationImageToImageMetricv4Test(int, char *[]) movingImage->SetRegions(region); movingImage->Allocate(); - /* Fill images */ - itk::ImageRegionIterator itFixed(fixedImage, region); - itFixed.GoToBegin(); - unsigned int count = 1; - while (!itFixed.IsAtEnd()) - { - itFixed.Set(count * count); - count++; - ++itFixed; - } - itk::ImageRegionIteratorWithIndex itMoving(movingImage, region); - itMoving.GoToBegin(); - count = 1; - while (!itMoving.IsAtEnd()) + /* Fill with a smooth blob plus gradients; moving is a monotone nonlinear + * remap of fixed so the marginals differ and MI varies with alignment. */ + itk::ImageRegionIteratorWithIndex itFixed(fixedImage, region); + itk::ImageRegionIterator itMoving(movingImage, region); + constexpr double center = (imageSize - 1) / 2.0; + for (itFixed.GoToBegin(), itMoving.GoToBegin(); !itFixed.IsAtEnd(); ++itFixed, ++itMoving) { - itMoving.Set((count)); - count++; - ++itMoving; + const ImageType::IndexType index = itFixed.GetIndex(); + double radius2 = 0.0; + for (unsigned int dim = 0; dim < imageDimensionality; ++dim) + { + radius2 += ((index[dim] - center) / (imageSize / 3.0)) * ((index[dim] - center) / (imageSize / 3.0)); + } + const double value = std::exp(-radius2) + 0.3 * index[0] / imageSize + 0.15 * index[1] / imageSize; + itFixed.Set(value); + itMoving.Set(value * value * value); } /* Transforms */ @@ -81,7 +77,8 @@ itkJointHistogramMutualInformationImageToImageMetricv4Test(int, char *[]) ITK_EXERCISE_BASIC_OBJECT_METHODS(metric, JointHistogramMutualInformationImageToImageMetricv4, ImageToImageMetricv4); - constexpr itk::SizeValueType numberOfHistogramBins{ 6 }; + // Fewer than 7 bins degenerates the padded-PDF derivative stencils to zero. + constexpr itk::SizeValueType numberOfHistogramBins{ 20 }; metric->SetNumberOfHistogramBins(numberOfHistogramBins); ITK_TEST_SET_GET_VALUE(numberOfHistogramBins, metric->GetNumberOfHistogramBins()); @@ -110,6 +107,37 @@ itkJointHistogramMutualInformationImageToImageMetricv4Test(int, char *[]) if (itk::Math::NotExactlyEquals(valueReturn1, valueReturn2)) { std::cerr << "Value return results are not identical: " << valueReturn1 << ", " << valueReturn2 << std::endl; + return EXIT_FAILURE; + } + + // Analytic derivative must point downhill and match a central finite difference. + { + constexpr double detune = 1.0; + constexpr double h = 0.5; + MovingTransformType::ParametersType offsetParameters(movingTransform->GetNumberOfParameters(), 0.0); + offsetParameters[0] = detune; + movingTransform->SetParameters(offsetParameters); + MetricType::MeasureType detunedValue = NAN; + MetricType::DerivativeType detunedDerivative; + metric->GetValueAndDerivative(detunedValue, detunedDerivative); + offsetParameters[0] = detune + h; + movingTransform->SetParameters(offsetParameters); + const MetricType::MeasureType valuePlus = metric->GetValue(); + offsetParameters[0] = detune - h; + movingTransform->SetParameters(offsetParameters); + const MetricType::MeasureType valueMinus = metric->GetValue(); + movingTransform->SetIdentity(); + const double finiteDifference = (valuePlus - valueMinus) / (2.0 * h); + const double analytic = detunedDerivative[0]; + if (itk::Math::abs(finiteDifference) < 1e-10 || analytic * finiteDifference >= 0.0 || + itk::Math::abs(analytic) < 0.3 * itk::Math::abs(finiteDifference) || + itk::Math::abs(analytic) > 3.0 * itk::Math::abs(finiteDifference)) + { + std::cerr << "Analytic derivative disagrees with finite difference of the value:" << std::endl + << " analytic[0]: " << analytic << std::endl + << " central finite difference: " << finiteDifference << std::endl; + return EXIT_FAILURE; + } } std::cout << "JointPDF: " << metric->GetJointPDF() << std::endl; @@ -128,6 +156,7 @@ itkJointHistogramMutualInformationImageToImageMetricv4Test(int, char *[]) << " Number of valid points: " << metric->GetNumberOfValidPoints() << std::endl << " Metric value: " << valueReturn2 << std::endl << " Expected metric max value: " << expectedMetricMax << std::endl; + return EXIT_FAILURE; } movingTransform->SetIdentity(); diff --git a/Modules/Registration/RegistrationMethodsv4/test/Baseline/itkSimpleImageRegistrationTest.1.nii.cid b/Modules/Registration/RegistrationMethodsv4/test/Baseline/itkSimpleImageRegistrationTest.1.nii.cid index cb9262735b2..0302bcf6794 100644 --- a/Modules/Registration/RegistrationMethodsv4/test/Baseline/itkSimpleImageRegistrationTest.1.nii.cid +++ b/Modules/Registration/RegistrationMethodsv4/test/Baseline/itkSimpleImageRegistrationTest.1.nii.cid @@ -1 +1 @@ -bafkreiagtlczegb2znsj5ohtooykgzqrdcquc3fqwncsxmvietu6eyea7e +bafkreidario4d4czuebxmx2eqahh7enromv22i4c76buhgol2h5ysccofe diff --git a/Modules/Registration/RegistrationMethodsv4/test/Baseline/itkSimpleImageRegistrationTest.2.nii.cid b/Modules/Registration/RegistrationMethodsv4/test/Baseline/itkSimpleImageRegistrationTest.2.nii.cid index 2ac2bddb3c1..cfa77126bea 100644 --- a/Modules/Registration/RegistrationMethodsv4/test/Baseline/itkSimpleImageRegistrationTest.2.nii.cid +++ b/Modules/Registration/RegistrationMethodsv4/test/Baseline/itkSimpleImageRegistrationTest.2.nii.cid @@ -1 +1 @@ -bafkreibuvr4ual4tvllt356minn6dikbeagwsggukxlp44zs7x7pjjat34 +bafkreiab3v6m4fy5jaypgfmamk727vc3fxiilvnocmotfnegv3zt5kedyu diff --git a/Modules/Registration/RegistrationMethodsv4/test/Baseline/itkSimpleImageRegistrationTest.3.nii.cid b/Modules/Registration/RegistrationMethodsv4/test/Baseline/itkSimpleImageRegistrationTest.3.nii.cid deleted file mode 100644 index e6680a21f33..00000000000 --- a/Modules/Registration/RegistrationMethodsv4/test/Baseline/itkSimpleImageRegistrationTest.3.nii.cid +++ /dev/null @@ -1 +0,0 @@ -bafkreic3fpf34kla3gduhp3cvam5jdhjjk4uyr3gyap6jekuoi5qxurhim diff --git a/Modules/Registration/RegistrationMethodsv4/test/Baseline/itkSimpleImageRegistrationTest.nii.cid b/Modules/Registration/RegistrationMethodsv4/test/Baseline/itkSimpleImageRegistrationTest.nii.cid index 1b06ec166dd..0e48aec29a8 100644 --- a/Modules/Registration/RegistrationMethodsv4/test/Baseline/itkSimpleImageRegistrationTest.nii.cid +++ b/Modules/Registration/RegistrationMethodsv4/test/Baseline/itkSimpleImageRegistrationTest.nii.cid @@ -1 +1 @@ -bafkreiafvfkt24cokys3irj4fzcobsa6ke54hooczbcuccmccz62loaoq4 +bafkreigqpynjmaslsa2ad5bvfvajnnpkbpzmjvwgb45ujq4tz6rokj2xo4 diff --git a/Testing/Data/content-links.manifest b/Testing/Data/content-links.manifest index 6c7bd546938..eb1bf4661d4 100644 --- a/Testing/Data/content-links.manifest +++ b/Testing/Data/content-links.manifest @@ -1355,10 +1355,9 @@ bafkreiaam33ofytsf3wmp4gakpnmwgtgtrxiyn4no7gxz7ueg4s2bbvh6e Modules/Registration bafkreieheksshyog3ddid5dj4jnda4j5pppllybwfrvzukcedeihrpka6m Modules/Registration/Montage/test/Input/Synthetic/Rect2N.png bafkreig5bchpphdtty3g2mmvmte7y6p2toygbwnk2e3b7pk4dd34awgj5a Modules/Registration/Montage/test/Input/VisibleHumanRGB/TileConfiguration.registered.txt bafkreig5bchpphdtty3g2mmvmte7y6p2toygbwnk2e3b7pk4dd34awgj5a Modules/Registration/Montage/test/Input/VisibleHumanRGB/TileConfiguration.txt -bafkreiagtlczegb2znsj5ohtooykgzqrdcquc3fqwncsxmvietu6eyea7e Modules/Registration/RegistrationMethodsv4/test/Baseline/itkSimpleImageRegistrationTest.1.nii -bafkreibuvr4ual4tvllt356minn6dikbeagwsggukxlp44zs7x7pjjat34 Modules/Registration/RegistrationMethodsv4/test/Baseline/itkSimpleImageRegistrationTest.2.nii -bafkreic3fpf34kla3gduhp3cvam5jdhjjk4uyr3gyap6jekuoi5qxurhim Modules/Registration/RegistrationMethodsv4/test/Baseline/itkSimpleImageRegistrationTest.3.nii -bafkreiafvfkt24cokys3irj4fzcobsa6ke54hooczbcuccmccz62loaoq4 Modules/Registration/RegistrationMethodsv4/test/Baseline/itkSimpleImageRegistrationTest.nii +bafkreidario4d4czuebxmx2eqahh7enromv22i4c76buhgol2h5ysccofe Modules/Registration/RegistrationMethodsv4/test/Baseline/itkSimpleImageRegistrationTest.1.nii +bafkreiab3v6m4fy5jaypgfmamk727vc3fxiilvnocmotfnegv3zt5kedyu Modules/Registration/RegistrationMethodsv4/test/Baseline/itkSimpleImageRegistrationTest.2.nii +bafkreigqpynjmaslsa2ad5bvfvajnnpkbpzmjvwgb45ujq4tz6rokj2xo4 Modules/Registration/RegistrationMethodsv4/test/Baseline/itkSimpleImageRegistrationTest.nii bafkreiblhixcduhuia4lfjykffuwvahnauwb54y6rr3lne4ksfukcrwqki Modules/Registration/RegistrationMethodsv4/test/Baseline/r16slice_rigid.nii.gz bafkreibqg6nybsvkxnnpgjxeuhscitfrj7e4t4kttamprrpvd6hm4o7bu4 Modules/Registration/RegistrationMethodsv4/test/Baseline/r64slice.nii.gz bafkreibue3a2zelpvucixsnxg53ft6mk7eu75qz5tlbqqtutsilec52qzu Modules/Registration/RegistrationMethodsv4/test/Input/Blob1.mha