james-cole / brainageR

Software for generating a brain-predicted age value, using Gaussian Processes regression, implemented in R
GNU Lesser General Public License v3.0
71 stars 25 forks source link

SPM12 version #3

Closed samanemati89 closed 2 years ago

samanemati89 commented 2 years ago

Hi There,

May I ask what was the version of SPM12 you used when you developed this brainageR toolbox? Does it matter which version of SPM12 we are using to run this toolbox?

james-cole commented 2 years ago

brainageR v2.1 used SPM12 release 7219 (https://www.fil.ion.ucl.ac.uk/spm/software/spm12/), though my understanding is that it should be compatible with the subsequent releases of SPM12. The key SPM12 tools used are Segment and DARTEL, so unless these have changed, it should be fine.

samanemati89 commented 2 years ago

Thank you so much for your response. I ran brainage on T1 scans of 60 healthy people. When I run this using the most updated version of SPM12 (7771), results are way different compared to when using SPM12 version 7487. For example, for a person with chronological age of 24, predicted brain age is 30 when using SPM12 version 7487. For the same person, predicted brain age is 51 when using the version 7771. The results from the previous version of SPM makes more sense.

james-cole commented 2 years ago

Thanks for sharing this feedback. Useful to know. I recommend that you check the parameters of the SPM12 Segment and DARTEL and how the defaults differ between versions. It's possible that some of these settings have altered, causing the inconsistency. Also, did you visually QC the segmentation results? If the segmentation has failed for one version, this would also cause a discrepancy.

neurolabusc commented 2 years ago

@james-cole thanks for the quick response. In our own dataset, we noticed that SPM12 r7487 generated much more accurate brain age estimates than r7771. The latter tends to estimate the individual is older than their chronological age. For example, we had a 24 year old where r7487 estimated 30.7 while r7771 estimated 51.8. Our data was acquired on a Siemens Prisma and was defaced. Therefore, we decided to check if this was specific to our scans. We evaluated data from a Philips Ingenia Elition X and observed the same influence.

We isolated the problem to the SPM12 spatial.preproc function. Here is the sample Philips data set and minimal script. The script is based on brainageR, and changes a few of the SPM12 segment defaults for this stage. In particular the number of Gaussians used for each tissue type.

Stepping through the patches we observe that r7505 changes the results a tiny bit, but the big change is r7593:

r7593 | john | 2019-05-20 19:58:16 +0100 (Mon, 20 May 2019) | 28 lines
Changed paths:
   M /trunk/spm_diffeo.m
   M /trunk/spm_preproc8.m
   M /trunk/src/shoot_boundary.h
   M /trunk/src/shoot_dartel.c
   M /trunk/src/shoot_diffeo3d.c
   M /trunk/src/shoot_diffeo3d.h
   M /trunk/src/shoot_optim3d.c
   M /trunk/src/shoot_regularisers.c
   M /trunk/src/spm_diffeo.c
   M /trunk/toolbox/Shoot/spm_shoot_greens.m

Extensive changes, ready to break things during the SPM course ;o)
* Fixed a bug that affected the full-multigrid optimisation used by many of the
  nonlinear registration algorithms.  The bug most strongly effected registration of
  2D images, but is likely have made 3D registration a tiny bit worse (because of its
  effect of the lower resolution stages of the multigrid solver). This bug was in
  shoot_regularisers.c, so recompilation of mex files is needed.
* Increased stability throughout shoot_regularisers.c by multiplying values on the
  diagonal of the differential operators by 1.000001.  This ensures things are more
  likely to remain positive definite, and makes spm_shoot_greens more stable.
* Re-done spm_shoot_greens so that things are a little bit more stable.  Greens
  functions of differential operators are now computed using double precision. Only
  the centre (DC) element of the FFT is handled differently in situations where the
  differential operator may be singular. The changes to shoot_regularisers.c
  (multiplying the diagonal by 1.000001) should mean that the DC element of the fft
  of the differential operator is now the only one that might approach zero.
* Changed the way that compositions of deformations (by spm_diffeo) are done so that
  unwrapping should no longer be needed.
* The deformation unwrapping procedure (and sampn_vox) has been moved to shoot_dartel.c
  as this is the only place it is used.  Also fixed the bug that caused seg faults
  when composing warps of different sizes (passed the wrog dimensions to unwrap).
* Removed the mexFunction for the Lie bracket stuff as it is not used (although the
  code is retained for possible future reference.
* Added a function for computing image gradients, which handles boundary conditions
  better. Gradients computed by bsplins were maybe not as suitable because the mirror
  boundary conditions it uses are not quite the same as those used by push and pull.
* Included a slight change to how the wp parameters are regularised in spm_preproc8.m.

r7487 r7771

While the more recent release may be internally consistent, it does not seem to match the training dataset used for brainageR.

@gllmflndn I looked on jiscmail and could not find a clear example of a report on this. Is this a well known issue? Is it advised that brainageR be re-trained on output from the latest version?

neurolabusc commented 2 years ago

r7593 seems to have an adverse outcome when applied to cropped images. Here is a link to two copies of the same Philips T1 scan, the only difference is that one is cropped to remove excess signal in the head-foot direction. Here I run SPM12 R7771 segment command with the default parameters for SPM (e.g. not the brainageR modifications). Cropping the image from 170x256x256 to 170x256x177 in the LR, AP, IS dimensions dramatically impacts the quality of the images. Our team often crops images: some FSL tools use the center of brightness as a starting estimate. Since obesity is a risk factor in our population, many individuals have excess neck fat. The issue is introduced by r7593 and is consistent across operating systems (macOS and Linux) and Matlab versions (2020a, 2021b).

@JohnAshburner do you have any thoughts on this?

Images show SPM's estimate for bone tissue as red superimposed on the T1 anatomical scan. The two images are identical, except one is cropped. Note the cropped image shows dramatic errors.

crop nocrop

james-cole commented 2 years ago

Thanks all for the interesting discussion and update. Unfortunately, re-processing and re-training the brainageR model based on the latest version isn't feasible at the moment, though I will amend the guidance to recommed the specific SPM12 version.

neurolabusc commented 2 years ago

@james-cole this is clearly an issue that is general to SPM12 and not specific to brainageR. My informal test of the restricted development branch suggests this issue has already been resolved by the SPM team. Give our team a few days to check this: we can ensure that the age estimates from the latest patches (8220) are as good as our previous findings with r7487. Perhaps @JohnAshburner or @gllmflndn can provide some insight.

Assuming I am correct that the latest change resolves the issue, I would suggest the brainageR script generate an error if the version of SPM12 is not supported:

[v,r] = spm('Ver','',1);
r = str2double(r);
if (r > 7592) && (r < 8220)
    error('brainageR may be unreliable with this release of SPM12 (issue 3). Either upgrade or downgrade SPM12')
end 

Result for r8220: r8220

james-cole commented 2 years ago

Thanks @neurolabusc, much appreciated. I'll wait to hear that outcome of your further tests, but can certainly add that line as you suggest.

JohnAshburner commented 2 years ago

There does indeed seem to be an issue with 7771, which is something that is fixed in the developmental version (spm_preproc8.m 7965 2020-09-28). The fix is simple, and involves changing line 800, from:

            prm(6)       = param(6)*scal^2;

to:

            prm(6)       = param(6)*scal;

My svn log entry was:

* Multigrid unstable if param(6) too large relative to the Alpha (likelihood Hessian)