ENH: Add Multimodal Surface Matching

.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}

.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}

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 \

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():

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__ = [

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,

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

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}
+}

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

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