docs/changes/newsfragments/311.change

@@ -0,0 +1 @@
+Add a positional argument ``using`` for Markers and Preprocessors having implementation-based variations, in particular :class:`.ReHoParcels`, :class:`.ReHoSpheres`, :class:`.ALFFParcels`, :class:`.ALFFSpheres` and :class:`.BOLDWarper` by `Synchon Mandal`_

docs/changes/newsfragments/311.enh

@@ -0,0 +1 @@
+Improve external dependency handling for :class:`.PipelineStepMixin`-derived objects having implementation-based variations by `Synchon Mandal`_

junifer/markers/falff/_afni_falff.py

@@ -0,0 +1,153 @@
+"""Provide class for computing ALFF using AFNI."""
+
+# Authors: Synchon Mandal <s.mandal@fz-juelich.de>
+# License: AGPL
+
+from functools import lru_cache
+from pathlib import Path
+from typing import (
+    TYPE_CHECKING,
+    ClassVar,
+    Dict,
+    List,
+    Optional,
+    Tuple,
+    Union,
+)
+
+import nibabel as nib
+
+from ...pipeline import WorkDirManager
+from ...pipeline.singleton import singleton
+from ...utils import logger, run_ext_cmd
+
+
+if TYPE_CHECKING:
+    from nibabel import Nifti1Image
+
+
+@singleton
+class AFNIALFF:
+    """Class for computing ALFF using AFNI.
+
+    This class uses AFNI's 3dRSFC to compute ALFF. It's designed as a singleton
+    with caching for efficient computation.
+
+    """
+
+    _EXT_DEPENDENCIES: ClassVar[List[Dict[str, Union[str, List[str]]]]] = [
+        {
+            "name": "afni",
+            "commands": ["3dRSFC", "3dAFNItoNIFTI"],
+        },
+    ]
+
+    def __del__(self) -> None:
+        """Terminate the class."""
+        # Clear the computation cache
+        logger.debug("Clearing cache for ALFF computation via AFNI")
+        self.compute.cache_clear()
+
+    @lru_cache(maxsize=None, typed=True)
+    def compute(
+        self,
+        data: "Nifti1Image",
+        highpass: float,
+        lowpass: float,
+        tr: Optional[float],
+    ) -> Tuple["Nifti1Image", "Nifti1Image", Path, Path]:
+        """Compute ALFF + fALFF map.
+
+        Parameters
+        ----------
+        data : 4D Niimg-like object
+            Images to process.
+        highpass : positive float
+            Highpass cutoff frequency.
+        lowpass : positive float
+            Lowpass cutoff frequency.
+        tr : positive float, optional
+            The Repetition Time of the BOLD data.
+
+        Returns
+        -------
+        Niimg-like object
+            ALFF map.
+        Niimg-like object
+            fALFF map.
+        pathlib.Path
+            The path to the ALFF map as NIfTI.
+        pathlib.Path
+            The path to the fALFF map as NIfTI.
+
+        """
+        logger.debug("Creating cache for ALFF computation via AFNI")
+
+        # Create component-scoped tempdir
+        tempdir = WorkDirManager().get_tempdir(prefix="afni_alff+falff")
+
+        # Save target data to a component-scoped tempfile
+        nifti_in_file_path = tempdir / "input.nii"  # needs to be .nii
+        nib.save(data, nifti_in_file_path)
+
+        # Set 3dRSFC command
+        alff_falff_out_path_prefix = tempdir / "alff_falff"
+        bp_cmd = [
+            "3dRSFC",
+            f"-prefix {alff_falff_out_path_prefix.resolve()}",
+            f"-input {nifti_in_file_path.resolve()}",
+            f"-band {highpass} {lowpass}",
+            "-no_rsfa -nosat -nodetrend",
+        ]
+        # Check tr
+        if tr is not None:
+            bp_cmd.append(f"-dt {tr}")
+        # Call 3dRSFC
+        run_ext_cmd(name="3dRSFC", cmd=bp_cmd)
+
+        # Create element-scoped tempdir so that the ALFF and fALFF maps are
+        # available later as nibabel stores file path reference for
+        # loading on computation
+        element_tempdir = WorkDirManager().get_element_tempdir(
+            prefix="afni_alff_falff"
+        )
+
+        params_suffix = f"_{highpass}_{lowpass}_{tr}"
+
+        # Convert alff afni to nifti
+        alff_afni_to_nifti_out_path = (
+            element_tempdir / f"alff{params_suffix}_output.nii"
+        )  # needs to be .nii
+        convert_alff_cmd = [
+            "3dAFNItoNIFTI",
+            f"-prefix {alff_afni_to_nifti_out_path.resolve()}",
+            f"{alff_falff_out_path_prefix}_ALFF+tlrc.BRIK",
+        ]
+        # Call 3dAFNItoNIFTI
+        run_ext_cmd(name="3dAFNItoNIFTI", cmd=convert_alff_cmd)
+
+        # Convert falff afni to nifti
+        falff_afni_to_nifti_out_path = (
+            element_tempdir / f"falff{params_suffix}_output.nii"
+        )  # needs to be .nii
+        convert_falff_cmd = [
+            "3dAFNItoNIFTI",
+            f"-prefix {falff_afni_to_nifti_out_path.resolve()}",
+            f"{alff_falff_out_path_prefix}_fALFF+tlrc.BRIK",
+        ]
+        # Call 3dAFNItoNIFTI
+        run_ext_cmd(name="3dAFNItoNIFTI", cmd=convert_falff_cmd)
+
+        # Load nifti
+        alff_data = nib.load(alff_afni_to_nifti_out_path)
+        falff_data = nib.load(falff_afni_to_nifti_out_path)
+
+        # Delete tempdir
+        WorkDirManager().delete_tempdir(tempdir)
+
+        return (
+            alff_data,
+            falff_data,
+            alff_afni_to_nifti_out_path,
+            falff_afni_to_nifti_out_path,
+        )  # type: ignore

junifer/markers/falff/_junifer_falff.py

@@ -0,0 +1,142 @@
+"""Provide class for computing ALFF using junifer."""
+
+# Authors: Synchon Mandal <s.mandal@fz-juelich.de>
+# License: AGPL
+
+from functools import lru_cache
+from pathlib import Path
+from typing import (
+    TYPE_CHECKING,
+    ClassVar,
+    Optional,
+    Set,
+    Tuple,
+)
+
+import nibabel as nib
+import numpy as np
+import scipy as sp
+from nilearn import image as nimg
+
+from ...pipeline import WorkDirManager
+from ...pipeline.singleton import singleton
+from ...utils import logger
+
+
+if TYPE_CHECKING:
+    from nibabel import Nifti1Image
+
+
+@singleton
+class JuniferALFF:
+    """Class for computing ALFF using junifer.
+
+    It's designed as a singleton with caching for efficient computation.
+
+    """
+
+    _DEPENDENCIES: ClassVar[Set[str]] = {"numpy", "nilearn", "scipy"}
+
+    def __del__(self) -> None:
+        """Terminate the class."""
+        # Clear the computation cache
+        logger.debug("Clearing cache for ALFF computation via junifer")
+        self.compute.cache_clear()
+
+    @lru_cache(maxsize=None, typed=True)
+    def compute(
+        self,
+        data: "Nifti1Image",
+        highpass: float,
+        lowpass: float,
+        tr: Optional[float],
+    ) -> Tuple["Nifti1Image", "Nifti1Image", Path, Path]:
+        """Compute ALFF + fALFF map.
+
+        Parameters
+        ----------
+        data : 4D Niimg-like object
+            Images to process.
+        highpass : positive float
+            Highpass cutoff frequency.
+        lowpass : positive float
+            Lowpass cutoff frequency.
+        tr : positive float, optional
+            The Repetition Time of the BOLD data.
+
+        Returns
+        -------
+        Niimg-like object
+            ALFF map.
+        Niimg-like object
+            fALFF map.
+        pathlib.Path
+            The path to the ALFF map as NIfTI.
+        pathlib.Path
+            The path to the fALFF map as NIfTI.
+
+        """
+        logger.debug("Creating cache for ALFF computation via junifer")
+
+        # Get scan data
+        niimg_data = data.get_fdata().copy()
+        if tr is None:
+            tr = float(data.header["pixdim"][4])  # type: ignore
+            logger.info(f"`tr` not provided, using `tr` from header: {tr}")
+
+        # Bandpass the data within the lowpass and highpass cutoff freqs
+        fft_data = sp.fft.fft(niimg_data, axis=-1)
+        fft_freqs = np.abs(sp.fft.fftfreq(niimg_data.shape[-1], tr))
+        # Frequency difference
+        fft_freqs_diff = fft_freqs[1] - fft_freqs[0]
+        # Nyquist frequency
+        nyquist = np.max(fft_freqs)
+        # FFT sample frequency count
+        n_fft = len(fft_freqs)
+        logger.info(
+            f"FFT: nfft = {n_fft}, dFreq = {fft_freqs_diff}, "
+            f"nyquist = {nyquist}"
+        )
+
+        # Compute the denominator on the broadband signal
+        all_freq_mask = fft_freqs > 0
+        denominator = np.sum(np.abs(fft_data[..., all_freq_mask]), axis=-1)
+
+        # Compute the numerator on the bandpassed signal
+        freq_mask = np.logical_and(fft_freqs > highpass, fft_freqs < lowpass)
+
+        # Compute ALFF
+        numerator = np.sum(np.abs(fft_data[..., freq_mask]), axis=-1)
+
+        # Compute fALFF, but avoid division by zero
+        denom_mask = denominator <= 0.000001
+        denominator[denom_mask] = 1  # set to 1 to avoid division by zero
+        # Calculate fALFF
+        falff = np.divide(numerator, denominator)
+        # Set the values where denominator is zero to zero
+        falff[denom_mask] = 0
+
+        # Calculate ALFF
+        alff = numerator / np.sqrt(niimg_data.shape[-1])
+        alff_data = nimg.new_img_like(
+            ref_niimg=data,
+            data=alff,
+        )
+        falff_data = nimg.new_img_like(
+            ref_niimg=data,
+            data=falff,
+        )
+
+        # Create element-scoped tempdir so that the ALFF and fALFF maps are
+        # available later as nibabel stores file path reference for
+        # loading on computation
+        element_tempdir = WorkDirManager().get_element_tempdir(
+            prefix="junifer_alff+falff"
+        )
+        output_alff_path = element_tempdir / "output_alff.nii.gz"
+        output_falff_path = element_tempdir / "output_falff.nii.gz"
+        # Save computed data to file
+        nib.save(alff_data, output_alff_path)
+        nib.save(falff_data, output_falff_path)
+
+        return alff_data, falff_data, output_alff_path, output_falff_path  # type: ignore

junifer/markers/falff/falff_base.py

@@ -1,16 +1,32 @@
-"""Provide abstract class for computing fALFF."""
+"""Provide base class for ALFF / fALFF."""
 
 # Authors: Federico Raimondo <f.raimondo@fz-juelich.de>
 #          Amir Omidvarnia <a.omidvarnia@fz-juelich.de>
 #          Kaustubh R. Patil <k.patil@fz-juelich.de>
+#          Synchon Mandal <s.mandal@fz-juelich.de>
 # License: AGPL
 
-from abc import abstractmethod
+from pathlib import Path
-from typing import ClassVar, Dict, List, Optional, Union
+from typing import (
+    TYPE_CHECKING,
+    Any,
+    ClassVar,
+    Dict,
+    List,
+    Optional,
+    Tuple,
+    Type,
+    Union,
+)
 
-from ...utils.logging import raise_error
+from ...utils.logging import logger, raise_error
 from ..base import BaseMarker
-from .falff_estimator import ALFFEstimator
+from ._afni_falff import AFNIALFF
+from ._junifer_falff import JuniferALFF
+
+
+if TYPE_CHECKING:
+    from nibabel import Nifti1Image
 
 
 class ALFFBase(BaseMarker):
@@ -24,33 +40,45 @@ class ALFFBase(BaseMarker):
         Highpass cutoff frequency.
     lowpass : positive float
         Lowpass cutoff frequency.
+    using : {"junifer", "afni"}
+        Implementation to use for computing ALFF:
+
+        * "junifer" : Use ``junifer``'s own ALFF implementation
+        * "afni" : Use AFNI's ``3dRSFC``
+
     tr : positive float, optional
         The Repetition Time of the BOLD data. If None, will extract
-        the TR from NIFTI header (default None).
+        the TR from NIfTI header (default None).
-    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).
 
     Notes
     -----
-        The `tr` parameter is crucial for the correctness of fALFF/ALFF
+    The ``tr`` parameter is crucial for the correctness of fALFF/ALFF
-        computation. If a dataset is correctly preprocessed, the TR should be
+    computation. If a dataset is correctly preprocessed, the ``tr`` should be
-        extracted from the NIFTI without any issue. However, it has been
+    extracted from the NIfTI without any issue. However, it has been
-        reported that some preprocessed data might not have the correct TR in
+    reported that some preprocessed data might not have the correct ``tr`` in
-        the NIFTI header.
+    the NIfTI header.
+
+    Raises
+    ------
+    ValueError
+        If ``highpass`` is not positive or zero or
+        if ``lowpass`` is not positive or
+        if ``highpass`` is higher than ``lowpass`` or
+        if ``using`` is invalid.
 
     """
 
-    _EXT_DEPENDENCIES: ClassVar[
+    _CONDITIONAL_DEPENDENCIES: ClassVar[List[Dict[str, Union[str, Type]]]] = [
-        List[Dict[str, Union[str, bool, List[str]]]]
-    ] = [
         {
-            "name": "afni",
+            "using": "afni",
-            "optional": True,
+            "depends_on": AFNIALFF,
-            "commands": ["3dRSFC", "3dAFNItoNIFTI"],
+        },
+        {
+            "using": "junifer",
+            "depends_on": JuniferALFF,
         },
     ]
 
@@ -59,8 +87,8 @@ class ALFFBase(BaseMarker):
         fractional: bool,
         highpass: float,
         lowpass: float,
+        using: str,
         tr: Optional[float] = None,
-        use_afni: Optional[bool] = None,
         name: Optional[str] = None,
     ) -> None:
         if highpass < 0:
@@ -71,8 +99,14 @@ class ALFFBase(BaseMarker):
             raise_error("Highpass must be lower than lowpass")
         self.highpass = highpass
         self.lowpass = lowpass
+        # Validate `using` parameter
+        valid_using = [dep["using"] for dep in self._CONDITIONAL_DEPENDENCIES]
+        if using not in valid_using:
+            raise_error(
+                f"Invalid value for `using`, should be one of: {valid_using}"
+            )
+        self.using = using
         self.tr = tr
-        self.use_afni = use_afni
         self.fractional = fractional
 
         # Create a name based on the class name if none is provided
@@ -108,84 +142,52 @@ class ALFFBase(BaseMarker):
         """
         return "vector"
 
-    def compute(
+    def _compute(
         self,
-        input: Dict[str, Dict],
+        input_data: Dict[str, Any],
-        extra_input: Optional[Dict] = None,
+    ) -> Tuple["Nifti1Image", Path]:
-    ) -> Dict:
+        """Compute ALFF and fALFF.
-        """Compute.
 
         Parameters
         ----------
-        input : dict
+        input_data : dict
-            A single input from the pipeline data object in which to compute
+            The input to the marker.
-            the marker.
         extra_input : dict, optional
-            The other fields in the pipeline data object. Useful for accessing
+            The other fields in the pipeline data object (default None).
-            other data kind that needs to be used in the computation. For
-            example, the functional connectivity markers can make use of the
-            confounds if available (default None).
 
         Returns
         -------
-        dict
+        Niimg-like object
-            The computed result as dictionary. This will be either returned
+            The ALFF / fALFF as NIfTI.
-            to the user or stored in the storage by calling the store method
+        pathlib.Path
-            with this as a parameter. The dictionary has the following keys:
+            The path to the ALFF / fALFF as NIfTI.
-
-            * ``data`` : the actual computed values as a numpy.ndarray
-            * ``col_names`` : the column labels for the computed values as list
 
         """
-        if self.use_afni is None:
+        logger.debug("Calculating ALFF and fALFF")
-            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."
-            )
 
-        estimator = ALFFEstimator()
+        # Conditional estimator
-
+        if self.using == "afni":
-        # If the input data space is "native", then alff_path and falff_path
+            estimator = AFNIALFF()
-        # both point to the input data path as it might be required to use
+        elif self.using == "junifer":
-        # in get_corrdinates() for transforming coordinates to native space.
+            estimator = JuniferALFF()
-        alff, falff, alff_path, falff_path = estimator.fit_transform(
+        # Compute ALFF + fALFF
-            use_afni=self.use_afni,
+        alff, falff, alff_path, falff_path = estimator.compute(  # type: ignore
-            input_data=input,
+            data=input_data["data"],
             highpass=self.highpass,
             lowpass=self.lowpass,
             tr=self.tr,
         )
-        post_data = falff if self.fractional else alff
-        post_path = falff_path if self.fractional else alff_path
 
-        post_input = dict(input.items())
+        # If the input data space is native already, the original path should
-        post_input["data"] = post_data
+        # be propagated down as it might be required for transforming
-        post_input["path"] = post_path
+        # parcellation / coordinates to native space, else the
-
+        # path should be passed for use later if required.
-        out = self._postprocess(post_input, extra_input=extra_input)
+        # TODO(synchon): will be taken care in #292
-
+        if input_data["space"] == "native" and self.fractional:
-        return out
+            return falff, input_data["path"]
-
+        elif input_data["space"] == "native" and not self.fractional:
-    @abstractmethod
+            return alff, input_data["path"]
-    def _postprocess(
+        elif input_data["space"] != "native" and self.fractional:
-        self, input: Dict, extra_input: Optional[Dict] = None
+            return falff, falff_path
-    ) -> Dict:
+        else:
-        """Postprocess the output of the estimator.
+            return alff, alff_path
-
-        Parameters
-        ----------
-        input : dict
-            The output of the estimator. It must have the following
-        extra_input : dict, optional
-            The other fields in the pipeline data object. Useful for accessing
-            other data kind that needs to be used in the computation. For
-            example, the functional connectivity markers can make use of the
-            confounds if available (default None).
-
-        """
-        raise_error(
-            "_postprocess must be implemented", klass=NotImplementedError
-        )

junifer/markers/falff/falff_estimator.py

@@ -1,349 +0,0 @@
-"""Provide estimator class for (f)ALFF."""
-
-# Authors: Synchon Mandal <s.mandal@fz-juelich.de>
-#          Federico Raimondo <f.raimondo@fz-juelich.de>
-# License: AGPL
-
-import typing
-from functools import lru_cache
-from pathlib import Path
-from typing import TYPE_CHECKING, Any, Dict, Optional, Tuple, Union
-
-import nibabel as nib
-import numpy as np
-from nilearn import image as nimg
-from scipy.fft import fft, fftfreq
-
-from ...pipeline import WorkDirManager
-from ...pipeline.singleton import singleton
-from ...utils import logger, run_ext_cmd
-
-
-if TYPE_CHECKING:
-    from nibabel import Nifti1Image, Nifti2Image
-
-
-@singleton
-class ALFFEstimator:
-    """Estimator class for (fractional) Amplitude Low Frequency Fluctuation.
-
-    This class is a singleton and is used for efficient computation of fALFF,
-    by caching the voxel-wise ALFF map for a given set of file path and
-    computation parameters.
-
-    .. warning:: This class can only be used via :class:`.ALFFBase` 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 3dRSFC
-        self.temp_dir_path = None
-
-    def __del__(self) -> None:
-        """Cleanup."""
-        # Delete temporary directory and ignore errors for read-only files
-        if self.temp_dir_path is not None:
-            WorkDirManager().delete_tempdir(self.temp_dir_path)
-
-    def _compute_alff_afni(
-        self,
-        data: Union["Nifti1Image", "Nifti2Image"],
-        highpass: float,
-        lowpass: float,
-        tr: Optional[float],
-    ) -> Tuple["Nifti1Image", "Nifti1Image", Path, Path]:
-        """Compute ALFF map via afni's commands.
-
-        Parameters
-        ----------
-        data : 4D Niimg-like object
-            Images to process.
-        highpass : positive float
-            Highpass cutoff frequency.
-        lowpass : positive float
-            Lowpass cutoff frequency.
-        tr : positive float, optional
-            The Repetition Time of the BOLD data.
-
-        Returns
-        -------
-        Niimg-like object
-            ALFF map.
-        Niimg-like object
-            fALFF map.
-        pathlib.Path
-            The path to the ALFF map as NIfTI.
-        pathlib.Path
-            The path to the fALFF map as NIfTI.
-
-        """
-        # Note: self.temp_dir_path is sure to exist before proceeding, so
-        #       types checks are ignored further on.
-
-        # Save niimg to nii.gz
-        nifti_in_file_path = self.temp_dir_path / "input.nii"  # type: ignore
-        nib.save(data, nifti_in_file_path)
-
-        params_suffix = f"_{highpass}_{lowpass}_{tr}"
-        # Create element-scoped tempdir so that the ALFF and fALFF maps are
-        # available later as get_coordinates and the like need it
-        # in fALFFSpheres and the like to transform to other template
-        # spaces
-        element_tempdir = WorkDirManager().get_element_tempdir(
-            prefix="alff_falff_afni"
-        )
-        alff_fname = element_tempdir / f"alff{params_suffix}.nii"
-        falff_fname = element_tempdir / f"falff{params_suffix}.nii"
-
-        # Use afni's 3dRSFC to compute ALFF and fALFF
-        falff_afni_out_path_prefix = (
-            self.temp_dir_path / "temp_falff"  # type: ignore
-        )
-
-        # Set 3dRSFC command
-        bp_cmd = [
-            "3dRSFC",
-            f"-prefix {falff_afni_out_path_prefix.resolve()}",
-            f"-input {nifti_in_file_path.resolve()}",
-            f"-band {highpass} {lowpass}",
-            "-no_rsfa -nosat -nodetrend",
-        ]
-        if tr is not None:
-            bp_cmd.append(f"-dt {tr}")
-        # Call 3dRSFC
-        run_ext_cmd(name="3dRSFC", cmd=bp_cmd)
-
-        # Convert alff output to nifti
-        convert_alff_cmd = [
-            "3dAFNItoNIFTI",
-            f"-prefix {alff_fname.resolve()}",
-            f"{falff_afni_out_path_prefix}_ALFF+tlrc.BRIK",
-        ]
-        # Call 3dAFNItoNIFTI
-        run_ext_cmd(name="3dAFNItoNIFTI", cmd=convert_alff_cmd)
-
-        # Convert falff output to nifti
-        convert_falff_cmd = [
-            "3dAFNItoNIFTI",
-            f"-prefix {falff_fname.resolve()}",
-            f"{falff_afni_out_path_prefix}_fALFF+tlrc.BRIK",
-        ]
-        # Call 3dAFNItoNIFTI
-        run_ext_cmd(name="3dAFNItoNIFTI", cmd=convert_falff_cmd)
-
-        # Cleanup intermediate files
-        for fname in self.temp_dir_path.glob("temp_*"):  # type: ignore
-            fname.unlink()
-
-        # Load niftis
-        alff_img = nib.load(alff_fname)
-        falff_img = nib.load(falff_fname)
-
-        return alff_img, falff_img, alff_fname, falff_fname  # type: ignore
-
-    def _compute_alff_python(
-        self,
-        data: Union["Nifti1Image", "Nifti2Image"],
-        highpass: float,
-        lowpass: float,
-        tr: Optional[float],
-    ) -> Tuple["Nifti1Image", "Nifti1Image", Path, Path]:
-        """Compute (f)ALFF map.
-
-        Parameters
-        ----------
-        data : 4D Niimg-like object
-            Images to process.
-        highpass : positive float
-            Highpass cutoff frequency.
-        lowpass : positive float
-            Lowpass cutoff frequency.
-        tr : positive float, optional
-            The Repetition Time of the BOLD data.
-
-        Returns
-        -------
-        Niimg-like object
-            ALFF map.
-        Niimg-like object
-            fALFF map.
-        pathlib.Path
-            The path to the ALFF map as NIfTI.
-        pathlib.Path
-            The path to the fALFF map as NIfTI.
-
-        """
-        timeseries = data.get_fdata().copy()
-        if tr is None:
-            tr = float(data.header["pixdim"][4])  # type: ignore
-            logger.info(f"TR Not provided, using TR from header = {tr}")
-        # bandpass the data within the lowpass and highpass cutoff freqs
-
-        ts_fft = fft(timeseries, axis=-1)
-        ts_fft = typing.cast(np.ndarray, ts_fft)
-        fft_freqs = np.abs(fftfreq(timeseries.shape[-1], tr))
-
-        dFreq = fft_freqs[1] - fft_freqs[0]
-        nyquist = np.max(fft_freqs)
-        nfft = len(fft_freqs)
-        logger.info(
-            f"FFT: nfft = {nfft}, dFreq = {dFreq}, nyquist = {nyquist}"
-        )
-
-        # First compute the denominator on the broadband signal
-        all_freq_mask = fft_freqs > 0
-        denominator = np.sum(np.abs(ts_fft[..., all_freq_mask]), axis=-1)
-
-        # Compute the numerator on the bandpassed signal
-        freq_mask = np.logical_and(fft_freqs > highpass, fft_freqs < lowpass)
-        # Compute ALFF
-        numerator = np.sum(np.abs(ts_fft[..., freq_mask]), axis=-1)
-
-        # Compute fALFF, but avoid division by zero
-        denom_mask = denominator <= 0.000001
-        denominator[denom_mask] = 1  # set to 1 to avoid division by zero
-        python_falff = np.divide(numerator, denominator)
-        # Set the values where denominator is zero to zero
-        python_falff[denom_mask] = 0
-
-        python_alff = numerator / np.sqrt(timeseries.shape[-1])
-        alff_img = nimg.new_img_like(data, python_alff)
-        falff_img = nimg.new_img_like(data, python_falff)
-        # Create element-scoped tempdir so that the ALFF and fALFF maps are
-        # available later as get_coordinates and the like need it
-        # in fALFFSpheres and the like to transform to other template
-        # spaces
-        element_tempdir = WorkDirManager().get_element_tempdir(
-            prefix="alff_falff_python"
-        )
-        alff_output_path = element_tempdir / "alff_map_python.nii.gz"
-        falff_output_path = element_tempdir / "falff_map_python.nii.gz"
-        nib.save(alff_img, alff_output_path)
-        nib.save(falff_img, falff_output_path)
-        return alff_img, falff_img, alff_output_path, falff_output_path  # type: ignore
-
-    @lru_cache(maxsize=None, typed=True)
-    def _compute(
-        self,
-        use_afni: bool,
-        data: Union["Nifti1Image", "Nifti2Image"],
-        highpass: float,
-        lowpass: float,
-        tr: Optional[float],
-    ) -> Tuple["Nifti1Image", "Nifti1Image", Path, Path]:
-        """Compute the ALFF map with memorization.
-
-        Parameters
-        ----------
-        use_afni : bool
-            Whether to use AFNI for computing.
-        data : 4D Niimg-like object
-            Images to process.
-        highpass : positive float
-            Highpass cutoff frequency.
-        lowpass : positive float
-            Lowpass cutoff frequency.
-        tr : positive float, optional
-            The Repetition Time of the BOLD data.
-
-        Returns
-        -------
-        Niimg-like object
-            ALFF map.
-        Niimg-like object
-            fALFF map.
-        pathlib.Path
-            The path to the ALFF map as NIfTI.
-        pathlib.Path
-            The path to the fALFF map as NIfTI.
-
-        """
-        if use_afni:
-            # Create new temporary directory before using AFNI
-            self.temp_dir_path = WorkDirManager().get_tempdir(prefix="falff")
-            return self._compute_alff_afni(
-                data=data,
-                highpass=highpass,
-                lowpass=lowpass,
-                tr=tr,
-            )
-
-        return self._compute_alff_python(
-            data, highpass=highpass, lowpass=lowpass, tr=tr
-        )
-
-    def fit_transform(
-        self,
-        use_afni: bool,
-        input_data: Dict[str, Any],
-        highpass: float,
-        lowpass: float,
-        tr: Optional[float],
-    ) -> Tuple["Nifti1Image", "Nifti1Image", Path, Path]:
-        """Fit and transform for the estimator.
-
-        Parameters
-        ----------
-        use_afni : bool
-            Whether to use AFNI for computing.
-        input_data : dict
-            The BOLD data as dictionary.
-        highpass : positive float
-            Highpass cutoff frequency.
-        lowpass : positive float
-            Lowpass cutoff frequency.
-        tr : positive float, optional
-            The Repetition Time of the BOLD data.
-
-        Returns
-        -------
-        Niimg-like object
-            ALFF map.
-        Niimg-like object
-            fALFF map.
-        pathlib.Path
-            The path to the ALFF map as NIfTI or the input data path if the
-            input data space is "native".
-        pathlib.Path
-            The path to the fALFF map as NIfTI or the input data path if the
-            input data space is "native".
-
-        """
-        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 fALFF map cache for {self._file_path}.")
-            # Clear the cache
-            self._compute.cache_clear()
-            # Clear temporary directory files
-            if self.temp_dir_path is not None:
-                WorkDirManager().delete_tempdir(self.temp_dir_path)
-            # Set the new file path
-            self._file_path = bold_path
-        else:
-            logger.info(f"Using fALFF map cache for {self._file_path}.")
-        # Compute
-        alff_map, falff_map, alff_map_path, falff_map_path = self._compute(
-            use_afni=use_afni,
-            data=bold_data,
-            highpass=highpass,
-            lowpass=lowpass,
-            tr=tr,
-        )
-        # If the input data space is native already, the original path should
-        # be propagated down as it might be required for transforming
-        # coordinates to native space via get_coordinates(), else the alff
-        # / falff map path should be passed for use later if required.
-        if input_data["space"] == "native":
-            return alff_map, falff_map, input_data["path"], input_data["path"]
-
-        return alff_map, falff_map, alff_map_path, falff_map_path

junifer/markers/falff/falff_parcels.py

@@ -1,20 +1,22 @@
-"""Provide class for computing fALFF on parcels."""
+"""Provide class for ALFF / fALFF on parcels."""
 
 # Authors: Federico Raimondo <f.raimondo@fz-juelich.de>
 #          Amir Omidvarnia <a.omidvarnia@fz-juelich.de>
 #          Kaustubh R. Patil <k.patil@fz-juelich.de>
+#          Synchon Mandal <s.mandal@fz-juelich.de>
 # License: AGPL
 
-from typing import Dict, List, Optional, Union
+from typing import Any, Dict, List, Optional, Union
 
 from ...api.decorators import register_marker
-from .. import ParcelAggregation
+from ...utils import logger
+from ..parcel_aggregation import ParcelAggregation
 from .falff_base import ALFFBase
 
 
 @register_marker
 class ALFFParcels(ALFFBase):
-    """Class for computing fALFF/ALFF on parcels.
+    """Class for ALFF / fALFF on parcels.
 
     Parameters
     ----------
@@ -23,6 +25,12 @@ class ALFFParcels(ALFFBase):
         :func:`.list_parcellations`.
     fractional : bool
         Whether to compute fractional ALFF.
+    using : {"junifer", "afni"}
+        Implementation to use for computing ALFF:
+
+        * "junifer" : Use ``junifer``'s own ALFF implementation
+        * "afni" : Use AFNI's ``3dRSFC``
+
     highpass : positive float, optional
         The highpass cutoff frequency for the bandpass filter. If 0,
         it will not apply a highpass filter (default 0.01).
@@ -30,20 +38,17 @@ class ALFFParcels(ALFFBase):
         The lowpass cutoff frequency for the bandpass filter (default 0.1).
     tr : positive float, optional
         The Repetition Time of the BOLD data. If None, will extract
-        the TR from NIFTI header (default None).
+        the TR from NIfTI header (default None).
-    use_afni : bool, optional
+    agg_method : str, optional
-        Whether to use AFNI for computing. If None, will use AFNI only
+        The method to perform aggregation using. Check valid options in
-        if available (default None).
+        :func:`.get_aggfunc_by_name` (default "mean").
+    agg_method_params : dict, optional
+        Parameters to pass to the aggregation function. Check valid options in
+        :func:`.get_aggfunc_by_name` (default None).
     masks : str, dict or list of dict or str, optional
         The specification of the masks to apply to regions before extracting
         signals. Check :ref:`Using Masks <using_masks>` for more details.
         If None, will not apply any mask (default None).
-    method : str, optional
-        The method to perform aggregation using. Check valid options in
-        :func:`.get_aggfunc_by_name` (default "mean").
-    method_params : dict, optional
-        Parameters to pass to the aggregation function. Check valid options in
-        :func:`.get_aggfunc_by_name`.
     name : str, optional
         The name of the marker. If None, will use the class name (default
         None).
@@ -51,10 +56,10 @@ class ALFFParcels(ALFFBase):
     Notes
     -----
     The ``tr`` parameter is crucial for the correctness of fALFF/ALFF
-    computation. If a dataset is correctly preprocessed, the TR should be
+    computation. If a dataset is correctly preprocessed, the ``tr`` should be
-    extracted from the NIFTI without any issue. However, it has been
+    extracted from the NIfTI without any issue. However, it has been
-    reported that some preprocessed data might not have the correct TR in
+    reported that some preprocessed data might not have the correct ``tr`` in
-    the NIFTI header.
+    the NIfTI header.
 
     ALFF/fALFF are computed using a bandpass butterworth filter. See
     :func:`scipy.signal.butter` and :func:`scipy.signal.filtfilt` for more
@@ -66,63 +71,73 @@ class ALFFParcels(ALFFBase):
         self,
         parcellation: Union[str, List[str]],
         fractional: bool,
+        using: str,
         highpass: float = 0.01,
         lowpass: float = 0.1,
         tr: Optional[float] = None,
-        use_afni: Optional[bool] = None,
+        agg_method: str = "mean",
+        agg_method_params: Optional[Dict] = None,
         masks: Union[str, Dict, List[Union[Dict, str]], None] = None,
-        method: str = "mean",
-        method_params: Optional[Dict] = None,
         name: Optional[str] = None,
     ) -> None:
-        self.parcellation = parcellation
+        # Superclass init first to validate `using` parameter
-        self.masks = masks
-        self.method = method
-        self.method_params = method_params
         super().__init__(
             fractional=fractional,
             highpass=highpass,
             lowpass=lowpass,
+            using=using,
             tr=tr,
             name=name,
-            use_afni=use_afni,
         )
+        self.parcellation = parcellation
+        self.agg_method = agg_method
+        self.agg_method_params = agg_method_params
+        self.masks = masks
 
-    def _postprocess(
+    def compute(
-        self, input: Dict, extra_input: Optional[Dict] = None
+        self,
-    ) -> Dict:
+        input: Dict[str, Any],
-        """Compute ALFF and fALFF.
+        extra_input: Optional[Dict[str, Any]] = None,
+    ) -> Dict[str, Any]:
+        """Compute.
 
         Parameters
         ----------
         input : dict
-            A single input from the pipeline data object in which to compute
+            The BOLD data as dictionary.
-            the marker.
         extra_input : dict, optional
-            The other fields in the pipeline data object. Useful for accessing
+            The other fields in the pipeline data object (default None).
-            other data kind that needs to be used in the computation. For
-            example, the functional connectivity markers can make use of the
-            confounds if available (default None).
 
         Returns
         -------
         dict
-            The computed ALFF as dictionary. The dictionary has the following
+            The computed result as dictionary. The dictionary has the following
             keys:
 
             * ``data`` : the actual computed values as a numpy.ndarray
             * ``col_names`` : the column labels for the computed values as list
 
         """
-        pa = ParcelAggregation(
+        logger.info("Calculating ALFF / fALFF for parcels")
+
+        # Compute ALFF / fALFF
+        output_data, output_file_path = self._compute(input_data=input)
+
+        # Initialize parcel aggregation
+        parcel_aggregation = ParcelAggregation(
             parcellation=self.parcellation,
-            method=self.method,
+            method=self.agg_method,
-            method_params=self.method_params,
+            method_params=self.agg_method_params,
             masks=self.masks,
-            on="fALFF",
+            on="BOLD",
+        )
+        # Perform aggregation on ALFF / fALFF
+        parcel_aggregation_input = dict(input.items())
+        parcel_aggregation_input["data"] = output_data
+        parcel_aggregation_input["path"] = output_file_path
+        output = parcel_aggregation.compute(
+            input=parcel_aggregation_input,
+            extra_input=extra_input,
         )
 
-        # get the 2D timeseries after parcel aggregation
+        return output
-        out = pa.compute(input, extra_input=extra_input)
-
-        return out

junifer/markers/falff/falff_spheres.py

@@ -1,26 +1,36 @@
-"""Provide class for computing fALFF on spheres."""
+"""Provide class for ALFF / fALFF on spheres."""
 
 # Authors: Federico Raimondo <f.raimondo@fz-juelich.de>
 #          Amir Omidvarnia <a.omidvarnia@fz-juelich.de>
 #          Kaustubh R. Patil <k.patil@fz-juelich.de>
+#          Synchon Mandal <s.mandal@fz-juelich.de>
 # License: AGPL
 
-from typing import Dict, List, Optional, Union
+from typing import Any, Dict, List, Optional, Union
 
 from ...api.decorators import register_marker
-from .. import SphereAggregation
+from ...utils import logger
+from ..sphere_aggregation import SphereAggregation
 from .falff_base import ALFFBase
 
 
 @register_marker
 class ALFFSpheres(ALFFBase):
-    """Class for computing fALFF/ALFF on spheres.
+    """Class for computing ALFF / fALFF on spheres.
 
     Parameters
     ----------
     coords : str
         The name of the coordinates list to use. See
         :func:`.list_coordinates` for options.
+    fractional : bool
+        Whether to compute fractional ALFF.
+    using : {"junifer", "afni"}
+        Implementation to use for computing ALFF:
+
+        * "junifer" : Use ``junifer``'s own ALFF implementation
+        * "afni" : Use AFNI's ``3dRSFC``
+
     radius : float, optional
         The radius of the sphere in mm. If None, the signal will be extracted
         from a single voxel. See :class:`nilearn.maskers.NiftiSpheresMasker`
@@ -28,8 +38,6 @@ class ALFFSpheres(ALFFBase):
     allow_overlap : bool, optional
         Whether to allow overlapping spheres. If False, an error is raised if
         the spheres overlap (default is False).
-    fractional : bool
-        Whether to compute fractional ALFF.
     highpass : positive float, optional
         The highpass cutoff frequency for the bandpass filter. If 0,
         it will not apply a highpass filter (default 0.01).
@@ -37,20 +45,17 @@ class ALFFSpheres(ALFFBase):
         The lowpass cutoff frequency for the bandpass filter (default 0.1).
     tr : positive float, optional
         The Repetition Time of the BOLD data. If None, will extract
-        the TR from NIFTI header (default None).
+        the TR from NIfTI header (default None).
-    use_afni : bool, optional
+    agg_method : str, optional
-        Whether to use AFNI for computing. If None, will use AFNI only
+        The method to perform aggregation using. Check valid options in
-        if available (default None).
+        :func:`.get_aggfunc_by_name` (default "mean").
+    agg_method_params : dict, optional
+        Parameters to pass to the aggregation function. Check valid options in
+        :func:`.get_aggfunc_by_name`.
     masks : str, dict or list of dict or str, optional
         The specification of the masks to apply to regions before extracting
         signals. Check :ref:`Using Masks <using_masks>` for more details.
         If None, will not apply any mask (default None).
-    method : str, optional
-        The method to perform aggregation using. Check valid options in
-        :func:`.get_aggfunc_by_name` (default "mean").
-    method_params : dict, optional
-        Parameters to pass to the aggregation function. Check valid options in
-        :func:`.get_aggfunc_by_name`.
     name : str, optional
         The name of the marker. If None, will use the class name (default
         None).
@@ -58,9 +63,9 @@ class ALFFSpheres(ALFFBase):
     Notes
     -----
     The ``tr`` parameter is crucial for the correctness of fALFF/ALFF
-    computation. If a dataset is correctly preprocessed, the TR should be
+    computation. If a dataset is correctly preprocessed, the ``tr`` should be
-    extracted from the NIFTI without any issue. However, it has been
+    extracted from the NIfTI without any issue. However, it has been
-    reported that some preprocessed data might not have the correct TR in
+    reported that some preprocessed data might not have the correct ``tr`` in
     the NIFTI header.
 
     ALFF/fALFF are computed using a bandpass butterworth filter. See
@@ -73,69 +78,79 @@ class ALFFSpheres(ALFFBase):
         self,
         coords: str,
         fractional: bool,
+        using: str,
         radius: Optional[float] = None,
         allow_overlap: bool = False,
         highpass: float = 0.01,
         lowpass: float = 0.1,
         tr: Optional[float] = None,
-        use_afni: Optional[bool] = None,
+        agg_method: str = "mean",
+        agg_method_params: Optional[Dict] = None,
         masks: Union[str, Dict, List[Union[Dict, str]], None] = None,
-        method: str = "mean",
-        method_params: Optional[Dict] = None,
         name: Optional[str] = None,
     ) -> None:
-        self.coords = coords
+        # Superclass init first to validate `using` parameter
-        self.radius = radius
-        self.allow_overlap = allow_overlap
-        self.masks = masks
-        self.method = method
-        self.method_params = method_params
         super().__init__(
             fractional=fractional,
             highpass=highpass,
             lowpass=lowpass,
+            using=using,
             tr=tr,
             name=name,
-            use_afni=use_afni,
         )
+        self.coords = coords
+        self.radius = radius
+        self.allow_overlap = allow_overlap
+        self.agg_method = agg_method
+        self.agg_method_params = agg_method_params
+        self.masks = masks
 
-    def _postprocess(
+    def compute(
-        self, input: Dict, extra_input: Optional[Dict] = None
+        self,
-    ) -> Dict:
+        input: Dict[str, Any],
-        """Compute ALFF and fALFF.
+        extra_input: Optional[Dict[str, Any]] = None,
+    ) -> Dict[str, Any]:
+        """Compute.
 
         Parameters
         ----------
         input : dict
-            A single input from the pipeline data object in which to compute
+            The BOLD data as dictionary.
-            the marker.
         extra_input : dict, optional
-            The other fields in the pipeline data object. Useful for accessing
+            The other fields in the pipeline data object (default None).
-            other data kind that needs to be used in the computation. For
-            example, the functional connectivity markers can make use of the
-            confounds if available (default None).
 
         Returns
         -------
         dict
-            The computed ALFF as dictionary. The dictionary has the following
+            The computed result as dictionary. The dictionary has the following
             keys:
 
             * ``data`` : the actual computed values as a numpy.ndarray
             * ``col_names`` : the column labels for the computed values as list
 
         """
-        pa = SphereAggregation(
+        logger.info("Calculating ALFF / fALFF for spheres")
+
+        # Compute ALFF / fALFF
+        output_data, output_file_path = self._compute(input_data=input)
+
+        # Initialize sphere aggregation
+        sphere_aggregation = SphereAggregation(
             coords=self.coords,
             radius=self.radius,
             allow_overlap=self.allow_overlap,
-            method=self.method,
+            method=self.agg_method,
-            method_params=self.method_params,
+            method_params=self.agg_method_params,
             masks=self.masks,
-            on="fALFF",
+            on="BOLD",
+        )
+        # Perform aggregation on ALFF / fALFF
+        sphere_aggregation_input = dict(input.items())
+        sphere_aggregation_input["data"] = output_data
+        sphere_aggregation_input["path"] = output_file_path
+        output = sphere_aggregation.compute(
+            input=sphere_aggregation_input,
+            extra_input=extra_input,
         )
 
-        # get the 2D timeseries after sphere aggregation
+        return output
-        out = pa.compute(input, extra_input=extra_input)
-
-        return out

junifer/markers/falff/tests/test_falff_estimator.py

@@ -1,267 +0,0 @@
-"""Provide test for (f)ALFF estimator."""
-
-# Authors: Federico Raimondo <f.raimondo@fz-juelich.de>
-# License: AGPL
-
-import time
-from pathlib import Path
-
-import pytest
-from nibabel import Nifti1Image
-from scipy.stats import pearsonr
-
-from junifer.datareader import DefaultDataReader
-from junifer.markers.falff.falff_estimator import ALFFEstimator
-from junifer.pipeline import WorkDirManager
-from junifer.pipeline.utils import _check_afni
-from junifer.testing.datagrabbers import PartlyCloudyTestingDataGrabber
-from junifer.utils import logger
-
-
-def test_ALFFEstimator_cache_python(tmp_path: Path) -> None:
-    """Test that the cache works properly when using python.
-
-    Parameters
-    ----------
-    tmp_path : pathlib.Path
-        The path to the test directory.
-
-    """
-    # Get subject from datagrabber
-    with PartlyCloudyTestingDataGrabber() as dg:
-        subject = dg["sub-01"]
-    # Read data for subject
-    subject_data = DefaultDataReader().fit_transform(subject)
-    # Update workdir to current test's tmp_path
-    WorkDirManager().workdir = tmp_path
-    # Setup estimator
-    estimator = ALFFEstimator()
-
-    # Compute without cache
-    start_time = time.time()
-    alff, falff, alff_path, falff_path = estimator.fit_transform(
-        use_afni=False,
-        input_data=subject_data["BOLD"],
-        highpass=0.01,
-        lowpass=0.1,
-        tr=None,
-    )
-    first_time = time.time() - start_time
-    logger.info(f"ALFF Estimator First time: {first_time}")
-    assert isinstance(alff, Nifti1Image)
-    assert isinstance(falff, Nifti1Image)
-    assert isinstance(alff_path, Path)
-    assert isinstance(falff_path, Path)
-
-    # Compute again with cache, should be faster
-    start_time = time.time()
-    alff, falff, _, _ = estimator.fit_transform(
-        use_afni=False,
-        input_data=subject_data["BOLD"],
-        highpass=0.01,
-        lowpass=0.1,
-        tr=None,
-    )
-    second_time = time.time() - start_time
-    logger.info(f"ALFF Estimator Second time: {second_time}")
-    assert second_time < (first_time / 1000)
-
-    # Change a parameter and compute again without cache
-    start_time = time.time()
-    alff, falff, _, _ = estimator.fit_transform(
-        use_afni=False,
-        input_data=subject_data["BOLD"],
-        highpass=0.01,
-        lowpass=0.11,
-        tr=None,
-    )
-    third_time = time.time() - start_time
-    logger.info(f"ALFF Estimator Third time: {third_time}")
-    assert third_time > (first_time / 10)
-
-    # Compute again with cache, should be faster
-    start_time = time.time()
-    alff, falff, _, _ = estimator.fit_transform(
-        use_afni=False,
-        input_data=subject_data["BOLD"],
-        highpass=0.01,
-        lowpass=0.1,
-        tr=None,
-    )
-    fourth = time.time() - start_time
-    logger.info(f"ALFF Estimator Fourth time: {fourth}")
-    assert fourth < (first_time / 1000)
-
-    # Change the data and it should clear the cache
-    with PartlyCloudyTestingDataGrabber() as dg:
-        subject = dg["sub-02"]
-    # Read data for new subject
-    subject_data = DefaultDataReader().fit_transform(subject)
-
-    start_time = time.time()
-    alff, falff, _, _ = estimator.fit_transform(
-        use_afni=False,
-        input_data=subject_data["BOLD"],
-        highpass=0.01,
-        lowpass=0.1,
-        tr=None,
-    )
-    fifth = time.time() - start_time
-    logger.info(f"ALFF Estimator Fifth time: {fifth}")
-    assert fifth > (first_time / 10)
-
-
-@pytest.mark.skipif(
-    _check_afni() is False, reason="requires afni to be in PATH"
-)
-def test_ALFFEstimator_cache_afni(tmp_path: Path) -> None:
-    """Test that the cache works properly when using afni.
-
-    Parameters
-    ----------
-    tmp_path : pathlib.Path
-        The path to the test directory.
-
-    """
-    # Get subject from datagrabber
-    with PartlyCloudyTestingDataGrabber() as dg:
-        subject = dg["sub-01"]
-    # Read data for subject
-    subject_data = DefaultDataReader().fit_transform(subject)
-    # Update workdir to current test's tmp_path
-    WorkDirManager().workdir = tmp_path
-    # Setup estimator
-    estimator = ALFFEstimator()
-
-    # Compute with cache
-    start_time = time.time()
-    alff, falff, alff_path, falff_path = estimator.fit_transform(
-        use_afni=True,
-        input_data=subject_data["BOLD"],
-        highpass=0.01,
-        lowpass=0.1,
-        tr=None,
-    )
-    first_time = time.time() - start_time
-    logger.info(f"ALFF Estimator First time: {first_time}")
-    assert isinstance(alff, Nifti1Image)
-    assert isinstance(falff, Nifti1Image)
-    assert isinstance(alff_path, Path)
-    assert isinstance(falff_path, Path)
-    n_files = len(list(estimator.temp_dir_path.glob("*")))
-    assert n_files == 1  # only input file
-
-    # Compute again with cache, should be faster
-    start_time = time.time()
-    alff, falff, _, _ = estimator.fit_transform(
-        use_afni=True,
-        input_data=subject_data["BOLD"],
-        highpass=0.01,
-        lowpass=0.1,
-        tr=None,
-    )
-    second_time = time.time() - start_time
-    logger.info(f"ALFF Estimator Second time: {second_time}")
-    assert second_time < (first_time / 1000)
-    n_files = len(list(estimator.temp_dir_path.glob("*")))
-    assert n_files == 1  # only input file
-
-    # Change a parameter and compute again without cache
-    start_time = time.time()
-    alff, falff, _, _ = estimator.fit_transform(
-        use_afni=True,
-        input_data=subject_data["BOLD"],
-        highpass=0.01,
-        lowpass=0.11,
-        tr=None,
-    )
-    third_time = time.time() - start_time
-    logger.info(f"ALFF Estimator Third time: {third_time}")
-    assert third_time > (first_time / 10)
-    n_files = len(list(estimator.temp_dir_path.glob("*")))
-    assert n_files == 1  # only input file
-
-    # Compute with cache, should be faster
-    start_time = time.time()
-    alff, falff, _, _ = estimator.fit_transform(
-        use_afni=True,
-        input_data=subject_data["BOLD"],
-        highpass=0.01,
-        lowpass=0.1,
-        tr=None,
-    )
-    fourth = time.time() - start_time
-    logger.info(f"ALFF Estimator Fourth time: {fourth}")
-    assert fourth < (first_time / 1000)
-    n_files = len(list(estimator.temp_dir_path.glob("*")))
-    assert n_files == 1  # only input file
-
-    # Change the data and it should clear the cache
-    with PartlyCloudyTestingDataGrabber() as dg:
-        subject = dg["sub-02"]
-    # Read data for new subject
-    subject_data = DefaultDataReader().fit_transform(subject)
-
-    start_time = time.time()
-    alff, falff, _, _ = estimator.fit_transform(
-        use_afni=True,
-        input_data=subject_data["BOLD"],
-        highpass=0.01,
-        lowpass=0.1,
-        tr=None,
-    )
-    fifth = time.time() - start_time
-    logger.info(f"ALFF Estimator Fifth time: {fifth}")
-    assert fifth > (first_time / 10)
-    n_files = len(list(estimator.temp_dir_path.glob("*")))
-    assert n_files == 1  # only input file
-
-
-@pytest.mark.skipif(
-    _check_afni() is False, reason="requires afni to be in PATH"
-)
-def test_ALFFEstimator_afni_vs_python(tmp_path: Path) -> None:
-    """Test that the cache works properly when using afni.
-
-    Parameters
-    ----------
-    tmp_path : pathlib.Path
-        The path to the test directory.
-
-    """
-    # Get subject from datagrabber
-    with PartlyCloudyTestingDataGrabber() as dg:
-        subject = dg["sub-01"]
-    # Read data for subject
-    subject_data = DefaultDataReader().fit_transform(subject)
-    # Update workdir to current test's tmp_path
-    WorkDirManager().workdir = tmp_path
-    # Setup estimator
-    estimator = ALFFEstimator()
-
-    # Use an arbitrary TR to test the AFNI vs Python implementation
-    afni_alff, afni_falff, _, _ = estimator.fit_transform(
-        use_afni=True,
-        input_data=subject_data["BOLD"],
-        highpass=0.01,
-        lowpass=0.1,
-        tr=2.5,
-    )
-
-    python_alff, python_falff, _, _ = estimator.fit_transform(
-        use_afni=False,
-        input_data=subject_data["BOLD"],
-        highpass=0.01,
-        lowpass=0.1,
-        tr=2.5,
-    )
-
-    r, _ = pearsonr(
-        afni_alff.get_fdata().flatten(), python_alff.get_fdata().flatten()
-    )
-    assert r > 0.99
-
-    r, _ = pearsonr(
-        afni_falff.get_fdata().flatten(), python_falff.get_fdata().flatten()
-    )
-    assert r > 0.99

junifer/markers/falff/tests/test_falff_parcels.py

@@ -1,14 +1,14 @@
-"""Provide test for parcel-aggregated (f)ALFF."""
+"""Provide tests for ALFFParcels."""
 
 # Authors: Federico Raimondo <f.raimondo@fz-juelich.de>
 #          Synchon Mandal <s.mandal@fz-juelich.de>
 # License: AGPL
 
+import logging
 from pathlib import Path
 
 import pytest
-from numpy.testing import assert_array_equal
+import scipy as sp
-from scipy.stats import pearsonr
 
 from junifer.datareader import DefaultDataReader
 from junifer.markers.falff import ALFFParcels
@@ -16,89 +16,72 @@ from junifer.pipeline import WorkDirManager
 from junifer.pipeline.utils import _check_afni
 from junifer.storage import SQLiteFeatureStorage
 from junifer.testing.datagrabbers import PartlyCloudyTestingDataGrabber
-from junifer.utils import logger
 
 
-_PARCELLATION = "Schaefer100x7"
+PARCELLATION = "Schaefer100x7"
 
 
-def test_ALFFParcels_python(tmp_path: Path) -> None:
+def test_ALFFParcels(caplog: pytest.LogCaptureFixture, tmp_path: Path) -> None:
-    """Test ALFFParcels using python.
+    """Test ALFFParcels.
 
     Parameters
     ----------
+    caplog : pytest.LogCaptureFixture
+        The pytest.LogCaptureFixture object.
     tmp_path : pathlib.Path
         The path to the test directory.
 
     """
-    with PartlyCloudyTestingDataGrabber() as dg:
+    with caplog.at_level(logging.DEBUG):
-        input_ = dg["sub-01"]
+        with PartlyCloudyTestingDataGrabber() as dg:
+            element_data = DefaultDataReader().fit_transform(dg["sub-01"])
+            # Update workdir to current test's tmp_path
+            WorkDirManager().workdir = tmp_path
 
-    input_ = DefaultDataReader().fit_transform(input_)
+            # Initialize marker
-    WorkDirManager().workdir = tmp_path
+            marker = ALFFParcels(
-    marker = ALFFParcels(
+                parcellation=PARCELLATION,
-        parcellation=_PARCELLATION,
+                fractional=False,
-        method="mean",
+                using="junifer",
-        use_afni=False,
+            )
-        fractional=False,
+            # Fit transform marker on data
-    )
+            output = marker.fit_transform(element_data)
-    python_values = marker.fit_transform(input_)["BOLD"]["data"]
 
-    assert marker.use_afni is False
+            assert "Creating cache" in caplog.text
-    assert python_values.ndim == 2
+
-    assert python_values.shape == (1, 100)
+            # Get BOLD output
+            assert "BOLD" in output
+            output_bold = output["BOLD"]
+            # Assert BOLD output keys
+            assert "data" in output_bold
+            assert "col_names" in output_bold
+
+            output_bold_data = output_bold["data"]
+            # Assert BOLD output data dimension
+            assert output_bold_data.ndim == 2
+            assert output_bold_data.shape == (1, 100)
+
+            # Reset log capture
+            caplog.clear()
+            # Initialize storage
+            storage = SQLiteFeatureStorage(tmp_path / "falff_parcels.sqlite")
+            # Fit transform marker on data with storage
+            marker.fit_transform(
+                input=element_data,
+                storage=storage,
+            )
+            # Cache working correctly
+            assert "Creating cache" not in caplog.text
 
 
 @pytest.mark.skipif(
-    _check_afni() is False, reason="requires afni to be in PATH"
+    _check_afni() is False, reason="requires AFNI to be in PATH"
-)
-def test_ALFFParcels_afni(tmp_path: Path) -> None:
-    """Test ALFFParcels using afni.
-
-    Parameters
-    ----------
-    tmp_path : pathlib.Path
-        The path to the test directory.
-
-    """
-    with PartlyCloudyTestingDataGrabber() as dg:
-        input_ = dg["sub-01"]
-
-    input_ = DefaultDataReader().fit_transform(input_)
-    WorkDirManager().workdir = tmp_path
-    marker = ALFFParcels(
-        parcellation=_PARCELLATION,
-        method="mean",
-        use_afni=True,
-        fractional=False,
-    )
-    assert marker.use_afni is True
-    afni_values = marker.fit_transform(input_)["BOLD"]["data"]
-
-    assert afni_values.ndim == 2
-    assert afni_values.shape == (1, 100)
-
-    # Again, should be blazing fast
-    marker = ALFFParcels(
-        parcellation=_PARCELLATION, method="mean", fractional=False
-    )
-    assert marker.use_afni is None
-    afni_values2 = marker.fit_transform(input_)["BOLD"]["data"]
-    assert marker.use_afni is True
-    assert_array_equal(afni_values, afni_values2)
-
-
-@pytest.mark.skipif(
-    _check_afni() is False, reason="requires afni to be in PATH"
 )
 @pytest.mark.parametrize(
     "fractional", [True, False], ids=["fractional", "non-fractional"]
 )
-def test_ALFFParcels_python_vs_afni(
+def test_ALFFParcels_comparison(tmp_path: Path, fractional: bool) -> None:
-    tmp_path: Path,
+    """Test ALFFParcels implementation comparison.
-    fractional: bool,
-) -> None:
-    """Test ALFFParcels using python.
 
     Parameters
     ----------
@@ -109,65 +92,35 @@ def test_ALFFParcels_python_vs_afni(
 
     """
     with PartlyCloudyTestingDataGrabber() as dg:
-        input_ = dg["sub-01"]
+        element_data = DefaultDataReader().fit_transform(dg["sub-01"])
-
+        # Update workdir to current test's tmp_path
-    input_ = DefaultDataReader().fit_transform(input_)
-    WorkDirManager().workdir = tmp_path
-    marker_python = ALFFParcels(
-        parcellation=_PARCELLATION,
-        method="mean",
-        use_afni=False,
-        fractional=fractional,
-    )
-    python_values = marker_python.fit_transform(input_)["BOLD"]["data"]
-
-    assert marker_python.use_afni is False
-    assert python_values.ndim == 2
-    assert python_values.shape == (1, 100)
-
-    marker_afni = ALFFParcels(
-        parcellation=_PARCELLATION,
-        method="mean",
-        use_afni=True,
-        fractional=fractional,
-    )
-    afni_values = marker_afni.fit_transform(input_)["BOLD"]["data"]
-
-    assert marker_afni.use_afni is True
-    assert afni_values.ndim == 2
-    assert afni_values.shape == (1, 100)
-
-    r, p = pearsonr(python_values[0], afni_values[0])
-    logger.info(f"Correlation between python and afni: {r} (p={p})")
-    assert r > 0.99
-
-
-def test_ALFFParcels_storage(
-    tmp_path: Path,
-) -> None:
-    """Test ALFFParcels storage.
-
-    Parameters
-    ----------
-    tmp_path : pathlib.Path
-        The path to the test directory.
-
-    """
-    with PartlyCloudyTestingDataGrabber() as dg:
-        # Use first subject
-        input_ = dg["sub-01"]
-        input_ = DefaultDataReader().fit_transform(input_)
         WorkDirManager().workdir = tmp_path
-        marker = ALFFParcels(
-            parcellation=_PARCELLATION,
-            method="mean",
-            use_afni=False,
-            fractional=True,
-        )
-        storage = SQLiteFeatureStorage(tmp_path / "alff_parcels.sqlite")
 
-        # Fit transform marker on data with storage
+        # Initialize marker
-        marker.fit_transform(
+        junifer_marker = ALFFParcels(
-            input=input_,
+            parcellation=PARCELLATION,
-            storage=storage,
+            fractional=fractional,
+            using="junifer",
         )
+        # Fit transform marker on data
+        junifer_output = junifer_marker.fit_transform(element_data)
+        # Get BOLD output
+        junifer_output_bold = junifer_output["BOLD"]
+
+        # Initialize marker
+        afni_marker = ALFFParcels(
+            parcellation=PARCELLATION,
+            fractional=fractional,
+            using="afni",
+        )
+        # Fit transform marker on data
+        afni_output = afni_marker.fit_transform(element_data)
+        # Get BOLD output
+        afni_output_bold = afni_output["BOLD"]
+
+        # Check for Pearson correlation coefficient
+        r, _ = sp.stats.pearsonr(
+            junifer_output_bold["data"][0],
+            afni_output_bold["data"][0],
+        )
+        assert r > 0.99

junifer/markers/falff/tests/test_falff_spheres.py

@@ -1,14 +1,14 @@
-"""Provide test for sphere-aggregated (f)ALFF."""
+"""Provide tests for ALFFSpheres."""
 
 # Authors: Federico Raimondo <f.raimondo@fz-juelich.de>
 #          Synchon Mandal <s.mandal@fz-juelich.de>
 # License: AGPL
 
+import logging
 from pathlib import Path
 
 import pytest
-from numpy.testing import assert_array_equal
+import scipy as sp
-from scipy.stats import pearsonr
 
 from junifer.datareader import DefaultDataReader
 from junifer.markers.falff import ALFFSpheres
@@ -16,169 +16,114 @@ from junifer.pipeline import WorkDirManager
 from junifer.pipeline.utils import _check_afni
 from junifer.storage import SQLiteFeatureStorage
 from junifer.testing.datagrabbers import PartlyCloudyTestingDataGrabber
-from junifer.utils import logger
 
 
-_COORDINATES = "DMNBuckner"
+COORDINATES = "DMNBuckner"
 
 
-def test_ALFFSpheres_python(tmp_path: Path) -> None:
+def test_ALFFSpheres(caplog: pytest.LogCaptureFixture, tmp_path: Path) -> None:
-    """Test ALFFSpheres using python.
+    """Test ALFFSpheres.
 
     Parameters
     ----------
+    caplog : pytest.LogCaptureFixture
+        The pytest.LogCaptureFixture object.
     tmp_path : pathlib.Path
         The path to the test directory.
 
     """
-    with PartlyCloudyTestingDataGrabber() as dg:
+    with caplog.at_level(logging.DEBUG):
-        input_ = dg["sub-01"]
+        with PartlyCloudyTestingDataGrabber() as dg:
+            element_data = DefaultDataReader().fit_transform(dg["sub-01"])
+            # Update workdir to current test's tmp_path
+            WorkDirManager().workdir = tmp_path
 
-    input_ = DefaultDataReader().fit_transform(input_)
+            # Initialize marker
+            marker = ALFFSpheres(
+                coords=COORDINATES,
+                fractional=False,
+                using="junifer",
+                radius=5.0,
+            )
+            # Fit transform marker on data
+            output = marker.fit_transform(element_data)
 
-    WorkDirManager().workdir = tmp_path
+            assert "Creating cache" in caplog.text
-    marker = ALFFSpheres(
-        coords=_COORDINATES,
-        radius=5,
-        method="mean",
-        use_afni=False,
-        fractional=False,
-    )
-    python_values = marker.fit_transform(input_)["BOLD"]["data"]
 
-    assert marker.use_afni is False
+            # Get BOLD output
-    assert python_values.ndim == 2
+            assert "BOLD" in output
-    assert python_values.shape == (1, 6)
+            output_bold = output["BOLD"]
+            # Assert BOLD output keys
+            assert "data" in output_bold
+            assert "col_names" in output_bold
+
+            output_bold_data = output_bold["data"]
+            # Assert BOLD output data dimension
+            assert output_bold_data.ndim == 2
+            assert output_bold_data.shape == (1, 6)
+
+            # Reset log capture
+            caplog.clear()
+            # Initialize storage
+            storage = SQLiteFeatureStorage(tmp_path / "falff_spheres.sqlite")
+            # Fit transform marker on data with storage
+            marker.fit_transform(
+                input=element_data,
+                storage=storage,
+            )
+            # Cache working correctly
+            assert "Creating cache" not in caplog.text
 
 
 @pytest.mark.skipif(
-    _check_afni() is False, reason="requires afni to be in PATH"
+    _check_afni() is False, reason="requires AFNI to be in PATH"
 )
-def test_ALFFSpheres_afni(tmp_path: Path) -> None:
+@pytest.mark.parametrize(
+    "fractional", [True, False], ids=["fractional", "non-fractional"]
+)
+def test_ALFFSpheres_comparison(tmp_path: Path, fractional: bool) -> None:
     """Test ALFFSpheres using afni.
 
     Parameters
     ----------
     tmp_path : pathlib.Path
-        The path to the test directory.
+         The path to the test directory.
-
-    """
-    with PartlyCloudyTestingDataGrabber() as dg:
-        input = dg["sub-01"]
-
-    input = DefaultDataReader().fit_transform(input)
-
-    WorkDirManager().workdir = tmp_path
-    marker = ALFFSpheres(
-        coords=_COORDINATES,
-        radius=5,
-        method="mean",
-        use_afni=True,
-        fractional=False,
-    )
-    assert marker.use_afni is True
-    afni_values = marker.fit_transform(input)["BOLD"]["data"]
-
-    assert afni_values.ndim == 2
-    assert afni_values.shape == (1, 6)
-
-    # Again, should be blazing fast
-    marker = ALFFSpheres(
-        coords=_COORDINATES,
-        radius=5,
-        method="mean",
-        fractional=False,
-    )
-    assert marker.use_afni is None
-    afni_values2 = marker.fit_transform(input)["BOLD"]["data"]
-    assert marker.use_afni is True
-    assert_array_equal(afni_values, afni_values2)
-
-
-@pytest.mark.skipif(
-    _check_afni() is False, reason="requires afni to be in PATH"
-)
-@pytest.mark.parametrize(
-    "fractional", [True, False], ids=["fractional", "non-fractional"]
-)
-def test_ALFFSpheres_python_vs_afni(
-    tmp_path: Path,
-    fractional: bool,
-) -> None:
-    """Test ALFFSpheres python vs afni results.
-
-    Parameters
-    ----------
-    tmp_path : pathlib.Path
-        The path to the test directory.
     fractional : bool
-        Whether to compute fractional ALFF or not.
+         Whether to compute fractional ALFF or not.
 
     """
     with PartlyCloudyTestingDataGrabber() as dg:
-        input_ = dg["sub-01"]
+        element_data = DefaultDataReader().fit_transform(dg["sub-01"])
-
+        # Update workdir to current test's tmp_path
-    input_ = DefaultDataReader().fit_transform(input_)
-
-    WorkDirManager().workdir = tmp_path
-    marker_python = ALFFSpheres(
-        coords=_COORDINATES,
-        radius=5,
-        method="mean",
-        use_afni=False,
-        fractional=fractional,
-    )
-    python_values = marker_python.fit_transform(input_)["BOLD"]["data"]
-
-    assert marker_python.use_afni is False
-    assert python_values.ndim == 2
-    assert python_values.shape == (1, 6)
-
-    marker_afni = ALFFSpheres(
-        coords=_COORDINATES,
-        radius=5,
-        method="mean",
-        use_afni=True,
-        fractional=fractional,
-    )
-    afni_values = marker_afni.fit_transform(input_)["BOLD"]["data"]
-
-    assert marker_afni.use_afni is True
-    assert afni_values.ndim == 2
-    assert afni_values.shape == (1, 6)
-
-    r, p = pearsonr(python_values[0], afni_values[0])
-    logger.info(f"Correlation between python and afni: {r} (p={p})")
-    assert r > 0.99
-
-
-def test_ALFFSpheres_storage(
-    tmp_path: Path,
-) -> None:
-    """Test ALFFSpheres storage.
-
-    Parameters
-    ----------
-    tmp_path : pathlib.Path
-        The path to the test directory.
-
-    """
-    with PartlyCloudyTestingDataGrabber() as dg:
-        # Use first subject
-        input_ = dg["sub-01"]
-        input_ = DefaultDataReader().fit_transform(input_)
         WorkDirManager().workdir = tmp_path
-        marker = ALFFSpheres(
-            coords=_COORDINATES,
-            radius=5,
-            method="mean",
-            use_afni=False,
-            fractional=True,
-        )
-        storage = SQLiteFeatureStorage(tmp_path / "alff_parcels.sqlite")
 
-        # Fit transform marker on data with storage
+        # Initialize marker
-        marker.fit_transform(
+        junifer_marker = ALFFSpheres(
-            input=input_,
+            coords=COORDINATES,
-            storage=storage,
+            fractional=fractional,
+            using="junifer",
+            radius=5.0,
         )
+        # Fit transform marker on data
+        junifer_output = junifer_marker.fit_transform(element_data)
+        # Get BOLD output
+        junifer_output_bold = junifer_output["BOLD"]
+
+        # Initialize marker
+        afni_marker = ALFFSpheres(
+            coords=COORDINATES,
+            fractional=fractional,
+            using="afni",
+            radius=5.0,
+        )
+        # Fit transform marker on data
+        afni_output = afni_marker.fit_transform(element_data)
+        # Get BOLD output
+        afni_output_bold = afni_output["BOLD"]
+
+        # Check for Pearson correlation coefficient
+        r, _ = sp.stats.pearsonr(
+            junifer_output_bold["data"][0],
+            afni_output_bold["data"][0],
+        )
+        assert r > 0.99

junifer/markers/reho/_afni_reho.py

@@ -0,0 +1,192 @@
+"""Provide class for computing regional homogeneity (ReHo) using AFNI."""
+
+# Authors: Synchon Mandal <s.mandal@fz-juelich.de>
+# License: AGPL
+
+from functools import lru_cache
+from pathlib import Path
+from typing import (
+    TYPE_CHECKING,
+    ClassVar,
+    Dict,
+    List,
+    Optional,
+    Tuple,
+    Union,
+)
+
+import nibabel as nib
+
+from ...pipeline import WorkDirManager
+from ...pipeline.singleton import singleton
+from ...utils import logger, run_ext_cmd
+
+
+if TYPE_CHECKING:
+    from nibabel import Nifti1Image
+
+
+@singleton
+class AFNIReHo:
+    """Class for computing ReHo using AFNI.
+
+    This class uses AFNI's 3dReHo to compute ReHo. It's designed as a singleton
+    with caching for efficient computation.
+
+    """
+
+    _EXT_DEPENDENCIES: ClassVar[List[Dict[str, Union[str, List[str]]]]] = [
+        {
+            "name": "afni",
+            "commands": ["3dReHo", "3dAFNItoNIFTI"],
+        },
+    ]
+
+    def __del__(self) -> None:
+        """Terminate the class."""
+        # Clear the computation cache
+        logger.debug("Clearing cache for ReHo computation via AFNI")
+        self.compute.cache_clear()
+
+    @lru_cache(maxsize=None, typed=True)
+    def compute(
+        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,
+    ) -> Tuple["Nifti1Image", Path]:
+        """Compute ReHo map.
+
+        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
+            The ReHo map as NIfTI.
+        pathlib.Path
+            The path to the ReHo map as NIfTI.
+
+        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
+
+        """
+        logger.debug("Creating cache for ReHo computation via AFNI")
+
+        # Create component-scoped tempdir
+        tempdir = WorkDirManager().get_tempdir(prefix="afni_reho")
+
+        # Save target data to a component-scoped tempfile
+        nifti_in_file_path = tempdir / "input.nii"  # needs to be .nii
+        nib.save(data, nifti_in_file_path)
+
+        # Set 3dReHo command
+        reho_out_path_prefix = tempdir / "reho"
+        reho_cmd = [
+            "3dReHo",
+            f"-prefix {reho_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
+        run_ext_cmd(name="3dReHo", cmd=reho_cmd)
+
+        # Create element-scoped tempdir so that the ReHo map is
+        # available later as nibabel stores file path reference for
+        # loading on computation
+        element_tempdir = WorkDirManager().get_element_tempdir(
+            prefix="afni_reho"
+        )
+        # Convert afni to nifti
+        reho_afni_to_nifti_out_path = (
+            element_tempdir / "output.nii"  # needs to be .nii
+        )
+        convert_cmd = [
+            "3dAFNItoNIFTI",
+            f"-prefix {reho_afni_to_nifti_out_path.resolve()}",
+            f"{reho_out_path_prefix}+tlrc.BRIK",
+        ]
+        # Call 3dAFNItoNIFTI
+        run_ext_cmd(name="3dAFNItoNIFTI", cmd=convert_cmd)
+
+        # Load nifti
+        output_data = nib.load(reho_afni_to_nifti_out_path)
+
+        # Delete tempdir
+        WorkDirManager().delete_tempdir(tempdir)
+
+        return output_data, reho_afni_to_nifti_out_path  # type: ignore

junifer/markers/reho/_junifer_reho.py

@@ -0,0 +1,281 @@
+"""Provide class for computing regional homogeneity (ReHo) using junifer."""
+
+# Authors: Synchon Mandal <s.mandal@fz-juelich.de>
+# License: AGPL
+
+from functools import lru_cache
+from itertools import product
+from pathlib import Path
+from typing import (
+    TYPE_CHECKING,
+    ClassVar,
+    Set,
+    Tuple,
+)
+
+import nibabel as nib
+import numpy as np
+import scipy as sp
+from nilearn import image as nimg
+from nilearn import masking as nmask
+
+from ...pipeline import WorkDirManager
+from ...pipeline.singleton import singleton
+from ...utils import logger, raise_error
+
+
+if TYPE_CHECKING:
+    from nibabel import Nifti1Image
+
+
+@singleton
+class JuniferReHo:
+    """Class for computing ReHo using junifer.
+
+    It's designed as a singleton with caching for efficient computation.
+
+    """
+
+    _DEPENDENCIES: ClassVar[Set[str]] = {"numpy", "nilearn", "scipy"}
+
+    def __del__(self) -> None:
+        """Terminate the class."""
+        # Clear the computation cache
+        logger.debug("Clearing cache for ReHo computation via junifer")
+        self.compute.cache_clear()
+
+    @lru_cache(maxsize=None, typed=True)
+    def compute(
+        self,
+        data: "Nifti1Image",
+        nneigh: int = 27,
+    ) -> Tuple["Nifti1Image", Path]:
+        """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
+            The ReHo map as NIfTI.
+        pathlib.Path
+            The path to the ReHo map as NIfTI.
+
+        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.debug("Creating cache for ReHo computation via junifer")
+
+        # 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 = sp.stats.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)
+
+        # TODO(synchon): this will give incorrect results if
+        # template doesn't match, hence needs to be changed
+        # after #299 is merged
+        # Calculate whole brain mask
+        mni152_whole_brain_mask = nmask.compute_brain_mask(
+            target_img=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,
+            )
+
+        # Create new image like target image
+        output_data = nimg.new_img_like(
+            ref_niimg=data,
+            data=reho_map,
+            copy_header=False,
+        )
+
+        # Create element-scoped tempdir so that the ReHo map is
+        # available later as nibabel stores file path reference for
+        # loading on computation
+        element_tempdir = WorkDirManager().get_element_tempdir(
+            prefix="junifer_reho"
+        )
+        output_path = element_tempdir / "output.nii.gz"
+        # Save computed data to file
+        nib.save(output_data, output_path)
+
+        return output_data, output_path  # type: ignore
+
+
+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_ranks : 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_base.py

@@ -13,12 +13,14 @@ from typing import (
     List,
     Optional,
     Tuple,
+    Type,
     Union,
 )
 
 from ...utils import logger, raise_error
 from ..base import BaseMarker
-from .reho_estimator import ReHoEstimator
+from ._afni_reho import AFNIReHo
+from ._junifer_reho import JuniferReHo
 
 
 if TYPE_CHECKING:
@@ -30,32 +32,47 @@ class ReHoBase(BaseMarker):
 
     Parameters
     ----------
-    use_afni : bool, optional
+    using : {"junifer", "afni"}
-        Whether to use AFNI for computing. If None, will use AFNI only
+        Implementation to use for computing ReHo:
-        if available (default None).
+
+        * "junifer" : Use ``junifer``'s own ReHo implementation
+        * "afni" : Use AFNI's ``3dReHo``
+
     name : str, optional
         The name of the marker. If None, it will use the class name
         (default None).
 
+    Raises
+    ------
+    ValueError
+        If ``using`` is invalid.
+
     """
 
-    _EXT_DEPENDENCIES: ClassVar[
+    _CONDITIONAL_DEPENDENCIES: ClassVar[List[Dict[str, Union[str, Type]]]] = [
-        List[Dict[str, Union[str, bool, List[str]]]]
-    ] = [
         {
-            "name": "afni",
+            "using": "afni",
-            "optional": True,
+            "depends_on": AFNIReHo,
-            "commands": ["3dReHo", "3dAFNItoNIFTI"],
+        },
+        {
+            "using": "junifer",
+            "depends_on": JuniferReHo,
         },
     ]
 
     def __init__(
         self,
-        use_afni: Optional[bool] = None,
+        using: str,
         name: Optional[str] = None,
     ) -> None:
+        # Validate `using` parameter
+        valid_using = [dep["using"] for dep in self._CONDITIONAL_DEPENDENCIES]
+        if using not in valid_using:
+            raise_error(
+                f"Invalid value for `using`, should be one of: {valid_using}"
+            )
+        self.using = using
         super().__init__(on="BOLD", name=name)
-        self.use_afni = use_afni
 
     def get_valid_inputs(self) -> List[str]:
         """Get valid data types for input.
@@ -84,12 +101,12 @@ class ReHoBase(BaseMarker):
         """
         return "vector"
 
-    def compute_reho_map(
+    def _compute(
         self,
-        input: Dict[str, Any],
+        input_data: Dict[str, Any],
         **reho_params: Any,
     ) -> Tuple["Nifti1Image", Path]:
-        """Compute.
+        """Compute voxel-wise ReHo.
 
         Calculates Kendall's W per voxel using neighborhood voxels.
         Instead of the time series values themselves, Kendall's W uses the
@@ -100,7 +117,7 @@ class ReHoBase(BaseMarker):
 
         Parameters
         ----------
-        input : dict
+        input_data : dict
             The BOLD data as dictionary.
         **reho_params : dict
             Extra keyword arguments for ReHo.
@@ -122,20 +139,25 @@ class ReHoBase(BaseMarker):
                https://doi.org/10.1177/1073858415595004
 
         """
-        if self.use_afni is None:
+        logger.debug("Calculating voxel-wise ReHo")
-            raise_error(
+
-                "Parameter `use_afni` must be set to True or False in order "
+        # Conditional estimator
-                "to compute this marker. It is currently set to None (default "
+        if self.using == "afni":
-                "behaviour). This is intended to be for auto-detection. In "
+            estimator = AFNIReHo()
-                "order for that to happen, please call the `validate` method "
+        elif self.using == "junifer":
-                "before calling the `compute` method."
+            estimator = JuniferReHo()
-            )
+        # Compute reho
-        logger.info("Calculating ReHO map.")
+        reho_map, reho_map_path = estimator.compute(  # type: ignore
-        # Initialize reho estimator
+            data=input_data["data"],
-        reho_estimator = ReHoEstimator()
-        # Fit-transform reho estimator
-        return reho_estimator.fit_transform(
-            use_afni=self.use_afni,
-            input_data=input,
             **reho_params,
         )
+
+        # If the input data space is native already, the original path should
+        # be propagated down as it might be required for transforming
+        # parcellation / coordinates to native space, else the reho map
+        # path should be passed for use later if required.
+        # TODO(synchon): will be taken care in #292
+        if input_data["space"] == "native":
+            return reho_map, input_data["path"]
+
+        return reho_map, reho_map_path

junifer/markers/reho/reho_estimator.py

@@ -1,510 +0,0 @@
-"""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
-from functools import lru_cache
-from itertools import product
-from pathlib import Path
-from typing import TYPE_CHECKING, Any, Dict, List, Optional, Tuple
-
-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 ...pipeline import WorkDirManager
-from ...pipeline.singleton import singleton
-from ...utils import logger, raise_error, run_ext_cmd
-
-
-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 = None
-
-    def __del__(self) -> None:
-        """Cleanup."""
-        # Delete temporary directory and ignore errors for read-only files
-        if self.temp_dir_path is not None:
-            WorkDirManager().delete_tempdir(self.temp_dir_path)
-
-    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,
-    ) -> Tuple["Nifti1Image", Path]:
-        """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
-            The ReHo map as NIfTI.
-        pathlib.Path
-            The path to the ReHo map as NIfTI.
-
-        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
-
-        """
-        # Note: self.temp_dir_path is sure to exist before proceeding, so
-        #       types checks are ignored further on.
-
-        # Save niimg to nii.gz
-        nifti_in_file_path = self.temp_dir_path / "input.nii"  # type: ignore
-        nib.save(data, nifti_in_file_path)
-
-        # Set 3dReHo command
-        reho_afni_out_path_prefix = self.temp_dir_path / "reho"  # type: ignore
-        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
-        run_ext_cmd(name="3dReHo", cmd=reho_cmd)
-
-        # SHA256 for bypassing memmap
-        sha256_params = hashlib.sha256(bytes(" ".join(reho_cmd), "utf-8"))
-        # Create element-scoped tempdir so that the ReHo map is
-        # available later as get_coordinates and the like need it
-        # in ReHoSpheres and the like to transform to other template
-        # spaces
-        element_tempdir = WorkDirManager().get_element_tempdir(
-            prefix="reho_afni"
-        )
-        # Convert afni to nifti
-        reho_afni_to_nifti_out_path = (
-            element_tempdir
-            / f"output_{sha256_params.hexdigest()}.nii"  # type: ignore
-        )
-        convert_cmd: List[str] = [
-            "3dAFNItoNIFTI",
-            f"-prefix {reho_afni_to_nifti_out_path.resolve()}",
-            f"{reho_afni_out_path_prefix}+tlrc.BRIK",
-        ]
-        # Call 3dAFNItoNIFTI
-        run_ext_cmd(name="3dAFNItoNIFTI", cmd=convert_cmd)
-
-        # Cleanup intermediate files
-        for fname in self.temp_dir_path.glob("reho*"):  # type: ignore
-            fname.unlink()
-
-        # Load nifti
-        output_data = nib.load(reho_afni_to_nifti_out_path)
-
-        return output_data, reho_afni_to_nifti_out_path  # type: ignore
-
-    def _compute_reho_python(
-        self,
-        data: "Nifti1Image",
-        nneigh: int = 27,
-    ) -> Tuple["Nifti1Image", Path]:
-        """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
-            The ReHo map as NIfTI.
-        pathlib.Path
-            The path to the ReHo map as NIfTI.
-
-        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)
-        # Create element-scoped tempdir so that the ReHo map is
-        # available later as get_coordinates and the like need it
-        # in ReHoSpheres and the like to transform to other template
-        # spaces
-        element_tempdir = WorkDirManager().get_element_tempdir(
-            prefix="reho_python"
-        )
-        output_path = element_tempdir / "reho_map_python.nii.gz"
-        nib.save(output, output_path)
-        return output, output_path  # type: ignore
-
-    @lru_cache(maxsize=None, typed=True)
-    def _compute(
-        self,
-        use_afni: bool,
-        data: "Nifti1Image",
-        **reho_params: Any,
-    ) -> Tuple["Nifti1Image", Path]:
-        """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
-            The ReHo map as NIfTI.
-        pathlib.Path
-            The path to the ReHo map as NIfTI.
-
-        """
-        if use_afni:
-            # Create new temporary directory before using AFNI
-            self.temp_dir_path = WorkDirManager().get_tempdir(prefix="reho")
-            return self._compute_reho_afni(data, **reho_params)
-
-        return self._compute_reho_python(data, **reho_params)
-
-    def fit_transform(
-        self,
-        use_afni: bool,
-        input_data: Dict[str, Any],
-        **reho_params: Any,
-    ) -> Tuple["Nifti1Image", Path]:
-        """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
-            The ReHo map as NIfTI.
-        pathlib.Path
-            The path to the ReHo map as NIfTI or the input data path if the
-            input data space is "native".
-
-        """
-        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 for {self._file_path}.")
-            # Clear the cache
-            self._compute.cache_clear()
-            # Clear temporary directory files
-            if self.temp_dir_path is not None:
-                WorkDirManager().delete_tempdir(self.temp_dir_path)
-            # Set the new file path
-            self._file_path = bold_path
-        else:
-            logger.info(f"Using ReHo map cache for {self._file_path}.")
-        # Compute
-        reho_map, reho_map_path = self._compute(
-            use_afni, bold_data, **reho_params
-        )
-        # If the input data space is native already, the original path should
-        # be propagated down as it might be required for transforming
-        # coordinates to native space via get_coordinates(), else the reho map
-        # path should be passed for use later if required.
-        if input_data["space"] == "native":
-            return reho_map, input_data["path"]
-
-        return reho_map, reho_map_path
-
-
-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_ranks : 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

@@ -20,15 +20,18 @@ class ReHoParcels(ReHoBase):
 
     Parameters
     ----------
-    parcellation : str
+    parcellation : str or list of str
-        The name of the parcellation. Check valid options by calling
+        The name(s) of the parcellation(s). Check valid options by calling
         :func:`.list_parcellations`.
-    use_afni : bool, optional
+    using : {"junifer", "afni"}
-        Whether to use AFNI for computing. If None, will use AFNI only
+        Implementation to use for computing ReHo:
-        if available (default None).
+
+        * "junifer" : Use ``junifer``'s own ReHo implementation
+        * "afni" : Use AFNI's ``3dReHo``
+
     reho_params : dict, optional
         Extra parameters for computing ReHo map as a dictionary (default None).
-        If ``use_afni = True``, then the valid keys are:
+        If ``using="afni"``, then the valid keys are:
 
         * ``nneigh`` : {7, 19, 27}, optional (default 27)
             Number of voxels in the neighbourhood, inclusive. Can be:
@@ -58,7 +61,7 @@ class ReHoParcels(ReHoBase):
             The number of voxels for +/- z-axis of cuboidal volumes
             (default None).
 
-        else if ``use_afni = False``, then the valid keys are:
+        else if ``using="junifer"``, then the valid keys are:
 
         * ``nneigh`` : {7, 19, 27, 125}, optional (default 27)
             Number of voxels in the neighbourhood, inclusive. Can be:
@@ -86,20 +89,21 @@ class ReHoParcels(ReHoBase):
 
     def __init__(
         self,
-        parcellation: str,
+        parcellation: Union[str, List[str]],
-        use_afni: Optional[bool] = None,
+        using: str,
         reho_params: Optional[Dict] = None,
         agg_method: str = "mean",
         agg_method_params: Optional[Dict] = None,
         masks: Union[str, Dict, List[Union[Dict, str]], None] = None,
         name: Optional[str] = None,
     ) -> None:
+        # Superclass init first to validate `using` parameter
+        super().__init__(using=using, name=name)
         self.parcellation = parcellation
         self.reho_params = reho_params
         self.agg_method = agg_method
         self.agg_method_params = agg_method_params
         self.masks = masks
-        super().__init__(use_afni=use_afni, name=name)
 
     def compute(
         self,
@@ -125,17 +129,19 @@ class ReHoParcels(ReHoBase):
             * ``col_names`` : the column labels for the parcels as a list
 
         """
-        logger.info("Calculating ReHo for parcels.")
+        logger.info("Calculating ReHo for parcels")
-        # Calculate reho map
+
+        # Compute voxelwise reho
         # If the input data space is "native", then reho_file_path points to
-        # the input data path as it might be required to use in
+        # the input data path as it might be required for parcellation
-        # get_coordinates() for transforming coordinates to native space.
+        # transformation to native space.
         if self.reho_params is not None:
-            reho_map, reho_file_path = self.compute_reho_map(
+            reho_map, reho_file_path = self._compute(
-                input=input, **self.reho_params
+                input_data=input, **self.reho_params
             )
         else:
-            reho_map, reho_file_path = self.compute_reho_map(input=input)
+            reho_map, reho_file_path = self._compute(input_data=input)
+
         # Initialize parcel aggregation
         parcel_aggregation = ParcelAggregation(
             parcellation=self.parcellation,
@@ -148,7 +154,6 @@ class ReHoParcels(ReHoBase):
         parcel_aggregation_input = dict(input.items())
         parcel_aggregation_input["data"] = reho_map
         parcel_aggregation_input["path"] = reho_file_path
-
         output = parcel_aggregation.compute(
             input=parcel_aggregation_input,
             extra_input=extra_input,

junifer/markers/reho/reho_spheres.py

@@ -23,6 +23,12 @@ class ReHoSpheres(ReHoBase):
     coords : str
         The name of the coordinates list to use. See
         :func:`.list_coordinates` for options.
+    using : {"junifer", "afni"}
+        Implementation to use for computing ReHo:
+
+        * "junifer" : Use ``junifer``'s own ReHo implementation
+        * "afni" : Use AFNI's ``3dReHo``
+
     radius : float, optional
         The radius of the sphere in millimeters. If None, the signal will be
         extracted from a single voxel. See
@@ -36,7 +42,7 @@ class ReHoSpheres(ReHoBase):
         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:
+        If ``using="afni"``, then the valid keys are:
 
         * ``nneigh`` : {7, 19, 27}, optional (default 27)
             Number of voxels in the neighbourhood, inclusive. Can be:
@@ -66,7 +72,7 @@ class ReHoSpheres(ReHoBase):
             The number of voxels for +/- z-axis of cuboidal volumes
             (default None).
 
-        else if ``use_afni = False``, then the valid keys are:
+        else if ``using="junifer"``, then the valid keys are:
 
         * ``nneigh`` : {7, 19, 27, 125}, optional (default 27)
             Number of voxels in the neighbourhood, inclusive. Can be:
@@ -95,15 +101,17 @@ class ReHoSpheres(ReHoBase):
     def __init__(
         self,
         coords: str,
+        using: str,
         radius: Optional[float] = None,
         allow_overlap: bool = False,
-        use_afni: Optional[bool] = None,
         reho_params: Optional[Dict] = None,
         agg_method: str = "mean",
         agg_method_params: Optional[Dict] = None,
         masks: Union[str, Dict, List[Union[Dict, str]], None] = None,
         name: Optional[str] = None,
     ) -> None:
+        # Superclass init first to validate `using` parameter
+        super().__init__(using=using, name=name)
         self.coords = coords
         self.radius = radius
         self.allow_overlap = allow_overlap
@@ -111,7 +119,6 @@ class ReHoSpheres(ReHoBase):
         self.agg_method = agg_method
         self.agg_method_params = agg_method_params
         self.masks = masks
-        super().__init__(use_afni=use_afni, name=name)
 
     def compute(
         self,
@@ -137,17 +144,19 @@ class ReHoSpheres(ReHoBase):
             * ``col_names`` : the column labels for the spheres as a list
 
         """
-        logger.info("Calculating ReHo for spheres.")
+        logger.info("Calculating ReHo for spheres")
-        # Calculate reho map
+
+        # Compute voxelwise reho
         # If the input data space is "native", then reho_file_path points to
-        # the input data path as it might be required to use in
+        # the input data path as it might be required for coordinates
-        # get_coordinates() for transforming coordinates to native space.
+        # transformation to native space.
         if self.reho_params is not None:
-            reho_map, reho_file_path = self.compute_reho_map(
+            reho_map, reho_file_path = self._compute(
-                input=input, **self.reho_params
+                input_data=input, **self.reho_params
             )
         else:
-            reho_map, reho_file_path = self.compute_reho_map(input=input)
+            reho_map, reho_file_path = self._compute(input_data=input)
+
         # Initialize sphere aggregation
         sphere_aggregation = SphereAggregation(
             coords=self.coords,

junifer/markers/reho/tests/test_reho_estimator.py

@@ -1,314 +0,0 @@
-"""Provide tests for ReHo map compute comparison."""
-
-# Authors: Synchon Mandal <s.mandal@fz-juelich.de>
-# License: AGPL
-
-import time
-from pathlib import Path
-
-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 import WorkDirManager
-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(tmp_path: Path) -> None:
-    """Test that the cache works properly when using Python implementation.
-
-    Parameters
-    ----------
-    tmp_path : pathlib.Path
-        The path to the test directory.
-
-    """
-    # Get subject from datagrabber
-    with PartlyCloudyTestingDataGrabber() as dg:
-        subject = dg["sub-01"]
-    # Read data for subject
-    subject_data = DefaultDataReader().fit_transform(subject)
-    # Update workdir to current test's tmp_path
-    WorkDirManager().workdir = tmp_path
-    # Setup estimator
-    reho_estimator = ReHoEstimator()
-
-    # Compute without cache
-    first_tic = time.time()
-    (
-        reho_map_without_cache,
-        reho_map_without_cache_path,
-    ) = 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)
-    assert isinstance(reho_map_without_cache_path, Path)
-
-    # Compute again with cache, should be faster
-    second_tic = time.time()
-    (
-        reho_map_with_cache,
-        reho_map_with_cache_path,
-    ) = 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)
-    assert isinstance(reho_map_with_cache_path, Path)
-    # Check that cache is being used
-    assert (second_toc - second_tic) < ((first_toc - first_tic) / 1000)
-
-    # Change a parameter and compute again without cache
-    third_tic = time.time()
-    (
-        reho_map_with_partial_cache,
-        reho_map_with_partial_cache_path,
-    ) = 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)
-    assert isinstance(reho_map_with_partial_cache_path, Path)
-    # Should require more time
-    assert (third_toc - third_tic) > ((first_toc - first_tic) / 10)
-
-    # Compute again with cache, should be faster
-    fourth_tic = time.time()
-    (
-        reho_map_with_new_cache,
-        reho_map_with_new_cache_path,
-    ) = 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)
-    assert isinstance(reho_map_with_new_cache_path, Path)
-    # Should require less time
-    assert (fourth_toc - fourth_tic) < ((first_toc - first_tic) / 1000)
-
-    # Change the data and 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_map_with_different_cache_path,
-    ) = 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)
-    assert isinstance(reho_map_with_different_cache_path, Path)
-    # Should take less time
-    assert (fifth_toc - fifth_tic) > ((first_toc - first_tic) / 10)
-
-
-@pytest.mark.skipif(
-    _check_afni() is False, reason="requires afni to be in PATH"
-)
-def test_reho_estimator_cache_afni(tmp_path: Path) -> None:
-    """Test that the cache works properly when using afni.
-
-    Parameters
-    ----------
-    tmp_path : pathlib.Path
-        The path to the test directory.
-
-    """
-    # Get subject from datagrabber
-    with PartlyCloudyTestingDataGrabber() as dg:
-        subject = dg["sub-01"]
-    # Read data for subject
-    subject_data = DefaultDataReader().fit_transform(subject)
-    # Update workdir to current test's tmp_path
-    WorkDirManager().workdir = tmp_path
-    # Setup estimator
-    reho_estimator = ReHoEstimator()
-
-    # Compute without cache
-    first_tic = time.time()
-    (
-        reho_map_without_cache,
-        reho_map_without_cache_path,
-    ) = 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)
-    assert isinstance(reho_map_without_cache_path, Path)
-    # Count intermediate files
-    n_files = len(list(reho_estimator.temp_dir_path.glob("*")))
-    assert n_files == 1  # only input file
-
-    # Compute again with cache, should be faster
-    second_tic = time.time()
-    (
-        reho_map_with_cache,
-        reho_map_with_cache_path,
-    ) = 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 isinstance(reho_map_with_cache_path, Path)
-    assert (second_toc - second_tic) < ((first_toc - first_tic) / 1000)
-    # Count intermediate files
-    n_files = len(list(reho_estimator.temp_dir_path.glob("*")))
-    assert n_files == 1  # only input file
-
-    # Change a parameter and compute again without cache
-    third_tic = time.time()
-    (
-        reho_map_with_partial_cache,
-        reho_map_with_partial_cache_path,
-    ) = 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)
-    assert isinstance(reho_map_with_partial_cache_path, Path)
-    # Should require more time
-    assert (third_toc - third_tic) > ((first_toc - first_tic) / 10)
-    # Count intermediate files
-    n_files = len(list(reho_estimator.temp_dir_path.glob("*")))
-    assert n_files == 1  # only input file
-
-    # Compute with cache, should be faster
-    fourth_tic = time.time()
-    (
-        reho_map_with_new_cache,
-        reho_map_with_new_cache_path,
-    ) = 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)
-    assert isinstance(reho_map_with_new_cache_path, Path)
-    # Should require less time
-    assert (fourth_toc - fourth_tic) < ((first_toc - first_tic) / 1000)
-    n_files = len(list(reho_estimator.temp_dir_path.glob("*")))
-    assert n_files == 1  # only input file
-
-    # Change the data and 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_map_with_different_cache_path,
-    ) = 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)
-    assert isinstance(reho_map_with_different_cache_path, Path)
-    # Should take less time
-    assert (fifth_toc - fifth_tic) > ((first_toc - first_tic) / 10)
-    # Count intermediate files
-    n_files = len(list(reho_estimator.temp_dir_path.glob("*")))
-    assert n_files == 1  # only input file
-
-
-@pytest.mark.skipif(
-    _check_afni() is False, reason="requires afni to be in PATH"
-)
-def test_reho_estimator_afni_vs_python(tmp_path: Path) -> None:
-    """Compare afni and Python implementations.
-
-    Parameters
-    ----------
-    tmp_path : pathlib.Path
-        The path to the test directory.
-
-    """
-    # Get subject from datagrabber
-    with PartlyCloudyTestingDataGrabber() as dg:
-        subject = dg["sub-01"]
-    # Read data for subject
-    subject_data = DefaultDataReader().fit_transform(subject)
-    # Update workdir to current test's tmp_path
-    WorkDirManager().workdir = tmp_path
-    # 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

@@ -1,73 +1,82 @@
-"""Provide tests for ReHo on parcels."""
+"""Provide tests for ReHoParcels."""
 
 # Authors: Synchon Mandal <s.mandal@fz-juelich.de>
 # License: AGPL
 
+import logging
 from pathlib import Path
 
 import pytest
-from nilearn import image as nimg
+import scipy as sp
-from scipy.stats import pearsonr
 
-from junifer.markers.reho.reho_parcels import ReHoParcels
+from junifer.datareader import DefaultDataReader
+from junifer.markers import ReHoParcels
 from junifer.pipeline import WorkDirManager
 from junifer.pipeline.utils import _check_afni
-from junifer.storage.sqlite import SQLiteFeatureStorage
+from junifer.storage import SQLiteFeatureStorage
 from junifer.testing.datagrabbers import SPMAuditoryTestingDataGrabber
 
 
 PARCELLATION = "Schaefer100x7"
 
 
-def test_reho_parcels_computation(tmp_path: Path) -> None:
+def test_ReHoParcels(caplog: pytest.LogCaptureFixture, tmp_path: Path) -> None:
-    """Test ReHoParcels fit-transform.
+    """Test ReHoParcels.
 
     Parameters
     ----------
+    caplog : pytest.LogCaptureFixture
+        The pytest.LogCaptureFixture object.
     tmp_path : pathlib.Path
         The path to the test directory.
 
     """
-    with SPMAuditoryTestingDataGrabber() as dg:
+    with caplog.at_level(logging.DEBUG):
-        # Use first subject
+        with SPMAuditoryTestingDataGrabber() as dg:
-        subject_data = dg["sub001"]
+            element_data = DefaultDataReader().fit_transform(dg["sub001"])
-        # Load image to memory
+            # Update workdir to current test's tmp_path
-        fmri_img = nimg.load_img(subject_data["BOLD"]["path"])
+            WorkDirManager().workdir = tmp_path
-        # Update workdir to current test's tmp_path
-        WorkDirManager().workdir = tmp_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": {},
-                    "space": "MNI",
-                }
-            }
-        )
-        # Get BOLD output
-        reho_parcels_output_bold = reho_parcels_output["BOLD"]
-        # Assert BOLD output keys
-        assert "data" in reho_parcels_output_bold
-        assert "col_names" in reho_parcels_output_bold
 
-        reho_parcels_output_bold_data = reho_parcels_output_bold["data"]
+            # Initialize marker
-        # Assert BOLD output data dimension
+            marker = ReHoParcels(parcellation=PARCELLATION, using="junifer")
-        assert reho_parcels_output_bold_data.ndim == 2
+            # Fit transform marker on data
-        # Assert BOLD output data is normalized
+            output = marker.fit_transform(element_data)
-        assert (reho_parcels_output_bold_data > 0).all() and (
+
-            reho_parcels_output_bold_data < 1
+            assert "Creating cache" in caplog.text
-        ).all()
+
+            # Get BOLD output
+            assert "BOLD" in output
+            output_bold = output["BOLD"]
+            # Assert BOLD output keys
+            assert "data" in output_bold
+            assert "col_names" in output_bold
+
+            output_bold_data = output_bold["data"]
+            # Assert BOLD output data dimension
+            assert output_bold_data.ndim == 2
+            # Assert BOLD output data is normalized
+            assert (output_bold_data > 0).all() and (
+                output_bold_data < 1
+            ).all()
+
+            # Reset log capture
+            caplog.clear()
+            # Initialize storage
+            storage = SQLiteFeatureStorage(tmp_path / "reho_parcels.sqlite")
+            # Fit transform marker on data with storage
+            marker.fit_transform(
+                input=element_data,
+                storage=storage,
+            )
+            # Cache working correctly
+            assert "Creating cache" not in caplog.text
 
 
 @pytest.mark.skipif(
-    _check_afni() is False, reason="requires afni to be in PATH"
+    _check_afni() is False, reason="requires AFNI to be in PATH"
 )
-def test_reho_parcels_computation_comparison(tmp_path: Path) -> None:
+def test_ReHoParcels_comparison(tmp_path: Path) -> None:
-    """Test ReHoParcels fit-transform implementation comparison.
+    """Test ReHoParcels implementation comparison.
 
     Parameters
     ----------
@@ -76,92 +85,29 @@ def test_reho_parcels_computation_comparison(tmp_path: Path) -> None:
 
     """
     with SPMAuditoryTestingDataGrabber() as dg:
-        # Use first subject
+        element_data = DefaultDataReader().fit_transform(dg["sub001"])
-        subject_data = dg["sub001"]
-        # Load image to memory
-        fmri_img = nimg.load_img(subject_data["BOLD"]["path"])
-
         # Update workdir to current test's tmp_path
         WorkDirManager().workdir = tmp_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": {},
-                    "space": "MNI",
-                }
-            }
-        )
-        # Get BOLD output
-        reho_parcels_output_bold_python = reho_parcels_output_python["BOLD"]
 
-        # Initialize marker with use_afni=True
+        # Initialize marker
-        reho_parcels_marker_afni = ReHoParcels(
+        junifer_marker = ReHoParcels(
-            parcellation=PARCELLATION, use_afni=True
+            parcellation=PARCELLATION, using="junifer"
         )
         # Fit transform marker on data
-        reho_parcels_output_afni = reho_parcels_marker_afni.fit_transform(
+        junifer_output = junifer_marker.fit_transform(element_data)
-            {
-                "BOLD": {
-                    "path": "/tmp",
-                    "data": fmri_img,
-                    "meta": {},
-                    "space": "MNI",
-                }
-            }
-        )
         # Get BOLD output
-        reho_parcels_output_bold_afni = reho_parcels_output_afni["BOLD"]
+        junifer_output_bold = junifer_output["BOLD"]
+
+        # Initialize marker
+        afni_marker = ReHoParcels(parcellation=PARCELLATION, using="afni")
+        # Fit transform marker on data
+        afni_output = afni_marker.fit_transform(element_data)
+        # Get BOLD output
+        afni_output_bold = afni_output["BOLD"]
 
         # Check for Pearson correlation coefficient
-        r, _ = pearsonr(
+        r, _ = sp.stats.pearsonr(
-            reho_parcels_output_bold_python["data"].flatten(),
+            junifer_output_bold["data"].flatten(),
-            reho_parcels_output_bold_afni["data"].flatten(),
+            afni_output_bold["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"])
-        # Update workdir to current test's tmp_path
-        WorkDirManager().workdir = tmp_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,
-                    "space": "MNI",
-                }
-            },
-            storage=reho_parcels_storage,
-        )

junifer/markers/reho/tests/test_reho_spheres.py

@@ -1,73 +1,83 @@
-"""Provide tests for ReHo on spheres."""
+"""Provide tests for ReHoSpheres."""
 
 # Authors: Synchon Mandal <s.mandal@fz-juelich.de>
 # License: AGPL
 
+import logging
 from pathlib import Path
 
 import pytest
-from nilearn import image as nimg
+import scipy as sp
-from scipy.stats import pearsonr
 
-from junifer.markers.reho.reho_spheres import ReHoSpheres
+from junifer.datareader import DefaultDataReader
+from junifer.markers import ReHoSpheres
 from junifer.pipeline import WorkDirManager
 from junifer.pipeline.utils import _check_afni
-from junifer.storage.sqlite import SQLiteFeatureStorage
+from junifer.storage import SQLiteFeatureStorage
 from junifer.testing.datagrabbers import SPMAuditoryTestingDataGrabber
 
 
 COORDINATES = "DMNBuckner"
 
 
-def test_reho_spheres_computation(tmp_path: Path) -> None:
+def test_ReHoSpheres(caplog: pytest.LogCaptureFixture, tmp_path: Path) -> None:
-    """Test ReHoSpheres fit-transform.
+    """Test ReHoSpheres.
 
     Parameters
     ----------
+    caplog : pytest.LogCaptureFixture
+        The pytest.LogCaptureFixture object.
     tmp_path : pathlib.Path
         The path to the test directory.
 
     """
-    with SPMAuditoryTestingDataGrabber() as dg:
+    with caplog.at_level(logging.DEBUG):
-        # Use first subject
+        with SPMAuditoryTestingDataGrabber() as dg:
-        subject_data = dg["sub001"]
+            element_data = DefaultDataReader().fit_transform(dg["sub001"])
-        # Load image to memory
+            # Update workdir to current test's tmp_path
-        fmri_img = nimg.load_img(subject_data["BOLD"]["path"])
+            WorkDirManager().workdir = tmp_path
-        # Update workdir to current test's tmp_path
+            # Initialize marker
-        WorkDirManager().workdir = tmp_path
+            marker = ReHoSpheres(
-        # Initialize marker
+                coords=COORDINATES, using="junifer", radius=10.0
-        reho_spheres_marker = ReHoSpheres(coords=COORDINATES, radius=10.0)
+            )
-        # Fit transform marker on data
+            # Fit transform marker on data
-        reho_spheres_output = reho_spheres_marker.fit_transform(
+            output = marker.fit_transform(element_data)
-            {
-                "BOLD": {
-                    "path": "/tmp",
-                    "data": fmri_img,
-                    "meta": {},
-                    "space": "MNI",
-                }
-            }
-        )
-        # Get BOLD output
-        reho_spheres_output_bold = reho_spheres_output["BOLD"]
-        # Assert BOLD output keys
-        assert "data" in reho_spheres_output_bold
-        assert "col_names" in reho_spheres_output_bold
 
-        reho_spheres_output_bold_data = reho_spheres_output_bold["data"]
+            assert "Creating cache" in caplog.text
-        # Assert BOLD output data dimension
+
-        assert reho_spheres_output_bold_data.ndim == 2
+            # Get BOLD output
-        # Assert BOLD output data is normalized
+            assert "BOLD" in output
-        assert (reho_spheres_output_bold_data > 0).all() and (
+            output_bold = output["BOLD"]
-            reho_spheres_output_bold_data < 1
+            # Assert BOLD output keys
-        ).all()
+            assert "data" in output_bold
+            assert "col_names" in output_bold
+
+            output_bold_data = output_bold["data"]
+            # Assert BOLD output data dimension
+            assert output_bold_data.ndim == 2
+            # Assert BOLD output data is normalized
+            assert (output_bold_data > 0).all() and (
+                output_bold_data < 1
+            ).all()
+
+            # Reset log capture
+            caplog.clear()
+            # Initialize storage
+            storage = SQLiteFeatureStorage(tmp_path / "reho_spheres.sqlite")
+            # Fit transform marker on data with storage
+            marker.fit_transform(
+                input=element_data,
+                storage=storage,
+            )
+            # Cache working correctly
+            assert "Creating cache" not in caplog.text
 
 
 @pytest.mark.skipif(
-    _check_afni() is False, reason="requires afni to be in PATH"
+    _check_afni() is False, reason="requires AFNI to be in PATH"
 )
-def test_reho_spheres_computation_comparison(tmp_path: Path) -> None:
+def test_ReHoSpheres_comparison(tmp_path: Path) -> None:
-    """Test ReHoSpheres fit-transform implementation comparison.
+    """Test ReHoSpheres implementation comparison.
 
     Parameters
     ----------
@@ -76,92 +86,35 @@ def test_reho_spheres_computation_comparison(tmp_path: Path) -> None:
 
     """
     with SPMAuditoryTestingDataGrabber() as dg:
-        # Use first subject
+        element_data = DefaultDataReader().fit_transform(dg["sub001"])
-        subject_data = dg["sub001"]
-        # Load image to memory
-        fmri_img = nimg.load_img(subject_data["BOLD"]["path"])
-
         # Update workdir to current test's tmp_path
         WorkDirManager().workdir = tmp_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": {},
-                    "space": "MNI",
-                }
-            }
-        )
-        # Get BOLD output
-        reho_spheres_output_bold_python = reho_spheres_output_python["BOLD"]
 
-        # Initialize marker with use_afni=True
+        # Initialize marker
-        reho_spheres_marker_afni = ReHoSpheres(
+        junifer_marker = ReHoSpheres(
-            coords=COORDINATES, radius=10.0, use_afni=True
+            coords=COORDINATES,
+            using="junifer",
+            radius=10.0,
         )
         # Fit transform marker on data
-        reho_spheres_output_afni = reho_spheres_marker_afni.fit_transform(
+        junifer_output = junifer_marker.fit_transform(element_data)
-            {
-                "BOLD": {
-                    "path": "/tmp",
-                    "data": fmri_img,
-                    "meta": {},
-                    "space": "MNI",
-                }
-            }
-        )
         # Get BOLD output
-        reho_spheres_output_bold_afni = reho_spheres_output_afni["BOLD"]
+        junifer_output_bold = junifer_output["BOLD"]
+
+        # Initialize marker
+        afni_marker = ReHoSpheres(
+            coords=COORDINATES,
+            using="afni",
+            radius=10.0,
+        )
+        # Fit transform marker on data
+        afni_output = afni_marker.fit_transform(element_data)
+        # Get BOLD output
+        afni_output_bold = afni_output["BOLD"]
 
         # Check for Pearson correlation coefficient
-        r, _ = pearsonr(
+        r, _ = sp.stats.pearsonr(
-            reho_spheres_output_bold_python["data"].flatten(),
+            junifer_output_bold["data"].flatten(),
-            reho_spheres_output_bold_afni["data"].flatten(),
+            afni_output_bold["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"])
-        # Update workdir to current test's tmp_path
-        WorkDirManager().workdir = tmp_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,
-                    "space": "MNI",
-                }
-            },
-            storage=reho_spheres_storage,
-        )

junifer/pipeline/pipeline_step_mixin.py

@@ -44,7 +44,7 @@ class PipelineStepMixin:
         raise_error(
             msg="Concrete classes need to implement validate_input().",
             klass=NotImplementedError,
-        )
+        )  # pragma: no cover
 
     def get_output_type(self, input_type: str) -> str:
         """Get output type.
@@ -63,10 +63,11 @@ class PipelineStepMixin:
         raise_error(
             msg="Concrete classes need to implement get_output_type().",
             klass=NotImplementedError,
-        )
+        )  # pragma: no cover
 
     def _fit_transform(
-        self, input: Dict[str, Dict], **kwargs: Any
+        self,
+        input: Dict[str, Dict],
     ) -> Dict[str, Dict]:
         """Fit and transform.
 
@@ -74,8 +75,6 @@ class PipelineStepMixin:
         ----------
         input : dict
             The Junifer Data object.
-        **kwargs : dict
-            Extra keyword arguments.
 
         Returns
         -------
@@ -86,7 +85,7 @@ class PipelineStepMixin:
         raise_error(
             msg="Concrete classes need to implement _fit_transform().",
             klass=NotImplementedError,
-        )
+        )  # pragma: no cover
 
     def validate(self, input: List[str]) -> List[str]:
         """Validate the the pipeline step.
@@ -101,44 +100,109 @@ class PipelineStepMixin:
         list of str
             The output of the pipeline step.
 
-        Raises
-        ------
-        ValueError
-            If the pipeline step object is missing dependencies required for
-            its working or if the input does not have the required data.
-
         """
-        # Check if _DEPENDENCIES attribute is found;
-        # (markers and preprocessors will have them but not datareaders
-        # as of now)
-        dependencies_not_found = []
-        if hasattr(self, "_DEPENDENCIES"):
-            # Check if dependencies are importable
-            for dependency in self._DEPENDENCIES:  # type: ignore
-                # First perform an easy check
-                if find_spec(dependency) is None:
-                    # Then check mapped names
-                    if dependency not in list(
-                        chain.from_iterable(packages_distributions().values())
-                    ):
-                        dependencies_not_found.append(dependency)
-        # Raise error if any dependency is not found
-        if dependencies_not_found:
-            raise_error(
-                msg=f"{dependencies_not_found} are not installed but are "
-                "required for using {self.name}.",
-                klass=ImportError,
-            )
-        # Check if _EXT_DEPENDENCIES attribute is found;
-        # (some markers might have them like ReHo-family)
-        if hasattr(self, "_EXT_DEPENDENCIES"):
-            for dependency in self._EXT_DEPENDENCIES:  # type: ignore
-                out = check_ext_dependencies(**dependency)
-                if getattr(self, f"use_{dependency['name']}", None) is None:
-                    # Set attribute for using external tools
-                    setattr(self, f"use_{dependency['name']}", out)
 
+        def _check_dependencies(obj) -> None:
+            """Check obj._DEPENDENCIES.
+
+            Parameters
+            ----------
+            obj : object
+                Object to check _DEPENDENCIES of.
+
+            Raises
+            ------
+            ImportError
+                If the pipeline step object is missing dependencies required
+                for its working.
+
+            """
+            # Check if _DEPENDENCIES attribute is found;
+            # (markers and preprocessors will have them but not datareaders
+            # as of now)
+            dependencies_not_found = []
+            if hasattr(obj, "_DEPENDENCIES"):
+                # Check if dependencies are importable
+                for dependency in obj._DEPENDENCIES:
+                    # First perform an easy check
+                    if find_spec(dependency) is None:
+                        # Then check mapped names
+                        if dependency not in list(
+                            chain.from_iterable(
+                                packages_distributions().values()
+                            )
+                        ):
+                            dependencies_not_found.append(dependency)
+            # Raise error if any dependency is not found
+            if dependencies_not_found:
+                raise_error(
+                    msg=(
+                        f"{dependencies_not_found} are not installed but are "
+                        f"required for using {obj.__class__.__name__}."
+                    ),
+                    klass=ImportError,
+                )
+
+        def _check_ext_dependencies(obj) -> None:
+            """Check obj._EXT_DEPENDENCIES.
+
+            Parameters
+            ----------
+            obj : object
+                Object to check _EXT_DEPENDENCIES of.
+
+            """
+            # Check if _EXT_DEPENDENCIES attribute is found;
+            # (some markers and preprocessors might have them)
+            if hasattr(obj, "_EXT_DEPENDENCIES"):
+                for dependency in obj._EXT_DEPENDENCIES:
+                    check_ext_dependencies(**dependency)
+
+        def _check_conditional_dependencies(obj) -> None:
+            """Check obj._CONDITIONAL_DEPENDENCIES.
+
+            Parameters
+            ----------
+            obj : object
+                Object to check _CONDITIONAL_DEPENDENCIES of.
+
+            Raises
+            ------
+            AttributeError
+                If the pipeline step object does not have `using` as a
+                constructor parameter.
+
+            """
+            # Check if _CONDITIONAL_DEPENDENCIES attribute is found;
+            # (some markers and preprocessors might have them)
+            if hasattr(obj, "_CONDITIONAL_DEPENDENCIES"):
+                if not hasattr(obj, "using"):
+                    raise_error(
+                        msg=(
+                            f"The pipeline step: {obj.__class__.__name__} has "
+                            "`_CONDITIONAL_DEPENDENCIES` but does not have "
+                            "`using` as a constructor parameter"
+                        ),
+                        klass=AttributeError,
+                    )
+                else:
+                    for dependency in obj._CONDITIONAL_DEPENDENCIES:
+                        if dependency["using"] == obj.using:
+                            depends_on = dependency["depends_on"]
+                            # Check dependencies
+                            _check_dependencies(depends_on)
+                            # Check external dependencies
+                            _check_ext_dependencies(depends_on)
+
+        # Check dependencies
+        _check_dependencies(self)
+        # Check external dependencies
+        _check_ext_dependencies(self)
+        # Check conditional dependencies
+        _check_conditional_dependencies(self)
+        # Validate input
         fit_input = self.validate_input(input=input)
+        # Validate output type
         outputs = [self.get_output_type(t_input) for t_input in fit_input]
         return outputs
 

junifer/pipeline/tests/test_pipeline_step_mixin.py

@@ -1,11 +1,11 @@
-"""Provide tests for pipeline mixin."""
+"""Provide tests for PipelineStepMixin."""
 
 # Authors: Federico Raimondo <f.raimondo@fz-juelich.de>
 #          Synchon Mandal <s.mandal@fz-juelich.de>
 # License: AGPL
 
 import warnings
-from typing import ClassVar, Dict, List, Set, Union
+from typing import ClassVar, Dict, List, Set, Type, Union
 
 import pytest
 
@@ -13,24 +13,13 @@ from junifer.pipeline.pipeline_step_mixin import PipelineStepMixin
 from junifer.pipeline.utils import _check_afni
 
 
-def test_PipelineStepMixin() -> None:
+def test_PipelineStepMixin_correct_dependencies() -> None:
-    """Test PipelineStepMixin."""
+    """Test fit-transform with correct dependencies."""
-    mixin = PipelineStepMixin()
-    with pytest.raises(NotImplementedError):
-        mixin.validate_input([])
-    with pytest.raises(NotImplementedError):
-        mixin.get_output_type("")
-    with pytest.raises(NotImplementedError):
-        mixin._fit_transform({})
-
-
-def test_pipeline_step_mixin_validate_correct_dependencies() -> None:
-    """Test validate with correct dependencies."""
 
     class CorrectMixer(PipelineStepMixin):
         """Test class for validation."""
 
-        _DEPENDENCIES: ClassVar[Set[str]] = {"setuptools"}
+        _DEPENDENCIES: ClassVar[Set[str]] = {"math"}
 
         def validate_input(self, input: List[str]) -> List[str]:
             return input
@@ -38,15 +27,15 @@ def test_pipeline_step_mixin_validate_correct_dependencies() -> None:
         def get_output_type(self, input_type: str) -> str:
             return input_type
 
-        def fit_transform(self, input: Dict[str, Dict]) -> Dict[str, Dict]:
+        def _fit_transform(self, input: Dict[str, Dict]) -> Dict[str, Dict]:
             return {"input": input}
 
     mixer = CorrectMixer()
-    mixer.validate([])
+    mixer.fit_transform({})
 
 
-def test_pipeline_step_mixin_validate_incorrect_dependencies() -> None:
+def test_PipelineStepMixin_incorrect_dependencies() -> None:
-    """Test validate with incorrect dependencies."""
+    """Test fit-transform with incorrect dependencies."""
 
     class IncorrectMixer(PipelineStepMixin):
         """Test class for validation."""
@@ -59,25 +48,49 @@ def test_pipeline_step_mixin_validate_incorrect_dependencies() -> None:
         def get_output_type(self, input_type: str) -> str:
             return input_type
 
-        def fit_transform(self, input: Dict[str, Dict]) -> Dict[str, Dict]:
+        def _fit_transform(self, input: Dict[str, Dict]) -> Dict[str, Dict]:
             return {"input": input}
 
     mixer = IncorrectMixer()
     with pytest.raises(ImportError, match="not installed"):
-        mixer.validate([])
+        mixer.fit_transform({})
 
 
 @pytest.mark.skipif(
-    _check_afni() is False, reason="requires afni to be in PATH"
+    _check_afni() is False, reason="requires AFNI to be in PATH"
 )
-def test_pipeline_step_mixin_validate_correct_ext_dependencies() -> None:
+def test_PipelineStepMixin_correct_ext_dependencies() -> None:
-    """Test validate with correct external dependencies."""
+    """Test fit-transform with correct external dependencies."""
 
     class CorrectMixer(PipelineStepMixin):
         """Test class for validation."""
 
-        _EXT_DEPENDENCIES: ClassVar[List[Dict[str, Union[str, bool]]]] = [
+        _EXT_DEPENDENCIES: ClassVar[List[Dict[str, str]]] = [{"name": "afni"}]
-            {"name": "afni", "optional": False}
+
+        def validate_input(self, input: List[str]) -> List[str]:
+            return input
+
+        def get_output_type(self, input_type: str) -> str:
+            return input_type
+
+        def _fit_transform(self, input: Dict[str, Dict]) -> Dict[str, Dict]:
+            return {"input": input}
+
+    mixer = CorrectMixer()
+    mixer.fit_transform({})
+
+
+@pytest.mark.skipif(
+    _check_afni() is False, reason="requires AFNI to be in PATH"
+)
+def test_PipelineStepMixin_ext_deps_correct_commands() -> None:
+    """Test fit-transform with correct external dependency commands."""
+
+    class CorrectMixer(PipelineStepMixin):
+        """Test class for validation."""
+
+        _EXT_DEPENDENCIES: ClassVar[List[Dict[str, Union[str, List[str]]]]] = [
+            {"name": "afni", "commands": ["3dReHo"]}
         ]
 
         def validate_input(self, input: List[str]) -> List[str]:
@@ -86,53 +99,27 @@ def test_pipeline_step_mixin_validate_correct_ext_dependencies() -> None:
         def get_output_type(self, input_type: str) -> str:
             return input_type
 
-        def fit_transform(self, input: Dict[str, Dict]) -> Dict[str, Dict]:
+        def _fit_transform(self, input: Dict[str, Dict]) -> Dict[str, Dict]:
-            return {"input": input}
-
-    mixer = CorrectMixer()
-    mixer.validate([])
-
-
-@pytest.mark.skipif(
-    _check_afni() is False, reason="requires afni to be in PATH"
-)
-def test_pipeline_step_mixin_validate_ext_deps_correct_commands() -> None:
-    """Test validate with correct external dependencies' correct commands."""
-
-    class CorrectMixer(PipelineStepMixin):
-        """Test class for validation."""
-
-        _EXT_DEPENDENCIES: ClassVar[
-            List[Dict[str, Union[str, bool, List[str]]]]
-        ] = [{"name": "afni", "optional": False, "commands": ["3dReHo"]}]
-
-        def validate_input(self, input: List[str]) -> List[str]:
-            return input
-
-        def get_output_type(self, input_type: str) -> str:
-            return input_type
-
-        def fit_transform(self, input: Dict[str, Dict]) -> Dict[str, Dict]:
             return {"input": input}
 
     mixer = CorrectMixer()
     with warnings.catch_warnings():
         warnings.simplefilter("error")
-        mixer.validate([])
+        mixer.fit_transform({})
 
 
 @pytest.mark.skipif(
-    _check_afni() is False, reason="requires afni to be in PATH"
+    _check_afni() is False, reason="requires AFNI to be in PATH"
 )
-def test_pipeline_step_mixin_validate_ext_deps_incorrect_commands() -> None:
+def test_PipelineStepMixin_ext_deps_incorrect_commands() -> None:
-    """Test validate with correct external dependencies' incorrect commands."""
+    """Test fit-transform with inccorrect external dependency commands."""
 
     class CorrectMixer(PipelineStepMixin):
         """Test class for validation."""
 
-        _EXT_DEPENDENCIES: ClassVar[
+        _EXT_DEPENDENCIES: ClassVar[List[Dict[str, Union[str, List[str]]]]] = [
-            List[Dict[str, Union[str, bool, List[str]]]]
+            {"name": "afni", "commands": ["3d"]}
-        ] = [{"name": "afni", "optional": False, "commands": ["3d"]}]
+        ]
 
         def validate_input(self, input: List[str]) -> List[str]:
             return input
@@ -140,16 +127,16 @@ def test_pipeline_step_mixin_validate_ext_deps_incorrect_commands() -> None:
         def get_output_type(self, input_type: str) -> str:
             return input_type
 
-        def fit_transform(self, input: Dict[str, Dict]) -> Dict[str, Dict]:
+        def _fit_transform(self, input: Dict[str, Dict]) -> Dict[str, Dict]:
             return {"input": input}
 
     mixer = CorrectMixer()
     with pytest.warns(RuntimeWarning, match="AFNI is installed"):
-        mixer.validate([])
+        mixer.fit_transform({})
 
 
-def test_pipeline_step_mixin_validate_incorrect_ext_dependencies() -> None:
+def test_PipelineStepMixin_incorrect_ext_dependencies() -> None:
-    """Test validate with incorrect external dependencies."""
+    """Test fit-transform with incorrect external dependencies."""
 
     class IncorrectMixer(PipelineStepMixin):
         """Test class for validation."""
@@ -164,9 +151,110 @@ def test_pipeline_step_mixin_validate_incorrect_ext_dependencies() -> None:
         def get_output_type(self, input_type: str) -> str:
             return input_type
 
-        def fit_transform(self, input: Dict[str, Dict]) -> Dict[str, Dict]:
+        def _fit_transform(self, input: Dict[str, Dict]) -> Dict[str, Dict]:
             return {"input": input}
 
     mixer = IncorrectMixer()
-    with pytest.raises(ValueError, match="too adventurous"):
+    with pytest.raises(ValueError, match="Invalid value"):
-        mixer.validate([])
+        mixer.fit_transform({})
+
+
+def test_PipelineStepMixin_correct_conditional_dependencies() -> None:
+    """Test fit-transform with correct conditional dependencies."""
+
+    class Dependency:
+        _DEPENDENCIES: ClassVar[Set[str]] = {"math"}
+
+    class CorrectMixer(PipelineStepMixin):
+        """Test class for validation."""
+
+        _CONDITIONAL_DEPENDENCIES: ClassVar[
+            List[Dict[str, Union[str, Type]]]
+        ] = [
+            {
+                "using": "math",
+                "depends_on": Dependency,
+            },
+        ]
+
+        using = "math"
+
+        def validate_input(self, input: List[str]) -> List[str]:
+            return input
+
+        def get_output_type(self, input_type: str) -> str:
+            return input_type
+
+        def _fit_transform(self, input: Dict[str, Dict]) -> Dict[str, Dict]:
+            return {"input": input}
+
+    mixer = CorrectMixer()
+    mixer.fit_transform({})
+
+
+def test_PipelineStepMixin_incorrect_conditional_dependencies() -> None:
+    """Test fit-transform with incorrect conditional dependencies."""
+
+    class Dependency:
+        _DEPENDENCIES: ClassVar[Set[str]] = {"math"}
+
+    class IncorrectMixer(PipelineStepMixin):
+        """Test class for validation."""
+
+        _CONDITIONAL_DEPENDENCIES: ClassVar[
+            List[Dict[str, Union[str, Type]]]
+        ] = [
+            {
+                "using": "math",
+                "depends_on": Dependency,
+            },
+        ]
+
+        def validate_input(self, input: List[str]) -> List[str]:
+            return input
+
+        def get_output_type(self, input_type: str) -> str:
+            return input_type
+
+        def _fit_transform(self, input: Dict[str, Dict]) -> Dict[str, Dict]:
+            return {"input": input}
+
+    mixer = IncorrectMixer()
+    with pytest.raises(AttributeError, match="`using`"):
+        mixer.fit_transform({})
+
+
+@pytest.mark.skipif(
+    _check_afni() is False, reason="requires AFNI to be in PATH"
+)
+def test_PipelineStepMixin_correct_conditional_ext_dependencies() -> None:
+    """Test fit-transform with correct conditional external dependencies."""
+
+    class ExternalDependency:
+        _EXT_DEPENDENCIES: ClassVar[List[Dict[str, str]]] = [{"name": "afni"}]
+
+    class CorrectMixer(PipelineStepMixin):
+        """Test class for validation."""
+
+        _CONDITIONAL_DEPENDENCIES: ClassVar[
+            List[Dict[str, Union[str, Type]]]
+        ] = [
+            {
+                "using": "afni",
+                "depends_on": ExternalDependency,
+            },
+        ]
+
+        using = "afni"
+
+        def validate_input(self, input: List[str]) -> List[str]:
+            return input
+
+        def get_output_type(self, input_type: str) -> str:
+            return input_type
+
+        def _fit_transform(self, input: Dict[str, Dict]) -> Dict[str, Dict]:
+            return {"input": input}
+
+    mixer = CorrectMixer()
+    mixer.fit_transform({})

junifer/pipeline/utils.py

@@ -10,16 +10,17 @@ from typing import Any, List, Optional
 from junifer.utils.logging import raise_error, warn_with_log
 
 
-def check_ext_dependencies(name: str, optional: bool, **kwargs: Any) -> bool:
+def check_ext_dependencies(
+    name: str, optional: bool = False, **kwargs: Any
+) -> bool:
     """Check if external dependency `name` is found if mandatory.
 
     Parameters
     ----------
     name : str
         The name of the dependency.
-    optional : bool
+    optional : bool, optional
-        Whether the dependency is optional. For external dependencies marked
+        Whether the dependency is optional (default False).
-        as optional, there should be an implementation provided with junfier.
     **kwargs : dict
         Extra keyword arguments.
 
@@ -28,7 +29,20 @@ def check_ext_dependencies(name: str, optional: bool, **kwargs: Any) -> bool:
     bool
         Whether the external dependency was found.
 
+    Raises
+    ------
+    ValueError
+        If ``name`` is invalid.
+    RuntimeError
+        If ``name`` is mandatory and is not found.
+
     """
+    valid_ext_dependencies = ("afni", "fsl", "ants")
+    if name not in valid_ext_dependencies:
+        raise_error(
+            "Invalid value for `name`, should be one of: "
+            f"{valid_ext_dependencies}"
+        )
     # Check for afni
     if name == "afni":
         found = _check_afni(**kwargs)
@@ -38,29 +52,26 @@ def check_ext_dependencies(name: str, optional: bool, **kwargs: Any) -> bool:
     # Check for ants
     elif name == "ants":
         found = _check_ants(**kwargs)
-    # Went off the rails
+
-    else:
-        raise_error(
-            f"The external dependency {name} has no check. "
-            f"Either the name '{name}' is incorrect or you were too "
-            "adventurous. Raise an issue if it's the latter ;-)."
-        )
     # Check if the dependency is mandatory in case it's not found
     if not found and not optional:
         raise_error(
-            f"{name} is not installed but is "
+            msg=(
-            "required by one of the pipeline steps."
+                f"{name} is not installed but is "
+                "required by one of the pipeline steps"
+            ),
+            klass=RuntimeError,
         )
     return found
 
 
 def _check_afni(commands: Optional[List[str]] = None) -> bool:
-    """Check if afni is present in the system.
+    """Check if AFNI is present in the system.
 
     Parameters
     ----------
     commands : list of str, optional
-        The commands to specifically check for from afni. If None, only
+        The commands to specifically check for from AFNI. If None, only
         the basic afni version would be looked up, else, would also
         check for specific commands (default None).
 
@@ -115,19 +126,19 @@ def _check_afni(commands: Optional[List[str]] = None) -> bool:
 
 
 def _check_fsl(commands: Optional[List[str]] = None) -> bool:
-    """Check if fsl is present in the system.
+    """Check if FSL is present in the system.
 
     Parameters
     ----------
     commands : list of str, optional
-        The commands to specifically check for from fsl. If None, only
+        The commands to specifically check for from FSL. If None, only
-        the basic fsl flirt version would be looked up, else, would also
+        the basic FSL flirt version would be looked up, else, would also
         check for specific commands (default None).
 
     Returns
     -------
     bool
-        Whether fsl is found or not.
+        Whether FSL is found or not.
 
     """
     completed_process = subprocess.run(

junifer/preprocess/ants/ants_apply_transforms_warper.py

@@ -20,14 +20,13 @@ import numpy as np
 
 from ...pipeline import WorkDirManager
 from ...utils import logger, raise_error, run_ext_cmd
-from ..base import BasePreprocessor
 
 
 if TYPE_CHECKING:
     from nibabel import Nifti1Image
 
 
-class _AntsApplyTransformsWarper(BasePreprocessor):
+class _AntsApplyTransformsWarper:
     """Class for warping NIfTI images via ANTs antsApplyTransforms.
 
     Warps ANTs ``antsApplyTransforms``.
@@ -46,12 +45,9 @@ class _AntsApplyTransformsWarper(BasePreprocessor):
 
     """
 
-    _EXT_DEPENDENCIES: ClassVar[
+    _EXT_DEPENDENCIES: ClassVar[List[Dict[str, Union[str, List[str]]]]] = [
-        List[Dict[str, Union[str, bool, List[str]]]]
-    ] = [
         {
             "name": "ants",
-            "optional": False,
             "commands": ["ResampleImage", "antsApplyTransforms"],
         },
     ]
@@ -62,42 +58,7 @@ class _AntsApplyTransformsWarper(BasePreprocessor):
         # Check only single data type is passed
         if isinstance(on, list):
             raise_error("Can only work on single data type, list was passed.")
-        self.on = on  # needed for the base validation to work
+        self.on = on
-        super().__init__(
-            on=self.on, required_data_types=[self.on, self.ref, "Warp"]
-        )
-
-    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
-            preprocessor.
-
-        """
-        # Constructed dynamically
-        return [self.on]
-
-    def get_output_type(self, input: List[str]) -> List[str]:
-        """Get output type.
-
-        Parameters
-        ----------
-        input : list of str
-            The input to the preprocessor. The list must contain the
-            available Junifer Data dictionary keys.
-
-        Returns
-        -------
-        list of str
-            The updated list of available Junifer Data object keys after
-            the pipeline step.
-
-        """
-        # Does not add any new keys
-        return input
 
     def _run_apply_transforms(
         self,

junifer/preprocess/ants/tests/test_ants_apply_transforms_warper.py

@@ -4,10 +4,7 @@
 # License: AGPL
 
 import socket
-from pathlib import Path
-from typing import List
 
-import nibabel as nib
 import pytest
 
 from junifer.datagrabber import DMCC13Benchmark
@@ -25,71 +22,6 @@ def test_AntsApplyTransformsWarper_init() -> None:
     )
     assert ants_apply_transforms_warper.ref == "T1w"
     assert ants_apply_transforms_warper.on == "BOLD"
-    assert ants_apply_transforms_warper._on == ["BOLD"]
-
-
-def test_AntsApplyTransformsWarper_get_valid_inputs() -> None:
-    """Test AntsApplyTransformsWarper get_valid_inputs."""
-    ants_apply_transforms_warper = _AntsApplyTransformsWarper(
-        reference="T1w", on="BOLD"
-    )
-    assert ants_apply_transforms_warper.get_valid_inputs() == ["BOLD"]
-
-
-@pytest.mark.parametrize(
-    "input_",
-    [
-        ["BOLD", "T1w", "Warp"],
-        ["BOLD", "T1w"],
-        ["BOLD"],
-    ],
-)
-def test_AntsApplyTransformsWarper_get_output_type(input_: List[str]) -> None:
-    """Test AntsApplyTransformsWarper get_output_type.
-
-    Parameters
-    ----------
-    input_ : list of str
-        The input data types.
-
-    """
-    ants_apply_transforms_warper = _AntsApplyTransformsWarper(
-        reference="T1w", on="BOLD"
-    )
-    assert ants_apply_transforms_warper.get_output_type(input_) == input_
-
-
-@pytest.mark.skipif(
-    _check_ants() is False, reason="requires ANTs to be in PATH"
-)
-@pytest.mark.skipif(
-    socket.gethostname() != "juseless",
-    reason="only for juseless",
-)
-def test_AntsApplyTransformsWarper__run_apply_transform() -> None:
-    """Test AntsApplyTransformsWarper _run_apply_transform."""
-    with DMCC13Benchmark(
-        types=["BOLD", "T1w", "Warp"],
-        sessions=["wave1bas"],
-        tasks=["Rest"],
-        phase_encodings=["AP"],
-        runs=["1"],
-        native_t1w=True,
-    ) as dg:
-        # Read data
-        element_data = DefaultDataReader().fit_transform(
-            dg[("f9057kp", "wave1bas", "Rest", "AP", "1")]
-        )
-        # Preprocess data
-        warped_data, resampled_ref_path = _AntsApplyTransformsWarper(
-            reference="T1w", on="BOLD"
-        )._run_apply_transforms(
-            input_data=element_data["BOLD"],
-            ref_path=element_data["T1w"]["path"],
-            warp_path=element_data["Warp"]["path"],
-        )
-        assert isinstance(warped_data, nib.Nifti1Image)
-        assert isinstance(resampled_ref_path, Path)
 
 
 @pytest.mark.skipif(

junifer/preprocess/bold_warper.py

@@ -10,6 +10,7 @@ from typing import (
     List,
     Optional,
     Tuple,
+    Type,
     Union,
 )
 
@@ -31,6 +32,12 @@ class BOLDWarper(BasePreprocessor):
 
     Parameters
     ----------
+    using : {"fsl", "ants"}
+        Implementation to use for warping:
+
+        * "fsl" : Use FSL's ``applywarp``
+        * "afni" : Use ANTs' ``antsApplyTransforms``
+
     reference : str
         The data type to use as reference for warping, can be either a data
         type like "T1w" or a template space like "MNI152NLin2009cAsym".
@@ -38,27 +45,39 @@ class BOLDWarper(BasePreprocessor):
     Raises
     ------
     ValueError
-        If ``reference`` is invalid.
+        If ``using`` is invalid or
+        if ``reference`` is invalid.
+
+    Notes
+    -----
+    If you are setting ``reference`` to a template space like
+    "MNI152NLin2009cAsym", make sure ANTs is available for the
+    transformation else it will fail during runtime. It is tricky to validate
+    this beforehand and difficult to enforce this as a requirement, hence the
+    heads-up.
 
     """
 
-    _EXT_DEPENDENCIES: ClassVar[
+    _CONDITIONAL_DEPENDENCIES: ClassVar[List[Dict[str, Union[str, Type]]]] = [
-        List[Dict[str, Union[str, bool, List[str]]]]
-    ] = [
         {
-            "name": "fsl",
+            "using": "fsl",
-            "optional": True,
+            "depends_on": _ApplyWarper,
-            "commands": ["flirt", "applywarp"],
         },
         {
-            "name": "ants",
+            "using": "ants",
-            "optional": True,
+            "depends_on": _AntsApplyTransformsWarper,
-            "commands": ["ResampleImage", "antsApplyTransforms"],
         },
     ]
 
-    def __init__(self, reference: str) -> None:
+    def __init__(self, using: str, reference: str) -> None:
         """Initialize the class."""
+        # Validate `using` parameter
+        valid_using = [dep["using"] for dep in self._CONDITIONAL_DEPENDENCIES]
+        if using not in valid_using:
+            raise_error(
+                f"Invalid value for `using`, should be one of: {valid_using}"
+            )
+        self.using = using
         self.ref = reference
         # Initialize superclass based on reference
         if self.ref == "T1w":

junifer/preprocess/fsl/apply_warper.py

@@ -20,14 +20,13 @@ import numpy as np
 
 from ...pipeline import WorkDirManager
 from ...utils import logger, raise_error, run_ext_cmd
-from ..base import BasePreprocessor
 
 
 if TYPE_CHECKING:
     from nibabel import Nifti1Image
 
 
-class _ApplyWarper(BasePreprocessor):
+class _ApplyWarper:
     """Class for warping NIfTI images via FSL FLIRT.
 
     Wraps FSL FLIRT ``applywarp``.
@@ -46,12 +45,9 @@ class _ApplyWarper(BasePreprocessor):
 
     """
 
-    _EXT_DEPENDENCIES: ClassVar[
+    _EXT_DEPENDENCIES: ClassVar[List[Dict[str, Union[str, List[str]]]]] = [
-        List[Dict[str, Union[str, bool, List[str]]]]
-    ] = [
         {
             "name": "fsl",
-            "optional": False,
             "commands": ["flirt", "applywarp"],
         },
     ]
@@ -62,42 +58,7 @@ class _ApplyWarper(BasePreprocessor):
         # Check only single data type is passed
         if isinstance(on, list):
             raise_error("Can only work on single data type, list was passed.")
-        self.on = on  # needed for the base validation to work
+        self.on = on
-        super().__init__(
-            on=self.on, required_data_types=[self.on, self.ref, "Warp"]
-        )
-
-    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
-            preprocessor.
-
-        """
-        # Constructed dynamically
-        return [self.on]
-
-    def get_output_type(self, input: List[str]) -> List[str]:
-        """Get output type.
-
-        Parameters
-        ----------
-        input : list of str
-            The input to the preprocessor. The list must contain the
-            available Junifer Data dictionary keys.
-
-        Returns
-        -------
-        list of str
-            The updated list of available Junifer Data object keys after
-            the pipeline step.
-
-        """
-        # Does not add any new keys
-        return input
 
     def _run_applywarp(
         self,

junifer/preprocess/fsl/tests/test_apply_warper.py

@@ -4,10 +4,7 @@
 # License: AGPL
 
 import socket
-from pathlib import Path
-from typing import List
 
-import nibabel as nib
 import pytest
 
 from junifer.datagrabber import DataladHCP1200
@@ -21,62 +18,6 @@ def test_ApplyWarper_init() -> None:
     apply_warper = _ApplyWarper(reference="T1w", on="BOLD")
     assert apply_warper.ref == "T1w"
     assert apply_warper.on == "BOLD"
-    assert apply_warper._on == ["BOLD"]
-
-
-def test_ApplyWarper_get_valid_inputs() -> None:
-    """Test ApplyWarper get_valid_inputs."""
-    apply_warper = _ApplyWarper(reference="T1w", on="BOLD")
-    assert apply_warper.get_valid_inputs() == ["BOLD"]
-
-
-@pytest.mark.parametrize(
-    "input_",
-    [
-        ["BOLD", "T1w", "Warp"],
-        ["BOLD", "T1w"],
-        ["BOLD"],
-    ],
-)
-def test_ApplyWarper_get_output_type(input_: List[str]) -> None:
-    """Test ApplyWarper get_output_type.
-
-    Parameters
-    ----------
-    input_ : list of str
-        The input data types.
-
-    """
-    apply_warper = _ApplyWarper(reference="T1w", on="BOLD")
-    assert apply_warper.get_output_type(input_) == input_
-
-
-@pytest.mark.skipif(_check_fsl() is False, reason="requires FSL to be in PATH")
-@pytest.mark.skipif(
-    socket.gethostname() != "juseless",
-    reason="only for juseless",
-)
-def test_ApplyWarper__run_applywarp() -> None:
-    """Test ApplyWarper _run_applywarp."""
-    with DataladHCP1200(
-        tasks=["REST1"],
-        phase_encodings=["LR"],
-        ica_fix=True,
-    ) as dg:
-        # Read data
-        element_data = DefaultDataReader().fit_transform(
-            dg[("100206", "REST1", "LR")]
-        )
-        # Preprocess data
-        warped_data, resampled_ref_path = _ApplyWarper(
-            reference="T1w", on="BOLD"
-        )._run_applywarp(
-            input_data=element_data["BOLD"],
-            ref_path=element_data["T1w"]["path"],
-            warp_path=element_data["Warp"]["path"],
-        )
-        assert isinstance(warped_data, nib.Nifti1Image)
-        assert isinstance(resampled_ref_path, Path)
 
 
 @pytest.mark.skipif(_check_fsl() is False, reason="requires FSL to be in PATH")

junifer/preprocess/tests/test_bold_warper.py

@@ -21,19 +21,19 @@ if TYPE_CHECKING:
 
 def test_BOLDWarper_init() -> None:
     """Test BOLDWarper init."""
-    bold_warper = BOLDWarper(reference="T1w")
+    bold_warper = BOLDWarper(using="ants", reference="T1w")
     assert bold_warper._on == ["BOLD"]
 
 
 def test_BOLDWarper_get_valid_inputs() -> None:
     """Test BOLDWarper get_valid_inputs."""
-    bold_warper = BOLDWarper(reference="T1w")
+    bold_warper = BOLDWarper(using="ants", reference="T1w")
     assert bold_warper.get_valid_inputs() == ["BOLD"]
 
 
 def test_BOLDWarper_get_output_type() -> None:
     """Test BOLDWarper get_output_type."""
-    bold_warper = BOLDWarper(reference="T1w")
+    bold_warper = BOLDWarper(using="ants", reference="T1w")
     assert bold_warper.get_output_type("BOLD") == "BOLD"