ENH: compute grey and white matter masks

[alpinho]

This PR tries to address the issue #2487.

[alpinho]

I got 2 errors (1 error duplicated) in my local machine with pytest in lines 632 and 724 of nilearn/masking.py, which might be causing many of these failing tests. The problem narrows down to line 223 nilearn/datasets/struct.py with the following log:

*** nibabel.filebasedimages.ImageFileError: Empty file: '/tmp/pytest-of-analu/pytest-35/temp_nilearn_home4/nilearn_shared_data/icbm152_2009/mni_icbm152_nlin_sym_09a/mni_icbm152_gm_tal_nlin_sym_09a.nii.gz'

I don't understand because at the end of the function check_niimg in file nilearn/_utils/niimg_conversions.py, the niimg seems to have been successfully created (it is not empty!), but then there is a problem with the assignment to gm_img in line 223 nilearn/datasets/struct.py. I have the impression that the issue might come from line 134 of nilearn/_utils/niimg.py where dtype is None. Any advice? Thx in advance.

[NicolasGensollen]

Thanks for working on this @alpinho ! 👍
I left a few comments below.

[NicolasGensollen]

@alpinho I looked at the errors you're getting with the empty niimg, and I think they come from the fact that you are using datasets fetchers in compute_brain_mask(), which are mocked when you run the tests.

The previous behaviour of compute_brain_mask() was to systematically rely on load_mni152_brain_mask() which doesn't download anything, but uses the template nilearn/nilearn/datasets/data/avg152T1_brain.nii.gz that is included in nilearn.

[alpinho]

Thanks @NicolasGensollen for the feedback! I am working on your requests. Regarding your last comment about my error, what should I do? Should I just discard this error, since this only happens when running the local tests? Note that the default behavior of compute_brain_mask() is still meant to work with the whole-brain mask. I am just adding two more options: the possibility to also compute the gray-matter or 'white-matter' masks.

[alpinho]

@NicolasGensollen I did several updates in order to both fix my previous local pytest errors and address your comments. To fix the local pytest errors, I decided to use the load functions instead of the fetch functions. Therefore, I had to download and store in nilearn/nilearn/datasets/data/ the gm and wm masks, similarly to what have been done for the whole-brain mask. I didn't abide for the moment with the DRY principles, as suggested by you, since I want to check first whether you agree with this approach. Besides, I am still getting some failing tests. Let me thus know what you think and how I shall proceed.

[NicolasGensollen]

Thanks for addressing my comments @alpinho ! 👍
I don't think adding the templates for gm and wm is a problem (@bthirion could confirm this or not). However, I'd be in favor of having something consistent in terms of API and data. I made a few comments below regarding this.

Maybe we should update the brain template also to have the same resolutions:

from nilearn import plotting
from nilearn.datasets.struct import (MNI152_FILE_PATH, 
                                     GM_MNI152_FILE_PATH, 
                                     WM_MNI152_FILE_PATH)
plotting.plot_img(MNI152_FILE_PATH, colorbar=True)
plotting.plot_img(GM_MNI152_FILE_PATH, colorbar=True)
plotting.plot_img(WM_MNI152_FILE_PATH, colorbar=True)

WDYT?

[alpinho]

WDYT?

I also prefer to wait for @bthirion comments. I am actually surprised that they don't share already the same resolution... I downloaded these masks from https://osf.io/7pj92/download, i.e. one of the links used by fetch_icbm152_2009 and, consequently, by fetch_icbm152_brain_gm_mask.

[alpinho]

@NicolasGensollen Could you advice what else I should fix?, since I still get many failing tests. I don't get more errors in my local machine with pytest. Thanks in advance!

[NicolasGensollen]

@alpinho I left a couple comments to solve the CircleCI failure and to clarify a point I don't understand.
I'll take a look at the test errors you're getting to see what might be causing them.

[NicolasGensollen]

@alpinho I had a look at the test failure that you get.
First of all, there is only one test failing: input_data/tests/test_multi_nifti_masker.py::test_compute_multi_gray_matter_mask.
In this test, a MultiNiftiMasker is fitted on some random data with a mask_strategy set to template, which triggers a call to compute_multi_gray_matter_mask in fit():

nilearn/nilearn/input_data/multi_nifti_masker.py

Lines 209 to 211 in d066c23

	elif self.mask_strategy == 'template':
	compute_mask = masking.compute_multi_gray_matter_mask
	else:

compute_multi_gray_matter_mask then calls compute_brain_mask with mask_type set to gm.
Within this function, the template is loaded and resampled to the provided images before being thresholded.
I initially thought that the thresholding was zeroing the template out, but it turns out that it is the resampling step which produces an empty template image (all zeros) for gm and wm.

Here is a little snippet that reproduces the issue:

import numpy as np
from nilearn.image import get_data
from numpy.testing import assert_array_equal
from nilearn._utils import check_niimg
from nibabel import Nifti1Image
from nilearn.masking import compute_brain_mask

rng = np.random.RandomState(42) 
imgs = [Nifti1Image(rng.uniform(size=(9, 9, 5)), np.eye(4)),
        Nifti1Image(rng.uniform(size=(9, 9, 5)), np.eye(4))] 
item_gen = check_niimg(imgs, return_iterator=True)
for _ in item_gen:
    pass
for mask_type in ['whole-brain', 'gm', 'wm']:
    print('-'*4 + ' ' + mask_type + ' ' + '-'*4)
    mask = compute_brain_mask(imgs[0], 
                              threshold=0.5, # This doen't change the results for gm and wm
                              connected=True, 
                              opening=2, 
                              mask_type=mask_type)
    assert mask.shape == (9, 9, 5)
    assert_array_equal(mask.affine, imgs[0].affine)
    print(np.sum(get_data(mask)))

---- whole-brain ----
25
---- gm ----
0
---- wm ----
0

Does this make sense?
I'll try to understand whether the zeroing of the template by the resampling step is normal or not. 🤔
If it is, we will most likely need to update the test.

[NicolasGensollen]

@alpinho I think the problem happens here (all selected slices are empty...):

nilearn/nilearn/image/resampling.py

Lines 579 to 580 in d066c23

	subset_indices = tuple(slice(0, s.stop-s.start) for s in slices)
	resampled_data[slices] = _get_data(cropped_img)[subset_indices]

We get to this line because we satisfy this condition:

nilearn/nilearn/image/resampling.py

Lines 552 to 555 in d066c23

	# if (A == I OR some combination of permutation(I) and sign-flipped(I)) AND
	# all(b == integers):
	if (np.all(np.eye(3) == A) and all(bt == np.round(bt) for bt in b) and
	not force_resample):

I haven't thought too much about it yet, but I think a quick fix would be to force the resampling in compute_brain_mask:

resampled_template = cache(resampling.resample_to_img, memory)(template, target_img, force_resample=True)

I tried it briefly and it seems to work as expected.

WDYT?

[codecov]

Codecov Report

Merging #2738 (373f03a) into main (7878897) will increase coverage by 0.07%.
The diff coverage is 99.05%.

@@            Coverage Diff             @@
##             main    #2738      +/-   ##
==========================================
+ Coverage   88.47%   88.54%   +0.07%     
==========================================
  Files         100      100              
  Lines       13588    13667      +79     
  Branches     2652     2670      +18     
==========================================
+ Hits        12022    12102      +80     
+ Misses        975      974       -1     
  Partials      591      591              

Impacted Files	Coverage Δ
nilearn/image/image.py	96.20% <ø> (ø)
nilearn/image/resampling.py	92.92% <ø> (ø)
nilearn/datasets/struct.py	91.09% <98.83%> (+4.65%)	⬆️
nilearn/datasets/__init__.py	100.00% <100.00%> (ø)
nilearn/input_data/multi_nifti_masker.py	89.87% <100.00%> (ø)
nilearn/masking.py	88.88% <100.00%> (+1.16%)	⬆️

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 7878897...373f03a. Read the comment docs.

[alpinho]

Thank you very much for your thorough and clear explanation @NicolasGensollen. Your advice also makes sense to me and I have implemented it in my last commit. All checks have passed now. Yet, "Codecov Report" still reports some red segments. Is this problematic?

[NicolasGensollen]

Thanks for adressing my reviews @alpinho !
I think this is getting close now! 🎉
I left a few comments below. In addition, we also need to decide whether we should update the brain template to have the same resolutions (see one of my previous comments).
Maybe @bthirion has an opinion on that?

[alpinho]

@jeromedockes But that solution does not work on my side. This is my error log:

---------------------------------------------------------------------------
UnboundLocalError                         Traceback (most recent call last)
<ipython-input-3-aa844c641ab4> in <module>
----> 1 load_mni152_gm_template()

~/mygit/nilearn/nilearn/datasets/struct.py in load_mni152_gm_template(resolution)
    210 
    211     if resolution is None:
--> 212         if not _MNI_RES_WARNING_ALREADY_SHOWN:
    213             warnings.warn("Default resolution of the MNI template will change "
    214                           "from 2mm to 1mm in version 0.10.0", FutureWarning)

UnboundLocalError: local variable '_MNI_RES_WARNING_ALREADY_SHOWN' referenced before assignment

That's why I decided to do a class. Is there any other workaround?

[jeromedockes]

UnboundLocalError: local variable '_MNI_RES_WARNING_ALREADY_SHOWN' referenced before assignment ``` That's why I decided to do a class. Is there any other workaround?

right, sorry I forgot you need to add `global _MNI_RES_WARNING_ALREADY_SHOWN` in the function before using this variable

[alpinho]

UnboundLocalError: local variable '_MNI_RES_WARNING_ALREADY_SHOWN' referenced before assignment ``` That's why I decided to do a class. Is there any other workaround?
right, sorry I forgot you need to add global _MNI_RES_WARNING_ALREADY_SHOWN in the function before using this variable

Ah! Great! Yes, now it seems it is working as expected. Well, I guess I have addressed all of your comments. Let us know whether it looks good from your side. Thanks!

[jeromedockes]

2 small details

[jeromedockes]

Well, I guess I have addressed all of your comments. Let us know whether it looks good from your side. Thanks!

Yes! thanks again for all this work. when having a final look I saw 2 small issues that will be easily fixed (see other comments). otherwise it LGTM! Let's give @bthirion and @NicolasGensollen time to do another review if they want then I'll merge it

[alpinho]

@jeromedockes One last question: codecov is complaining because it seems that the lines of the warning messages are not being tested now. I am not sure how to fix this now with the global var.

[jeromedockes]

One last question: codecov is complaining because it seems that the lines of the warning messages are not being tested now

you can add a test like this one in test_struct.py:

@pytest.mark.parametrize("part", ["_brain", "_gm", "_wm"])
@pytest.mark.parametrize("kind", ["template", "mask"])
def test_mni152_resolution_warnings(part, kind):
    struct._MNI_RES_WARNING_ALREADY_SHOWN = False
    if kind == "template" and part == "_brain":
        part = ""
    loader = getattr(struct, f"load_mni152{part}_{kind}")
    try:
        loader.cache_clear()
    except AttributeError:
        pass
    with warnings.catch_warnings(record=True) as w:
        loader(resolution=1)
    assert len(w) == 0
    with warnings.catch_warnings(record=True) as w:
        loader()
        loader()
    assert len(w) == 1

[alpinho]

I think functools is throwing some errors for Python <= 3.7

[jeromedockes]

Ok. But this last snippet of code is already written in the function, exactly after defining the `loader`, no?

yes, I edited it after posting this comment

[jeromedockes]

I think functools is throwing some errors for Python <= 3.7

replace ``` functools.lru_cache ``` with ``` @functools.lru_cache(maxsize=3) ```

[alpinho]

Hello! @jeromedockes , is there anything else that should be fixed? @NicolasGensollen and @bthirion , let us know what you think. Thx!

[jeromedockes]

Hello! @jeromedockes , is there anything else that should be fixed? @NicolasGensollen and @bthirion , let us know what you think. Thx!

everything LGTM! I'm just waiting to see if there are more comments from the other reviewers before merging :)

[NicolasGensollen]

LGTM too! 👍
Thanks a lot!

[jeromedockes]

thanks a lot @alpinho !