Skip to content

Commit

Permalink
Load spheres for anatomical workflow from TemplateFlow and BIDS deriv…
Browse files Browse the repository at this point in the history 
…atives (#1207)
  • Loading branch information
@tsalo
tsalo authored Jul 19, 2024
1 parent 4878b5b commit ad49779
Show file tree
Hide file tree
Showing 8 changed files with 187 additions and 203 deletions.
88 changes: 29 additions & 59 deletions xcp_d/interfaces/bids.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -50,11 +50,6 @@ class DerivativesDataSink(BaseDerivativesDataSink):


class _CollectRegistrationFilesInputSpec(BaseInterfaceInputSpec):
segmentation_dir = Directory(
exists=True,
required=True,
desc="Path to FreeSurfer or MCRIBS derivatives.",
)
software = traits.Enum(
"FreeSurfer",
"MCRIBS",
Expand All @@ -67,17 +62,9 @@ class _CollectRegistrationFilesInputSpec(BaseInterfaceInputSpec):
required=True,
desc="The hemisphere being used.",
)
participant_id = traits.Str(
required=True,
desc="Participant ID. Used to select the subdirectory of the FreeSurfer derivatives.",
)


class _CollectRegistrationFilesOutputSpec(TraitedSpec):
subject_sphere = File(
exists=True,
desc="Subject-space sphere.",
)
source_sphere = File(
exists=True,
desc="Source-space sphere (namely, fsaverage).",
Expand All @@ -99,29 +86,13 @@ class CollectRegistrationFiles(SimpleInterface):
output_spec = _CollectRegistrationFilesOutputSpec

def _run_interface(self, runtime):
import os

from templateflow.api import get as get_template

from xcp_d.data import load as load_data

hemisphere = self.inputs.hemisphere
hstr = f"{hemisphere.lower()}h"
participant_id = self.inputs.participant_id
if not participant_id.startswith("sub-"):
participant_id = f"sub-{participant_id}"

if self.inputs.software == "FreeSurfer":
# Find the subject's sphere in the FreeSurfer derivatives.
# TODO: Collect from the preprocessing derivatives if they're a compliant version.
# Namely, fMRIPrep >= 23.1.2, Nibabies >= 24.0.0a1.
self._results["subject_sphere"] = os.path.join(
self.inputs.segmentation_dir,
participant_id,
"surf",
f"{hstr}.sphere.reg",
)

# Load the fsaverage-164k sphere
# FreeSurfer: tpl-fsaverage_hemi-?_den-164k_sphere.surf.gii
self._results["source_sphere"] = str(
Expand Down Expand Up @@ -158,40 +129,39 @@ def _run_interface(self, runtime):
)

elif self.inputs.software == "MCRIBS":
# Find the subject's sphere in the MCRIBS derivatives.
# TODO: Collect from the preprocessing derivatives if they're a compliant version.
# Namely, fMRIPrep >= 23.1.2, Nibabies >= 24.0.0a1.
self._results["subject_sphere"] = os.path.join(
self.inputs.segmentation_dir,
participant_id,
"freesurfer",
participant_id,
"surf",
f"{hstr}.sphere.reg2",
)

# TODO: Collect from templateflow once it's uploaded.
# MCRIBS: tpl-fsaverage_hemi-?_den-41k_desc-reg_sphere.surf.gii
self._results["source_sphere"] = os.path.join(
self.inputs.segmentation_dir,
"templates_fsLR",
f"tpl-fsaverage_hemi-{hemisphere}_den-41k_desc-reg_sphere.surf.gii",
self._results["source_sphere"] = str(
get_template(
template="fsaverage",
space=None,
hemi=hemisphere,
density="41k",
desc=None,
suffix="sphere",
),
)

# TODO: Collect from templateflow once it's uploaded.
# MCRIBS: tpl-dHCP_space-fsaverage_hemi-?_den-41k_desc-reg_sphere.surf.gii
self._results["sphere_to_sphere"] = os.path.join(
self.inputs.segmentation_dir,
"templates_fsLR",
f"tpl-dHCP_space-fsaverage_hemi-{hemisphere}_den-41k_desc-reg_sphere.surf.gii",
self._results["sphere_to_sphere"] = str(
get_template(
template="dhcpAsym",
cohort="42",
space="fsaverage",
hemi=hemisphere,
density="41k",
desc="reg",
suffix="sphere",
),
)

# TODO: Collect from templateflow once it's uploaded.
# MCRIBS: tpl-dHCP_space-fsLR_hemi-?_den-32k_desc-week42_sphere.surf.gii
self._results["target_sphere"] = os.path.join(
self.inputs.segmentation_dir,
"templates_fsLR",
f"tpl-dHCP_space-fsLR_hemi-{hemisphere}_den-32k_desc-week42_sphere.surf.gii",
self._results["target_sphere"] = str(
get_template(
template="dhcpAsym",
cohort="42",
space=None,
hemi=hemisphere,
density="32k",
desc=None,
suffix="sphere",
),
)

return runtime
Expand Down
2 changes: 1 addition & 1 deletion xcp_d/interfaces/connectivity.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -231,6 +231,7 @@ class _TSVConnectOutputSpec(TraitedSpec):
def correlate_timeseries(timeseries, temporal_mask):
"""Correlate timeseries stored in a TSV file."""
timeseries_df = pd.read_table(timeseries)
correlations_exact = {}
if isdefined(temporal_mask):
censoring_df = pd.read_table(temporal_mask)

Expand All @@ -243,7 +244,6 @@ def correlate_timeseries(timeseries, temporal_mask):
censored_censoring_df = censoring_df.loc[censoring_df["framewise_displacement"] == 0]
censored_censoring_df.reset_index(drop=True, inplace=True)
exact_columns = [c for c in censoring_df.columns if c.startswith("exact_")]
correlations_exact = {}
for exact_column in exact_columns:
exact_timeseries_df = timeseries_df.loc[censored_censoring_df[exact_column] == 0]
exact_correlations_df = exact_timeseries_df.corr()
Expand Down
26 changes: 16 additions & 10 deletions xcp_d/tests/conftest.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -148,38 +148,44 @@ def pnc_data(datasets):
)
files["brain_mask_file"] = os.path.join(
func_dir,
"sub-1648798153_task-rest_space-MNI152NLin2009cAsym_res-2_desc-brain_mask.nii.gz",
"sub-1648798153_ses-PNC1_task-rest_space-MNI152NLin2009cAsym_res-2_desc-brain_mask.nii.gz",
)
files["confounds_file"] = os.path.join(
func_dir,
"sub-1648798153_task-rest_desc-confounds_timeseries.tsv",
"sub-1648798153_ses-PNC1_task-rest_desc-confounds_timeseries.tsv",
)
files["confounds_json"] = os.path.join(
func_dir,
"sub-1648798153_task-rest_desc-confounds_timeseries.json",
"sub-1648798153_ses-PNC1_task-rest_desc-confounds_timeseries.json",
)
files["anat_to_template_xfm"] = os.path.join(
anat_dir,
"sub-1648798153_from-T1w_to-MNI152NLin2009cAsym_mode-image_xfm.h5",
"sub-1648798153_ses-PNC1_acq-refaced_from-T1w_to-MNI152NLin2009cAsym_mode-image_xfm.h5",
)
files["template_to_anat_xfm"] = os.path.join(
anat_dir,
"sub-1648798153_from-MNI152NLin2009cAsym_to-T1w_mode-image_xfm.h5",
"sub-1648798153_ses-PNC1_acq-refaced_from-MNI152NLin2009cAsym_to-T1w_mode-image_xfm.h5",
)
files["boldref"] = os.path.join(
func_dir,
"sub-1648798153_task-rest_space-MNI152NLin2009cAsym_res-2_boldref.nii.gz",
"sub-1648798153_ses-PNC1_task-rest_space-MNI152NLin2009cAsym_res-2_boldref.nii.gz",
)
files["boldref_t1w"] = os.path.join(
func_dir,
"sub-1648798153_task-rest_space-T1w_boldref.nii.gz",
"sub-1648798153_ses-PNC1_task-rest_space-T1w_boldref.nii.gz",
)
files["t1w"] = os.path.join(anat_dir, "sub-1648798153_desc-preproc_T1w.nii.gz")
files["t1w"] = os.path.join(anat_dir, "sub-1648798153_acq-refaced_desc-preproc_T1w.nii.gz")
files["t1w_mni"] = os.path.join(
anat_dir,
"sub-1648798153_space-MNI152NLin2009cAsym_res-2_desc-preproc_T1w.nii.gz",
(
"sub-1648798153_ses-PNC1_acq-refaced_space-MNI152NLin2009cAsym_"
"res-2_desc-preproc_T1w.nii.gz"
),
)
files["anat_dseg"] = os.path.join(
anat_dir,
"sub-1648798153_ses-PNC1_acq-refaced_desc-aseg_dseg.nii.gz",
)
files["anat_dseg"] = os.path.join(anat_dir, "sub-1648798153_desc-aseg_dseg.nii.gz")

return files

Expand Down
68 changes: 23 additions & 45 deletions xcp_d/tests/test_utils_bids.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -126,60 +126,60 @@ def test_collect_mesh_data(datasets, tmp_path_factory):
"""Test collect_mesh_data."""
# Dataset without mesh files
layout = BIDSLayout(datasets["fmriprep_without_freesurfer"], validate=False)
mesh_available, standard_space_mesh, _ = xbids.collect_mesh_data(layout, "01")
mesh_available, standard_space_mesh, _, _ = xbids.collect_mesh_data(layout, "1648798153")
assert mesh_available is False
assert standard_space_mesh is False

# Dataset with native-space mesh files (one file matching each query)
layout = BIDSLayout(datasets["ds001419"], validate=False)
mesh_available, standard_space_mesh, _ = xbids.collect_mesh_data(layout, "01")
layout = BIDSLayout(datasets["pnc"], validate=False)
mesh_available, standard_space_mesh, _, _ = xbids.collect_mesh_data(layout, "1648798153")
assert mesh_available is True
assert standard_space_mesh is False

# Dataset with standard-space mesh files (one file matching each query)
std_mesh_dir = tmp_path_factory.mktemp("standard_mesh")
shutil.copyfile(
os.path.join(datasets["ds001419"], "dataset_description.json"),
os.path.join(datasets["pnc"], "dataset_description.json"),
std_mesh_dir / "dataset_description.json",
)
os.makedirs(std_mesh_dir / "sub-01/anat", exist_ok=True)
os.makedirs(std_mesh_dir / "sub-1648798153/ses-PNC1/anat", exist_ok=True)
files = [
"sub-01_space-fsLR_den-32k_hemi-L_pial.surf.gii",
"sub-01_space-fsLR_den-32k_hemi-L_white.surf.gii",
"sub-01_space-fsLR_den-32k_hemi-R_pial.surf.gii",
"sub-01_space-fsLR_den-32k_hemi-R_white.surf.gii",
"sub-1648798153_ses-PNC1_space-fsLR_den-32k_hemi-L_pial.surf.gii",
"sub-1648798153_ses-PNC1_space-fsLR_den-32k_hemi-L_white.surf.gii",
"sub-1648798153_ses-PNC1_space-fsLR_den-32k_hemi-R_pial.surf.gii",
"sub-1648798153_ses-PNC1_space-fsLR_den-32k_hemi-R_white.surf.gii",
]
for f in files:
(std_mesh_dir / "sub-01/anat").joinpath(f).touch()
(std_mesh_dir / "sub-1648798153/ses-PNC1/anat").joinpath(f).touch()

layout = BIDSLayout(std_mesh_dir, validate=False)
mesh_available, standard_space_mesh, _ = xbids.collect_mesh_data(layout, "01")
mesh_available, standard_space_mesh, _, _ = xbids.collect_mesh_data(layout, "1648798153")
assert mesh_available is True
assert standard_space_mesh is True

# Dataset with multiple files matching each query (raises an error)
bad_mesh_dir = tmp_path_factory.mktemp("standard_mesh")
shutil.copyfile(
os.path.join(datasets["ds001419"], "dataset_description.json"),
os.path.join(datasets["pnc"], "dataset_description.json"),
bad_mesh_dir / "dataset_description.json",
)
os.makedirs(bad_mesh_dir / "sub-01/anat", exist_ok=True)
os.makedirs(bad_mesh_dir / "sub-1648798153/ses-PNC1/anat", exist_ok=True)
files = [
"sub-01_space-fsLR_den-32k_hemi-L_pial.surf.gii",
"sub-01_space-fsLR_den-32k_hemi-L_white.surf.gii",
"sub-01_space-fsLR_den-32k_hemi-R_pial.surf.gii",
"sub-01_space-fsLR_den-32k_hemi-R_white.surf.gii",
"sub-01_acq-test_space-fsLR_den-32k_hemi-L_pial.surf.gii",
"sub-01_acq-test_space-fsLR_den-32k_hemi-L_white.surf.gii",
"sub-01_acq-test_space-fsLR_den-32k_hemi-R_pial.surf.gii",
"sub-01_acq-test_space-fsLR_den-32k_hemi-R_white.surf.gii",
"sub-1648798153_ses-PNC1_space-fsLR_den-32k_hemi-L_pial.surf.gii",
"sub-1648798153_ses-PNC1_space-fsLR_den-32k_hemi-L_white.surf.gii",
"sub-1648798153_ses-PNC1_space-fsLR_den-32k_hemi-R_pial.surf.gii",
"sub-1648798153_ses-PNC1_space-fsLR_den-32k_hemi-R_white.surf.gii",
"sub-1648798153_ses-PNC1_acq-test_space-fsLR_den-32k_hemi-L_pial.surf.gii",
"sub-1648798153_ses-PNC1_acq-test_space-fsLR_den-32k_hemi-L_white.surf.gii",
"sub-1648798153_ses-PNC1_acq-test_space-fsLR_den-32k_hemi-R_pial.surf.gii",
"sub-1648798153_ses-PNC1_acq-test_space-fsLR_den-32k_hemi-R_white.surf.gii",
]
for f in files:
(std_mesh_dir / "sub-01/anat").joinpath(f).touch()
(std_mesh_dir / "sub-1648798153/ses-PNC1/anat").joinpath(f).touch()

layout = BIDSLayout(std_mesh_dir, validate=False)
with pytest.raises(ValueError, match="More than one surface found"):
xbids.collect_mesh_data(layout, "01")
xbids.collect_mesh_data(layout, "1648798153")


def test_write_dataset_description(datasets, tmp_path_factory, caplog):
Expand Down Expand Up @@ -306,28 +306,6 @@ def test_get_tr(ds001419_data):
assert t_r == 3.0


def test_get_freesurfer_dir(datasets):
"""Test get_freesurfer_dir."""
with pytest.raises(NotADirectoryError, match="No FreeSurfer/MCRIBS derivatives found"):
xbids.get_freesurfer_dir(".")

fs_dir, software = xbids.get_freesurfer_dir(datasets["nibabies"])
assert os.path.isdir(fs_dir)
assert software == "FreeSurfer"

# Create fake FreeSurfer folder so there are two possible folders and it grabs the closest
tmp_fs_dir = os.path.join(datasets["nibabies"], "sourcedata/mcribs")
os.makedirs(tmp_fs_dir, exist_ok=True)
fs_dir, software = xbids.get_freesurfer_dir(datasets["nibabies"])
assert os.path.isdir(fs_dir)
assert software == "MCRIBS"
os.rmdir(tmp_fs_dir)

fs_dir, software = xbids.get_freesurfer_dir(datasets["pnc"])
assert os.path.isdir(fs_dir)
assert software == "FreeSurfer"


def test_get_entity(datasets):
"""Test get_entity."""
fname = os.path.join(datasets["ds001419"], "sub-01", "anat", "sub-01_desc-preproc_T1w.nii.gz")
Expand Down
Loading

0 comments on commit ad49779

Please sign in to comment.