From ad4977989910589dda5e2caa517bc2d1b06700d0 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Fri, 19 Jul 2024 16:48:10 -0400 Subject: [PATCH] Load spheres for anatomical workflow from TemplateFlow and BIDS derivatives (#1207) --- xcp_d/interfaces/bids.py | 88 ++++++++---------------- xcp_d/interfaces/connectivity.py | 2 +- xcp_d/tests/conftest.py | 26 ++++--- xcp_d/tests/test_utils_bids.py | 68 +++++++----------- xcp_d/tests/test_workflows_anatomical.py | 48 +++++++++---- xcp_d/utils/bids.py | 67 ++++++++++++------ xcp_d/workflows/anatomical.py | 81 ++++++++-------------- xcp_d/workflows/base.py | 10 ++- 8 files changed, 187 insertions(+), 203 deletions(-) diff --git a/xcp_d/interfaces/bids.py b/xcp_d/interfaces/bids.py index efb621acd..433e5d34e 100644 --- a/xcp_d/interfaces/bids.py +++ b/xcp_d/interfaces/bids.py @@ -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", @@ -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).", @@ -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( @@ -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 diff --git a/xcp_d/interfaces/connectivity.py b/xcp_d/interfaces/connectivity.py index d1f2eeba7..b52bc9754 100644 --- a/xcp_d/interfaces/connectivity.py +++ b/xcp_d/interfaces/connectivity.py @@ -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) @@ -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() diff --git a/xcp_d/tests/conftest.py b/xcp_d/tests/conftest.py index 4075fbbf1..d6d3f7922 100644 --- a/xcp_d/tests/conftest.py +++ b/xcp_d/tests/conftest.py @@ -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 diff --git a/xcp_d/tests/test_utils_bids.py b/xcp_d/tests/test_utils_bids.py index c8732f716..cc66431ea 100644 --- a/xcp_d/tests/test_utils_bids.py +++ b/xcp_d/tests/test_utils_bids.py @@ -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): @@ -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") diff --git a/xcp_d/tests/test_workflows_anatomical.py b/xcp_d/tests/test_workflows_anatomical.py index 6127a55ba..41869d394 100644 --- a/xcp_d/tests/test_workflows_anatomical.py +++ b/xcp_d/tests/test_workflows_anatomical.py @@ -15,29 +15,48 @@ def surface_files(datasets, tmp_path_factory): """Collect real and fake surface files to test the anatomical workflow.""" tmpdir = tmp_path_factory.mktemp("surface_files") - anat_dir = os.path.join(datasets["ds001419"], "sub-01", "anat") + anat_dir = os.path.join(datasets["pnc"], "sub-1648798153", "ses-PNC1", "anat") files = { - "native_lh_pial": os.path.join(anat_dir, "sub-01_hemi-L_pial.surf.gii"), - "native_lh_wm": os.path.join(anat_dir, "sub-01_hemi-L_smoothwm.surf.gii"), - "native_rh_pial": os.path.join(anat_dir, "sub-01_hemi-R_pial.surf.gii"), - "native_rh_wm": os.path.join(anat_dir, "sub-01_hemi-R_smoothwm.surf.gii"), + "native_lh_pial": os.path.join( + anat_dir, "sub-1648798153_ses-PNC1_acq-refaced_hemi-L_pial.surf.gii" + ), + "native_lh_wm": os.path.join( + anat_dir, "sub-1648798153_ses-PNC1_acq-refaced_hemi-L_white.surf.gii" + ), + "native_rh_pial": os.path.join( + anat_dir, "sub-1648798153_ses-PNC1_acq-refaced_hemi-R_pial.surf.gii" + ), + "native_rh_wm": os.path.join( + anat_dir, "sub-1648798153_ses-PNC1_acq-refaced_hemi-R_white.surf.gii" + ), } final_files = files.copy() for fref, fpath in files.items(): std_fref = fref.replace("native_", "fsLR_") std_fname = os.path.basename(fpath) - std_fname = std_fname.replace("sub-01_", "sub-01_space-fsLR_den-32k_") + std_fname = std_fname.replace( + "sub-1648798153_ses-PNC1_acq-refaced_", + "sub-1648798153_ses-PNC1_acq-refaced_space-fsLR_den-32k_", + ) std_fpath = os.path.join(tmpdir, std_fname) shutil.copyfile(fpath, std_fpath) final_files[std_fref] = std_fpath + final_files["lh_subject_sphere"] = os.path.join( + anat_dir, + "sub-1648798153_ses-PNC1_acq-refaced_hemi-L_desc-reg_sphere.surf.gii", + ) + final_files["rh_subject_sphere"] = os.path.join( + anat_dir, + "sub-1648798153_ses-PNC1_acq-refaced_hemi-R_desc-reg_sphere.surf.gii", + ) + return final_files def test_warp_surfaces_to_template_wf( - datasets, - ds001419_data, + pnc_data, surface_files, tmp_path_factory, ): @@ -47,12 +66,9 @@ def test_warp_surfaces_to_template_wf( """ tmpdir = tmp_path_factory.mktemp("test_warp_surfaces_to_template_wf") - subject_id = "01" - wf = anatomical.init_warp_surfaces_to_template_wf( - fmri_dir=datasets["ds001419"], - subject_id=subject_id, output_dir=tmpdir, + software="FreeSurfer", omp_nthreads=1, ) @@ -60,15 +76,17 @@ def test_warp_surfaces_to_template_wf( wf.inputs.inputnode.rh_pial_surf = surface_files["native_rh_pial"] wf.inputs.inputnode.lh_wm_surf = surface_files["native_lh_wm"] wf.inputs.inputnode.rh_wm_surf = surface_files["native_rh_wm"] + wf.inputs.inputnode.lh_subject_sphere = surface_files["lh_subject_sphere"] + wf.inputs.inputnode.rh_subject_sphere = surface_files["rh_subject_sphere"] # transforms (only used if warp_to_standard is True) - wf.inputs.inputnode.anat_to_template_xfm = ds001419_data["anat_to_template_xfm"] - wf.inputs.inputnode.template_to_anat_xfm = ds001419_data["template_to_anat_xfm"] + wf.inputs.inputnode.anat_to_template_xfm = pnc_data["anat_to_template_xfm"] + wf.inputs.inputnode.template_to_anat_xfm = pnc_data["template_to_anat_xfm"] wf.base_dir = tmpdir wf.run() # All of the possible fsLR surfaces should be available. - out_anat_dir = os.path.join(tmpdir, "sub-01", "anat") + out_anat_dir = os.path.join(tmpdir, "sub-1648798153", "ses-PNC1", "anat") for key, filename in surface_files.items(): if "fsLR" in key: out_fname = os.path.basename(filename) diff --git a/xcp_d/utils/bids.py b/xcp_d/utils/bids.py index 26820c022..904a06044 100644 --- a/xcp_d/utils/bids.py +++ b/xcp_d/utils/bids.py @@ -318,7 +318,6 @@ def collect_data( queries["anat_to_template_xfm"]["from"] = "T2w" # Nibabies may include space-T2w for some derivatives queries["anat_dseg"]["space"] = ["T2w", None] - queries["anat_brainmask"]["space"] = ["T2w", None] # Search for the files. subj_data = { @@ -373,6 +372,8 @@ def collect_mesh_data(layout, participant_label): True if surface mesh files (pial and smoothwm) were found. False if they were not. standard_space_mesh : :obj:`bool` True if standard-space (fsLR) surface mesh files were found. False if they were not. + software : {"MCRIBS", "FreeSurfer"} + The software used to generate the surfaces. mesh_files : :obj:`dict` Dictionary of surface file identifiers and their paths. If the surface files weren't found, then the paths will be Nones. @@ -381,16 +382,18 @@ def collect_mesh_data(layout, participant_label): # The base surfaces can be used to generate the derived surfaces. # The base surfaces may be in native or standard space. base_queries = { - "pial_surf": "pial", - "wm_surf": ["smoothwm", "white"], - } - query_extras = { - "space": "fsLR", - "den": "32k", + "pial_surf": { + "desc": None, + "suffix": "pial", + }, + "wm_surf": { + "desc": None, + "suffix": ["smoothwm", "white"], + }, } standard_space_mesh = True - for name, suffixes in base_queries.items(): + for name, query in base_queries.items(): # First, try to grab the first base surface file in standard space. # If it's not available, switch to native T1w-space data. for hemisphere in ["L", "R"]: @@ -399,11 +402,12 @@ def collect_mesh_data(layout, participant_label): subject=participant_label, datatype="anat", hemi=hemisphere, - desc=None, - suffix=suffixes, + space="fsLR", + den="32k", extension=".surf.gii", - **query_extras, + **query, ) + if len(temp_files) == 0: LOGGER.info("No standard-space surfaces found.") standard_space_mesh = False @@ -411,23 +415,29 @@ def collect_mesh_data(layout, participant_label): LOGGER.warning(f"{name}: More than one standard-space surface found.") # Now that we know if there are standard-space surfaces available, we can grab the files. + query_extras = {} if not standard_space_mesh: query_extras = { "space": None, } + # We need the subject spheres if we're using native-space surfaces. + base_queries["subject_sphere"] = { + "space": None, + "desc": "reg", + "suffix": "sphere", + } initial_mesh_files = {} queries = {} - for name, suffixes in base_queries.items(): + for name, query in base_queries.items(): for hemisphere in ["L", "R"]: key = f"{hemisphere.lower()}h_{name}" queries[key] = { "subject": participant_label, "datatype": "anat", "hemi": hemisphere, - "desc": None, - "suffix": suffixes, "extension": ".surf.gii", + **query, **query_extras, } initial_mesh_files[key] = layout.get(return_type="file", **queries[key]) @@ -451,12 +461,31 @@ def collect_mesh_data(layout, participant_label): f"Query: {queries[dtype]}" ) + mesh_files["lh_subject_sphere"] = mesh_files.get("lh_subject_sphere", None) + mesh_files["rh_subject_sphere"] = mesh_files.get("rh_subject_sphere", None) + + # Check for *_space-dhcpAsym_desc-reg_sphere.surf.gii + # If we find it, we assume segmentation was done with MCRIBS. Otherwise, assume FreeSurfer. + dhcp_file = layout.get( + return_type="file", + datatype="anat", + subject=participant_label, + hemi="L", + space="dhcpAsym", + desc="reg", + suffix="sphere", + extension=".surf.gii", + ) + software = "MCRIBS" if bool(len(dhcp_file)) else "FreeSurfer" + LOGGER.log( 25, f"Collected mesh files:\n{yaml.dump(mesh_files, default_flow_style=False, indent=4)}", ) + if mesh_available: + LOGGER.log(25, f"Assuming segmentation was performed with {software}.") - return mesh_available, standard_space_mesh, mesh_files + return mesh_available, standard_space_mesh, software, mesh_files @fill_doc @@ -866,7 +895,7 @@ def _get_tr(img): return img.header.get_zooms()[-1] -def get_freesurfer_dir(fmri_dir): +def get_segmentation_software(fmri_dir): """Find FreeSurfer or MCRIBS derivatives associated with preprocessing pipeline. NOTE: This is a Node function. @@ -878,9 +907,7 @@ def get_freesurfer_dir(fmri_dir): Returns ------- - seg_path : :obj:`str` - Path to FreeSurfer or MCRIBS derivatives. - seg + software Raises ------ @@ -925,7 +952,7 @@ def get_freesurfer_dir(fmri_dir): f"{software} derivatives associated with {desc} preprocessing derivatives found " f"at {pattern}" ) - return pattern, software + return software # Otherwise, continue to the next pattern diff --git a/xcp_d/workflows/anatomical.py b/xcp_d/workflows/anatomical.py index 38a8b2494..77f3c1c59 100644 --- a/xcp_d/workflows/anatomical.py +++ b/xcp_d/workflows/anatomical.py @@ -1,7 +1,7 @@ # emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- # vi: set ft=python sts=4 ts=4 sw=4 et: """Anatomical post-processing workflows.""" -from nipype import Function, logging +from nipype import logging from nipype.interfaces import utility as niu from nipype.interfaces.ants import CompositeTransformUtil # MB from nipype.interfaces.freesurfer import MRIsConvert @@ -28,7 +28,6 @@ SurfaceGenerateInflated, SurfaceSphereProjectUnproject, ) -from xcp_d.utils.bids import get_freesurfer_dir from xcp_d.utils.doc import fill_doc from xcp_d.utils.utils import list_to_str from xcp_d.workflows.execsummary import ( @@ -287,12 +286,12 @@ def init_postprocess_anat_wf( @fill_doc def init_postprocess_surfaces_wf( - subject_id, mesh_available, standard_space_mesh, morphometry_files, t1w_available, t2w_available, + software, name="postprocess_surfaces_wf", ): """Postprocess surfaces. @@ -324,18 +323,17 @@ def init_postprocess_surfaces_wf( with mock_config(): wf = init_postprocess_surfaces_wf( - subject_id="01", mesh_available=True, standard_space_mesh=False, morphometry_files=[], t1w_available=True, t2w_available=True, + software="FreeSurfer", name="postprocess_surfaces_wf", ) Parameters ---------- - subject_id mesh_available : bool standard_space_mesh : bool morphometry_files : list of str @@ -343,6 +341,8 @@ def init_postprocess_surfaces_wf( True if a T1w image is available. t2w_available : bool True if a T2w image is available. + software : {"MCRIBS", "FreeSurfer"} + The software used to generate the surfaces. %(name)s Default is "postprocess_surfaces_wf". @@ -356,6 +356,7 @@ def init_postprocess_surfaces_wf( %(template_to_anat_xfm)s lh_pial_surf, rh_pial_surf lh_wm_surf, rh_wm_surf + lh_subject_sphere, rh_subject_sphere sulcal_depth sulcal_curv cortical_thickness @@ -365,7 +366,6 @@ def init_postprocess_surfaces_wf( """ workflow = Workflow(name=name) - fmri_dir = config.execution.fmri_dir abcc_qc = config.workflow.abcc_qc process_surfaces = config.workflow.process_surfaces output_dir = config.execution.xcp_d_dir @@ -378,6 +378,8 @@ def init_postprocess_surfaces_wf( "t2w", "anat_to_template_xfm", "template_to_anat_xfm", + "lh_subject_sphere", + "rh_subject_sphere", "lh_pial_surf", "rh_pial_surf", "lh_wm_surf", @@ -485,14 +487,15 @@ def init_postprocess_surfaces_wf( workflow.__desc__ += " fsnative-space surfaces were then warped to fsLR space." # Mesh files are in fsnative and must be warped to fsLR. warp_surfaces_to_template_wf = init_warp_surfaces_to_template_wf( - fmri_dir=fmri_dir, - subject_id=subject_id, output_dir=output_dir, + software=software, omp_nthreads=omp_nthreads, name="warp_surfaces_to_template_wf", ) workflow.connect([ (inputnode, warp_surfaces_to_template_wf, [ + ("lh_subject_sphere", "inputnode.lh_subject_sphere"), + ("rh_subject_sphere", "inputnode.rh_subject_sphere"), ("lh_pial_surf", "inputnode.lh_pial_surf"), ("rh_pial_surf", "inputnode.rh_pial_surf"), ("lh_wm_surf", "inputnode.lh_wm_surf"), @@ -531,9 +534,8 @@ def init_postprocess_surfaces_wf( @fill_doc def init_warp_surfaces_to_template_wf( - fmri_dir, - subject_id, output_dir, + software, omp_nthreads, name="warp_surfaces_to_template_wf", ): @@ -547,18 +549,17 @@ def init_warp_surfaces_to_template_wf( from xcp_d.workflows.anatomical import init_warp_surfaces_to_template_wf wf = init_warp_surfaces_to_template_wf( - fmri_dir=".", - subject_id="01", output_dir=".", + software="FreeSurfer", omp_nthreads=1, name="warp_surfaces_to_template_wf", ) Parameters ---------- - %(fmri_dir)s - %(subject_id)s %(output_dir)s + software : {"MCRIBS", "FreeSurfer"} + The software used to generate the surfaces. %(omp_nthreads)s %(name)s Default is "warp_surfaces_to_template_wf". @@ -573,6 +574,8 @@ def init_warp_surfaces_to_template_wf( The template in question should match the volumetric space of the BOLD CIFTI files being processed by the main xcpd workflow. For example, MNI152NLin6Asym for fsLR-space CIFTIs. + lh_subject_sphere, rh_subject_sphere : :obj:`str` + Left- and right-hemisphere sphere registration files. lh_pial_surf, rh_pial_surf : :obj:`str` Left- and right-hemisphere pial surface files in fsnative space. lh_wm_surf, rh_wm_surf : :obj:`str` @@ -594,6 +597,8 @@ def init_warp_surfaces_to_template_wf( "anat_to_template_xfm", "template_to_anat_xfm", # surfaces + "lh_subject_sphere", + "rh_subject_sphere", "lh_pial_surf", "rh_pial_surf", "lh_wm_surf", @@ -617,16 +622,6 @@ def init_warp_surfaces_to_template_wf( ) # Warp the surfaces to space-fsLR, den-32k. - get_freesurfer_dir_node = pe.Node( - Function( - function=get_freesurfer_dir, - input_names=["fmri_dir"], - output_names=["freesurfer_path", "segmentation_software"], - ), - name="get_freesurfer_dir_node", - ) - get_freesurfer_dir_node.inputs.fmri_dir = fmri_dir - # First, we create the Connectome WorkBench-compatible transform files. update_xfm_wf = init_ants_xfm_to_fsl_wf( mem_gb=1, @@ -658,16 +653,15 @@ def init_warp_surfaces_to_template_wf( ]) # fmt:skip apply_transforms_wf = init_warp_one_hemisphere_wf( - participant_id=subject_id, hemisphere=hemi, + software=software, mem_gb=2, omp_nthreads=omp_nthreads, name=f"{hemi_label}_apply_transforms_wf", ) workflow.connect([ - (get_freesurfer_dir_node, apply_transforms_wf, [ - ("freesurfer_path", "inputnode.freesurfer_path"), - ("segmentation_software", "inputnode.segmentation_software"), + (inputnode, apply_transforms_wf, [ + (f"{hemi_label}_subject_sphere", "inputnode.subject_sphere"), ]), (update_xfm_wf, apply_transforms_wf, [ ("outputnode.merged_warpfield", "inputnode.merged_warpfield"), @@ -1099,8 +1093,8 @@ def init_ants_xfm_to_fsl_wf(mem_gb, omp_nthreads, name="ants_xfm_to_fsl_wf"): @fill_doc def init_warp_one_hemisphere_wf( - participant_id, hemisphere, + software, mem_gb, omp_nthreads, name="warp_one_hemisphere_wf", @@ -1115,8 +1109,8 @@ def init_warp_one_hemisphere_wf( from xcp_d.workflows.anatomical import init_warp_one_hemisphere_wf wf = init_warp_one_hemisphere_wf( - participant_id="01", hemisphere="L", + software="FreeSurfer", mem_gb=0.1, omp_nthreads=1, name="warp_one_hemisphere_wf", @@ -1125,6 +1119,8 @@ def init_warp_one_hemisphere_wf( Parameters ---------- hemisphere : {"L", "R"} + software : {"MCRIBS", "FreeSurfer"} + The software used for the segmentation. %(mem_gb)s %(omp_nthreads)s %(name)s @@ -1137,12 +1133,7 @@ def init_warp_one_hemisphere_wf( world_xfm merged_warpfield merged_inv_warpfield - freesurfer_path - Path to FreeSurfer derivatives. Used to load the subject's sphere file. - segmentation_software : {"FreeSurfer", "MCRIBS"} - The software used for the segmentation. - participant_id - Set from parameters. + subject_sphere Outputs ------- @@ -1157,28 +1148,18 @@ def init_warp_one_hemisphere_wf( "world_xfm", "merged_warpfield", "merged_inv_warpfield", - "freesurfer_path", - "segmentation_software", - "participant_id", + "subject_sphere", ], ), name="inputnode", ) - inputnode.inputs.participant_id = participant_id collect_registration_files = pe.Node( - CollectRegistrationFiles(hemisphere=hemisphere), + CollectRegistrationFiles(hemisphere=hemisphere, software=software), name="collect_registration_files", mem_gb=0.1, n_procs=1, ) - workflow.connect([ - (inputnode, collect_registration_files, [ - ("freesurfer_path", "segmentation_dir"), - ("participant_id", "participant_id"), - ("segmentation_software", "software"), - ]), - ]) # fmt:skip # NOTE: What does this step do? sphere_to_surf_gii = pe.Node( @@ -1187,9 +1168,7 @@ def init_warp_one_hemisphere_wf( mem_gb=mem_gb, n_procs=omp_nthreads, ) - workflow.connect([ - (collect_registration_files, sphere_to_surf_gii, [("subject_sphere", "in_file")]), - ]) # fmt:skip + workflow.connect([(inputnode, sphere_to_surf_gii, [("subject_sphere", "in_file")])]) # NOTE: What does this step do? surface_sphere_project_unproject = pe.Node( diff --git a/xcp_d/workflows/base.py b/xcp_d/workflows/base.py index 0933f8d43..bdd2d9f9c 100644 --- a/xcp_d/workflows/base.py +++ b/xcp_d/workflows/base.py @@ -135,7 +135,7 @@ def init_single_subject_wf(subject_id: str): t2w_available = subj_data["t2w"] is not None anat_mod = "t1w" if t1w_available else "t2w" - mesh_available, standard_space_mesh, mesh_files = collect_mesh_data( + mesh_available, standard_space_mesh, software, mesh_files = collect_mesh_data( layout=config.execution.layout, participant_label=subject_id, ) @@ -167,6 +167,8 @@ def init_single_subject_wf(subject_id: str): "rh_pial_surf", "lh_wm_surf", "rh_wm_surf", + "lh_subject_sphere", + "rh_subject_sphere", # morphometry files "sulcal_depth", "sulcal_curv", @@ -191,6 +193,8 @@ def init_single_subject_wf(subject_id: str): inputnode.inputs.rh_pial_surf = mesh_files["rh_pial_surf"] inputnode.inputs.lh_wm_surf = mesh_files["lh_wm_surf"] inputnode.inputs.rh_wm_surf = mesh_files["rh_wm_surf"] + inputnode.inputs.lh_subject_sphere = mesh_files["lh_subject_sphere"] + inputnode.inputs.rh_subject_sphere = mesh_files["rh_subject_sphere"] # optional surface shape files (used by surface-warping workflow) inputnode.inputs.sulcal_depth = morphometry_files["sulcal_depth"] @@ -305,12 +309,12 @@ def init_single_subject_wf(subject_id: str): # Run surface post-processing workflow if we want to warp meshes to standard space *or* # generate brainsprite. postprocess_surfaces_wf = init_postprocess_surfaces_wf( - subject_id=subject_id, mesh_available=mesh_available, standard_space_mesh=standard_space_mesh, morphometry_files=morph_file_types, t1w_available=t1w_available, t2w_available=t2w_available, + software=software, ) workflow.connect([ @@ -319,6 +323,8 @@ def init_single_subject_wf(subject_id: str): ("rh_pial_surf", "inputnode.rh_pial_surf"), ("lh_wm_surf", "inputnode.lh_wm_surf"), ("rh_wm_surf", "inputnode.rh_wm_surf"), + ("lh_subject_sphere", "inputnode.lh_subject_sphere"), + ("rh_subject_sphere", "inputnode.rh_subject_sphere"), ("anat_to_template_xfm", "inputnode.anat_to_template_xfm"), ("template_to_anat_xfm", "inputnode.template_to_anat_xfm"), ]),