diff --git a/bdpy/mri/fmriprep.py b/bdpy/mri/fmriprep.py index 6567a69..7e786d0 100644 --- a/bdpy/mri/fmriprep.py +++ b/bdpy/mri/fmriprep.py @@ -11,7 +11,6 @@ from collections import OrderedDict import nibabel -import nipy import numpy as np import pandas as pd @@ -406,17 +405,17 @@ def __load_volume(self): - Returns voxel ijk indexes (3 x voxel) - Data, xyz, and ijk are flattened by Fortran-like index order """ - img = nipy.load_image(self.__dpath) + img = nibabel.load(self.__dpath) data = img.get_fdata() if data.ndim == 4: data = data.reshape(-1, data.shape[-1], order='F').T i_len, j_len, k_len, t = img.shape - affine = np.delete(np.delete(img.coordmap.affine, 3, axis=0), 3, axis=1) + affine = img.affine elif data.ndim == 3: data = data.flatten(order='F') i_len, j_len, k_len = img.shape - affine = img.coordmap.affine + affine = img.affine else: raise ValueError('Invalid shape.') @@ -838,17 +837,17 @@ def __load_mri(fpath): - Returns voxel ijk indexes (3 x voxel) - Data, xyz, and ijk are flattened by Fortran-like index order """ - img = nipy.load_image(fpath) + img = nibabel.load(fpath) data = img.get_fdata() if data.ndim == 4: data = data.reshape(-1, data.shape[-1], order='F').T i_len, j_len, k_len, t = img.shape - affine = np.delete(np.delete(img.coordmap.affine, 3, axis=0), 3, axis=1) + affine = np.delete(np.delete(img.affine, 3, axis=0), 3, axis=1) elif data.ndim == 3: data = data.flatten(order='F') i_len, j_len, k_len = img.shape - affine = img.coordmap.affine + affine = img.affine else: raise ValueError('Invalid shape.') diff --git a/bdpy/mri/glm.py b/bdpy/mri/glm.py index 2916cac..c775c96 100644 --- a/bdpy/mri/glm.py +++ b/bdpy/mri/glm.py @@ -2,6 +2,7 @@ import csv +import warnings import numpy as np from nipy.modalities.fmri.experimental_paradigm import ( @@ -45,6 +46,14 @@ def make_paradigm(event_files, num_vols, tr=2., cond_col=2, label_col=None, regr run_regressors : nuisance regressors for runs run_regressors_label : labels for the run regressors """ + warnings.warn( + "make_paradigm depends on nipy, which is unmaintained. nipy support " + "will be removed in a future release of bdpy. Consider migrating to " + "nilearn (e.g. nilearn.glm.first_level.make_first_level_design_matrix).", + DeprecationWarning, + stacklevel=2, + ) + onset = [] duration = [] conds = [] diff --git a/bdpy/mri/load_epi.py b/bdpy/mri/load_epi.py index 7e1d87b..95ec1d4 100644 --- a/bdpy/mri/load_epi.py +++ b/bdpy/mri/load_epi.py @@ -6,7 +6,7 @@ import itertools as itr -import nipy +import nibabel import numpy as np @@ -35,7 +35,7 @@ def load_epi(datafiles): print("Loading %s" % df) # Load an EPI image - img = nipy.load_image(df) + img = nibabel.load(df) xyz = _check_xyz(xyz, img) data_list.append(np.array(img.get_fdata().flatten(), dtype=np.float64)) @@ -49,7 +49,7 @@ def _check_xyz(xyz, img): """Check voxel xyz consistency.""" # Use the image shape metadata; avoid get_fdata() here, which would read the # whole volume just to check consistency (load_epi reads the data later). - xyz_current = _get_xyz(img.coordmap.affine, img.shape[:3]) + xyz_current = _get_xyz(img.affine, img.shape[:3]) if xyz.size == 0: xyz = xyz_current diff --git a/bdpy/mri/load_mri.py b/bdpy/mri/load_mri.py index 772d6b0..94c31b1 100644 --- a/bdpy/mri/load_mri.py +++ b/bdpy/mri/load_mri.py @@ -1,7 +1,7 @@ """load_mri""" -import nipy +import nibabel import numpy as np @@ -13,17 +13,17 @@ def load_mri(fpath): - Returns voxel ijk indexes (3 x voxel) - Data, xyz, and ijk are flattened by Fortran-like index order """ - img = nipy.load_image(fpath) + img = nibabel.load(fpath) data = img.get_fdata() if data.ndim == 4: data = data.reshape(-1, data.shape[-1], order='F').T i_len, j_len, k_len, t = img.shape - affine = np.delete(np.delete(img.coordmap.affine, 3, axis=0), 3, axis=1) + affine = img.affine elif data.ndim == 3: data = data.flatten(order='F') i_len, j_len, k_len = img.shape - affine = img.coordmap.affine + affine = img.affine else: raise ValueError('Invalid shape.') diff --git a/tests/test_mri.py b/tests/test_mri.py index bee53ae..e1e27fb 100644 --- a/tests/test_mri.py +++ b/tests/test_mri.py @@ -47,6 +47,94 @@ def test_load_mri_3d(self) -> None: self.assertEqual(xyz.shape, (3, arr.size)) self.assertEqual(ijk.shape, (3, arr.size)) + def test_load_mri_4d(self) -> None: + """Test load_mri on a 4D volume. + + Guards a bug in the 4D path where the affine was reduced to 3x3 before + being applied to homogeneous (4 x N) voxel indices, raising a shape + mismatch. The 4D path must use the full 4x4 affine, just like the 3D + path, and return one row per time point. + """ + import os + import tempfile + + import nibabel + + i_len, j_len, k_len, t_len = 2, 3, 4, 5 + n_vox = i_len * j_len * k_len + arr = np.arange(n_vox * t_len, dtype=np.float32).reshape( + i_len, j_len, k_len, t_len) + affine = np.array([[2., 0., 0., -10.], + [0., 2., 0., -12.], + [0., 0., 2., -8.], + [0., 0., 0., 1.]]) + + with tempfile.TemporaryDirectory() as tmpdir: + fpath = os.path.join(tmpdir, 'test4d.nii.gz') + nibabel.save(nibabel.Nifti1Image(arr, affine=affine), fpath) + data, xyz, ijk = bmr.load_mri(fpath) # type: ignore + + # Data is (sample x voxel): one sample per time point, voxels flattened + # in Fortran order. + self.assertEqual(data.shape, (t_len, n_vox)) + np.testing.assert_array_equal( + data, arr.reshape(-1, t_len, order='F').T) + + # ijk are the voxel indices in Fortran order. + expected_ijk = np.array(np.unravel_index( + np.arange(n_vox), (i_len, j_len, k_len), order='F')) + np.testing.assert_array_equal(ijk, expected_ijk) + + # xyz are the world coordinates from the full 4x4 affine. + ijk_b = np.vstack([expected_ijk, np.ones((1, n_vox))]) + expected_xyz = np.dot(affine, ijk_b)[:3] + np.testing.assert_allclose(xyz, expected_xyz) + self.assertEqual(xyz.shape, (3, n_vox)) + + def test_braindata_load_volume_4d(self) -> None: + """Test BrainData volume loading on a 4D volume. + + Guards the same 4D affine bug as test_load_mri_4d in the duplicated + loader inside fmriprep.BrainData.__load_volume. The surrounding + create_bdata_fmriprep code treats ``data.shape[0]`` as the number of + volumes (time points), so the 4D path must return ``(T, N)`` data and + spatial xyz from the full 4x4 affine. + """ + import os + import tempfile + + import nibabel + + from bdpy.mri.fmriprep import BrainData + + i_len, j_len, k_len, t_len = 2, 3, 4, 5 + n_vox = i_len * j_len * k_len + arr = np.arange(n_vox * t_len, dtype=np.float32).reshape( + i_len, j_len, k_len, t_len) + affine = np.array([[2., 0., 0., -10.], + [0., 2., 0., -12.], + [0., 0., 2., -8.], + [0., 0., 0., 1.]]) + + with tempfile.TemporaryDirectory() as tmpdir: + fpath = os.path.join(tmpdir, 'bold4d.nii.gz') + nibabel.save(nibabel.Nifti1Image(arr, affine=affine), fpath) + brain = BrainData(fpath, dtype='volume') + + # data.shape[0] must be the number of time points (volumes). + self.assertEqual(brain.data.shape, (t_len, n_vox)) + np.testing.assert_array_equal( + brain.data, arr.reshape(-1, t_len, order='F').T) + + expected_ijk = np.array(np.unravel_index( + np.arange(n_vox), (i_len, j_len, k_len), order='F')) + np.testing.assert_array_equal(brain.index, expected_ijk) + + ijk_b = np.vstack([expected_ijk, np.ones((1, n_vox))]) + expected_xyz = np.dot(affine, ijk_b)[:3] + np.testing.assert_allclose(brain.xyz, expected_xyz) + self.assertEqual(brain.xyz.shape, (3, n_vox)) + def test_get_roiflag_pass0002(self) -> None: """Test for get_roiflag (pass case 0002).""" roi_xyz = [np.array([[1, 2, 3],