BSplineSyN registration non-deterministic with sampling strategy != NONE

[cookpa]

Describe the problem

If a sampling strategy is used, BSplineSyN registrations are non-deterministic even when single-threaded and with a fixed seed.

To Reproduce

Using data from templateflow:

export ITK_GLOBAL_DEFAULT_NUMBER_OF_THREADS=1
export ANTS_RANDOM_SEED=12345

# data from templateflow
moving=tpl-OASIS30ANTs_res-01_T1w.nii.gz
fixed=tpl-MNI152NLin2009cAsym_res-01_T1w.nii.gz

# consistent warps between runs
for i in 1 2; do 
  antsRegistration \
    -d 3 -r [ ${fixed}, ${moving}, 1 ] \
    -t BSplineSyN[ 0.1, 26, 0 ] \
    -m CC[ ${fixed}, ${moving}, 1, 3 ] \
    -f 5x4x1 -s 3x2x0vox -c 10x5x0 \
    -o bsplineSyN_None_${i}_ --verbose | tee bsplineSyN_None_${i}.txt
done

# variable warps
for i in 1 2; do 
  antsRegistration \
    -d 3 -r [ ${fixed}, ${moving}, 1 ] \
    -t BSplineSyN[ 0.1, 26, 0 ] \
    -m CC[ ${fixed}, ${moving}, 1, 3, Regular, 0.5 ] \
    -f 5x4x1 -s 3x2x0vox -c 10x5x0 \
    -o bsplineSyN_Regular_${i}_ --verbose | tee bsplineSyN_Regular_${i}.txt
done

If I replace the transform with -t GaussianDisplacementField[ 0.1, 3, 0 ], the results are repeatable.

System information (please complete the following information)

• OS: Mac OS
• OS version: 12.6.2
• Type of system: Laptop

ANTs version information

• ANTs code version: v2.4.2.post11-g066bebc
• ANTs installation type: Compiled from source

[ntustison]

Thanks @cookpa . This is actually a two-part bug. The first one is in ANTs and here's the fix. The problem was that when the SyN method was refactored for special treatment, setting the metric sampling strategy and percentage was never included in the refactoring. So anytime SyN was called, "None" was used for the sampling strategy regardless of what the user set. So with this fix, both BSplineSyN and SyN behave similarly which is good in that sampling is now actually employed for both but now neither one are reproducible regardless of single-threading and explicit random seed setting. This latter problem is corrected by changing this line to

m_CurrentRandomSeed = m_RandomSeed = seed;

I'll post an ITK pull request tomorrow.

[ntustison]

InsightSoftwareConsortium/ITK#3847

[ntustison]

Okay, once we update the ITK version, I think everything should be fixed now.

[cookpa]

Thanks for looking at this @ntustison, I'll try this out.

BTW, would you mind making future changes via PR? That will let them be automatically added to the release notes. If desirable, you can bypass running the CI tests by putting "[skip actions]" (no quotes needed) in your commit message and merge right away.

For things like this I often expand the automated description from the PR, but it's helpful to have the prompt so I don't forget to mention them.

[ntustison]

Sure thing---I'll try to remember.

[cookpa]

@ntustison I'm building the patched ITK now and looking at the code. I'm having trouble seeing how the patch affects the BSplineSyN registration in particular. For example, Affine or GaussianDisplacementField both worked before, and it seems they are prepared the same way by the ants registration helper. Does BSplineSyN do anything special with random numbers?

I see that greedy SyN is a special case, which explains why your other commit (to ANTs) was needed.

Thanks

[ntustison]

I actually played around with the GaussianDisplacementField transform as well and I thought I noticed it behaving similarly to BSplineSyN. So I don't know why it was working for you. If you take a look at the function SetMetricSamplePoints, you'll see that m_CurrentRandomSeed is used to set the sample points at perturbation values. However, that value was previously only initialized in the constructor. So when the ANTs random seed was set in itkantsRegistrationHelper, m_RandomSeed was set to e.g. 12345 with m_CurrentRandomSeed having the initial random value set in the constructor.

[cookpa]

I see the initialization in the constructor, but I thought that would be overwritten in GenerateData

https://github.com/InsightSoftwareConsortium/ITK/blob/de142a0dbde44fa918225e574fb4d1e83e626100/Modules/Registration/RegistrationMethodsv4/include/itkImageRegistrationMethodv4.hxx#L828-L829

  // Ensure the same seed is used for each update
  this->m_CurrentRandomSeed = this->m_RandomSeed;

after the above, GenerateData calls InitializeRegistrationAtEachLevel, which in turn calls SetMetricSamplePoints.

At first I thought it was because itkBSplineSyNImageRegistrationMethod overrides InitializeRegistrationAtEachLevel, but the overridden method calls the superclass method

This is what I don't get, if SetMetricSamplePoints was using a different seed at each run, it wouldn't produce the same output for any transform run to run.

I actually played around with the GaussianDisplacementField transform as well and I thought I noticed it behaving similarly to BSplineSyN.

This is strange! I'll re-check my experiments

[ntustison]

Yeah, I noticed that, too, but one other thing I did to run down the issue was to print out m_CurrentRandomSeed in SetMetricSamplePoints after setting the random seed. I could see where the m_RandomSeed was being set with ANTS_RANDOM_SEED but then m_CurrentRandomSeed wouldn't even be close to that specified value (taking into account that it's incremented in the function).

[cookpa]

I still don't get how, but the patch to ImageRegistrationMethodv4<TFixedImage, TMovingImage, TTransform, TVirtualImage, TPointSet>:: MetricSamplingReinitializeSeed(int seed) works, so I guess we can bump to that once it's merged

[ntustison]

Yeah, I don't disagree with you--it's weird. But I'm glad both of those issues are fixed now. Thanks.