BUG: Fix JointHistogramMutualInformation metric marginals and derivative

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
 ```

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).

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.

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);

Modules/Registration/Metricsv4/include/itkJointHistogramMutualInformationGetValueAndDerivativeThreader.hxx

@@ -20,6 +20,9 @@
 
 #include "itkMakeUniqueForOverwrite.h"
 
+#include <algorithm> // For clamp.
+#include <cmath>     // For log.
+
 namespace itk
 {
 
@@ -60,13 +63,6 @@ JointHistogramMutualInformationGetValueAndDerivativeThreader<TDomainPartitioner,
     }
     this->m_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 TDomainPartitioner, typename TImageToImageMetric, typename TJointHistogramMetric>
-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 TDomainPartitioner, typename TImageToImageMetric, typename TJointHistogramMetric>
 typename JointHistogramMutualInformationGetValueAndDerivativeThreader<
   TDomainPartitioner,

Modules/Registration/Metricsv4/include/itkJointHistogramMutualInformationImageToImageMetricv4.hxx

@@ -251,11 +251,11 @@ JointHistogramMutualInformationImageToImageMetricv4<TFixedImage,
     this->m_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<TInternalComputationValueType> sum;
   while (!linearIter.IsAtEnd())
   {
@@ -266,15 +266,16 @@ JointHistogramMutualInformationImageToImageMetricv4<TFixedImage,
       ++linearIter;
     }
     MarginalPDFIndexType mind;
-    mind[0] = fixedIndex;
-    m_FixedImageMarginalPDF->SetPixel(mind, static_cast<PDFValueType>(sum.GetSum()));
+    mind[0] = movingIndex;
+    m_MovingImageMarginalPDF->SetPixel(mind, static_cast<PDFValueType>(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 @@ JointHistogramMutualInformationImageToImageMetricv4<TFixedImage,
       ++linearIter;
     }
     MarginalPDFIndexType mind;
-    mind[0] = movingIndex;
-    m_MovingImageMarginalPDF->SetPixel(mind, static_cast<PDFValueType>(sum.GetSum()));
+    mind[0] = fixedIndex;
+    m_FixedImageMarginalPDF->SetPixel(mind, static_cast<PDFValueType>(sum.GetSum()));
     linearIter.NextLine();
-    ++movingIndex;
+    ++fixedIndex;
   }
 }