diff --git a/.circleci/ds005_run.sh b/.circleci/ds005_run.sh index 864ea87f52..6981fe7478 100644 --- a/.circleci/ds005_run.sh +++ b/.circleci/ds005_run.sh @@ -17,6 +17,7 @@ docker run -it -e FMRIPREP_DEV=1 -u $(id -u) \ --skull-strip-template MNI152NLin2009cAsym:res-2 \ --sloppy --mem-gb 4 \ --ncpus 2 --omp-nthreads 2 -vv \ + --no-msm \ --fs-license-file /tmp/fslicense/license.txt \ --fs-subjects-dir /tmp/ds005/freesurfer \ ${@:1} diff --git a/.circleci/ds054_run.sh b/.circleci/ds054_run.sh index b984a44b21..366479cc92 100644 --- a/.circleci/ds054_run.sh +++ b/.circleci/ds054_run.sh @@ -17,5 +17,6 @@ docker run --rm -it -e FMRIPREP_DEV=1 -u $(id -u) \ --skull-strip-template OASIS30ANTs:res-1 \ --output-spaces MNI152Lin MNI152NLin2009cAsym:res-2:res-native \ --mem-gb 4 --ncpus 2 --omp-nthreads 2 -vv \ + --no-msm \ --fs-license-file /tmp/fslicense/license.txt \ ${@:1} diff --git a/Dockerfile b/Dockerfile index a74b8c39b8..1c7f53c3ad 100644 --- a/Dockerfile +++ b/Dockerfile @@ -241,6 +241,10 @@ ENV LANG="C.UTF-8" \ FSLREMOTECALL="" \ FSLGECUDAQ="cuda.q" +# MSM +RUN curl -L -H "Accept: application/octet-stream" https://api.github.com/repos/ecr05/MSM_HOCR/releases/assets/16253707 -o /usr/local/bin/msm \ + && chmod +x /usr/local/bin/msm + # Unless otherwise specified each process should only use one thread - nipype # will handle parallelization ENV MKL_NUM_THREADS=1 \ diff --git a/scripts/fetch_templates.py b/scripts/fetch_templates.py index 65ed6c8a82..5987cced57 100755 --- a/scripts/fetch_templates.py +++ b/scripts/fetch_templates.py @@ -79,10 +79,14 @@ def fetch_fsaverage(): tpl-fsaverage/tpl-fsaverage_hemi-R_den-164k_desc-std_sphere.surf.gii tpl-fsaverage/tpl-fsaverage_hemi-L_den-164k_desc-vaavg_midthickness.shape.gii tpl-fsaverage/tpl-fsaverage_hemi-R_den-164k_desc-vaavg_midthickness.shape.gii + tpl-fsaverage/tpl-fsaverage_hemi-L_den-164k_sulc.shape.gii + tpl-fsaverage/tpl-fsaverage_hemi-R_den-164k_sulc.shape.gii """ template = "fsaverage" tf.get(template, density="164k", suffix="dseg", extension=".tsv") + tf.get(template, density='164k', desc='std', suffix='sphere', extension='.surf.gii') + tf.get(template, density='164k', suffix='sulc', extension='.shape.gii') def fetch_all(): diff --git a/smriprep/__about__.py b/smriprep/__about__.py index 12e0aeb140..41d821a637 100644 --- a/smriprep/__about__.py +++ b/smriprep/__about__.py @@ -24,7 +24,10 @@ Base module variables """ -from ._version import __version__ +try: + from ._version import __version__ +except ImportError: # pragma: no cover + __version__ = "0+unknown" __copyright__ = "Copyright 2019, Center for Reproducible Neuroscience, Stanford University" __credits__ = [ diff --git a/smriprep/cli/run.py b/smriprep/cli/run.py index f4cf25be6b..ffff602d21 100644 --- a/smriprep/cli/run.py +++ b/smriprep/cli/run.py @@ -214,6 +214,12 @@ def get_parser(): dest="hires", help="disable sub-millimeter (hires) reconstruction", ) + g_surfs.add_argument( + "--no-msm", + action="store_false", + dest="msm_sulc", + help="Disable Multimodal Surface Matching surface registration." + ) g_surfs_xor = g_surfs.add_mutually_exclusive_group() g_surfs_xor.add_argument( @@ -567,6 +573,7 @@ def build_workflow(opts, retval): layout=layout, longitudinal=opts.longitudinal, low_mem=opts.low_mem, + msm_sulc=opts.msm_sulc, omp_nthreads=omp_nthreads, output_dir=str(output_dir), run_uuid=run_uuid, diff --git a/smriprep/conftest.py b/smriprep/conftest.py index 138e5b9db7..839237af98 100644 --- a/smriprep/conftest.py +++ b/smriprep/conftest.py @@ -1,3 +1,16 @@ import os +import pytest +import numpy + +from smriprep.data import load_resource + os.environ['NO_ET'] = '1' + + +@pytest.fixture(autouse=True) +def populate_namespace(doctest_namespace, tmp_path): + doctest_namespace['os'] = os + doctest_namespace['np'] = numpy + doctest_namespace['load'] = load_resource + doctest_namespace['testdir'] = tmp_path diff --git a/smriprep/data/boilerplate.bib b/smriprep/data/boilerplate.bib index 9acffc48dc..93a5db1b7f 100644 --- a/smriprep/data/boilerplate.bib +++ b/smriprep/data/boilerplate.bib @@ -322,3 +322,15 @@ @article{posse_t2s volume = 42, year = 1999 } + +@article{msm, + author = {Emma C. Robinson and Saad Jbabdi and Matthew F. Glasser and Jesper Andersson and Gregory C. Burgess and Michael P. Harms and Stephen M. Smith and David C. Van Essen and Mark Jenkinson}, + doi = {10.1016/j.neuroimage.2014.05.069}, + year = 2014, + month = {oct}, + volume = {100}, + pages = {414-426}, + title = {{MSM}: A new flexible framework for Multimodal Surface Matching}, + url = {https://doi.org/10.1016/j.neuroimage.2014.05.069}, + journal = {NeuroImage} +} diff --git a/smriprep/data/msm/MSMSulcStrainFinalconf b/smriprep/data/msm/MSMSulcStrainFinalconf new file mode 100644 index 0000000000..f1b2f5b806 --- /dev/null +++ b/smriprep/data/msm/MSMSulcStrainFinalconf @@ -0,0 +1,18 @@ +--simval=3,2,2,2 +--sigma_in=0,0,0,0 +--sigma_ref=0,0,0,0 +--lambda=0,10,7.5,7.5 +--it=50,10,15,15 +--opt=AFFINE,DISCRETE,DISCRETE,DISCRETE +--CPgrid=6,2,3,4 +--SGgrid=6,4,5,6 +--datagrid=6,4,5,6 +--regoption=3 +--regexp=2 +--dopt=HOCR +--VN +--rescaleL +--triclique +--k_exponent=2 +--bulkmod=1.6 +--shearmod=0.4 diff --git a/smriprep/interfaces/msm.py b/smriprep/interfaces/msm.py new file mode 100644 index 0000000000..b9148a7f0c --- /dev/null +++ b/smriprep/interfaces/msm.py @@ -0,0 +1,159 @@ +from pathlib import Path + +from nipype.interfaces.base import ( + CommandLine, + CommandLineInputSpec, + File, + TraitedSpec, + traits, +) + + +class MSMInputSpec(CommandLineInputSpec): + in_mesh = File( + exists=True, + mandatory=True, + argstr="--inmesh=%s", + desc="input mesh (available formats: VTK, ASCII, GIFTI). Needs to be a sphere", + ) + out_base = File( + name_source=["in_mesh"], + name_template="%s_msm", + argstr="--out=%s", + desc="output basename", + ) + reference_mesh = File( + exists=True, + argstr="--refmesh=%s", + desc="reference mesh (available formats: VTK, ASCII, GIFTI). Needs to be a sphere." + "If not included algorithm assumes reference mesh is equivalent input", + ) + in_data = File( + exists=True, + argstr="--indata=%s", + desc="scalar or multivariate data for input - can be ASCII (.asc,.dpv,.txt) " + "or GIFTI (.func.gii or .shape.gii)", + ) + reference_data = File( + exists=True, + argstr="--refdata=%s", + desc="scalar or multivariate data for reference - can be ASCII (.asc,.dpv,.txt) " + "or GIFTI (.func.gii or .shape.gii)", + ) + transformed_mesh = File( + exists=True, + argstr="--trans=%s", + desc="Transformed source mesh (output of a previous registration). " + "Use this to initiliase the current registration.", + ) + in_register = File( + exists=True, + argstr="--in_register=%s", + desc="Input mesh at data resolution. Used to resample data onto input mesh if data " + "is supplied at a different resolution. Note this mesh HAS to be in alignment with " + "either the input_mesh of (if supplied) the transformed source mesh. " + "Use with supreme caution.", + ) + in_weight = File( + exists=True, + argstr="--inweight=%s", + desc="cost function weighting for input - weights data in these vertices when calculating " + "similarity (ASCII or GIFTI). Can be multivariate provided dimension equals that of data", + ) + reference_weight = File( + exists=True, + argstr="--refweight=%s", + desc="cost function weighting for reference - weights data in these vertices when " + "calculating similarity (ASCII or GIFTI). Can be multivariate provided dimension " + "equals that of data", + ) + output_format = traits.Enum( + "GIFTI", + "VTK", + "ASCII", + "ASCII_MAT", + argstr="--format=%s", + desc="format of output files", + ) + config_file = File( + exists=True, + argstr="--conf=%s", + desc="configuration file", + ) + levels = traits.Int( + argstr="--levels=%d", + desc="number of resolution levels (default = number of resolution levels specified " + "by --opt in config file)", + ) + smooth_output_sigma = traits.Int( + argstr="--smoothout=%d", + desc="smooth tranformed output with this sigma (default=0)", + ) + verbose = traits.Bool( + argstr="--verbose", + desc="switch on diagnostic messages", + ) + + +class MSMOutputSpec(TraitedSpec): + warped_mesh = File( + desc="the warped input mesh (i.e., new vertex locations - this capture the warp field, " + "much like a _warp.nii.gz file would for volumetric warps created by FNIRT)." + ) + downsampled_warped_mesh = File( + desc="a downsampled version of the warped_mesh where the resolution of this mesh will " + "be equivalent to the resolution of the final datamesh" + ) + warped_data = File( + desc="the input data passed through the MSM warp and projected onto the target surface" + ) + + +class MSM(CommandLine): + """ + MSM (Multimodal Surface Matching) is a tool for registering cortical surfaces. + The tool has been developed and tested using FreeSurfer extracted surfaces. + However, in principle the tool with work with any cortical surface extraction method provided + the surfaces can be mapped to the sphere. + The key advantage of the method is that alignment may be driven using a wide variety of + univariate (sulcal depth, curvature, myelin), multivariate (Task fMRI, or Resting State + Networks) or multimodal (combinations of folding, myelin and fMRI) feature sets. + + The main MSM tool is currently run from the command line using the program ``msm``. + This enables fast alignment of spherical cortical surfaces by utilising a fast discrete + optimisation framework (FastPD Komodakis 2007), which significantly reduces the search + space of possible deformations for each vertex, and allows flexibility with regards to the + choice of similarity metric used to match the images. + + >>> msm = MSM( + ... config_file=load('msm/MSMSulcStrainFinalconf'), + ... in_mesh='sub-01_hemi-L_sphere.surf.gii', + ... reference_mesh='tpl-fsaverage_hemi-L_den-164k_desc-std_sphere.surf.gii', + ... in_data='sub-01_hemi-L_sulc.shape.gii', + ... reference_data='tpl-fsaverage_hemi-L_den-164k_sulc.shape.gii', + ... out_base='L.', + ... ) + >>> msm.cmdline # doctest: +ELLIPSIS +NORMALIZE_WHITESPACE + 'msm --conf=.../MSMSulcStrainFinalconf \ + --indata=sub-01_hemi-L_sulc.shape.gii \ + --inmesh=sub-01_hemi-L_sphere.surf.gii \ + --out=L. \ + --refdata=tpl-fsaverage_hemi-L_den-164k_sulc.shape.gii \ + --refmesh=tpl-fsaverage_hemi-L_den-164k_desc-std_sphere.surf.gii' + + """ + + input_spec = MSMInputSpec + output_spec = MSMOutputSpec + _cmd = "msm" + + def _list_outputs(self): + from nipype.utils.filemanip import split_filename + + outputs = self._outputs().get() + out_base = self.inputs.out_base or split_filename(self.inputs.in_mesh)[1] + cwd = Path.cwd() + outputs['warped_mesh'] = str(cwd / (out_base + 'sphere.reg.surf.gii')) + outputs['downsampled_warped_mesh'] = str(cwd / (out_base + 'sphere.LR.reg.surf.gii')) + outputs['warped_data'] = str(cwd / (out_base + 'transformed_and_reprojected.func.gii')) + return outputs diff --git a/smriprep/interfaces/tests/data/sub-01_desc-warped_T1w.nii.gz b/smriprep/interfaces/tests/data/sub-01_desc-warped_T1w.nii.gz new file mode 100644 index 0000000000..e69de29bb2 diff --git a/smriprep/interfaces/tests/data/sub-01_hemi-L_sulc.shape.gii b/smriprep/interfaces/tests/data/sub-01_hemi-L_sulc.shape.gii new file mode 100644 index 0000000000..e69de29bb2 diff --git a/smriprep/interfaces/tests/data/tpl-fsaverage_hemi-L_den-164k_desc-std_sphere.surf.gii b/smriprep/interfaces/tests/data/tpl-fsaverage_hemi-L_den-164k_desc-std_sphere.surf.gii new file mode 100644 index 0000000000..e69de29bb2 diff --git a/smriprep/interfaces/tests/data/tpl-fsaverage_hemi-L_den-164k_sulc.shape.gii b/smriprep/interfaces/tests/data/tpl-fsaverage_hemi-L_den-164k_sulc.shape.gii new file mode 100644 index 0000000000..e69de29bb2 diff --git a/smriprep/interfaces/workbench.py b/smriprep/interfaces/workbench.py index 9e4c59caa0..1916a52396 100644 --- a/smriprep/interfaces/workbench.py +++ b/smriprep/interfaces/workbench.py @@ -144,6 +144,214 @@ class CreateSignedDistanceVolume(WBCommand): _cmd = "wb_command -create-signed-distance-volume" +class SurfaceAffineRegressionInputSpec(CommandLineInputSpec): + in_surface = File( + exists=True, + mandatory=True, + argstr="%s", + position=0, + desc="Surface to warp", + ) + target_surface = File( + exists=True, + mandatory=True, + argstr="%s", + position=1, + desc="Surface to match the coordinates of", + ) + out_affine = File( + name_template="%s_xfm", + name_source=["in_surface"], + argstr="%s", + position=2, + desc="the output affine file", + ) + + +class SurfaceAffineRegressionOutputSpec(TraitedSpec): + out_affine = File(desc="The output affine file") + + +class SurfaceAffineRegression(WBCommand): + """REGRESS THE AFFINE TRANSFORM BETWEEN SURFACES ON THE SAME MESH + + wb_command -surface-affine-regression + - the surface to warp + - the surface to match the coordinates of + - output - the output affine file + + Use linear regression to compute an affine that minimizes the sum of + squares of the coordinate differences between the target surface and the + warped source surface. Note that this has a bias to shrink the surface + that is being warped. The output is written as a NIFTI 'world' matrix, + see -convert-affine to convert it for use in other software. + + >>> sar = SurfaceAffineRegression() + >>> sar.inputs.in_surface = 'sub-01_hemi-L_sulc.shape.gii' + >>> sar.inputs.target_surface = 'tpl-fsaverage_hemi-L_den-164k_sulc.shape.gii' + >>> sar.cmdline # doctest: +ELLIPSIS +NORMALIZE_WHITESPACE + 'wb_command -surface-affine-regression \ + sub-01_hemi-L_sulc.shape.gii \ + tpl-fsaverage_hemi-L_den-164k_sulc.shape.gii \ + sub-01_hemi-L_sulc.shape_xfm' + """ + input_spec = SurfaceAffineRegressionInputSpec + output_spec = SurfaceAffineRegressionOutputSpec + _cmd = "wb_command -surface-affine-regression" + + +class SurfaceApplyAffineInputSpec(CommandLineInputSpec): + in_surface = File( + exists=True, + mandatory=True, + argstr="%s", + position=0, + desc="the surface to transform", + ) + in_affine = File( + exists=True, + mandatory=True, + argstr="%s", + position=1, + desc="the affine file", + ) + out_surface = File( + name_template="%s_xformed.surf.gii", + name_source=["in_surface"], + argstr="%s", + position=2, + desc="the output transformed surface", + ) + flirt_source = File( + exists=True, + requires=["flirt_target"], + argstr="-flirt %s", + position=3, + desc="Source volume (must be used if affine is a flirt affine)", + ) + flirt_target = File( + exists=True, + requires=["flirt_source"], + argstr="%s", + position=4, + desc="Target volume (must be used if affine is a flirt affine)", + ) + + +class SurfaceApplyAffineOutputSpec(TraitedSpec): + out_surface = File(desc="the output transformed surface") + + +class SurfaceApplyAffine(WBCommand): + """APPLY AFFINE TRANSFORM TO SURFACE FILE + + wb_command -surface-apply-affine + - the surface to transform + - the affine file + - output - the output transformed surface + + [-flirt] - MUST be used if affine is a flirt affine + - the source volume used when generating the affine + - the target volume used when generating the affine + + For flirt matrices, you must use the -flirt option, because flirt + matrices are not a complete description of the coordinate transform they + represent. If the -flirt option is not present, the affine must be a + nifti 'world' affine, which can be obtained with the -convert-affine + command, or aff_conv from the 4dfp suite. + + .. testsetup:: + + >>> np.savetxt('affine.txt', np.eye(4), delimiter='\t') + + .. doctest:: + + >>> saa = SurfaceApplyAffine() + >>> saa.inputs.in_surface = 'sub-01_hemi-L_sphere.surf.gii' + >>> saa.inputs.in_affine = 'affine.txt' + >>> saa.cmdline # doctest: +ELLIPSIS +NORMALIZE_WHITESPACE + 'wb_command -surface-apply-affine \ + sub-01_hemi-L_sphere.surf.gii \ + affine.txt \ + sub-01_hemi-L_sphere.surf_xformed.surf.gii' + + .. testcleanup:: + + >>> os.unlink('affine.txt') + """ + input_spec = SurfaceApplyAffineInputSpec + output_spec = SurfaceApplyAffineOutputSpec + _cmd = "wb_command -surface-apply-affine" + + +class SurfaceApplyWarpfieldInputSpec(CommandLineInputSpec): + in_surface = File( + exists=True, + mandatory=True, + argstr="%s", + position=0, + desc="the surface to transform", + ) + warpfield = File( + exists=True, + mandatory=True, + argstr="%s", + position=1, + desc="the INVERSE warpfield", + ) + out_surface = File( + name_template="%s_warped.surf.gii", + name_source=["in_surface"], + argstr="%s", + position=2, + desc="the output transformed surface", + ) + fnirt_forward_warp = File( + exists=True, + argstr="-fnirt %s", + position=3, + desc="the forward warpfield (must be used if fnirt warpfield)", + ) + + +class SurfaceApplyWarpfieldOutputSpec(TraitedSpec): + out_surface = File(desc="the output transformed surface") + + +class SurfaceApplyWarpfield(WBCommand): + """APPLY WARPFIELD TO SURFACE FILE + + wb_command -surface-apply-warpfield + - the surface to transform + - the INVERSE warpfield + - output - the output transformed surface + + [-fnirt] - MUST be used if using a fnirt warpfield + - the forward warpfield + + NOTE: warping a surface requires the INVERSE of the warpfield used to + warp the volume it lines up with. The header of the forward warp is + needed by the -fnirt option in order to correctly interpret the + displacements in the fnirt warpfield. + + If the -fnirt option is not present, the warpfield must be a nifti + 'world' warpfield, which can be obtained with the -convert-warpfield + command. + + >>> saw = SurfaceApplyWarpfield() + >>> saw.inputs.in_surface = 'sub-01_hemi-L_sphere.surf.gii' + >>> saw.inputs.warpfield = 'sub-01_desc-warped_T1w.nii.gz' + >>> saw.cmdline # doctest: +ELLIPSIS +NORMALIZE_WHITESPACE + 'wb_command -surface-apply-warpfield \ + sub-01_hemi-L_sphere.surf.gii \ + sub-01_desc-warped_T1w.nii.gz \ + sub-01_hemi-L_sphere.surf_warped.surf.gii' + """ + input_spec = SurfaceApplyWarpfieldInputSpec + output_spec = SurfaceApplyWarpfieldOutputSpec + _cmd = "wb_command -surface-apply-warpfield" + + class SurfaceSphereProjectUnprojectInputSpec(TraitedSpec): """COPY REGISTRATION DEFORMATIONS TO DIFFERENT SPHERE. diff --git a/smriprep/workflows/anatomical.py b/smriprep/workflows/anatomical.py index d458d95e90..6187044895 100644 --- a/smriprep/workflows/anatomical.py +++ b/smriprep/workflows/anatomical.py @@ -77,6 +77,7 @@ def init_anat_preproc_wf( cifti_output=False, debug=False, sloppy=False, + msm_sulc=False, existing_derivatives=None, name="anat_preproc_wf", skull_strip_fixed_seed=False, @@ -529,7 +530,7 @@ def _check_img(img): name="surface_recon_wf", omp_nthreads=omp_nthreads, hires=hires ) applyrefined = pe.Node(fsl.ApplyMask(), name="applyrefined") - sphere_reg_wf = init_sphere_reg_wf(name="sphere_reg_wf") + sphere_reg_wf = init_sphere_reg_wf(msm_sulc=msm_sulc, name="sphere_reg_wf") if t2w: t2w_template_wf = init_anat_template_wf( @@ -619,6 +620,7 @@ def _check_img(img): (surface_recon_wf, sphere_reg_wf, [ ('outputnode.subject_id', 'inputnode.subject_id'), ('outputnode.subjects_dir', 'inputnode.subjects_dir'), + ('outputnode.sulc', 'inputnode.sulc'), ]), (sphere_reg_wf, outputnode, [ ('outputnode.sphere_reg', 'sphere_reg'), diff --git a/smriprep/workflows/base.py b/smriprep/workflows/base.py index e39c3a2181..9ef3884b2f 100644 --- a/smriprep/workflows/base.py +++ b/smriprep/workflows/base.py @@ -52,6 +52,7 @@ def init_smriprep_wf( layout, longitudinal, low_mem, + msm_sulc, omp_nthreads, output_dir, run_uuid, @@ -90,6 +91,7 @@ def init_smriprep_wf( layout=BIDSLayout('.'), longitudinal=False, low_mem=False, + msm_sulc=False, omp_nthreads=1, output_dir='.', run_uuid='testrun', @@ -123,6 +125,8 @@ def init_smriprep_wf( See sub-workflows for specific differences low_mem : :obj:`bool` Write uncompressed .nii files in some cases to reduce memory usage + msm_sulc : :obj:`bool` + Run Multimodal Surface Matching with sulcal depth maps omp_nthreads : :obj:`int` Maximum number of threads an individual process may use output_dir : :obj:`str` @@ -176,6 +180,7 @@ def init_smriprep_wf( layout=layout, longitudinal=longitudinal, low_mem=low_mem, + msm_sulc=msm_sulc, name="single_subject_%s_wf" % subject_id, omp_nthreads=omp_nthreads, output_dir=output_dir, @@ -210,6 +215,7 @@ def init_single_subject_wf( layout, longitudinal, low_mem, + msm_sulc, name, omp_nthreads, output_dir, @@ -250,6 +256,7 @@ def init_single_subject_wf( layout=BIDSLayout('.'), longitudinal=False, low_mem=False, + msm_sulc=False, name='single_subject_wf', omp_nthreads=1, output_dir='.', @@ -409,6 +416,7 @@ def init_single_subject_wf( freesurfer=freesurfer, hires=hires, longitudinal=longitudinal, + msm_sulc=msm_sulc, name="anat_preproc_wf", t1w=subject_data["t1w"], t2w=subject_data["t2w"], diff --git a/smriprep/workflows/surfaces.py b/smriprep/workflows/surfaces.py index 71b4dbecad..50c0ad8ae3 100644 --- a/smriprep/workflows/surfaces.py +++ b/smriprep/workflows/surfaces.py @@ -49,6 +49,7 @@ PatchedRobustRegister as RobustRegister, RefineBrainMask, ) +import templateflow.api as tf from ..interfaces.workbench import CreateSignedDistanceVolume @@ -493,7 +494,7 @@ def _dedup(in_list): return workflow -def init_sphere_reg_wf(*, name="sphere_reg_wf"): +def init_sphere_reg_wf(*, msm_sulc: bool = False, name: str = "sphere_reg_wf"): """Generate GIFTI registration files to fsLR space""" from ..interfaces.surf import FixGiftiMetadata from ..interfaces.workbench import SurfaceSphereProjectUnproject @@ -501,7 +502,7 @@ def init_sphere_reg_wf(*, name="sphere_reg_wf"): workflow = Workflow(name=name) inputnode = pe.Node( - niu.IdentityInterface(["subjects_dir", "subject_id"]), + niu.IdentityInterface(["subjects_dir", "subject_id", "sulc"]), name="inputnode", ) outputnode = pe.Node( @@ -513,11 +514,11 @@ def init_sphere_reg_wf(*, name="sphere_reg_wf"): # Via FreeSurfer2CaretConvertAndRegisterNonlinear.sh#L270-L273 # # See https://github.com/DCAN-Labs/DCAN-HCP/tree/9291324 - sphere_gii = pe.MapNode( - fs.MRIsConvert(out_datatype="gii"), iterfield="in_file", name="sphere_gii" + sphere_reg_gii = pe.MapNode( + fs.MRIsConvert(out_datatype="gii"), iterfield="in_file", name="sphere_reg_gii" ) - fix_meta = pe.MapNode(FixGiftiMetadata(), iterfield="in_file", name="fix_meta") + fix_reg_meta = pe.MapNode(FixGiftiMetadata(), iterfield="in_file", name="fix_reg_meta") # Via # ${CARET7DIR}/wb_command -surface-sphere-project-unproject @@ -546,14 +547,120 @@ def init_sphere_reg_wf(*, name="sphere_reg_wf"): ('subjects_dir', 'subjects_dir'), ('subject_id', 'subject_id'), ]), - (get_surfaces, sphere_gii, [(('sphere_reg', _sorted_by_basename), 'in_file')]), - (sphere_gii, fix_meta, [('converted', 'in_file')]), - (fix_meta, project_unproject, [('out_file', 'sphere_in')]), - (sphere_gii, outputnode, [('converted', 'sphere_reg')]), - (project_unproject, outputnode, [('sphere_out', 'sphere_reg_fsLR')]), + (get_surfaces, sphere_reg_gii, [(('sphere_reg', _sorted_by_basename), 'in_file')]), + (sphere_reg_gii, fix_reg_meta, [('converted', 'in_file')]), + (fix_reg_meta, project_unproject, [('out_file', 'sphere_in')]), + (fix_reg_meta, outputnode, [('out_file', 'sphere_reg')]), ]) # fmt:on + if not msm_sulc: + workflow.connect(project_unproject, 'sphere_out', outputnode, 'sphere_reg_fsLR') + return workflow + + sphere_gii = pe.MapNode( + fs.MRIsConvert(out_datatype='gii'), iterfield='in_file', name='sphere_gii', + ) + fix_sphere_meta = pe.MapNode( + FixGiftiMetadata(), iterfield='in_file', name='fix_sphere_meta', + ) + msm_sulc_wf = init_msm_sulc_wf() + # fmt:off + workflow.connect([ + (get_surfaces, sphere_gii, [(('sphere', _sorted_by_basename), 'in_file')]), + (sphere_gii, fix_sphere_meta, [('converted', 'in_file')]), + (fix_sphere_meta, msm_sulc_wf, [('out_file', 'inputnode.sphere')]), + (inputnode, msm_sulc_wf, [('sulc', 'inputnode.sulc')]), + (project_unproject, msm_sulc_wf, [('sphere_out', 'inputnode.sphere_reg_fsLR')]), + (msm_sulc_wf, outputnode, [('outputnode.sphere_reg_fsLR', 'sphere_reg_fsLR')]), + ]) + # fmt:on + return workflow + + +def init_msm_sulc_wf(*, name: str = 'msm_sulc_wf'): + """Run MSMSulc registration to fsLR surfaces, per hemisphere.""" + from ..interfaces.msm import MSM + from ..interfaces.workbench import SurfaceAffineRegression, SurfaceApplyAffine + + workflow = Workflow(name=name) + inputnode = pe.Node( + niu.IdentityInterface(fields=['sulc', 'sphere', 'sphere_reg_fsLR']), + name='inputnode', + ) + outputnode = pe.Node(niu.IdentityInterface(fields=['sphere_reg_fsLR']), name='outputnode') + + # 0) Calculate affine + # ${CARET7DIR}/wb_command -surface-affine-regression \ + # $SUB.L.sphere.native.surf.gii \ + # $SUB.sphere.reg.reg_LR.native.surf.gii \ + # "$AtlasSpaceFolder"/"$NativeFolder"/MSMSulc/${Hemisphere}.mat + regress_affine = pe.MapNode( + SurfaceAffineRegression(), + iterfield=['in_surface', 'target_surface'], + name='regress_affine', + ) + + # 1) Apply affine to native sphere: + # wb_command -surface-apply-affine \ + # ${SUB}.L.sphere.native.surf.gii \ + # L.mat \ + # ${SUB}.L.sphere_rot.native.surf.gii + apply_surface_affine = pe.MapNode( + SurfaceApplyAffine(), + iterfield=['in_surface', 'in_affine'], + name='apply_surface_affine', + ) + # 2) Run MSMSulc + # ./msm_centos_v3 --conf=MSMSulcStrainFinalconf \ + # --inmesh=${SUB}.${HEMI}.sphere_rot.native.surf.gii + # --refmesh=fsaverage.${HEMI}_LR.spherical_std.164k_fs_LR.surf.gii + # --indata=sub-${SUB}_ses-${SES}_hemi-${HEMI)_sulc.shape.gii \ + # --refdata=tpl-fsaverage_hemi-${HEMI}_den-164k_sulc.shape.gii \ + # --out=${HEMI}. --verbose + msmsulc = pe.MapNode( + MSM(verbose=True, config_file=load_resource('msm/MSMSulcStrainFinalconf')), + iterfield=['in_mesh', 'reference_mesh', 'in_data', 'reference_data', 'out_base'], + name='msmsulc', + mem_gb=2, + ) + msmsulc.inputs.out_base = ['lh.', 'rh.'] # To placate Path2BIDS + msmsulc.inputs.reference_mesh = [ + str( + tf.get( + 'fsaverage', + hemi=hemi, + density='164k', + desc='std', + suffix='sphere', + extension='.surf.gii', + ) + ) + for hemi in 'LR' + ] + msmsulc.inputs.reference_data = [ + str( + tf.get( + 'fsaverage', + hemi=hemi, + density='164k', + suffix='sulc', + extension='.shape.gii', + ) + ) + for hemi in 'LR' + ] + # fmt:off + workflow.connect([ + (inputnode, regress_affine, [('sphere', 'in_surface'), + ('sphere_reg_fsLR', 'target_surface')]), + (inputnode, apply_surface_affine, [('sphere', 'in_surface')]), + (regress_affine, apply_surface_affine, [('out_affine', 'in_affine')]), + (inputnode, msmsulc, [('sulc', 'in_data')]), + (apply_surface_affine, msmsulc, [('out_surface', 'in_mesh')]), + (msmsulc, outputnode, [('warped_mesh', 'sphere_reg_fsLR')]), + ]) + # fmt:on return workflow @@ -985,9 +1092,9 @@ def init_morph_grayords_wf( Paths to JSON files containing metadata corresponding to ``cifti_morph`` """ - import templateflow.api as tf from niworkflows.engine.workflows import LiterateWorkflow as Workflow from smriprep.interfaces.cifti import GenerateDScalar + import templateflow.api as tf workflow = Workflow(name=name) workflow.__desc__ = f"""\ diff --git a/wrapper/smriprep_docker.py b/wrapper/smriprep_docker.py index 5edc2b30d6..4ccf0d12f9 100755 --- a/wrapper/smriprep_docker.py +++ b/wrapper/smriprep_docker.py @@ -34,7 +34,7 @@ MISSING = """ Image '{}' is missing Would you like to download? [Y/n] """ -PKG_PATH = '/opt/conda/lib/python3.9/site-packages' +PKG_PATH = '/opt/conda/envs/smriprep/lib/python3.10/site-packages' # Monkey-patch Py2 subprocess if not hasattr(subprocess, "DEVNULL"):