Mattes metric fails for very small voxels

[cookpa]

Related to ANTsPy issue ANTsX/ANTsPy#362

@ntustison , I found that antsRegistration and antsMotionCorr will fail with very small voxels. In the attached archive, I used RNifti to edit the pixel spacing of the MNI152NLin6Asym template from 2mm to 0.2mm, 0.02mm, 0.002mm, etc. Once the voxel spacing gets sufficiently small, the registration starts to fail with Mattes, but still works with CC or MSQ.

With mask (I included the mask because antspy's motion correction uses a mask by default)

z="02mm" # works with 0.2mm spacing
# z="0002mm" # fails
antsRegistration -d 3 -m Mattes\[ t1w_${z}.nii.gz, t1w_${z}.nii.gz, 1, 32 \] -f 1x1 -s 0x0 -c 1x1 -t Rigid\[ 0.1 \] -x mask_${z}.nii.gz -o reg/result_${z}_ --verbose

Without the mask, it goes one level smaller before failing

z="0002mm"
# z="00002mm" # fails
antsRegistration -d 3 -m Mattes\[ t1w_${z}.nii.gz, t1w_${z}.nii.gz, 1, 32 \] -f 1x1 -s 0x0 -c 1x1 -t Rigid\[ 0.1 \]  -o reg/result_${z}_ --verbose

Data here: innerSpace.zip

ANTs Version: v2.3.5.post101-g1c5ae7c
Compiled: Mar 23 2022 11:49:51

[cookpa]

This is not a perfect test because the images are identical, and the original issue in antspy is intermittent...but I ran it 100 times with the smallest voxels and Mattes fails consistently, and CC does not fail at all.

[cookpa]

Reproduced on Russ' data:

ImageMath 4 series-tdtomato_size-mm_units-mm_vols/series-tdtomato_size-mm_units-mm_vol_.nii.gz TimeSeriesDisassemble ants_example_data/series-tdtomato_size-mm_units-mm.nii

for i in `seq 1000 1056`; do echo $i; antsRegistration -d 3 -m CC[ ants_example_data/series-tdtomato_size-mm_units-mm_mean.nii , series-tdtomato_size-mm_units-mm_vols/series-tdtomato_size-mm_units-mm_vol_${i}.nii.gz , 1, 2 ] -t Rigid[ 0.1 ] -c 5x0 -f 2x1 -s 1x0vox -o [ pairwise/vol${i}ToMean , pairwise/vol${i}ToMeanDeformed.nii.gz ] -x mask.nii.gz

Pairwise calls to antsRegistration work with CC metric, with no failures (antsMotionCorr, not so much, will post more in the other thread). Mattes errors out consistently. This is with dense sampling so not a random initialization thing.

A sample of the error, there's a whole bunch of warnings like this

WARNING: In /Users/pcook/tmp/NOT_BACKED_UP/antsMaster/build/staging/include/ITK-5.3/itkImageBase.hxx, line 87
Image (0x7feda965d350): Negative spacing is not supported and may result in undefined behavior.
Proceeding to set spacing to [-inf, 0.0356921]

/Users/pcook/tmp/NOT_BACKED_UP/antsMaster/build/ITKv5/Modules/ThirdParty/VNL/src/vxl/core/vnl/algo/vnl_svd.hxx: suspicious return value (2) from SVDC
/Users/pcook/tmp/NOT_BACKED_UP/antsMaster/build/ITKv5/Modules/ThirdParty/VNL/src/vxl/core/vnl/algo/vnl_svd.hxx: M is 2x2
M = [ ...
            -inf                0 
             nan  0.0356921405224  ]

Then

Exception caught: 
itk::ExceptionObject (0x7feda9407a00)
Location: "unknown" 
File: /Users/pcook/tmp/NOT_BACKED_UP/antsMaster/build/staging/include/ITK-5.3/itkMattesMutualInformationImageToImageMetricv4.hxx
Line: 311
Description: ITK ERROR: MattesMutualInformationImageToImageMetricv4(0x7feda9657ae0): All samples map outside moving image buffer. The images do not sufficiently overlap. They need to be initialized to have more overlap before this metric will work. For instance, you can align the image centers by translation.

Wait a minute - M is 2x2 - I remember this #756

[cookpa]

And without a mask, it runs without any error:

ImageMath 4 series-tdtomato_size-mm_units-mm_vols/series-tdtomato_size-mm_units-mm_vol_.nii.gz \
  TimeSeriesDisassemble ants_example_data/series-tdtomato_size-mm_units-mm.nii

for i in `seq 1000 1056`; do 
  antsRegistration -d 3 
  -m Mattes[ ants_example_data/series-tdtomato_size-mm_units-mm_mean.nii , \
                     series-tdtomato_size-mm_units-mm_vols/series-tdtomato_size-mm_units-mm_vol_${i}.nii.gz , 1, 32 ] \
  -t Rigid[ 0.1 ] -c 5x0 -f 2x1 -s 1x0vox \
  -o [ pairwise/vol${i}ToMean , pairwise/vol${i}ToMeanDeformed.nii.gz ]
done

[ntustison]

As I mentioned before, I think it's the scale parameter estimator for the optimizer. When the voxel size gets small, the perturbation blows up transforming the image domain outside of the domain. Still trying to figure out the source of this behavior.

[ntustison]

Okay, so it's the Mattes gradient that is used to estimate the scale parameter which blows up. Note that the measurement seems stable across voxel size. Here's the gradient for the different image pairs:

02
gradient =  [-0.00007291265663351622, 0.00004265629618109531, -0.00004060390151223091, -0.002249708770099317, -0.01573422042220483, -0.0051309019203129805]
step size = 0.0306068
 
002
gradient = [-0.0008968057751967737, 0.0013755250886906616, 0.0003417147076551544, -0.1020198617887443, -0.0506777312583488, -0.12274299980452325]
step size = 0.442716

0002
gradient = [0.008129704388570034, -0.0012587711492438214, 0.007701400336923561, -0.20802612665886497, 1.7621225187161993, 0.35251437736267455]
step size = 3.55652

00002
gradient = [0.22635166283827757, -0.23517012913934726, 0.020511407844790164, -13.876707745938255, -13.963464469888887, 45.81445385262064]
step size = 95.5934

So I think this should be brought up over at the ITK forum regardless of what we do here. Perhaps @hjmjohnson or @thewtex could weigh in as they have more experience with the ITK development of Mattes (if I'm remembering correctly).
Also, perhaps this issue has already been brought up which they would know.

One thing we could do here is expose the gradient filter option within all the metrics to the ANTs user which is currently hard-coded as false (for time considerations) although, as noted, that can have a significant impact on performance. For example, when gradientfilter = true, the metric gradient seems fine:

02
gradient = [6.812373208743657e-7, 0.000009766801576907105, 0.000010615970500715887, -0.0014030361562300047, 0.0012085649563486683, -0.0005154792818747656]
step size = 0.00429771
  

002
gradient = [0.000036171825083973186, -0.00003891936402654846, 0.0000019300638316771268, 0.000897978938187557, 0.0017966492291090058, 0.005073424707967914]
step size = 0.014331

0002
gradient = [0.000020454224819336065, -0.000007224404640241085, 0.00001545085394110889, -0.001058836622265709, 0.0026392894829883554, 0.0019142092628396242]
step size = 0.00793043

00002
gradient = [-0.0000015401530948162888, -0.0000273019305472448, -0.000028163501132950842, 0.004592840495396785, -0.0010733920383553086, 0.0003557576159049231]
step size = 0.0110835

[cookpa]

Thanks for the additional option @ntustison, I'll check it out.

I'll update the Wiki when I get a chance to document how to deal with these kinds of errors with gradient filtering / workaround with unit change as suggested in ANTsX/ANTsPy#362 .

I had a quick look around the ITK forum and didn't see this exact issue. I might try to make a minimal ITK example to highlight the issue.

It would be nice to add gradient filtering option to antsMotionCorr and see if it improves stability. But I got pretty decent results on my tests after fixing the smoothing units.

[ntustison]

Sure, give me a sec and I'll make a pull request.

[gdevenyi]

So I think this should be brought up over at the ITK forum regardless of what we do here. Perhaps @hjmjohnson or @thewtex could weigh in as they have more experience with the ITK development of Mattes (if I'm remembering correctly).

Is there a forum thread I can follow?

[cookpa]

I'll cross post when I get to it, but I haven't had time to make an example that only depends on ITK, not including all of ANTs. Maybe in a week or so. If you or anyone would care to attempt it, I was going to base off this

https://examples.itk.org/src/registration/common/perform2dtranslationregistrationwithmeansquares/documentation

and see if I can reproduce the metric gradient problem.

[gdevenyi]

Following up here, I can't tell if antsRegistration needs this adjustment as well? It only got added to antsMotionCorr.

[cookpa]

It does affect antsRegistration as well. I haven't been able to make a minimal ITK example using the v4 registration framework with the scales estimators.

[gdevenyi]

I've now done some very rough tests using my custom registration script I use for everything, and a typical registration (multiple scales of rigid/similarity/linear with Mattes + SyNCC) from a 1mm T1 to the MNI 09c model has zero effect on computational time (~10sec faster on 13min with 40 cores)

[cookpa]

I'm thinking to close this now as it seems using the gradient filter fixes these issues? Is there still a need to investigate at the ITK level?

[ntustison]

Sounds good to me.