juaml/junifer
Watch
3
0
Fork
You've already forked junifer
 0
Code Issues 14 Pull requests 2 Packages Activity Actions

Improve external dependency handling of PipelineStepMixin #311

Merged
synchon merged 22 commits from refactor/ext-dependency into main 2024-03-14 14:26:24 +00:00
Conversation 5 Commits 22 Files changed 29 +1696 -2524
29 changed files with 1696 additions and 2524 deletions
Download patch file Download diff file Expand all files Collapse all files

1
docs/changes/newsfragments/311.change Normal file
View file

@ -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`_

1
docs/changes/newsfragments/311.enh Normal file
View file

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

153
junifer/markers/falff/_afni_falff.py Normal file
View file

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

142
junifer/markers/falff/_junifer_falff.py Normal file
View file

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

176
junifer/markers/falff/falff_base.py
View file

@ -1,16 +1,32 @@
"""Provide abstract class for computing fALFF."""
"""Provide base class for ALFF / fALFF."""
# Authors: Federico Raimondo <f.raimondo@fz-juelich.de>
# Amir Omidvarnia <a.omidvarnia@fz-juelich.de>
# Kaustubh R. Patil <k.patil@fz-juelich.de>
# Synchon Mandal <s.mandal@fz-juelich.de>
# License: AGPL
from abc import abstractmethod
from 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

349
junifer/markers/falff/falff_estimator.py
View file

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

105
junifer/markers/falff/falff_parcels.py
View file

@ -1,20 +1,22 @@
"""Provide class for computing fALFF on parcels."""
"""Provide class for ALFF / fALFF on parcels."""
# Authors: Federico Raimondo <f.raimondo@fz-juelich.de>
# Amir Omidvarnia <a.omidvarnia@fz-juelich.de>
# Kaustubh R. Patil <k.patil@fz-juelich.de>
# Synchon Mandal <s.mandal@fz-juelich.de>
# License: AGPL
from typing import Dict, List, Optional, Union
from typing import Any, Dict, List, Optional, Union
from ...api.decorators import register_marker
from .. import ParcelAggregation
from ...utils import logger
from ..parcel_aggregation import ParcelAggregation
from .falff_base import ALFFBase
@register_marker
class ALFFParcels(ALFFBase):
"""Class for computing fALFF/ALFF on parcels.
"""Class for ALFF / fALFF on parcels.
Parameters
----------
@ -23,6 +25,12 @@ class ALFFParcels(ALFFBase):
:func:`.list_parcellations`.
fractional : bool
Whether to compute fractional ALFF.
using : {"junifer", "afni"}
Implementation to use for computing ALFF:
* "junifer" : Use ``junifer``'s own ALFF implementation
* "afni" : Use AFNI's ``3dRSFC``
highpass : positive float, optional
The highpass cutoff frequency for the bandpass filter. If 0,
it will not apply a highpass filter (default 0.01).
@ -30,20 +38,17 @@ class ALFFParcels(ALFFBase):
The lowpass cutoff frequency for the bandpass filter (default 0.1).
tr : positive float, optional
The Repetition Time of the BOLD data. If None, will extract
the TR from NIFTI header (default None).
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 <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

111
junifer/markers/falff/falff_spheres.py
View file

@ -1,26 +1,36 @@
"""Provide class for computing fALFF on spheres."""
"""Provide class for ALFF / fALFF on spheres."""
# Authors: Federico Raimondo <f.raimondo@fz-juelich.de>
# Amir Omidvarnia <a.omidvarnia@fz-juelich.de>
# Kaustubh R. Patil <k.patil@fz-juelich.de>
# Synchon Mandal <s.mandal@fz-juelich.de>
# License: AGPL
from typing import Dict, List, Optional, Union
from typing import Any, Dict, List, Optional, Union
from ...api.decorators import register_marker
from .. import SphereAggregation
from ...utils import logger
from ..sphere_aggregation import SphereAggregation
from .falff_base import ALFFBase
@register_marker
class ALFFSpheres(ALFFBase):
"""Class for computing fALFF/ALFF on spheres.
"""Class for computing ALFF / fALFF on spheres.
Parameters
----------
coords : str
The name of the coordinates list to use. See
:func:`.list_coordinates` for options.
fractional : bool
Whether to compute fractional ALFF.
using : {"junifer", "afni"}
Implementation to use for computing ALFF:
* "junifer" : Use ``junifer``'s own ALFF implementation
* "afni" : Use AFNI's ``3dRSFC``
radius : float, optional
The radius of the sphere in mm. If None, the signal will be extracted
from a single voxel. See :class:`nilearn.maskers.NiftiSpheresMasker`
@ -28,8 +38,6 @@ class ALFFSpheres(ALFFBase):
allow_overlap : bool, optional
Whether to allow overlapping spheres. If False, an error is raised if
the spheres overlap (default is False).
fractional : bool
Whether to compute fractional ALFF.
highpass : positive float, optional
The highpass cutoff frequency for the bandpass filter. If 0,
it will not apply a highpass filter (default 0.01).
@ -37,20 +45,17 @@ class ALFFSpheres(ALFFBase):
The lowpass cutoff frequency for the bandpass filter (default 0.1).
tr : positive float, optional
The Repetition Time of the BOLD data. If None, will extract
the TR from NIFTI header (default None).
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 <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

267
junifer/markers/falff/tests/test_falff_estimator.py
View file

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

203
junifer/markers/falff/tests/test_falff_parcels.py
View file

@ -1,14 +1,14 @@
"""Provide test for parcel-aggregated (f)ALFF."""
"""Provide tests for ALFFParcels."""
# Authors: Federico Raimondo <f.raimondo@fz-juelich.de>
# Synchon Mandal <s.mandal@fz-juelich.de>
# License: AGPL
import logging
from pathlib import Path
import pytest
from numpy.testing import assert_array_equal
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

223
junifer/markers/falff/tests/test_falff_spheres.py
View file

@ -1,14 +1,14 @@
"""Provide test for sphere-aggregated (f)ALFF."""
"""Provide tests for ALFFSpheres."""
# Authors: Federico Raimondo <f.raimondo@fz-juelich.de>
# Synchon Mandal <s.mandal@fz-juelich.de>
# License: AGPL
import logging
from pathlib import Path
import pytest
from numpy.testing import assert_array_equal
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

192
junifer/markers/reho/_afni_reho.py Normal file
View file

@ -0,0 +1,192 @@
"""Provide class for computing regional homogeneity (ReHo) using AFNI."""
# Authors: Synchon Mandal <s.mandal@fz-juelich.de>
# License: AGPL
from functools import lru_cache
from pathlib import Path
from typing import (
TYPE_CHECKING,
ClassVar,
Dict,
List,
Optional,
Tuple,
Union,
)
import nibabel as nib
from ...pipeline import WorkDirManager
from ...pipeline.singleton import singleton
from ...utils import logger, run_ext_cmd
if TYPE_CHECKING:
from nibabel import Nifti1Image
@singleton
class AFNIReHo:
"""Class for computing ReHo using AFNI.
This class uses AFNI's 3dReHo to compute ReHo. It's designed as a singleton
with caching for efficient computation.
"""
_EXT_DEPENDENCIES: ClassVar[List[Dict[str, Union[str, List[str]]]]] = [
{
"name": "afni",
"commands": ["3dReHo", "3dAFNItoNIFTI"],
},
]
def __del__(self) -> None:
"""Terminate the class."""
# Clear the computation cache
logger.debug("Clearing cache for ReHo computation via AFNI")
self.compute.cache_clear()
@lru_cache(maxsize=None, typed=True)
def compute(
self,
data: "Nifti1Image",
nneigh: int = 27,
neigh_rad: Optional[float] = None,
neigh_x: Optional[float] = None,
neigh_y: Optional[float] = None,
neigh_z: Optional[float] = None,
box_rad: Optional[int] = None,
box_x: Optional[int] = None,
box_y: Optional[int] = None,
box_z: Optional[int] = None,
) -> Tuple["Nifti1Image", Path]:
"""Compute ReHo map.
Parameters
----------
data : 4D Niimg-like object
Images to process.
nneigh : {7, 19, 27}, optional
Number of voxels in the neighbourhood, inclusive. Can be:
* 7 : for facewise neighbours only
* 19 : for face- and edge-wise nieghbours
* 27 : for face-, edge-, and node-wise neighbors
(default 27).
neigh_rad : positive float, optional
The radius of a desired neighbourhood (default None).
neigh_x : positive float, optional
The semi-radius for x-axis of ellipsoidal volumes (default None).
neigh_y : positive float, optional
The semi-radius for y-axis of ellipsoidal volumes (default None).
neigh_z : positive float, optional
The semi-radius for z-axis of ellipsoidal volumes (default None).
box_rad : positive int, optional
The number of voxels outward in a given cardinal direction for a
cubic box centered on a given voxel (default None).
box_x : positive int, optional
The number of voxels for +/- x-axis of cuboidal volumes
(default None).
box_y : positive int, optional
The number of voxels for +/- y-axis of cuboidal volumes
(default None).
box_z : positive int, optional
The number of voxels for +/- z-axis of cuboidal volumes
(default None).
Returns
-------
Niimg-like object
The ReHo map as NIfTI.
pathlib.Path
The path to the ReHo map as NIfTI.
Notes
-----
For more information on the publication, please check [1]_ , and for
3dReHo help check:
https://afni.nimh.nih.gov/pub/dist/doc/program_help/3dReHo.html
Please note that that you cannot mix ``box_*`` and ``neigh_*``
arguments. The arguments are prioritized by their order in the function
signature.
As the process also depends on the conversion of AFNI files to NIFTI
via afni's 3dAFNItoNIFTI, the help for that can be found at:
https://afni.nimh.nih.gov/pub/dist/doc/program_help/3dAFNItoNIFTI.html
References
----------
.. [1] Taylor, P.A., & Saad, Z.S. (2013).
FATCAT: (An Efficient) Functional And Tractographic Connectivity
Analysis Toolbox.
Brain connectivity, Volume 3(5), Pages 523-35.
https://doi.org/10.1089/brain.2013.0154
"""
logger.debug("Creating cache for ReHo computation via AFNI")
# Create component-scoped tempdir
tempdir = WorkDirManager().get_tempdir(prefix="afni_reho")
# Save target data to a component-scoped tempfile
nifti_in_file_path = tempdir / "input.nii" # needs to be .nii
nib.save(data, nifti_in_file_path)
# Set 3dReHo command
reho_out_path_prefix = tempdir / "reho"
reho_cmd = [
"3dReHo",
f"-prefix {reho_out_path_prefix.resolve()}",
f"-inset {nifti_in_file_path.resolve()}",
]
# Check ellipsoidal / cuboidal volume arguments
if neigh_rad:
reho_cmd.append(f"-neigh_RAD {neigh_rad}")
elif neigh_x and neigh_y and neigh_z:
reho_cmd.extend(
[
f"-neigh_X {neigh_x}",
f"-neigh_Y {neigh_y}",
f"-neigh_Z {neigh_z}",
]
)
elif box_rad:
reho_cmd.append(f"-box_RAD {box_rad}")
elif box_x and box_y and box_z:
reho_cmd.extend(
[f"-box_X {box_x}", f"-box_Y {box_y}", f"-box_Z {box_z}"]
)
else:
reho_cmd.append(f"-nneigh {nneigh}")
# Call 3dReHo
run_ext_cmd(name="3dReHo", cmd=reho_cmd)
# Create element-scoped tempdir so that the ReHo map is
# available later as nibabel stores file path reference for
# loading on computation
element_tempdir = WorkDirManager().get_element_tempdir(
prefix="afni_reho"
)
# Convert afni to nifti
reho_afni_to_nifti_out_path = (
element_tempdir / "output.nii" # needs to be .nii
)
convert_cmd = [
"3dAFNItoNIFTI",
f"-prefix {reho_afni_to_nifti_out_path.resolve()}",
f"{reho_out_path_prefix}+tlrc.BRIK",
]
# Call 3dAFNItoNIFTI
run_ext_cmd(name="3dAFNItoNIFTI", cmd=convert_cmd)
# Load nifti
output_data = nib.load(reho_afni_to_nifti_out_path)
# Delete tempdir
WorkDirManager().delete_tempdir(tempdir)
return output_data, reho_afni_to_nifti_out_path # type: ignore

281
junifer/markers/reho/_junifer_reho.py Normal file
View file

@ -0,0 +1,281 @@
"""Provide class for computing regional homogeneity (ReHo) using junifer."""
# Authors: Synchon Mandal <s.mandal@fz-juelich.de>
# License: AGPL
from functools import lru_cache
from itertools import product
from pathlib import Path
from typing import (
TYPE_CHECKING,
ClassVar,
Set,
Tuple,
)
import nibabel as nib
import numpy as np
import scipy as sp
from nilearn import image as nimg
from nilearn import masking as nmask
from ...pipeline import WorkDirManager
from ...pipeline.singleton import singleton
from ...utils import logger, raise_error
if TYPE_CHECKING:
from nibabel import Nifti1Image
@singleton
class JuniferReHo:
"""Class for computing ReHo using junifer.
It's designed as a singleton with caching for efficient computation.
"""
_DEPENDENCIES: ClassVar[Set[str]] = {"numpy", "nilearn", "scipy"}
def __del__(self) -> None:
"""Terminate the class."""
# Clear the computation cache
logger.debug("Clearing cache for ReHo computation via junifer")
self.compute.cache_clear()
@lru_cache(maxsize=None, typed=True)
def compute(
self,
data: "Nifti1Image",
nneigh: int = 27,
) -> Tuple["Nifti1Image", Path]:
"""Compute ReHo map.
Parameters
----------
data : 4D Niimg-like object
Images to process.
nneigh : {7, 19, 27, 125}, optional
Number of voxels in the neighbourhood, inclusive. Can be:
* 7 : for facewise neighbours only
* 19 : for face- and edge-wise nieghbours
* 27 : for face-, edge-, and node-wise neighbors
* 125 : for 5x5 cuboidal volume
(default 27).
Returns
-------
Niimg-like object
The ReHo map as NIfTI.
pathlib.Path
The path to the ReHo map as NIfTI.
Raises
------
ValueError
If ``nneigh`` is invalid.
"""
valid_nneigh = (7, 19, 27, 125)
if nneigh not in valid_nneigh:
raise_error(
f"Invalid value for `nneigh`, should be one of: {valid_nneigh}"
)
logger.debug("Creating cache for ReHo computation via junifer")
# Get scan data
niimg_data = data.get_fdata()
# Get scan dimensions
n_x, n_y, n_z, _ = niimg_data.shape
# Get rank of every voxel across time series
ranks_niimg_data = sp.stats.rankdata(niimg_data, axis=-1)
# Initialize 3D array to store tied rank correction for every voxel
tied_rank_corrections = np.zeros((n_x, n_y, n_z), dtype=np.float64)
# Calculate tied rank correction for every voxel
for i_x, i_y, i_z in product(range(n_x), range(n_y), range(n_z)):
# Calculate tied rank count for every voxel across time series
_, tie_count = np.unique(
ranks_niimg_data[i_x, i_y, i_z, :],
return_counts=True,
)
# Calculate and store tied rank correction for every voxel across
# timeseries
tied_rank_corrections[i_x, i_y, i_z] = np.sum(
tie_count**3 - tie_count
)
# Initialize 3D array to store reho map
reho_map = np.ones((n_x, n_y, n_z), dtype=np.float32)
# TODO(synchon): this will give incorrect results if
# template doesn't match, hence needs to be changed
# after #299 is merged
# Calculate whole brain mask
mni152_whole_brain_mask = nmask.compute_brain_mask(
target_img=data,
threshold=0.5,
mask_type="whole-brain",
)
# Convert 0 / 1 array to bool
logical_mni152_whole_brain_mask = (
mni152_whole_brain_mask.get_fdata().astype(bool)
)
# Create mask cluster and set start and end indices
if nneigh in (7, 19, 27):
mask_cluster = np.ones((3, 3, 3))
if nneigh == 7:
mask_cluster[0, 0, 0] = 0
mask_cluster[0, 1, 0] = 0
mask_cluster[0, 2, 0] = 0
mask_cluster[0, 0, 1] = 0
mask_cluster[0, 2, 1] = 0
mask_cluster[0, 0, 2] = 0
mask_cluster[0, 1, 2] = 0
mask_cluster[0, 2, 2] = 0
mask_cluster[1, 0, 0] = 0
mask_cluster[1, 2, 0] = 0
mask_cluster[1, 0, 2] = 0
mask_cluster[1, 2, 2] = 0
mask_cluster[2, 0, 0] = 0
mask_cluster[2, 1, 0] = 0
mask_cluster[2, 2, 0] = 0
mask_cluster[2, 0, 1] = 0
mask_cluster[2, 2, 1] = 0
mask_cluster[2, 0, 2] = 0
mask_cluster[2, 1, 2] = 0
mask_cluster[2, 2, 2] = 0
elif nneigh == 19:
mask_cluster[0, 0, 0] = 0
mask_cluster[0, 2, 0] = 0
mask_cluster[2, 0, 0] = 0
mask_cluster[2, 2, 0] = 0
mask_cluster[0, 0, 2] = 0
mask_cluster[0, 2, 2] = 0
mask_cluster[2, 0, 2] = 0
mask_cluster[2, 2, 2] = 0
start_idx = 1
end_idx = 2
elif nneigh == 125:
mask_cluster = np.ones((5, 5, 5))
start_idx = 2
end_idx = 3
# Convert 0 / 1 array to bool
logical_mask_cluster = mask_cluster.astype(bool)
for i, j, k in product(
range(start_idx, n_x - (end_idx - 1)),
range(start_idx, n_y - (end_idx - 1)),
range(start_idx, n_z - (end_idx - 1)),
):
# Get mask only for neighbourhood
logical_neighbourhood_mni152_whole_brain_mask = (
logical_mni152_whole_brain_mask[
i - start_idx : i + end_idx,
j - start_idx : j + end_idx,
k - start_idx : k + end_idx,
]
)
# Perform logical AND to get neighbourhood mask;
# done to take care of brain boundaries
neighbourhood_mask = (
logical_mask_cluster
& logical_neighbourhood_mni152_whole_brain_mask
)
# Continue if voxel is restricted by mask
if neighbourhood_mask[1, 1, 1] == 0:
continue
# Get ranks for the neighbourhood
neighbourhood_ranks = ranks_niimg_data[
i - start_idx : i + end_idx,
j - start_idx : j + end_idx,
k - start_idx : k + end_idx,
:,
]
# Get tied ranks corrections for the neighbourhood
neighbourhood_tied_ranks_corrections = tied_rank_corrections[
i - start_idx : i + end_idx,
j - start_idx : j + end_idx,
k - start_idx : k + end_idx,
]
# Mask neighbourhood ranks
masked_neighbourhood_ranks = neighbourhood_ranks[
logical_mask_cluster, :
]
# Mask tied ranks corrections for the neighbourhood
masked_tied_rank_corrections = (
neighbourhood_tied_ranks_corrections[logical_mask_cluster]
)
# Calculate KCC
reho_map[i, j, k] = _kendall_w_reho(
timeseries_ranks=masked_neighbourhood_ranks,
tied_rank_corrections=masked_tied_rank_corrections,
)
# Create new image like target image
output_data = nimg.new_img_like(
ref_niimg=data,
data=reho_map,
copy_header=False,
)
# Create element-scoped tempdir so that the ReHo map is
# available later as nibabel stores file path reference for
# loading on computation
element_tempdir = WorkDirManager().get_element_tempdir(
prefix="junifer_reho"
)
output_path = element_tempdir / "output.nii.gz"
# Save computed data to file
nib.save(output_data, output_path)
return output_data, output_path # type: ignore
def _kendall_w_reho(
timeseries_ranks: np.ndarray, tied_rank_corrections: np.ndarray
) -> float:
"""Calculate Kendall's coefficient of concordance (KCC) for ReHo map.
..note:: This function should only be used to calculate KCC for a ReHo map.
For general use, check out ``junifer.stats.kendall_w``.
Parameters
----------
timeseries_ranks : 2D numpy.ndarray
A matrix of ranks of a subset subject's brain voxels.
tied_rank_corrections : 3D numpy.ndarray
A 3D array consisting of the tied rank corrections for the ranks
of a subset subject's brain voxels.
Returns
-------
float
Kendall's W (KCC) of the given timeseries matrix.
"""
m, n = timeseries_ranks.shape # annotators X items
numerator = (12 * np.sum(np.square(np.sum(timeseries_ranks, axis=0)))) - (
3 * m**2 * n * (n + 1) ** 2
)
denominator = (m**2 * n * (n**2 - 1)) - (m * np.sum(tied_rank_corrections))
if denominator == 0:
kcc = 1.0
else:
kcc = numerator / denominator
return kcc

84
junifer/markers/reho/reho_base.py
View file

@ -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

510
junifer/markers/reho/reho_estimator.py
View file

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

41
junifer/markers/reho/reho_parcels.py
View file

@ -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,

31
junifer/markers/reho/reho_spheres.py
View file

@ -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,

314
junifer/markers/reho/tests/test_reho_estimator.py
View file

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

188
junifer/markers/reho/tests/test_reho_parcels.py
View file

@ -1,73 +1,82 @@
"""Provide tests for ReHo on parcels."""
"""Provide tests for ReHoParcels."""
# Authors: Synchon Mandal <s.mandal@fz-juelich.de>
# License: AGPL
import logging
from pathlib import Path
import pytest
from nilearn import image as nimg
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,
)

195
junifer/markers/reho/tests/test_reho_spheres.py
View file

@ -1,73 +1,83 @@
"""Provide tests for ReHo on spheres."""
"""Provide tests for ReHoSpheres."""
# Authors: Synchon Mandal <s.mandal@fz-juelich.de>
# License: AGPL
import logging
from pathlib import Path
import pytest
from nilearn import image as nimg
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,
)

146
junifer/pipeline/pipeline_step_mixin.py
View file

@ -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

224
junifer/pipeline/tests/test_pipeline_step_mixin.py
View file

@ -1,11 +1,11 @@
"""Provide tests for pipeline mixin."""
"""Provide tests for PipelineStepMixin."""
# Authors: Federico Raimondo <f.raimondo@fz-juelich.de>
# Synchon Mandal <s.mandal@fz-juelich.de>
# License: AGPL
import warnings
from typing import ClassVar, Dict, List, Set, Union
from typing import ClassVar, Dict, List, Set, Type, Union
import pytest
@ -13,24 +13,13 @@ from junifer.pipeline.pipeline_step_mixin import PipelineStepMixin
from junifer.pipeline.utils import _check_afni
def test_PipelineStepMixin() -> None:
"""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({})

49
junifer/pipeline/utils.py
View file

@ -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(

45
junifer/preprocess/ants/ants_apply_transforms_warper.py
View file

@ -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,

68
junifer/preprocess/ants/tests/test_ants_apply_transforms_warper.py
View file

@ -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(

41
junifer/preprocess/bold_warper.py
View file

@ -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":

45
junifer/preprocess/fsl/apply_warper.py
View file

@ -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,

59
junifer/preprocess/fsl/tests/test_apply_warper.py
View file

@ -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")

6
junifer/preprocess/tests/test_bold_warper.py
View file

@ -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"