ENH: Make Mattes metric derivatives deterministic with multiple threads

[cookpa]

Description of changes from Codex:

Implemented deterministic multithreaded Mattes mutual information derivatives.

Cause: derivative buffers were merged into shared storage according to mutex acquisition and thread completion order, producing schedule-dependent floating-point summation.

Changes:

• Accumulate derivative PDFs independently per work unit.
• Merge them after threading in ascending work-unit ID order.
• Remove the mutex-based DerivativeBufferManager and FinalizeThread() reduction.
• Add test that requires exact equality across multiple runs.
• Preserve existing deterministic joint-PDF merging behavior.

Description by cookpa:

Mattes mutual information was non-deterministic with multiple threads, even with dense sampling and a fixed random seed.

These changes produce consistent run-to-run results with the same number of threads > 1, and random seed.

Performance seems to be similar, I have not closely profiled memory use but CPU time is similar.

PR Checklist

• No API changes were made (or the changes have been approved)
• No major design changes were made (or the changes have been approved)
• Added test (or behavior not changed)
• Updated API documentation (or API not changed)
• Added license to new files (if any)
• Added Python wrapping to new files (if any) as described in ITK Software Guide Section 9.5
• Added ITK examples for all new major features (if any)

Refer to the ITK Software Guide for
further development details if necessary.

[cookpa]

@ntustison I was curious what Codex would come up with for resolving the non-determinism in Mattes metrics with multiple threads. I tried building ANTs against this as system ITK and it seems to work for ordinary brain affine registration at least.

The ITK tests run on my machine but I started a PR just to see if that reproduces, or if any of the experts here can see an obvious flaw I'm missing.

[hjmjohnson]

@cookpa As we review this, can you do a performance analysis (perhaps using the PerformanceBenchmarks), to get a handle on how much of a performance penalty this will introduce?

[cookpa]

Great idea @hjmjohnson , can I follow the instructions in the README here?

https://github.com/InsightSoftwareConsortium/ITKPerformanceBenchmarking

[dzenanz]

Looks good on a glance.

Interestingly, despite adding a new test, this PR substantially reduced amount of code 😄

[blowekamp]

We recently had a long discussion about the performance of small images with registration when there were a large number of threads:
https://discourse.itk.org/t/8x-slower-registration-with-itk-elastix-python-api-vs-elastix-cli-minimal-reproducible-example/7736/38

The underlying cause was that with a large number of work units or threads, the single threaded histogram merged became a time bottleneck.

Interesting we have an experimental PR in Elastix which goes from the single threaded method to a multi-threaded method here: SuperElastix/elastix#1451

What is needed here a a parallel concurrent "reduce" algorithm to combine all the join PDFs together in a deterministic order. This can be done we a binary tree-like merge up a heap data structure. There are other methods in this class of algorithms. I think an idea solution would be have some type of "Reduce" classes/interface in ITK to provide difference algorithm options.

I may be able to find some time to look into the generic interface this week.

[cookpa]

I built the Benchmarking module but couldn't see any use of Mattes in its code, so I made a small test script using the 2D data from the discourse link above

Results here show real time elapsed, final metric differences, and differences under changes in the number of threads

./affine_mattes_benchmark blobs.tif blobs-rot15deg.tif 10 1,2,4,8 

Summary (mean excludes run 1):

Threads   Mean time (s)   Max metric delta      Max within-thread param delta Param delta vs 1 thread
----------------------------------------------------------------------------------------------------
1         0.347578        0.000000e+00          0.000000e+00                  0.000000e+00
2         0.199482        0.000000e+00          0.000000e+00                  8.483832e-01
4         0.126126        0.000000e+00          0.000000e+00                  4.192933e-01
8         0.126703        0.000000e+00          0.000000e+00                  2.145228e-02

The last column is non-zero because the results still change with different numbers of threads.

Compared to the baseline (from main):

Summary (mean excludes run 1):

Threads   Mean time (s)   Max metric delta      Max within-thread param delta Param delta vs 1 thread
----------------------------------------------------------------------------------------------------
1         0.381328        0.000000e+00          0.000000e+00                  0.000000e+00
2         0.182197        5.764668e-04          1.039299e+00                  8.495317e-01
4         0.163378        6.872430e-04          7.888613e-01                  2.056043e-01
8         0.149558        6.764406e-04          7.425776e-01                  2.193141e-01

test_mattes_perf.zip

This is on a machine with 4 physical cores (8 virtual).

[hjmjohnson]

I'm supportive of this as is to making consistent results. It appears that the performance degredation of forcing an ordering for accumulation appears to be minimal.

[cookpa]

I think this could be ready for review now, but happy to hold off on this if @blowekamp is about to supersede with higher level thread reduction changes.

[greptile-apps]

Greptile Summary

This PR makes multithreaded Mattes metric derivatives deterministic. The main changes are:

• Accumulate derivative PDFs independently for each work unit.
• Merge derivative buffers in ascending work-unit order.
• Remove the mutex-based buffer manager and final thread reduction.
• Add exact repeated-value and repeated-derivative comparisons for the multithreaded affine metric.

Confidence Score: 5/5

This looks safe to merge, with a memory-scaling concern for large global transforms.

The ordered reduction preserves the derivative buffer layout and accumulation semantics. Reused buffers are cleared before each evaluation.

The full derivative image allocated per work unit can substantially increase peak memory.

T-Rex Logs

What T-Rex did

• T-Rex attempted to configure and build the project, but essential tools were not found, including cmake, ctest, and ninja.
• It was also observed that Python had no CMake module, there was no existing build tree, and no ITKMetricsv4TestDriver was located to run tests.
• Because of these missing prerequisites, the configure, build, CTest, and direct execution steps exited with status 127.
• Log artifacts were collected documenting the failure for review.

_{Ran code and verified through T-Rex}

Important Files Changed

Filename	Overview
Modules/Registration/Metricsv4/include/itkMattesMutualInformationImageToImageMetricv4.h	Removes the mutex-backed derivative buffer manager and adds per-work-unit derivative images.
Modules/Registration/Metricsv4/include/itkMattesMutualInformationImageToImageMetricv4.hxx	Removes the obsolete final-thread flush and buffered reduction implementation.
Modules/Registration/Metricsv4/include/itkMattesMutualInformationImageToImageMetricv4GetValueAndDerivativeThreader.hxx	Adds direct per-work-unit accumulation and ordered merging, with increased memory use proportional to the work-unit count.
Modules/Registration/Metricsv4/test/itkMattesMutualInformationImageToImageMetricv4Test.cxx	Adds exact repeated-value and repeated-derivative comparisons for the multithreaded affine metric.

_{Reviews (1): Last reviewed commit: "ENH: Make Mattes metric derivatives dete..." | Re-trigger Greptile}

[greptile-apps]

Dense Buffer Per Work Unit

Each work unit now owns a full parameters × histogram bins² derivative image. Global transforms with many parameters and many work units can therefore multiply the dominant allocation by the thread count and fail allocation even when the previous bounded staging buffers fit in memory.

[cookpa]

The memory use issue pointed out by greptile is concerning, but I'm having trouble seeing it.

I ran an ANTs registration on two 1mm T1w brain volumes

ITK_GLOBAL_DEFAULT_NUMBER_OF_THREADS=8 \
antsRegistration \
  --dimensionality 3 \
  --float 0 \
  --random-seed 12345 \
  --output '[bspline_,bspline_Warped.nii.gz]' \
  --initial-moving-transform '[fixed.nii.gz,moving.nii.gz,1]' \
  --transform 'BSpline[0.1,4x4x4]' \
  --metric 'Mattes[fixed.nii.gz,moving.nii.gz,1,50,Random,0.2]' \
  --convergence '[40x40x20,1e-6,5]' \
  --shrink-factors 8x4x2 \
  --smoothing-sigmas 4x2x1vox \
  --verbose 1

with both the main and PR code. The registration isn't very good, BTW, I don't use BSpline much and I didn't do any real affine initialization first, I was just testing memory use. But the PR appears to use less memory (about 13Gb vs 19Gb), perhaps because of the lack of dynamic buffering?