docs/builtin.rst

@@ -174,7 +174,14 @@ Available
      - Calculate regional homogeneity over spheres placed on coordinates
      - Done
      - 0.0.1
-
+   * - :class:`junifer.markers.AmplitudeLowFrequencyFluctuationParcels`
+     - Calculate (f)ALFF and aggregate using parcellations
+     - Done
+     - 0.0.1
+   * - :class:`junifer.markers.AmplitudeLowFrequencyFluctuationSpheres`
+     - Calculate (f)ALFF and aggregate using spheres placed on coordinates
+     - Done
+     - 0.0.1
 
 Planned
 ~~~~~~~
@@ -189,9 +196,6 @@ Planned
    * - Connectedness
      - Compute connectedness
      - :gh:`34`
-   * - ALFF and (f)ALFF
-     - Detect amplitude of low-frequency fluctuation (ALFF) for resting-state fMRI
-     - :gh:`35`
    * - Permutation entropy, Range entropy, Multiscale entropy and Hurst exponent
      - Calculate Permutation entropy, Range entropy, Multiscale entropy and Hurst exponent
      - :gh:`61`

docs/changes/latest.inc

@@ -98,6 +98,8 @@ Enhancements
 
 - Implement :class:`junifer.markers.ReHoParcels` and :class:`junifer.markers.ReHoSpheres` markers (:gh:`36` by `Synchon Mandal`_).
 
+- Implement :class:`junifer.markers.AmplitudeLowFrequencyFluctuationParcels` and :class:`junifer.markers.AmplitudeLowFrequencyFluctuationSpheres` markers (:gh:`35` by `Fede Raimondo`_).
+
 Bugs
 ~~~~
 

docs/conf.py

@@ -155,6 +155,8 @@ numpydoc_xref_ignore = {
     "class",
     "objects",
     "Engine",
+    "positive",
+    "negative",
 }
 # numpydoc_validation_checks = {
 #     "all",

junifer/api/res/afni/3dRSFC

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

junifer/markers/__init__.py

@@ -13,3 +13,7 @@ from .functional_connectivity_spheres import FunctionalConnectivitySpheres
 from .parcel_aggregation import ParcelAggregation
 from .sphere_aggregation import SphereAggregation
 from .reho import ReHoParcels, ReHoSpheres
+from .falff import (
+    AmplitudeLowFrequencyFluctuationParcels,
+    AmplitudeLowFrequencyFluctuationSpheres,
+)

junifer/markers/falff/__init__.py

@@ -0,0 +1,7 @@
+"""Provide imports for falff sub-package."""
+
+# Authors: Federico Raimondo <f.raimondo@fz-juelich.de>
+# License: AGPL
+
+from .falff_parcels import AmplitudeLowFrequencyFluctuationParcels
+from .falff_spheres import AmplitudeLowFrequencyFluctuationSpheres

junifer/markers/falff/falff_base.py

@@ -0,0 +1,180 @@
+"""Provide abstract class for computing fALFF."""
+
+# Authors: Federico Raimondo <f.raimondo@fz-juelich.de>
+#          Amir Omidvarnia <a.omidvarnia@fz-juelich.de>
+#          Kaustubh R. Patil <k.patil@fz-juelich.de>
+# License: AGPL
+
+from typing import Dict, List, Optional
+
+from abc import abstractmethod
+
+from ..base import BaseMarker
+from .falff_estimator import AmplitudeLowFrequencyFluctuationEstimator
+from ...utils.logging import raise_error
+
+
+class AmplitudeLowFrequencyFluctuationBase(BaseMarker):
+    """Base class for (fractional) Amplitude Low Frequency Fluctuation.
+
+    Parameters
+    ----------
+    fractional : bool
+        Whether to compute fractional ALFF.
+    highpass : positive float
+        Highpass cutoff frequency.
+    lowpass : positive float
+        Lowpass cutoff frequency.
+    tr : positive float, optional
+        The Repetition Time of the BOLD data. If None, will extract
+        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
+        computation. If a dataset is correctly preprocessed, the TR should be
+        extracted from the NIFTI without any issue. However, it has been
+        reported that some preprocessed data might not have the correct TR in
+        the NIFTI header.
+
+    """
+
+    _EXT_DEPENDENCIES = [
+        {
+            "name": "afni",
+            "optional": True,
+            "commands": ["3dRSFC", "3dAFNItoNIFTI"],
+        },
+    ]
+
+    def __init__(
+        self,
+        fractional: bool,
+        highpass: float,
+        lowpass: float,
+        tr: Optional[float] = None,
+        use_afni: Optional[bool] = None,
+        name: Optional[str] = None,
+    ) -> None:
+        if highpass <= 0:
+            raise_error("Highpass must be positive")
+        if lowpass <= 0:
+            raise_error("Lowpass must be positive")
+        if highpass >= lowpass:
+            raise_error("Highpass must be lower than lowpass")
+        self.highpass = highpass
+        self.lowpass = lowpass
+        self.tr = tr
+        self.use_afni = use_afni
+        self.fractional = fractional
+
+        # Create a name based on the class name if none is provided
+        if name is None:
+            suffix = "_fractional" if fractional else ""
+            name = f"{self.__class__.__name__}{suffix}"
+        super().__init__(on="BOLD", name=name)
+
+    def get_valid_inputs(self) -> List[str]:
+        """Get valid data types for input.
+
+        Returns
+        -------
+        list of str
+            The list of data types that can be used as input for this marker.
+
+        """
+        return ["BOLD"]
+
+    def get_output_type(self, input_type: str) -> str:
+        """Get output type.
+
+        Parameters
+        ----------
+        input_type : str
+            The data type input to the marker.
+
+        Returns
+        -------
+        str
+            The storage type output by the marker.
+
+        """
+        return "table"
+
+    def compute(
+        self,
+        input: Dict[str, Dict],
+        extra_input: Optional[Dict] = None,
+    ) -> Dict:
+        """Compute.
+
+        Parameters
+        ----------
+        input : dict
+            A single input from the pipeline data object in which to compute
+            the marker.
+        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).
+
+        Returns
+        -------
+        dict
+            The computed result as dictionary. This will be either returned
+            to the user or stored in the storage by calling the store method
+            with this as a parameter. The dictionary has the following keys:
+
+            * ``data`` : the actual computed values as a numpy.ndarray
+            * ``columns`` : the column labels for the computed values as a list
+            * ``row_names`` (if more than one row is present in data): "scan"
+
+        """
+        if self.use_afni is None:
+            raise_error(
+                "Parameter `use_afni` must be set to True or False in order "
+                "to compute this marker. It is currently set to None (default "
+                "behaviour). This is intended to be for auto-detection. In "
+                "order for that to happen, please call the `validate` method "
+                "before calling the `compute` method."
+            )
+
+        estimator = AmplitudeLowFrequencyFluctuationEstimator()
+
+        alff, falff = estimator.fit_transform(
+            use_afni=self.use_afni,
+            input_data=input,
+            highpass=self.highpass,
+            lowpass=self.lowpass,
+            tr=self.tr,
+        )
+        post_data = falff if self.fractional else alff
+
+        post_input = {
+            "data": post_data,
+            "path": None,
+        }
+
+        out = self._postprocess(post_input)
+
+        return out
+
+    @abstractmethod
+    def _postprocess(self, input: Dict) -> Dict:
+        """Postprocess the output of the estimator.
+
+        Parameters
+        ----------
+        input : dict
+            The output of the estimator. It must have the following
+        """
+        raise_error(
+            "_postprocess must be implemented", klass=NotImplementedError
+        )

junifer/markers/falff/falff_estimator.py

@@ -0,0 +1,338 @@
+"""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 typing import TYPE_CHECKING, Any, Dict, Tuple, Union, Optional
+
+import shutil
+import subprocess
+import tempfile
+from functools import lru_cache
+from pathlib import Path
+
+import nibabel as nib
+import numpy as np
+from scipy.fft import fft, fftfreq
+
+from nilearn import image as nimg
+
+from ...utils import logger, raise_error
+from ..utils import singleton
+
+
+if TYPE_CHECKING:
+    from nibabel import Nifti1Image, Nifti2Image
+
+
+@singleton
+class AmplitudeLowFrequencyFluctuationEstimator:
+    """Estimator class for AmplitudeLowFrequencyFluctuationBase.
+
+    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:`junifer.markers.falff.AmplitudeLowFrequencyFluctuationBase`
+    as it serves a specific purpose.
+
+    Parameters
+    ----------
+    use_afni : bool
+        Whether to use afni for computation. If False, will use python.
+
+    """
+
+    def __init__(self) -> None:
+        self._file_path = None
+        # Create temporary directory for intermittent storage of assets during
+        # computation via afni's 3dReHo
+        self.temp_dir_path = Path(tempfile.mkdtemp())
+
+    def __del__(self) -> None:
+        """Cleanup."""
+        print("Cleaning up temporary directory...")
+        # Delete temporary directory and ignore errors for read-only files
+        shutil.rmtree(self.temp_dir_path, ignore_errors=True)
+
+    @staticmethod
+    def _run_afni_cmd(cmd: str) -> None:
+        """Run AFNI command.
+
+        Parameters
+        ----------
+        cmd : str
+            AFNI command to be executed.
+
+        Raises
+        ------
+        RuntimeError
+            If AFNI command fails.
+        """
+        logger.info(f"AFNI command to be executed: {cmd}")
+        # TODO: Figure out how to capture stdout and stderr
+        process = subprocess.run(
+            cmd,
+            stdin=subprocess.DEVNULL,
+            # stdout=subprocess.STDOUT,
+            # stderr=subprocess.STDOUT,
+            shell=True,
+            check=False,
+        )
+        if process.returncode == 0:
+            logger.info(
+                "AFNI command succeeded with the following output: "
+                f"{process.stdout}"
+            )
+        else:
+            raise_error(
+                msg="AFNI command failed with the following error: "
+                f"{process.stdout}",
+                klass=RuntimeError,
+            )
+
+    def _compute_alff_afni(
+        self,
+        data: Union["Nifti1Image", "Nifti2Image"],
+        highpass: float,
+        lowpass: float,
+        tr: Optional[float],
+    ) -> Tuple["Nifti1Image", "Nifti1Image"]:
+        """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
+        -------
+        alff: Niimg-like object
+            ALFF map.
+        falff: Niimg-like object
+            fALFF map.
+
+        Raises
+        ------
+        RuntimeError
+            If the AFNI commands fails due to some issues
+
+        """
+
+        # Save niimg to nii.gz
+        nifti_in_file_path = self.temp_dir_path / "input.nii"
+        nib.save(data, nifti_in_file_path)
+
+        params_suffix = f"_{highpass}_{lowpass}_{tr}"
+        alff_fname = self.temp_dir_path / f"alff{params_suffix}.nii"
+        falff_fname = self.temp_dir_path / 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"
+
+        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 += f"-dt {tr} "
+        self._run_afni_cmd(bp_cmd)
+
+        # Convert afni's output to nifti
+        convert_cmd = (
+            "3dAFNItoNIFTI "
+            f"-prefix {alff_fname.resolve()} "
+            f"{falff_afni_out_path_prefix}_ALFF+tlrc.BRIK "
+        )
+        self._run_afni_cmd(convert_cmd)
+
+        convert_cmd = (
+            "3dAFNItoNIFTI "
+            f"-prefix {falff_fname.resolve()} "
+            f"{falff_afni_out_path_prefix}_fALFF+tlrc.BRIK "
+        )
+        self._run_afni_cmd(convert_cmd)
+
+        # Cleanup intermediate files
+        for fname in self.temp_dir_path.glob("temp_*"):
+            fname.unlink()
+
+        # Load niftis
+        alff_img = nib.load(alff_fname)
+        falff_img = nib.load(falff_fname)
+
+        return alff_img, falff_img
+
+    def _compute_alff_python(
+        self,
+        data: Union["Nifti1Image", "Nifti2Image"],
+        highpass: float,
+        lowpass: float,
+        tr: Optional[float],
+    ) -> Tuple["Nifti1Image", "Nifti1Image"]:
+        """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
+        -------
+        alff: Niimg-like object
+            ALFF map.
+        falff: Niimg-like object
+            fALFF map.
+        """
+        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)
+        return alff_img, falff_img
+
+    @lru_cache(maxsize=None, typed=True)  # noqa: B019
+    def _compute(
+        self,
+        use_afni: bool,
+        data: Union["Nifti1Image", "Nifti2Image"],
+        highpass: float,
+        lowpass: float,
+        tr: Optional[float],
+    ) -> Tuple["Nifti1Image", "Nifti1Image"]:
+        """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
+        -------
+        alff: Niimg-like object
+            ALFF map.
+        falff: Niimg-like object
+            fALFF map.
+        """
+        if use_afni:
+            output = self._compute_alff_afni(
+                data=data,
+                highpass=highpass,
+                lowpass=lowpass,
+                tr=tr,
+            )
+        else:
+            output = self._compute_alff_python(
+                data, highpass=highpass, lowpass=lowpass, tr=tr
+            )
+        return output
+
+    def fit_transform(
+        self,
+        use_afni: bool,
+        input_data: Dict[str, Any],
+        highpass: float,
+        lowpass: float,
+        tr: Optional[float],
+    ) -> Tuple["Nifti1Image", "Nifti1Image"]:
+        """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
+        -------
+        alff: Niimg-like object
+            ALFF map.
+        falff: Niimg-like object
+            fALFF map.
+        """
+        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 at {self._file_path}.")
+            # Clear the cache
+            self._compute.cache_clear()
+            # Clear temporary directory files
+            for file_ in self.temp_dir_path.iterdir():
+                file_.unlink(missing_ok=True)
+            # Set the new file path
+            self._file_path = bold_path
+        else:
+            logger.info(f"Using fALFF map cache at {self._file_path}.")
+        # Compute
+        return self._compute(
+            use_afni=use_afni,
+            data=bold_data,
+            highpass=highpass,
+            lowpass=lowpass,
+            tr=tr,
+        )

junifer/markers/falff/falff_parcels.py

@@ -0,0 +1,126 @@
+"""Provide class for computing 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>
+# License: AGPL
+
+from typing import Dict, List, Optional, Union
+
+
+from ...api.decorators import register_marker
+from .falff_base import AmplitudeLowFrequencyFluctuationBase
+from .. import ParcelAggregation
+
+
+@register_marker
+class AmplitudeLowFrequencyFluctuationParcels(
+    AmplitudeLowFrequencyFluctuationBase
+):
+    """Class for computing fALFF/ALFF on parcels.
+
+    Parameters
+    ----------
+    parcellation : str or list of str
+        The name(s) of the parcellation(s). Check valid options by calling
+        :func:`junifer.data.parcellations.list_parcellations`.
+    fractional : bool
+        Whether to compute fractional ALFF.
+    highpass : positive float, optional
+        The highpass cutoff frequency for the bandpass filter (default 0.01).
+    lowpass : positive float, optional
+        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).
+    use_afni : bool, optional
+        Whether to use AFNI for computing. If None, will use AFNI only
+        if available (default None).
+    mask : str, optional
+        The name of the mask to apply to regions before extracting signals.
+        Check valid options by calling :func:`junifer.data.masks.list_masks`
+        (default None).
+    method : str, optional
+        The method to perform aggregation using. Check valid options in
+        :func:`junifer.stats.get_aggfunc_by_name` (default "mean").
+    method_params : dict, optional
+        Parameters to pass to the aggregation function. Check valid options in
+        :func:`junifer.stats.get_aggfunc_by_name`.
+    name : str, optional
+        The name of the marker. If None, will use the class name (default
+        None).
+
+    Notes
+    -----
+    The ``tr`` parameter is crucial for the correctness of fALFF/ALFF
+    computation. If a dataset is correctly preprocessed, the TR should be
+    extracted from the NIFTI without any issue. However, it has been
+    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
+    :func:`scipy.signal.butter` and :func:`scipy.signal.filtfilt` for more
+    details.
+    """
+
+    def __init__(
+        self,
+        parcellation: Union[str, List[str]],
+        fractional: bool,
+        highpass: float = 0.01,
+        lowpass: float = 0.1,
+        tr: Optional[float] = None,
+        use_afni: Optional[bool] = None,
+        mask: Optional[str] = None,
+        method: str = "mean",
+        method_params: Optional[Dict] = None,
+        name: Optional[str] = None,
+    ) -> None:
+        self.parcellation = parcellation
+        self.mask = mask
+        self.method = method
+        self.method_params = method_params
+        super().__init__(
+            fractional=fractional,
+            highpass=highpass,
+            lowpass=lowpass,
+            tr=tr,
+            name=name,
+            use_afni=use_afni,
+        )
+
+    def _postprocess(self, input: Dict) -> Dict:
+        """Compute ALFF and fALFF.
+
+        Parameters
+        ----------
+        input : dict
+            A single input from the pipeline data object in which to compute
+            the marker.
+        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).
+
+        Returns
+        -------
+        dict
+            The computed ALFF as dictionary. The dictionary has the following
+            keys:
+
+            * ``data`` : the actual computed values as a numpy.ndarray
+            * ``columns`` : the column labels for the computed values as a list
+        """
+        pa = ParcelAggregation(
+            parcellation=self.parcellation,
+            method=self.method,
+            method_params=self.method_params,
+            mask=self.mask,
+            on="fALFF",
+        )
+
+        # get the 2D timeseries after parcel aggregation
+        out = pa.compute(input)
+
+        return out

junifer/markers/falff/falff_spheres.py

@@ -0,0 +1,134 @@
+"""Provide class for computing 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>
+# License: AGPL
+
+from typing import Dict, Optional
+
+
+from ...api.decorators import register_marker
+from .falff_base import AmplitudeLowFrequencyFluctuationBase
+from .. import SphereAggregation
+
+
+@register_marker
+class AmplitudeLowFrequencyFluctuationSpheres(
+    AmplitudeLowFrequencyFluctuationBase
+):
+    """Class for computing fALFF/ALFF on spheres.
+
+    Parameters
+    ----------
+    coords : str
+        The name of the coordinates list to use. See
+        :func:`junifer.data.coordinates.list_coordinates` for options.
+    radius : float, optional
+        The radius of the sphere in mm. If None, the signal will be extracted
+        from a single voxel. See :class:`nilearn.maskers.NiftiSpheresMasker`
+        for more information (default None).
+    fractional : bool
+        Whether to compute fractional ALFF.
+    highpass : positive float, optional
+        The highpass cutoff frequency for the bandpass filter (default 0.01).
+    lowpass : positive float, optional
+        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).
+    use_afni : bool, optional
+        Whether to use AFNI for computing. If None, will use AFNI only
+        if available (default None).
+    mask : str, optional
+        The name of the mask to apply to regions before extracting signals.
+        Check valid options by calling :func:`junifer.data.masks.list_masks`
+        (default None).
+    method : str, optional
+        The method to perform aggregation using. Check valid options in
+        :func:`junifer.stats.get_aggfunc_by_name` (default "mean").
+    method_params : dict, optional
+        Parameters to pass to the aggregation function. Check valid options in
+        :func:`junifer.stats.get_aggfunc_by_name`.
+    name : str, optional
+        The name of the marker. If None, will use the class name (default
+        None).
+
+    Notes
+    -----
+    The ``tr`` parameter is crucial for the correctness of fALFF/ALFF
+    computation. If a dataset is correctly preprocessed, the TR should be
+    extracted from the NIFTI without any issue. However, it has been
+    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
+    :func:`scipy.signal.butter` and :func:`scipy.signal.filtfilt` for more
+    details.
+    """
+
+    def __init__(
+        self,
+        coords: str,
+        fractional: bool,
+        radius: Optional[float] = None,
+        highpass: float = 0.01,
+        lowpass: float = 0.1,
+        tr: Optional[float] = None,
+        use_afni: Optional[bool] = None,
+        mask: Optional[str] = None,
+        method: str = "mean",
+        method_params: Optional[Dict] = None,
+        name: Optional[str] = None,
+    ) -> None:
+        self.coords = coords
+        self.radius = radius
+        self.mask = mask
+        self.method = method
+        self.method_params = method_params
+        super().__init__(
+            fractional=fractional,
+            highpass=highpass,
+            lowpass=lowpass,
+            tr=tr,
+            name=name,
+            use_afni=use_afni,
+        )
+
+    def _postprocess(self, input: Dict) -> Dict:
+        """Compute ALFF and fALFF.
+
+        Parameters
+        ----------
+        input : dict
+            A single input from the pipeline data object in which to compute
+            the marker.
+        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).
+
+        Returns
+        -------
+        dict
+            The computed ALFF as dictionary. The dictionary has the following
+            keys:
+
+            * ``data`` : the actual computed values as a numpy.ndarray
+            * ``columns`` : the column labels for the computed values as a list
+
+        """
+        pa = SphereAggregation(
+            coords=self.coords,
+            radius=self.radius,
+            method=self.method,
+            method_params=self.method_params,
+            mask=self.mask,
+            on="fALFF",
+        )
+
+        # get the 2D timeseries after parcel aggregation
+        out = pa.compute(input)
+
+        return out

junifer/markers/falff/tests/test_falff_estimator.py

@@ -0,0 +1,241 @@
+"""Provide test for (f)ALFF estimator."""
+
+# Authors: Federico Raimondo <f.raimondo@fz-juelich.de>
+# License: AGPL
+
+import pytest
+import time
+from scipy.stats import pearsonr
+from nibabel import Nifti1Image
+
+from junifer.datareader import DefaultDataReader
+from junifer.markers.falff.falff_estimator import (
+    AmplitudeLowFrequencyFluctuationEstimator,
+)
+from junifer.testing.datagrabbers import (
+    PartlyCloudyTestingDataGrabber,
+)
+from junifer.pipeline.utils import _check_afni
+from junifer.utils import logger
+
+
+def test_AmplitudeLowFrequencyFluctuationEstimator_cache_python() -> None:
+    """Test that the cache works properly when using python."""
+    with PartlyCloudyTestingDataGrabber() as dg:
+        input = dg["sub-01"]
+
+    input = DefaultDataReader().fit_transform(input)
+
+    estimator = AmplitudeLowFrequencyFluctuationEstimator()
+    start_time = time.time()
+    alff, falff = estimator.fit_transform(
+        use_afni=False,
+        input_data=input["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)
+    n_files = len([x for x in estimator.temp_dir_path.glob("*")])
+    assert n_files == 0  # no files in python
+
+    # Now fit again, should be faster
+    start_time = time.time()
+    alff, falff = estimator.fit_transform(
+        use_afni=False,
+        input_data=input["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([x for x in estimator.temp_dir_path.glob("*")])
+    assert n_files == 0  # no files in python
+
+    # Now change a parameter, should compute again, without clearing the
+    # cache
+    start_time = time.time()
+    alff, falff = estimator.fit_transform(
+        use_afni=False,
+        input_data=input["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([x for x in estimator.temp_dir_path.glob("*")])
+    assert n_files == 0  # no files in python
+
+    # Now fit again with the previous params, should be fast
+    start_time = time.time()
+    alff, falff = estimator.fit_transform(
+        use_afni=False,
+        input_data=input["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([x for x in estimator.temp_dir_path.glob("*")])
+    assert n_files == 0  # no files in python
+
+    # Now change the data, it should clear the cache
+    with PartlyCloudyTestingDataGrabber() as dg:
+        input = dg["sub-02"]
+
+    input = DefaultDataReader().fit_transform(input)
+
+    start_time = time.time()
+    alff, falff = estimator.fit_transform(
+        use_afni=False,
+        input_data=input["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([x for x in estimator.temp_dir_path.glob("*")])
+    assert n_files == 0  # no files in python
+
+
+@pytest.mark.skipif(
+    _check_afni() is False, reason="requires afni to be in PATH"
+)
+def test_AmplitudeLowFrequencyFluctuationEstimator_cache_afni() -> None:
+    """Test that the cache works properly when using afni."""
+    with PartlyCloudyTestingDataGrabber() as dg:
+        input = dg["sub-01"]
+
+    input = DefaultDataReader().fit_transform(input)
+
+    estimator = AmplitudeLowFrequencyFluctuationEstimator()
+    start_time = time.time()
+    alff, falff = estimator.fit_transform(
+        use_afni=True,
+        input_data=input["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)
+    n_files = len([x for x in estimator.temp_dir_path.glob("*")])
+    assert n_files == 3  # input + alff + falff
+
+    # Now fit again, should be faster
+    start_time = time.time()
+    alff, falff = estimator.fit_transform(
+        use_afni=True,
+        input_data=input["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([x for x in estimator.temp_dir_path.glob("*")])
+    assert n_files == 3  # input + alff + falff
+
+    # Now change a parameter, should compute again, without clearing the
+    # cache
+    start_time = time.time()
+    alff, falff = estimator.fit_transform(
+        use_afni=True,
+        input_data=input["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([x for x in estimator.temp_dir_path.glob("*")])
+    assert n_files == 5  # input + 2 * alff + 2 * falff
+
+    # Now fit again with the previous params, should be fast
+    start_time = time.time()
+    alff, falff = estimator.fit_transform(
+        use_afni=True,
+        input_data=input["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([x for x in estimator.temp_dir_path.glob("*")])
+    assert n_files == 5  # input + 2 * alff + 2 * falff
+
+    # Now change the data, it should clear the cache
+    with PartlyCloudyTestingDataGrabber() as dg:
+        input = dg["sub-02"]
+
+    input = DefaultDataReader().fit_transform(input)
+
+    start_time = time.time()
+    alff, falff = estimator.fit_transform(
+        use_afni=True,
+        input_data=input["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([x for x in estimator.temp_dir_path.glob("*")])
+    assert n_files == 3  # input + alff + falff
+
+
+@pytest.mark.skipif(
+    _check_afni() is False, reason="requires afni to be in PATH"
+)
+def test_AmplitudeLowFrequencyFluctuationEstimator_afni_vs_python() -> None:
+    """Test that the cache works properly when using afni."""
+    with PartlyCloudyTestingDataGrabber() as dg:
+        input = dg["sub-01"]
+
+    input = DefaultDataReader().fit_transform(input)
+    estimator = AmplitudeLowFrequencyFluctuationEstimator()
+
+    # Use an arbitrary TR to test the AFNI vs Python implementation
+    afni_alff, afni_falff = estimator.fit_transform(
+        use_afni=True,
+        input_data=input["BOLD"],
+        highpass=0.01,
+        lowpass=0.1,
+        tr=2.5,
+    )
+
+    python_alff, python_falff = estimator.fit_transform(
+        use_afni=False,
+        input_data=input["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

@@ -0,0 +1,157 @@
+"""Provide test for parcel-aggregated (f)ALFF."""
+
+# Authors: Federico Raimondo <f.raimondo@fz-juelich.de>
+#          Synchon Mandal <s.mandal@fz-juelich.de>
+# License: AGPL
+
+import pytest
+
+from pathlib import Path
+from numpy.testing import assert_array_equal
+from scipy.stats import pearsonr
+
+from junifer.datareader import DefaultDataReader
+from junifer.markers.falff import AmplitudeLowFrequencyFluctuationParcels
+from junifer.testing.datagrabbers import PartlyCloudyTestingDataGrabber
+from junifer.pipeline.utils import _check_afni
+from junifer.storage import SQLiteFeatureStorage
+from junifer.utils import logger
+
+
+_PARCELLATION = "Schaefer100x7"
+
+
+def test_AmplitudeLowFrequencyFluctuationParcels_python() -> None:
+    """Test AmplitudeLowFrequencyFluctuationParcels using python."""
+    # Get the SPM auditory data:
+
+    with PartlyCloudyTestingDataGrabber() as dg:
+        input = dg["sub-01"]
+
+    input = DefaultDataReader().fit_transform(input)
+    # Create ParcelAggregation object
+    marker = AmplitudeLowFrequencyFluctuationParcels(
+        parcellation=_PARCELLATION,
+        method="mean",
+        use_afni=False,
+        fractional=False,
+    )
+    python_values = marker.fit_transform(input)["BOLD"]["data"]
+
+    assert marker.use_afni is False
+    assert python_values.ndim == 2
+    assert python_values.shape == (1, 100)
+
+
+@pytest.mark.skipif(
+    _check_afni() is False, reason="requires afni to be in PATH"
+)
+def test_AmplitudeLowFrequencyFluctuationParcels_afni() -> None:
+    """Test AmplitudeLowFrequencyFluctuationParcels using afni."""
+    # Get the SPM auditory data:
+    with PartlyCloudyTestingDataGrabber() as dg:
+        input = dg["sub-01"]
+
+    input = DefaultDataReader().fit_transform(input)
+    # Create ParcelAggregation object
+    marker = AmplitudeLowFrequencyFluctuationParcels(
+        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 = AmplitudeLowFrequencyFluctuationParcels(
+        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_AmplitudeLowFrequencyFluctuationParcels_python_vs_afni(
+    fractional: bool,
+) -> None:
+    """Test AmplitudeLowFrequencyFluctuationParcels using python.
+
+    Parameters
+    ----------
+    factional : bool
+        Whether to compute fractional ALFF or not.
+    """
+
+    with PartlyCloudyTestingDataGrabber() as dg:
+        input = dg["sub-01"]
+
+    input = DefaultDataReader().fit_transform(input)
+    # Create ParcelAggregation object
+    marker_python = AmplitudeLowFrequencyFluctuationParcels(
+        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 = AmplitudeLowFrequencyFluctuationParcels(
+        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_AmplitudeLowFrequencyFluctuationParcels_storage(
+    tmp_path: Path,
+) -> None:
+    """Test AmplitudeLowFrequencyFluctuationParcels 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)
+        # Create ParcelAggregation object
+        marker = AmplitudeLowFrequencyFluctuationParcels(
+            parcellation=_PARCELLATION,
+            method="mean",
+            use_afni=False,
+            fractional=True,
+        )
+        storage = SQLiteFeatureStorage(tmp_path / "alff_parcels.sqlite")
+
+        # Fit transform marker on data with storage
+        marker.fit_transform(
+            input=input,
+            storage=storage,
+        )

junifer/markers/falff/tests/test_falff_spheres.py

@@ -0,0 +1,165 @@
+"""Provide test for sphere-aggregated (f)ALFF."""
+
+# Authors: Federico Raimondo <f.raimondo@fz-juelich.de>
+#          Synchon Mandal <s.mandal@fz-juelich.de>
+# License: AGPL
+
+import pytest
+
+from pathlib import Path
+
+from numpy.testing import assert_array_equal
+from scipy.stats import pearsonr
+
+from junifer.datareader import DefaultDataReader
+from junifer.markers.falff import AmplitudeLowFrequencyFluctuationSpheres
+from junifer.testing.datagrabbers import PartlyCloudyTestingDataGrabber
+from junifer.pipeline.utils import _check_afni
+from junifer.storage import SQLiteFeatureStorage
+from junifer.utils import logger
+
+
+_COORDINATES = "DMNBuckner"
+
+
+def test_AmplitudeLowFrequencyFluctuationSpheres_python() -> None:
+    """Test AmplitudeLowFrequencyFluctuationSpheres using python."""
+    # Get the SPM auditory data:
+
+    with PartlyCloudyTestingDataGrabber() as dg:
+        input = dg["sub-01"]
+
+    input = DefaultDataReader().fit_transform(input)
+    # Create ParcelAggregation object
+    marker = AmplitudeLowFrequencyFluctuationSpheres(
+        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
+    assert python_values.ndim == 2
+    assert python_values.shape == (1, 6)
+
+
+@pytest.mark.skipif(
+    _check_afni() is False, reason="requires afni to be in PATH"
+)
+def test_AmplitudeLowFrequencyFluctuationSpheres_afni() -> None:
+    """Test AmplitudeLowFrequencyFluctuationSpheres using afni."""
+    # Get the SPM auditory data:
+    with PartlyCloudyTestingDataGrabber() as dg:
+        input = dg["sub-01"]
+
+    input = DefaultDataReader().fit_transform(input)
+    # Create ParcelAggregation object
+    marker = AmplitudeLowFrequencyFluctuationSpheres(
+        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 = AmplitudeLowFrequencyFluctuationSpheres(
+        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_AmplitudeLowFrequencyFluctuationSpheres_python_vs_afni(
+    fractional: bool,
+) -> None:
+    """Test AmplitudeLowFrequencyFluctuationSpheres python vs afni results.
+
+    Parameters
+    ----------
+    fractional : bool
+        Whether to compute fractional ALFF or not.
+    """
+    with PartlyCloudyTestingDataGrabber() as dg:
+        input = dg["sub-01"]
+
+    input = DefaultDataReader().fit_transform(input)
+    # Create ParcelAggregation object
+    marker_python = AmplitudeLowFrequencyFluctuationSpheres(
+        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 = AmplitudeLowFrequencyFluctuationSpheres(
+        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_AmplitudeLowFrequencyFluctuationSpheres_storage(
+    tmp_path: Path,
+) -> None:
+    """Test AmplitudeLowFrequencyFluctuationSpheres 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)
+        # Create ParcelAggregation object
+        marker = AmplitudeLowFrequencyFluctuationSpheres(
+            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
+        marker.fit_transform(
+            input=input,
+            storage=storage,
+        )