diff --git a/docs/builtin.rst b/docs/builtin.rst index df91e9924..a814a751d 100644 --- a/docs/builtin.rst +++ b/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` diff --git a/docs/changes/latest.inc b/docs/changes/latest.inc index 196986ab5..f8c3d88ed 100644 --- a/docs/changes/latest.inc +++ b/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 ~~~~ diff --git a/docs/conf.py b/docs/conf.py index f738b1e68..a11f4ed53 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -155,6 +155,8 @@ numpydoc_xref_ignore = { "class", "objects", "Engine", + "positive", + "negative", } # numpydoc_validation_checks = { # "all", diff --git a/junifer/api/res/afni/3dRSFC b/junifer/api/res/afni/3dRSFC new file mode 100755 index 000000000..f3036ac16 --- /dev/null +++ b/junifer/api/res/afni/3dRSFC @@ -0,0 +1,3 @@ +#!/bin/bash + +run_afni_docker.sh 3dRSFC "$@" \ No newline at end of file diff --git a/junifer/markers/__init__.py b/junifer/markers/__init__.py index 2c853aad0..bc069e3db 100644 --- a/junifer/markers/__init__.py +++ b/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, +) diff --git a/junifer/markers/falff/__init__.py b/junifer/markers/falff/__init__.py new file mode 100644 index 000000000..6e1023665 --- /dev/null +++ b/junifer/markers/falff/__init__.py @@ -0,0 +1,7 @@ +"""Provide imports for falff sub-package.""" + +# Authors: Federico Raimondo +# License: AGPL + +from .falff_parcels import AmplitudeLowFrequencyFluctuationParcels +from .falff_spheres import AmplitudeLowFrequencyFluctuationSpheres diff --git a/junifer/markers/falff/falff_base.py b/junifer/markers/falff/falff_base.py new file mode 100644 index 000000000..dcf028bc8 --- /dev/null +++ b/junifer/markers/falff/falff_base.py @@ -0,0 +1,180 @@ +"""Provide abstract class for computing fALFF.""" + +# Authors: Federico Raimondo +# Amir Omidvarnia +# Kaustubh R. Patil +# 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 + ) diff --git a/junifer/markers/falff/falff_estimator.py b/junifer/markers/falff/falff_estimator.py new file mode 100644 index 000000000..49d0f39d5 --- /dev/null +++ b/junifer/markers/falff/falff_estimator.py @@ -0,0 +1,338 @@ +"""Provide estimator class for (f)ALFF.""" + +# Authors: Synchon Mandal +# Federico Raimondo +# 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, + ) diff --git a/junifer/markers/falff/falff_parcels.py b/junifer/markers/falff/falff_parcels.py new file mode 100644 index 000000000..0d68ca34b --- /dev/null +++ b/junifer/markers/falff/falff_parcels.py @@ -0,0 +1,126 @@ +"""Provide class for computing fALFF on parcels.""" + +# Authors: Federico Raimondo +# Amir Omidvarnia +# Kaustubh R. Patil +# 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 diff --git a/junifer/markers/falff/falff_spheres.py b/junifer/markers/falff/falff_spheres.py new file mode 100644 index 000000000..38c59f26e --- /dev/null +++ b/junifer/markers/falff/falff_spheres.py @@ -0,0 +1,134 @@ +"""Provide class for computing fALFF on spheres.""" + +# Authors: Federico Raimondo +# Amir Omidvarnia +# Kaustubh R. Patil +# 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 diff --git a/junifer/markers/falff/tests/test_falff_estimator.py b/junifer/markers/falff/tests/test_falff_estimator.py new file mode 100644 index 000000000..64c9d4360 --- /dev/null +++ b/junifer/markers/falff/tests/test_falff_estimator.py @@ -0,0 +1,241 @@ +"""Provide test for (f)ALFF estimator.""" + +# Authors: Federico Raimondo +# 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 diff --git a/junifer/markers/falff/tests/test_falff_parcels.py b/junifer/markers/falff/tests/test_falff_parcels.py new file mode 100644 index 000000000..aa12c0cdf --- /dev/null +++ b/junifer/markers/falff/tests/test_falff_parcels.py @@ -0,0 +1,157 @@ +"""Provide test for parcel-aggregated (f)ALFF.""" + +# Authors: Federico Raimondo +# Synchon Mandal +# 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, + ) diff --git a/junifer/markers/falff/tests/test_falff_spheres.py b/junifer/markers/falff/tests/test_falff_spheres.py new file mode 100644 index 000000000..6d10ae463 --- /dev/null +++ b/junifer/markers/falff/tests/test_falff_spheres.py @@ -0,0 +1,165 @@ +"""Provide test for sphere-aggregated (f)ALFF.""" + +# Authors: Federico Raimondo +# Synchon Mandal +# 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, + )