.github/workflows/ci.yml

@@ -22,8 +22,8 @@ jobs:
     - name: Set up system
       run: |
         bash -c "$(curl -fsSL http://neuro.debian.net/_files/neurodebian-travis.sh)"
-        sudo apt-get update -qq
+        sudo apt-get -qq update
-        sudo apt-get install git-annex-standalone
+        sudo apt-get -qq install git-annex-standalone
     - name: Configure git for datalad
       run: |
         git config --global user.email "runner@github.com"
@@ -39,8 +39,29 @@ jobs:
         python -m pip install tox tox-gh-actions
     - name: Install AFNI
       run: |
-        docker pull afni/afni_make_build
+        echo "++ Distro information"
-        echo "$(pwd)/junifer/api/res/afni" >> $GITHUB_PATH
+        lsb_release -a
+
+        echo "++ Add universe repo"
+        sudo add-apt-repository -y universe
+
+        echo "++ Add PPAs for R programs"
+        sudo add-apt-repository -y "ppa:marutter/rrutter4.0"
+        sudo add-apt-repository -y "ppa:c2d4u.team/c2d4u4.0+"
+
+        echo "++ Update package manager info"
+        sudo apt-get -qq update
+
+        echo "++ Get main dependencies"
+        sudo apt-get -qq install -y tcsh libssl-dev gsl-bin libjpeg62 vim curl \
+        build-essential libcurl4-openssl-dev libxml2-dev libgfortran-11-dev \
+        libgomp1 r-base cmake rsync libxm4
+
+        sudo ln -s /usr/lib/x86_64-linux-gnu/libgsl.so.27 /usr/lib/x86_64-linux-gnu/libgsl.so.19
+
+        curl -O https://afni.nimh.nih.gov/pub/dist/bin/misc/@update.afni.binaries
+        sudo tcsh @update.afni.binaries -package linux_ubuntu_16_64 -bindir /afni
+        echo "/afni" >> $GITHUB_PATH
       if: matrix.python-version == 3.10
     - name: Check AFNI
       run: |

docs/builtin.rst

@@ -166,6 +166,14 @@ Available
      - Compute root sum of squares of edgewise timeseries
      - Done
      - 0.0.1
+   * - :class:`junifer.markers.ReHoParcels`
+     - Calculate regional homogeneity over parcellation
+     - Done
+     - 0.0.1
+   * - :class:`junifer.markers.ReHoSpheres`
+     - Calculate regional homogeneity over spheres placed on coordinates
+     - Done
+     - 0.0.1
 
 
 Planned
@@ -184,9 +192,6 @@ Planned
    * - ALFF and (f)ALFF
      - Detect amplitude of low-frequency fluctuation (ALFF) for resting-state fMRI
      - :gh:`35`
-   * - ReHo
-     - Calculate regional homogeneity
-     - :gh:`36`
    * - Permutation entropy, Range entropy, Multiscale entropy and Hurst exponent
      - Calculate Permutation entropy, Range entropy, Multiscale entropy and Hurst exponent
      - :gh:`61`

docs/changes/latest.inc

@@ -96,6 +96,8 @@ Enhancements
 - Refactor :class:`junifer.pipeline.PipelineStepMixin` to improve its implementation and validation for pipeline steps
   (:gh:`152` by `Synchon Mandal`_).
 
+- Implement :class:`junifer.markers.ReHoParcels` and :class:`junifer.markers.ReHoSpheres` markers (:gh:`36` by `Synchon Mandal`_).
+
 Bugs
 ~~~~
 

junifer/api/cli.py

@@ -308,6 +308,7 @@ def setup() -> None:  # pragma: no cover
 def afni_docker() -> None:  # pragma: no cover
     """Configure AFNI-Docker wrappers."""
     import junifer
+
     pkg_path = Path(junifer.__path__[0])  # type: ignore
     afni_wrappers_path = pkg_path / "api" / "res" / "afni"
     msg = f"""

junifer/api/res/afni/3dAFNItoNIFTI

@@ -0,0 +1,3 @@
+#!/bin/bash
+
+run_afni_docker.sh 3dAFNItoNIFTI "$@"

junifer/markers/__init__.py

@@ -12,3 +12,4 @@ from .functional_connectivity_parcels import FunctionalConnectivityParcels
 from .functional_connectivity_spheres import FunctionalConnectivitySpheres
 from .parcel_aggregation import ParcelAggregation
 from .sphere_aggregation import SphereAggregation
+from .reho import ReHoParcels, ReHoSpheres

junifer/markers/reho/__init__.py

@@ -0,0 +1,7 @@
+"""Provide imports for reho sub-package."""
+
+# Authors: Synchon Mandal <s.mandal@fz-juelich.de>
+# License: AGPL
+
+from .reho_parcels import ReHoParcels
+from .reho_spheres import ReHoSpheres

junifer/markers/reho/reho_base.py

@@ -0,0 +1,126 @@
+"""Provide base class for regional homogeneity (ReHo)."""
+
+# Authors: Synchon Mandal <s.mandal@fz-juelich.de>
+# License: AGPL
+
+
+from typing import TYPE_CHECKING, Any, Dict, List, Optional
+
+from ...utils import logger, raise_error
+from ..base import BaseMarker
+from .reho_estimator import ReHoEstimator
+
+
+if TYPE_CHECKING:
+    from nibabel import Nifti1Image
+
+
+class ReHoBase(BaseMarker):
+    """Base class for regional homogeneity computation.
+
+    Parameters
+    ----------
+    use_afni : bool, optional
+        Whether to use AFNI for computing. If None, will use AFNI only
+        if available (default None).
+    name : str, optional
+        The name of the marker. If None, it will use the class name
+        (default None).
+
+    """
+
+    _EXT_DEPENDENCIES = [
+        {
+            "name": "afni",
+            "optional": True,
+            "commands": ["3dReHo", "3dAFNItoNIFTI"],
+        },
+    ]
+
+    def __init__(
+        self,
+        use_afni: Optional[bool] = None,
+        name: Optional[str] = None,
+    ) -> None:
+        super().__init__(on="BOLD", name=name)
+        self.use_afni = use_afni
+
+    def get_valid_inputs(self) -> List[str]:
+        """Get valid data types for input.
+
+        Returns
+        -------
+        list of str
+            The list of data types that can be used as input for this marker.
+
+        """
+        return ["BOLD"]
+
+    def get_output_type(self, input_type: str) -> str:
+        """Get output type.
+
+        Parameters
+        ----------
+        input_type : str
+            The data type input to the marker.
+
+        Returns
+        -------
+        str
+            The storage type output by the marker.
+
+        """
+        return "table"
+
+    def compute_reho_map(
+        self,
+        input: Dict[str, Any],
+        **reho_params: Any,
+    ) -> "Nifti1Image":
+        """Compute.
+
+        Calculates Kendall's W per voxel using neighborhood voxels.
+        Instead of the time series values themselves, Kendall's W uses the
+        relative rank ordering of a hood over all time points to evaluate
+        a parameter W in range 0-1, with 0 reflecting no trend of agreement
+        between time series and 1 reflecting perfect agreement. For more
+        information about the method, please check [1]_.
+
+        Parameters
+        ----------
+        input : dict
+            The BOLD data as dictionary.
+        **reho_params : dict
+            Extra keyword arguments for ReHo.
+
+        Returns
+        -------
+        Niimg-like object
+
+        References
+        ----------
+        .. [1] Jiang, L., & Zuo, X. N. (2016).
+               Regional Homogeneity: A Multimodal, Multiscale Neuroimaging
+               Marker of the Human Connectome.
+               The Neuroscientist, Volume 22(5), Pages 486–505.
+               https://doi.org/10.1177/1073858415595004
+
+        """
+        if self.use_afni is None:
+            raise_error(
+                "Parameter `use_afni` must be set to True or False in order "
+                "to compute this marker. It is currently set to None (default "
+                "behaviour). This is intended to be for auto-detection. In "
+                "order for that to happen, please call the `validate` method "
+                "before calling the `compute` method."
+            )
+        logger.info("Calculating ReHO map.")
+        # Initialize reho estimator
+        reho_estimator = ReHoEstimator()
+        # Fit-transform reho estimator
+        reho_map = reho_estimator.fit_transform(
+            use_afni=self.use_afni,
+            input_data=input,
+            **reho_params,
+        )
+        return reho_map

junifer/markers/reho/reho_estimator.py

@@ -0,0 +1,510 @@
+"""Provide estimator class for regional homogeneity (ReHo)."""
+
+# Authors: Synchon Mandal <s.mandal@fz-juelich.de>
+#          Federico Raimondo <f.raimondo@fz-juelich.de>
+# License: AGPL
+
+
+import hashlib
+import shutil
+import subprocess
+import tempfile
+from functools import lru_cache
+from itertools import product
+from pathlib import Path
+from typing import TYPE_CHECKING, Any, Dict, List, Optional, cast
+
+import nibabel as nib
+import numpy as np
+from nilearn import image as nimg
+from nilearn import masking as nmask
+from scipy.stats import rankdata
+
+from ...utils import logger, raise_error
+from ..utils import singleton
+
+
+if TYPE_CHECKING:
+    from nibabel import Nifti1Image
+
+
+@singleton
+class ReHoEstimator:
+    """Estimator class for regional homogeneity.
+
+    This class is a singleton and is used for efficient computation of ReHo,
+    by caching the ReHo map for a given set of file path and computation
+    parameters.
+
+    .. warning:: This class can only be used via ReHoBase() and is a deliberate
+                 decision as it serves a specific purpose.
+
+    Attributes
+    ----------
+    temp_dir_path : pathlib.Path
+        Path to the temporary directory for assets storage.
+
+    """
+
+    def __init__(self) -> None:
+        self._file_path = None
+        # Create temporary directory for intermittent storage of assets during
+        # computation via afni's 3dReHo
+        self.temp_dir_path = Path(tempfile.mkdtemp())
+
+    def __del__(self) -> None:
+        """Cleanup."""
+        # Delete temporary directory and ignore errors for read-only files
+        shutil.rmtree(self.temp_dir_path, ignore_errors=True)
+
+    def _compute_reho_afni(
+        self,
+        data: "Nifti1Image",
+        nneigh: int = 27,
+        neigh_rad: Optional[float] = None,
+        neigh_x: Optional[float] = None,
+        neigh_y: Optional[float] = None,
+        neigh_z: Optional[float] = None,
+        box_rad: Optional[int] = None,
+        box_x: Optional[int] = None,
+        box_y: Optional[int] = None,
+        box_z: Optional[int] = None,
+    ) -> "Nifti1Image":
+        """Compute ReHo map via afni's 3dReHo.
+
+        Parameters
+        ----------
+        data : 4D Niimg-like object
+            Images to process.
+        nneigh : {7, 19, 27}, optional
+            Number of voxels in the neighbourhood, inclusive. Can be:
+
+            * 7 : for facewise neighbours only
+            * 19 : for face- and edge-wise nieghbours
+            * 27 : for face-, edge-, and node-wise neighbors
+
+            (default 27).
+        neigh_rad : positive float, optional
+            The radius of a desired neighbourhood (default None).
+        neigh_x : positive float, optional
+            The semi-radius for x-axis of ellipsoidal volumes (default None).
+        neigh_y : positive float, optional
+            The semi-radius for y-axis of ellipsoidal volumes (default None).
+        neigh_z : positive float, optional
+            The semi-radius for z-axis of ellipsoidal volumes (default None).
+        box_rad : positive int, optional
+            The number of voxels outward in a given cardinal direction for a
+            cubic box centered on a given voxel (default None).
+        box_x : positive int, optional
+            The number of voxels for +/- x-axis of cuboidal volumes
+            (default None).
+        box_y : positive int, optional
+            The number of voxels for +/- y-axis of cuboidal volumes
+            (default None).
+        box_z : positive int, optional
+            The number of voxels for +/- z-axis of cuboidal volumes
+            (default None).
+
+        Returns
+        -------
+        Niimg-like object
+
+        Raises
+        ------
+        RuntimeError
+            If the 3dReHo command fails due to some issue.
+
+        Notes
+        -----
+        For more information on the publication, please check [1]_ , and for
+        3dReHo help check:
+        https://afni.nimh.nih.gov/pub/dist/doc/program_help/3dReHo.html
+
+        Please note that that you cannot mix ``box_*`` and ``neigh_*``
+        arguments. The arguments are prioritized by their order in the function
+        signature.
+
+        As the process also depends on the conversion of AFNI files to NIFTI
+        via afni's 3dAFNItoNIFTI, the help for that can be found at:
+        https://afni.nimh.nih.gov/pub/dist/doc/program_help/3dAFNItoNIFTI.html
+
+        References
+        ----------
+        .. [1] Taylor, P.A., & Saad, Z.S. (2013).
+               FATCAT: (An Efficient) Functional And Tractographic Connectivity
+               Analysis Toolbox.
+               Brain connectivity, Volume 3(5), Pages 523-35.
+               https://doi.org/10.1089/brain.2013.0154
+
+        """
+        # Save niimg to nii.gz
+        nifti_in_file_path = self.temp_dir_path / "input.nii"
+        nib.save(data, nifti_in_file_path)
+
+        # Set 3dReHo command
+        reho_afni_out_path_prefix = self.temp_dir_path / "reho"
+        reho_cmd: List[str] = [
+            "3dReHo",
+            f"-prefix {reho_afni_out_path_prefix.resolve()}",
+            f"-inset {nifti_in_file_path.resolve()}",
+        ]
+        # Check ellipsoidal / cuboidal volume arguments
+        if neigh_rad:
+            reho_cmd.append(f"-neigh_RAD {neigh_rad}")
+        elif neigh_x and neigh_y and neigh_z:
+            reho_cmd.extend(
+                [
+                    f"-neigh_X {neigh_x}",
+                    f"-neigh_Y {neigh_y}",
+                    f"-neigh_Z {neigh_z}",
+                ]
+            )
+        elif box_rad:
+            reho_cmd.append(f"-box_RAD {box_rad}")
+        elif box_x and box_y and box_z:
+            reho_cmd.extend(
+                [f"-box_X {box_x}", f"-box_Y {box_y}", f"-box_Z {box_z}"]
+            )
+        else:
+            reho_cmd.append(f"-nneigh {nneigh}")
+        # Call 3dReHo
+        reho_cmd_str = " ".join(reho_cmd)
+        logger.info(f"3dReHo command to be executed: {reho_cmd_str}")
+        reho_process = subprocess.run(
+            reho_cmd_str,  # string needed with shell=True
+            stdin=subprocess.DEVNULL,
+            shell=True,  # needed for respecting $PATH
+            check=False,
+        )
+        if reho_process.returncode == 0:
+            logger.info(
+                "3dReHo succeeded with the following output: "
+                f"{reho_process.stdout}"
+            )
+        else:
+            raise_error(
+                msg="3dReHo failed with the following error: "
+                f"{reho_process.stdout}",
+                klass=RuntimeError,
+            )
+
+        # SHA256 for bypassing memmap
+        sha256_params = hashlib.sha256(bytes(reho_cmd_str, "utf-8"))
+        # Convert afni to nifti
+        reho_afni_to_nifti_out_path = (
+            self.temp_dir_path / f"output_{sha256_params.hexdigest()}.nii"
+        )
+        convert_cmd: List[str] = [
+            "3dAFNItoNIFTI",
+            f"-prefix {reho_afni_to_nifti_out_path.resolve()}",
+            f"{reho_afni_out_path_prefix}+tlrc.BRIK",
+        ]
+        # Call 3dAFNItoNIFTI
+        convert_cmd_str = " ".join(convert_cmd)
+        logger.info(f"3dAFNItoNIFTI command to be executed: {convert_cmd_str}")
+        convert_process = subprocess.run(
+            convert_cmd_str,  # string needed with shell=True
+            stdin=subprocess.DEVNULL,
+            shell=True,  # needed for respecting $PATH
+            check=False,
+        )
+        if convert_process.returncode == 0:
+            logger.info(
+                "3dAFNItoNIFTI succeeded with the following output: "
+                f"{convert_process.stdout}"
+            )
+        else:
+            raise_error(
+                msg="3dAFNItoNIFTI failed with the following error: "
+                f"{convert_process.stdout}",
+                klass=RuntimeError,
+            )
+
+        # Cleanup intermediate files
+        for fname in self.temp_dir_path.glob("reho*"):
+            fname.unlink()
+
+        # Load nifti
+        output_data = nib.load(reho_afni_to_nifti_out_path)
+        # Stupid casting
+        output_data = cast("Nifti1Image", output_data)
+        return output_data
+
+    def _compute_reho_python(
+        self,
+        data: "Nifti1Image",
+        nneigh: int = 27,
+    ) -> "Nifti1Image":
+        """Compute ReHo map.
+
+        Parameters
+        ----------
+        data : 4D Niimg-like object
+            Images to process.
+        nneigh : {7, 19, 27, 125}, optional
+            Number of voxels in the neighbourhood, inclusive. Can be:
+
+            * 7 : for facewise neighbours only
+            * 19 : for face- and edge-wise nieghbours
+            * 27 : for face-, edge-, and node-wise neighbors
+            * 125 : for 5x5 cuboidal volume
+
+            (default 27).
+        Returns
+        -------
+        Niimg-like object
+
+        Raises
+        ------
+        ValueError
+            If ``nneigh`` is invalid.
+
+        """
+        valid_nneigh = (7, 19, 27, 125)
+        if nneigh not in valid_nneigh:
+            raise_error(
+                f"Invalid value for `nneigh`, should be one of {valid_nneigh}."
+            )
+
+        logger.info(f"Computing ReHo map using {nneigh} neighbours.")
+        # Get scan data
+        niimg_data = data.get_fdata()
+        # Get scan dimensions
+        n_x, n_y, n_z, _ = niimg_data.shape
+
+        # Get rank of every voxel across time series
+        ranks_niimg_data = rankdata(niimg_data, axis=-1)
+
+        # Initialize 3D array to store tied rank correction for every voxel
+        tied_rank_corrections = np.zeros((n_x, n_y, n_z), dtype=np.float64)
+        # Calculate tied rank correction for every voxel
+        for i_x, i_y, i_z in product(range(n_x), range(n_y), range(n_z)):
+            # Calculate tied rank count for every voxel across time series
+            _, tie_count = np.unique(
+                ranks_niimg_data[i_x, i_y, i_z, :],
+                return_counts=True,
+            )
+            # Calculate and store tied rank correction for every voxel across
+            # timeseries
+            tied_rank_corrections[i_x, i_y, i_z] = np.sum(
+                tie_count**3 - tie_count
+            )
+
+        # Initialize 3D array to store reho map
+        reho_map = np.ones((n_x, n_y, n_z), dtype=np.float32)
+
+        # Calculate whole brain mask
+        mni152_whole_brain_mask = nmask.compute_brain_mask(
+            data, threshold=0.5, mask_type="whole-brain"
+        )
+        # Convert 0 / 1 array to bool
+        logical_mni152_whole_brain_mask = (
+            mni152_whole_brain_mask.get_fdata().astype(bool)
+        )
+
+        # Create mask cluster and set start and end indices
+        if nneigh in (7, 19, 27):
+            mask_cluster = np.ones((3, 3, 3))
+
+            if nneigh == 7:
+                mask_cluster[0, 0, 0] = 0
+                mask_cluster[0, 1, 0] = 0
+                mask_cluster[0, 2, 0] = 0
+                mask_cluster[0, 0, 1] = 0
+                mask_cluster[0, 2, 1] = 0
+                mask_cluster[0, 0, 2] = 0
+                mask_cluster[0, 1, 2] = 0
+                mask_cluster[0, 2, 2] = 0
+                mask_cluster[1, 0, 0] = 0
+                mask_cluster[1, 2, 0] = 0
+                mask_cluster[1, 0, 2] = 0
+                mask_cluster[1, 2, 2] = 0
+                mask_cluster[2, 0, 0] = 0
+                mask_cluster[2, 1, 0] = 0
+                mask_cluster[2, 2, 0] = 0
+                mask_cluster[2, 0, 1] = 0
+                mask_cluster[2, 2, 1] = 0
+                mask_cluster[2, 0, 2] = 0
+                mask_cluster[2, 1, 2] = 0
+                mask_cluster[2, 2, 2] = 0
+
+            elif nneigh == 19:
+                mask_cluster[0, 0, 0] = 0
+                mask_cluster[0, 2, 0] = 0
+                mask_cluster[2, 0, 0] = 0
+                mask_cluster[2, 2, 0] = 0
+                mask_cluster[0, 0, 2] = 0
+                mask_cluster[0, 2, 2] = 0
+                mask_cluster[2, 0, 2] = 0
+                mask_cluster[2, 2, 2] = 0
+
+            start_idx = 1
+            end_idx = 2
+
+        elif nneigh == 125:
+            mask_cluster = np.ones((5, 5, 5))
+            start_idx = 2
+            end_idx = 3
+
+        # Convert 0 / 1 array to bool
+        logical_mask_cluster = mask_cluster.astype(bool)
+
+        for i, j, k in product(
+            range(start_idx, n_x - (end_idx - 1)),
+            range(start_idx, n_y - (end_idx - 1)),
+            range(start_idx, n_z - (end_idx - 1)),
+        ):
+            # Get mask only for neighbourhood
+            logical_neighbourhood_mni152_whole_brain_mask = (
+                logical_mni152_whole_brain_mask[
+                    i - start_idx : i + end_idx,
+                    j - start_idx : j + end_idx,
+                    k - start_idx : k + end_idx,
+                ]
+            )
+            # Perform logical AND to get neighbourhood mask;
+            # done to take care of brain boundaries
+            neighbourhood_mask = (
+                logical_mask_cluster
+                & logical_neighbourhood_mni152_whole_brain_mask
+            )
+            # Continue if voxel is restricted by mask
+            if neighbourhood_mask[1, 1, 1] == 0:
+                continue
+
+            # Get ranks for the neighbourhood
+            neighbourhood_ranks = ranks_niimg_data[
+                i - start_idx : i + end_idx,
+                j - start_idx : j + end_idx,
+                k - start_idx : k + end_idx,
+                :,
+            ]
+            # Get tied ranks corrections for the neighbourhood
+            neighbourhood_tied_ranks_corrections = tied_rank_corrections[
+                i - start_idx : i + end_idx,
+                j - start_idx : j + end_idx,
+                k - start_idx : k + end_idx,
+            ]
+            # Mask neighbourhood ranks
+            masked_neighbourhood_ranks = neighbourhood_ranks[
+                logical_mask_cluster, :
+            ]
+            # Mask tied ranks corrections for the neighbourhood
+            masked_tied_rank_corrections = (
+                neighbourhood_tied_ranks_corrections[logical_mask_cluster]
+            )
+            # Calculate KCC
+            reho_map[i, j, k] = _kendall_w_reho(
+                timeseries_ranks=masked_neighbourhood_ranks,
+                tied_rank_corrections=masked_tied_rank_corrections,
+            )
+
+        output = nimg.new_img_like(data, reho_map, copy_header=False)
+        return output
+
+    @lru_cache(maxsize=None, typed=True)
+    def _compute(
+        self,
+        use_afni: bool,
+        data: "Nifti1Image",
+        **reho_params: Any,
+    ) -> "Nifti1Image":
+        """Compute the ReHo map with memoization.
+
+        Parameters
+        ----------
+        use_afni : bool
+            Whether to use afni or not.
+        data : 4D Niimg-like object
+            Images to process.
+        **reho_params : dict
+            Extra keyword arguments for ReHo.
+
+        Returns
+        -------
+        Niimg-like object
+
+        """
+        if use_afni:
+            output = self._compute_reho_afni(data, **reho_params)
+        else:
+            output = self._compute_reho_python(data, **reho_params)
+        return output
+
+    def fit_transform(
+        self,
+        use_afni: bool,
+        input_data: Dict[str, Any],
+        **reho_params: Any,
+    ) -> "Nifti1Image":
+        """Fit and transform for the estimator.
+
+        Parameters
+        ----------
+        use_afni : bool
+            Whether to use afni or not.
+        input_data : dict
+            The BOLD data as dictionary.
+        **reho_params : dict
+            Extra keyword arguments for ReHo.
+
+        Returns
+        -------
+        Niimg-like object
+
+        """
+        bold_path = input_data["path"]
+        bold_data = input_data["data"]
+        # Clear cache if file path is different from when caching was done
+        if self._file_path != bold_path:
+            logger.info(f"Removing ReHo map cache at {self._file_path}.")
+            # Clear the cache
+            self._compute.cache_clear()
+            # Clear temporary directory files
+            for file_ in self.temp_dir_path.iterdir():
+                file_.unlink(missing_ok=True)
+            # Set the new file path
+            self._file_path = bold_path
+        else:
+            logger.info(f"Using ReHo map cache at {self._file_path}.")
+        # Compute
+        return self._compute(use_afni, bold_data, **reho_params)
+
+
+def _kendall_w_reho(
+    timeseries_ranks: np.ndarray, tied_rank_corrections: np.ndarray
+) -> float:
+    """Calculate Kendall's coefficient of concordance (KCC) for ReHo map.
+
+    ..note:: This function should only be used to calculate KCC for a ReHo map.
+             For general use, check out ``junifer.stats.kendall_w``.
+
+    Parameters
+    ----------
+    timeseries_matrix : 2D numpy.ndarray
+        A matrix of ranks of a subset subject's brain voxels.
+    tied_rank_corrections : 3D numpy.ndarray
+        A 3D array consisting of the tied rank corrections for the ranks
+        of a subset subject's brain voxels.
+
+    Returns
+    -------
+    float
+        Kendall's W (KCC) of the given timeseries matrix.
+
+    """
+    m, n = timeseries_ranks.shape  # annotators X items
+
+    numerator = (12 * np.sum(np.square(np.sum(timeseries_ranks, axis=0)))) - (
+        3 * m**2 * n * (n + 1) ** 2
+    )
+    denominator = (m**2 * n * (n**2 - 1)) - (
+        m * np.sum(tied_rank_corrections)
+    )
+
+    if denominator == 0:
+        kcc = 1.0
+    else:
+        kcc = numerator / denominator
+
+    return kcc

junifer/markers/reho/reho_parcels.py

@@ -0,0 +1,148 @@
+"""Provide class for regional homogeneity (ReHo) on parcels."""
+
+# Authors: Synchon Mandal <s.mandal@fz-juelich.de>
+# License: AGPL
+
+
+from typing import Any, Dict, Optional
+
+import numpy as np
+
+from ...api.decorators import register_marker
+from ...utils import logger
+from ..parcel_aggregation import ParcelAggregation
+from .reho_base import ReHoBase
+
+
+@register_marker
+class ReHoParcels(ReHoBase):
+    """Class for regional homogeneity on parcels.
+
+    Parameters
+    ----------
+    parcellation : str
+        The name of the parcellation. Check valid options by calling
+        :func:`junifer.data.parcellations.list_parcellations`.
+    use_afni : bool, optional
+        Whether to use AFNI for computing. If None, will use AFNI only
+        if available (default None).
+    reho_params : dict, optional
+        Extra parameters for computing ReHo map as a dictionary (default None).
+        If ``use_afni = True``, then the valid keys are:
+
+        * ``nneigh`` : {7, 19, 27}, optional (default 27)
+            Number of voxels in the neighbourhood, inclusive. Can be:
+
+            - 7 : for facewise neighbours only
+            - 19 : for face- and edge-wise nieghbours
+            - 27 : for face-, edge-, and node-wise neighbors
+
+        * ``neigh_rad`` : positive float, optional
+            The radius of a desired neighbourhood (default None).
+        * ``neigh_x`` : positive float, optional
+            The semi-radius for x-axis of ellipsoidal volumes (default None).
+        * ``neigh_y`` : positive float, optional
+            The semi-radius for y-axis of ellipsoidal volumes (default None).
+        * ``neigh_z`` : positive float, optional
+            The semi-radius for z-axis of ellipsoidal volumes (default None).
+        * ``box_rad`` : positive int, optional
+            The number of voxels outward in a given cardinal direction for a
+            cubic box centered on a given voxel (default None).
+        * ``box_x`` : positive int, optional
+            The number of voxels for +/- x-axis of cuboidal volumes
+            (default None).
+        * ``box_y`` : positive int, optional
+            The number of voxels for +/- y-axis of cuboidal volumes
+            (default None).
+        * ``box_z`` : positive int, optional
+            The number of voxels for +/- z-axis of cuboidal volumes
+            (default None).
+
+        else if ``use_afni = False``, then the valid keys are:
+
+        * ``nneigh`` : {7, 19, 27, 125}, optional (default 27)
+            Number of voxels in the neighbourhood, inclusive. Can be:
+
+            * 7 : for facewise neighbours only
+            * 19 : for face- and edge-wise nieghbours
+            * 27 : for face-, edge-, and node-wise neighbors
+            * 125 : for 5x5 cuboidal volume
+
+    agg_method : str, optional
+        The method to perform aggregation using. Check valid options in
+        :func:`junifer.stats.get_aggfunc_by_name` (default "mean").
+    agg_method_params : dict, optional
+        Parameters to pass to the aggregation function. Check valid options in
+        :func:`junifer.stats.get_aggfunc_by_name` (default None).
+    mask : str, optional
+        The name of the mask to apply to regions before extracting signals.
+        Check valid options by calling :func:`junifer.data.masks.list_masks`
+        (default None).
+    name : str, optional
+        The name of the marker. If None, it will use the class name
+        (default None).
+
+    """
+
+    def __init__(
+        self,
+        parcellation: str,
+        use_afni: Optional[bool] = None,
+        reho_params: Optional[Dict] = None,
+        agg_method: str = "mean",
+        agg_method_params: Optional[Dict] = None,
+        mask: Optional[str] = None,
+        name: Optional[str] = None,
+    ) -> None:
+        self.parcellation = parcellation
+        self.reho_params = reho_params
+        self.agg_method = agg_method
+        self.agg_method_params = agg_method_params
+        self.mask = mask
+        super().__init__(use_afni=use_afni, name=name)
+
+    def compute(
+        self,
+        input: Dict[str, Any],
+        extra_input: Optional[Dict[str, Any]] = None,
+    ) -> Dict[str, Any]:
+        """Compute.
+
+        Parameters
+        ----------
+        input : dict
+            The BOLD data as dictionary.
+        extra_input : dict, optional
+            The other fields in the pipeline data object (default None).
+
+        Returns
+        -------
+        dict
+            The computed result as dictionary. The dictionary has the following
+            keys:
+
+            * ``data`` : the actual computed values as a 1D numpy.ndarray
+            * ``columns`` : the column labels for the parcels as a list
+            * ``row_names`` : ``None``
+
+        """
+        logger.info("Calculating ReHo for parcels.")
+        # Calculate reho map
+        if self.reho_params is not None:
+            reho_map = self.compute_reho_map(input=input, **self.reho_params)
+        else:
+            reho_map = self.compute_reho_map(input=input)
+        # Initialize parcel aggregation
+        parcel_aggregation = ParcelAggregation(
+            parcellation=self.parcellation,
+            method=self.agg_method,
+            method_params=self.agg_method_params,
+            mask=self.mask,
+            on="BOLD",
+        )
+        # Perform aggregation on reho map
+        parcel_aggregation_input = {"data": reho_map}
+        output = parcel_aggregation.compute(input=parcel_aggregation_input)
+        # Only use the first row and expand row dimension
+        output["data"] = output["data"][0][np.newaxis, :]
+        return output

junifer/markers/reho/reho_spheres.py

@@ -0,0 +1,156 @@
+"""Provide class for regional homogeneity (ReHo) on spheres."""
+
+# Authors: Synchon Mandal <s.mandal@fz-juelich.de>
+# License: AGPL
+
+
+from typing import Any, Dict, Optional
+
+import numpy as np
+
+from ...api.decorators import register_marker
+from ...utils import logger
+from ..sphere_aggregation import SphereAggregation
+from .reho_base import ReHoBase
+
+
+@register_marker
+class ReHoSpheres(ReHoBase):
+    """Class for regional homogeneity on spheres.
+
+    Parameters
+    ----------
+    coords : str
+        The name of the coordinates list to use. See
+        :func:`junifer.data.coordinates.list_coordinates` for options.
+    radius : float, optional
+        The radius of the sphere in millimeters. If None, the signal will be
+        extracted from a single voxel. See
+        :class:`nilearn.maskers.NiftiSpheresMasker` for more information
+        (default None).
+    use_afni : bool, optional
+        Whether to use AFNI for computing. If None, will use AFNI only
+        if available (default None).
+    reho_params : dict, optional
+        Extra parameters for computing ReHo map as a dictionary (default None).
+        If ``use_afni = True``, then the valid keys are:
+
+        * ``nneigh`` : {7, 19, 27}, optional (default 27)
+            Number of voxels in the neighbourhood, inclusive. Can be:
+
+            - 7 : for facewise neighbours only
+            - 19 : for face- and edge-wise nieghbours
+            - 27 : for face-, edge-, and node-wise neighbors
+
+        * ``neigh_rad`` : positive float, optional
+            The radius of a desired neighbourhood (default None).
+        * ``neigh_x`` : positive float, optional
+            The semi-radius for x-axis of ellipsoidal volumes (default None).
+        * ``neigh_y`` : positive float, optional
+            The semi-radius for y-axis of ellipsoidal volumes (default None).
+        * ``neigh_z`` : positive float, optional
+            The semi-radius for z-axis of ellipsoidal volumes (default None).
+        * ``box_rad`` : positive int, optional
+            The number of voxels outward in a given cardinal direction for a
+            cubic box centered on a given voxel (default None).
+        * ``box_x`` : positive int, optional
+            The number of voxels for +/- x-axis of cuboidal volumes
+            (default None).
+        * ``box_y`` : positive int, optional
+            The number of voxels for +/- y-axis of cuboidal volumes
+            (default None).
+        * ``box_z`` : positive int, optional
+            The number of voxels for +/- z-axis of cuboidal volumes
+            (default None).
+
+        else if ``use_afni = False``, then the valid keys are:
+
+        * ``nneigh`` : {7, 19, 27, 125}, optional (default 27)
+            Number of voxels in the neighbourhood, inclusive. Can be:
+
+            * 7 : for facewise neighbours only
+            * 19 : for face- and edge-wise nieghbours
+            * 27 : for face-, edge-, and node-wise neighbors
+            * 125 : for 5x5 cuboidal volume
+
+    agg_method : str, optional
+        The aggregation method to use.
+        See :func:`junifer.stats.get_aggfunc_by_name` for more information
+        (default None).
+    agg_method_params : dict, optional
+        The parameters to pass to the aggregation method (default None).
+    mask : str, optional
+        The name of the mask to apply to regions before extracting signals.
+        Check valid options by calling :func:`junifer.data.masks.list_masks`
+        (default None).
+    name : str, optional
+        The name of the marker. If None, it will use the class name
+        (default None).
+
+    """
+
+    def __init__(
+        self,
+        coords: str,
+        radius: Optional[float] = None,
+        use_afni: Optional[bool] = None,
+        reho_params: Optional[Dict] = None,
+        agg_method: str = "mean",
+        agg_method_params: Optional[Dict] = None,
+        mask: Optional[str] = None,
+        name: Optional[str] = None,
+    ) -> None:
+        self.coords = coords
+        self.radius = radius
+        self.reho_params = reho_params
+        self.agg_method = agg_method
+        self.agg_method_params = agg_method_params
+        self.mask = mask
+        super().__init__(use_afni=use_afni, name=name)
+
+    def compute(
+        self,
+        input: Dict[str, Any],
+        extra_input: Optional[Dict[str, Any]] = None,
+    ) -> Dict[str, Any]:
+        """Compute.
+
+        Parameters
+        ----------
+        input : dict
+            The BOLD data as dictionary.
+        extra_input : dict, optional
+            The other fields in the pipeline data object (default None).
+
+        Returns
+        -------
+        dict
+            The computed result as dictionary. The dictionary has the following
+            keys:
+
+            * ``data`` : the actual computed values as a 1D numpy.ndarray
+            * ``columns`` : the column labels for the spheres as a list
+            * ``rows_col_name`` : ``None``
+
+        """
+        logger.info("Calculating ReHo for spheres.")
+        # Calculate reho map
+        if self.reho_params is not None:
+            reho_map = self.compute_reho_map(input=input, **self.reho_params)
+        else:
+            reho_map = self.compute_reho_map(input=input)
+        # Initialize sphere aggregation
+        sphere_aggregation = SphereAggregation(
+            coords=self.coords,
+            radius=self.radius,
+            method=self.agg_method,
+            method_params=self.agg_method_params,
+            mask=self.mask,
+            on="BOLD",
+        )
+        # Perform aggregation on reho map
+        sphere_aggregation_input = {"data": reho_map}
+        output = sphere_aggregation.compute(input=sphere_aggregation_input)
+        # Only use the first row and expand row dimension
+        output["data"] = output["data"][0][np.newaxis, :]
+        return output

junifer/markers/reho/tests/test_reho_estimator.py

@@ -0,0 +1,260 @@
+"""Provide tests for ReHo map compute comparison."""
+
+# Authors: Synchon Mandal <s.mandal@fz-juelich.de>
+# License: AGPL
+
+import time
+
+import nibabel as nib
+import pytest
+from scipy.stats import pearsonr
+
+from junifer.datareader.default import DefaultDataReader
+from junifer.markers.reho.reho_estimator import ReHoEstimator
+from junifer.pipeline.utils import _check_afni
+from junifer.testing.datagrabbers import PartlyCloudyTestingDataGrabber
+from junifer.utils.logging import logger
+
+
+def test_reho_estimator_cache_python() -> None:
+    """Test that the cache works properly when using Python implementation."""
+    # Get subject from datagrabber
+    with PartlyCloudyTestingDataGrabber() as dg:
+        subject = dg["sub-01"]
+    # Read data for subject
+    subject_data = DefaultDataReader().fit_transform(subject)
+    # Setup estimator
+    reho_estimator = ReHoEstimator()
+
+    first_tic = time.time()
+    reho_map_without_cache = reho_estimator.fit_transform(
+        use_afni=False,
+        input_data=subject_data["BOLD"],
+        nneigh=27,
+    )
+    first_toc = time.time()
+    logger.info(
+        f"ReHo estimator in Python without cache: {first_toc - first_tic}"
+    )
+    assert isinstance(reho_map_without_cache, nib.Nifti1Image)
+    # Count intermediate files
+    n_files = len([x for x in reho_estimator.temp_dir_path.glob("*")])
+    assert n_files == 0  # no files in python
+
+    # Now fit again, should be faster
+    second_tic = time.time()
+    reho_map_with_cache = reho_estimator.fit_transform(
+        use_afni=False,
+        input_data=subject_data["BOLD"],
+        nneigh=27,
+    )
+    second_toc = time.time()
+    logger.info(
+        f"ReHo estimator in Python with cache: {second_toc - second_tic}"
+    )
+    assert isinstance(reho_map_with_cache, nib.Nifti1Image)
+    # Check that cache is being used
+    assert (second_toc - second_tic) < ((first_toc - first_tic) / 1000)
+    # Count intermediate files
+    n_files = len([x for x in reho_estimator.temp_dir_path.glob("*")])
+    assert n_files == 0  # no files in python
+
+    # Now change a parameter, should compute again, without clearing the
+    # cache
+    third_tic = time.time()
+    reho_map_with_partial_cache = reho_estimator.fit_transform(
+        use_afni=False,
+        input_data=subject_data["BOLD"],
+        nneigh=125,
+    )
+    third_toc = time.time()
+    logger.info(
+        f"ReHo estimator in Python with partial cache: {third_toc - third_tic}"
+    )
+    assert isinstance(reho_map_with_partial_cache, nib.Nifti1Image)
+    # Should require more time
+    assert (third_toc - third_tic) > ((first_toc - first_tic) / 10)
+    # Count intermediate files
+    n_files = len([x for x in reho_estimator.temp_dir_path.glob("*")])
+    assert n_files == 0  # no files in python
+
+    # Now fit again with the previous params, should be fast
+    fourth_tic = time.time()
+    reho_map_with_new_cache = reho_estimator.fit_transform(
+        use_afni=False,
+        input_data=subject_data["BOLD"],
+        nneigh=125,
+    )
+    fourth_toc = time.time()
+    logger.info(
+        f"ReHo estimator in Python with new cache: {fourth_toc - fourth_tic}"
+    )
+    assert isinstance(reho_map_with_new_cache, nib.Nifti1Image)
+    # Should require less time
+    assert (fourth_toc - fourth_tic) < ((first_toc - first_tic) / 1000)
+    # Count intermediate files
+    n_files = len([x for x in reho_estimator.temp_dir_path.glob("*")])
+    assert n_files == 0  # no files in python
+
+    # Now change the data, it should clear the cache
+    with PartlyCloudyTestingDataGrabber() as dg:
+        subject = dg["sub-02"]
+    # Read data for new subject
+    subject_data = DefaultDataReader().fit_transform(subject)
+
+    fifth_tic = time.time()
+    reho_map_with_different_cache = reho_estimator.fit_transform(
+        use_afni=False,
+        input_data=subject_data["BOLD"],
+        nneigh=27,
+    )
+    fifth_toc = time.time()
+    logger.info(
+        "ReHo estimator in Python with different cache: "
+        f"{fifth_toc - fifth_tic}"
+    )
+    assert isinstance(reho_map_with_different_cache, nib.Nifti1Image)
+    # Should take less time
+    assert (fifth_toc - fifth_tic) > ((first_toc - first_tic) / 10)
+    # Count intermediate files
+    n_files = len([x for x in reho_estimator.temp_dir_path.glob("*")])
+    assert n_files == 0  # no files in python
+
+
+@pytest.mark.skipif(
+    _check_afni() is False, reason="requires afni to be in PATH"
+)
+def test_reho_estimator_cache_afni() -> None:
+    """Test that the cache works properly when using afni."""
+    # Get subject from datagrabber
+    with PartlyCloudyTestingDataGrabber() as dg:
+        subject = dg["sub-01"]
+    # Read data for subject
+    subject_data = DefaultDataReader().fit_transform(subject)
+    # Setup estimator
+    reho_estimator = ReHoEstimator()
+
+    first_tic = time.time()
+    reho_map_without_cache = reho_estimator.fit_transform(
+        use_afni=True,
+        input_data=subject_data["BOLD"],
+        nneigh=19,
+    )
+    first_toc = time.time()
+    logger.info(
+        f"ReHo estimator in AFNI without cache: {first_toc - first_tic}"
+    )
+    assert isinstance(reho_map_without_cache, nib.Nifti1Image)
+    # Count intermediate files
+    n_files = len([x for x in reho_estimator.temp_dir_path.glob("*")])
+    assert n_files == 2  # input + reho
+
+    # Now fit again, should be faster
+    second_tic = time.time()
+    reho_map_with_cache = reho_estimator.fit_transform(
+        use_afni=True,
+        input_data=subject_data["BOLD"],
+        nneigh=19,
+    )
+    second_toc = time.time()
+    logger.info(
+        f"ReHo estimator in AFNI with cache: {second_toc - second_tic}"
+    )
+    assert isinstance(reho_map_with_cache, nib.Nifti1Image)
+    assert (second_toc - second_tic) < ((first_toc - first_tic) / 1000)
+    # Count intermediate files
+    n_files = len([x for x in reho_estimator.temp_dir_path.glob("*")])
+    assert n_files == 2  # input + reho
+
+    # Now change a parameter, should compute again, without clearing the
+    # cache
+    third_tic = time.time()
+    reho_map_with_partial_cache = reho_estimator.fit_transform(
+        use_afni=True,
+        input_data=subject_data["BOLD"],
+        nneigh=27,
+    )
+    third_toc = time.time()
+    logger.info(
+        f"ReHo estimator in AFNI with partial cache: {third_toc - third_tic}"
+    )
+    assert isinstance(reho_map_with_partial_cache, nib.Nifti1Image)
+    # Should require more time
+    assert (third_toc - third_tic) > ((first_toc - first_tic) / 10)
+    # Count intermediate files
+    n_files = len([x for x in reho_estimator.temp_dir_path.glob("*")])
+    assert n_files == 3  # input + 2 * reho
+
+    # Now fit again with the previous params, should be fast
+    fourth_tic = time.time()
+    reho_map_with_new_cache = reho_estimator.fit_transform(
+        use_afni=True,
+        input_data=subject_data["BOLD"],
+        nneigh=27,
+    )
+    fourth_toc = time.time()
+    logger.info(
+        f"ReHo estimator in AFNI with new cache: {fourth_toc - fourth_tic}"
+    )
+    assert isinstance(reho_map_with_new_cache, nib.Nifti1Image)
+    # Should require less time
+    assert (fourth_toc - fourth_tic) < ((first_toc - first_tic) / 1000)
+    n_files = len([x for x in reho_estimator.temp_dir_path.glob("*")])
+    assert n_files == 3  # input + 2 * reho
+
+    # Now change the data, it should clear the cache
+    with PartlyCloudyTestingDataGrabber() as dg:
+        subject = dg["sub-02"]
+    # Read data for new subject
+    subject_data = DefaultDataReader().fit_transform(subject)
+
+    fifth_tic = time.time()
+    reho_map_with_different_cache = reho_estimator.fit_transform(
+        use_afni=True,
+        input_data=subject_data["BOLD"],
+        nneigh=27,
+    )
+    fifth_toc = time.time()
+    logger.info(
+        f"ReHo estimator in AFNI with different cache: {fifth_toc - fifth_tic}"
+    )
+    assert isinstance(reho_map_with_different_cache, nib.Nifti1Image)
+    # Should take less time
+    assert (fifth_toc - fifth_tic) > ((first_toc - first_tic) / 10)
+    # Count intermediate files
+    n_files = len([x for x in reho_estimator.temp_dir_path.glob("*")])
+    assert n_files == 2  # input + reho
+
+
+@pytest.mark.skipif(
+    _check_afni() is False, reason="requires afni to be in PATH"
+)
+def test_reho_estimator_afni_vs_python() -> None:
+    """Compare afni and Python implementations."""
+    # Get subject from datagrabber
+    with PartlyCloudyTestingDataGrabber() as dg:
+        subject = dg["sub-01"]
+    # Read data for subject
+    subject_data = DefaultDataReader().fit_transform(subject)
+    # Setup estimator
+    reho_estimator = ReHoEstimator()
+
+    # Compare using 27 neighbours
+    reho_map_afni = reho_estimator.fit_transform(
+        use_afni=True,
+        input_data=subject_data["BOLD"],
+        nneigh=27,
+    )
+    reho_map_python = reho_estimator.fit_transform(
+        use_afni=False,
+        input_data=subject_data["BOLD"],
+        nneigh=27,
+    )
+
+    # Calculate Pearson correlation coefficient
+    r, _ = pearsonr(
+        reho_map_afni.get_fdata().flatten(),
+        reho_map_python.get_fdata().flatten(),
+    )
+    # Assert good correlation
+    assert r > 0.70

junifer/markers/reho/tests/test_reho_parcels.py

@@ -0,0 +1,118 @@
+"""Provide tests for ReHo on parcels."""
+
+# Authors: Synchon Mandal <s.mandal@fz-juelich.de>
+# License: AGPL
+
+from pathlib import Path
+
+import pytest
+from nilearn import image as nimg
+from scipy.stats import pearsonr
+
+from junifer.markers.reho.reho_parcels import ReHoParcels
+from junifer.pipeline.utils import _check_afni
+from junifer.storage.sqlite import SQLiteFeatureStorage
+from junifer.testing.datagrabbers import SPMAuditoryTestingDatagrabber
+
+
+PARCELLATION = "Schaefer100x7"
+
+
+def test_reho_parcels_computation() -> None:
+    """Test ReHoParcels fit-transform."""
+    with SPMAuditoryTestingDatagrabber() as dg:
+        # Use first subject
+        subject_data = dg["sub001"]
+        # Load image to memory
+        fmri_img = nimg.load_img(subject_data["BOLD"]["path"])
+        # Initialize marker
+        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": {}}}
+        )
+        # Get BOLD output
+        reho_parcels_output_bold = reho_parcels_output["BOLD"]
+        # Assert BOLD output keys
+        assert "data" in reho_parcels_output_bold
+        assert "columns" in reho_parcels_output_bold
+
+        reho_parcels_output_bold_data = reho_parcels_output_bold["data"]
+        # Assert BOLD output data dimension
+        assert reho_parcels_output_bold_data.ndim == 2
+        # Assert BOLD output data is normalized
+        assert (reho_parcels_output_bold_data > 0).all() and (
+            reho_parcels_output_bold_data < 1
+        ).all()
+
+
+@pytest.mark.skipif(
+    _check_afni() is False, reason="requires afni to be in PATH"
+)
+def test_reho_parcels_computation_comparison() -> None:
+    """Test ReHoParcels fit-transform implementation comparison.."""
+    with SPMAuditoryTestingDatagrabber() as dg:
+        # Use first subject
+        subject_data = dg["sub001"]
+        # Load image to memory
+        fmri_img = nimg.load_img(subject_data["BOLD"]["path"])
+
+        # Initialize marker with use_afni=False
+        reho_parcels_marker_python = ReHoParcels(
+            parcellation=PARCELLATION, use_afni=False
+        )
+        # Fit transform marker on data
+        reho_parcels_output_python = reho_parcels_marker_python.fit_transform(
+            {"BOLD": {"path": "/tmp", "data": fmri_img, "meta": {}}}
+        )
+        # Get BOLD output
+        reho_parcels_output_bold_python = reho_parcels_output_python["BOLD"]
+
+        # Initialize marker with use_afni=True
+        reho_parcels_marker_afni = ReHoParcels(
+            parcellation=PARCELLATION, use_afni=True
+        )
+        # Fit transform marker on data
+        reho_parcels_output_afni = reho_parcels_marker_afni.fit_transform(
+            {"BOLD": {"path": "/tmp", "data": fmri_img, "meta": {}}}
+        )
+        # Get BOLD output
+        reho_parcels_output_bold_afni = reho_parcels_output_afni["BOLD"]
+
+        # Check for Pearson correlation coefficient
+        r, _ = pearsonr(
+            reho_parcels_output_bold_python["data"].flatten(),
+            reho_parcels_output_bold_afni["data"].flatten(),
+        )
+        assert r >= 0.3  # this is very bad, but they differ...
+
+
+def test_reho_parcels_storage(tmp_path: Path) -> None:
+    """Test ReHoParcels storage.
+
+    Parameters
+    ----------
+    tmp_path : pathlib.Path
+        The path to the test directory.
+
+    """
+    with SPMAuditoryTestingDatagrabber() as dg:
+        # Use first subject
+        subject_data = dg["sub001"]
+        # Load image to memory
+        fmri_img = nimg.load_img(subject_data["BOLD"]["path"])
+        # Initialize marker
+        reho_parcels_marker = ReHoParcels(parcellation=PARCELLATION)
+        # Initialize storage
+        reho_parcels_storage = SQLiteFeatureStorage(
+            tmp_path / "reho_parcels.sqlite"
+        )
+        # Generate meta
+        meta = {
+            "element": {"subject": "sub001"}
+        }  # 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}},
+            storage=reho_parcels_storage,
+        )

junifer/markers/reho/tests/test_reho_spheres.py

@@ -0,0 +1,118 @@
+"""Provide tests for ReHo on spheres."""
+
+# Authors: Synchon Mandal <s.mandal@fz-juelich.de>
+# License: AGPL
+
+from pathlib import Path
+
+import pytest
+from nilearn import image as nimg
+from scipy.stats import pearsonr
+
+from junifer.markers.reho.reho_spheres import ReHoSpheres
+from junifer.pipeline.utils import _check_afni
+from junifer.storage.sqlite import SQLiteFeatureStorage
+from junifer.testing.datagrabbers import SPMAuditoryTestingDatagrabber
+
+
+COORDINATES = "DMNBuckner"
+
+
+def test_reho_spheres_computation() -> None:
+    """Test ReHoSpheres fit-transform."""
+    with SPMAuditoryTestingDatagrabber() as dg:
+        # Use first subject
+        subject_data = dg["sub001"]
+        # Load image to memory
+        fmri_img = nimg.load_img(subject_data["BOLD"]["path"])
+        # Initialize marker
+        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": {}}}
+        )
+        # Get BOLD output
+        reho_spheres_output_bold = reho_spheres_output["BOLD"]
+        # Assert BOLD output keys
+        assert "data" in reho_spheres_output_bold
+        assert "columns" in reho_spheres_output_bold
+
+        reho_spheres_output_bold_data = reho_spheres_output_bold["data"]
+        # Assert BOLD output data dimension
+        assert reho_spheres_output_bold_data.ndim == 2
+        # Assert BOLD output data is normalized
+        assert (reho_spheres_output_bold_data > 0).all() and (
+            reho_spheres_output_bold_data < 1
+        ).all()
+
+
+@pytest.mark.skipif(
+    _check_afni() is False, reason="requires afni to be in PATH"
+)
+def test_reho_spheres_computation_comparison() -> None:
+    """Test ReHoSpheres fit-transform implementation comparison.."""
+    with SPMAuditoryTestingDatagrabber() as dg:
+        # Use first subject
+        subject_data = dg["sub001"]
+        # Load image to memory
+        fmri_img = nimg.load_img(subject_data["BOLD"]["path"])
+
+        # Initialize marker with use_afni=False
+        reho_spheres_marker_python = ReHoSpheres(
+            coords=COORDINATES, radius=10.0, use_afni=False
+        )
+        # Fit transform marker on data
+        reho_spheres_output_python = reho_spheres_marker_python.fit_transform(
+            {"BOLD": {"path": "/tmp", "data": fmri_img, "meta": {}}}
+        )
+        # Get BOLD output
+        reho_spheres_output_bold_python = reho_spheres_output_python["BOLD"]
+
+        # Initialize marker with use_afni=True
+        reho_spheres_marker_afni = ReHoSpheres(
+            coords=COORDINATES, radius=10.0, use_afni=True
+        )
+        # Fit transform marker on data
+        reho_spheres_output_afni = reho_spheres_marker_afni.fit_transform(
+            {"BOLD": {"path": "/tmp", "data": fmri_img, "meta": {}}}
+        )
+        # Get BOLD output
+        reho_spheres_output_bold_afni = reho_spheres_output_afni["BOLD"]
+
+        # Check for Pearson correlation coefficient
+        r, _ = pearsonr(
+            reho_spheres_output_bold_python["data"].flatten(),
+            reho_spheres_output_bold_afni["data"].flatten(),
+        )
+        assert r >= 0.8  # 0.8 is a loose threshold
+
+
+def test_reho_spheres_storage(tmp_path: Path) -> None:
+    """Test ReHoSpheres storage.
+
+    Parameters
+    ----------
+    tmp_path : pathlib.Path
+        The path to the test directory.
+
+    """
+    with SPMAuditoryTestingDatagrabber() as dg:
+        # Use first subject
+        subject_data = dg["sub001"]
+        # Load image to memory
+        fmri_img = nimg.load_img(subject_data["BOLD"]["path"])
+        # Initialize marker
+        reho_spheres_marker = ReHoSpheres(coords=COORDINATES, radius=10.0)
+        # Initialize storage
+        reho_spheres_storage = SQLiteFeatureStorage(
+            tmp_path / "reho_spheres.sqlite"
+        )
+        # Generate meta
+        meta = {
+            "element": {"subject": "sub001"}
+        }  # 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}},
+            storage=reho_spheres_storage,
+        )

junifer/markers/tests/test_crossparcellation_functional_connectivity.py

@@ -71,8 +71,7 @@ def test_store(tmp_path: Path) -> None:
         crossparcellation.fit_transform(input_dict, storage=storage)
         features = storage.list_features()
         assert any(
-            x["name"] == "BOLD_CrossParcellationFC"
+            x["name"] == "BOLD_CrossParcellationFC" for x in features.values()
-            for x in features.values()
         )
 
 

junifer/markers/tests/test_sphere_aggregation.py

@@ -4,8 +4,8 @@
 # License: AGPL
 
 import typing
-from typing import Dict
 from pathlib import Path
+from typing import Dict
 
 import nibabel as nib
 import pytest

junifer/markers/utils.py

@@ -29,7 +29,7 @@ def singleton(cls: Type) -> Type:
     class
         The only instance of the class.
 
-    """ ""
+    """
     instances: Dict = {}
 
     def get_instance(*args: Any, **kwargs: Any) -> Type:
@@ -47,7 +47,7 @@ def singleton(cls: Type) -> Type:
         class
             The only instance of the class.
 
-        """ ""
+        """
         if cls not in instances:
             instances[cls] = cls(*args, **kwargs)
         return instances[cls]

junifer/pipeline/tests/test_pipeline_step_mixin.py

@@ -4,10 +4,10 @@
 #          Synchon Mandal <s.mandal@fz-juelich.de>
 # License: AGPL
 
+import warnings
 from typing import Dict, List
 
 import pytest
-import warnings
 
 from junifer.pipeline.pipeline_step_mixin import PipelineStepMixin
 from junifer.pipeline.utils import _check_afni

junifer/storage/sqlite.py

@@ -4,10 +4,10 @@
 #          Synchon Mandal <s.mandal@fz-juelich.de>
 # License: AGPL
 
+import json
 from pathlib import Path
 from typing import TYPE_CHECKING, Dict, List, Optional, Union
 
-import json
 import numpy as np
 import pandas as pd
 from pandas.core.base import NoNewAttributesMixin

tox.ini

@@ -93,6 +93,10 @@ exclude =
     __init__.py
 max-line-length = 79
 extend-ignore =
+    ; Use of `functools.lru_cache` or `functools.cache` on methods can lead to
+    ; memory leaks. The cache may retain instance references, preventing garbage
+    ; collection.
+    B019
     ; abstract class with no abstract methods
     B024
     D202