refactor(mri): replace nipy image I/O with nibabel and deprecate nipy GLM paradigm helper

bdpy/mri/fmriprep.py

@@ -11,7 +11,6 @@
 from collections import OrderedDict
 
 import nibabel
-import nipy
 import numpy as np
 import pandas as pd
 
@@ -95,13 +94,13 @@
         # File name patterns
         # FIXME
         if self.__fmriprep_version == '1.2':
             file_pattern = {'volume_native'   : '.*_space-T1w_desc-preproc_bold\.nii\.gz$',
                             'volume_standard' : '.*_space-MNI152NLin2009cAsym_desc-preproc_bold\.nii\.gz$',
                             'surf_native_left'  : '.*_space-fsnative_hemi-L\.func\.gii$',
                             'surf_native_right' : '.*_space-fsnative_hemi-R\.func\.gii$',
                             'surf_standard_left'  : '.*_space-fsaverage_hemi-L\.func\.gii$',
                             'surf_standard_right' : '.*_space-fsaverage_hemi-R\.func\.gii$',
                             'surf_standard_41k_left'  : '.*_space-fsaverage6_hemi-L\.func\.gii$',
                             'surf_standard_41k_right' : '.*_space-fsaverage6_hemi-R\.func\.gii$',
                             'surf_standard_10k_left'  : '.*_space-fsaverage5_hemi-L\.func\.gii$',
                             'surf_standard_10k_right' : '.*_space-fsaverage5_hemi-R\.func\.gii$',
@@ -406,17 +405,17 @@
         - 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 @@
     - 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.')
 

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 = []

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

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.')
 

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],