diff --git a/docs/builtin.rst b/docs/builtin.rst index 574b645f3..9c865f0c7 100644 --- a/docs/builtin.rst +++ b/docs/builtin.rst @@ -247,6 +247,7 @@ Available * - Name - Options - Keys + - Spaces - Version added - Publication * - Schaefer @@ -255,6 +256,7 @@ Available | ``Schaefer200x17``, ``Schaefer300x17``, ``Schaefer400x17``, | ``Schaefer500x17``, ``Schaefer600x17``, ``Schaefer700x17``, | ``Schaefer800x17``, ``Schaefer900x17``, ``Schaefer1000x17`` + - ``MNI152NLin6Asym`` - 0.0.1 - | Schaefer, A., Kong, R., Gordon, E.M. et al. | Local-Global Parcellation of the Human Cerebral Cortex from @@ -264,6 +266,7 @@ Available * - SUIT - ``space`` - ``SUITxMNI``, ``SUITxSUIT`` + - ``SUIT``, ``MNI152Lin6Asym`` - 0.0.1 - | Diedrichsen, J. | A spatially unbiased atlas template of the human cerebellum. @@ -279,6 +282,7 @@ Available | ``TianxS2x3TxMNInonlinear2009cAsym``, | ``TianxS3x3TxMNInonlinear2009cAsym``, | ``TianxS4x3TxMNInonlinear2009cAsym`` + - ``MNI152NLin6Asym``, ``MNI152NLin2009cAsym`` - 0.0.1 - | Tian, Y., Margulies, D.S., Breakspear, M. et al. | Topographic organization of the human subcortex @@ -288,6 +292,7 @@ Available * - AICHA - ``version`` - ``AICHA_v1``, ``AICHA_v2`` + - ``MNI152Lin6Asym`` - 0.0.3 - | Joliot, M., Jobard, G., Naveau, M. et al. | AICHA: An atlas of intrinsic connectivity of homotopic areas. @@ -297,6 +302,7 @@ Available - ``year``, ``n_rois`` - | ``Shen_2013_50``, ``Shen_2013_100``, ``Shen_2013_150``, | ``Shen_2015_268``, ``Shen_2019_368`` + - ``MNI152NLin2009cAsym`` - 0.0.3 - | Shen, X., Tokoglu, F., Papademetris, X., Constable, R.T. | Groupwise whole-brain parcellation from resting-state fMRI data @@ -322,6 +328,7 @@ Available | ``Yan400xKong17``, ``Yan500xKong17``, ``Yan600xKong17``, | ``Yan700xKong17``, ``Yan800xKong17``, ``Yan900xKong17``, | ``Yan1000xKong17`` + - ``MNI152NLin6Asym`` - 0.0.3 - | Yan, X., Kong, R., Xue, A., et al. | Homotopic local-global parcellation of the human cerebral cortex from @@ -590,28 +597,33 @@ Available * - Name - Keys + - Spaces - Version added - Description - Publication * - Vickery-Patil (Gray Matter) - | ``GM_prob0.2`` + - ``MNI152Lin6Asym`` - 0.0.1 - | Vickery, Sam, & Patil, Kaustubh. (2022). | Chimpanzee and Human Gray Matter Masks [Data set]. Zenodo. | https://doi.org/10.5281/zenodo.6463123 * - Vickery-Patil (Cortex + Basal Ganglia) - | ``GM_prob0.2_cortex`` + - ``MNI152Lin6Asym`` - 0.0.1 - | Vickery, Sam, & Patil, Kaustubh. (2022). | Chimpanzee and Human Gray Matter Masks [Data set]. Zenodo. | https://doi.org/10.5281/zenodo.6463123 * - Nilearn's MNI152 1mm-resolution mask - | ``compute_brain_mask`` + - Adapts to the target data - 0.0.2 - | Compute the whole-brain mask. This mask is calculated using | MNI152 1mm-resolution template mask onto the target image. | See :func:`nilearn.masking.compute_brain_mask` * - Nilearn's mask computed from FMRI data - | ``compute_epi_mask`` + - Adapts to the target data - 0.0.2 - | Compute a brain mask from fMRI data. This is based on an heuristic | proposed by T.Nichols: find the least dense point of the histogram, @@ -619,6 +631,7 @@ Available | image histogram. See :func:`nilearn.masking.compute_epi_mask` * - Nilearn's background mask - | ``compute_background_mask`` + - Adapts to the target data - 0.0.2 - | Compute a brain mask for the images by guessing the value of the | background from the border of the image. @@ -626,6 +639,7 @@ Available * - Nilearn's ICBM152 template gray-matter mask - | ``fetch_icbm152_brain_gm_mask`` + - ``MNI152NLin2009aAsym`` - 0.0.2 - | Compute a gray-matter mask from the asymmetrical ICBM152 2009 template, | release a. diff --git a/docs/changes/newsfragments/268.change b/docs/changes/newsfragments/268.change new file mode 100644 index 000000000..78a2f2030 --- /dev/null +++ b/docs/changes/newsfragments/268.change @@ -0,0 +1 @@ +Add ``space`` parameter to :func:`.register_coordinates`, :func:`.register_parcellation` and :func:`.register_mask` and return space from :func:`.load_coordinates`, :func:`.load_parcellation` and :func:`.load_mask` by `Synchon Mandal`_ and `Fede Raimondo`_ diff --git a/docs/changes/newsfragments/268.enh b/docs/changes/newsfragments/268.enh new file mode 100644 index 000000000..f45544395 --- /dev/null +++ b/docs/changes/newsfragments/268.enh @@ -0,0 +1 @@ +Add ``space`` information to existing datagrabbers, masks, parcellations and coordinates by `Synchon Mandal`_ and `Fede Raimondo`_ diff --git a/docs/extending/coordinates.rst b/docs/extending/coordinates.rst index 7f1580ee4..e9ca8b38a 100644 --- a/docs/extending/coordinates.rst +++ b/docs/extending/coordinates.rst @@ -20,8 +20,8 @@ example if you came up with your own set of coordinates), then junifer provides an easy way for you to register them using the :func:`.register_coordinates` function, so you can use your own set of coordinates within a junifer pipeline. -From the API reference, we can see that it has 3 positional arguments -(``name``, ``coordinates``, and ``voi_names``) as well as one +From the API reference, we can see that it has 4 positional arguments +(``name``, ``coordinates``, ``voi_names`` and ``space``) as well as one optional keyword argument (``overwrite``). The ``name`` argument takes a string indicating the name you want to give to @@ -40,12 +40,15 @@ belong to this set. Note, that junifer (as of yet) only works in MNI space, and so therefore these coordinates should always be real-world coordinates of the MNI space. -Lastly, the ``voi_names`` argument takes a list of strings +The ``voi_names`` argument takes a list of strings indicating the names of each coordinate (i.e. volume-of-interest) in the ``coordinates`` array. Therefore, the length of this list should correspond to the number of rows in the coordinates array. Now, we know everything we need to know to register a set of coordinates. +Lastly, we specify the ``space`` that the coordinates are in, for example, +``"MNI"`` or ``"Native"`` (scanner-native space). + Step 1: Prepare code to register a set of coordinates ----------------------------------------------------- @@ -90,7 +93,8 @@ simply use this to register our coordinates: register_coordinates( name="DMNCustom", coordinates=dmn_coords, - voi_names=voi_names + voi_names=voi_names, + space="MNI" ) Now, when we run this script, junifer registers these coordinates and we can diff --git a/docs/extending/masks.rst b/docs/extending/masks.rst index 1950941eb..02657b098 100644 --- a/docs/extending/masks.rst +++ b/docs/extending/masks.rst @@ -15,9 +15,9 @@ learn how to use your own masks. The principle is fairly simple and quite similar to :ref:`adding_parcellations` and :ref:`adding_coordinates`. junifer provides a :func:`.register_mask` -function that lets you register your own custom masks. It consists of two -positional arguments (``name`` and ``mask_path``) and one optional keyword -argument (``overwrite``). +function that lets you register your own custom masks. It consists of three +positional arguments (``name``, ``mask_path`` and ``space``) and one optional +keyword argument (``overwrite``). The ``name`` argument is a string indicating the name of the mask. This name is used to refer to that mask in junifer internally in order to obtain the @@ -29,6 +29,9 @@ The ``mask_path`` should contain the path to a valid NIfTI image with binary voxel values (i.e. 0 or 1). This data can then be used by junifer to mask other MR images. +Lastly, we specify the ``space`` that the coordinates are in, for example, +``"MNI"`` or ``"Native"`` (scanner-native space). + Step 1: Prepare code to register a mask --------------------------------------- @@ -46,7 +49,7 @@ look as follows: # on your system: mask_path = Path("..") / ".." / "my_custom_mask.nii.gz" - register_mask(name="my_custom_mask", mask_path=mask_path) + register_mask(name="my_custom_mask", mask_path=mask_path, space="Native") Simple, right? Now we just have to configure a YAML file to register this mask so we can use it for :ref:`codeless configuration of junifer `. diff --git a/docs/extending/parcellations.rst b/docs/extending/parcellations.rst index 888022097..b059aa781 100644 --- a/docs/extending/parcellations.rst +++ b/docs/extending/parcellations.rst @@ -20,9 +20,9 @@ the easy-to-use :func:`.register_parcellation` function to do just that. Let's try to understand the API reference and then use this function to register our own parcellation. -From the API reference, we can see that it has 3 positional arguments -(``name``, ``parcellation_path``, and ``parcels_labels``) as well as one -optional keyword argument (``overwrite``). +From the API reference, we can see that it has 4 positional arguments +(``name``, ``parcellation_path``, ``parcels_labels`` and ``space``) as well as +one optional keyword argument (``overwrite``). The ``name`` of the parcellation is up to you and will be the name that junifer will use to refer to this particular parcellation. You can think of this as @@ -40,7 +40,7 @@ this parcellation should be indicated by 0, and the labels of ROIs should go from 1 to N (where N is the total number of ROIs in your parcellation). Now, we nearly have everything we need. -Lastly, we also want to make sure that we can associate each integer label with +We also want to make sure that we can associate each integer label with a human readable name (i.e. the name for each ROI). This serves naming the features that parcellation-based markers produce in an unambiguous way, such that a user can easily identify which ROIs were used to produce a specific @@ -51,6 +51,9 @@ list, the label at the i-th position indicates the i-th integer label (i.e. the first label in this list corresponds to the first integer label in the parcellation and so on). +Lastly, we specify the ``space`` that the parcellation is in, for example, +``"MNI"`` or ``"Native"`` (scanner-native space). + Step 1: Prepare code to register a parcellation ----------------------------------------------- @@ -78,7 +81,8 @@ a simple example could look like this: register_parcellation( name="my_custom_parcellation", parcellation_path=path_to_parcellation, - parcels_labels=my_labels + parcels_labels=my_labels, + space="MNI" ) We can run this code and it seems to work, however, how can we actually diff --git a/examples/run_compute_parcel_mean.py b/examples/run_compute_parcel_mean.py index d4d00f046..0a82dab96 100644 --- a/examples/run_compute_parcel_mean.py +++ b/examples/run_compute_parcel_mean.py @@ -2,18 +2,21 @@ Computer Parcel Aggregation. ============================ -This example uses a ParcelAggregation marker to compute the mean of each parcel -using the Schaefer parcellations (100 rois, 7 Yeo networks) for both a 3D and -4D nifti +This example uses the ``ParcelAggregation`` marker to compute the mean of each +parcel using the Schaefer parcellations (100 rois, 7 Yeo networks) for both 3D +and 4D NIfTI. -Authors: Federico Raimondo +Authors: Federico Raimondo, Synchon Mandal License: BSD 3 clause """ -import nilearn - -from junifer.markers.parcel_aggregation import ParcelAggregation +from junifer.testing.datagrabbers import ( + OasisVBMTestingDataGrabber, + SPMAuditoryTestingDataGrabber, +) +from junifer.datareader import DefaultDataReader +from junifer.markers import ParcelAggregation from junifer.utils import configure_logging @@ -21,41 +24,34 @@ from junifer.utils import configure_logging # Set the logging level to info to see extra information configure_logging(level="INFO") +############################################################################### +# Perform parcel aggregation on VBM GM data (3D) from OASIS dataset +with OasisVBMTestingDataGrabber() as dg: + # Get the first element + element = dg.get_elements()[0] + # Read the element + element_data = DefaultDataReader().fit_transform(dg[element]) + # Initialize marker + marker = ParcelAggregation(parcellation="Schaefer100x7", method="mean") + # Compute feature + feature = marker.fit_transform(element_data) + # Print the output + print(feature.keys()) + print(feature["VBM_GM"]["data"].shape) # Shape is (1 x parcels) ############################################################################### -# Load the VBM GM data (3d): -# - Fetch the Oasis dataset -oasis_dataset = nilearn.datasets.fetch_oasis_vbm(n_subjects=1) -vbm_fname = oasis_dataset.gray_matter_maps[0] -vbm_img = nilearn.image.load_img(vbm_fname) - -############################################################################### -# Load the functional data (4d): -# - Fetch the SPM auditory dataset -# - Concatenate the functional data into one 4D image -s_func_data = nilearn.datasets.fetch_spm_auditory() -fmri_img = nilearn.image.concat_imgs(s_func_data.func) - -############################################################################### -# Define the marker -marker = ParcelAggregation(parcellation="Schaefer100x7", method="mean") - -############################################################################### -# Prepare the input -input = { - "BOLD": {"data": fmri_img, "meta": {"element": "subject1"}}, - "VBM_GM": {"data": vbm_img, "meta": {"element": "subject1"}}, -} - -############################################################################### -# Fit transform the data -out = marker.fit_transform(input) - -############################################################################### -# Check the results - -print(out.keys()) -print(out["VBM_GM"]["data"].shape) # Shape is (1 x parcels) - -print(out.keys()) -print(out["BOLD"]["data"].shape) # Shape is (timepoints x parcels) +# Perform parcel aggregation on BOLD data (4D) from SPM Auditory dataset +with SPMAuditoryTestingDataGrabber() as dg: + # Get the first element + element = dg.get_elements()[0] + # Read the element + element_data = DefaultDataReader().fit_transform(dg[element]) + # Initialize marker + marker = ParcelAggregation( + parcellation="Schaefer100x7", method="mean", on="BOLD" + ) + # Compute feature + feature = marker.fit_transform(element_data) + # Print the output + print(feature.keys()) + print(feature["BOLD"]["data"].shape) # Shape is (timepoints x parcels) diff --git a/junifer/data/coordinates.py b/junifer/data/coordinates.py index c9acc55e9..9bd903fa2 100644 --- a/junifer/data/coordinates.py +++ b/junifer/data/coordinates.py @@ -4,6 +4,7 @@ # Synchon Mandal # License: AGPL +import subprocess import typing from pathlib import Path from typing import Any, Dict, List, Optional, Tuple, Union @@ -12,6 +13,7 @@ import numpy as np import pandas as pd from numpy.typing import ArrayLike +from ..pipeline import WorkDirManager from ..utils.logging import logger, raise_error, warn_with_log @@ -24,31 +26,83 @@ _vois_meta_path = _vois_path / "meta" # A dictionary containing all supported coordinates and their respective file # or data. +# Each entry is a dictionary that must contain at least the following keys: +# * 'space': the coordinates' space (e.g., 'MNI') + # The built-in coordinates are files that are shipped with the package in the # data/VOIs directory. The user can also register their own coordinates, which # will be stored as numpy arrays in the dictionary. _available_coordinates: Dict[ - str, Union[Path, Dict[str, Union[ArrayLike, List[str]]]] + str, Dict[str, Union[Path, ArrayLike, List[str]]] ] = { - "CogAC": _vois_meta_path / "CogAC_VOIs.txt", - "CogAR": _vois_meta_path / "CogAR_VOIs.txt", - "DMNBuckner": _vois_meta_path / "DMNBuckner_VOIs.txt", - "eMDN": _vois_meta_path / "eMDN_VOIs.txt", - "Empathy": _vois_meta_path / "Empathy_VOIs.txt", - "eSAD": _vois_meta_path / "eSAD_VOIs.txt", - "extDMN": _vois_meta_path / "extDMN_VOIs.txt", - "Motor": _vois_meta_path / "Motor_VOIs.txt", - "MultiTask": _vois_meta_path / "MultiTask_VOIs.txt", - "PhysioStress": _vois_meta_path / "PhysioStress_VOIs.txt", - "Rew": _vois_meta_path / "Rew_VOIs.txt", - "Somatosensory": _vois_meta_path / "Somatosensory_VOIs.txt", - "ToM": _vois_meta_path / "ToM_VOIs.txt", - "VigAtt": _vois_meta_path / "VigAtt_VOIs.txt", - "WM": _vois_meta_path / "WM_VOIs.txt", - "Power": _vois_meta_path / "Power2011_MNI_VOIs.txt", - "Power2011": _vois_meta_path / "Power2011_MNI_VOIs.txt", - "Dosenbach": _vois_meta_path / "Dosenbach2010_MNI_VOIs.txt", - "Power2013": _vois_meta_path / "Power2013_MNI_VOIs.tsv", + "CogAC": { + "path": _vois_meta_path / "CogAC_VOIs.txt", + "space": "MNI", + }, + "CogAR": { + "path": _vois_meta_path / "CogAR_VOIs.txt", + "space": "MNI", + }, + "DMNBuckner": { + "path": _vois_meta_path / "DMNBuckner_VOIs.txt", + "space": "MNI", + }, + "eMDN": { + "path": _vois_meta_path / "eMDN_VOIs.txt", + "space": "MNI", + }, + "Empathy": { + "path": _vois_meta_path / "Empathy_VOIs.txt", + "space": "MNI", + }, + "eSAD": { + "path": _vois_meta_path / "eSAD_VOIs.txt", + "space": "MNI", + }, + "extDMN": { + "path": _vois_meta_path / "extDMN_VOIs.txt", + "space": "MNI", + }, + "Motor": { + "path": _vois_meta_path / "Motor_VOIs.txt", + "space": "MNI", + }, + "MultiTask": { + "path": _vois_meta_path / "MultiTask_VOIs.txt", + "space": "MNI", + }, + "PhysioStress": { + "path": _vois_meta_path / "PhysioStress_VOIs.txt", + "space": "MNI", + }, + "Rew": { + "path": _vois_meta_path / "Rew_VOIs.txt", + "space": "MNI", + }, + "Somatosensory": { + "path": _vois_meta_path / "Somatosensory_VOIs.txt", + "space": "MNI", + }, + "ToM": { + "path": _vois_meta_path / "ToM_VOIs.txt", + "space": "MNI", + }, + "VigAtt": { + "path": _vois_meta_path / "VigAtt_VOIs.txt", + "space": "MNI", + }, + "WM": { + "path": _vois_meta_path / "WM_VOIs.txt", + "space": "MNI", + }, + "Power": { + "path": _vois_meta_path / "Power2011_MNI_VOIs.txt", + "space": "MNI", + }, + "Dosenbach": { + "path": _vois_meta_path / "Dosenbach2010_MNI_VOIs.txt", + "space": "MNI", + }, } @@ -56,6 +110,7 @@ def register_coordinates( name: str, coordinates: ArrayLike, voi_names: List[str], + space: str, overwrite: Optional[bool] = False, ) -> None: """Register a custom user coordinates. @@ -71,6 +126,8 @@ def register_coordinates( z-coordinates). voi_names : list of str The names of the VOIs. + space : str + The space of the coordinates. overwrite : bool, optional If True, overwrite an existing list of coordinates with the same name. Does not apply to built-in coordinates (default False). @@ -88,7 +145,7 @@ def register_coordinates( """ if name in _available_coordinates: - if isinstance(_available_coordinates[name], Path): + if isinstance(_available_coordinates[name].get("path"), Path): raise_error( f"Coordinates {name} already registered as built-in " "coordinates." @@ -122,6 +179,7 @@ def register_coordinates( _available_coordinates[name] = { "coords": coordinates, "voi_names": voi_names, + "space": space, } @@ -166,16 +224,79 @@ def get_coordinates( Raises ------ ValueError - If ``extra_input`` is None when ``target_data``'s space is not MNI. + If ``extra_input`` is None when ``target_data``'s space is native. """ # Load the coordinates - seeds, labels = load_coordinates(name=coords) + seeds, labels, _ = load_coordinates(name=coords) + + # Transform coordinate if target data is native + if target_data["space"] == "native": + # Check for extra inputs + if extra_input is None: + raise_error( + "No extra input provided, requires `Warp` and `T1w` " + "data types in particular for transformation to " + f"{target_data['space']} space for further computation." + ) + + # Create tempdir + tempdir = WorkDirManager().get_tempdir(prefix="coordinates") + + # Save existing coordinates + pretransform_coordinates_path = ( + tempdir / "pretransform_coordinates.txt" + ) + np.savetxt(pretransform_coordinates_path, seeds) + + # Create a tempfile for transformed coordinates output + std2imgcoord_out_path = tempdir / "coordinates_transformed.txt" + # Set std2imgcoord command + std2imgcoord_cmd = [ + "std2imgcoord", + f"-img {target_data['reference_path'].resolve()}", + f"-warp {extra_input['Warp']['path'].resolve()}", + f"{pretransform_coordinates_path}", + f"> {std2imgcoord_out_path}", + ] + # Call std2imgcoord + std2imgcoord_cmd_str = " ".join(std2imgcoord_cmd) + logger.info( + f"std2imgcoord command to be executed: {std2imgcoord_cmd_str}" + ) + std2imgcoord_process = subprocess.run( + std2imgcoord_cmd_str, # string needed with shell=True + stdin=subprocess.DEVNULL, + stdout=subprocess.PIPE, + stderr=subprocess.STDOUT, + shell=True, # needed for respecting $PATH + check=False, + ) + # Delete saved coordinates file + pretransform_coordinates_path.unlink() + # Check for success or failure + if std2imgcoord_process.returncode == 0: + logger.info( + "std2imgcoord succeeded with the following output: " + f"{std2imgcoord_process.stdout}" + ) + else: + raise_error( + msg="std2imgcoord failed with the following error: " + f"{std2imgcoord_process.stdout}", + klass=RuntimeError, + ) + + # Load coordinates + seeds = np.loadtxt(std2imgcoord_out_path) + + # Delete tempdir + WorkDirManager().delete_tempdir(tempdir) return seeds, labels -def load_coordinates(name: str) -> Tuple[ArrayLike, List[str]]: +def load_coordinates(name: str) -> Tuple[ArrayLike, List[str], str]: """Load coordinates. Parameters @@ -189,6 +310,8 @@ def load_coordinates(name: str) -> Tuple[ArrayLike, List[str]]: The coordinates. list of str The names of the VOIs. + str + The space of the coordinates. Raises ------ @@ -221,9 +344,9 @@ def load_coordinates(name: str) -> Tuple[ArrayLike, List[str]]: # Load coordinates t_coord = _available_coordinates[name] - if isinstance(t_coord, Path): + if isinstance(t_coord.get("path"), Path): # Load via pandas - df_coords = pd.read_csv(t_coord, sep="\t", header=None) + df_coords = pd.read_csv(t_coord["path"], sep="\t", header=None) coords = df_coords.iloc[:, [0, 1, 2]].to_numpy() names = list(df_coords.iloc[:, [3]].values[:, 0]) else: @@ -232,4 +355,4 @@ def load_coordinates(name: str) -> Tuple[ArrayLike, List[str]]: names = t_coord["voi_names"] names = typing.cast(List[str], names) - return coords, names + return coords, names, t_coord["space"] diff --git a/junifer/data/masks.py b/junifer/data/masks.py index 7f22e907b..2837130cc 100644 --- a/junifer/data/masks.py +++ b/junifer/data/masks.py @@ -3,6 +3,8 @@ # Authors: Federico Raimondo # License: AGPL +import subprocess +import typing from pathlib import Path from typing import ( TYPE_CHECKING, @@ -26,6 +28,7 @@ from nilearn.masking import ( intersect_masks, ) +from ..pipeline import WorkDirManager from ..utils.logging import logger, raise_error from .utils import closest_resolution @@ -33,7 +36,7 @@ from .utils import closest_resolution if TYPE_CHECKING: from nibabel import Nifti1Image -# Path to the VOIs +# Path to the masks _masks_path = Path(__file__).parent / "masks" @@ -55,6 +58,7 @@ def _fetch_icbm152_brain_gm_mask( ------- nibabel.Nifti1Image The resampled mask. + """ mask = fetch_icbm152_brain_gm_mask(**kwargs) mask = resample_to_img( @@ -66,29 +70,40 @@ def _fetch_icbm152_brain_gm_mask( # A dictionary containing all supported masks and their respective file or # data. +# Each entry is a dictionary that must contain at least the following keys: +# * 'family': the mask's family name (e.g., 'Vickery-Patil', 'Callable') +# * 'space': the mask's space (e.g., 'MNI', 'inherit') + # The built-in masks are files that are shipped with the package in the # data/masks directory. The user can also register their own masks. # Callable masks should be functions that take at least one parameter: # * `target_img`: the image to which the mask will be applied. _available_masks: Dict[str, Dict[str, Any]] = { - "GM_prob0.2": {"family": "Vickery-Patil"}, - "GM_prob0.2_cortex": {"family": "Vickery-Patil"}, + "GM_prob0.2": {"family": "Vickery-Patil", "space": "IXI549Space"}, + "GM_prob0.2_cortex": { + "family": "Vickery-Patil", + "space": "IXI549Space", + }, "compute_brain_mask": { "family": "Callable", "func": compute_brain_mask, + "space": "inherit", }, "compute_background_mask": { "family": "Callable", "func": compute_background_mask, + "space": "inherit", }, "compute_epi_mask": { "family": "Callable", "func": compute_epi_mask, + "space": "inherit", }, "fetch_icbm152_brain_gm_mask": { "family": "Callable", "func": _fetch_icbm152_brain_gm_mask, + "space": "MNI152NLin2009aAsym", }, } @@ -96,6 +111,7 @@ _available_masks: Dict[str, Dict[str, Any]] = { def register_mask( name: str, mask_path: Union[str, Path], + space: str, overwrite: bool = False, ) -> None: """Register a custom user mask. @@ -106,6 +122,8 @@ def register_mask( The name of the mask. mask_path : str or pathlib.Path The path to the mask file. + space : str + The space of the mask. overwrite : bool, optional If True, overwrite an existing mask with the same name. Does not apply to built-in mask (default False). @@ -115,6 +133,7 @@ def register_mask( ValueError If the mask name is already registered and overwrite is set to False or if the mask name is a built-in mask. + """ # Check for attempt of overwriting built-in parcellations if name in _available_masks: @@ -122,11 +141,11 @@ def register_mask( logger.info(f"Overwriting {name} mask") if _available_masks[name]["family"] != "CustomUserMask": raise_error( - f"Cannot overwrite {name} mask. " "It is a built-in mask." + f"Cannot overwrite {name} mask. It is a built-in mask." ) else: raise_error( - f"Mask {name} already registered. Set `overwrite=True`" + f"Mask {name} already registered. Set `overwrite=True` " "to update its value." ) # Convert str to Path @@ -136,6 +155,7 @@ def register_mask( _available_masks[name] = { "path": str(mask_path.absolute()), "family": "CustomUserMask", + "space": space, } @@ -146,11 +166,12 @@ def list_masks() -> List[str]: ------- list of str A list with all available masks names. + """ return sorted(_available_masks.keys()) -def get_mask( +def get_mask( # noqa: C901 masks: Union[str, Dict, List[Union[Dict, str]]], target_data: Dict[str, Any], extra_input: Optional[Dict[str, Any]] = None, @@ -173,6 +194,23 @@ def get_mask( ------- Nifti1Image The mask image. + + Raises + ------ + RuntimeError + If masks are in different spaces and they need to be intersected / + unionized. + ValueError + If extra key is provided in addition to mask name in ``masks`` or + if no mask is provided or + if ``masks = "inherit"`` but ``extra_input`` is None or ``mask_item`` + is None or ``mask_items``'s value is not in ``extra_input`` or + if callable parameters are passed to non-callable mask or + if multiple masks are provided and their spaces do not match or + if parameters are passed to :func:`nilearn.masking.intersect_masks` + when there is only one mask or + if ``extra_input`` is None when ``target_data``'s space is native. + """ # Get the min of the voxels sizes and use it as the resolution target_img = target_data["data"] @@ -211,6 +249,7 @@ def get_mask( raise_error("No mask was passed. At least one mask is required.") # Get all the masks all_masks = [] + all_spaces = [] for t_mask in true_masks: if isinstance(t_mask, dict): mask_name = next(iter(t_mask.keys())) @@ -219,46 +258,73 @@ def get_mask( mask_name = t_mask mask_params = None + # If mask is being inherited from previous steps like preprocessing if mask_name == "inherit": + # Requires extra input to be passed if extra_input is None: raise_error( "Cannot inherit mask from another data item " "because no extra data was passed." ) + # Missing inherited mask item if inherited_mask_item is None: raise_error( "Cannot inherit mask from another data item " "because no mask item was specified " "(missing `mask_item` key in the data object)." ) + # Missing inherited mask item in extra input if inherited_mask_item not in extra_input: raise_error( "Cannot inherit mask from another data item " f"because the item ({inherited_mask_item}) does not exist." ) mask_img = extra_input[inherited_mask_item]["data"] + # Starting with new mask else: - mask_object, _ = load_mask( + # Load mask + mask_object, _, mask_space = load_mask( mask_name, path_only=False, resolution=resolution ) + # If mask is callable like from nilearn if callable(mask_object): if mask_params is None: mask_params = {} mask_img = mask_object(target_img, **mask_params) - else: # Mask is a Nifti1Image + # Mask is a Nifti1Image + else: + # Mask params provided if mask_params is not None: + # Unused params raise_error( "Cannot pass callable params to a non-callable mask." ) + # Resample mask to target image mask_img = resample_to_img( mask_object, target_img, interpolation="nearest", copy=True, ) + all_spaces.append(mask_space) all_masks.append(mask_img) + + # Multiple masks, need intersection / union if len(all_masks) > 1: - mask_img = intersect_masks(all_masks, **intersect_params) + # Filter out "inherit" and make a set for spaces + filtered_spaces = set(filter(lambda x: x != "inherit", all_spaces)) + # Intersect / union of masks only if all masks are in the same space + if len(filtered_spaces) == 1: + mask_img = intersect_masks(all_masks, **intersect_params) + else: + raise_error( + msg=( + f"Masks are in different spaces: {filtered_spaces}, " + "unable to merge." + ), + klass=RuntimeError, + ) + # Single mask else: if len(intersect_params) > 0: # Yes, I'm this strict! @@ -268,6 +334,67 @@ def get_mask( ) mask_img = all_masks[0] + # Warp mask if target data is native + if target_data["space"] == "native": + # Check for extra inputs + if extra_input is None: + raise_error( + "No extra input provided, requires `Warp` and `T1w` " + "data types in particular for transformation to " + f"{target_data['space']} space for further computation." + ) + + # Create tempdir + tempdir = WorkDirManager().get_tempdir(prefix="masks") + + # Save mask image + prewarp_mask_path = tempdir / "prewarp_mask.nii.gz" + nib.save(mask_img, prewarp_mask_path) + + # Create a tempfile for warped output + applywarp_out_path = tempdir / "mask_warped.nii.gz" + # Set applywarp command + applywarp_cmd = [ + "applywarp", + "--interp=spline", + f"-i {prewarp_mask_path.resolve()}", + # use resampled reference + f"-r {target_data['reference_path'].resolve()}", + f"-w {extra_input['Warp']['path'].resolve()}", + f"-o {applywarp_out_path.resolve()}", + ] + # Call applywarp + applywarp_cmd_str = " ".join(applywarp_cmd) + logger.info(f"applywarp command to be executed: {applywarp_cmd_str}") + applywarp_process = subprocess.run( + applywarp_cmd_str, # string needed with shell=True + stdin=subprocess.DEVNULL, + stdout=subprocess.PIPE, + stderr=subprocess.STDOUT, + shell=True, # needed for respecting $PATH + check=False, + ) + # Check for success or failure + if applywarp_process.returncode == 0: + logger.info( + "applywarp succeeded with the following output: " + f"{applywarp_process.stdout}" + ) + else: + raise_error( + msg="applywarp failed with the following error: " + f"{applywarp_process.stdout}", + klass=RuntimeError, + ) + + # Load nifti + mask_img = nib.load(applywarp_out_path) + + # Delete tempdir + WorkDirManager().delete_tempdir(tempdir) + + # Type-cast to remove errors + mask_img = typing.cast("Nifti1Image", mask_img) return mask_img @@ -275,13 +402,14 @@ def load_mask( name: str, resolution: Optional[float] = None, path_only: bool = False, -) -> Tuple[Optional[Union["Nifti1Image", Callable]], Optional[Path]]: - """Load mask. +) -> Tuple[Optional[Union["Nifti1Image", Callable]], Optional[Path], str]: + """Load a mask. Parameters ---------- name : str - The name of the mask. + The name of the mask. Check valid options by calling + :func:`.list_masks`. resolution : float, optional The desired resolution of the mask to load. If it is not available, the closest resolution will be loaded. Preferably, use a @@ -296,16 +424,27 @@ def load_mask( Loaded mask image. pathlib.Path or None File path to the mask image. + str + The space of the mask. + + Raises + ------ + ValueError + If the ``name`` is invalid of if the mask family is invalid. + """ - mask_img = None + # Check for valid mask name if name not in _available_masks: raise_error( - f"Mask {name} not found. " f"Valid options are: {list_masks()}" + f"Mask {name} not found. Valid options are: {list_masks()}" ) + # Copy mask definition to avoid edits in original object mask_definition = _available_masks[name].copy() t_family = mask_definition.pop("family") + # Check if the mask family is custom or built-in + mask_img = None if t_family == "CustomUserMask": mask_fname = Path(mask_definition["path"]) elif t_family == "Vickery-Patil": @@ -316,12 +455,16 @@ def load_mask( else: raise_error(f"I don't know about the {t_family} mask family.") + # Load mask if mask_fname is not None: - logger.info(f"Loading mask {mask_fname.absolute()}") + logger.info(f"Loading mask {mask_fname.absolute()!s}") if path_only is False: + # Load via nibabel mask_img = nib.load(mask_fname) - return mask_img, mask_fname + # Type-cast to remove error + mask_img = typing.cast("Nifti1Image", mask_img) + return mask_img, mask_fname, mask_definition["space"] def _load_vickery_patil_mask( @@ -332,7 +475,7 @@ def _load_vickery_patil_mask( Parameters ---------- - name : str + name : {"GM_prob0.2", "GM_prob0.2_cortex"} The name of the mask. resolution : float, optional The desired resolution of the mask to load. If it is not @@ -344,6 +487,13 @@ def _load_vickery_patil_mask( ------- pathlib.Path File path to the mask image. + + Raises + ------ + ValueError + If ``name`` is invalid or if ``resolution`` is invalid for + ``name = "GM_prob0.2"``. + """ if name == "GM_prob0.2": available_resolutions = [1.5, 3.0] @@ -356,12 +506,14 @@ def _load_vickery_patil_mask( mask_fname = "CAT12_IXI555_MNI152_TMP_GS_GMprob0.2_clean.nii.gz" else: raise_error( - f"Cannot find a GM_prob0.2 mask for resolution {resolution}" + f"Cannot find a GM_prob0.2 mask of resolution {resolution}" ) elif name == "GM_prob0.2_cortex": mask_fname = "GMprob0.2_cortex_3mm_NA_rm.nii.gz" else: raise_error(f"Cannot find a Vickery-Patil mask called {name}") + + # Set path for masks mask_fname = _masks_path / "vickery-patil" / mask_fname return mask_fname diff --git a/junifer/data/parcellations.py b/junifer/data/parcellations.py index dc5b9a4c4..25b690691 100644 --- a/junifer/data/parcellations.py +++ b/junifer/data/parcellations.py @@ -7,6 +7,7 @@ import io import shutil +import subprocess import tarfile import tempfile import typing @@ -21,6 +22,7 @@ import requests from nilearn import datasets, image from requests.exceptions import ConnectionError, HTTPError, ReadTimeout +from ..pipeline import WorkDirManager from ..utils.logging import logger, raise_error, warn_with_log from .utils import closest_resolution @@ -34,6 +36,7 @@ if TYPE_CHECKING: # Each entry is a dictionary that must contain at least the following keys: # * 'family': the parcellation's family name (e.g. 'Schaefer', 'SUIT') +# * 'space': the parcellation's space (e.g., 'MNI', 'SUIT') # Optional keys: # * 'valid_resolutions': a list of valid resolutions for the parcellation @@ -42,7 +45,7 @@ if TYPE_CHECKING: # TODO: have separate dictionary for built-in _available_parcellations: Dict[str, Dict[Any, Any]] = { "SUITxSUIT": {"family": "SUIT", "space": "SUIT"}, - "SUITxMNI": {"family": "SUIT", "space": "MNI"}, + "SUITxMNI": {"family": "SUIT", "space": "MNI152NLin6Asym"}, } # Add Schaefer parcellation info @@ -53,6 +56,7 @@ for n_rois in range(100, 1001, 100): "family": "Schaefer", "n_rois": n_rois, "yeo_networks": t_net, + "space": "MNI152NLin6Asym", } # Add Tian parcellation info for scale in range(1, 5): @@ -61,27 +65,28 @@ for scale in range(1, 5): "family": "Tian", "scale": scale, "magneticfield": "7T", - "space": "MNI6thgeneration", + "space": "MNI152NLin6Asym", } t_name = f"TianxS{scale}x3TxMNI6thgeneration" _available_parcellations[t_name] = { "family": "Tian", "scale": scale, "magneticfield": "3T", - "space": "MNI6thgeneration", + "space": "MNI152NLin6Asym", } t_name = f"TianxS{scale}x3TxMNInonlinear2009cAsym" _available_parcellations[t_name] = { "family": "Tian", "scale": scale, "magneticfield": "3T", - "space": "MNInonlinear2009cAsym", + "space": "MNI152NLin2009cAsym", } # Add AICHA parcellation info for version in (1, 2): _available_parcellations[f"AICHA_v{version}"] = { "family": "AICHA", "version": version, + "space": "IXI549Space", } # Add Shen parcellation info for year in (2013, 2015, 2019): @@ -91,18 +96,21 @@ for year in (2013, 2015, 2019): "family": "Shen", "year": 2013, "n_rois": n_rois, + "space": "MNI152NLin2009cAsym", } elif year == 2015: _available_parcellations["Shen_2015_268"] = { "family": "Shen", "year": 2015, "n_rois": 268, + "space": "MNI152NLin2009cAsym", } elif year == 2019: _available_parcellations["Shen_2019_368"] = { "family": "Shen", "year": 2019, "n_rois": 368, + "space": "MNI152NLin2009cAsym", } # Add Yan parcellation info for n_rois in range(100, 1001, 100): @@ -112,12 +120,14 @@ for n_rois in range(100, 1001, 100): "family": "Yan", "n_rois": n_rois, "yeo_networks": yeo_network, + "space": "MNI152NLin6Asym", } # Add Kong networks _available_parcellations[f"Yan{n_rois}xKong17"] = { "family": "Yan", "n_rois": n_rois, "kong_networks": 17, + "space": "MNI152NLin6Asym", } @@ -125,6 +135,7 @@ def register_parcellation( name: str, parcellation_path: Union[str, Path], parcels_labels: List[str], + space: str, overwrite: bool = False, ) -> None: """Register a custom user parcellation. @@ -137,6 +148,8 @@ def register_parcellation( The path to the parcellation file. parcels_labels : list of str The list of labels for the parcellation. + space : str + The space of the parcellation. overwrite : bool, optional If True, overwrite an existing parcellation with the same name. Does not apply to built-in parcellations (default False). @@ -173,6 +186,7 @@ def register_parcellation( "path": str(parcellation_path.absolute()), "labels": parcels_labels, "family": "CustomUserParcellation", + "space": space, } @@ -214,6 +228,13 @@ def get_parcellation( list of str Parcellation labels. + Raises + ------ + RuntimeError + If parcellations are in different spaces and they need to be merged. + ValueError + If ``extra_input`` is None when ``target_data``'s space is native. + """ # Get the min of the voxels sizes and use it as the resolution target_img = target_data["data"] @@ -222,8 +243,9 @@ def get_parcellation( # Load the parcellations all_parcellations = [] all_labels = [] + all_spaces = [] for name in parcellation: - img, labels, _ = load_parcellation( + img, labels, _, space = load_parcellation( name=name, resolution=resolution, ) @@ -236,17 +258,89 @@ def get_parcellation( ) all_parcellations.append(resampled_img) all_labels.append(labels) + all_spaces.append(space) # Avoid merging if there is only one parcellation if len(all_parcellations) == 1: resampled_parcellation_img = all_parcellations[0] labels = all_labels[0] else: - # Merge the parcellations - resampled_parcellation_img, labels = merge_parcellations( - parcellations_list=all_parcellations, - parcellations_names=parcellation, - labels_lists=all_labels, + # Merge the parcellations only if all parcellations are in the same + # space + if len(set(all_spaces)) == 1: + resampled_parcellation_img, labels = merge_parcellations( + parcellations_list=all_parcellations, + parcellations_names=parcellation, + labels_lists=all_labels, + ) + else: + raise_error( + msg="Parcellations are in different spaces, unable to merge.", + klass=RuntimeError, + ) + + # Warp parcellation if target data is native + if target_data["space"] == "native": + # Check for extra inputs + if extra_input is None: + raise_error( + "No extra input provided, requires `Warp` and `T1w` " + "data types in particular for transformation to " + f"{target_data['space']} space for further computation." + ) + + # Create tempdir + tempdir = WorkDirManager().get_tempdir(prefix="parcellations") + + # Save parcellation image + prewarp_parcellation_path = tempdir / "prewarp_parcellation.nii.gz" + nib.save(resampled_parcellation_img, prewarp_parcellation_path) + + # Create a tempfile for warped output + applywarp_out_path = tempdir / "parcellation_warped.nii.gz" + # Set applywarp command + applywarp_cmd = [ + "applywarp", + "--interp=spline", + f"-i {prewarp_parcellation_path.resolve()}", + # use resampled reference + f"-r {target_data['reference_path'].resolve()}", + f"-w {extra_input['Warp']['path'].resolve()}", + f"-o {applywarp_out_path.resolve()}", + ] + # Call applywarp + applywarp_cmd_str = " ".join(applywarp_cmd) + logger.info(f"applywarp command to be executed: {applywarp_cmd_str}") + applywarp_process = subprocess.run( + applywarp_cmd_str, # string needed with shell=True + stdin=subprocess.DEVNULL, + stdout=subprocess.PIPE, + stderr=subprocess.STDOUT, + shell=True, # needed for respecting $PATH + check=False, + ) + # Check for success or failure + if applywarp_process.returncode == 0: + logger.info( + "applywarp succeeded with the following output: " + f"{applywarp_process.stdout}" + ) + else: + raise_error( + msg="applywarp failed with the following error: " + f"{applywarp_process.stdout}", + klass=RuntimeError, + ) + + # Load nifti + resampled_parcellation_img = nib.load(applywarp_out_path) + + # Delete tempdir + WorkDirManager().delete_tempdir(tempdir) + + # Stupid casting + resampled_parcellation_img = typing.cast( + "Nifti1Image", resampled_parcellation_img ) return resampled_parcellation_img, labels @@ -257,7 +351,7 @@ def load_parcellation( parcellations_dir: Union[str, Path, None] = None, resolution: Optional[float] = None, path_only: bool = False, -) -> Tuple[Optional["Nifti1Image"], List[str], Path]: +) -> Tuple[Optional["Nifti1Image"], List[str], Path, str]: """Load a brain parcellation (including a label file). If it is a built-in parcellation and the file is not present in the @@ -287,6 +381,8 @@ def load_parcellation( Parcellation labels. pathlib.Path File path to the parcellation image. + str + The space of the parcellation. Raises ------ @@ -305,6 +401,11 @@ def load_parcellation( # Copy parcellation definition to avoid edits in original object parcellation_definition = _available_parcellations[name].copy() t_family = parcellation_definition.pop("family") + # Remove space conditionally + if t_family not in ["SUIT", "Tian"]: + space = parcellation_definition.pop("space") + else: + space = parcellation_definition["space"] # Check if the parcellation family is custom or built-in if t_family == "CustomUserParcellation": @@ -343,7 +444,7 @@ def load_parcellation( # Type-cast to remove errors parcellation_img = typing.cast("Nifti1Image", parcellation_img) - return parcellation_img, parcellation_labels, parcellation_fname + return parcellation_img, parcellation_labels, parcellation_fname, space def _retrieve_parcellation( @@ -384,13 +485,13 @@ def _retrieve_parcellation( * Tian : ``scale`` : {1, 2, 3, 4} Scale of parcellation (defines granularity). - ``space`` : {"MNI6thgeneration", "MNInonlinear2009cAsym"}, optional - Space of parcellation (default "MNI6thgeneration"). (For more + ``space`` : {"MNI152NLin6Asym", "MNI152NLin2009cAsym"}, optional + Space of parcellation (default "MNI152NLin6Asym"). (For more information see https://github.com/yetianmed/subcortex) ``magneticfield`` : {"3T", "7T"}, optional Magnetic field (default "3T"). * SUIT : - ``space`` : {"MNI", "SUIT"}, optional + ``space`` : {"MNI152NLin6Asym", "SUIT"}, optional Space of parcellation (default "MNI"). (For more information see http://www.diedrichsenlab.org/imaging/suit.htm). * AICHA : @@ -588,7 +689,7 @@ def _retrieve_tian( parcellations_dir: Path, resolution: Optional[float] = None, scale: Optional[int] = None, - space: str = "MNI6thgeneration", + space: str = "MNI152NLin6Asym", magneticfield: str = "3T", ) -> Tuple[Path, List[str]]: """Retrieve Tian parcellation. @@ -605,8 +706,8 @@ def _retrieve_tian( parcellation depend on the space and magnetic field. scale : {1, 2, 3, 4}, optional Scale of parcellation (defines granularity) (default None). - space : {"MNI6thgeneration", "MNInonlinear2009cAsym"}, optional - Space of parcellation (default "MNI6thgeneration"). (For more + space : {"MNI152NLin6Asym", "MNI152NLin2009cAsym"}, optional + Space of parcellation (default "MNI152NLin6Asym"). (For more information see https://github.com/yetianmed/subcortex) magneticfield : {"3T", "7T"}, optional Magnetic field (default "3T"). @@ -642,10 +743,10 @@ def _retrieve_tian( _valid_resolutions = [] # avoid pylance error if magneticfield == "3T": - _valid_spaces = ["MNI6thgeneration", "MNInonlinear2009cAsym"] - if space == "MNI6thgeneration": + _valid_spaces = ["MNI152NLin6Asym", "MNI152NLin2009cAsym"] + if space == "MNI152NLin6Asym": _valid_resolutions = [1, 2] - elif space == "MNInonlinear2009cAsym": + elif space == "MNI152NLin2009cAsym": _valid_resolutions = [2] else: raise_error( @@ -654,10 +755,10 @@ def _retrieve_tian( ) elif magneticfield == "7T": _valid_resolutions = [1.6] - if space != "MNI6thgeneration": + if space != "MNI152NLin6Asym": raise_error( f"The parameter `space` ({space}) for 7T needs to be " - f"MNI6thgeneration" + f"MNI152NLin6Asym" ) else: raise_error( @@ -675,7 +776,7 @@ def _retrieve_tian( parcellation_lname = parcellation_fname_base_3T / ( f"Tian_Subcortex_S{scale}_3T_label.txt" ) - if space == "MNI6thgeneration": + if space == "MNI152NLin6Asym": parcellation_fname = parcellation_fname_base_3T / ( f"Tian_Subcortex_S{scale}_{magneticfield}.nii.gz" ) @@ -684,7 +785,7 @@ def _retrieve_tian( parcellation_fname_base_3T / f"Tian_Subcortex_S{scale}_{magneticfield}_1mm.nii.gz" ) - elif space == "MNInonlinear2009cAsym": + elif space == "MNI152NLin2009cAsym": space = "2009cAsym" parcellation_fname = parcellation_fname_base_3T / ( f"Tian_Subcortex_S{scale}_{magneticfield}_{space}.nii.gz" @@ -754,7 +855,9 @@ def _retrieve_tian( def _retrieve_suit( - parcellations_dir: Path, resolution: Optional[float], space: str = "MNI" + parcellations_dir: Path, + resolution: Optional[float], + space: str = "MNI152NLin6Asym", ) -> Tuple[Path, List[str]]: """Retrieve SUIT parcellation. @@ -768,9 +871,9 @@ def _retrieve_suit( resolution higher than the desired one. By default, will load the highest one (default None). Available resolutions for this parcellation are 1mm and 2mm. - space : {"MNI", "SUIT"}, optional - Space of parcellation (default "MNI"). (For more information - see http://www.diedrichsenlab.org/imaging/suit.htm). + space : {"MNI152NLin6Asym", "SUIT"}, optional + Space of parcellation (default "MNI152NLin6Asym"). (For more + information see http://www.diedrichsenlab.org/imaging/suit.htm). Returns ------- @@ -790,7 +893,7 @@ def _retrieve_suit( logger.info(f"\tresolution: {resolution}") logger.info(f"\tspace: {space}") - _valid_spaces = ["MNI", "SUIT"] + _valid_spaces = ["MNI152NLin6Asym", "SUIT"] # check validity of parameters if space not in _valid_spaces: @@ -804,7 +907,11 @@ def _retrieve_suit( resolution = closest_resolution(resolution, _valid_resolutions) - # define file names + # Format space if MNI; required for the file name + if space == "MNI152NLin6Asym": + space = "MNI" + + # Define parcellation and label file names parcellation_fname = ( parcellations_dir / "SUIT" / (f"SUIT_{space}Space_{resolution}mm.nii") ) @@ -812,7 +919,7 @@ def _retrieve_suit( parcellations_dir / "SUIT" / (f"SUIT_{space}Space_{resolution}mm.tsv") ) - # check existence of parcellation + # Check existence of parcellation if not (parcellation_fname.exists() and parcellation_lname.exists()): parcellation_fname.parent.mkdir(exist_ok=True, parents=True) logger.info( @@ -827,6 +934,7 @@ def _retrieve_suit( url_SUIT = url_basis + "atl-Anatom_space-SUIT_dseg.nii" url_labels = url_basis + "atl-Anatom.tsv" + # TODO: improve HTTP request / response handling if space == "MNI": logger.info(f"Downloading {url_MNI}") parcellation_download = requests.get(url_MNI) @@ -891,12 +999,23 @@ def _retrieve_aicha( If invalid value is provided for ``version`` or if there is a problem fetching the parcellation. + Warns + ----- + RuntimeWarning + Until the authors confirm the space, the warning will be issued. + Notes ----- The resolution of the parcellation is 2mm and although v2 provides 1mm, it is only for display purpose as noted in the release document. """ + # Issue warning until space is confirmed by authors + warn_with_log( + "The current space for AICHA parcellations are IXI549Space, but are " + "not confirmed by authors, until that this warning will be issued." + ) + # show parameters to user logger.info("Parcellation parameters:") logger.info(f"\tresolution: {resolution}") diff --git a/junifer/data/tests/test_coordinates.py b/junifer/data/tests/test_coordinates.py index f762c2f79..8a1928a12 100644 --- a/junifer/data/tests/test_coordinates.py +++ b/junifer/data/tests/test_coordinates.py @@ -24,6 +24,7 @@ def test_register_coordinates_built_in_check() -> None: name="DMNBuckner", coordinates=np.zeros(2), voi_names=["1", "2"], + space="MNI", overwrite=True, ) @@ -34,6 +35,7 @@ def test_register_coordinates_overwrite() -> None: name="MyList", coordinates=np.zeros((2, 3)), voi_names=["roi1", "roi2"], + space="MNI", overwrite=True, ) with pytest.raises(ValueError, match=r"already registered"): @@ -41,18 +43,21 @@ def test_register_coordinates_overwrite() -> None: name="MyList", coordinates=np.ones((2, 3)), voi_names=["roi2", "roi3"], + space="MNI", ) register_coordinates( name="MyList", coordinates=np.ones((2, 3)), voi_names=["roi2", "roi3"], + space="MNI", overwrite=True, ) - coord, names = load_coordinates("MyList") + coord, names, space = load_coordinates("MyList") assert_array_equal(coord, np.ones((2, 3))) assert names == ["roi2", "roi3"] + assert space == "MNI" def test_register_coordinates_valid_input() -> None: @@ -62,6 +67,7 @@ def test_register_coordinates_valid_input() -> None: name="MyList", coordinates=[1, 2], voi_names=["roi1", "roi2"], + space="MNI", overwrite=True, ) with pytest.raises(ValueError, match=r"2D array"): @@ -69,6 +75,7 @@ def test_register_coordinates_valid_input() -> None: name="MyList", coordinates=np.zeros((2, 3, 4)), voi_names=["roi1", "roi2"], + space="MNI", overwrite=True, ) @@ -77,6 +84,7 @@ def test_register_coordinates_valid_input() -> None: name="MyList", coordinates=np.zeros((2, 4)), voi_names=["roi1", "roi2"], + space="MNI", overwrite=True, ) with pytest.raises(ValueError, match=r"voi_names"): @@ -84,6 +92,7 @@ def test_register_coordinates_valid_input() -> None: name="MyList", coordinates=np.zeros((2, 3)), voi_names=["roi1", "roi2", "roi3"], + space="MNI", overwrite=True, ) @@ -99,9 +108,10 @@ def test_list_coordinates() -> None: def test_load_coordinates() -> None: """Test loading coordinates from file.""" - coord, names = load_coordinates("DMNBuckner") + coord, names, space = load_coordinates("DMNBuckner") assert coord.shape == (6, 3) # type: ignore assert names == ["PCC", "MPFC", "lAG", "rAG", "lHF", "rHF"] + assert space == "MNI" def test_load_coordinates_nonexisting() -> None: @@ -122,7 +132,7 @@ def test_get_coordinates() -> None: coords="DMNBuckner", target_data=vbm_gm ) # Get raw coordinates - raw_coords, raw_labels = load_coordinates("DMNBuckner") + raw_coords, raw_labels, _ = load_coordinates("DMNBuckner") # Both tailored and raw should be same for now assert_array_equal(tailored_coords, raw_coords) assert tailored_labels == raw_labels diff --git a/junifer/data/tests/test_masks.py b/junifer/data/tests/test_masks.py index 85c4e252a..1a530f7ae 100644 --- a/junifer/data/tests/test_masks.py +++ b/junifer/data/tests/test_masks.py @@ -6,7 +6,7 @@ # License: AGPL from pathlib import Path -from typing import Callable, Dict, List, Union +from typing import Callable, Dict, List, Optional, Union import numpy as np import pytest @@ -41,6 +41,7 @@ def test_register_mask_built_in_check() -> None: register_mask( name="GM_prob0.2", mask_path="testmask.nii.gz", + space="MNI", overwrite=True, ) @@ -57,6 +58,7 @@ def test_register_mask_already_registered() -> None: register_mask( name="testmask", mask_path="testmask.nii.gz", + space="MNI", ) out = load_mask("testmask", path_only=True) assert out[1] is not None @@ -67,10 +69,12 @@ def test_register_mask_already_registered() -> None: register_mask( name="testmask", mask_path="testmask.nii.gz", + space="MNI", ) register_mask( name="testmask", mask_path="testmask2.nii.gz", + space="MNI", overwrite=True, ) @@ -80,16 +84,17 @@ def test_register_mask_already_registered() -> None: @pytest.mark.parametrize( - "name, mask_path, overwrite", + "name, mask_path, space, overwrite", [ - ("testmask_1", "testmask_1.nii.gz", True), - ("testmask_2", "testmask_2.nii.gz", True), - ("testmask_3", Path("testmask_3.nii.gz"), True), + ("testmask_1", "testmask_1.nii.gz", "MNI", True), + ("testmask_2", "testmask_2.nii.gz", "MNI", True), + ("testmask_3", Path("testmask_3.nii.gz"), "MNI", True), ], ) def test_register_mask( name: str, mask_path: str, + space: str, overwrite: bool, ) -> None: """Test mask registration. @@ -100,6 +105,8 @@ def test_register_mask( The parametrized mask name. mask_path : str or pathlib.Path The parametrized mask path. + space : str + The parametrized mask space. overwrite : bool The parametrized mask overwrite value. @@ -108,16 +115,18 @@ def test_register_mask( register_mask( name=name, mask_path=mask_path, + space=space, overwrite=overwrite, ) # List available mask and check registration masks = list_masks() assert name in masks # Load registered mask - _, fname = load_mask(name=name, path_only=True) + _, fname, mask_space = load_mask(name=name, path_only=True) # Check values for registered mask assert fname is not None assert fname.name == f"{name}.nii.gz" + assert space == mask_space @pytest.mark.parametrize( @@ -146,36 +155,62 @@ def test_load_mask_incorrect() -> None: load_mask("wrongmask") -def test_vickery_patil() -> None: - """Test Vickery-Patil mask.""" - mask, fname = load_mask("GM_prob0.2") +@pytest.mark.parametrize( + "name, resolution, pixdim, fname", + [ + ( + "GM_prob0.2", + None, + [1.5, 1.5, 1.5], + "CAT12_IXI555_MNI152_TMP_GS_GMprob0.2_clean.nii.gz", + ), + ( + "GM_prob0.2", + 3.0, + [3.0, 3.0, 3.0], + "CAT12_IXI555_MNI152_TMP_GS_GMprob0.2_clean_3mm.nii.gz", + ), + ( + "GM_prob0.2_cortex", + None, + [3.0, 3.0, 3.0], + "GMprob0.2_cortex_3mm_NA_rm.nii.gz", + ), + ], +) +def test_vickery_patil( + name: str, + resolution: Optional[float], + pixdim: List[float], + fname: str, +) -> None: + """Test Vickery-Patil mask. + + Parameters + ---------- + name : str + The parametrized name of the mask. + resolution : float or None + The parametrized resolution of the mask. + pixdim : list of float + The parametrized pixel dimensions of the mask. + fname : str + The parametrized name of the mask file. + + """ + mask, mask_fname, space = load_mask(name, resolution=resolution) assert_array_almost_equal( - mask.header["pixdim"][1:4], [1.5, 1.5, 1.5] # type: ignore + mask.header["pixdim"][1:4], pixdim # type: ignore ) + assert space == "IXI549Space" + assert mask_fname is not None + assert mask_fname.name == fname - assert fname is not None - assert fname.name == "CAT12_IXI555_MNI152_TMP_GS_GMprob0.2_clean.nii.gz" - - mask, fname = load_mask("GM_prob0.2", resolution=3) - assert_array_almost_equal( - mask.header["pixdim"][1:4], [3.0, 3.0, 3.0] # type: ignore - ) - - assert fname is not None - assert ( - fname.name == "CAT12_IXI555_MNI152_TMP_GS_GMprob0.2_clean_3mm.nii.gz" - ) - - mask, fname = load_mask("GM_prob0.2_cortex") - assert_array_almost_equal( - mask.header["pixdim"][1:4], [3.0, 3.0, 3.0] # type: ignore - ) - - assert fname is not None - assert fname.name == "GMprob0.2_cortex_3mm_NA_rm.nii.gz" +def test_vickery_patil_error() -> None: + """Test error for Vickery-Patil mask.""" with pytest.raises(ValueError, match=r"find a Vickery-Patil mask "): - _load_vickery_patil_mask("wrong", resolution=2) + _load_vickery_patil_mask(name="wrong", resolution=2.0) def test_get_mask() -> None: @@ -191,7 +226,7 @@ def test_get_mask() -> None: assert mask.shape == vbm_gm_img.shape assert_array_equal(mask.affine, vbm_gm_img.affine) - raw_mask_img, _ = load_mask("GM_prob0.2", resolution=1.5) + raw_mask_img, _, _ = load_mask("GM_prob0.2", resolution=1.5) res_mask_img = resample_to_img( raw_mask_img, vbm_gm_img, @@ -207,7 +242,11 @@ def test_mask_callable() -> None: def ident(x): return x - _available_masks["identity"] = {"family": "Callable", "func": ident} + _available_masks["identity"] = { + "family": "Callable", + "func": ident, + "space": "MNI", + } reader = DefaultDataReader() with OasisVBMTestingDataGrabber() as dg: input = dg["sub-01"] @@ -371,14 +410,6 @@ def test_get_mask_inherit() -> None: ["GM_prob0.2", "compute_brain_mask"], {"threshold": 0.2}, ), - ( - [ - "GM_prob0.2", - "compute_brain_mask", - "fetch_icbm152_brain_gm_mask", - ], - {"threshold": 1, "connected": True}, - ), ], ) def test_get_mask_multiple( @@ -445,3 +476,21 @@ def test_get_mask_multiple( expected = intersect_masks(mask_imgs, **params) assert_array_equal(computed.get_fdata(), expected.get_fdata()) + + +def test_get_mask_multiple_incorrect_space() -> None: + """Test incorrect space error for getting multiple masks.""" + reader = DefaultDataReader() + with SPMAuditoryTestingDataGrabber() as dg: + input = dg["sub001"] + input = reader.fit_transform(input) + + with pytest.raises(RuntimeError, match="unable to merge."): + get_mask( + masks=[ + "GM_prob0.2", + "compute_brain_mask", + "fetch_icbm152_brain_gm_mask", + ], + target_data=input["BOLD"], + ) diff --git a/junifer/data/tests/test_parcellations.py b/junifer/data/tests/test_parcellations.py index 7eae10aa1..be190f1f9 100644 --- a/junifer/data/tests/test_parcellations.py +++ b/junifer/data/tests/test_parcellations.py @@ -39,6 +39,7 @@ def test_register_parcellation_built_in_check() -> None: name="SUITxSUIT", parcellation_path="testparc.nii.gz", parcels_labels=["1", "2", "3"], + space="SUIT", overwrite=True, ) @@ -56,6 +57,7 @@ def test_register_parcellation_already_registered() -> None: name="testparc", parcellation_path="testparc.nii.gz", parcels_labels=["1", "2", "3"], + space="MNI", ) assert ( load_parcellation("testparc", path_only=True)[2].name @@ -68,11 +70,13 @@ def test_register_parcellation_already_registered() -> None: name="testparc", parcellation_path="testparc.nii.gz", parcels_labels=["1", "2", "3"], + space="MNI", ) register_parcellation( name="testparc", parcellation_path="testparc2.nii.gz", parcels_labels=["1", "2", "3"], + space="MNI", overwrite=True, ) @@ -90,17 +94,19 @@ def test_parcellation_wrong_labels_values(tmp_path: Path) -> None: tmp_path : pathlib.Path The path to the test directory. """ - schaefer, labels, schaefer_path = load_parcellation("Schaefer100x7") + schaefer, labels, schaefer_path, _ = load_parcellation("Schaefer100x7") assert schaefer is not None # Test wrong number of labels - register_parcellation("WrongLabels", schaefer_path, labels[:10]) + register_parcellation("WrongLabels", schaefer_path, labels[:10], "MNI") with pytest.raises(ValueError, match=r"has 100 parcels but 10"): load_parcellation("WrongLabels") # Test wrong number of labels - register_parcellation("WrongLabels2", schaefer_path, [*labels, "wrong"]) + register_parcellation( + "WrongLabels2", schaefer_path, [*labels, "wrong"], "MNI" + ) with pytest.raises(ValueError, match=r"has 100 parcels but 101"): load_parcellation("WrongLabels2") @@ -111,7 +117,7 @@ def test_parcellation_wrong_labels_values(tmp_path: Path) -> None: new_schaefer_img = new_img_like(schaefer, schaefer_data) nib.save(new_schaefer_img, new_schaefer_path) - register_parcellation("WrongValues", new_schaefer_path, labels[:-1]) + register_parcellation("WrongValues", new_schaefer_path, labels[:-1], "MNI") with pytest.raises(ValueError, match=r"the range [0, 99]"): load_parcellation("WrongValues") @@ -121,23 +127,30 @@ def test_parcellation_wrong_labels_values(tmp_path: Path) -> None: new_schaefer_img = new_img_like(schaefer, schaefer_data) nib.save(new_schaefer_img, new_schaefer_path) - register_parcellation("WrongValues2", new_schaefer_path, labels) + register_parcellation("WrongValues2", new_schaefer_path, labels, "MNI") with pytest.raises(ValueError, match=r"the range [0, 100]"): load_parcellation("WrongValues2") @pytest.mark.parametrize( - "name, parcellation_path, parcels_labels, overwrite", + "name, parcellation_path, parcels_labels, space, overwrite", [ - ("testparc_1", "testparc_1.nii.gz", ["1", "2", "3"], True), - ("testparc_2", "testparc_2.nii.gz", ["1", "2", "6"], True), - ("testparc_3", Path("testparc_3.nii.gz"), ["1", "2", "6"], True), + ("testparc_1", "testparc_1.nii.gz", ["1", "2", "3"], "MNI", True), + ("testparc_2", "testparc_2.nii.gz", ["1", "2", "6"], "MNI", True), + ( + "testparc_3", + Path("testparc_3.nii.gz"), + ["1", "2", "6"], + "MNI", + True, + ), ], ) def test_register_parcellation( name: str, parcellation_path: str, parcels_labels: List[str], + space: str, overwrite: bool, ) -> None: """Test parcellation registration. @@ -150,6 +163,8 @@ def test_register_parcellation( The parametrized parcellation path. parcels_labels : list of str The parametrized parcellation labels. + space : str + The parametrized parcellation space. overwrite : bool The parametrized parcellation overwrite value. @@ -159,16 +174,20 @@ def test_register_parcellation( name=name, parcellation_path=parcellation_path, parcels_labels=parcels_labels, + space=space, overwrite=overwrite, ) # List available parcellation and check registration parcellations = list_parcellations() assert name in parcellations # Load registered parcellation - _, lbl, fname = load_parcellation(name=name, path_only=True) + _, lbl, fname, parcellation_space = load_parcellation( + name=name, path_only=True + ) # Check values for registered parcellation assert lbl == parcels_labels assert fname.name == f"{name}.nii.gz" + assert parcellation_space == space @pytest.mark.parametrize( @@ -284,7 +303,7 @@ def test_schaefer( f"{int(resolution)}mm.nii.gz" ) # Load parcellation - img, label, img_path = load_parcellation( + img, label, img_path, space = load_parcellation( name=parcellation_name, parcellations_dir=tmp_path, resolution=resolution, @@ -292,6 +311,7 @@ def test_schaefer( assert img is not None assert img_path.name == parcellation_file assert len(label) == n_rois + assert space == "MNI152NLin6Asym" assert_array_equal( img.header["pixdim"][1:4], 3 * [resolution] # type: ignore ) @@ -334,29 +354,32 @@ def test_retrieve_schaefer_incorrect_yeo_networks(tmp_path: Path) -> None: @pytest.mark.parametrize( - "space", - ["SUIT", "MNI"], + "space_key, space", + [("SUIT", "SUIT"), ("MNI", "MNI152NLin6Asym")], ) -def test_suit(tmp_path: Path, space: str) -> None: +def test_suit(tmp_path: Path, space_key: str, space: str) -> None: """Test SUIT parcellation. Parameters ---------- tmp_path : pathlib.Path The path to the test directory. + space_key : str + The parametrized space values for the key. space : str The parametrized space values. """ parcellations = list_parcellations() - assert f"SUITx{space}" in parcellations + assert f"SUITx{space_key}" in parcellations # Load parcellation - img, label, img_path = load_parcellation( - name=f"SUITx{space}", + img, label, img_path, parcellation_space = load_parcellation( + name=f"SUITx{space_key}", parcellations_dir=tmp_path, ) assert img is not None - assert img_path.name == f"SUIT_{space}Space_1mm.nii" + assert img_path.name == f"SUIT_{space_key}Space_1mm.nii" + assert parcellation_space == space assert len(label) == 34 assert_array_equal(img.header["pixdim"][1:4], [1, 1, 1]) # type: ignore @@ -400,16 +423,17 @@ def test_tian_3T_6thgeneration( assert "TianxS3x3TxMNI6thgeneration" in parcellations assert "TianxS4x3TxMNI6thgeneration" in parcellations # Load parcellation - img, lbl, fname = load_parcellation( + img, lbl, fname, parcellation_space_1 = load_parcellation( name=f"TianxS{scale}x3TxMNI6thgeneration", parcellations_dir=tmp_path ) fname1 = f"Tian_Subcortex_S{scale}_3T_1mm.nii.gz" assert img is not None assert fname.name == fname1 + assert parcellation_space_1 == "MNI152NLin6Asym" assert len(lbl) == n_label assert_array_equal(img.header["pixdim"][1:4], [1, 1, 1]) # type: ignore # Load parcellation - img, lbl, fname = load_parcellation( + img, lbl, fname, parcellation_space_2 = load_parcellation( name=f"TianxS{scale}x3TxMNI6thgeneration", parcellations_dir=tmp_path, resolution=2, @@ -417,6 +441,7 @@ def test_tian_3T_6thgeneration( fname1 = f"Tian_Subcortex_S{scale}_3T.nii.gz" assert img is not None assert fname.name == fname1 + assert parcellation_space_2 == "MNI152NLin6Asym" assert len(lbl) == n_label assert_array_equal(img.header["pixdim"][1:4], [2, 2, 2]) # type: ignore @@ -445,13 +470,14 @@ def test_tian_3T_nonlinear2009cAsym( assert "TianxS3x3TxMNInonlinear2009cAsym" in parcellations assert "TianxS4x3TxMNInonlinear2009cAsym" in parcellations # Load parcellation - img, lbl, fname = load_parcellation( + img, lbl, fname, space = load_parcellation( name=f"TianxS{scale}x3TxMNInonlinear2009cAsym", parcellations_dir=tmp_path, ) fname1 = f"Tian_Subcortex_S{scale}_3T_2009cAsym.nii.gz" assert img is not None assert fname.name == fname1 + assert space == "MNI152NLin2009cAsym" assert len(lbl) == n_label assert_array_equal(img.header["pixdim"][1:4], [2, 2, 2]) # type: ignore @@ -480,12 +506,13 @@ def test_tian_7T_6thgeneration( assert "TianxS3x7TxMNI6thgeneration" in parcellations assert "TianxS4x7TxMNI6thgeneration" in parcellations # Load parcellation - img, lbl, fname = load_parcellation( + img, lbl, fname, space = load_parcellation( name=f"TianxS{scale}x7TxMNI6thgeneration", parcellations_dir=tmp_path ) fname1 = f"Tian_Subcortex_S{scale}_7T.nii.gz" assert img is not None assert fname.name == fname1 + assert space == "MNI152NLin6Asym" assert len(lbl) == n_label assert_array_almost_equal( img.header["pixdim"][1:4], [1.6, 1.6, 1.6] # type: ignore @@ -506,13 +533,13 @@ def test_retrieve_tian_incorrect_space(tmp_path: Path) -> None: parcellations_dir=tmp_path, resolution=1, scale=1, space="wrong" ) - with pytest.raises(ValueError, match=r"MNI6thgeneration"): + with pytest.raises(ValueError, match=r"MNI152NLin6Asym"): _retrieve_tian( parcellations_dir=tmp_path, resolution=1, scale=1, magneticfield="7T", - space="MNInonlinear2009cAsym", + space="MNI152NLin2009cAsym", ) @@ -548,7 +575,7 @@ def test_retrieve_tian_incorrect_scale(tmp_path: Path) -> None: parcellations_dir=tmp_path, resolution=1, scale=5, - space="MNI6thgeneration", + space="MNI152NLin6Asym", ) @@ -567,11 +594,12 @@ def test_aicha(tmp_path: Path, version: int) -> None: parcellations = list_parcellations() assert f"AICHA_v{version}" in parcellations # Load parcellation - img, label, img_path = load_parcellation( + img, label, img_path, space = load_parcellation( name=f"AICHA_v{version}", parcellations_dir=tmp_path ) assert img is not None assert img_path.name == "AICHA.nii" + assert space == "IXI549Space" assert len(label) == 384 assert_array_equal(img.header["pixdim"][1:4], [2, 2, 2]) # type: ignore @@ -635,13 +663,14 @@ def test_shen( parcellations = list_parcellations() assert f"Shen_{year}_{n_rois}" in parcellations # Load parcellation - img, label, img_path = load_parcellation( + img, label, img_path, space = load_parcellation( name=f"Shen_{year}_{n_rois}", parcellations_dir=tmp_path, resolution=resolution, ) assert img is not None assert img_name in img_path.name + assert space == "MNI152NLin2009cAsym" assert len(label) == n_labels assert_array_equal( img.header["pixdim"][1:4], 3 * [resolution] # type: ignore @@ -830,13 +859,14 @@ def test_yan( f"{int(resolution)}mm.nii.gz" ) # Load parcellation - img, label, img_path = load_parcellation( + img, label, img_path, space = load_parcellation( name=parcellation_name, # type: ignore parcellations_dir=tmp_path, resolution=resolution, ) assert img is not None assert img_path.name == parcellation_file # type: ignore + assert space == "MNI152NLin6Asym" assert len(label) == n_rois assert_array_equal( img.header["pixdim"][1:4], 3 * [resolution] # type: ignore @@ -927,10 +957,10 @@ def test_retrieve_yan_incorrect_kong_networks(tmp_path: Path) -> None: def test_merge_parcellations() -> None: """Test merging parcellations.""" # load some parcellations for testing - schaefer_parcellation, schaefer_labels, _ = load_parcellation( + schaefer_parcellation, schaefer_labels, _, _ = load_parcellation( "Schaefer100x17" ) - tian_parcellation, tian_labels, _ = load_parcellation( + tian_parcellation, tian_labels, _, _ = load_parcellation( "TianxS2x3TxMNInonlinear2009cAsym" ) # prepare the list of the actual parcellations @@ -963,7 +993,7 @@ def test_merge_parcellations_3D_multiple_non_overlapping( """ # Get the testing parcellation - parcellation, labels, _ = load_parcellation("Schaefer100x7") + parcellation, labels, _, _ = load_parcellation("Schaefer100x7") assert parcellation is not None @@ -984,7 +1014,7 @@ def test_merge_parcellations_3D_multiple_non_overlapping( names = ["high", "low"] labels_lists = [labels1, labels2] - merged_parc, merged_labels = merge_parcellations( + merged_parc, _ = merge_parcellations( parcellation_list, names, labels_lists ) @@ -994,18 +1024,11 @@ def test_merge_parcellations_3D_multiple_non_overlapping( assert len(np.unique(parc_data)) == 101 # 100 + 1 because background 0 -def test_merge_parcellations_3D_multiple_overlapping(tmp_path: Path) -> None: - """Test merge_parcellations with multiple overlapping parcellations. - - Parameters - ---------- - tmp_path : pathlib.Path - The path to the test directory. - - """ +def test_merge_parcellations_3D_multiple_overlapping() -> None: + """Test merge_parcellations with multiple overlapping parcellations.""" # Get the testing parcellation - parcellation, labels, _ = load_parcellation("Schaefer100x7") + parcellation, labels, _, _ = load_parcellation("Schaefer100x7") assert parcellation is not None @@ -1038,20 +1061,11 @@ def test_merge_parcellations_3D_multiple_overlapping(tmp_path: Path) -> None: assert len(np.unique(parc_data)) == 101 # 100 + 1 because background 0 -def test_merge_parcellations_3D_multiple_duplicated_labels( - tmp_path: Path, -) -> None: - """Test merge_parcellations with two parcellations with duplicated labels. - - Parameters - ---------- - tmp_path : pathlib.Path - The path to the test directory. - - """ +def test_merge_parcellations_3D_multiple_duplicated_labels() -> None: + """Test merge_parcellations with duplicated labels.""" # Get the testing parcellation - parcellation, labels, _ = load_parcellation("Schaefer100x7") + parcellation, labels, _, _ = load_parcellation("Schaefer100x7") assert parcellation is not None @@ -1100,7 +1114,7 @@ def test_get_parcellation_single() -> None: assert tailored_parcellation.shape == vbm_gm_img.shape assert_array_equal(tailored_parcellation.affine, vbm_gm_img.affine) # Get raw parcellation - raw_parcellation, raw_labels, _ = load_parcellation( + raw_parcellation, raw_labels, _, _ = load_parcellation( "Schaefer100x7", resolution=1.5, ) @@ -1118,8 +1132,8 @@ def test_get_parcellation_single() -> None: assert tailored_labels == raw_labels -def test_get_parcellation_multi() -> None: - """Test tailored multi parcellation fetch.""" +def test_get_parcellation_multi_same_space() -> None: + """Test tailored multi parcellation fetch in same space.""" reader = DefaultDataReader() with OasisVBMTestingDataGrabber() as dg: element = dg["sub-01"] @@ -1130,7 +1144,7 @@ def test_get_parcellation_multi() -> None: tailored_parcellation, tailored_labels = get_parcellation( parcellation=[ "Schaefer100x7", - "TianxS2x3TxMNInonlinear2009cAsym", + "TianxS2x3TxMNI6thgeneration", ], target_data=vbm_gm, ) @@ -1142,10 +1156,10 @@ def test_get_parcellation_multi() -> None: raw_labels = [] parcellations_names = [ "Schaefer100x7", - "TianxS2x3TxMNInonlinear2009cAsym", + "TianxS2x3TxMNI6thgeneration", ] for name in parcellations_names: - img, labels, _ = load_parcellation(name=name, resolution=1.5) + img, labels, _, _ = load_parcellation(name=name, resolution=1.5) # Resample raw parcellations resampled_img = resample_to_img( source_img=img, @@ -1167,3 +1181,21 @@ def test_get_parcellation_multi() -> None: merged_resampled_parcellations.get_fdata(), ) assert tailored_labels == merged_labels + + +def test_get_parcellation_multi_different_space() -> None: + """Test tailored multi parcellation fetch in different space.""" + reader = DefaultDataReader() + with OasisVBMTestingDataGrabber() as dg: + element = dg["sub-01"] + element_data = reader.fit_transform(element) + vbm_gm = element_data["VBM_GM"] + # Get tailored parcellation + with pytest.raises(RuntimeError, match="unable to merge."): + get_parcellation( + parcellation=[ + "Schaefer100x7", + "SUITxSUIT", + ], + target_data=vbm_gm, + ) diff --git a/junifer/datagrabber/aomic/id1000.py b/junifer/datagrabber/aomic/id1000.py index 4e0859ff6..a955f6b8a 100644 --- a/junifer/datagrabber/aomic/id1000.py +++ b/junifer/datagrabber/aomic/id1000.py @@ -121,6 +121,10 @@ class DataladAOMICID1000(PatternDataladDataGrabber): out = super().get_item(subject=subject) if out.get("BOLD"): out["BOLD"]["mask_item"] = "BOLD_mask" + # Add space information + out["BOLD"].update({"space": "MNI152NLin2009cAsym"}) if out.get("T1w"): out["T1w"]["mask_item"] = "T1w_mask" + # Add space information + out["T1w"].update({"space": "native"}) return out diff --git a/junifer/datagrabber/aomic/piop1.py b/junifer/datagrabber/aomic/piop1.py index 0a3c1ecf8..88a59005c 100644 --- a/junifer/datagrabber/aomic/piop1.py +++ b/junifer/datagrabber/aomic/piop1.py @@ -166,8 +166,12 @@ class DataladAOMICPIOP1(PatternDataladDataGrabber): out = super().get_item(subject=subject, task=new_task) if out.get("BOLD"): out["BOLD"]["mask_item"] = "BOLD_mask" + # Add space information + out["BOLD"].update({"space": "MNI152NLin2009cAsym"}) if out.get("T1w"): out["T1w"]["mask_item"] = "T1w_mask" + # Add space information + out["T1w"].update({"space": "native"}) return out def get_elements(self) -> List: diff --git a/junifer/datagrabber/aomic/piop2.py b/junifer/datagrabber/aomic/piop2.py index 934cd48fb..a973d93e2 100644 --- a/junifer/datagrabber/aomic/piop2.py +++ b/junifer/datagrabber/aomic/piop2.py @@ -166,6 +166,10 @@ class DataladAOMICPIOP2(PatternDataladDataGrabber): out = super().get_item(subject=subject, task=f"{task}_acq-seq") if out.get("BOLD"): out["BOLD"]["mask_item"] = "BOLD_mask" + # Add space information + out["BOLD"].update({"space": "MNI152NLin2009cAsym"}) if out.get("T1w"): out["T1w"]["mask_item"] = "T1w_mask" + # Add space information + out["T1w"].update({"space": "native"}) return out diff --git a/junifer/datagrabber/base.py b/junifer/datagrabber/base.py index a34db6bf2..76384d48f 100644 --- a/junifer/datagrabber/base.py +++ b/junifer/datagrabber/base.py @@ -74,12 +74,16 @@ class BaseDataGrabber(ABC, UpdateMetaMixin): """ logger.info(f"Getting element {element}") + # Convert element to tuple if not already if not isinstance(element, tuple): element = (element,) - named_element = dict(zip(self.get_element_keys(), element)) + # Zip through element keys and actual values to construct element + # access dictionary + named_element: Dict = dict(zip(self.get_element_keys(), element)) logger.debug(f"Named element: {named_element}") + # Fetch element out = self.get_item(**named_element) - + # Update metadata for _, t_val in out.items(): self.update_meta(t_val, "datagrabber") t_val["meta"]["element"] = named_element diff --git a/junifer/datagrabber/pattern.py b/junifer/datagrabber/pattern.py index eba2ebd66..7c60d476a 100644 --- a/junifer/datagrabber/pattern.py +++ b/junifer/datagrabber/pattern.py @@ -42,6 +42,11 @@ class PatternDataGrabber(BaseDataGrabber): confounds_format : {"fmriprep", "adhoc"} or None, optional The format of the confounds for the dataset (default None). + Raises + ------ + ValueError + If ``confounds_format`` is invalid. + """ def __init__( @@ -55,6 +60,7 @@ class PatternDataGrabber(BaseDataGrabber): # Validate patterns validate_patterns(types=types, patterns=patterns) + # Convert replacements to list if not already if not isinstance(replacements, list): replacements = [replacements] # Validate replacements diff --git a/junifer/datareader/default.py b/junifer/datareader/default.py index 2956e3aec..c2a486e93 100644 --- a/junifer/datareader/default.py +++ b/junifer/datareader/default.py @@ -91,42 +91,65 @@ class DefaultDataReader(PipelineStepMixin, UpdateMetaMixin): The processed output as dictionary. The "data" key is added to each data type dictionary. + Warns + ----- + RuntimeWarning + If input data type has no key called ``"path"``. + """ - # For each type of data, try to read it + # Copy input to not modify the original out = input.copy() + # Set default extra parameters if params is None: params = {} + # For each type of data, try to read it for type_ in input.keys(): + # Check for malformed datagrabber specification if "path" not in input[type_]: warn_with_log( f"Input type {type_} does not provide a path. Skipping." ) continue + + # Retrieve actual path t_path = input[type_]["path"] + # Retrieve loading params for the data type t_params = params.get(type_, {}) - # Convert to Path if datareader is not well done + # Convert str to Path if not isinstance(t_path, Path): t_path = Path(t_path) out[type_]["path"] = t_path - logger.info(f"Reading {type_} from {t_path.as_posix()}") - fread = None + logger.info(f"Reading {type_} from {t_path.as_posix()}") + # Initialize variable for file data + fread = None + # Lowercase path fname = t_path.name.lower() + # Loop through extensions to find the correct one for ext, ftype in _extensions.items(): if fname.endswith(ext): logger.info(f"{type_} is type {ftype}") + # Retrieve reader function reader_func = _readers[ftype]["func"] + # Retrieve reader function params reader_params = _readers[ftype]["params"] + # Update reader function params if reader_params is not None: t_params.update(reader_params) logger.debug(f"Calling {reader_func} with {t_params}") + # Read data fread = reader_func(t_path, **t_params) break + # If no file data is found due to unknown extension if fread is None: logger.info( f"Unknown file type {t_path.as_posix()}, skipping reading" ) + + # Set file data for output out[type_]["data"] = fread + # Update metadata for step self.update_meta(out[type_], "datareader") + return out diff --git a/junifer/markers/base.py b/junifer/markers/base.py index 87553273e..6a652d200 100644 --- a/junifer/markers/base.py +++ b/junifer/markers/base.py @@ -34,12 +34,15 @@ class BaseMarker(ABC, PipelineStepMixin, UpdateMetaMixin): on: Optional[Union[List[str], str]] = None, name: Optional[str] = None, ) -> None: + # Use all data types if not provided if on is None: on = self.get_valid_inputs() + # Convert data types to list if not isinstance(on, list): on = [on] + # Set default name if not provided self.name = self.__class__.__name__ if name is None else name - + # Check if required inputs are found if any(x not in self.get_valid_inputs() for x in on): wrong_on = [x for x in on if x not in self.get_valid_inputs()] raise ValueError(f"{self.name} cannot be computed on {wrong_on}") diff --git a/junifer/markers/collection.py b/junifer/markers/collection.py index bf40c8bb2..01c9194a8 100644 --- a/junifer/markers/collection.py +++ b/junifer/markers/collection.py @@ -9,7 +9,7 @@ from typing import TYPE_CHECKING, Dict, List, Optional from ..datareader.default import DefaultDataReader from ..markers.base import BaseMarker -from ..pipeline import PipelineStepMixin +from ..pipeline import PipelineStepMixin, WorkDirManager from ..preprocess.base import BasePreprocessor from ..storage.base import BaseFeatureStorage from ..utils import logger @@ -98,6 +98,9 @@ class MarkerCollection: out[marker.name] = m_value logger.info("Marker collection fitting done") + # Cleanup element directory + WorkDirManager().cleanup_elementdir() + return None if self._storage else out def validate(self, datagrabber: "BaseDataGrabber") -> None: diff --git a/junifer/markers/functional_connectivity/tests/test_crossparcellation_functional_connectivity.py b/junifer/markers/functional_connectivity/tests/test_crossparcellation_functional_connectivity.py index bf2787ac9..720b817fa 100644 --- a/junifer/markers/functional_connectivity/tests/test_crossparcellation_functional_connectivity.py +++ b/junifer/markers/functional_connectivity/tests/test_crossparcellation_functional_connectivity.py @@ -29,6 +29,7 @@ def test_compute() -> None: "data": niimg, "path": out["BOLD"]["path"], "meta": {"element": "sub001"}, + "space": "MNI", } } diff --git a/junifer/markers/functional_connectivity/tests/test_edge_functional_connectivity_parcels.py b/junifer/markers/functional_connectivity/tests/test_edge_functional_connectivity_parcels.py index 5257980a5..7822652fc 100644 --- a/junifer/markers/functional_connectivity/tests/test_edge_functional_connectivity_parcels.py +++ b/junifer/markers/functional_connectivity/tests/test_edge_functional_connectivity_parcels.py @@ -30,7 +30,9 @@ def test_EdgeCentricFCParcels(tmp_path: Path) -> None: parcellation="TianxS1x3TxMNInonlinear2009cAsym", cor_method_params={"empirical": True}, ) - all_out = efc.fit_transform({"BOLD": {"data": fmri_img, "meta": {}}}) + all_out = efc.fit_transform( + {"BOLD": {"data": fmri_img, "meta": {}, "space": "MNI"}} + ) out = all_out["BOLD"] @@ -51,7 +53,7 @@ def test_EdgeCentricFCParcels(tmp_path: Path) -> None: # Single storage, must be the uri storage = SQLiteFeatureStorage(uri=uri, upsert="ignore") meta = {"element": {"subject": "test"}, "dependencies": {"numpy"}} - input = {"BOLD": {"data": fmri_img, "meta": meta}} + input = {"BOLD": {"data": fmri_img, "meta": meta, "space": "MNI"}} all_out = efc.fit_transform(input, storage=storage) features = storage.list_features() diff --git a/junifer/markers/functional_connectivity/tests/test_edge_functional_connectivity_spheres.py b/junifer/markers/functional_connectivity/tests/test_edge_functional_connectivity_spheres.py index 6b9ad17f6..8695e55f5 100644 --- a/junifer/markers/functional_connectivity/tests/test_edge_functional_connectivity_spheres.py +++ b/junifer/markers/functional_connectivity/tests/test_edge_functional_connectivity_spheres.py @@ -28,7 +28,9 @@ def test_EdgeCentricFCSpheres(tmp_path: Path) -> None: efc = EdgeCentricFCSpheres( coords="DMNBuckner", radius=5.0, cor_method="correlation" ) - all_out = efc.fit_transform({"BOLD": {"data": fmri_img, "meta": {}}}) + all_out = efc.fit_transform( + {"BOLD": {"data": fmri_img, "meta": {}, "space": "MNI"}} + ) out = all_out["BOLD"] @@ -58,13 +60,15 @@ def test_EdgeCentricFCSpheres(tmp_path: Path) -> None: "element": {"subject": "sub001"}, "dependencies": {"nilearn"}, } - all_out = efc.fit_transform({"BOLD": {"data": fmri_img, "meta": meta}}) + all_out = efc.fit_transform( + {"BOLD": {"data": fmri_img, "meta": meta, "space": "MNI"}} + ) uri = tmp_path / "test_fc_parcellation.sqlite" # Single storage, must be the uri storage = SQLiteFeatureStorage(uri=uri, upsert="ignore") meta = {"element": {"subject": "test"}, "dependencies": {"numpy"}} - input = {"BOLD": {"data": fmri_img, "meta": meta}} + input = {"BOLD": {"data": fmri_img, "meta": meta, "space": "MNI"}} all_out = efc.fit_transform(input, storage=storage) features = storage.list_features() diff --git a/junifer/markers/functional_connectivity/tests/test_functional_connectivity_parcels.py b/junifer/markers/functional_connectivity/tests/test_functional_connectivity_parcels.py index aa438d8a5..560635f26 100644 --- a/junifer/markers/functional_connectivity/tests/test_functional_connectivity_parcels.py +++ b/junifer/markers/functional_connectivity/tests/test_functional_connectivity_parcels.py @@ -32,7 +32,9 @@ def test_FunctionalConnectivityParcels(tmp_path: Path) -> None: fmri_img = image.concat_imgs(ni_data.func) # type: ignore fc = FunctionalConnectivityParcels(parcellation="Schaefer100x7") - all_out = fc.fit_transform({"BOLD": {"data": fmri_img, "meta": {}}}) + all_out = fc.fit_transform( + {"BOLD": {"data": fmri_img, "meta": {}, "space": "MNI"}} + ) out = all_out["BOLD"] @@ -52,7 +54,7 @@ def test_FunctionalConnectivityParcels(tmp_path: Path) -> None: "element": {"subject": "sub001"}, "dependencies": {"nilearn"}, } - ts = pa.compute({"data": fmri_img, "meta": meta}) + ts = pa.compute({"data": fmri_img, "meta": meta, "space": "MNI"}) # compare with nilearn # Get the testing parcellation (for nilearn) @@ -80,13 +82,15 @@ def test_FunctionalConnectivityParcels(tmp_path: Path) -> None: parcellation="Schaefer100x7", cor_method_params={"empirical": True} ) - all_out = fc.fit_transform({"BOLD": {"data": fmri_img, "meta": meta}}) + all_out = fc.fit_transform( + {"BOLD": {"data": fmri_img, "meta": meta, "space": "MNI"}} + ) uri = tmp_path / "test_fc_parcellation.sqlite" # Single storage, must be the uri storage = SQLiteFeatureStorage(uri=uri, upsert="ignore") meta = {"element": {"subject": "test"}, "dependencies": {"numpy"}} - input = {"BOLD": {"data": fmri_img, "meta": meta}} + input = {"BOLD": {"data": fmri_img, "meta": meta, "space": "MNI"}} all_out = fc.fit_transform(input, storage=storage) features = storage.list_features() diff --git a/junifer/markers/functional_connectivity/tests/test_functional_connectivity_spheres.py b/junifer/markers/functional_connectivity/tests/test_functional_connectivity_spheres.py index ddc4c842b..6f695c787 100644 --- a/junifer/markers/functional_connectivity/tests/test_functional_connectivity_spheres.py +++ b/junifer/markers/functional_connectivity/tests/test_functional_connectivity_spheres.py @@ -36,7 +36,9 @@ def test_FunctionalConnectivitySpheres(tmp_path: Path) -> None: fc = FunctionalConnectivitySpheres( coords="DMNBuckner", radius=5.0, cor_method="correlation" ) - all_out = fc.fit_transform({"BOLD": {"data": fmri_img, "meta": {}}}) + all_out = fc.fit_transform( + {"BOLD": {"data": fmri_img, "meta": {}, "space": "MNI"}} + ) out = all_out["BOLD"] @@ -52,7 +54,7 @@ def test_FunctionalConnectivitySpheres(tmp_path: Path) -> None: sa = SphereAggregation( coords="DMNBuckner", radius=5.0, method="mean", on="BOLD" ) - ts = sa.compute({"data": fmri_img, "meta": {}}) + ts = sa.compute({"data": fmri_img, "meta": {}, "space": "MNI"}) # Check that FC are almost equal when using nileran cm = ConnectivityMeasure(kind="correlation") @@ -69,7 +71,7 @@ def test_FunctionalConnectivitySpheres(tmp_path: Path) -> None: "element": {"subject": "test"}, "dependencies": {"numpy", "nilearn"}, } - input = {"BOLD": {"data": fmri_img, "meta": meta}} + input = {"BOLD": {"data": fmri_img, "meta": meta, "space": "MNI"}} all_out = fc.fit_transform(input, storage=storage) features = storage.list_features() @@ -99,7 +101,9 @@ def test_FunctionalConnectivitySpheres_empirical(tmp_path: Path) -> None: cor_method="correlation", cor_method_params={"empirical": True}, ) - all_out = fc.fit_transform({"BOLD": {"data": fmri_img, "meta": {}}}) + all_out = fc.fit_transform( + {"BOLD": {"data": fmri_img, "meta": {}, "space": "MNI"}} + ) out = all_out["BOLD"] @@ -115,7 +119,7 @@ def test_FunctionalConnectivitySpheres_empirical(tmp_path: Path) -> None: sa = SphereAggregation( coords="DMNBuckner", radius=5.0, method="mean", on="BOLD" ) - ts = sa.compute({"data": fmri_img}) + ts = sa.compute({"data": fmri_img, "space": "MNI"}) # Check that FC are almost equal when using nileran cm = ConnectivityMeasure( diff --git a/junifer/markers/reho/tests/test_reho_parcels.py b/junifer/markers/reho/tests/test_reho_parcels.py index 488eda53a..6bd7564bd 100644 --- a/junifer/markers/reho/tests/test_reho_parcels.py +++ b/junifer/markers/reho/tests/test_reho_parcels.py @@ -39,7 +39,14 @@ def test_reho_parcels_computation(tmp_path: Path) -> None: reho_parcels_marker = ReHoParcels(parcellation=PARCELLATION) # Fit transform marker on data reho_parcels_output = reho_parcels_marker.fit_transform( - {"BOLD": {"path": "/tmp", "data": fmri_img, "meta": {}}} + { + "BOLD": { + "path": "/tmp", + "data": fmri_img, + "meta": {}, + "space": "MNI", + } + } ) # Get BOLD output reho_parcels_output_bold = reho_parcels_output["BOLD"] @@ -82,7 +89,14 @@ def test_reho_parcels_computation_comparison(tmp_path: Path) -> None: ) # Fit transform marker on data reho_parcels_output_python = reho_parcels_marker_python.fit_transform( - {"BOLD": {"path": "/tmp", "data": fmri_img, "meta": {}}} + { + "BOLD": { + "path": "/tmp", + "data": fmri_img, + "meta": {}, + "space": "MNI", + } + } ) # Get BOLD output reho_parcels_output_bold_python = reho_parcels_output_python["BOLD"] @@ -93,7 +107,14 @@ def test_reho_parcels_computation_comparison(tmp_path: Path) -> None: ) # Fit transform marker on data reho_parcels_output_afni = reho_parcels_marker_afni.fit_transform( - {"BOLD": {"path": "/tmp", "data": fmri_img, "meta": {}}} + { + "BOLD": { + "path": "/tmp", + "data": fmri_img, + "meta": {}, + "space": "MNI", + } + } ) # Get BOLD output reho_parcels_output_bold_afni = reho_parcels_output_afni["BOLD"] @@ -134,6 +155,13 @@ def test_reho_parcels_storage(tmp_path: Path) -> None: } # only requires element key for storing # Fit transform marker on data with storage reho_parcels_marker.fit_transform( - input={"BOLD": {"path": "/tmp", "data": fmri_img, "meta": meta}}, + input={ + "BOLD": { + "path": "/tmp", + "data": fmri_img, + "meta": meta, + "space": "MNI", + } + }, storage=reho_parcels_storage, ) diff --git a/junifer/markers/reho/tests/test_reho_spheres.py b/junifer/markers/reho/tests/test_reho_spheres.py index c662980be..e3eba63e8 100644 --- a/junifer/markers/reho/tests/test_reho_spheres.py +++ b/junifer/markers/reho/tests/test_reho_spheres.py @@ -39,7 +39,14 @@ def test_reho_spheres_computation(tmp_path: Path) -> None: reho_spheres_marker = ReHoSpheres(coords=COORDINATES, radius=10.0) # Fit transform marker on data reho_spheres_output = reho_spheres_marker.fit_transform( - {"BOLD": {"path": "/tmp", "data": fmri_img, "meta": {}}} + { + "BOLD": { + "path": "/tmp", + "data": fmri_img, + "meta": {}, + "space": "MNI", + } + } ) # Get BOLD output reho_spheres_output_bold = reho_spheres_output["BOLD"] @@ -82,7 +89,14 @@ def test_reho_spheres_computation_comparison(tmp_path: Path) -> None: ) # Fit transform marker on data reho_spheres_output_python = reho_spheres_marker_python.fit_transform( - {"BOLD": {"path": "/tmp", "data": fmri_img, "meta": {}}} + { + "BOLD": { + "path": "/tmp", + "data": fmri_img, + "meta": {}, + "space": "MNI", + } + } ) # Get BOLD output reho_spheres_output_bold_python = reho_spheres_output_python["BOLD"] @@ -93,7 +107,14 @@ def test_reho_spheres_computation_comparison(tmp_path: Path) -> None: ) # Fit transform marker on data reho_spheres_output_afni = reho_spheres_marker_afni.fit_transform( - {"BOLD": {"path": "/tmp", "data": fmri_img, "meta": {}}} + { + "BOLD": { + "path": "/tmp", + "data": fmri_img, + "meta": {}, + "space": "MNI", + } + } ) # Get BOLD output reho_spheres_output_bold_afni = reho_spheres_output_afni["BOLD"] @@ -134,6 +155,13 @@ def test_reho_spheres_storage(tmp_path: Path) -> None: } # only requires element key for storing # Fit transform marker on data with storage reho_spheres_marker.fit_transform( - input={"BOLD": {"path": "/tmp", "data": fmri_img, "meta": meta}}, + input={ + "BOLD": { + "path": "/tmp", + "data": fmri_img, + "meta": meta, + "space": "MNI", + } + }, storage=reho_spheres_storage, ) diff --git a/junifer/markers/temporal_snr/tests/test_temporal_snr_parcels.py b/junifer/markers/temporal_snr/tests/test_temporal_snr_parcels.py index e0ef668f4..e888a74ff 100644 --- a/junifer/markers/temporal_snr/tests/test_temporal_snr_parcels.py +++ b/junifer/markers/temporal_snr/tests/test_temporal_snr_parcels.py @@ -26,7 +26,7 @@ def test_TemporalSNRParcels(tmp_path: Path) -> None: tsnr_parcels = TemporalSNRParcels(parcellation="Schaefer100x7") all_out = tsnr_parcels.fit_transform( - {"BOLD": {"data": fmri_img, "meta": {}}} + {"BOLD": {"data": fmri_img, "meta": {}, "space": "MNI"}} ) out = all_out["BOLD"] @@ -45,7 +45,7 @@ def test_TemporalSNRParcels(tmp_path: Path) -> None: # Single storage, must be the uri storage = SQLiteFeatureStorage(uri=uri, upsert="ignore") meta = {"element": {"subject": "test"}, "dependencies": {"numpy"}} - input = {"BOLD": {"data": fmri_img, "meta": meta}} + input = {"BOLD": {"data": fmri_img, "meta": meta, "space": "MNI"}} all_out = tsnr_parcels.fit_transform(input, storage=storage) features = storage.list_features() diff --git a/junifer/markers/temporal_snr/tests/test_temporal_snr_spheres.py b/junifer/markers/temporal_snr/tests/test_temporal_snr_spheres.py index 448de5e54..70ad1752e 100644 --- a/junifer/markers/temporal_snr/tests/test_temporal_snr_spheres.py +++ b/junifer/markers/temporal_snr/tests/test_temporal_snr_spheres.py @@ -27,7 +27,7 @@ def test_TemporalSNRSpheres(tmp_path: Path) -> None: tsnr_spheres = TemporalSNRSpheres(coords="DMNBuckner", radius=5.0) all_out = tsnr_spheres.fit_transform( - {"BOLD": {"data": fmri_img, "meta": {}}} + {"BOLD": {"data": fmri_img, "meta": {}, "space": "MNI"}} ) out = all_out["BOLD"] @@ -48,7 +48,7 @@ def test_TemporalSNRSpheres(tmp_path: Path) -> None: "element": {"subject": "test"}, "dependencies": {"numpy", "nilearn"}, } - input = {"BOLD": {"data": fmri_img, "meta": meta}} + input = {"BOLD": {"data": fmri_img, "meta": meta, "space": "MNI"}} all_out = tsnr_spheres.fit_transform(input, storage=storage) features = storage.list_features() diff --git a/junifer/markers/tests/test_ets_rss.py b/junifer/markers/tests/test_ets_rss.py index 0c6c872c4..04c102893 100644 --- a/junifer/markers/tests/test_ets_rss.py +++ b/junifer/markers/tests/test_ets_rss.py @@ -29,13 +29,17 @@ def test_compute() -> None: # Load BOLD image niimg = image.load_img(str(out["BOLD"]["path"].absolute())) # Create input data - input_dict = {"data": niimg, "path": out["BOLD"]["path"]} + input_dict = { + "data": niimg, + "path": out["BOLD"]["path"], + "space": "MNI", + } # Compute the RSSETSMarker ets_rss_marker = RSSETSMarker(parcellation=PARCELLATION) new_out = ets_rss_marker.compute(input_dict) # Load parcellation - test_parcellation, _, _ = load_parcellation(PARCELLATION) + test_parcellation, _, _, _ = load_parcellation(PARCELLATION) # Compute the NiftiLabelsMasker test_masker = NiftiLabelsMasker(test_parcellation) test_ts = test_masker.fit_transform(niimg) diff --git a/junifer/markers/tests/test_parcel_aggregation.py b/junifer/markers/tests/test_parcel_aggregation.py index 84d5d4c58..12bd01c3f 100644 --- a/junifer/markers/tests/test_parcel_aggregation.py +++ b/junifer/markers/tests/test_parcel_aggregation.py @@ -82,7 +82,7 @@ def test_ParcelAggregation_3D() -> None: name="gmd_schaefer100x7_mean", on="VBM_GM", ) # Test passing "on" as a keyword argument - input = {"VBM_GM": {"data": img, "meta": {}}} + input = {"VBM_GM": {"data": img, "meta": {}, "space": "MNI"}} jun_values3d_mean = marker.fit_transform(input)["VBM_GM"]["data"] assert jun_values3d_mean.ndim == 2 @@ -98,7 +98,7 @@ def test_ParcelAggregation_3D() -> None: # Use the ParcelAggregation object marker = ParcelAggregation(parcellation="Schaefer100x7", method="std") - input = {"VBM_GM": {"data": img, "meta": {}}} + input = {"VBM_GM": {"data": img, "meta": {}, "space": "MNI"}} jun_values3d_std = marker.fit_transform(input)["VBM_GM"]["data"] assert jun_values3d_std.ndim == 2 @@ -122,7 +122,7 @@ def test_ParcelAggregation_3D() -> None: method="trim_mean", method_params={"proportiontocut": 0.1}, ) - input = {"VBM_GM": {"data": img, "meta": {}}} + input = {"VBM_GM": {"data": img, "meta": {}, "space": "MNI"}} jun_values3d_tm = marker.fit_transform(input)["VBM_GM"]["data"] assert jun_values3d_tm.ndim == 2 @@ -147,7 +147,7 @@ def test_ParcelAggregation_4D(): # Create ParcelAggregation object marker = ParcelAggregation(parcellation="Schaefer100x7", method="mean") - input = {"BOLD": {"data": fmri_img, "meta": {}}} + input = {"BOLD": {"data": fmri_img, "meta": {}, "space": "MNI"}} jun_values4d = marker.fit_transform(input)["BOLD"]["data"] assert jun_values4d.ndim == 2 @@ -175,7 +175,7 @@ def test_ParcelAggregation_storage(tmp_path: Path) -> None: "element": {"subject": "sub-01", "session": "ses-01"}, "dependencies": {"nilearn", "nibabel"}, } - input = {"VBM_GM": {"data": img, "meta": meta}} + input = {"VBM_GM": {"data": img, "meta": meta, "space": "MNI"}} marker = ParcelAggregation( parcellation="Schaefer100x7", method="mean", on="VBM_GM" ) @@ -194,7 +194,7 @@ def test_ParcelAggregation_storage(tmp_path: Path) -> None: # Get the SPM auditory data subject_data = datasets.fetch_spm_auditory() fmri_img = concat_imgs(subject_data.func) # type: ignore - input = {"BOLD": {"data": fmri_img, "meta": meta}} + input = {"BOLD": {"data": fmri_img, "meta": meta, "space": "MNI"}} marker = ParcelAggregation( parcellation="Schaefer100x7", method="mean", on="BOLD" ) @@ -213,7 +213,7 @@ def test_ParcelAggregation_3D_mask() -> None: parcellation = datasets.fetch_atlas_schaefer_2018(n_rois=100) # Get one mask - mask_img, _ = load_mask("GM_prob0.2") + mask_img, _, _ = load_mask("GM_prob0.2") # Get the oasis VBM data oasis_dataset = datasets.fetch_oasis_vbm(n_subjects=1) @@ -234,7 +234,7 @@ def test_ParcelAggregation_3D_mask() -> None: name="gmd_schaefer100x7_mean", on="VBM_GM", ) # Test passing "on" as a keyword argument - input = {"VBM_GM": {"data": img, "meta": {}}} + input = {"VBM_GM": {"data": img, "meta": {}, "space": "MNI"}} jun_values3d_mean = marker.fit_transform(input)["VBM_GM"]["data"] assert jun_values3d_mean.ndim == 2 @@ -279,7 +279,7 @@ def test_ParcelAggregation_3D_mask_computed() -> None: name="gmd_schaefer100x7_mean", on="VBM_GM", ) # Test passing "on" as a keyword argument - input = {"VBM_GM": {"data": img, "meta": {}}} + input = {"VBM_GM": {"data": img, "meta": {}, "space": "MNI"}} jun_values3d_mean = marker.fit_transform(input)["VBM_GM"]["data"] assert jun_values3d_mean.ndim == 2 @@ -301,7 +301,7 @@ def test_ParcelAggregation_3D_multiple_non_overlapping(tmp_path: Path) -> None: """ # Get the testing parcellation - parcellation, labels, _ = load_parcellation("Schaefer100x7") + parcellation, labels, _, _ = load_parcellation("Schaefer100x7") assert parcellation is not None @@ -330,10 +330,18 @@ def test_ParcelAggregation_3D_multiple_non_overlapping(tmp_path: Path) -> None: nib.save(parcellation2_img, parcellation2_path) register_parcellation( - "Schaefer100x7_low", parcellation1_path, labels1, overwrite=True + name="Schaefer100x7_low", + parcellation_path=parcellation1_path, + parcels_labels=labels1, + space="MNI", + overwrite=True, ) register_parcellation( - "Schaefer100x7_high", parcellation2_path, labels2, overwrite=True + name="Schaefer100x7_high", + parcellation_path=parcellation2_path, + parcels_labels=labels2, + space="MNI", + overwrite=True, ) # Use the ParcelAggregation object on the original parcellation @@ -343,7 +351,7 @@ def test_ParcelAggregation_3D_multiple_non_overlapping(tmp_path: Path) -> None: name="gmd_schaefer100x7_mean", on="VBM_GM", ) # Test passing "on" as a keyword argument - input = {"VBM_GM": {"data": img, "meta": {}}} + input = {"VBM_GM": {"data": img, "meta": {}, "space": "MNI"}} orig_mean = marker_original.fit_transform(input)["VBM_GM"] orig_mean_data = orig_mean["data"] @@ -359,7 +367,7 @@ def test_ParcelAggregation_3D_multiple_non_overlapping(tmp_path: Path) -> None: name="gmd_schaefer100x7_mean", on="VBM_GM", ) # Test passing "on" as a keyword argument - input = {"VBM_GM": {"data": img, "meta": {}}} + input = {"VBM_GM": {"data": img, "meta": {}, "space": "MNI"}} # No warnings should be raised with warnings.catch_warnings(): @@ -387,7 +395,7 @@ def test_ParcelAggregation_3D_multiple_overlapping(tmp_path: Path) -> None: """ # Get the testing parcellation - parcellation, labels, _ = load_parcellation("Schaefer100x7") + parcellation, labels, _, _ = load_parcellation("Schaefer100x7") assert parcellation is not None @@ -418,10 +426,18 @@ def test_ParcelAggregation_3D_multiple_overlapping(tmp_path: Path) -> None: nib.save(parcellation2_img, parcellation2_path) register_parcellation( - "Schaefer100x7_low2", parcellation1_path, labels1, overwrite=True + name="Schaefer100x7_low2", + parcellation_path=parcellation1_path, + parcels_labels=labels1, + space="MNI", + overwrite=True, ) register_parcellation( - "Schaefer100x7_high2", parcellation2_path, labels2, overwrite=True + name="Schaefer100x7_high2", + parcellation_path=parcellation2_path, + parcels_labels=labels2, + space="MNI", + overwrite=True, ) # Use the ParcelAggregation object on the original parcellation @@ -431,7 +447,7 @@ def test_ParcelAggregation_3D_multiple_overlapping(tmp_path: Path) -> None: name="gmd_schaefer100x7_mean", on="VBM_GM", ) # Test passing "on" as a keyword argument - input = {"VBM_GM": {"data": img, "meta": {}}} + input = {"VBM_GM": {"data": img, "meta": {}, "space": "MNI"}} orig_mean = marker_original.fit_transform(input)["VBM_GM"] orig_mean_data = orig_mean["data"] @@ -447,7 +463,7 @@ def test_ParcelAggregation_3D_multiple_overlapping(tmp_path: Path) -> None: name="gmd_schaefer100x7_mean", on="VBM_GM", ) # Test passing "on" as a keyword argument - input = {"VBM_GM": {"data": img, "meta": {}}} + input = {"VBM_GM": {"data": img, "meta": {}, "space": "MNI"}} with pytest.warns(RuntimeWarning, match="overlapping voxels"): split_mean = marker_split.fit_transform(input)["VBM_GM"] split_mean_data = split_mean["data"] @@ -481,7 +497,7 @@ def test_ParcelAggregation_3D_multiple_duplicated_labels( """ # Get the testing parcellation - parcellation, labels, _ = load_parcellation("Schaefer100x7") + parcellation, labels, _, _ = load_parcellation("Schaefer100x7") assert parcellation is not None @@ -510,10 +526,18 @@ def test_ParcelAggregation_3D_multiple_duplicated_labels( nib.save(parcellation2_img, parcellation2_path) register_parcellation( - "Schaefer100x7_low", parcellation1_path, labels1, overwrite=True + name="Schaefer100x7_low", + parcellation_path=parcellation1_path, + parcels_labels=labels1, + space="MNI", + overwrite=True, ) register_parcellation( - "Schaefer100x7_high", parcellation2_path, labels2, overwrite=True + name="Schaefer100x7_high", + parcellation_path=parcellation2_path, + parcels_labels=labels2, + space="MNI", + overwrite=True, ) # Use the ParcelAggregation object on the original parcellation @@ -523,7 +547,7 @@ def test_ParcelAggregation_3D_multiple_duplicated_labels( name="gmd_schaefer100x7_mean", on="VBM_GM", ) # Test passing "on" as a keyword argument - input = {"VBM_GM": {"data": img, "meta": {}}} + input = {"VBM_GM": {"data": img, "meta": {}, "space": "MNI"}} orig_mean = marker_original.fit_transform(input)["VBM_GM"] orig_mean_data = orig_mean["data"] @@ -539,7 +563,7 @@ def test_ParcelAggregation_3D_multiple_duplicated_labels( name="gmd_schaefer100x7_mean", on="VBM_GM", ) # Test passing "on" as a keyword argument - input = {"VBM_GM": {"data": img, "meta": {}}} + input = {"VBM_GM": {"data": img, "meta": {}, "space": "MNI"}} with pytest.warns(RuntimeWarning, match="duplicated labels."): split_mean = marker_split.fit_transform(input)["VBM_GM"] @@ -578,7 +602,7 @@ def test_ParcelAggregation_4D_agg_time(): marker = ParcelAggregation( parcellation="Schaefer100x7", method="mean", time_method="mean" ) - input = {"BOLD": {"data": fmri_img, "meta": {}}} + input = {"BOLD": {"data": fmri_img, "meta": {}, "space": "MNI"}} jun_values4d = marker.fit_transform(input)["BOLD"]["data"] assert jun_values4d.ndim == 1 @@ -593,7 +617,7 @@ def test_ParcelAggregation_4D_agg_time(): time_method_params={"pick": [0]}, ) - input = {"BOLD": {"data": fmri_img, "meta": {}}} + input = {"BOLD": {"data": fmri_img, "meta": {}, "space": "MNI"}} jun_values4d = marker.fit_transform(input)["BOLD"]["data"] assert jun_values4d.ndim == 2 @@ -620,5 +644,11 @@ def test_ParcelAggregation_4D_agg_time(): ) with pytest.warns(RuntimeWarning, match="No time dimension to aggregate"): - input = {"BOLD": {"data": fmri_img.slicer[..., 0:1], "meta": {}}} + input = { + "BOLD": { + "data": fmri_img.slicer[..., 0:1], + "meta": {}, + "space": "MNI", + } + } marker.fit_transform(input) diff --git a/junifer/markers/tests/test_sphere_aggregation.py b/junifer/markers/tests/test_sphere_aggregation.py index 68ea96714..04b9de474 100644 --- a/junifer/markers/tests/test_sphere_aggregation.py +++ b/junifer/markers/tests/test_sphere_aggregation.py @@ -37,7 +37,7 @@ def test_SphereAggregation_input_output() -> None: def test_SphereAggregation_3D() -> None: """Test SphereAggregation object on 3D images.""" # Get the testing coordinates (for nilearn) - coordinates, labels = load_coordinates(COORDS) + coordinates, _, _ = load_coordinates(COORDS) # Get the oasis VBM data oasis_dataset = datasets.fetch_oasis_vbm(n_subjects=1) @@ -52,7 +52,7 @@ def test_SphereAggregation_3D() -> None: marker = SphereAggregation( coords=COORDS, method="mean", radius=RADIUS, on="VBM_GM" ) - input = {"VBM_GM": {"data": img, "meta": {}}} + input = {"VBM_GM": {"data": img, "meta": {}, "space": "MNI"}} jun_values4d = marker.fit_transform(input)["VBM_GM"]["data"] assert jun_values4d.ndim == 2 @@ -63,7 +63,7 @@ def test_SphereAggregation_3D() -> None: def test_SphereAggregation_4D() -> None: """Test SphereAggregation object on 4D images.""" # Get the testing coordinates (for nilearn) - coordinates, _ = load_coordinates(COORDS) + coordinates, _, _ = load_coordinates(COORDS) # Get the SPM auditory data subject_data = datasets.fetch_spm_auditory() @@ -75,7 +75,7 @@ def test_SphereAggregation_4D() -> None: # Create SphereAggregation object marker = SphereAggregation(coords=COORDS, method="mean", radius=RADIUS) - input = {"BOLD": {"data": fmri_img, "meta": {}}} + input = {"BOLD": {"data": fmri_img, "meta": {}, "space": "MNI"}} jun_values4d = marker.fit_transform(input)["BOLD"]["data"] assert jun_values4d.ndim == 2 @@ -103,7 +103,7 @@ def test_SphereAggregation_storage(tmp_path: Path) -> None: "element": {"subject": "sub-01", "session": "ses-01"}, "dependencies": {"nilearn", "nibabel"}, } - input = {"VBM_GM": {"data": img, "meta": meta}} + input = {"VBM_GM": {"data": img, "meta": meta, "space": "MNI"}} marker = SphereAggregation( coords=COORDS, method="mean", radius=RADIUS, on="VBM_GM" ) @@ -122,7 +122,7 @@ def test_SphereAggregation_storage(tmp_path: Path) -> None: # Get the SPM auditory data subject_data = datasets.fetch_spm_auditory() fmri_img = concat_imgs(subject_data.func) # type: ignore - input = {"BOLD": {"data": fmri_img, "meta": meta}} + input = {"BOLD": {"data": fmri_img, "meta": meta, "space": "MNI"}} marker = SphereAggregation( coords=COORDS, method="mean", radius=RADIUS, on="BOLD" ) @@ -137,10 +137,10 @@ def test_SphereAggregation_storage(tmp_path: Path) -> None: def test_SphereAggregation_3D_mask() -> None: """Test SphereAggregation object on 3D images using mask.""" # Get the testing coordinates (for nilearn) - coordinates, _ = load_coordinates(COORDS) + coordinates, _, _ = load_coordinates(COORDS) # Get one mask - mask_img, _ = load_mask("GM_prob0.2") + mask_img, _, _ = load_mask("GM_prob0.2") # Get the oasis VBM data oasis_dataset = datasets.fetch_oasis_vbm(n_subjects=1) @@ -161,7 +161,7 @@ def test_SphereAggregation_3D_mask() -> None: on="VBM_GM", masks="GM_prob0.2", ) - input = {"VBM_GM": {"data": img, "meta": {}}} + input = {"VBM_GM": {"data": img, "meta": {}, "space": "MNI"}} jun_values4d = marker.fit_transform(input)["VBM_GM"]["data"] assert jun_values4d.ndim == 2 @@ -172,7 +172,7 @@ def test_SphereAggregation_3D_mask() -> None: def test_SphereAggregation_4D_agg_time() -> None: """Test SphereAggregation object on 4D images, aggregating time.""" # Get the testing coordinates (for nilearn) - coordinates, _ = load_coordinates(COORDS) + coordinates, _, _ = load_coordinates(COORDS) # Get the SPM auditory data subject_data = datasets.fetch_spm_auditory() @@ -187,7 +187,7 @@ def test_SphereAggregation_4D_agg_time() -> None: marker = SphereAggregation( coords=COORDS, method="mean", radius=RADIUS, time_method="mean" ) - input = {"BOLD": {"data": fmri_img, "meta": {}}} + input = {"BOLD": {"data": fmri_img, "meta": {}, "space": "MNI"}} jun_values4d = marker.fit_transform(input)["BOLD"]["data"] assert jun_values4d.ndim == 1 @@ -203,7 +203,7 @@ def test_SphereAggregation_4D_agg_time() -> None: time_method_params={"pick": [0]}, ) - input = {"BOLD": {"data": fmri_img, "meta": {}}} + input = {"BOLD": {"data": fmri_img, "meta": {}, "space": "MNI"}} jun_values4d = marker.fit_transform(input)["BOLD"]["data"] assert jun_values4d.ndim == 2 @@ -232,5 +232,11 @@ def test_SphereAggregation_4D_agg_time() -> None: ) with pytest.warns(RuntimeWarning, match="No time dimension to aggregate"): - input = {"BOLD": {"data": fmri_img.slicer[..., 0:1], "meta": {}}} + input = { + "BOLD": { + "data": fmri_img.slicer[..., 0:1], + "meta": {}, + "space": "MNI", + } + } marker.fit_transform(input) diff --git a/junifer/testing/datagrabbers.py b/junifer/testing/datagrabbers.py index 280ac3530..ffb744cfd 100644 --- a/junifer/testing/datagrabbers.py +++ b/junifer/testing/datagrabbers.py @@ -55,7 +55,10 @@ class OasisVBMTestingDataGrabber(BaseDataGrabber): """ out = {} i_sub = int(subject.split("-")[1]) - 1 - out["VBM_GM"] = {"path": Path(self._dataset.gray_matter_maps[i_sub])} + out["VBM_GM"] = { + "path": Path(self._dataset.gray_matter_maps[i_sub]), + "space": "MNI152Lin", + } return out @@ -141,8 +144,8 @@ class SPMAuditoryTestingDataGrabber(BaseDataGrabber): anat_fname = self.datadir / f"{subject}_T1w.nii.gz" nib.save(fmri_img, fmri_fname) nib.save(anat_img, anat_fname) - out["BOLD"] = {"path": fmri_fname} - out["T1w"] = {"path": anat_fname} + out["BOLD"] = {"path": fmri_fname, "space": "MNI152Lin"} + out["T1w"] = {"path": anat_fname, "space": "native"} return out @@ -236,12 +239,13 @@ class PartlyCloudyTestingDataGrabber(BaseDataGrabber): """ out = {} i_sub = int(subject.split("-")[1]) - 1 - out["BOLD"] = {"path": Path(self._dataset["func"][i_sub])} - conf_format = "fmriprep" - + out["BOLD"] = { + "path": Path(self._dataset["func"][i_sub]), + "space": "MNI152NLin2009cAsym", + } out["BOLD_confounds"] = { "path": Path(self._dataset["confounds"][i_sub]), - "format": conf_format, + "format": "fmriprep", } return out