From 9fe66035990564ef9b3fa5d95aa979b3c2273a16 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Tue, 12 Mar 2024 12:19:59 -0400 Subject: [PATCH] Raise error if atlas affines don't match across participants (#1075) --- xcp_d/interfaces/bids.py | 120 ++++++++++++++++++++++++++++ xcp_d/tests/test_interfaces_bids.py | 84 +++++++++++++++++++ xcp_d/tests/test_utils_atlas.py | 54 ------------- xcp_d/utils/atlas.py | 65 --------------- xcp_d/utils/bids.py | 9 ++- xcp_d/workflows/connectivity.py | 45 ++--------- 6 files changed, 214 insertions(+), 163 deletions(-) create mode 100644 xcp_d/tests/test_interfaces_bids.py diff --git a/xcp_d/interfaces/bids.py b/xcp_d/interfaces/bids.py index 06e1c18d9..bfdc4bcf1 100644 --- a/xcp_d/interfaces/bids.py +++ b/xcp_d/interfaces/bids.py @@ -1,8 +1,12 @@ """Adapted interfaces from Niworkflows.""" +import os +import shutil from json import loads from pathlib import Path +import nibabel as nb +import numpy as np from bids.layout import Config from nipype import logging from nipype.interfaces.base import ( @@ -16,6 +20,8 @@ from niworkflows.interfaces.bids import DerivativesDataSink as BaseDerivativesDataSink from pkg_resources import resource_filename as pkgrf +from xcp_d.utils.bids import get_entity + # NOTE: Modified for xcpd's purposes xcp_d_spec = loads(Path(pkgrf("xcp_d", "data/xcp_d_bids_config.json")).read_text()) bids_config = Config.load("bids") @@ -190,3 +196,117 @@ def _run_interface(self, runtime): ) return runtime + + +class _CopyAtlasInputSpec(BaseInterfaceInputSpec): + name_source = traits.Str( + desc="The source file's name.", + mandatory=False, + ) + in_file = File( + exists=True, + desc="The atlas file to copy.", + mandatory=True, + ) + output_dir = Directory( + exists=True, + desc="The output directory.", + mandatory=True, + ) + atlas = traits.Str( + desc="The atlas name.", + mandatory=True, + ) + + +class _CopyAtlasOutputSpec(TraitedSpec): + out_file = File( + exists=True, + desc="The copied atlas file.", + ) + + +class CopyAtlas(SimpleInterface): + """Copy atlas file to output directory. + + Parameters + ---------- + name_source : :obj:`str` + The source name of the atlas file. + in_file : :obj:`str` + The atlas file to copy. + output_dir : :obj:`str` + The output directory. + atlas : :obj:`str` + The name of the atlas. + + Returns + ------- + out_file : :obj:`str` + The path to the copied atlas file. + + Notes + ----- + I can't use DerivativesDataSink because it has a problem with dlabel CIFTI files. + It gives the following error: + "AttributeError: 'Cifti2Header' object has no attribute 'set_data_dtype'" + + I can't override the CIFTI atlas's data dtype ahead of time because setting it to int8 or int16 + somehow converts all of the values in the data array to weird floats. + This could be a version-specific nibabel issue. + + I've also updated this function to handle JSON and TSV files as well. + """ + + input_spec = _CopyAtlasInputSpec + output_spec = _CopyAtlasOutputSpec + + def _run_interface(self, runtime): + output_dir = self.inputs.output_dir + in_file = self.inputs.in_file + name_source = self.inputs.name_source + atlas = self.inputs.atlas + + atlas_out_dir = os.path.join(output_dir, f"xcp_d/atlases/atlas-{atlas}") + + if in_file.endswith(".json"): + out_basename = f"atlas-{atlas}_dseg.json" + elif in_file.endswith(".tsv"): + out_basename = f"atlas-{atlas}_dseg.tsv" + else: + extension = ".nii.gz" if name_source.endswith(".nii.gz") else ".dlabel.nii" + space = get_entity(name_source, "space") + res = get_entity(name_source, "res") + den = get_entity(name_source, "den") + cohort = get_entity(name_source, "cohort") + + cohort_str = f"_cohort-{cohort}" if cohort else "" + res_str = f"_res-{res}" if res else "" + den_str = f"_den-{den}" if den else "" + if extension == ".dlabel.nii": + out_basename = f"space-{space}_atlas-{atlas}{den_str}{cohort_str}_dseg{extension}" + elif extension == ".nii.gz": + out_basename = f"space-{space}_atlas-{atlas}{res_str}{cohort_str}_dseg{extension}" + + os.makedirs(atlas_out_dir, exist_ok=True) + out_file = os.path.join(atlas_out_dir, out_basename) + + if out_file.endswith(".nii.gz") and os.path.isfile(out_file): + # Check that native-resolution atlas doesn't have a different resolution from the last + # run's atlas. + old_img = nb.load(out_file) + new_img = nb.load(in_file) + if not np.allclose(old_img.affine, new_img.affine): + raise ValueError( + f"Existing '{atlas}' atlas affine ({out_file}) is different from the input " + f"file affine ({in_file})." + ) + + # Don't copy the file if it exists, to prevent any race conditions between parallel + # processes. + if not os.path.isfile(out_file): + shutil.copyfile(in_file, out_file) + + self._results["out_file"] = out_file + + return runtime diff --git a/xcp_d/tests/test_interfaces_bids.py b/xcp_d/tests/test_interfaces_bids.py new file mode 100644 index 000000000..89bf26e09 --- /dev/null +++ b/xcp_d/tests/test_interfaces_bids.py @@ -0,0 +1,84 @@ +"""Tests for xcp_d.interfaces.bids.""" + +import os + +import pytest + +from xcp_d.interfaces import bids +from xcp_d.utils import atlas + + +def test_copy_atlas(tmp_path_factory): + """Test xcp_d.interfaces.bids.CopyAtlas.""" + tmpdir = tmp_path_factory.mktemp("test_copy_atlas") + os.makedirs(os.path.join(tmpdir, "xcp_d"), exist_ok=True) + + # NIfTI + atlas_file, _, _ = atlas.get_atlas_nifti("Gordon") + name_source = "sub-01_task-A_run-01_space-MNI152NLin2009cAsym_res-2_desc-z_bold.nii.gz" + copyatlas = bids.CopyAtlas( + name_source=name_source, in_file=atlas_file, output_dir=tmpdir, atlas="Y" + ) + result = copyatlas.run(cwd=tmpdir) + assert os.path.isfile(result.outputs.out_file) + assert ( + os.path.basename(result.outputs.out_file) + == "space-MNI152NLin2009cAsym_atlas-Y_res-2_dseg.nii.gz" + ) + + # Check that the NIfTI file raises an error if the resolution varies + # Gordon atlas is 1mm, HCP is 2mm + atlas_file_diff_affine, _, _ = atlas.get_atlas_nifti("HCP") + with pytest.raises(ValueError, match="is different from the input file affine"): + copyatlas = bids.CopyAtlas( + name_source=name_source, + in_file=atlas_file_diff_affine, + output_dir=tmpdir, + atlas="Y", + ) + copyatlas.run(cwd=tmpdir) + + # CIFTI + atlas_file, atlas_labels_file, atlas_metadata_file = atlas.get_atlas_cifti("Gordon") + name_source = "sub-01_task-imagery_run-01_space-fsLR_den-91k_desc-denoised_bold.dtseries.nii" + copyatlas = bids.CopyAtlas( + name_source=name_source, in_file=atlas_file, output_dir=tmpdir, atlas="Y" + ) + result = copyatlas.run(cwd=tmpdir) + assert os.path.isfile(result.outputs.out_file) + assert ( + os.path.basename(result.outputs.out_file) == "space-fsLR_atlas-Y_den-91k_dseg.dlabel.nii" + ) + + # TSV + name_source = "sub-01_task-imagery_run-01_space-fsLR_den-91k_desc-denoised_bold.dtseries.nii" + copyatlas = bids.CopyAtlas( + name_source=name_source, in_file=atlas_labels_file, output_dir=tmpdir, atlas="Y" + ) + result = copyatlas.run(cwd=tmpdir) + assert os.path.isfile(result.outputs.out_file) + assert os.path.basename(result.outputs.out_file) == "atlas-Y_dseg.tsv" + + # JSON + name_source = "sub-01_task-imagery_run-01_space-fsLR_den-91k_desc-denoised_bold.dtseries.nii" + copyatlas = bids.CopyAtlas( + name_source=name_source, in_file=atlas_metadata_file, output_dir=tmpdir, atlas="Y" + ) + result = copyatlas.run(cwd=tmpdir) + assert os.path.isfile(result.outputs.out_file) + assert os.path.basename(result.outputs.out_file) == "atlas-Y_dseg.json" + + # Ensure that out_file isn't overwritten if it already exists + fake_in_file = os.path.join(tmpdir, "fake.json") + with open(fake_in_file, "w") as fo: + fo.write("fake") + + copyatlas = bids.CopyAtlas( + name_source=name_source, in_file=fake_in_file, output_dir=tmpdir, atlas="Y" + ) + result = copyatlas.run(cwd=tmpdir) + assert os.path.isfile(result.outputs.out_file) + assert os.path.basename(result.outputs.out_file) == "atlas-Y_dseg.json" + # The file should not be overwritten, so the contents shouldn't be "fake" + with open(result.outputs.out_file, "r") as fo: + assert fo.read() != "fake" diff --git a/xcp_d/tests/test_utils_atlas.py b/xcp_d/tests/test_utils_atlas.py index 43ffdf765..f877fd709 100644 --- a/xcp_d/tests/test_utils_atlas.py +++ b/xcp_d/tests/test_utils_atlas.py @@ -45,57 +45,3 @@ def test_get_atlas_cifti(): with pytest.raises(FileNotFoundError, match="DNE"): atlas.get_atlas_cifti("tofail") - - -def test_copy_atlas(tmp_path_factory): - """Test xcp_d.utils.atlas.copy_atlas.""" - tmpdir = tmp_path_factory.mktemp("test_copy_atlas") - os.makedirs(os.path.join(tmpdir, "xcp_d"), exist_ok=True) - - # NIfTI - atlas_file, _, _ = atlas.get_atlas_nifti("Gordon") - name_source = "sub-01_task-A_run-01_space-MNI152NLin2009cAsym_res-2_desc-z_bold.nii.gz" - out_file = atlas.copy_atlas( - name_source=name_source, in_file=atlas_file, output_dir=tmpdir, atlas="Y" - ) - assert os.path.isfile(out_file) - assert os.path.basename(out_file) == "space-MNI152NLin2009cAsym_atlas-Y_res-2_dseg.nii.gz" - - # CIFTI - atlas_file, atlas_labels_file, atlas_metadata_file = atlas.get_atlas_cifti("Gordon") - name_source = "sub-01_task-imagery_run-01_space-fsLR_den-91k_desc-denoised_bold.dtseries.nii" - out_file = atlas.copy_atlas( - name_source=name_source, in_file=atlas_file, output_dir=tmpdir, atlas="Y" - ) - assert os.path.isfile(out_file) - assert os.path.basename(out_file) == "space-fsLR_atlas-Y_den-91k_dseg.dlabel.nii" - - # TSV - name_source = "sub-01_task-imagery_run-01_space-fsLR_den-91k_desc-denoised_bold.dtseries.nii" - out_file = atlas.copy_atlas( - name_source=name_source, in_file=atlas_labels_file, output_dir=tmpdir, atlas="Y" - ) - assert os.path.isfile(out_file) - assert os.path.basename(out_file) == "atlas-Y_dseg.tsv" - - # JSON - name_source = "sub-01_task-imagery_run-01_space-fsLR_den-91k_desc-denoised_bold.dtseries.nii" - out_file = atlas.copy_atlas( - name_source=name_source, in_file=atlas_metadata_file, output_dir=tmpdir, atlas="Y" - ) - assert os.path.isfile(out_file) - assert os.path.basename(out_file) == "atlas-Y_dseg.json" - - # Ensure that out_file isn't overwritten if it already exists - fake_in_file = os.path.join(tmpdir, "fake.json") - with open(fake_in_file, "w") as fo: - fo.write("fake") - - out_file = atlas.copy_atlas( - name_source=name_source, in_file=fake_in_file, output_dir=tmpdir, atlas="Y" - ) - assert os.path.isfile(out_file) - assert os.path.basename(out_file) == "atlas-Y_dseg.json" - # The file should not be overwritten, so the contents shouldn't be "fake" - with open(out_file, "r") as fo: - assert fo.read() != "fake" diff --git a/xcp_d/utils/atlas.py b/xcp_d/utils/atlas.py index 561973b83..ac1b70d60 100644 --- a/xcp_d/utils/atlas.py +++ b/xcp_d/utils/atlas.py @@ -155,68 +155,3 @@ def get_atlas_cifti(atlas): ) return atlas_file, atlas_labels_file, atlas_metadata_file - - -def copy_atlas(name_source, in_file, output_dir, atlas): - """Copy atlas file to output directory. - - Parameters - ---------- - name_source : :obj:`str` - The source name of the atlas file. - in_file : :obj:`str` - The atlas file to copy. - output_dir : :obj:`str` - The output directory. - atlas : :obj:`str` - The name of the atlas. - - Returns - ------- - out_file : :obj:`str` - The path to the copied atlas file. - - Notes - ----- - I can't use DerivativesDataSink because it has a problem with dlabel CIFTI files. - It gives the following error: - "AttributeError: 'Cifti2Header' object has no attribute 'set_data_dtype'" - - I can't override the CIFTI atlas's data dtype ahead of time because setting it to int8 or int16 - somehow converts all of the values in the data array to weird floats. - This could be a version-specific nibabel issue. - - I've also updated this function to handle JSON and TSV files as well. - """ - import os - import shutil - - from xcp_d.utils.bids import get_entity - - if in_file.endswith(".json"): - out_basename = f"atlas-{atlas}_dseg.json" - elif in_file.endswith(".tsv"): - out_basename = f"atlas-{atlas}_dseg.tsv" - else: - extension = ".nii.gz" if name_source.endswith(".nii.gz") else ".dlabel.nii" - space = get_entity(name_source, "space") - res = get_entity(name_source, "res") - den = get_entity(name_source, "den") - cohort = get_entity(name_source, "cohort") - - cohort_str = f"_cohort-{cohort}" if cohort else "" - res_str = f"_res-{res}" if res else "" - den_str = f"_den-{den}" if den else "" - if extension == ".dlabel.nii": - out_basename = f"space-{space}_atlas-{atlas}{den_str}{cohort_str}_dseg{extension}" - elif extension == ".nii.gz": - out_basename = f"space-{space}_atlas-{atlas}{res_str}{cohort_str}_dseg{extension}" - - atlas_out_dir = os.path.join(output_dir, f"xcp_d/atlases/atlas-{atlas}") - os.makedirs(atlas_out_dir, exist_ok=True) - out_file = os.path.join(atlas_out_dir, out_basename) - # Don't copy the file if it exists, to prevent any race conditions between parallel processes. - if not os.path.isfile(out_file): - shutil.copyfile(in_file, out_file) - - return out_file diff --git a/xcp_d/utils/bids.py b/xcp_d/utils/bids.py index 8230b3c37..0938f7024 100644 --- a/xcp_d/utils/bids.py +++ b/xcp_d/utils/bids.py @@ -317,11 +317,12 @@ def collect_data( # This probably works well for resolution (1 typically means 1x1x1, # 2 typically means 2x2x2, etc.), but probably doesn't work well for density. resolutions = layout.get_res(**queries["bold"]) - densities = layout.get_den(**queries["bold"]) - if len(resolutions) > 1: - queries["bold"]["resolution"] = resolutions[0] + if len(resolutions) >= 1: + # This will also select res-* when there are both res-* and native-resolution files. + queries["bold"]["res"] = resolutions[0] - if len(densities) > 1: + densities = layout.get_den(**queries["bold"]) + if len(densities) >= 1: queries["bold"]["den"] = densities[0] # Check for anatomical images, and determine if T2w xfms must be used. diff --git a/xcp_d/workflows/connectivity.py b/xcp_d/workflows/connectivity.py index 0535d4935..643651a0f 100644 --- a/xcp_d/workflows/connectivity.py +++ b/xcp_d/workflows/connectivity.py @@ -7,7 +7,7 @@ from niworkflows.engine.workflows import LiterateWorkflow as Workflow from xcp_d.interfaces.ants import ApplyTransforms -from xcp_d.interfaces.bids import DerivativesDataSink +from xcp_d.interfaces.bids import CopyAtlas, DerivativesDataSink from xcp_d.interfaces.connectivity import ( CiftiConnect, CiftiParcellate, @@ -21,12 +21,7 @@ CiftiCreateDenseFromTemplate, CiftiParcellateWorkbench, ) -from xcp_d.utils.atlas import ( - copy_atlas, - get_atlas_cifti, - get_atlas_nifti, - select_atlases, -) +from xcp_d.utils.atlas import get_atlas_cifti, get_atlas_nifti, select_atlases from xcp_d.utils.doc import fill_doc from xcp_d.utils.utils import get_std2bold_xfms @@ -214,21 +209,11 @@ def init_load_atlases_wf( ]) # fmt:skip ds_atlas = pe.MapNode( - Function( - function=copy_atlas, - input_names=[ - "name_source", - "in_file", - "output_dir", - "atlas", - ], - output_names=["out_file"], - ), + CopyAtlas(output_dir=output_dir), name="ds_atlas", iterfield=["in_file", "atlas"], run_without_submitting=True, ) - ds_atlas.inputs.output_dir = output_dir ds_atlas.inputs.atlas = atlases workflow.connect([ @@ -238,21 +223,11 @@ def init_load_atlases_wf( ]) # fmt:skip ds_atlas_labels_file = pe.MapNode( - Function( - function=copy_atlas, - input_names=[ - "name_source", - "in_file", - "output_dir", - "atlas", - ], - output_names=["out_file"], - ), + CopyAtlas(output_dir=output_dir), name="ds_atlas_labels_file", iterfield=["in_file", "atlas"], run_without_submitting=True, ) - ds_atlas_labels_file.inputs.output_dir = output_dir ds_atlas_labels_file.inputs.atlas = atlases workflow.connect([ @@ -262,21 +237,11 @@ def init_load_atlases_wf( ]) # fmt:skip ds_atlas_metadata = pe.MapNode( - Function( - function=copy_atlas, - input_names=[ - "name_source", - "in_file", - "output_dir", - "atlas", - ], - output_names=["out_file"], - ), + CopyAtlas(output_dir=output_dir), name="ds_atlas_metadata", iterfield=["in_file", "atlas"], run_without_submitting=True, ) - ds_atlas_metadata.inputs.output_dir = output_dir ds_atlas_metadata.inputs.atlas = atlases workflow.connect([