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.

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`_

docs/changes/newsfragments/268.enh

@@ -0,0 +1 @@
+Add ``space`` information to existing datagrabbers, masks, parcellations and coordinates by `Synchon Mandal`_ and `Fede Raimondo`_

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

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 <codeless>`.

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

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")
 
-
 ###############################################################################
-# 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
+# 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)
 
 ###############################################################################
-# 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)

junifer/data/coordinates.py

@@ -4,6 +4,7 @@
 #          Synchon Mandal <s.mandal@fz-juelich.de>
 # 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"]

junifer/data/masks.py

@@ -3,6 +3,8 @@
 # Authors: Federico Raimondo <f.raimondo@fz-juelich.de>
 # 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,7 +141,7 @@ 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(
@@ -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:
+        # 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

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,18 +258,90 @@ 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
+        # 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}")

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

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"],
+            )

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,
+            )

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

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:

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

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

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

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

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}")

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:

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",
             }
         }
 

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()

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()

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()

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(

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

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

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()

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()

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)

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)

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)

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