From 47356717c5155c9dd99efae69eebffe1901eb283 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Mon, 29 Jul 2024 09:08:29 -0400 Subject: [PATCH] Continue to improve documentation of anatomical workflow (#1216) --- xcp_d/interfaces/bids.py | 3 +- xcp_d/interfaces/workbench.py | 86 ++++++++++++++++++----- xcp_d/workflows/anatomical.py | 128 +++++++++++++++++++++++----------- 3 files changed, 158 insertions(+), 59 deletions(-) diff --git a/xcp_d/interfaces/bids.py b/xcp_d/interfaces/bids.py index 433e5d34e..5fd4b8fe2 100644 --- a/xcp_d/interfaces/bids.py +++ b/xcp_d/interfaces/bids.py @@ -71,7 +71,7 @@ class _CollectRegistrationFilesOutputSpec(TraitedSpec): ) target_sphere = File( exists=True, - desc="Target-space sphere (fsLR for FreeSurfer, dHCP-in-fsLR for MCRIBS).", + desc="Target-space sphere (fsLR for FreeSurfer, dhcpAsym for MCRIBS).", ) sphere_to_sphere = File( exists=True, @@ -108,6 +108,7 @@ def _run_interface(self, runtime): # TODO: Collect from templateflow once it's uploaded. # FreeSurfer: fs_?/fs_?-to-fs_LR_fsaverage.?_LR.spherical_std.164k_fs_?.surf.gii + # Should be tpl-fsLR_hemi-?_space-fsaverage_den-164k_sphere.surf.gii on TemplateFlow self._results["sphere_to_sphere"] = str( load_data( f"standard_mesh_atlases/fs_{hemisphere}/" diff --git a/xcp_d/interfaces/workbench.py b/xcp_d/interfaces/workbench.py index df879037e..b31151073 100644 --- a/xcp_d/interfaces/workbench.py +++ b/xcp_d/interfaces/workbench.py @@ -294,6 +294,40 @@ class _SurfaceSphereProjectUnprojectOutputSpec(TraitedSpec): class SurfaceSphereProjectUnproject(WBCommand): """Copy registration deformations to different sphere. + A surface registration starts with an input sphere, + and moves its vertices around on the sphere until it matches the template data. + This means that the registration deformation is actually represented as the difference + between two separate files - the starting sphere, and the registered sphere. + Since the starting sphere of the registration may not have vertex correspondence to any + other sphere (often, it is a native sphere), + it can be inconvenient to manipulate or compare these deformations across subjects, etc. + + The purpose of this command is to be able to apply these deformations onto a new sphere + of the user's choice, + to make it easier to compare or manipulate them. + Common uses are to concatenate two successive separate registrations + (e.g. Human to Chimpanzee, and then Chimpanzee to Macaque) + or inversion (for dedrifting or symmetric registration schemes). + + must already be considered to be in alignment with one of the two ends of the + registration (if your registration is Human to Chimpanzee, + must be in register with either Human or Chimpanzee). + + The 'project-to' sphere must be the side of the registration that is aligned with + (if your registration is Human to Chimpanzee, and is aligned with Human, then + 'project-to' should be the original Human sphere). + + The 'unproject-from' sphere must be the remaining sphere of the registration + (original vs deformed/registered). + The output is as if you had run the same registration with as the starting sphere, + in the direction of deforming the 'project-to' sphere to create the 'unproject-from' sphere. + + Note that this command cannot check for you what spheres are aligned with other spheres, + and using the wrong spheres or in the incorrect order will not necessarily cause an error + message. + In some cases, it may be useful to use a new, arbitrary sphere as the input, + which can be created with the -surface-create-sphere command. + .. code-block:: wb_command -surface-sphere-project-unproject @@ -615,37 +649,55 @@ class CiftiParcellateWorkbench(WBCommand): class _CiftiSurfaceResampleInputSpec(CommandLineInputSpec): - """Input specification for the CiftiSurfaceResample command.""" + """Input specification for the CiftiSurfaceResample command. + + Resamples a surface file, given two spherical surfaces that are in register. + If ADAP_BARY_AREA is used, exactly one of -area-surfs or -area-metrics must be specified. + This method is not generally recommended for surface resampling, + but is provided for completeness. + + For cut surfaces (including flatmaps), use -surface-cut-resample. + + Instead of resampling a spherical surface, + the -surface-sphere-project-unproject command is recommended. + """ in_file = File( exists=True, mandatory=True, - argstr="%s ", + argstr="%s", position=0, - desc="the gifti file", + desc="the surface file to resample", ) current_sphere = File( exists=True, position=1, - argstr=" %s", - desc=" the current sphere surface in gifti for in_file", + argstr="%s", + desc="a sphere surface with the mesh that the input surface is currently on", ) new_sphere = File( exists=True, position=2, - argstr=" %s", - desc=" the new sphere surface to be resample the in_file to, eg fsaverag5 or fsl32k", - ) - - metric = traits.Str( - argstr=" %s ", + argstr="%s", + desc=( + "a sphere surface that is in register with " + "and has the desired output mesh" + ), + ) + method = traits.Enum( + "BARYCENTRIC", + "ADAP_BARY_AREA", + argstr="%s", position=3, - desc=" fixed for anatomic", - default=" BARYCENTRIC ", + desc=( + "the method name. " + "The BARYCENTRIC method is generally recommended for anatomical surfaces, " + "in order to minimize smoothing." + ), + default="BARYCENTRIC", ) - out_file = File( name_source=["in_file"], name_template="resampled_%s.gii", @@ -653,7 +705,7 @@ class _CiftiSurfaceResampleInputSpec(CommandLineInputSpec): extensions=[".shape.gii", ".surf.gii"], argstr="%s", position=4, - desc="The gifti output, either left and right", + desc="The output surface file.", ) @@ -664,14 +716,14 @@ class _CiftiSurfaceResampleOutputSpec(TraitedSpec): class CiftiSurfaceResample(WBCommand): - """Resample a surface from one sphere to another. + """Resample a surface to a different mesh. TODO: Improve documentation. """ input_spec = _CiftiSurfaceResampleInputSpec output_spec = _CiftiSurfaceResampleOutputSpec - _cmd = "wb_command -surface-resample" + _cmd = "wb_command -surface-resample" class _CiftiSeparateMetricInputSpec(CommandLineInputSpec): diff --git a/xcp_d/workflows/anatomical.py b/xcp_d/workflows/anatomical.py index 0f0cc5275..04e081295 100644 --- a/xcp_d/workflows/anatomical.py +++ b/xcp_d/workflows/anatomical.py @@ -870,6 +870,8 @@ def init_generate_hcp_surfaces_wf(name="generate_hcp_surfaces_wf"): def init_ants_xfm_to_fsl_wf(mem_gb, omp_nthreads, name="ants_xfm_to_fsl_wf"): """Modify ANTS-style fMRIPrep transforms to work with Connectome Workbench/FSL FNIRT. + XXX: Does this only work if the template is MNI152NLin6Asym? + Workflow Graph .. workflow:: :graph2use: orig @@ -893,18 +895,21 @@ def init_ants_xfm_to_fsl_wf(mem_gb, omp_nthreads, name="ants_xfm_to_fsl_wf"): Inputs ------ anat_to_template_xfm - ANTS/fMRIPrep-style H5 transform from T1w image to template. + ANTS/fMRIPrep-style H5 transform from anatomical image to template. template_to_anat_xfm - ANTS/fMRIPrep-style H5 transform from template to T1w image. + ANTS/fMRIPrep-style H5 transform from template to anatomical image. Outputs ------- world_xfm - TODO: Add description. + The affine portion of the volumetric anatomical-to-template transform, + in NIfTI (world) format. merged_warpfield - TODO: Add description. + The warpfield portion of the volumetric anatomical-to-template transform, + in FSL (FNIRT) format. merged_inv_warpfield - TODO: Add description. + The warpfield portion of the volumetric template-to-anatomical transform, + in FSL (FNIRT) format. """ workflow = Workflow(name=name) @@ -919,7 +924,7 @@ def init_ants_xfm_to_fsl_wf(mem_gb, omp_nthreads, name="ants_xfm_to_fsl_wf"): ) # Now we can start the actual workflow. - # use ANTs CompositeTransformUtil to separate the .h5 into affine and warpfield xfms + # Use ANTs CompositeTransformUtil to separate the .h5 into affine and warpfield xfms. disassemble_h5 = pe.Node( CompositeTransformUtil( process="disassemble", @@ -928,7 +933,7 @@ def init_ants_xfm_to_fsl_wf(mem_gb, omp_nthreads, name="ants_xfm_to_fsl_wf"): name="disassemble_h5", mem_gb=mem_gb, n_procs=omp_nthreads, - ) # MB + ) workflow.connect([(inputnode, disassemble_h5, [("anat_to_template_xfm", "in_file")])]) # Nipype's CompositeTransformUtil assumes a certain file naming and @@ -945,32 +950,27 @@ def init_ants_xfm_to_fsl_wf(mem_gb, omp_nthreads, name="ants_xfm_to_fsl_wf"): ) workflow.connect([(inputnode, disassemble_h5_inv, [("template_to_anat_xfm", "in_file")])]) - # convert affine from ITK binary to txt - convert_ants_transform = pe.Node( + # Convert anat-to-template affine from ITK binary to txt + convert_ants_xfm = pe.Node( ConvertTransformFile(dimension=3), - name="convert_ants_transform", + name="convert_ants_xfm", ) - workflow.connect([ - (disassemble_h5, convert_ants_transform, [("affine_transform", "in_transform")]), - ]) # fmt:skip + workflow.connect([(disassemble_h5, convert_ants_xfm, [("affine_transform", "in_transform")])]) - # change xfm type from "AffineTransform" to "MatrixOffsetTransformBase" + # Change xfm type from "AffineTransform" to "MatrixOffsetTransformBase" # since wb_command doesn't recognize "AffineTransform" - # (AffineTransform is a subclass of MatrixOffsetTransformBase - # which makes this okay to do AFAIK) + # (AffineTransform is a subclass of MatrixOffsetTransformBase which prob makes this okay to do) change_xfm_type = pe.Node(ChangeXfmType(), name="change_xfm_type") - workflow.connect([ - (convert_ants_transform, change_xfm_type, [("out_transform", "in_transform")]), - ]) # fmt:skip + workflow.connect([(convert_ants_xfm, change_xfm_type, [("out_transform", "in_transform")])]) - # convert affine xfm to "world" so it works with -surface-apply-affine + # Convert affine xfm to "world" so it works with -surface-apply-affine convert_xfm2world = pe.Node( ConvertAffine(fromwhat="itk", towhat="world"), name="convert_xfm2world", ) workflow.connect([(change_xfm_type, convert_xfm2world, [("out_transform", "in_file")])]) - # use C3d to separate the combined warpfield xfm into x, y, and z components + # Use C3d to separate the combined warpfield xfm into x, y, and z components get_xyz_components = pe.Node( C3d( is_4d=True, @@ -996,7 +996,7 @@ def init_ants_xfm_to_fsl_wf(mem_gb, omp_nthreads, name="ants_xfm_to_fsl_wf"): (disassemble_h5_inv, get_inv_xyz_components, [("displacement_field", "in_file")]), ]) # fmt:skip - # select x-component after separating warpfield above + # Select x-component after separating warpfield above select_x_component = pe.Node( niu.Select(index=[0]), name="select_x_component", @@ -1010,7 +1010,7 @@ def init_ants_xfm_to_fsl_wf(mem_gb, omp_nthreads, name="ants_xfm_to_fsl_wf"): n_procs=omp_nthreads, ) - # select y-component + # Select y-component select_y_component = pe.Node( niu.Select(index=[1]), name="select_y_component", @@ -1024,7 +1024,7 @@ def init_ants_xfm_to_fsl_wf(mem_gb, omp_nthreads, name="ants_xfm_to_fsl_wf"): n_procs=omp_nthreads, ) - # select z-component + # Select z-component select_z_component = pe.Node( niu.Select(index=[2]), name="select_z_component", @@ -1046,7 +1046,7 @@ def init_ants_xfm_to_fsl_wf(mem_gb, omp_nthreads, name="ants_xfm_to_fsl_wf"): (get_inv_xyz_components, select_inv_z_component, [("out_files", "inlist")]), ]) # fmt:skip - # reverse y-component of the warpfield + # Reverse y-component of the warpfield # (need to do this when converting a warpfield from ANTs to FNIRT format # for use with wb_command -surface-apply-warpfield) reverse_y_component = pe.Node( @@ -1122,6 +1122,12 @@ def init_warp_one_hemisphere_wf( ): """Apply transforms to warp one hemisphere's surface files into standard space. + Basically, the resulting surface files will have the same vertices as the standard-space + surfaces, but the coordinates/mesh of those vertices will be the subject's native-space + coordinates/mesh. + This way we can visualize surface statistical maps on the subject's unique morphology + (sulci, gyri, etc.). + Workflow Graph .. workflow:: :graph2use: orig @@ -1150,15 +1156,54 @@ def init_warp_one_hemisphere_wf( Inputs ------ hemi_files : list of str - A list of surface files for the requested hemisphere, in fsnative space. + A list of surface files (i.e., pial and white matter) for the requested hemisphere, + in fsnative space. world_xfm + The affine portion of the volumetric anatomical-to-template transform, + in NIfTI (world) format. merged_warpfield + The warpfield portion of the volumetric anatomical-to-template transform, + in FSL (FNIRT) format. merged_inv_warpfield + The warpfield portion of the volumetric template-to-anatomical transform, + in FSL (FNIRT) format. subject_sphere + The subject's fsnative sphere registration file to fsaverage + (sphere.reg in FreeSurfer parlance). + The file contains the vertices from the subject's fsnative sphere, + with coordinates that are aligned to the fsaverage sphere. Outputs ------- - warped_hemi_files + warped_hemi_files : list of str + The ``hemi_files`` warped from fsnative space to standard space. + + Notes + ----- + Steps: + + 1. Collect the registration files needed for the warp. + 2. Convert the subject's sphere to a GIFTI file. + - This step is unnecessary since fMRIPrep and Nibabies already write out a GIFTI file. + 3. Project the subject's fsnative-in-fsaverage sphere to a high-resolution + target-sphere-in-fsaverage-space. + This retains the subject's fsnative sphere's resolution and vertices + (e.g., 120079 vertices), but the coordinates are now aligned to the target sphere's space. + - For Freesurfer, this is the fsLR-164k-in-fsaverage sphere. + - For MCRIBS, this is the dhcpAsym-41k-in-fsaverage sphere. + - Nibabies and fMRIPrep do this already to produce the + space-_desc-reg_sphere files, so XCP-D could directly use those and skip + this step. + 4. Apply the warped sphere from the previous step to warp the pial and white matter surfaces + to the target space. This includes downsampling to 32k. + - For Freesurfer, this means the coordinates for these files are fsLR-32k. + - For MCRIBS, this means the coordinates for these files are dhcpAsym-32k. + 5. Apply the anatomical-to-template affine transform to the 32k surfaces. + 6. Apply the anatomical-to-template warpfield to the 32k surfaces. + This and the previous step make it so you can overlay the pial and white matter surfaces + on the associated volumetric template (e.g., for XCP-D's brainsprite). + - This important thing is that the volumetric template must match the template space + used here. """ workflow = Workflow(name=name) @@ -1174,6 +1219,10 @@ def init_warp_one_hemisphere_wf( ), name="inputnode", ) + outputnode = pe.Node( + niu.IdentityInterface(fields=["warped_hemi_files"]), + name="outputnode", + ) collect_registration_files = pe.Node( CollectRegistrationFiles(hemisphere=hemisphere, software=software), @@ -1182,7 +1231,8 @@ def init_warp_one_hemisphere_wf( n_procs=1, ) - # NOTE: What does this step do? + # XXX: Given that fMRIPrep and Nibabies write out the subject spheres as surf.gii, + # I think this is unnecessary. sphere_to_surf_gii = pe.Node( MRIsConvert(out_datatype="gii"), name="sphere_to_surf_gii", @@ -1192,6 +1242,10 @@ def init_warp_one_hemisphere_wf( workflow.connect([(inputnode, sphere_to_surf_gii, [("subject_sphere", "in_file")])]) # NOTE: What does this step do? + # Project the subject's sphere (fsnative) to the source-sphere (fsaverage) using the + # fsLR/dhcpAsym-in-fsaverage + # (fsLR or dhcpAsym vertices with coordinates on the fsaverage sphere) sphere? + # So what's the result? The fsLR or dhcpAsym vertices with coordinates on the fsnative sphere? surface_sphere_project_unproject = pe.Node( SurfaceSphereProjectUnproject(), name="surface_sphere_project_unproject", @@ -1204,10 +1258,9 @@ def init_warp_one_hemisphere_wf( (sphere_to_surf_gii, surface_sphere_project_unproject, [("converted", "in_file")]), ]) # fmt:skip - # resample the surfaces to fsLR-32k - # NOTE: Does that mean the data are in fsLR-164k before this? + # Resample the pial and white matter surfaces from fsnative to fsLR-32k or dhcpAsym-32k resample_to_fsLR32k = pe.MapNode( - CiftiSurfaceResample(metric="BARYCENTRIC"), + CiftiSurfaceResample(method="BARYCENTRIC"), name="resample_to_fsLR32k", mem_gb=mem_gb, n_procs=omp_nthreads, @@ -1219,8 +1272,8 @@ def init_warp_one_hemisphere_wf( (surface_sphere_project_unproject, resample_to_fsLR32k, [("out_file", "current_sphere")]), ]) # fmt:skip - # apply affine to 32k surfs - # NOTE: What does this step do? Aren't the data in fsLR-32k from resample_to_fsLR32k? + # Apply FLIRT-format anatomical-to-template affine transform to 32k surfs + # NOTE: What does this step do? Aren't the data in fsLR/dhcpAsym-32k from resample_to_fsLR32k? apply_affine_to_fsLR32k = pe.MapNode( ApplyAffine(), name="apply_affine_to_fsLR32k", @@ -1233,7 +1286,7 @@ def init_warp_one_hemisphere_wf( (resample_to_fsLR32k, apply_affine_to_fsLR32k, [("out_file", "in_file")]), ]) # fmt:skip - # apply FNIRT-format warpfield + # Apply FNIRT-format (forward) anatomical-to-template warpfield # NOTE: What does this step do? apply_warpfield_to_fsLR32k = pe.MapNode( ApplyWarpfield(), @@ -1248,13 +1301,6 @@ def init_warp_one_hemisphere_wf( ("merged_inv_warpfield", "warpfield"), ]), (apply_affine_to_fsLR32k, apply_warpfield_to_fsLR32k, [("out_file", "in_file")]), - ]) # fmt:skip - - outputnode = pe.Node( - niu.IdentityInterface(fields=["warped_hemi_files"]), - name="outputnode", - ) - workflow.connect([ (apply_warpfield_to_fsLR32k, outputnode, [("out_file", "warped_hemi_files")]), ]) # fmt:skip