Compute_contrast as second level produces inconsistent error

Software Support
1

Hi there! This is a clone of the github issue compute_contrast as second level produces inconsistent error · Issue #3562 · nilearn/nilearn · GitHub which, as said by bthirion in the first comment, should be here rather than in the nilearn issues.

Running compute_contrast for second level analysis crashes because it detects mismatching shapes of ref and image, but the actual shapes are just fine.

Nilearn version:

0.8.1

Expected behavior

compute correctly the contrasts

Actual behavior

produce a ValueError

Steps and code to reproduce bug

Unfortunately I was not able to reproduce this bug using nilearn data.
In the minimalistic version of my bug, I run it with only two sets of participants: for set 1, everything works just fine, while for set 2, I get the error. Note that the first level results have been derived using the same script and produced no error. Of course, all my participants underwent the same exact experimental paradigm, and fMRIPrep produced no errors when preprocessing these participants.

import numpy as np
import pandas as pd
from nilearn import plotting
from nilearn.glm.second_level import SecondLevelModel
from nilearn.plotting import plot_design_matrix

cmaps_cond1 = ['/path/to/sub1-con1', '/path/to/sub2-con1']
cmaps_cond2 = ['/path/to/sub1-con2', '/path/to/sub2-con2']

n_subjects = 2
second_level_input = cmaps_cond1 + cmaps_cond2
condition_effect = np.hstack(([1] * n_subjects, [- 1] * n_subjects))
subject_effect = np.vstack((np.eye(n_subjects), np.eye(n_subjects)))
subjects = [f'S{i:02d}' for i in range(1, n_subjects + 1)]
paired_design_matrix = pd.DataFrame(
    np.hstack((condition_effect[:, np.newaxis], subject_effect)),
    columns=['cond1 vs cond2'] + subjects)
plot_design_matrix(paired_design_matrix)
second_level_model_paired = SecondLevelModel().fit(
    second_level_input, design_matrix=paired_design_matrix)

This, for some sub1, works fine, but then if I use it for another subject (call it sub3, hence replace sub1 by sub3 everywhere), it produces the following error (full error message below!):

stat_maps_paired = second_level_model_paired.compute_contrast(  'cond1 vs cond2', output_type='all')

gives

ValueError: Field of view of image #2 is different from reference FOV.
Reference affine:
array([[   2.25,    0.  ,    0.  ,  -90.  ],
       [   0.  ,    2.25,    0.  , -126.  ],
       [   0.  ,    0.  ,    3.  ,  -72.  ],
       [   0.  ,    0.  ,    0.  ,    1.  ]])
Image affine:
array([[   2.25      ,    0.        ,    0.        ,  -90.        ],
       [   0.        ,    2.25      ,    0.        , -126.        ],
       [   0.        ,    0.        ,    3.00099993,  -72.        ],
       [   0.        ,    0.        ,    0.        ,    1.        ]])
Reference shape:
(81, 97, 61)
Image shape:
(81, 97, 61, 1)

What strikes me is that if I run

from nilearn import image
for fn in second_level_input:
    print(image.get_data(fn).shape)

I get always (81, 97, 61) (namely, for set1 or set2), which is indeed the shape of all my cmaps (themselves derived from fMRIPrep preproc bold files in MNI152NLin6Asym space). So the cmaps of sub1, sub2 and sub3 have the same exact shapes…

I’ve tried to dig a bit in nilearn code to see where this error came from, bug this is a bit too deep for me. Here is the complete error:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Input In [66], in <cell line: 1>()
----> 1 stat_maps_paired = second_level_model_paired.compute_contrast(
      2                                                 'cond1 vs cond2',
      3                                                 output_type='all')

File ~/.local/lib/python3.8/site-packages/nilearn/glm/second_level/second_level.py:494, in SecondLevelModel.compute_contrast(self, second_level_contrast, first_level_contrast, second_level_stat_type, output_type)
    491 _check_effect_maps(effect_maps, self.design_matrix_)
    493 # Fit an Ordinary Least Squares regression for parametric statistics
--> 494 Y = self.masker_.transform(effect_maps)
    495 if self.memory:
    496     mem_glm = self.memory.cache(run_glm, ignore=['n_jobs'])

File ~/.local/lib/python3.8/site-packages/nilearn/input_data/base_masker.py:200, in BaseMasker.transform(self, imgs, confounds, sample_mask)
    197 self._check_fitted()
    199 if confounds is None and not self.high_variance_confounds:
--> 200     return self.transform_single_imgs(imgs,
    201                                       confounds=confounds,
    202                                       sample_mask=sample_mask)
    204 # Compute high variance confounds if requested
    205 all_confounds = []

File ~/.local/lib/python3.8/site-packages/nilearn/input_data/nifti_masker.py:518, in NiftiMasker.transform_single_imgs(self, imgs, confounds, sample_mask, copy)
    515     else:
    516         sample_mask = self._sample_mask
--> 518 data = self._cache(filter_and_mask,
    519                    ignore=['verbose', 'memory', 'memory_level',
    520                            'copy'],
    521                    shelve=self._shelving)(
    522     imgs, self.mask_img_, params,
    523     memory_level=self.memory_level,
    524     memory=self.memory,
    525     verbose=self.verbose,
    526     confounds=confounds,
    527     sample_mask=sample_mask,
    528     copy=copy,
    529     dtype=self.dtype
    530 )
    532 return data

File ~/.local/lib/python3.8/site-packages/joblib/memory.py:349, in NotMemorizedFunc.__call__(self, *args, **kwargs)
    348 def __call__(self, *args, **kwargs):
--> 349     return self.func(*args, **kwargs)

File ~/.local/lib/python3.8/site-packages/nilearn/input_data/nifti_masker.py:85, in filter_and_mask(imgs, mask_img_, parameters, memory_level, memory, verbose, confounds, sample_mask, copy, dtype)
     62 def filter_and_mask(imgs, mask_img_, parameters,
     63                     memory_level=0, memory=Memory(location=None),
     64                     verbose=0,
   (...)
     67                     copy=True,
     68                     dtype=None):
     69     """Extract representative time series using given mask.
     70 
     71     Parameters
   (...)
     83 
     84     """
---> 85     imgs = _utils.check_niimg(imgs, atleast_4d=True, ensure_ndim=4)
     87     # Check whether resampling is truly necessary. If so, crop mask
     88     # as small as possible in order to speed up the process
     90     if not _check_same_fov(imgs, mask_img_):

File ~/.local/lib/python3.8/site-packages/nilearn/_utils/niimg_conversions.py:280, in check_niimg(niimg, ensure_ndim, atleast_4d, dtype, return_iterator, wildcards)
    277     if return_iterator:
    278         return _iter_check_niimg(niimg, ensure_ndim=ensure_ndim,
    279                                  dtype=dtype)
--> 280     return concat_niimgs(niimg, ensure_ndim=ensure_ndim, dtype=dtype)
    282 # Otherwise, it should be a filename or a SpatialImage, we load it
    283 niimg = load_niimg(niimg, dtype=dtype)

File ~/.local/lib/python3.8/site-packages/nilearn/_utils/niimg_conversions.py:484, in concat_niimgs(niimgs, dtype, ensure_ndim, memory, memory_level, auto_resample, verbose)
    481 data = np.ndarray(target_shape + (sum(lengths), ),
    482                   order="F", dtype=dtype)
    483 cur_4d_index = 0
--> 484 for index, (size, niimg) in enumerate(zip(lengths, _iter_check_niimg(
    485         iterator, atleast_4d=True, target_fov=target_fov,
    486         memory=memory, memory_level=memory_level))):
    488     if verbose > 0:
    489         if isinstance(niimg, str):

File ~/.local/lib/python3.8/site-packages/nilearn/_utils/niimg_conversions.py:159, in _iter_check_niimg(niimgs, ensure_ndim, atleast_4d, target_fov, dtype, memory, memory_level, verbose)
    153             niimg = cache(image.resample_img, memory,
    154                           func_memory_level=2,
    155                           memory_level=memory_level)(
    156                     niimg, target_affine=ref_fov[0],
    157                     target_shape=ref_fov[1])
    158         else:
--> 159             raise ValueError(
    160                 "Field of view of image #%d is different from "
    161                 "reference FOV.\n"
    162                 "Reference affine:\n%r\nImage affine:\n%r\n"
    163                 "Reference shape:\n%r\nImage shape:\n%r\n"
    164                 % (i, ref_fov[0], niimg.affine, ref_fov[1],
    165                    niimg.shape))
    166     yield niimg
    167 except DimensionError as exc:
    168     # Keep track of the additional dimension in the error

ValueError: Field of view of image #2 is different from reference FOV.
Reference affine:
array([[   2.25,    0.  ,    0.  ,  -90.  ],
       [   0.  ,    2.25,    0.  , -126.  ],
       [   0.  ,    0.  ,    3.  ,  -72.  ],
       [   0.  ,    0.  ,    0.  ,    1.  ]])
Image affine:
array([[   2.25      ,    0.        ,    0.        ,  -90.        ],
       [   0.        ,    2.25      ,    0.        , -126.        ],
       [   0.        ,    0.        ,    3.00099993,  -72.        ],
       [   0.        ,    0.        ,    0.        ,    1.        ]])
Reference shape:
(81, 97, 61)
Image shape:
(81, 97, 61, 1)

Any guess?

2

For the sake of clarity, I paste here the relevant part of bthirion here:

[…] the issue that seems to occur is a small glitch in the affine of sub3 (the [2,2] coefficient is 3.00099993 instead of 3), which is indeed extremely surprising given what you describe.
Any chance you can share an image from sub3 and another one from say, sub1 ?

3

Thanks bthirion for your reply!

I checked and this weird affine value of 3.00099993 is already there at the level of the preproc images by fMRIPrep:

preproc_fn = '/path/to/fmriprep/sub-3/func/sub-3_task-sometaskname_space-MNI152NLin6Asym_desc-preproc_bold.nii.gz'
nilearn.image.loa_img(preproc_fn).affine

gives

array([[   2.25      ,    0.        ,    0.        ,  -90.        ],
       [   0.        ,    2.25      ,    0.        , -126.        ],
       [   0.        ,    0.        ,    3.00099993,  -72.        ],
       [   0.        ,    0.        ,    0.        ,    1.        ]])

Of course, for the ‘normal’ subjects (say sub1), I get

array([[   2.25      ,    0.        ,    0.        ,  -90.        ],
       [   0.        ,    2.25      ,    0.        , -126.        ],
       [   0.        ,    0.        ,    3,  -72.        ],
       [   0.        ,    0.        ,    0.        ,    1.        ]])

So maybe this isn’t an issue with nilearn, maybe some fmriprep experts will have to say something about it?

Fmriprep version 21.0.4 was used.

PS1: by the way, running much more basic commands like

image.load_img(cmaps_cond2)

directly catches the inconsistency between the affines (essentially I get the same error).

PS2: a quick and dirty fix is as follows (basically loading the images and rounding the affines by hand)

import numpy as np
import pandas as pd
from nilearn import plotting
from nilearn.glm.second_level import SecondLevelModel
from nilearn.plotting import plot_design_matrix

cmaps_cond1 = ['/path/to/sub1-con1', '/path/to/sub3-con1']
cmaps_cond2 = ['/path/to/sub1-con2', '/path/to/sub3-con2']

n_subjects = 2

cmaps_c1_imgs_sub1 = image.load_img(cmaps_cond1[0])
cmaps_c1_imgs_sub3 = image.load_img(cmaps_cond1[1])
cmaps_c2_imgss_sub1 = image.load_img(cmaps_cond2[0])
cmaps_c2_imgs_sub3 = image.load_img(cmaps_cond2[1]) 

# say the only image with affine having this weird 3.00099993 is cmaps_c2_imgs_sub3

cmaps_c2_imgs_sub3._affine = np.around(cmaps_c2_imgs_sub3.affine, 2)

second_level_input = [cmaps_c1_imgs_sub1] + [cmaps_c1_imgs_sub3] + [cmaps_c2_imgs_sub1] + [cmaps_c2_imgs_sub3]
condition_effect = np.hstack(([1] * n_subjects, [- 1] * n_subjects))
subject_effect = np.vstack((np.eye(n_subjects), np.eye(n_subjects)))
subjects = [f'S{i:02d}' for i in range(1, n_subjects + 1)]
paired_design_matrix = pd.DataFrame(
    np.hstack((condition_effect[:, np.newaxis], subject_effect)),
    columns=['cond1 vs cond2'] + subjects)
plot_design_matrix(paired_design_matrix)
second_level_model_paired = SecondLevelModel().fit(
    second_level_input, design_matrix=paired_design_matrix)

stat_maps_paired = second_level_model_paired.compute_contrast(
                                                'cond1 vs cond2',
                                                output_type='all')
4

(This starts to looks like a monologue but whatever.)

There was other problematic subjects in my datasets. For these, the [2,2] entry of affine was not 3 nor 3.00099993, but rather 2.999-something. While this should be “fixed” by the quick and dirty trick mentioned above, it had another consequence: shape[2] was now 62 instead of 61.

To solve this in a (hopefully) clean way, I decided to resample all images to some ref image (taken to be the first one) before feeding it into the second level tools. In short, something like:

ref = cmaps_cond1[0]
cmaps_cond1_resampled = [image.resample_to_img(img, ref) for img in cmaps_cond1]
cmaps_cond2_resampled = [image.resample_to_img(img, ref) for img in cmaps_cond2]

(note that I purposely use cmaps_cond1[0] as the ref for the images from cond2).

5

Thx for posting.
As an alternative, you can use nilearn.image.resample_to_img solution, but this involves resampling your data, and the affine difference is so tiny that I think it is ride to simply override it.
Best,
Bertrand