Improve external dependency handling of PipelineStepMixin #311
29 changed files with 1696 additions and 2524 deletions
1
docs/changes/newsfragments/311.change
Normal file
1
docs/changes/newsfragments/311.change
Normal 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
1
docs/changes/newsfragments/311.enh
Normal 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
153
junifer/markers/falff/_afni_falff.py
Normal 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
142
junifer/markers/falff/_junifer_falff.py
Normal 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
|
||||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
@ -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 caplog.at_level(logging.DEBUG):
|
||||
with PartlyCloudyTestingDataGrabber() as dg:
|
||||
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=False,
|
||||
)
|
||||
python_values = marker.fit_transform(input_)["BOLD"]["data"]
|
||||
|
||||
assert marker.use_afni is False
|
||||
assert python_values.ndim == 2
|
||||
assert python_values.shape == (1, 100)
|
||||
# Initialize marker
|
||||
marker = ALFFParcels(
|
||||
parcellation=PARCELLATION,
|
||||
fractional=False,
|
||||
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 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_)
|
||||
element_data = DefaultDataReader().fit_transform(dg["sub-01"])
|
||||
# Update workdir to current test's tmp_path
|
||||
WorkDirManager().workdir = tmp_path
|
||||
marker_python = ALFFParcels(
|
||||
parcellation=_PARCELLATION,
|
||||
method="mean",
|
||||
use_afni=False,
|
||||
|
||||
# Initialize marker
|
||||
junifer_marker = ALFFParcels(
|
||||
parcellation=PARCELLATION,
|
||||
fractional=fractional,
|
||||
using="junifer",
|
||||
)
|
||||
python_values = marker_python.fit_transform(input_)["BOLD"]["data"]
|
||||
# Fit transform marker on data
|
||||
junifer_output = junifer_marker.fit_transform(element_data)
|
||||
# Get BOLD output
|
||||
junifer_output_bold = junifer_output["BOLD"]
|
||||
|
||||
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,
|
||||
# Initialize marker
|
||||
afni_marker = ALFFParcels(
|
||||
parcellation=PARCELLATION,
|
||||
fractional=fractional,
|
||||
using="afni",
|
||||
)
|
||||
afni_values = marker_afni.fit_transform(input_)["BOLD"]["data"]
|
||||
# Fit transform marker on data
|
||||
afni_output = afni_marker.fit_transform(element_data)
|
||||
# Get BOLD output
|
||||
afni_output_bold = afni_output["BOLD"]
|
||||
|
||||
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})")
|
||||
# Check for Pearson correlation coefficient
|
||||
r, _ = sp.stats.pearsonr(
|
||||
junifer_output_bold["data"][0],
|
||||
afni_output_bold["data"][0],
|
||||
)
|
||||
assert r > 0.99
|
||||
|
||||
|
||||
def test_ALFFParcels_storage(
|
||||
tmp_path: Path,
|
||||
) -> None:
|
||||
"""Test ALFFParcels storage.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
tmp_path : pathlib.Path
|
||||
The path to the test directory.
|
||||
|
||||
"""
|
||||
with PartlyCloudyTestingDataGrabber() as dg:
|
||||
# Use first subject
|
||||
input_ = dg["sub-01"]
|
||||
input_ = DefaultDataReader().fit_transform(input_)
|
||||
WorkDirManager().workdir = tmp_path
|
||||
marker = ALFFParcels(
|
||||
parcellation=_PARCELLATION,
|
||||
method="mean",
|
||||
use_afni=False,
|
||||
fractional=True,
|
||||
)
|
||||
storage = SQLiteFeatureStorage(tmp_path / "alff_parcels.sqlite")
|
||||
|
||||
# Fit transform marker on data with storage
|
||||
marker.fit_transform(
|
||||
input=input_,
|
||||
storage=storage,
|
||||
)
|
||||
|
|
|
|||
|
|
@ -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,96 +16,73 @@ 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 caplog.at_level(logging.DEBUG):
|
||||
with PartlyCloudyTestingDataGrabber() as dg:
|
||||
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=False,
|
||||
)
|
||||
python_values = marker.fit_transform(input_)["BOLD"]["data"]
|
||||
|
||||
assert marker.use_afni is False
|
||||
assert python_values.ndim == 2
|
||||
assert python_values.shape == (1, 6)
|
||||
# Initialize marker
|
||||
marker = ALFFSpheres(
|
||||
coords=COORDINATES,
|
||||
fractional=False,
|
||||
using="junifer",
|
||||
radius=5.0,
|
||||
)
|
||||
# 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 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"
|
||||
)
|
||||
def test_ALFFSpheres_afni(tmp_path: Path) -> 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"
|
||||
_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.
|
||||
def test_ALFFSpheres_comparison(tmp_path: Path, fractional: bool) -> None:
|
||||
"""Test ALFFSpheres using afni.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
|
|
@ -116,69 +93,37 @@ def test_ALFFSpheres_python_vs_afni(
|
|||
|
||||
"""
|
||||
with PartlyCloudyTestingDataGrabber() as dg:
|
||||
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_python = ALFFSpheres(
|
||||
coords=_COORDINATES,
|
||||
radius=5,
|
||||
method="mean",
|
||||
use_afni=False,
|
||||
|
||||
# Initialize marker
|
||||
junifer_marker = ALFFSpheres(
|
||||
coords=COORDINATES,
|
||||
fractional=fractional,
|
||||
using="junifer",
|
||||
radius=5.0,
|
||||
)
|
||||
python_values = marker_python.fit_transform(input_)["BOLD"]["data"]
|
||||
# Fit transform marker on data
|
||||
junifer_output = junifer_marker.fit_transform(element_data)
|
||||
# Get BOLD output
|
||||
junifer_output_bold = junifer_output["BOLD"]
|
||||
|
||||
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,
|
||||
# Initialize marker
|
||||
afni_marker = ALFFSpheres(
|
||||
coords=COORDINATES,
|
||||
fractional=fractional,
|
||||
using="afni",
|
||||
radius=5.0,
|
||||
)
|
||||
afni_values = marker_afni.fit_transform(input_)["BOLD"]["data"]
|
||||
# Fit transform marker on data
|
||||
afni_output = afni_marker.fit_transform(element_data)
|
||||
# Get BOLD output
|
||||
afni_output_bold = afni_output["BOLD"]
|
||||
|
||||
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})")
|
||||
# Check for Pearson correlation coefficient
|
||||
r, _ = sp.stats.pearsonr(
|
||||
junifer_output_bold["data"][0],
|
||||
afni_output_bold["data"][0],
|
||||
)
|
||||
assert r > 0.99
|
||||
|
||||
|
||||
def test_ALFFSpheres_storage(
|
||||
tmp_path: Path,
|
||||
) -> None:
|
||||
"""Test ALFFSpheres storage.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
tmp_path : pathlib.Path
|
||||
The path to the test directory.
|
||||
|
||||
"""
|
||||
with PartlyCloudyTestingDataGrabber() as dg:
|
||||
# Use first subject
|
||||
input_ = dg["sub-01"]
|
||||
input_ = DefaultDataReader().fit_transform(input_)
|
||||
WorkDirManager().workdir = tmp_path
|
||||
marker = ALFFSpheres(
|
||||
coords=_COORDINATES,
|
||||
radius=5,
|
||||
method="mean",
|
||||
use_afni=False,
|
||||
fractional=True,
|
||||
)
|
||||
storage = SQLiteFeatureStorage(tmp_path / "alff_parcels.sqlite")
|
||||
|
||||
# Fit transform marker on data with storage
|
||||
marker.fit_transform(
|
||||
input=input_,
|
||||
storage=storage,
|
||||
)
|
||||
|
|
|
|||
192
junifer/markers/reho/_afni_reho.py
Normal file
192
junifer/markers/reho/_afni_reho.py
Normal 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
281
junifer/markers/reho/_junifer_reho.py
Normal 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
|
||||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
@ -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,
|
||||
|
|
|
|||
|
|
@ -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,
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
@ -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 caplog.at_level(logging.DEBUG):
|
||||
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
|
||||
reho_parcels_marker = ReHoParcels(parcellation=PARCELLATION)
|
||||
# Fit transform marker on data
|
||||
reho_parcels_output = reho_parcels_marker.fit_transform(
|
||||
{
|
||||
"BOLD": {
|
||||
"path": "/tmp",
|
||||
"data": fmri_img,
|
||||
"meta": {},
|
||||
"space": "MNI",
|
||||
}
|
||||
}
|
||||
)
|
||||
# Get BOLD output
|
||||
reho_parcels_output_bold = reho_parcels_output["BOLD"]
|
||||
# Assert BOLD output keys
|
||||
assert "data" in reho_parcels_output_bold
|
||||
assert "col_names" in reho_parcels_output_bold
|
||||
|
||||
reho_parcels_output_bold_data = reho_parcels_output_bold["data"]
|
||||
# Initialize marker
|
||||
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 reho_parcels_output_bold_data.ndim == 2
|
||||
assert 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
|
||||
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,
|
||||
)
|
||||
|
|
|
|||
|
|
@ -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 caplog.at_level(logging.DEBUG):
|
||||
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
|
||||
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",
|
||||
}
|
||||
}
|
||||
marker = ReHoSpheres(
|
||||
coords=COORDINATES, using="junifer", radius=10.0
|
||||
)
|
||||
# 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
|
||||
# Fit transform marker on data
|
||||
output = marker.fit_transform(element_data)
|
||||
|
||||
reho_spheres_output_bold_data = reho_spheres_output_bold["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 reho_spheres_output_bold_data.ndim == 2
|
||||
assert 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
|
||||
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,
|
||||
)
|
||||
|
|
|
|||
|
|
@ -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.
|
||||
|
||||
"""
|
||||
|
||||
def _check_dependencies(obj) -> None:
|
||||
"""Check obj._DEPENDENCIES.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
obj : object
|
||||
Object to check _DEPENDENCIES of.
|
||||
|
||||
Raises
|
||||
------
|
||||
ValueError
|
||||
If the pipeline step object is missing dependencies required for
|
||||
its working or if the input does not have the required data.
|
||||
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(self, "_DEPENDENCIES"):
|
||||
if hasattr(obj, "_DEPENDENCIES"):
|
||||
# Check if dependencies are importable
|
||||
for dependency in self._DEPENDENCIES: # type: ignore
|
||||
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())
|
||||
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}.",
|
||||
msg=(
|
||||
f"{dependencies_not_found} are not installed but are "
|
||||
f"required for using {obj.__class__.__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_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
|
||||
|
||||
|
|
|
|||
|
|
@ -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({})
|
||||
|
|
|
|||
|
|
@ -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(
|
||||
msg=(
|
||||
f"{name} is not installed but is "
|
||||
"required by one of the pipeline steps."
|
||||
"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(
|
||||
|
|
|
|||
|
|
@ -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,
|
||||
|
|
|
|||
|
|
@ -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(
|
||||
|
|
|
|||
|
|
@ -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":
|
||||
|
|
|
|||
|
|
@ -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,
|
||||
|
|
|
|||
|
|
@ -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")
|
||||
|
|
|
|||
|
|
@ -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"
|
||||
|
||||
|
||||
|
|
|
|||
Loading…
Reference in a new issue