diff --git a/docs/changes/newsfragments/311.change b/docs/changes/newsfragments/311.change new file mode 100644 index 000000000..cf6117d92 --- /dev/null +++ b/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`_ diff --git a/docs/changes/newsfragments/311.enh b/docs/changes/newsfragments/311.enh new file mode 100644 index 000000000..9cfb9e82f --- /dev/null +++ b/docs/changes/newsfragments/311.enh @@ -0,0 +1 @@ +Improve external dependency handling for :class:`.PipelineStepMixin`-derived objects having implementation-based variations by `Synchon Mandal`_ diff --git a/junifer/markers/falff/_afni_falff.py b/junifer/markers/falff/_afni_falff.py new file mode 100644 index 000000000..81d1ec669 --- /dev/null +++ b/junifer/markers/falff/_afni_falff.py @@ -0,0 +1,153 @@ +"""Provide class for computing ALFF using AFNI.""" + +# Authors: Synchon Mandal +# 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 diff --git a/junifer/markers/falff/_junifer_falff.py b/junifer/markers/falff/_junifer_falff.py new file mode 100644 index 000000000..3dc6fd083 --- /dev/null +++ b/junifer/markers/falff/_junifer_falff.py @@ -0,0 +1,142 @@ +"""Provide class for computing ALFF using junifer.""" + +# Authors: Synchon Mandal +# 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 diff --git a/junifer/markers/falff/falff_base.py b/junifer/markers/falff/falff_base.py index f57411474..7bdd86c0a 100644 --- a/junifer/markers/falff/falff_base.py +++ b/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 # Amir Omidvarnia # Kaustubh R. Patil +# Synchon Mandal # License: AGPL -from abc import abstractmethod -from typing import ClassVar, Dict, List, Optional, Union +from pathlib import Path +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). - use_afni : bool, optional - Whether to use AFNI for computing. If None, will use AFNI only - if available (default None). + the TR from NIfTI header (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. + 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. + + 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[ - List[Dict[str, Union[str, bool, List[str]]]] - ] = [ + _CONDITIONAL_DEPENDENCIES: ClassVar[List[Dict[str, Union[str, Type]]]] = [ { - "name": "afni", - "optional": True, - "commands": ["3dRSFC", "3dAFNItoNIFTI"], + "using": "afni", + "depends_on": AFNIALFF, + }, + { + "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], - extra_input: Optional[Dict] = None, - ) -> Dict: - """Compute. + input_data: Dict[str, Any], + ) -> Tuple["Nifti1Image", Path]: + """Compute ALFF and fALFF. Parameters ---------- - input : dict - A single input from the pipeline data object in which to compute - the marker. + input_data : dict + The input to 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). + The other fields in the pipeline data object (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 - * ``col_names`` : the column labels for the computed values as list + Niimg-like object + The ALFF / fALFF as NIfTI. + pathlib.Path + The path to the ALFF / fALFF as NIfTI. """ - if self.use_afni is None: - raise_error( - "Parameter `use_afni` must be set to True or False in order " - "to compute this marker. It is currently set to None (default " - "behaviour). This is intended to be for auto-detection. In " - "order for that to happen, please call the `validate` method " - "before calling the `compute` method." - ) + logger.debug("Calculating ALFF and fALFF") - estimator = ALFFEstimator() - - # If the input data space is "native", then alff_path and falff_path - # both point to the input data path as it might be required to use - # in get_corrdinates() for transforming coordinates to native space. - alff, falff, alff_path, falff_path = estimator.fit_transform( - use_afni=self.use_afni, - input_data=input, + # Conditional estimator + if self.using == "afni": + estimator = AFNIALFF() + elif self.using == "junifer": + estimator = JuniferALFF() + # Compute ALFF + fALFF + alff, falff, alff_path, falff_path = estimator.compute( # type: ignore + 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()) - post_input["data"] = post_data - post_input["path"] = post_path - - out = self._postprocess(post_input, extra_input=extra_input) - - return out - - @abstractmethod - def _postprocess( - self, input: Dict, extra_input: Optional[Dict] = None - ) -> Dict: - """Postprocess the output of the estimator. - - 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 - ) + # 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 + # path should be passed for use later if required. + # TODO(synchon): will be taken care in #292 + if input_data["space"] == "native" and self.fractional: + return falff, input_data["path"] + elif input_data["space"] == "native" and not self.fractional: + return alff, input_data["path"] + elif input_data["space"] != "native" and self.fractional: + return falff, falff_path + else: + return alff, alff_path diff --git a/junifer/markers/falff/falff_estimator.py b/junifer/markers/falff/falff_estimator.py deleted file mode 100644 index f212710be..000000000 --- a/junifer/markers/falff/falff_estimator.py +++ /dev/null @@ -1,349 +0,0 @@ -"""Provide estimator class for (f)ALFF.""" - -# Authors: Synchon Mandal -# Federico Raimondo -# 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 diff --git a/junifer/markers/falff/falff_parcels.py b/junifer/markers/falff/falff_parcels.py index 25e8fcfc8..8741a219d 100644 --- a/junifer/markers/falff/falff_parcels.py +++ b/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 # Amir Omidvarnia # Kaustubh R. Patil +# Synchon Mandal # 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). - use_afni : bool, optional - Whether to use AFNI for computing. If None, will use AFNI only - if available (default None). + the TR from NIfTI header (default None). + agg_method : str, optional + The method to perform aggregation using. Check valid options in + :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 ` 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 - 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. + 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 @@ -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 - self.masks = masks - self.method = method - self.method_params = method_params + # Superclass init first to validate `using` parameter 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( - self, input: Dict, extra_input: Optional[Dict] = None - ) -> Dict: - """Compute ALFF and fALFF. + def compute( + self, + input: Dict[str, Any], + 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 marker. + The BOLD data as dictionary. 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). + The other fields in the pipeline data object (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_params=self.method_params, + method=self.agg_method, + 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 - out = pa.compute(input, extra_input=extra_input) - - return out + return output diff --git a/junifer/markers/falff/falff_spheres.py b/junifer/markers/falff/falff_spheres.py index 9768eadff..d60d3c2a2 100644 --- a/junifer/markers/falff/falff_spheres.py +++ b/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 # Amir Omidvarnia # Kaustubh R. Patil +# Synchon Mandal # 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). - use_afni : bool, optional - Whether to use AFNI for computing. If None, will use AFNI only - if available (default None). + the TR from NIfTI header (default None). + agg_method : str, optional + The method to perform aggregation using. Check valid options in + :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 ` 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 - extracted from the NIFTI without any issue. However, it has been - reported that some preprocessed data might not have the correct TR in + 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 @@ -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 - self.radius = radius - self.allow_overlap = allow_overlap - self.masks = masks - self.method = method - self.method_params = method_params + # Superclass init first to validate `using` parameter 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( - self, input: Dict, extra_input: Optional[Dict] = None - ) -> Dict: - """Compute ALFF and fALFF. + def compute( + self, + input: Dict[str, Any], + 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 marker. + The BOLD data as dictionary. 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). + The other fields in the pipeline data object (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_params=self.method_params, + method=self.agg_method, + 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 - out = pa.compute(input, extra_input=extra_input) - - return out + return output diff --git a/junifer/markers/falff/tests/test_falff_estimator.py b/junifer/markers/falff/tests/test_falff_estimator.py deleted file mode 100644 index 4025de807..000000000 --- a/junifer/markers/falff/tests/test_falff_estimator.py +++ /dev/null @@ -1,267 +0,0 @@ -"""Provide test for (f)ALFF estimator.""" - -# Authors: Federico Raimondo -# 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 diff --git a/junifer/markers/falff/tests/test_falff_parcels.py b/junifer/markers/falff/tests/test_falff_parcels.py index 38c907d7c..4356c4954 100644 --- a/junifer/markers/falff/tests/test_falff_parcels.py +++ b/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 # Synchon Mandal # License: AGPL +import logging from pathlib import Path import pytest -from numpy.testing import assert_array_equal -from scipy.stats import pearsonr +import scipy as sp 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: - """Test ALFFParcels using python. +def test_ALFFParcels(caplog: pytest.LogCaptureFixture, tmp_path: Path) -> None: + """Test ALFFParcels. Parameters ---------- + caplog : pytest.LogCaptureFixture + The pytest.LogCaptureFixture object. tmp_path : pathlib.Path The path to the test directory. """ - with PartlyCloudyTestingDataGrabber() as dg: - input_ = dg["sub-01"] + with caplog.at_level(logging.DEBUG): + 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_) - WorkDirManager().workdir = tmp_path - marker = ALFFParcels( - parcellation=_PARCELLATION, - method="mean", - use_afni=False, - fractional=False, - ) - python_values = marker.fit_transform(input_)["BOLD"]["data"] + # Initialize marker + marker = ALFFParcels( + parcellation=PARCELLATION, + fractional=False, + using="junifer", + ) + # Fit transform marker on data + output = marker.fit_transform(element_data) - assert marker.use_afni is False - assert python_values.ndim == 2 - assert python_values.shape == (1, 100) + assert "Creating cache" in caplog.text + + # 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" -) -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" + _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( - tmp_path: Path, - fractional: bool, -) -> None: - """Test ALFFParcels using python. +def test_ALFFParcels_comparison(tmp_path: Path, fractional: bool) -> None: + """Test ALFFParcels implementation comparison. Parameters ---------- @@ -109,65 +92,35 @@ def test_ALFFParcels_python_vs_afni( """ with PartlyCloudyTestingDataGrabber() as dg: - input_ = dg["sub-01"] - - 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_) + element_data = DefaultDataReader().fit_transform(dg["sub-01"]) + # Update workdir to current test's tmp_path 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 - marker.fit_transform( - input=input_, - storage=storage, + # Initialize marker + junifer_marker = ALFFParcels( + parcellation=PARCELLATION, + 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 diff --git a/junifer/markers/falff/tests/test_falff_spheres.py b/junifer/markers/falff/tests/test_falff_spheres.py index e4e7c1e19..115c6c205 100644 --- a/junifer/markers/falff/tests/test_falff_spheres.py +++ b/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 # Synchon Mandal # License: AGPL +import logging from pathlib import Path import pytest -from numpy.testing import assert_array_equal -from scipy.stats import pearsonr +import scipy as sp 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: - """Test ALFFSpheres using python. +def test_ALFFSpheres(caplog: pytest.LogCaptureFixture, tmp_path: Path) -> None: + """Test ALFFSpheres. Parameters ---------- + caplog : pytest.LogCaptureFixture + The pytest.LogCaptureFixture object. tmp_path : pathlib.Path The path to the test directory. """ - with PartlyCloudyTestingDataGrabber() as dg: - input_ = dg["sub-01"] + with caplog.at_level(logging.DEBUG): + 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 - marker = ALFFSpheres( - coords=_COORDINATES, - radius=5, - method="mean", - use_afni=False, - fractional=False, - ) - python_values = marker.fit_transform(input_)["BOLD"]["data"] + assert "Creating cache" in caplog.text - assert marker.use_afni is False - assert python_values.ndim == 2 - assert python_values.shape == (1, 6) + # 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, 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. - - """ - 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. + 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"] - - 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_) + element_data = DefaultDataReader().fit_transform(dg["sub-01"]) + # Update workdir to current test's tmp_path 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 - marker.fit_transform( - input=input_, - storage=storage, + # Initialize marker + junifer_marker = ALFFSpheres( + coords=COORDINATES, + 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 diff --git a/junifer/markers/reho/_afni_reho.py b/junifer/markers/reho/_afni_reho.py new file mode 100644 index 000000000..3cbaf8d12 --- /dev/null +++ b/junifer/markers/reho/_afni_reho.py @@ -0,0 +1,192 @@ +"""Provide class for computing regional homogeneity (ReHo) using AFNI.""" + +# Authors: Synchon Mandal +# 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 diff --git a/junifer/markers/reho/_junifer_reho.py b/junifer/markers/reho/_junifer_reho.py new file mode 100644 index 000000000..2f1d0d12f --- /dev/null +++ b/junifer/markers/reho/_junifer_reho.py @@ -0,0 +1,281 @@ +"""Provide class for computing regional homogeneity (ReHo) using junifer.""" + +# Authors: Synchon Mandal +# 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 diff --git a/junifer/markers/reho/reho_base.py b/junifer/markers/reho/reho_base.py index fc9b51ba4..1603754d3 100644 --- a/junifer/markers/reho/reho_base.py +++ b/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 - Whether to use AFNI for computing. If None, will use AFNI only - if available (default None). + using : {"junifer", "afni"} + Implementation to use for computing ReHo: + + * "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[ - List[Dict[str, Union[str, bool, List[str]]]] - ] = [ + _CONDITIONAL_DEPENDENCIES: ClassVar[List[Dict[str, Union[str, Type]]]] = [ { - "name": "afni", - "optional": True, - "commands": ["3dReHo", "3dAFNItoNIFTI"], + "using": "afni", + "depends_on": AFNIReHo, + }, + { + "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: - raise_error( - "Parameter `use_afni` must be set to True or False in order " - "to compute this marker. It is currently set to None (default " - "behaviour). This is intended to be for auto-detection. In " - "order for that to happen, please call the `validate` method " - "before calling the `compute` method." - ) - logger.info("Calculating ReHO map.") - # Initialize reho estimator - reho_estimator = ReHoEstimator() - # Fit-transform reho estimator - return reho_estimator.fit_transform( - use_afni=self.use_afni, - input_data=input, + logger.debug("Calculating voxel-wise ReHo") + + # Conditional estimator + if self.using == "afni": + estimator = AFNIReHo() + elif self.using == "junifer": + estimator = JuniferReHo() + # Compute reho + reho_map, reho_map_path = estimator.compute( # type: ignore + data=input_data["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 + # 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 diff --git a/junifer/markers/reho/reho_estimator.py b/junifer/markers/reho/reho_estimator.py deleted file mode 100644 index 9ebfdc2e9..000000000 --- a/junifer/markers/reho/reho_estimator.py +++ /dev/null @@ -1,510 +0,0 @@ -"""Provide estimator class for regional homogeneity (ReHo).""" - -# Authors: Synchon Mandal -# Federico Raimondo -# 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 diff --git a/junifer/markers/reho/reho_parcels.py b/junifer/markers/reho/reho_parcels.py index b3f1133d1..c00e8841d 100644 --- a/junifer/markers/reho/reho_parcels.py +++ b/junifer/markers/reho/reho_parcels.py @@ -20,15 +20,18 @@ class ReHoParcels(ReHoBase): Parameters ---------- - parcellation : str - The name of the parcellation. Check valid options by calling + parcellation : str or list of str + The name(s) of the parcellation(s). Check valid options by calling :func:`.list_parcellations`. - use_afni : bool, optional - Whether to use AFNI for computing. If None, will use AFNI only - if available (default None). + using : {"junifer", "afni"} + Implementation to use for computing ReHo: + + * "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, - use_afni: Optional[bool] = None, + parcellation: Union[str, List[str]], + 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.") - # Calculate reho map + logger.info("Calculating ReHo for parcels") + + # 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 - # get_coordinates() for transforming coordinates to native space. + # the input data path as it might be required for parcellation + # transformation to native space. if self.reho_params is not None: - reho_map, reho_file_path = self.compute_reho_map( - input=input, **self.reho_params + reho_map, reho_file_path = self._compute( + 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, diff --git a/junifer/markers/reho/reho_spheres.py b/junifer/markers/reho/reho_spheres.py index 73faaaaed..453004b71 100644 --- a/junifer/markers/reho/reho_spheres.py +++ b/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.") - # Calculate reho map + logger.info("Calculating ReHo for spheres") + + # 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 - # get_coordinates() for transforming coordinates to native space. + # the input data path as it might be required for coordinates + # transformation to native space. if self.reho_params is not None: - reho_map, reho_file_path = self.compute_reho_map( - input=input, **self.reho_params + reho_map, reho_file_path = self._compute( + 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, diff --git a/junifer/markers/reho/tests/test_reho_estimator.py b/junifer/markers/reho/tests/test_reho_estimator.py deleted file mode 100644 index f171ad62f..000000000 --- a/junifer/markers/reho/tests/test_reho_estimator.py +++ /dev/null @@ -1,314 +0,0 @@ -"""Provide tests for ReHo map compute comparison.""" - -# Authors: Synchon Mandal -# 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 diff --git a/junifer/markers/reho/tests/test_reho_parcels.py b/junifer/markers/reho/tests/test_reho_parcels.py index 6bd7564bd..531fe4fc8 100644 --- a/junifer/markers/reho/tests/test_reho_parcels.py +++ b/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 # License: AGPL +import logging from pathlib import Path import pytest -from nilearn import image as nimg -from scipy.stats import pearsonr +import scipy as sp -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: - """Test ReHoParcels fit-transform. +def test_ReHoParcels(caplog: pytest.LogCaptureFixture, tmp_path: Path) -> None: + """Test ReHoParcels. Parameters ---------- + caplog : pytest.LogCaptureFixture + The pytest.LogCaptureFixture object. 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) - # 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 + with caplog.at_level(logging.DEBUG): + with SPMAuditoryTestingDataGrabber() as dg: + element_data = DefaultDataReader().fit_transform(dg["sub001"]) + # Update workdir to current test's tmp_path + WorkDirManager().workdir = tmp_path - reho_parcels_output_bold_data = reho_parcels_output_bold["data"] - # Assert BOLD output data dimension - assert reho_parcels_output_bold_data.ndim == 2 - # Assert BOLD output data is normalized - assert (reho_parcels_output_bold_data > 0).all() and ( - reho_parcels_output_bold_data < 1 - ).all() + # Initialize marker + marker = ReHoParcels(parcellation=PARCELLATION, using="junifer") + # Fit transform marker on data + output = marker.fit_transform(element_data) + + assert "Creating cache" in caplog.text + + # 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: - """Test ReHoParcels fit-transform implementation comparison. +def test_ReHoParcels_comparison(tmp_path: Path) -> None: + """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 - subject_data = dg["sub001"] - # Load image to memory - fmri_img = nimg.load_img(subject_data["BOLD"]["path"]) - + element_data = DefaultDataReader().fit_transform(dg["sub001"]) # 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 - reho_parcels_marker_afni = ReHoParcels( - parcellation=PARCELLATION, use_afni=True + # Initialize marker + junifer_marker = ReHoParcels( + parcellation=PARCELLATION, using="junifer" ) # Fit transform marker on data - reho_parcels_output_afni = reho_parcels_marker_afni.fit_transform( - { - "BOLD": { - "path": "/tmp", - "data": fmri_img, - "meta": {}, - "space": "MNI", - } - } - ) + junifer_output = junifer_marker.fit_transform(element_data) # 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( - reho_parcels_output_bold_python["data"].flatten(), - reho_parcels_output_bold_afni["data"].flatten(), + r, _ = sp.stats.pearsonr( + junifer_output_bold["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, - ) diff --git a/junifer/markers/reho/tests/test_reho_spheres.py b/junifer/markers/reho/tests/test_reho_spheres.py index e3eba63e8..c094cf9eb 100644 --- a/junifer/markers/reho/tests/test_reho_spheres.py +++ b/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 # License: AGPL +import logging from pathlib import Path import pytest -from nilearn import image as nimg -from scipy.stats import pearsonr +import scipy as sp -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: - """Test ReHoSpheres fit-transform. +def test_ReHoSpheres(caplog: pytest.LogCaptureFixture, tmp_path: Path) -> None: + """Test ReHoSpheres. Parameters ---------- + caplog : pytest.LogCaptureFixture + The pytest.LogCaptureFixture object. 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) - # Fit transform marker on data - reho_spheres_output = reho_spheres_marker.fit_transform( - { - "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 + with caplog.at_level(logging.DEBUG): + with SPMAuditoryTestingDataGrabber() as dg: + element_data = DefaultDataReader().fit_transform(dg["sub001"]) + # Update workdir to current test's tmp_path + WorkDirManager().workdir = tmp_path + # Initialize marker + marker = ReHoSpheres( + coords=COORDINATES, using="junifer", radius=10.0 + ) + # Fit transform marker on data + output = marker.fit_transform(element_data) - reho_spheres_output_bold_data = reho_spheres_output_bold["data"] - # Assert BOLD output data dimension - assert reho_spheres_output_bold_data.ndim == 2 - # Assert BOLD output data is normalized - assert (reho_spheres_output_bold_data > 0).all() and ( - reho_spheres_output_bold_data < 1 - ).all() + assert "Creating cache" in caplog.text + + # 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_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: - """Test ReHoSpheres fit-transform implementation comparison. +def test_ReHoSpheres_comparison(tmp_path: Path) -> None: + """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 - subject_data = dg["sub001"] - # Load image to memory - fmri_img = nimg.load_img(subject_data["BOLD"]["path"]) - + element_data = DefaultDataReader().fit_transform(dg["sub001"]) # 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 - reho_spheres_marker_afni = ReHoSpheres( - coords=COORDINATES, radius=10.0, use_afni=True + # Initialize marker + junifer_marker = ReHoSpheres( + coords=COORDINATES, + using="junifer", + radius=10.0, ) # Fit transform marker on data - reho_spheres_output_afni = reho_spheres_marker_afni.fit_transform( - { - "BOLD": { - "path": "/tmp", - "data": fmri_img, - "meta": {}, - "space": "MNI", - } - } - ) + junifer_output = junifer_marker.fit_transform(element_data) # 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( - reho_spheres_output_bold_python["data"].flatten(), - reho_spheres_output_bold_afni["data"].flatten(), + r, _ = sp.stats.pearsonr( + junifer_output_bold["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, - ) diff --git a/junifer/pipeline/pipeline_step_mixin.py b/junifer/pipeline/pipeline_step_mixin.py index 81031bd79..39e701039 100644 --- a/junifer/pipeline/pipeline_step_mixin.py +++ b/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 diff --git a/junifer/pipeline/tests/test_pipeline_step_mixin.py b/junifer/pipeline/tests/test_pipeline_step_mixin.py index 6afccf113..659675ebc 100644 --- a/junifer/pipeline/tests/test_pipeline_step_mixin.py +++ b/junifer/pipeline/tests/test_pipeline_step_mixin.py @@ -1,11 +1,11 @@ -"""Provide tests for pipeline mixin.""" +"""Provide tests for PipelineStepMixin.""" # Authors: Federico Raimondo # Synchon Mandal # 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: - """Test PipelineStepMixin.""" - 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.""" +def test_PipelineStepMixin_correct_dependencies() -> None: + """Test fit-transform 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: - """Test validate with incorrect dependencies.""" +def test_PipelineStepMixin_incorrect_dependencies() -> None: + """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: - """Test validate with correct external dependencies.""" +def test_PipelineStepMixin_correct_ext_dependencies() -> None: + """Test fit-transform with correct external dependencies.""" class CorrectMixer(PipelineStepMixin): """Test class for validation.""" - _EXT_DEPENDENCIES: ClassVar[List[Dict[str, Union[str, bool]]]] = [ - {"name": "afni", "optional": False} + _EXT_DEPENDENCIES: ClassVar[List[Dict[str, str]]] = [{"name": "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({}) + + +@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]: - 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]: + 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: - """Test validate with correct external dependencies' incorrect commands.""" +def test_PipelineStepMixin_ext_deps_incorrect_commands() -> None: + """Test fit-transform with inccorrect external dependency commands.""" class CorrectMixer(PipelineStepMixin): """Test class for validation.""" - _EXT_DEPENDENCIES: ClassVar[ - List[Dict[str, Union[str, bool, List[str]]]] - ] = [{"name": "afni", "optional": False, "commands": ["3d"]}] + _EXT_DEPENDENCIES: ClassVar[List[Dict[str, Union[str, List[str]]]]] = [ + {"name": "afni", "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: - """Test validate with incorrect external dependencies.""" +def test_PipelineStepMixin_incorrect_ext_dependencies() -> None: + """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"): - mixer.validate([]) + with pytest.raises(ValueError, match="Invalid value"): + 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({}) diff --git a/junifer/pipeline/utils.py b/junifer/pipeline/utils.py index 156a100a2..0dfd8ed69 100644 --- a/junifer/pipeline/utils.py +++ b/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 - Whether the dependency is optional. For external dependencies marked - as optional, there should be an implementation provided with junfier. + optional : bool, optional + Whether the dependency is optional (default False). **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 " - "required by one of the pipeline steps." + msg=( + 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 basic fsl flirt version would be looked up, else, would also + The commands to specifically check for from FSL. If None, only + 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( diff --git a/junifer/preprocess/ants/ants_apply_transforms_warper.py b/junifer/preprocess/ants/ants_apply_transforms_warper.py index 66aa24420..9bfa0f379 100644 --- a/junifer/preprocess/ants/ants_apply_transforms_warper.py +++ b/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[ - List[Dict[str, Union[str, bool, List[str]]]] - ] = [ + _EXT_DEPENDENCIES: ClassVar[List[Dict[str, Union[str, 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 - 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 + self.on = on def _run_apply_transforms( self, diff --git a/junifer/preprocess/ants/tests/test_ants_apply_transforms_warper.py b/junifer/preprocess/ants/tests/test_ants_apply_transforms_warper.py index 68b43c648..9b66a327b 100644 --- a/junifer/preprocess/ants/tests/test_ants_apply_transforms_warper.py +++ b/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( diff --git a/junifer/preprocess/bold_warper.py b/junifer/preprocess/bold_warper.py index eb030374f..22aeaceb1 100644 --- a/junifer/preprocess/bold_warper.py +++ b/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[ - List[Dict[str, Union[str, bool, List[str]]]] - ] = [ + _CONDITIONAL_DEPENDENCIES: ClassVar[List[Dict[str, Union[str, Type]]]] = [ { - "name": "fsl", - "optional": True, - "commands": ["flirt", "applywarp"], + "using": "fsl", + "depends_on": _ApplyWarper, }, { - "name": "ants", - "optional": True, - "commands": ["ResampleImage", "antsApplyTransforms"], + "using": "ants", + "depends_on": _AntsApplyTransformsWarper, }, ] - 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": diff --git a/junifer/preprocess/fsl/apply_warper.py b/junifer/preprocess/fsl/apply_warper.py index 2ef437f94..4e7cb1ba5 100644 --- a/junifer/preprocess/fsl/apply_warper.py +++ b/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[ - List[Dict[str, Union[str, bool, List[str]]]] - ] = [ + _EXT_DEPENDENCIES: ClassVar[List[Dict[str, Union[str, 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 - 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 + self.on = on def _run_applywarp( self, diff --git a/junifer/preprocess/fsl/tests/test_apply_warper.py b/junifer/preprocess/fsl/tests/test_apply_warper.py index af2177628..e7a678a75 100644 --- a/junifer/preprocess/fsl/tests/test_apply_warper.py +++ b/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") diff --git a/junifer/preprocess/tests/test_bold_warper.py b/junifer/preprocess/tests/test_bold_warper.py index c36c45fe7..dd9cb1103 100644 --- a/junifer/preprocess/tests/test_bold_warper.py +++ b/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"