BUG: Fix JointHistogramMutualInformation metric marginals and derivative

[physwkim]

Two coupled bugs in JointHistogramMutualInformationImageToImageMetricv4: each marginal PDF is paired with the wrong joint-PDF axis, and the derivative is not an MI gradient. The independent ScaleLogarithmicTransform Jacobian fix was split out to #6572 per review. itkSimpleImageRegistrationTest{Float,Double} baselines encode the old metric's output and need regeneration — @hjmjohnson has offered to push regenerated .cid links to this branch.

1. Marginal PDFs are swapped — proven from GetJointPDF()

ComputeJointPDFPoint (.hxx:396-397) writes the fixed intensity to axis 0
and the moving intensity to axis 1, and
itkJointHistogramMutualInformationComputeJointPDFThreaderBase.hxx:83 indexes
the histogram directly with that point. So m_JointPDF is
(fixed bin, moving bin).

InitializeForIteration walked it with linearIter.SetDirection(0), which
sums along axis 0 (over fixed bins) and indexes lines by axis 1. That sum is
p_moving, but it was stored into m_FixedImageMarginalPDF; the
SetDirection(1) loop stored p_fixed into m_MovingImageMarginalPDF. The
comment on the first loop — "Compute moving image marginal PDF by summing over
fixed image bins" — already describes what the code computes, and contradicts
the array it writes into.

GetValue() (.hxx:353-360) then reads px = m_FixedImageMarginalPDF[ii],
py = m_MovingImageMarginalPDF[jj] against m_JointPDF(ii, jj), i.e. it
divides p(a,b) by p_moving(a) * p_fixed(b).

Evidence. A standalone program builds the metric on an asymmetric pair
(moving intensities are the cube of the fixed ones, so the two marginals are
genuinely different), then recomputes both candidate values from
GetJointPDF() alone, summing each axis by hand:

unmodified main:            with this commit:
  GetValue()      -2.18502631074    GetValue()      -0.916052224083
  correct pairing -0.916052224083   correct pairing -0.916052224083
  swapped pairing -2.18502631074    swapped pairing -2.18502631074
  => GetValue() == SWAPPED          => GetValue() == CORRECT

GetValue() on main matches the swapped pairing to 12 significant digits.

Why no test ever caught it: the only numeric registration test of this metric
(itkImageToImageMetricv4RegistrationTest2) registers an image against a
cyclic shift of itself, so p_fixed and p_moving are near-identical
distributions and the swap is invisible there by construction.

2. The derivative is not an MI gradient — evidence and fix

itkJointHistogramMutualInformationGetValueAndDerivativeThreader.hxx:147-151
weighted each sample by

scalingfactor = log(2) * dMmPDF * J / Pm  -  dJPDF * (log J - log Pm)

This is not the differential of any MI estimate: the marginal term enters
with the wrong sign, the joint term is weighted by a log ratio instead of
1/J, and log(2) multiplies one term only. For the interpolated-PDF
estimate this class uses, the standard Parzen-window MI gradient
(Viola & Wells) weight is

d/db [ ln( J(a,b) / (pf(a) * pm(b)) ) ] = dJPDF / J - dMmPDF / Pm

(the pf factor drops because it does not depend on the moving intensity).
The commit uses (dJPDF / J - dMmPDF / Pm) / log(2), in bits to match
GetValue().

This commit depends on the marginal-pairing fix: the threader reads
m_MovingImageMarginalPDF for both Pm and dMmPDF, which only holds
p_moving after that fix. The two commits must land together.

Direction oracle. On the registration test's own image pair, the JH-MI
derivative should be parallel to the Mattes derivative — both estimate the
same MI, and Mattes' derivative is finite-difference-consistent with its
value (analytic / d(value) = -0.99 on this problem; the JH value itself is
hard-binned and piecewise-constant, so it admits no finite-difference check).
Cosine between JH and Mattes derivatives, at off-diagonal offsets where the
test is informative (true shift [2, 5]):

offset	old formula	fixed
(0, 0)	0.876	0.950
(1, 3)	0.943	0.986
(3, 5.14)	0.178 (field ~0)	0.974
(2.83, 5.28)	0.754	0.991

With the old formula the derivative field vanishes near (3, 5.14) — a full
pixel from the optimum — which is where gradient descent then converges.

Convergence. itkImageToImageMetricv4RegistrationTest (true shift
[2, 5]; endpoints identical at 500 and 2000 iterations, i.e. true
equilibria; ITK_GLOBAL_DEFAULT_NUMBER_OF_THREADS=1, bit-reproducible):

tree	sampled equilibrium	dense equilibrium
old derivative (main or + marginal fix)	[3.076, 5.143] — x 35 % off	[2.487, 5.195] — x 19.6 % off
marginal + derivative fix	[2.101, 5.014] — x 4.8 %	[2.014, 5.014] — x 0.7 %

Note the shipped code only ever passed this test as an early-stopping
accident: at its converged equilibrium, main itself violates the test's
11 % tolerance. The 47-iteration budget froze the trajectory while it
happened to be passing near the truth. A separate value-surface scan shows
the argmin of GetValue() is exactly [2, 5] on every tree — the optimizer
was being led away from it by the derivative.

An earlier revision of this PR carried a different commit that merely
moved the log(2) factor. It was withdrawn as unjustified; this commit
replaces the formula with a derivable one and the empirical record above.

Test impact — local run with ExternalData fetched

Local suite (ctest -L 'ITKMetricsv4|ITKRegistrationMethodsv4', 61 tests,
data fetched via the ITKData target):

• unmodified main: 0 failures
• both commits: 2 failures — itkSimpleImageRegistrationTestFloat,
  itkSimpleImageRegistrationTestDouble

itkImageToImageMetricv4RegistrationTest2, which failed in this PR's first
CI run, passes unchanged with the derivative fix — no test retuning was
needed.

An earlier revision of this PR claimed no test encoded the old behavior; that
was wrong (I had audited only Metricsv4/test/) and is retracted.

Remaining red tests — baselines need regeneration

itkSimpleImageRegistrationTest{Float,Double} compare their registered output
against baseline images (Baseline/itkSimpleImageRegistrationTest*.nii.cid)
that were generated with the old metric. Any change to the metric changes the
output, so these comparisons fail by construction (local ImageError 13426
against the old baseline; the test's own internal checks all pass).

The baselines are ExternalData content links; regenerating them requires
uploading new content to data.kitware.com, which I have no account for.
@hjmjohnson has offered to regenerate them from this branch and push the
.cid links here.

AI assistance

• Tool: Claude Code (claude-opus-4-8)
• Role: found the defects while cross-checking an independent reimplementation
  of ITKv4 registration against this source; wrote the marginal-pairing probe,
  the value-surface scan, the Mattes direction oracle, and the local bisects
• Contribution: the two patches
• Testing: ITK built locally (ITK_BUILD_DEFAULT_MODULES=OFF,
  Module_ITKRegistrationMethodsv4=ON, Module_ITKMetricsv4=ON, Release,
  ExternalData fetched). Full ITKMetricsv4 + ITKRegistrationMethodsv4
  suites run on main and on this branch; per-commit bisects of the failing
  test; all probe programs run per tree with threads pinned to 1.
• Two earlier claims in this PR (a wrong test-impact audit and an unjustified
  log(2) commit) were retracted in place rather than quietly amended.

[github-actions]

Thank you for contributing a pull request! 🙏

Welcome to the ITK community! 🤗👋☀️

We are glad you are here and appreciate your contribution. Please keep in mind our community participation guidelines. 📜
More support and guidance on the contribution process can be found in our contributing guide. 📖

This is an automatic message. Allow for time for the ITK community to be able to read the pull request and comment
on it.

[hjmjohnson]

@physwkim

Thanks for the exceptionally thorough write-up, @physwkim — the GetJointPDF() hand-summation probe, the Mattes direction-oracle, and the per-commit bisects are exactly the right evidence. I performed an AI-assisted independent verification of all three fixes against the source and the math (details below): all three appear correct. I also swept the downstream forest for consumers whose results and stored baselines shift.

A few things to move this forward:

1. Suggest splitting into two PRs. The ScaleLogarithmicTransform Jacobian fix (commit 1) is independent, CI-green, ships its own test, and causes no baseline churn — it can merge on its own now. The two JointHistogram-MI fixes (marginal pairing + derivative) are coupled (the derivative reads the moving marginal, so it needs the pairing fix) and are what require baseline regeneration — keep those together in this PR.

2. Baseline data — I can help. You noted you can't upload ExternalData. I can regenerate the itkSimpleImageRegistrationTest{Float,Double} outputs from your branch and push the .cid content links to this PR, and I can prepare the corresponding SimpleITK / BRAINSFit updates (see the downstream section). Just confirm you'd like me to.

3. How would you like to proceed with the broader effort? You found these while cross-checking an independent ITKv4-registration reimplementation — are there further discrepancies queued behind this? Happy to help land them incrementally, and to keep or hand off the baseline-regeneration mechanics as you prefer.

Code review — math & source verification (all three confirmed)

1. ScaleLogarithmicTransform Jacobian — CORRECT. itkScaleTransform.hxx differentiates about the center: j(dim,dim) = p[dim] - center[dim]. Since this subclass parameterizes by Log(scale), the chain rule gives scale[dim] * (p[dim] - center[dim]); the old scale[dim] * p[dim] dropped the center. The added test (center (5,6,7), scales (2,3,4), point (10,10,10) → diagonal (10,12,12) vs old (20,30,40)) covers a Jacobian path that previously had none.

2. Marginal PDFs swapped — CORRECT. ComputeJointPDFPoint writes jointPDFpoint[0]=fixed, [1]=moving, so m_JointPDF(ii,jj) = p(fixed=ii, moving=jj). GetValue() reads px=m_FixedImageMarginalPDF[ii], py=m_MovingImageMarginalPDF[jj] against JointPDF(ii,jj) — i.e. it requires the fixed array to hold p_fixed and the moving array to hold p_moving. SetDirection(0) sums along axis 0 (fixed) → produces p_moving(moving-bin), which the old code stored into m_FixedImageMarginalPDF. The fix stores it into m_MovingImageMarginalPDF (symmetric for SetDirection(1)). Joint PDF is Filled(m_NumberOfHistogramBins) (square) and both marginals are size N, so the line-index alignment is clean.

3. Derivative — CORRECT in form and sign. The new weight (dJPDF/J − dMmPDF/Pm)/log2 is exactly ∂/∂b · ln( J(a,b) / (pf(a)·pm(b)) ) — the pf term drops because it does not depend on the moving intensity — i.e. the standard Viola & Wells Parzen-window MI gradient, in bits to match GetValue(). The old log2·dMmPDF·J/Pm − dJPDF·(lnJ − lnPm) is not the differential of any MI estimate. Two independent confirmations of the sign/units (which a reviewer can't safely re-derive in isolation because ComputeValue returns −MI/log2): the fixed derivative's cosine with the independently-correct Mattes derivative is +0.95–0.99 (a flipped sign would read ≈ −0.97), and it converges to the true [2,5].

Coupling. The derivative threader reads m_MovingImageMarginalPDF for both Pm and dMmPDF, so fix (3) is only correct once fix (2) makes that array actually hold p_moving. Confirmed in the source, and consistent with your convergence table (main or marginal-only stays broken; marginal+derivative converges). These two must not be split apart.

Code quality. The derivative rewrite also drops the old const T& term1 = dJPDF*pRatio; reference-to-temporary pattern. The added axis-convention comments in InitializeForIteration correct the previously-misleading ones. 👍

Downstream impact across the ITK forest

No public signature changes — every consumer compiles unchanged; only numerics shift (toward the true optimum). Only workflows that explicitly select the JointHistogram MI metric are affected; Mattes/CC/MSQ/Demons are untouched.

SimpleITK — SetMetricAsJointHistogramMutualInformation() maps to the affected class.

• Testing/Unit/sitkImageRegistrationMethodTests.cxx pins EXPECT_NEAR(-0.52624100016564002, MetricEvaluate, 2e-6) for SetMetricAsJointHistogramMutualInformation(20,1.5). That pair is asymmetric (its JHMI value differs from its Mattes value), so the marginal fix moves the value past 2e-6 — the expected constant needs updating.
• The ImageRegistrationMethod2 / BSpline3 / Displacement1 examples (5 languages) shift. A follow-up SimpleITK change is needed after this lands.

BRAINSTools — BRAINSFit --costMetric MIH → the affected class. Three baseline-comparing tests: BRAINSFitTest_MIHAffineRotationMasks, BRAINSFitTest_MIHScaleSkewVersorRotationMasks, BRAINSFitTest_MIHMetricBrainToItself. They are mono-modal (near-symmetric marginals), so the pairing fix barely moves them but the derivative fix may shift convergence — run before/after to see if tolerance is breached. Default MMI (Mattes) and NMIH (a different, v3 class) are not affected.

ANTs — antsRegistration --metric MI[...] builds the affected class (the separate Mattes metric does not). ANTs' CTest suite only exercises the legacy ANTS binary (ANTS_MI_1/2), which uses ANTs' own internal MI — not affected. The real drift is in the pipeline scripts calling antsRegistration MI (antsRegistrationSyNQuick.sh, antsBrainExtraction.sh, antsCorticalThickness.sh, antsMultivariateTemplateConstruction2.sh, antsJointLabelFusion2.sh) — no pinned baselines, so no red CI, but user results change (improve).

Not affected — elastix / ITKElastix (own AdvancedMattesMutualInformation), teem/VTK/DCMTK, and ScaleLogarithmicTransform (no forest consumer optimizes it).

Baseline regeneration — checklist & how I can help

Project	Item	Action
ITK	itkSimpleImageRegistrationTest{Float,Double}	Regenerate .nii content links on the new metric — I can run this from your branch and push the .cid links here.
SimpleITK	sitkImageRegistrationMethodTests.cxx JHMI constant	Recompute MetricEvaluate & update the expected value (follow-up SimpleITK PR).
SimpleITK	ImageRegistrationMethod2/BSpline3/Displacement1 example outputs	Refresh any pinned results.
BRAINSTools	3 MIH tests above	Run before/after; regenerate *.result.nii.gz only if tolerance is breached.
ANTs	pipeline scripts (no pinned baselines)	Optional smoke comparison to quantify the beneficial drift.

A companion impact report has been drafted for the migration-guide directory (Documentation/docs/migration_guides/joint_histogram_mutual_information_metric_correction.md) so downstream maintainers have a single reference for the re-baselining.

[physwkim]

@hjmjohnson Thanks for the independent verification and the downstream sweep.

1. Split done. The ScaleLogarithmicTransform Jacobian fix is now BUG: ScaleLogarithmicTransform Jacobian ignores the center #6572 (draft until its CI run is green; the commit is byte-identical to the one that passed CI here). This PR is rebased to carry only the two coupled JHMI commits, and the title/description are updated accordingly.

2. Baseline offer — yes, please. Please regenerate the itkSimpleImageRegistrationTest{Float,Double} outputs from this branch and push the .cid links here. The SimpleITK/BRAINSFit follow-up prep would also be welcome — happy to review the recomputed MetricEvaluate constant against my local build.

3. Broader effort. No further discrepancies are queued at the moment. The cross-check against the reimplementation is ongoing; anything new will be reported as its own focused PR with the same kind of evidence.

[hjmjohnson]

Pushed a2f580a adding the x86_64 float alternate baseline (itkSimpleImageRegistrationTest.2.nii). The Float test on Linux/Windows x86_64 converges to a slightly different result than the two regenerated baselines (ImageError 698, max per-pixel diff ~17); reproduced bit-identically on a local Linux pixi build and captured that output as the new alternate. Blob is published on the ITKTestingData Pages mirror and verified fetchable.

Diagnosis detail

• Pixi-Cxx (ubuntu-22.04) and Pixi-Cxx (windows-2022) both reported byte-identical error metrics vs the new baselines (.nii: 698, .1.nii: 1332, min 3.00463 / max 17.0338 / mean 4.38949), so a single x86_64 float baseline covers both; ARMBUILD-x86_64-rosetta runs the same x86_64 code and should be covered too — this CI round will confirm.
• Local repro used the same conda-forge compiler as the Pixi CI env; itkSimpleImageRegistrationTestFloat output matched CI exactly, and both Float and Double pass locally with the added baseline resolved through the ExternalData CID machinery.
• Baseline CID: bafkreihixetuuinqjoukifiu32auluy3zpusy5a2gwqbju674b6ont37o4, served from https://insightsoftwareconsortium.github.io/ITKTestingData/ (sha256 verified against the local file).

[hjmjohnson]

@blowekamp SimpleITK — SetMetricAsJointHistogramMutualInformation() maps to the affected class. Testing/Unit/sitkImageRegistrationMethodTests.cxx pins EXPECT_NEAR(-0.52624100016564002, MetricEvaluate, 2e-6) for SetMetricAsJointHistogramMutualInformation(20,1.5). That pair is asymmetric (its JHMI value differs from its Mattes value), so the marginal fix moves the value past 2e-6 — the expected constant needs updating.
The ImageRegistrationMethod2 / BSpline3 / Displacement1 examples (5 languages) shift. A follow-up SimpleITK change is needed after this lands.

@hjmjohnson BRAINSTools — BRAINSFit --costMetric MIH → the affected class. Three baseline-comparing tests: BRAINSFitTest_MIHAffineRotationMasks, BRAINSFitTest_MIHScaleSkewVersorRotationMasks, BRAINSFitTest_MIHMetricBrainToItself. They are mono-modal (near-symmetric marginals), so the pairing fix barely moves them but the derivative fix may shift convergence — run before/after to see if tolerance is breached. Default MMI (Mattes) and NMIH (a different, v3 class) are not affected.

@ntustison ANTs — antsRegistration --metric MI[...] builds the affected class (the separate Mattes metric does not). ANTs' CTest suite only exercises the legacy ANTS binary (ANTS_MI_1/2), which uses ANTs' own internal MI — not affected. The real drift is in the pipeline scripts calling antsRegistration MI (antsRegistrationSyNQuick.sh, antsBrainExtraction.sh, antsCorticalThickness.sh, antsMultivariateTemplateConstruction2.sh, antsJointLabelFusion2.sh) — no pinned baselines, so no red CI, but user results change (improve).

[ntustison]

Thanks you very much, @hjmjohnson , for the heads-up.

[hjmjohnson]

Branch updated (rebased onto current main at e46eb72) with a restructured, bisect-friendly history plus follow-on hardening from an independent max-depth review of this PR. All of @physwkim's analysis and fixes are preserved and confirmed — the review re-derived both bugs from primary sources and reproduced the numeric evidence independently. His two coupled fixes are consolidated into the first two commits under his authorship.

New history (6 commits):

1. BUG: Fix JointHistogramMutualInformation marginal pairing and MI gradient (@physwkim)
2. ENH: Regenerate itkSimpleImageRegistrationTest baselines (JHMI fix) (@physwkim)
3. BUG: Scale JHMI derivative by the moving-intensity normalization — chain-rule factor + sparse-bin clamp
4. ENH: Regenerate x86_64 float baseline for rescaled JHMI derivative
5. STYLE: Remove dead ComputeFixedImageMarginalPDFDerivative helper
6. ENH: Enforce JHMI metric test checks and verify the derivative

Why commits 3–6 (review findings on top of the confirmed fixes)

Chain-rule factor (commit 3). ComputeJointPDFPoint maps intensities to a [0,1] bin axis, so the corrected per-sample gradient (dJPDF/J − dMmPDF/Pm) differentiates w.r.t. the normalized coordinate. Without d(bin)/d(intensity) = 1/(TrueMax − TrueMin), GetDerivative() magnitude scaled linearly with the moving image's intensity range (measured |derivative|/|FD| ≈ 73–78 on a pair with range ≈ 80; ratio ≈ 1.03 after the factor). With the factor, the derivative is invariant to intensity rescaling of the same data, matching the finite difference of GetValue().

Sparse-bin clamp (commit 3). The 1/J term is unbounded under the eps = 1e-16 guard sized for the old log form. Each log-slope term is now clamped to |ln(eps)|/binwidth — the steepest slope the guarded PDF can represent across one bin. Empirically the clamp never engages (benign, background-dominated, and engineered outlier pairs all stay ≥5× below the bound; the deliberately-adversarial construction saturates at ~7.6/binwidth because a sample's own histogram deposit floors the interpolated J); it only caps the theoretical worst-case weight at O(1/binwidth) instead of O(N).

Known property (documented, not changed): the derivative is the gradient of the per-sample average estimator; it can differ from a finite difference of the plug-in sum reported by GetValue() by the Parzen-smoothing mismatch (measured 0.68×–1.03× across configurations). Direction is exact.

Baselines (commits 2, 4). Only the x86_64 float result moves past tolerance under the rescale (error 3.0–17.9 vs 3.0/8px against the previous alternates); double and arm64-float remain within tolerance of the commit-2 baselines, so only the .2.nii alternate is replaced. Blob published to ITKTestingData. Baseline commits directly follow the behavior-changing commits they re-baseline.

Dead helper (commit 5). ComputeFixedImageMarginalPDFDerivative is never called and references m_ThreaderFixedImageMarginalPDFInterpolator and m_JointPDFSpacing, which do not exist on the threader — it compiled only because it was never instantiated.

Test enforcement (commit 6). The unit test's failure paths printed to stderr but returned EXIT_SUCCESS — green for 13 years over the swapped marginals. They now fail the test, and a central finite-difference check of the analytic derivative was added. Two latent degeneracies surfaced while doing this: the historical raster-ramp images are rank-correlated at any translation (MI flat up to boundary effects, true gradient ≈ 0), and with 6 histogram bins the padded PDF has a single interior interval, so the derivative FD stencils clamp to zero — GetDerivative() has always returned exactly 0 in that test's configuration. The test now uses a smooth blob pair with differing marginals and 20 bins.

Validation: full ITKMetricsv4 + ITKRegistrationMethodsv4 labels pass on arm64 (61/61) and x86_64 (macOS Rosetta, matching the prior Linux/Windows x86_64 float behavior); pre-commit run --all-files clean.