Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

antsRegistration does not produce equivalent results #210

Open
dmirman-zz opened this issue Mar 29, 2018 · 49 comments
Open

antsRegistration does not produce equivalent results #210

dmirman-zz opened this issue Mar 29, 2018 · 49 comments

Comments

@dmirman-zz
Copy link

Multiple runs of antsRegistration() on the same images produces different results. Here is a simple example:

library(ANTsR)

template <- antsImageRead(getANTsRData("ch2")) #colins template
atlas <- antsImageRead(getANTsRData("mnia")) #aal atlas

#register template space to atlas
regRigid <- list()
regRigidcount <- NA
for(i in 1:10){
  regRigid[[i]] <- antsRegistration(fixed = atlas, moving = template, typeofTransform = c("Rigid"))
  regRigidcount[[i]] <- sum(regRigid[[i]]$warpedmovout)
  print(regRigidcount[[i]])
}
# the counts are not the same for the 10 runs of the registration

Similarly non-identical results show up for other transforms ("Affine", "SyNCC"). Although the differences seem somewhat small, it is not clear why they are different at all. More importantly, these small differences have substantive consequences down-stream: we ran a lesion-symptom mapping analysis and wanted to describe the results in terms of the AAL atlas, but the 10 somewhat different transformations produce somewhat different descriptions:

#read in data to cluster
# file is available from http://drive.google.com/file/d/1tm2PI6yxhEngB60jNjoP6YQ3T_H7u8JM"
data <- antsImageRead("LSM_result.nii") 

#Rigid Transformation
wrpRigid <- list()
clustRigid <- list()
statRigid <- list()
for(a in 1:length(regRigid)){
  wrpRigid[[a]] <- antsApplyTransforms(fixed = atlas, moving = data, 
                                      transformlist = regRigid[[a]]$fwdtransforms, interpolator = c("NearestNeighbor"))
  clustRigid[[a]] <- labelClusters(wrpRigid[[a]], maxThresh = 1, minClusterSize = 0, minThresh = 0.1)
  statRigid[[a]] <- labelStats(wrpRigid[[a]], clustRigid[[a]])
  print(statRigid[[a]])
}

In my test, I get anywhere from 7 to 11 clusters, with slightly different masses and (x,y,z) coordinates.

This was run most recently using R version 3.3.1, ANTsR version: 0.5.6, ANTsRCore version 0.3.4, but we've replicated this problem on two different Linux machines.

@ntustison
Copy link
Member

Please see this thread where these issues are discussed.

@stnava
Copy link
Member

stnava commented Mar 29, 2018 via email

@muschellij2
Copy link
Collaborator

muschellij2 commented Mar 29, 2018

It seems like that thread (ANTsX/ANTs#444 (comment)) indicated changing --float to --double may help, is that correct?

If so, what about changing
https://github.com/ANTsX/ANTsRCore/blob/master/R/antsRegistration.R#L494
to --double.

@dmirman-zz
Copy link
Author

That was also my understanding of that thread, but it doesn't seem like antsRegistration will accept a --double flag. When I try
regRigid[[i]] <- antsRegistration(fixed = atlas, moving = template, typeofTransform = c("Rigid"), "--double")
I get
ERROR: Invalid flag provided double

@stnava
Copy link
Member

stnava commented Mar 29, 2018 via email

@dmirman-zz
Copy link
Author

I'd like to test running this in single threaded mode without random sampling, but how do I specify that within ANTsR?

@cookpa
Copy link
Member

cookpa commented Mar 29, 2018

I think you would need dense sampling (strategy "NONE") to have completely reproducible results. With regular sampling, a perturbation is applied to the points to randomize their location within the voxel grid.

https://github.com/InsightSoftwareConsortium/ITK/blob/8a2a15f41218c925c0a89119e09419d48f83eb22/Modules/Registration/RegistrationMethodsv4/include/itkImageRegistrationMethodv4.hxx#L954-L984

@dmirman-zz
Copy link
Author

Please pardon my ignorance, but I do not know enough about how ANTsR communicates with ITK (or about how ITK works) in order to implement the solutions that have been described here. As I understand it, two strategies have been proposed:
(1) use --double precision, but antsRegistration doesn't seem to take that as a flag and there may be overly negative consequences for memory use.
(2) run in single threaded mode without random sampling and use dense sampling. It's not clear to me how to control those parameters from within ANTsR.

Any advice about implementing either or both of these would be most appreciated.

@ntustison
Copy link
Member

antsRegistration, as defined in ANTsR, is used more like a script which simplifies the underlying antsRegistration tool in ANTs. If you need to perform your experiments within ANTsR, you might want to access the underlying program using the antsRegistration( list_of_arguments ) format. For example, to see the short and long help menus, you can type

> antsRegistration( list( '-h', 1 ) )
> antsRegistration( list( '--help', 1 ) )

You could then use double by specifying "--float 0", i.e.,

> args <- list( "-d", ... "--float", "0", ... )
> antsRegistration( args )

See here for an example. In the same file you can see how all the different transforms are defined.

@dorianps
Copy link
Collaborator

dorianps commented Mar 30, 2018 via email

@stnava
Copy link
Member

stnava commented Mar 30, 2018 via email

@ntustison
Copy link
Member

Thanks Brian--informative as always. I added this wiki page as an initial attempt for addressing the relevant issues. Let me know if you'd like something different.

@dmirman-zz
Copy link
Author

I can confirm that using double precision is computationally impractical: I tried it for a Rigid registration and it still hadn't finished after several hours (float version takes ~20sec). Thanks everyone for an informative discussion.

@stnava
Copy link
Member

stnava commented Apr 2, 2018 via email

@stnava
Copy link
Member

stnava commented Jun 3, 2018

we can close this issue due to several recent changes in ANTs and ANTsR. you will need to define (in your environment) ANTS_RANDOM_SEED=1 (or whatever your favorite number is). e.g. my .Renviron file reads:

ANTS_RANDOM_SEED=1234
ITK_GLOBAL_DEFAULT_NUMBER_OF_THREADS=1

this should give you repeatable results.

ultimately, we may also allow a command line parameter to control the seed.

one note: histogram matching induced an additional unexpected source of random variability.

@ncullen93 may find this of interest as well.

@moskvich412
Copy link

In testing ANTs registration, we discovered it is not reproducible: running several times on the same inputs gives different morphisms. Implementing:

ANTS_RANDOM_SEED=1234
ITK_GLOBAL_DEFAULT_NUMBER_OF_THREADS=1

as suggested above, does not seem to completely resolve the issue. I see based on CPU usage that threads are not created. But it is impossible for me to check if fixing seed took effect. We are using version 2.3.1, should these variables have effect in this version?

Are there other methods to enforce reproducibility?

P.S. It is not obvious to me how randomness can help solving for brain morphism (random seed suggests it is done on purpose), but this is not the question at hand.

@ntustison
Copy link
Member

@moskvich412 I see that you have read a couple of the comments in this thread, specifically the last comment from @stnava dealing with ANTS_RANDOM_SEED and ITK_GLOBAL_DEFAULT_NUMBER_OF_THREADS but have you read some of the other comments pointing to some of the foundational problems with this issue that goes beyond ANTs? There are even some references pointing to discussion of these algorithmic stochastic elements and why they're important, including in the study of brain morphisms.

@moskvich412
Copy link

moskvich412 commented Sep 28, 2019 via email

@ntustison
Copy link
Member

You've probably done all you can to remove ANTs variability.

@cookpa
Copy link
Member

cookpa commented Sep 28, 2019 via email

@moskvich412
Copy link

moskvich412 commented Sep 28, 2019 via email

@cookpa
Copy link
Member

cookpa commented Sep 28, 2019 via email

@dorianps
Copy link
Collaborator

dorianps commented Sep 28, 2019 via email

@muschellij2
Copy link
Collaborator

muschellij2 commented Sep 28, 2019 via email

@moskvich412
Copy link

I use this bash script to compute totalTransform.nii.gz:

export ANTS_RANDOM_SEED=10
export ITK_GLOBAL_DEFAULT_NUMBER_OF_THREADS=1

N3BiasFieldCorrection 3 input.nii.gz input.nii.gz 4
N3BiasFieldCorrection 3 template.nii.gz template.nii.gz 4

ANTS 3 -m CC[template.nii.gz,input.nii.gz,1,4] -t SyN[0.25]
-r Gauss[3,0] -o prefix -i 100x90x30 --use-Histogram-Matching
--number-of-affine-iterations 10000x10000x10000x10000x10000 --MI-option 32x16000

ComposeMultiTransform 3 totalTransform.nii.gz -R template.nii.gz prefixWarp.nii.gz prefixAffine.txt

The output morphism, totalTransform.nii.gz, is then applied using:

WarpImageMultiTransform 3 input.nii.gz morphedInput.nii.gz totalTransform.nii.gz -R totalTransform.nii.gz

I see that both N3BiasFieldCorrection and ANTS are running in a single thread. We are using version 2.3.1.

Comments on the choice of metrics, iterations, etc are also welcome. Maybe there is a better set.

@stnava
Copy link
Member

stnava commented Sep 28, 2019

See an example here ... looks like Mattes is not reproducible. Comments welcome.

It would be good to add Python and bash examples to above. Please contribute via pull request.

We don't recommend using the program called ANTS which is not maintained except perhaps to "get it to compile." The current tool is called antsRegistration.

@stnava stnava reopened this Sep 28, 2019
@ntustison
Copy link
Member

@muschellij2 Can you post an input image where this occurs? There's no explicitly random element in N3BiasFieldCorrection so I'm curious as to where this is happening.

@muschellij2
Copy link
Collaborator

muschellij2 commented Oct 7, 2019 via email

@ntustison
Copy link
Member

Okay, thanks for the clarification. I was worried since N4 has gone through changes over in the ITK GitHub repo, some of which concern overlapping elements with N3 and I wanted to make sure that I wasn't behind in my upkeep.

@moskvich412
Copy link

I am sorry for not being clear. The set of commands I wrote were internal parts of my script, which first makes local copies of the input and template. Bias correction is not iteratively applied. I played with it separately and each iteration does change the image. There is no way to terminate this process because there is no good way to distinguish receive and transmit coils and T1 weighting in the context of MRI from the actual spatial variation in the image.

I also realized I was under influence of the older, version 1.9, manual which talked about the routines affected by randomness. I am glad the cause was my use of the deprecated commands which I will now fix and get reproducible results.

kousu added a commit to spinalcordtoolbox/spinalcordtoolbox that referenced this issue Apr 25, 2020
Upon upgrading our ANTS build we found that motion correction was no
longer deterministic (#2642 (comment)).

The ANTS maintainers suggested setting samplingStrategy=None
(ANTsX/ANTs#976 (comment))
to get deterministic output, but they also warned that this will induce
a variance-for-bias tradeoff:

* https://github.com/ANTsX/ANTs/wiki/antsRegistration-reproducibility-issues
* ANTsX/ANTsR#210 (comment)
* Thévenaz, P., Bierlaire, M., & Unser, M. (2008).
  Halt on sampling for image registration based on mutual information.
  Sampling Theory in Signal and Image Processing, 7(2).
  <http://bigwww.epfl.ch/preprints/thevenaz0602p.pdf>

Because we are hard-coding the sampling rate I removed it from ParamMoco
entirely. If we want to support it properly we would have to pass
.sampling_strategy and .sampling_percentage separately.
kousu added a commit to spinalcordtoolbox/spinalcordtoolbox that referenced this issue Apr 25, 2020
Upon upgrading our ANTS build we found that motion correction was no
longer deterministic (#2642 (comment)).

The ANTS maintainers suggested setting samplingStrategy=None
(ANTsX/ANTs#976 (comment))
to get deterministic output, but they also warned that this will induce
a variance-for-bias tradeoff:

* https://github.com/ANTsX/ANTs/wiki/antsRegistration-reproducibility-issues
* ANTsX/ANTsR#210 (comment)
* Thévenaz, P., Bierlaire, M., & Unser, M. (2008).
  Halt on sampling for image registration based on mutual information.
  Sampling Theory in Signal and Image Processing, 7(2).
  <http://bigwww.epfl.ch/preprints/thevenaz0602p.pdf>

Because we are hard-coding the sampling rate I removed it from ParamMoco
entirely. If we want to support it properly we would have to pass
.sampling_strategy and .sampling_percentage separately.
kousu added a commit to spinalcordtoolbox/spinalcordtoolbox that referenced this issue Apr 25, 2020
Upon upgrading our ANTS build we found that motion correction was no
longer deterministic (#2642 (comment)).

The ANTS maintainers suggested setting samplingStrategy=None
(ANTsX/ANTs#976 (comment))
to get deterministic output, but they also warned that this will induce
a variance-for-bias tradeoff:

* https://github.com/ANTsX/ANTs/wiki/antsRegistration-reproducibility-issues
* ANTsX/ANTsR#210 (comment)
* Thévenaz, P., Bierlaire, M., & Unser, M. (2008).
  Halt on sampling for image registration based on mutual information.
  Sampling Theory in Signal and Image Processing, 7(2).
  <http://bigwww.epfl.ch/preprints/thevenaz0602p.pdf>

Because we are hard-coding the sampling rate I removed it from ParamMoco
entirely. If we want to support it properly we would have to pass
.sampling_strategy and .sampling_percentage separately.
kousu added a commit to spinalcordtoolbox/spinalcordtoolbox that referenced this issue May 4, 2020
This makes the output deterministic, at the cost of running unnecessarily slowly.

The order of the floating point sums used internally is numerically unstable.

See

https://github.com/ANTsX/ANTs/wiki/antsRegistration-reproducibility-issues
ANTsX/ANTs#444 (comment)
ANTsX/ANTsR#210 (comment)
kousu added a commit to spinalcordtoolbox/spinalcordtoolbox that referenced this issue May 4, 2020
Upon upgrading our ANTS build we found that motion correction was no
longer deterministic (#2642 (comment)).

The ANTS maintainers suggested setting samplingStrategy=None
(ANTsX/ANTs#976 (comment))
to get deterministic output, but they also warned that this will induce
a variance-for-bias tradeoff:

* https://github.com/ANTsX/ANTs/wiki/antsRegistration-reproducibility-issues
* ANTsX/ANTsR#210 (comment)
* Thévenaz, P., Bierlaire, M., & Unser, M. (2008).
  Halt on sampling for image registration based on mutual information.
  Sampling Theory in Signal and Image Processing, 7(2).
  <http://bigwww.epfl.ch/preprints/thevenaz0602p.pdf>

Because we are hard-coding the sampling rate I removed it from ParamMoco
entirely. If we want to support it properly we would have to pass
.sampling_strategy and .sampling_percentage separately.
kousu added a commit to spinalcordtoolbox/spinalcordtoolbox that referenced this issue May 4, 2020
This is also to make the output deterministic, at the cost of running slow.

It turned out that using dense sampling wasn't enough; there was still some
numerical instability that came from the order of addition:

* https://github.com/ANTsX/ANTs/wiki/antsRegistration-reproducibility-issues#variance-due-to-floating-point-precision-errors

For some reason it only appeared on OS X, and only about 10% of the time, and never on Linux.
I [showed](#2642 (comment))
that the instability in isct_antsSliceRegularizedRegistration did exist on Linux, so something
still unknown about how we call it was hiding it there.

See:

* https://github.com/ANTsX/ANTs/wiki/antsRegistration-reproducibility-issues
* ANTsX/ANTs#444 (comment)
* ANTsX/ANTsR#210 (comment)
kousu added a commit to spinalcordtoolbox/spinalcordtoolbox that referenced this issue May 4, 2020
This is also to make the output deterministic, at the cost of running slow.

It turned out that using dense sampling wasn't enough; there was still some
numerical instability that came from the order of addition:

* https://github.com/ANTsX/ANTs/wiki/antsRegistration-reproducibility-issues#variance-due-to-floating-point-precision-errors

For some reason it only appeared on OS X, and only about 10% of the time, and never on Linux.
I [showed](#2642 (comment))
that the instability in isct_antsSliceRegularizedRegistration did exist on Linux, so something
still unknown about how we call it was hiding it there.

This was actually supposed to be in place already but the code had atrophied,
so all this does is fix it up.

See:

* https://github.com/ANTsX/ANTs/wiki/antsRegistration-reproducibility-issues
* ANTsX/ANTs#444 (comment)
* ANTsX/ANTsR#210 (comment)
kousu added a commit to spinalcordtoolbox/spinalcordtoolbox that referenced this issue May 9, 2020
Upon upgrading our ANTS build we found that motion correction was no
longer deterministic (#2642 (comment)).

The ANTS maintainers suggested setting samplingStrategy=None
(ANTsX/ANTs#976 (comment))
to get deterministic output, but they also warned that this will induce
a variance-for-bias tradeoff:

* https://github.com/ANTsX/ANTs/wiki/antsRegistration-reproducibility-issues
* ANTsX/ANTsR#210 (comment)
* Thévenaz, P., Bierlaire, M., & Unser, M. (2008).
  Halt on sampling for image registration based on mutual information.
  Sampling Theory in Signal and Image Processing, 7(2).
  <http://bigwww.epfl.ch/preprints/thevenaz0602p.pdf>

Because we are hard-coding the sampling rate I removed it from ParamMoco
entirely. If we want to support it properly we would have to pass
.sampling_strategy and .sampling_percentage separately.
kousu added a commit to spinalcordtoolbox/spinalcordtoolbox that referenced this issue May 9, 2020
This is also to make the output deterministic, at the cost of running slow.

It turned out that using dense sampling wasn't enough; there was still some
numerical instability that came from the order of addition:

* https://github.com/ANTsX/ANTs/wiki/antsRegistration-reproducibility-issues#variance-due-to-floating-point-precision-errors

For some reason it only appeared on OS X, and only about 10% of the time, and never on Linux.
I [showed](#2642 (comment))
that the instability in isct_antsSliceRegularizedRegistration did exist on Linux, so something
still unknown about how we call it was hiding it there.

This was actually supposed to be in place already but the code had atrophied,
so all this does is fix it up.

See:

* https://github.com/ANTsX/ANTs/wiki/antsRegistration-reproducibility-issues
* ANTsX/ANTs#444 (comment)
* ANTsX/ANTsR#210 (comment)
kousu added a commit to spinalcordtoolbox/spinalcordtoolbox that referenced this issue May 11, 2020
Upon upgrading our ANTS build we found that motion correction was no
longer deterministic (#2642 (comment)).

The ANTS maintainers suggested setting samplingStrategy=None
(ANTsX/ANTs#976 (comment))
to get deterministic output, but they also warned that this will induce
a variance-for-bias tradeoff:

* https://github.com/ANTsX/ANTs/wiki/antsRegistration-reproducibility-issues
* ANTsX/ANTsR#210 (comment)
* Thévenaz, P., Bierlaire, M., & Unser, M. (2008).
  Halt on sampling for image registration based on mutual information.
  Sampling Theory in Signal and Image Processing, 7(2).
  <http://bigwww.epfl.ch/preprints/thevenaz0602p.pdf>

Because we are hard-coding the sampling rate I removed it from ParamMoco
entirely. If we want to support it properly we would have to pass
.sampling_strategy and .sampling_percentage separately.
kousu added a commit to spinalcordtoolbox/spinalcordtoolbox that referenced this issue May 11, 2020
This is also to make the output deterministic, at the cost of running slow.

It turned out that using dense sampling wasn't enough; there was still some
numerical instability that came from the order of addition:

* https://github.com/ANTsX/ANTs/wiki/antsRegistration-reproducibility-issues#variance-due-to-floating-point-precision-errors

For some reason it only appeared on OS X, and only about 10% of the time, and never on Linux.
I [showed](#2642 (comment))
that the instability in isct_antsSliceRegularizedRegistration did exist on Linux, so something
still unknown about how we call it was hiding it there.

This was actually supposed to be in place already but the code had atrophied,
so all this does is fix it up.

See:

* https://github.com/ANTsX/ANTs/wiki/antsRegistration-reproducibility-issues
* ANTsX/ANTs#444 (comment)
* ANTsX/ANTsR#210 (comment)
kousu added a commit to spinalcordtoolbox/spinalcordtoolbox that referenced this issue May 17, 2020
Upon upgrading our ANTS build we found that motion correction was no
longer deterministic (#2642 (comment)).

The ANTS maintainers suggested setting samplingStrategy=None
(ANTsX/ANTs#976 (comment))
to get deterministic output, but they also warned that this will induce
a variance-for-bias tradeoff:

* https://github.com/ANTsX/ANTs/wiki/antsRegistration-reproducibility-issues
* ANTsX/ANTsR#210 (comment)
* Thévenaz, P., Bierlaire, M., & Unser, M. (2008).
  Halt on sampling for image registration based on mutual information.
  Sampling Theory in Signal and Image Processing, 7(2).
  <http://bigwww.epfl.ch/preprints/thevenaz0602p.pdf>

Because we are hard-coding the sampling rate I removed it from ParamMoco
entirely. If we want to support it properly we would have to pass
.sampling_strategy and .sampling_percentage separately.
kousu added a commit to spinalcordtoolbox/spinalcordtoolbox that referenced this issue May 17, 2020
This is also to make the output deterministic, at the cost of running slow.

It turned out that using dense sampling wasn't enough; there was still some
numerical instability that came from the order of addition:

* https://github.com/ANTsX/ANTs/wiki/antsRegistration-reproducibility-issues#variance-due-to-floating-point-precision-errors

For some reason it only appeared on OS X, and only about 10% of the time, and never on Linux.
I [showed](#2642 (comment))
that the instability in isct_antsSliceRegularizedRegistration did exist on Linux, so something
still unknown about how we call it was hiding it there.

This was actually supposed to be in place already but the code had atrophied,
so all this does is fix it up.

See:

* https://github.com/ANTsX/ANTs/wiki/antsRegistration-reproducibility-issues
* ANTsX/ANTs#444 (comment)
* ANTsX/ANTsR#210 (comment)
kousu added a commit to spinalcordtoolbox/spinalcordtoolbox that referenced this issue May 17, 2020
Upon upgrading our ANTS build we found that motion correction was no
longer deterministic (#2642 (comment)).

The ANTS maintainers suggested setting samplingStrategy=None
(ANTsX/ANTs#976 (comment))
to get deterministic output, but they also warned that this will induce
a variance-for-bias tradeoff:

* https://github.com/ANTsX/ANTs/wiki/antsRegistration-reproducibility-issues
* ANTsX/ANTsR#210 (comment)
* Thévenaz, P., Bierlaire, M., & Unser, M. (2008).
  Halt on sampling for image registration based on mutual information.
  Sampling Theory in Signal and Image Processing, 7(2).
  <http://bigwww.epfl.ch/preprints/thevenaz0602p.pdf>

I think ideally moco.py would thinly wrap antsRegistration, and
expose both .sampling_strategy and .sampling_percentage directly,
but until now we've always hardcoded 'samplingStrategy=Regular'
so in the interest of not rocking the boat too much all I did
was add a 'sampling=None' case, *and set it as the default*

Updated version re feedback: #2642 (comment)
jcohenadad pushed a commit to spinalcordtoolbox/spinalcordtoolbox that referenced this issue May 22, 2020
Upon upgrading our ANTS build we found that motion correction was no
longer deterministic (#2642 (comment)).

The ANTS maintainers suggested setting samplingStrategy=None
(ANTsX/ANTs#976 (comment))
to get deterministic output, but they also warned that this will induce
a variance-for-bias tradeoff:

* https://github.com/ANTsX/ANTs/wiki/antsRegistration-reproducibility-issues
* ANTsX/ANTsR#210 (comment)
* Thévenaz, P., Bierlaire, M., & Unser, M. (2008).
  Halt on sampling for image registration based on mutual information.
  Sampling Theory in Signal and Image Processing, 7(2).
  <http://bigwww.epfl.ch/preprints/thevenaz0602p.pdf>

Because we are hard-coding the sampling rate I removed it from ParamMoco
entirely. If we want to support it properly we would have to pass
.sampling_strategy and .sampling_percentage separately.
jcohenadad pushed a commit to spinalcordtoolbox/spinalcordtoolbox that referenced this issue May 22, 2020
This is also to make the output deterministic, at the cost of running slow.

It turned out that using dense sampling wasn't enough; there was still some
numerical instability that came from the order of addition:

* https://github.com/ANTsX/ANTs/wiki/antsRegistration-reproducibility-issues#variance-due-to-floating-point-precision-errors

For some reason it only appeared on OS X, and only about 10% of the time, and never on Linux.
I [showed](#2642 (comment))
that the instability in isct_antsSliceRegularizedRegistration did exist on Linux, so something
still unknown about how we call it was hiding it there.

This was actually supposed to be in place already but the code had atrophied,
so all this does is fix it up.

See:

* https://github.com/ANTsX/ANTs/wiki/antsRegistration-reproducibility-issues
* ANTsX/ANTs#444 (comment)
* ANTsX/ANTsR#210 (comment)
jcohenadad pushed a commit to spinalcordtoolbox/spinalcordtoolbox that referenced this issue May 22, 2020
Upon upgrading our ANTS build we found that motion correction was no
longer deterministic (#2642 (comment)).

The ANTS maintainers suggested setting samplingStrategy=None
(ANTsX/ANTs#976 (comment))
to get deterministic output, but they also warned that this will induce
a variance-for-bias tradeoff:

* https://github.com/ANTsX/ANTs/wiki/antsRegistration-reproducibility-issues
* ANTsX/ANTsR#210 (comment)
* Thévenaz, P., Bierlaire, M., & Unser, M. (2008).
  Halt on sampling for image registration based on mutual information.
  Sampling Theory in Signal and Image Processing, 7(2).
  <http://bigwww.epfl.ch/preprints/thevenaz0602p.pdf>

I think ideally moco.py would thinly wrap antsRegistration, and
expose both .sampling_strategy and .sampling_percentage directly,
but until now we've always hardcoded 'samplingStrategy=Regular'
so in the interest of not rocking the boat too much all I did
was add a 'sampling=None' case, *and set it as the default*

Updated version re feedback: #2642 (comment)
@muschellij2
Copy link
Collaborator

So revisiting this a bit again, I think there is a sneaky bug somewhere for this. I don't believe this is fixed.

Here are 3 different scenarios, setting environment variables before loading (like .Renviron would do), setting it after, and then setting it after some ANTsR operation has been completed.

Setting the Seed before ANTsR Loaded - Works

Sys.setenv(ANTS_RANDOM_SEED = 123)
Sys.setenv(ITK_GLOBAL_DEFAULT_NUMBER_OF_THREADS = 1)
library( ANTsR )
#> Loading required package: ANTsRCore
#> 
#> Attaching package: 'ANTsRCore'
#> The following objects are masked from 'package:stats':
#> 
#>     sd, var
#> The following objects are masked from 'package:base':
#> 
#>     all, any, apply, max, min, prod, range, sum

fi <- antsImageRead(getANTsRData("r16") )
mi <- antsImageRead(getANTsRData("r64") )

tx = "SyN"
print(Sys.getenv("ANTS_RANDOM_SEED"))
#> [1] "123"
mytx1 <- antsRegistration(fixed=fi, mi, typeofTransform = tx)
mytx2 <- antsRegistration(fixed=fi, mi, typeofTransform = tx)
print(sum(mytx1$warpedmovout - mytx2$warpedmovout))
#> [1] 0
max(abs(as.array(mytx1$warpedmovout) - as.array(mytx2$warpedmovout)))
#> [1] 0

Created on 2020-08-06 by the reprex package (v0.3.0)

Session info
devtools::session_info()
#> ─ Session info ───────────────────────────────────────────────────────────────
#>  setting  value                       
#>  version  R version 4.0.0 (2020-04-24)
#>  os       macOS Mojave 10.14.6        
#>  system   x86_64, darwin17.0          
#>  ui       X11                         
#>  language (EN)                        
#>  collate  en_US.UTF-8                 
#>  ctype    en_US.UTF-8                 
#>  tz       America/New_York            
#>  date     2020-08-06                  
#> 
#> ─ Packages ───────────────────────────────────────────────────────────────────
#>  package     * version    date       lib source                                
#>  ANTsR       * 0.5.6.1    2020-06-01 [1] Github (ANTsX/ANTsR@9c7c9b7)          
#>  ANTsRCore   * 0.7.4.6    2020-07-07 [1] Github (muschellij2/ANTsRCore@61c37a1)
#>  assertthat    0.2.1      2019-03-21 [1] CRAN (R 4.0.0)                        
#>  backports     1.1.8      2020-06-17 [1] CRAN (R 4.0.0)                        
#>  callr         3.4.3      2020-03-28 [1] CRAN (R 4.0.0)                        
#>  cli           2.0.2      2020-02-28 [1] CRAN (R 4.0.0)                        
#>  crayon        1.3.4      2017-09-16 [1] CRAN (R 4.0.0)                        
#>  desc          1.2.0      2020-06-01 [1] Github (muschellij2/desc@b0c374f)     
#>  devtools      2.3.0.9000 2020-06-01 [1] Github (hadley/devtools@202ea81)      
#>  digest        0.6.25     2020-02-23 [1] CRAN (R 4.0.0)                        
#>  ellipsis      0.3.1      2020-05-15 [1] CRAN (R 4.0.0)                        
#>  evaluate      0.14       2019-05-28 [1] CRAN (R 4.0.0)                        
#>  fansi         0.4.1      2020-01-08 [1] CRAN (R 4.0.0)                        
#>  fs            1.4.2      2020-06-30 [1] CRAN (R 4.0.0)                        
#>  glue          1.4.1      2020-05-13 [1] CRAN (R 4.0.0)                        
#>  highr         0.8        2019-03-20 [1] CRAN (R 4.0.0)                        
#>  htmltools     0.5.0      2020-06-16 [1] CRAN (R 4.0.0)                        
#>  ITKR          0.5.3.2.0  2020-06-01 [1] Github (stnava/ITKR@9bdd5f8)          
#>  knitr         1.29.3     2020-07-24 [1] Github (muschellij2/knitr@55195ad)    
#>  lattice       0.20-41    2020-04-02 [1] CRAN (R 4.0.0)                        
#>  magrittr      1.5        2014-11-22 [1] CRAN (R 4.0.0)                        
#>  Matrix        1.2-18     2019-11-27 [1] CRAN (R 4.0.0)                        
#>  memoise       1.1.0      2017-04-21 [1] CRAN (R 4.0.0)                        
#>  pkgbuild      1.0.8      2020-05-07 [1] CRAN (R 4.0.0)                        
#>  pkgload       1.1.0      2020-05-29 [1] CRAN (R 4.0.0)                        
#>  prettyunits   1.1.1      2020-01-24 [1] CRAN (R 4.0.0)                        
#>  processx      3.4.3      2020-07-05 [1] CRAN (R 4.0.0)                        
#>  ps            1.3.3      2020-05-08 [1] CRAN (R 4.0.0)                        
#>  R6            2.4.1      2019-11-12 [1] CRAN (R 4.0.0)                        
#>  Rcpp          1.0.5      2020-07-06 [1] CRAN (R 4.0.0)                        
#>  RcppEigen     0.3.3.7.0  2019-11-16 [1] CRAN (R 4.0.0)                        
#>  remotes       2.1.1      2020-02-15 [1] CRAN (R 4.0.0)                        
#>  rlang         0.4.6.9000 2020-07-08 [1] Github (r-lib/rlang@f3cdf54)          
#>  rmarkdown     2.3        2020-06-18 [1] CRAN (R 4.0.0)                        
#>  rprojroot     1.3-2      2018-01-03 [1] CRAN (R 4.0.0)                        
#>  sessioninfo   1.1.1      2018-11-05 [1] CRAN (R 4.0.0)                        
#>  stringi       1.4.6      2020-02-17 [1] CRAN (R 4.0.0)                        
#>  stringr       1.4.0      2019-02-10 [1] CRAN (R 4.0.0)                        
#>  testthat      2.3.2      2020-03-02 [1] CRAN (R 4.0.0)                        
#>  usethis       1.6.1.9000 2020-07-06 [1] Github (r-lib/usethis@c1cb6d9)        
#>  withr         2.2.0      2020-04-20 [1] CRAN (R 4.0.0)                        
#>  xfun          0.15       2020-06-21 [1] CRAN (R 4.0.0)                        
#>  yaml          2.2.1      2020-02-01 [1] CRAN (R 4.0.0)                        
#> 
#> [1] /Library/Frameworks/R.framework/Versions/4.0/Resources/library

Setting after Loading ANTsR - Works

library( ANTsR )
#> Loading required package: ANTsRCore
#> 
#> Attaching package: 'ANTsRCore'
#> The following objects are masked from 'package:stats':
#> 
#>     sd, var
#> The following objects are masked from 'package:base':
#> 
#>     all, any, apply, max, min, prod, range, sum

Sys.setenv(ANTS_RANDOM_SEED = 123)
Sys.setenv(ITK_GLOBAL_DEFAULT_NUMBER_OF_THREADS = 1)

fi <- antsImageRead(getANTsRData("r16") )
mi <- antsImageRead(getANTsRData("r64") )

tx = "SyN"
print(Sys.getenv("ANTS_RANDOM_SEED"))
#> [1] "123"
mytx1 <- antsRegistration(fixed=fi, mi, typeofTransform = tx)
mytx2 <- antsRegistration(fixed=fi, mi, typeofTransform = tx)
print(sum(mytx1$warpedmovout - mytx2$warpedmovout))
#> [1] 0
max(abs(as.array(mytx1$warpedmovout) - as.array(mytx2$warpedmovout)))
#> [1] 0

Created on 2020-08-06 by the reprex package (v0.3.0)

Session info
devtools::session_info()
#> ─ Session info ───────────────────────────────────────────────────────────────
#>  setting  value                       
#>  version  R version 4.0.0 (2020-04-24)
#>  os       macOS Mojave 10.14.6        
#>  system   x86_64, darwin17.0          
#>  ui       X11                         
#>  language (EN)                        
#>  collate  en_US.UTF-8                 
#>  ctype    en_US.UTF-8                 
#>  tz       America/New_York            
#>  date     2020-08-06                  
#> 
#> ─ Packages ───────────────────────────────────────────────────────────────────
#>  package     * version    date       lib source                                
#>  ANTsR       * 0.5.6.1    2020-06-01 [1] Github (ANTsX/ANTsR@9c7c9b7)          
#>  ANTsRCore   * 0.7.4.6    2020-07-07 [1] Github (muschellij2/ANTsRCore@61c37a1)
#>  assertthat    0.2.1      2019-03-21 [1] CRAN (R 4.0.0)                        
#>  backports     1.1.8      2020-06-17 [1] CRAN (R 4.0.0)                        
#>  callr         3.4.3      2020-03-28 [1] CRAN (R 4.0.0)                        
#>  cli           2.0.2      2020-02-28 [1] CRAN (R 4.0.0)                        
#>  crayon        1.3.4      2017-09-16 [1] CRAN (R 4.0.0)                        
#>  desc          1.2.0      2020-06-01 [1] Github (muschellij2/desc@b0c374f)     
#>  devtools      2.3.0.9000 2020-06-01 [1] Github (hadley/devtools@202ea81)      
#>  digest        0.6.25     2020-02-23 [1] CRAN (R 4.0.0)                        
#>  ellipsis      0.3.1      2020-05-15 [1] CRAN (R 4.0.0)                        
#>  evaluate      0.14       2019-05-28 [1] CRAN (R 4.0.0)                        
#>  fansi         0.4.1      2020-01-08 [1] CRAN (R 4.0.0)                        
#>  fs            1.4.2      2020-06-30 [1] CRAN (R 4.0.0)                        
#>  glue          1.4.1      2020-05-13 [1] CRAN (R 4.0.0)                        
#>  highr         0.8        2019-03-20 [1] CRAN (R 4.0.0)                        
#>  htmltools     0.5.0      2020-06-16 [1] CRAN (R 4.0.0)                        
#>  ITKR          0.5.3.2.0  2020-06-01 [1] Github (stnava/ITKR@9bdd5f8)          
#>  knitr         1.29.3     2020-07-24 [1] Github (muschellij2/knitr@55195ad)    
#>  lattice       0.20-41    2020-04-02 [1] CRAN (R 4.0.0)                        
#>  magrittr      1.5        2014-11-22 [1] CRAN (R 4.0.0)                        
#>  Matrix        1.2-18     2019-11-27 [1] CRAN (R 4.0.0)                        
#>  memoise       1.1.0      2017-04-21 [1] CRAN (R 4.0.0)                        
#>  pkgbuild      1.0.8      2020-05-07 [1] CRAN (R 4.0.0)                        
#>  pkgload       1.1.0      2020-05-29 [1] CRAN (R 4.0.0)                        
#>  prettyunits   1.1.1      2020-01-24 [1] CRAN (R 4.0.0)                        
#>  processx      3.4.3      2020-07-05 [1] CRAN (R 4.0.0)                        
#>  ps            1.3.3      2020-05-08 [1] CRAN (R 4.0.0)                        
#>  R6            2.4.1      2019-11-12 [1] CRAN (R 4.0.0)                        
#>  Rcpp          1.0.5      2020-07-06 [1] CRAN (R 4.0.0)                        
#>  RcppEigen     0.3.3.7.0  2019-11-16 [1] CRAN (R 4.0.0)                        
#>  remotes       2.1.1      2020-02-15 [1] CRAN (R 4.0.0)                        
#>  rlang         0.4.6.9000 2020-07-08 [1] Github (r-lib/rlang@f3cdf54)          
#>  rmarkdown     2.3        2020-06-18 [1] CRAN (R 4.0.0)                        
#>  rprojroot     1.3-2      2018-01-03 [1] CRAN (R 4.0.0)                        
#>  sessioninfo   1.1.1      2018-11-05 [1] CRAN (R 4.0.0)                        
#>  stringi       1.4.6      2020-02-17 [1] CRAN (R 4.0.0)                        
#>  stringr       1.4.0      2019-02-10 [1] CRAN (R 4.0.0)                        
#>  testthat      2.3.2      2020-03-02 [1] CRAN (R 4.0.0)                        
#>  usethis       1.6.1.9000 2020-07-06 [1] Github (r-lib/usethis@c1cb6d9)        
#>  withr         2.2.0      2020-04-20 [1] CRAN (R 4.0.0)                        
#>  xfun          0.15       2020-06-21 [1] CRAN (R 4.0.0)                        
#>  yaml          2.2.1      2020-02-01 [1] CRAN (R 4.0.0)                        
#> 
#> [1] /Library/Frameworks/R.framework/Versions/4.0/Resources/library

Setting Seed after reading data - FAILS

library( ANTsR )
#> Loading required package: ANTsRCore
#> 
#> Attaching package: 'ANTsRCore'
#> The following objects are masked from 'package:stats':
#> 
#>     sd, var
#> The following objects are masked from 'package:base':
#> 
#>     all, any, apply, max, min, prod, range, sum

fi <- antsImageRead(getANTsRData("r16") )
mi <- antsImageRead(getANTsRData("r64") )

Sys.setenv(ANTS_RANDOM_SEED = 123)
Sys.setenv(ITK_GLOBAL_DEFAULT_NUMBER_OF_THREADS = 1)
tx = "SyN"
print(Sys.getenv("ANTS_RANDOM_SEED"))
#> [1] "123"
mytx1 <- antsRegistration(fixed=fi, mi, typeofTransform = tx)
mytx2 <- antsRegistration(fixed=fi, mi, typeofTransform = tx)
print(sum(mytx1$warpedmovout - mytx2$warpedmovout))
#> [1] 117.2713
max(abs(as.array(mytx1$warpedmovout) - as.array(mytx2$warpedmovout)))
#> [1] 103.0909

Created on 2020-08-06 by the reprex package (v0.3.0)

Session info
devtools::session_info()
#> ─ Session info ───────────────────────────────────────────────────────────────
#>  setting  value                       
#>  version  R version 4.0.0 (2020-04-24)
#>  os       macOS Mojave 10.14.6        
#>  system   x86_64, darwin17.0          
#>  ui       X11                         
#>  language (EN)                        
#>  collate  en_US.UTF-8                 
#>  ctype    en_US.UTF-8                 
#>  tz       America/New_York            
#>  date     2020-08-06                  
#> 
#> ─ Packages ───────────────────────────────────────────────────────────────────
#>  package     * version    date       lib source                                
#>  ANTsR       * 0.5.6.1    2020-06-01 [1] Github (ANTsX/ANTsR@9c7c9b7)          
#>  ANTsRCore   * 0.7.4.6    2020-07-07 [1] Github (muschellij2/ANTsRCore@61c37a1)
#>  assertthat    0.2.1      2019-03-21 [1] CRAN (R 4.0.0)                        
#>  backports     1.1.8      2020-06-17 [1] CRAN (R 4.0.0)                        
#>  callr         3.4.3      2020-03-28 [1] CRAN (R 4.0.0)                        
#>  cli           2.0.2      2020-02-28 [1] CRAN (R 4.0.0)                        
#>  crayon        1.3.4      2017-09-16 [1] CRAN (R 4.0.0)                        
#>  desc          1.2.0      2020-06-01 [1] Github (muschellij2/desc@b0c374f)     
#>  devtools      2.3.0.9000 2020-06-01 [1] Github (hadley/devtools@202ea81)      
#>  digest        0.6.25     2020-02-23 [1] CRAN (R 4.0.0)                        
#>  ellipsis      0.3.1      2020-05-15 [1] CRAN (R 4.0.0)                        
#>  evaluate      0.14       2019-05-28 [1] CRAN (R 4.0.0)                        
#>  fansi         0.4.1      2020-01-08 [1] CRAN (R 4.0.0)                        
#>  fs            1.4.2      2020-06-30 [1] CRAN (R 4.0.0)                        
#>  glue          1.4.1      2020-05-13 [1] CRAN (R 4.0.0)                        
#>  highr         0.8        2019-03-20 [1] CRAN (R 4.0.0)                        
#>  htmltools     0.5.0      2020-06-16 [1] CRAN (R 4.0.0)                        
#>  ITKR          0.5.3.2.0  2020-06-01 [1] Github (stnava/ITKR@9bdd5f8)          
#>  knitr         1.29.3     2020-07-24 [1] Github (muschellij2/knitr@55195ad)    
#>  lattice       0.20-41    2020-04-02 [1] CRAN (R 4.0.0)                        
#>  magrittr      1.5        2014-11-22 [1] CRAN (R 4.0.0)                        
#>  Matrix        1.2-18     2019-11-27 [1] CRAN (R 4.0.0)                        
#>  memoise       1.1.0      2017-04-21 [1] CRAN (R 4.0.0)                        
#>  pkgbuild      1.0.8      2020-05-07 [1] CRAN (R 4.0.0)                        
#>  pkgload       1.1.0      2020-05-29 [1] CRAN (R 4.0.0)                        
#>  prettyunits   1.1.1      2020-01-24 [1] CRAN (R 4.0.0)                        
#>  processx      3.4.3      2020-07-05 [1] CRAN (R 4.0.0)                        
#>  ps            1.3.3      2020-05-08 [1] CRAN (R 4.0.0)                        
#>  R6            2.4.1      2019-11-12 [1] CRAN (R 4.0.0)                        
#>  Rcpp          1.0.5      2020-07-06 [1] CRAN (R 4.0.0)                        
#>  RcppEigen     0.3.3.7.0  2019-11-16 [1] CRAN (R 4.0.0)                        
#>  remotes       2.1.1      2020-02-15 [1] CRAN (R 4.0.0)                        
#>  rlang         0.4.6.9000 2020-07-08 [1] Github (r-lib/rlang@f3cdf54)          
#>  rmarkdown     2.3        2020-06-18 [1] CRAN (R 4.0.0)                        
#>  rprojroot     1.3-2      2018-01-03 [1] CRAN (R 4.0.0)                        
#>  sessioninfo   1.1.1      2018-11-05 [1] CRAN (R 4.0.0)                        
#>  stringi       1.4.6      2020-02-17 [1] CRAN (R 4.0.0)                        
#>  stringr       1.4.0      2019-02-10 [1] CRAN (R 4.0.0)                        
#>  testthat      2.3.2      2020-03-02 [1] CRAN (R 4.0.0)                        
#>  usethis       1.6.1.9000 2020-07-06 [1] Github (r-lib/usethis@c1cb6d9)        
#>  withr         2.2.0      2020-04-20 [1] CRAN (R 4.0.0)                        
#>  xfun          0.15       2020-06-21 [1] CRAN (R 4.0.0)                        
#>  yaml          2.2.1      2020-02-01 [1] CRAN (R 4.0.0)                        
#> 
#> [1] /Library/Frameworks/R.framework/Versions/4.0/Resources/library

Example video demonstrating this: https://www.youtube.com/watch?v=tyz5Pbp-35o&feature=youtu.be

@muschellij2
Copy link
Collaborator

Can anyone confirm this on their machine? @adigherman ? @cookpa @stnava

@cookpa
Copy link
Member

cookpa commented Sep 30, 2020

Will try but may take some days

@giadadam
Copy link

Hello! I would like to store the seed that was used in a specific run, so that we are able to reproduce that data in the future. How can this be done? Thanks in advance!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

8 participants