diff --git a/.vscode/settings.json b/.vscode/settings.json new file mode 100644 index 0000000000..5b94359216 --- /dev/null +++ b/.vscode/settings.json @@ -0,0 +1,3 @@ +{ + "esbonio.server.enabled": false +} \ No newline at end of file diff --git a/.zenodo.json b/.zenodo.json index eec1b4536c..b6c6e29c67 100644 --- a/.zenodo.json +++ b/.zenodo.json @@ -47,6 +47,11 @@ "orcid": "0000-0002-5312-6729", "type": "Researcher" }, + { + "affiliation": "Beckman Institute for Advanced Science & Technology, University of Illinois at Urbana-Champaign, IL, USA", + "name": "Camacho, Paul B.", + "orcid": "0000-0001-9048-7307", + }, { "affiliation": "Department of Pediatrics, University of Minnesota, MN, USA", "name": "Feczko, Eric", diff --git a/Dockerfile b/Dockerfile index b19b2deea6..54a0d7005c 100644 --- a/Dockerfile +++ b/Dockerfile @@ -33,6 +33,12 @@ RUN python -m build /src/fmriprep FROM ubuntu:jammy-20221130 # Prepare environment +ENV DEBIAN_FRONTEND="noninteractive" +ENV LANG=C.UTF-8 +ARG PYTHON_VERSION=3.8 +ARG CONDA_FILE=Miniconda3-py38_4.11.0-Linux-x86_64.sh +ENV DEBIAN_FRONTEND=noninteractive + RUN apt-get update && \ apt-get install -y --no-install-recommends \ apt-utils \ @@ -41,6 +47,9 @@ RUN apt-get update && \ bzip2 \ ca-certificates \ curl \ + wget \ + upx \ + file \ git \ gnupg \ libtool \ @@ -51,15 +60,45 @@ RUN apt-get update && \ xvfb && \ apt-get clean && rm -rf /var/lib/apt/lists/* /tmp/* /var/tmp/* -ENV DEBIAN_FRONTEND="noninteractive" \ - LANG="en_US.UTF-8" \ - LC_ALL="en_US.UTF-8" +# git clone stable branch of FastSurfer +RUN cd /opt && mkdir /fastsurfer \ + && git clone -b v2.0.4 https://github.com/Deep-MI/FastSurfer.git \ + && cp /opt/FastSurfer/fastsurfer_env_gpu.yml /fastsurfer/fastsurfer_env_gpu.yml + +# Copy fastsurfer from git folder +RUN cp -R /opt/FastSurfer/* /fastsurfer/ && rm -rf /opt/FastSurfer + +# Install conda +RUN wget --no-check-certificate -qO ~/miniconda.sh https://repo.continuum.io/miniconda/$CONDA_FILE && \ + chmod +x ~/miniconda.sh && \ + ~/miniconda.sh -b -p /opt/conda && \ + rm ~/miniconda.sh # Installing freesurfer COPY docker/files/freesurfer7.3.2-exclude.txt /usr/local/etc/freesurfer7.3.2-exclude.txt RUN curl -sSL https://surfer.nmr.mgh.harvard.edu/pub/dist/freesurfer/7.3.2/freesurfer-linux-ubuntu22_amd64-7.3.2.tar.gz \ | tar zxv --no-same-owner -C /opt --exclude-from=/usr/local/etc/freesurfer7.3.2-exclude.txt +# Install required packages for freesurfer to run +RUN apt-get update && apt-get install -y --no-install-recommends \ + tcsh \ + time \ + bc \ + gawk \ + libgomp1 && \ + apt clean && \ + rm -rf /var/lib/apt/lists/* /tmp/* /var/tmp/* + +# Add FreeSurfer Environment variables +ENV OS=Linux \ + FS_OVERRIDE=0 \ + FIX_VERTEX_AREA= \ + SUBJECTS_DIR=/opt/freesurfer/subjects \ + FSF_OUTPUT_FORMAT=nii.gz \ + FREESURFER_HOME=/opt/freesurfer \ + PYTHONUNBUFFERED=0 \ + PATH=/opt/freesurfer/bin:$PATH + # Simulate SetUpFreeSurfer.sh ENV FSL_DIR="/opt/fsl-6.0.5.1" \ OS="Linux" \ @@ -89,7 +128,6 @@ RUN apt-get update -qq \ libgl1-mesa-dev \ libgl1-mesa-dri \ libglu1-mesa-dev \ - libgomp1 \ libice6 \ libxcursor1 \ libxft2 \ @@ -98,6 +136,7 @@ RUN apt-get update -qq \ libxrender1 \ libxt6 \ sudo \ + curl \ wget \ && apt-get clean \ && rm -rf /var/lib/apt/lists/* \ @@ -167,21 +206,16 @@ RUN GNUPGHOME=/tmp gpg --keyserver hkps://keyserver.ubuntu.com --no-default-keyr # AFNI latest (neurodocker build) RUN apt-get update -qq \ && apt-get install -y -q --no-install-recommends \ - apt-utils \ ed \ gsl-bin \ libglib2.0-0 \ - libglu1-mesa-dev \ libglw1-mesa \ - libgomp1 \ libjpeg62 \ libpng12-0 \ libxm4 \ libxp6 \ netpbm \ - tcsh \ xfonts-base \ - xvfb \ && apt-get clean \ && rm -rf /var/lib/apt/lists/* \ && curl -sSL --retry 5 -o /tmp/multiarch.deb http://archive.ubuntu.com/ubuntu/pool/main/g/glibc/multiarch-support_2.27-3ubuntu1.5_amd64.deb \ @@ -281,10 +315,16 @@ RUN /opt/conda/bin/python fetch_templates.py && \ find $HOME/.cache/templateflow -type d -exec chmod go=u {} + && \ find $HOME/.cache/templateflow -type f -exec chmod go=u {} + +# Installing FastSurfer gpu dependencies +RUN conda env update -n base --file /fastsurfer/fastsurfer_env_gpu.yml + # Installing sMRIPREP COPY --from=src /src/fmriprep/dist/*.whl . RUN /opt/conda/bin/python -m pip install --no-cache-dir $( ls *.whl )[all] +# Installing nibabel version 4.0.2 +RUN /opt/conda/bin/python -m pip install --no-cache-dir nibabel==4.0.2 + RUN find $HOME -type d -exec chmod go=u {} + && \ find $HOME -type f -exec chmod go=u {} + && \ rm -rf $HOME/.npm $HOME/.conda $HOME/.empty @@ -294,6 +334,12 @@ RUN find $HOME -type d -exec chmod go=u {} + && \ ENV FREESURFER="/opt/freesurfer" ENV IS_DOCKER_8395080871=1 +ENV FASTSURFER_HOME=/fastsurfer +ENV PATH="$PATH:/fastsurfer" + +# Download all remote network checkpoints already +ENV PYTHONPATH=/fastsurfer:$PYTHONPATH +RUN cd /fastsurfer ; python3 FastSurferCNN/download_checkpoints.py --all RUN ldconfig WORKDIR /tmp diff --git a/docker/files/freesurfer7.3.2-exclude.txt b/docker/files/freesurfer7.3.2-exclude.txt index a817a062b1..09714a0bc1 100644 --- a/docker/files/freesurfer7.3.2-exclude.txt +++ b/docker/files/freesurfer7.3.2-exclude.txt @@ -4,20 +4,10 @@ freesurfer/average/3T18yoSchwartzReactN32_as_orig.4dfp.img freesurfer/average/3T18yoSchwartzReactN32_as_orig.4dfp.img.rec freesurfer/average/3T18yoSchwartzReactN32_as_orig.4dfp.mat freesurfer/average/3T18yoSchwartzReactN32_as_orig.lst -freesurfer/average/711-2B_as_mni_average_305.4dfp.hdr -freesurfer/average/711-2B_as_mni_average_305.4dfp.ifh -freesurfer/average/711-2B_as_mni_average_305.4dfp.img -freesurfer/average/711-2B_as_mni_average_305.4dfp.img.rec -freesurfer/average/711-2B_as_mni_average_305_mask.4dfp.hdr -freesurfer/average/711-2B_as_mni_average_305_mask.4dfp.img.rec -freesurfer/average/711-2C_as_mni_average_305.4dfp.hdr -freesurfer/average/711-2C_as_mni_average_305.4dfp.img.rec -freesurfer/average/711-2C_as_mni_average_305.4dfp.mat freesurfer/average/aseg+spmhead+vermis+pons.ixi.gca freesurfer/average/BrainstemSS freesurfer/average/Buckner_JNeurophysiol11_MNI152 freesurfer/average/Choi_JNeurophysiol12_MNI152 -freesurfer/average/colortable_desikan_killiany.txt freesurfer/average/face.gca freesurfer/average/HippoSF freesurfer/average/label_scales.dat @@ -41,7 +31,6 @@ freesurfer/average/mni305.cor.subfov1.mgz freesurfer/average/mni305.cor.subfov1.reg freesurfer/average/mni305.cor.subfov2.mgz freesurfer/average/mni305.cor.subfov2.reg -freesurfer/average/mni305.mask.cor.mgz freesurfer/average/mni_average_305.4dfp.hdr freesurfer/average/mni_average_305.4dfp.ifh freesurfer/average/mni_average_305.4dfp.img @@ -420,7 +409,6 @@ freesurfer/bin/mri_make_template freesurfer/bin/mri_map_cpdat freesurfer/bin/mri_maps2csd freesurfer/bin/mri_mark_temporal_lobe -freesurfer/bin/mri_mc freesurfer/bin/mri_mcsim freesurfer/bin/mri_mergelabels freesurfer/bin/mri_mi @@ -504,7 +492,6 @@ freesurfer/bin/mris_gradient freesurfer/bin/mris_hausdorff_dist freesurfer/bin/mris_image2vtk freesurfer/bin/mri_simulate_atrophy -freesurfer/bin/mris_info freesurfer/bin/mris_init_global_tractography freesurfer/bin/mris_intensity_profile freesurfer/bin/mris_interpolate_warp @@ -548,7 +535,6 @@ freesurfer/bin/mris_rf_label freesurfer/bin/mris_rf_train freesurfer/bin/mris_rotate freesurfer/bin/mris_sample_label -freesurfer/bin/mris_sample_parc freesurfer/bin/mris_seg2annot freesurfer/bin/mris_segment freesurfer/bin/mris_segmentation_stats @@ -749,11 +735,9 @@ freesurfer/bin/xhemireg freesurfer/bin/xhemi-tal freesurfer/bin/xsanatreg freesurfer/bin/zero_lt_4dfp -freesurfer/DefectLUT.txt freesurfer/diffusion freesurfer/docs/xml freesurfer/FreeSurferEnv.csh -freesurfer/FreeSurferEnv.sh freesurfer/fsfast freesurfer/lib/bem/ic0.tri freesurfer/lib/bem/ic1.tri @@ -861,11 +845,8 @@ freesurfer/python/lib/python3.8/site-packages/torch freesurfer/python/lib/python3.8/site-packages/**/tests freesurfer/python/**/__pycache__ freesurfer/python/share -freesurfer/SegmentNoLUT.txt freesurfer/sessions freesurfer/SetUpFreeSurfer.csh -freesurfer/SetUpFreeSurfer.sh -freesurfer/Simple_surface_labels2009.txt freesurfer/sources.sh freesurfer/subjects/bert freesurfer/subjects/cvs_avg35 diff --git a/docs/installation.rst b/docs/installation.rst index a4c6d204b1..148241d396 100644 --- a/docs/installation.rst +++ b/docs/installation.rst @@ -148,5 +148,6 @@ the ``smriprep`` package: - ANTs_ (version 2.2.0 - NeuroDocker build) - AFNI_ (version Debian-16.2.07) - `C3D `_ (version 1.0.0) -- FreeSurfer_ (version 6.0.1) +- FreeSurfer_ (version 7.2.0) - `bids-validator `_ (version 1.1.0) +- FastSurfer_ (version 2.0.1) diff --git a/docs/requirements.txt b/docs/requirements.txt index 20f6906786..fe5957b1cf 100644 --- a/docs/requirements.txt +++ b/docs/requirements.txt @@ -10,4 +10,4 @@ sphinx sphinx_rtd_theme sphinxcontrib-apidoc ~= 0.3.0 templateflow -networkx != 2.8.1 +networkx != 2.8.1 \ No newline at end of file diff --git a/smriprep/cli/run.py b/smriprep/cli/run.py index 4a588283f8..2db66800f7 100644 --- a/smriprep/cli/run.py +++ b/smriprep/cli/run.py @@ -222,14 +222,21 @@ def get_parser(): dest="hires", help="disable sub-millimeter (hires) reconstruction", ) - g_surfs_xor = g_surfs.add_mutually_exclusive_group() + g_surfs_xor = g_surfs.add_mutually_exclusive_group() g_surfs_xor.add_argument( "--fs-no-reconall", action="store_false", dest="run_reconall", help="disable FreeSurfer surface preprocessing.", ) + g_surfs_xor.add_argument( + "--fastsurfer", + action="store_true", + default=False, + dest="fastsurfer", + help="enable FastSurfer surface preprocessing.", + ) g_other = parser.add_argument_group("Other options") g_other.add_argument( @@ -394,6 +401,17 @@ def _warn_redirect(message, category, filename, lineno, file=None, line=None): from templateflow import api from niworkflows.utils.misc import _copy_any + dseg_tsv = str(api.get("fsaverage", suffix="dseg", extension=[".tsv"])) + _copy_any( + dseg_tsv, str(Path(output_dir) / "smriprep" / "desc-aseg_dseg.tsv") + ) + _copy_any( + dseg_tsv, str(Path(output_dir) / "smriprep" / "desc-aparcaseg_dseg.tsv") + ) + elif opts.fastsurfer: + from templateflow import api + from niworkflows.utils.misc import _copy_any + dseg_tsv = str(api.get("fsaverage", suffix="dseg", extension=[".tsv"])) _copy_any( dseg_tsv, str(Path(output_dir) / "smriprep" / "desc-aseg_dseg.tsv") @@ -584,6 +602,7 @@ def build_workflow(opts, retval): freesurfer=opts.run_reconall, fs_subjects_dir=opts.fs_subjects_dir, hires=opts.hires, + fastsurfer=opts.fastsurfer, layout=layout, longitudinal=opts.longitudinal, low_mem=opts.low_mem, diff --git a/smriprep/data/boilerplate.bib b/smriprep/data/boilerplate.bib index 9acffc48dc..730f47c74c 100644 --- a/smriprep/data/boilerplate.bib +++ b/smriprep/data/boilerplate.bib @@ -74,6 +74,34 @@ @article{fs_reconall year = 1999 } +@article{fastsurfer, +title = {FastSurfer - A fast and accurate deep learning based neuroimaging pipeline}, +journal = {NeuroImage}, +volume = {219}, +pages = {117012}, +year = {2020}, +issn = {1053-8119}, +doi = {https://doi.org/10.1016/j.neuroimage.2020.117012}, +url = {https://www.sciencedirect.com/science/article/pii/S1053811920304985}, +author = {Leonie Henschel and Sailesh Conjeti and Santiago Estrada and Kersten Diers and Bruce Fischl and Martin Reuter}, +keywords = {Freesurfer, Computational neuroimaging, Deep learning, Structural MRI, Artificial intelligence}, +abstract = {Traditional neuroimage analysis pipelines involve computationally intensive, time-consuming optimization steps, and thus, do not scale well to large cohort studies with thousands or tens of thousands of individuals. In this work we propose a fast and accurate deep learning based neuroimaging pipeline for the automated processing of structural human brain MRI scans, replicating FreeSurfer’s anatomical segmentation including surface reconstruction and cortical parcellation. To this end, we introduce an advanced deep learning architecture capable of whole-brain segmentation into 95 classes. The network architecture incorporates local and global competition via competitive dense blocks and competitive skip pathways, as well as multi-slice information aggregation that specifically tailor network performance towards accurate segmentation of both cortical and subcortical structures. Further, we perform fast cortical surface reconstruction and thickness analysis by introducing a spectral spherical embedding and by directly mapping the cortical labels from the image to the surface. This approach provides a full FreeSurfer alternative for volumetric analysis (in under 1 ​min) and surface-based thickness analysis (within only around 1 ​h runtime). For sustainability of this approach we perform extensive validation: we assert high segmentation accuracy on several unseen datasets, measure generalizability and demonstrate increased test-retest reliability, and high sensitivity to group differences in dementia.} +} + +@article{fastsufrer_hires, +title = {FastSurferVINN: Building resolution-independence into deep learning segmentation methods—A solution for HighRes brain MRI}, +journal = {NeuroImage}, +volume = {251}, +pages = {118933}, +year = {2022}, +issn = {1053-8119}, +doi = {https://doi.org/10.1016/j.neuroimage.2022.118933}, +url = {https://www.sciencedirect.com/science/article/pii/S1053811922000623}, +author = {Leonie Henschel and David Kügler and Martin Reuter}, +keywords = {Computational neuroimaging, Deep learning, Structural MRI, Artificial intelligence, High-resolution}, +abstract = {Leading neuroimaging studies have pushed 3T MRI acquisition resolutions below 1.0 mm for improved structure definition and morphometry. Yet, only few, time-intensive automated image analysis pipelines have been validated for high-resolution (HiRes) settings. Efficient deep learning approaches, on the other hand, rarely support more than one fixed resolution (usually 1.0 mm). Furthermore, the lack of a standard submillimeter resolution as well as limited availability of diverse HiRes data with sufficient coverage of scanner, age, diseases, or genetic variance poses additional, unsolved challenges for training HiRes networks. Incorporating resolution-independence into deep learning-based segmentation, i.e., the ability to segment images at their native resolution across a range of different voxel sizes, promises to overcome these challenges, yet no such approach currently exists. We now fill this gap by introducing a Voxel-size Independent Neural Network (VINN) for resolution-independent segmentation tasks and present FastSurferVINN, which (i) establishes and implements resolution-independence for deep learning as the first method simultaneously supporting 0.7–1.0 mm whole brain segmentation, (ii) significantly outperforms state-of-the-art methods across resolutions, and (iii) mitigates the data imbalance problem present in HiRes datasets. Overall, internal resolution-independence mutually benefits both HiRes and 1.0 mm MRI segmentation. With our rigorously validated FastSurferVINN we distribute a rapid tool for morphometric neuroimage analysis. The VINN architecture, furthermore, represents an efficient resolution-independent segmentation method for wider application.} +} + @article{mindboggle, author = {Klein, Arno and Ghosh, Satrajit S. and Bao, Forrest S. and Giard, Joachim and Häme, Yrjö and Stavsky, Eliezer and Lee, Noah and Rossa, Brian and Reuter, Martin and Neto, Elias Chaibub and Keshavan, Anisha}, doi = {10.1371/journal.pcbi.1005350}, diff --git a/smriprep/interfaces/fastsurfer.py b/smriprep/interfaces/fastsurfer.py new file mode 100644 index 0000000000..9833eaf483 --- /dev/null +++ b/smriprep/interfaces/fastsurfer.py @@ -0,0 +1,399 @@ +# -*- coding: utf-8 -*- +# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- +# vi: set ft=python sts=4 ts=4 sw=4 et: +""" +The FastSufer module provides basic functions +for running FastSurfer VINN and surface processing. + +""" +from nipype.interfaces.base import ( + CommandLine, + Directory, + CommandLineInputSpec, + OutputMultiPath, + BaseInterfaceInputSpec, + File, +) +from nipype.interfaces.base.traits_extension import traits +from nipype.interfaces.io import FSSourceOutputSpec, FreeSurferSource + + +class FastSInputSpec(CommandLineInputSpec): + r""" + Required arguments + ================== + + subjects_dir + Output directory + subject_id + Subject ID for directory inside ``sd`` to be created + t1 + T1 full head input, not bias corrected, global path. + The 'network was trained with conformed images + UCHAR, 256 x 256 x 256, 1 mm voxels and standard slice orientation. + These specifications are checked in the ``eval.py`` script and the image + is automatically conformed if it does not comply. + fs_license + Path to FreeSurfer license key file. + + + Optional arguments + ================== + + Network specific arguments + -------------------------- + seg + Global path with filename of segmentation + (where and under which name to store it). + Default location + ``$SUBJECTS_DIR/$sid/mri/aparc.DKTatlas+aseg.deep.mgz`` + weights_sag + Pretrained weights of sagittal network. + Default + ``../checkpoints/Sagittal_Weights_FastSurferCNN/ckpts/Epoch_30_training_state.pkl`` + weights_ax + Pretrained weights of axial network. + Default + ``../checkpoints/Axial_Weights_FastSurferCNN/ckpts/Epoch_30_training_state.pkl`` + weights_cor + Pretrained weights of coronal network. + Default + ``../checkpoints/Coronal_Weights_FastSurferCNN/ckpts/Epoch_30_training_state.pkl`` + seg_log + Name and location for the log-file for the segmentation (FastSurferCNN). + Default '$SUBJECTS_DIR/$sid/scripts/deep-seg.log' + clean_seg + Flag to clean up FastSurferCNN segmentation + run_viewagg_on + Define where the view aggregation should be run on. + By default, the program checks if you have enough memory to run + the view aggregation on the gpu. The total memory is considered for this decision. + If this fails, or you actively overwrote the check with setting + ``run_viewagg_on cpu``, view agg is run on the cpu. + Equivalently, if you define ``--run_viewagg_on gpu``, view agg will be run on the gpu + (no memory check will be done). + no_cuda + Flag to disable CUDA usage in FastSurferCNN (no GPU usage, inference on CPU). + batch + Batch size for inference. Default 16. Lower this to reduce memory requirement. + + + Surface pipeline arguments + -------------------------- + fstess + Use ``mri_tesselate`` instead of marching cube (default) for surface creation + fsqsphere + Use FreeSurfer default instead of + ndovel spectral spherical projection for qsphere + fsaparc + Use FS aparc segmentations in addition to DL prediction + (slower in this case and usually the mapped ones from the DL prediction are fine) + no_surfreg + Skip creating Surface-Atlas ``sphere.reg`` registration with FreeSurfer + (for cross-subject correspondence or other mappings) + parallel + Run both hemispheres in parallel + threads + Set openMP and ITK threads + + + Other + ---- + py + which python version to use. + Default ``python3.8`` + seg_only + only run FastSurferCNN + (generate segmentation, do not run the surface pipeline) + surf_only + only run the surface pipeline ``recon_surf``. + The segmentation created by FastSurferCNN must already exist in this case. + + + """ + subjects_dir = Directory( + exists=True, + argstr="--sd %s", + hash_files=False, + desc="Subjects directory", + genfile=True, + ) + subject_id = traits.String( + argstr="--sid %s", + mandatory=True, + desc="Subject ID" + ) + T1_files = File( + exists=True, + mandatory=False, + argstr="--t1 %s", + desc="T1 full head input (not bias corrected, global path)" + ) + fs_license = File( + exists=True, + argstr="--fs_license %s", + desc="Path to FreeSurfer license key file." + ) + seg = File( + exists=True, + argstr="--seg %s", + desc="Pre-computed segmentation file" + ) + weights_sag = File( + exists=True, + mandatory=False, + argstr="--weights_sag %s", + desc="Pretrained weights of sagittal network" + ) + weights_ax = File( + exists=True, + mandatory=False, + argstr="--weights_ax %s", + desc="Pretrained weights of axial network" + ) + weights_cor = File( + exists=True, + mandatory=False, + argstr="--weights_cor %s", + desc="Pretrained weights of coronal network" + ) + seg_log = File( + exists=True, + mandatory=False, + argstr="--seg_log %s", + desc="Name and location for the log-file for the segmentation (FastSurferCNN)." + ) + clean_seg = traits.Bool( + False, + mandatory=False, + argstr="--clean_seg", + desc="Flag to clean up FastSurferCNN segmentation" + ) + run_viewagg_on = File( + exists=True, + mandatory=False, + argstr="--run_viewagg_on %s", + desc="Define where the view aggregation should be run on." + ) + no_cuda = traits.Bool( + False, + mandatory=False, + argstr="--no_cuda", + desc="Flag to disable CUDA usage in FastSurferCNN (no GPU usage, inference on CPU)" + ) + batch = traits.Int( + 16, + usedefault=True, + mandatory=False, + argstr="--batch %d", + desc="Batch size for inference. default=16. Lower this to reduce memory requirement" + ) + fstess = traits.Bool( + False, + mandatory=False, + argstr="--fstess", + desc="Use mri_tesselate instead of marching cube (default) for surface creation" + ) + fsqsphere = traits.Bool( + False, + mandatory=False, + argstr="--fsqsphere", + desc="Use FreeSurfer default instead of novel spectral spherical projection for qsphere" + ) + fsaparc = traits.Bool( + False, + mandatory=False, + argstr="--fsaparc", + desc="Use FS aparc segmentations in addition to DL prediction" + ) + no_surfreg = traits.Bool( + False, + mandatory=False, + argstr="--no_surfreg", + desc="""Skip creating Surface-Atlas (sphere.reg) registration with FreeSurfer + (for cross-subject correspondence or other mappings)""" + ) + parallel = traits.Bool( + True, + usedefault=True, + mandatory=False, + argstr="--parallel", + desc="Run both hemispheres in parallel" + ) + threads = traits.Int( + 4, + usedefault=True, + mandatory=False, + argstr="--threads %d", + desc="Set openMP and ITK threads to" + ) + py = traits.String( + "python3.8", + usedefault=True, + mandatory=False, + argstr="--py %s", + desc="which python version to use. default=python3.6" + ) + seg_only = traits.Bool( + False, + mandatory=False, + argstr="--seg_only", + desc="only run FastSurferCNN (generate segmentation, do not surface)" + ) + surf_only = traits.Bool( + False, + mandatory=False, + argstr="--surf_only", + desc="only run the surface pipeline recon_surf." + ) + + +class FastSurfSourceInputSpec(BaseInterfaceInputSpec): + sd = Directory( + exists=True, + argstr="--sd %s", + mandatory=False, + desc="Subjects directory" + ) + sid = traits.String( + exists=True, + argstr="--sid %s", + mandatory=False, + desc="Subject ID" + ) + t1 = File( + exists=True, + mandatory=False, + argstr="--t1 %s", + desc="T1 full head input (not bias corrected, global path)" + ) + + +class FastSurferSourceOutputSpec(FSSourceOutputSpec): + orig_nu = File( + exists=True, + desc="Base image conformed to Fastsurfer space and nonuniformity corrected", + loc="mri" + ) + wm_asegedit = OutputMultiPath( + File(exists=True), + desc="Edited white matter volume post-aseg", + loc="mri", + altkey="mri_nu_correct.mni" + ) + wmparc_mapped = OutputMultiPath( + exists=True, + loc="mri", + desc="DKT atlas mapped Aparc into subcortical white matter", + altkey="wmparc.DKTatlas.mapped" + ) + mni_log = OutputMultiPath( + File(exists=True), + desc="Non-uniformity correction log", + loc="mri", + altkey="mri_nu_correct.mni" + ) + mni_log_bak = OutputMultiPath( + File(exists=True), + desc="Non-uniformity correction log bak", + loc="mri", + altkey="mri_nu_correct.mni" + ) + defects = OutputMultiPath( + File(exists=True), + desc="Defects", + loc="surf", + altkey=["defect_borders", "defect_chull", "defect_labels"], + ) + nofix = OutputMultiPath( + File(exists=True), + desc="Pre-tessellation original surface", + loc="surf", + altkey=["nofix"], + ) + aparc_ctab = OutputMultiPath( + File(exists=True), + loc="label", + altkey="aparc.annot.mapped.ctab", + desc="Aparc parcellation annotation ctab file", + ) + mapped_024 = OutputMultiPath( + File(exists=True), + loc="label", + altkey="mapped*024", + desc="Mapped label files", + ) + cortex = OutputMultiPath( + File(exists=True), + loc="label", + altkey="cortex", + desc="Cortex class label files", + ) + aparc_dkt_aseg = OutputMultiPath( + File(exists=True), + loc="mri", + altkey="aparc.DKTatlas*aseg*", + desc="Aparc parcellation from DKT atlas projected into aseg volume", + ) + segment_dat = OutputMultiPath( + File(exists=True), + loc="mri", + altkey="segment_dat", + desc="Segmentation .dat files", + ) + filled_pretess = OutputMultiPath( + File(exists=True), + loc="mri", + altkey="filled*pretess*", + desc="Pre-tessellation filled volume files", + ) + preaparc = OutputMultiPath( + File(exists=True), + loc="surf", + altkey="preaparc", + desc="Pre-Aparc files", + ) + w_g_stats = OutputMultiPath( + File(exists=True), + loc="stats", + altkey="w*g.pct", + desc="White minus gray statistics files" + ) + aseg_presurf_stats = OutputMultiPath( + File(exists=True), + loc="stats", + altkey="aseg.presurf.hypos", + desc="Automated segmentation pre-surface recon statistics files" + ) + subjects_dir = Directory( + exists=True, + argstr="--sd %s", + desc="Subjects directory" + ) + subject_id = traits.String( + exists=True, + argstr="--sid %s", + desc="Subject ID" + ) + + +class FastSurferSource(FreeSurferSource): + """Generates FastSurfer subject info from their directories. + + """ + + output_spec = FastSurferSourceOutputSpec + + +class FastSurfer(CommandLine): + """Wraps FastSurfer command for segmentation and surface processing. + + """ + + input_spec = FastSInputSpec + output_spec = FastSurferSourceOutputSpec + _cmd = 'run_fastsurfer.sh' + + def _list_outputs(self): + outputs = self.output_spec().get() + return outputs diff --git a/smriprep/interfaces/reports.py b/smriprep/interfaces/reports.py index 7db4f92c1b..0ba0a77e81 100644 --- a/smriprep/interfaces/reports.py +++ b/smriprep/interfaces/reports.py @@ -37,6 +37,8 @@ ) from nipype.interfaces import freesurfer as fs from nipype.interfaces.io import FSSourceInputSpec as _FSSourceInputSpec +from smriprep.interfaces import fastsurfer as fastsurf +from smriprep.utils.misc import check_fastsurfer from nipype.interfaces.mixins import reporting from niworkflows.interfaces.reportlets.base import _SVGReportCapableInputSpec @@ -48,6 +50,7 @@ \t\t
  • Structural images: {n_t1s:d} T1-weighted {t2w}
    • \t\t
  • Standard spaces: {output_spaces}
    • \t\t
  • FreeSurfer reconstruction: {freesurfer_status}
    • +\t\t
  • FreeSurfer reconstruction: {fastsurfer_status}
    • \t """ @@ -99,6 +102,7 @@ class SubjectSummary(SummaryInterface): input_spec = _SubjectSummaryInputSpec output_spec = _SubjectSummaryOutputSpec + fastsurfer_bool = False def _run_interface(self, runtime): if isdefined(self.inputs.subject_id): @@ -106,15 +110,28 @@ def _run_interface(self, runtime): return super(SubjectSummary, self)._run_interface(runtime) def _generate_segment(self): + fastsurfer_status = "Not run" if not isdefined(self.inputs.subjects_dir): freesurfer_status = "Not run" + fastsurfer_status = "Not run" else: - recon = fs.ReconAll( + fastsurfer_bool = check_fastsurfer( subjects_dir=self.inputs.subjects_dir, - subject_id=self.inputs.subject_id, - T1_files=self.inputs.t1w, - flags="-noskullstrip", - ) + subject_id=self.inputs.subject_id) + if fastsurfer_bool is True: + recon = fastsurf.FastSCommand( + sd=self.inputs.subjects_dir, + sid=self.inputs.subject_id, + t1=self.inputs.t1w, + ) + fastsurfer_status = "Run by sMRIPrep" + else: + recon = fs.ReconAll( + subjects_dir=self.inputs.subjects_dir, + subject_id=self.inputs.subject_id, + T1_files=self.inputs.t1w, + flags="-noskullstrip", + ) if recon.cmdline.startswith("echo"): freesurfer_status = "Pre-existing directory" else: @@ -136,6 +153,7 @@ def _generate_segment(self): t2w=t2w_seg, output_spaces=output_spaces, freesurfer_status=freesurfer_status, + fastsurfer_status=fastsurfer_status, ) diff --git a/smriprep/utils/misc.py b/smriprep/utils/misc.py index 99d1458866..4145517381 100644 --- a/smriprep/utils/misc.py +++ b/smriprep/utils/misc.py @@ -95,3 +95,40 @@ def fs_isRunning(subjects_dir, subject_id, mtime_tol=86400, logger=None): if logger: logger.warn(f'Removed "IsRunning*" files found under {subj_dir}') return subjects_dir + + +def check_fastsurfer(subjects_dir, subject_id, logger=None): + """ + Checks FreeSurfer subjects dir for presence of files in mri/ with names \ + indicating processing with FastSurfer. + + Parameters + ---------- + subjects_dir : os.PathLike or None + Existing FreeSurfer subjects directory + subject_id : str + Subject label + + Returns + ------- + fastsurfer_bool : Boolean + + """ + from pathlib import Path + + fastsurfer_bool = False + + if subjects_dir is None: + return fastsurfer_bool + subj_dir = Path(subjects_dir) / subject_id + if not subj_dir.exists(): + return fastsurfer_bool + + fastsurferfiles = tuple(subj_dir.glob("mri/*deep*mgz")) + if not fastsurferfiles: + return fastsurfer_bool + else: + if logger: + logger.warn(f'Evidence of FastSurfer processing found in {subj_dir}') + fastsurfer_bool = True + return fastsurfer_bool diff --git a/smriprep/workflows/anatomical.py b/smriprep/workflows/anatomical.py index 571e011dc1..78e4ecdb64 100644 --- a/smriprep/workflows/anatomical.py +++ b/smriprep/workflows/anatomical.py @@ -33,7 +33,7 @@ ) from nipype.interfaces.ants.base import Info as ANTsInfo -from nipype.interfaces.ants import DenoiseImage, N4BiasFieldCorrection +from nipype.interfaces.ants import N4BiasFieldCorrection from niworkflows.engine.workflows import LiterateWorkflow as Workflow from niworkflows.interfaces.fixes import FixHeaderApplyTransforms as ApplyTransforms @@ -48,10 +48,15 @@ from niworkflows.utils.misc import fix_multi_T1w_source_name, add_suffix from niworkflows.anat.ants import init_brain_extraction_wf, init_n4_only_wf from ..utils.bids import get_outputnode_spec -from ..utils.misc import apply_lut as _apply_bids_lut, fs_isRunning as _fs_isRunning +from ..utils.misc import ( + apply_lut as _apply_bids_lut, + fs_isRunning as _fs_isRunning, + check_fastsurfer as _check_fastsurfer +) from .norm import init_anat_norm_wf from .outputs import init_anat_reports_wf, init_anat_derivatives_wf -from .surfaces import init_anat_ribbon_wf, init_surface_recon_wf, init_morph_grayords_wf +from .surfaces import init_surface_recon_wf, init_fastsurf_recon_wf, init_anat_ribbon_wf, init_morph_grayords_wf +from .surfaces import init_anat_ribbon_wf, init_surface_recon_wf LOGGER = logging.getLogger("nipype.workflow") @@ -62,6 +67,7 @@ def init_anat_preproc_wf( freesurfer, hires, longitudinal, + fastsurfer, t1w, t2w, omp_nthreads, @@ -100,6 +106,7 @@ def init_anat_preproc_wf( freesurfer=True, hires=True, longitudinal=False, + fastsurfer=False, t1w=['t1w.nii.gz'], t2w=[], omp_nthreads=1, @@ -125,6 +132,10 @@ def init_anat_preproc_wf( longitudinal : :obj:`bool` Create unbiased structural template, regardless of number of inputs (may increase runtime) + fastsurfer : :obj:`bool` + Enable FastSurfer surface reconstruction (shorter runtime than + the FreeSurfer workflow). Uses similar input variables and output + file structure as FreeSurfer. t1w : :obj:`list` List of T1-weighted structural images. omp_nthreads : :obj:`int` @@ -241,6 +252,7 @@ def init_anat_preproc_wf( # Connect reportlets workflows anat_reports_wf = init_anat_reports_wf( freesurfer=freesurfer, + fastsurfer=fastsurfer, output_dir=output_dir, ) # fmt:off @@ -457,10 +469,11 @@ def _check_img(img): ]) # fmt:on - # Write outputs ############################################3 + # Write outputs # anat_derivatives_wf = init_anat_derivatives_wf( bids_root=bids_root, freesurfer=freesurfer, + fastsurfer=fastsurfer, num_t1w=num_t1w, t2w=t2w, output_dir=output_dir, @@ -513,7 +526,7 @@ def _check_img(img): (fast2bids, outputnode, [('out', 't1w_tpms')]), ]) # fmt:on - if not freesurfer: # Flag --fs-no-reconall is set - return + if not freesurfer and not fastsurfer: # Flag --fs-no-reconall is set - return # fmt:off workflow.connect([ (brain_extraction_wf, buffernode, [ @@ -529,10 +542,50 @@ def _check_img(img): ) fs_isrunning.inputs.logger = LOGGER - # 5. Surface reconstruction (--fs-no-reconall not set) - surface_recon_wf = init_surface_recon_wf( - name="surface_recon_wf", omp_nthreads=omp_nthreads, hires=hires + # check for FastSurfer .mgz files and + # touch mri/aseg.auto_noCCseg.label_intensities.txt + # to prevent failure in surfaces.py (temporary fix - DEPRECATED) + check_fastsurfer = pe.Node( + niu.Function(function=_check_fastsurfer), overwrite=True, name="check_fastsurfer" ) + check_fastsurfer.inputs.logger = LOGGER + + # Select which surface reconstruction workflow based on CLI arguments + if freesurfer and (not fastsurfer): + surface_recon_wf = init_surface_recon_wf( + name="surface_recon_wf", omp_nthreads=omp_nthreads, hires=hires) + # fmt:off + workflow.connect([ + (inputnode, fs_isrunning, [ + ('subjects_dir', 'subjects_dir'), + ('subject_id', 'subject_id')]), + (inputnode, surface_recon_wf, [ + ('t2w', 'inputnode.t2w'), + ('flair', 'inputnode.flair'), + ('subject_id', 'inputnode.subject_id')]), + (fs_isrunning, surface_recon_wf, [('out', 'inputnode.subjects_dir')]), + (surface_recon_wf, anat_reports_wf, [ + ('outputnode.subject_id', 'inputnode.subject_id'), + ('outputnode.subjects_dir', 'inputnode.subjects_dir'), + ]), + ]) + # fmt:on + elif fastsurfer: + surface_recon_wf = init_fastsurf_recon_wf( + name="surface_recon_wf", omp_nthreads=omp_nthreads, hires=hires) + # fmt:off + workflow.connect([ + (inputnode, surface_recon_wf, [ + ('subject_id', 'inputnode.subject_id'), + ('subjects_dir', 'inputnode.subjects_dir')]), + (surface_recon_wf, anat_reports_wf, [ + ('outputnode.subject_id', 'inputnode.subject_id'), + ('outputnode.subjects_dir', 'inputnode.subjects_dir'), + ]), + ]) + # fmt:on + + # Identical connections applyrefined = pe.Node(fsl.ApplyMask(), name="applyrefined") if t2w: @@ -568,9 +621,9 @@ def _check_img(img): workflow.connect([ (inputnode, t2w_template_wf, [('t2w', 'inputnode.anat_files')]), (t2w_template_wf, bbreg, [('outputnode.anat_ref', 'source_file')]), - (surface_recon_wf, bbreg, [ - ('outputnode.subject_id', 'subject_id'), - ('outputnode.subjects_dir', 'subjects_dir'), + (inputnode, bbreg, [ + ('subject_id', 'subject_id'), + ('subjects_dir', 'subjects_dir'), ]), (bbreg, coreg_xfms, [('out_lta_file', 'in1')]), (surface_recon_wf, coreg_xfms, [('outputnode.fsnative2t1w_xfm', 'in2')]), @@ -587,15 +640,8 @@ def _check_img(img): # Anatomical ribbon file using HCP signed-distance volume method anat_ribbon_wf = init_anat_ribbon_wf() # fmt:off + workflow.connect([ - (inputnode, fs_isrunning, [ - ('subjects_dir', 'subjects_dir'), - ('subject_id', 'subject_id')]), - (inputnode, surface_recon_wf, [ - ('t2w', 'inputnode.t2w'), - ('flair', 'inputnode.flair'), - ('subject_id', 'inputnode.subject_id')]), - (fs_isrunning, surface_recon_wf, [('out', 'inputnode.subjects_dir')]), (anat_validate, surface_recon_wf, [('out_file', 'inputnode.t1w')]), (brain_extraction_wf, surface_recon_wf, [ (('outputnode.out_file', _pop), 'inputnode.skullstripped_t1'), @@ -622,13 +668,9 @@ def _check_img(img): (applyrefined, buffernode, [('out_file', 't1w_brain')]), (surface_recon_wf, buffernode, [ ('outputnode.out_brainmask', 't1w_mask')]), - (surface_recon_wf, anat_reports_wf, [ - ('outputnode.subject_id', 'inputnode.subject_id'), - ('outputnode.subjects_dir', 'inputnode.subjects_dir')]), (surface_recon_wf, anat_derivatives_wf, [ ('outputnode.out_aseg', 'inputnode.t1w_fs_aseg'), - ('outputnode.out_aparc', 'inputnode.t1w_fs_aparc'), - ]), + ('outputnode.out_aparc', 'inputnode.t1w_fs_aparc')]), (outputnode, anat_derivatives_wf, [ ('t1w2fsnative_xfm', 'inputnode.t1w2fsnative_xfm'), ('fsnative2t1w_xfm', 'inputnode.fsnative2t1w_xfm'), @@ -639,7 +681,6 @@ def _check_img(img): ("outputnode.anat_ribbon", "inputnode.anat_ribbon"), ]), ]) - # fmt:on if cifti_output: morph_grayords_wf = init_morph_grayords_wf(grayord_density=cifti_output) @@ -730,28 +771,19 @@ def init_anat_template_wf( name="outputnode", ) - # 0. Denoise and reorient T1w image(s) to RAS and resample to common voxel space + # 0. Reorient T1w image(s) to RAS and resample to common voxel space anat_ref_dimensions = pe.Node(TemplateDimensions(), name="anat_ref_dimensions") - denoise = pe.MapNode( - DenoiseImage(noise_model="Rician", num_threads=omp_nthreads), - iterfield="input_image", - name="denoise", - ) anat_conform = pe.MapNode(Conform(), iterfield="in_file", name="anat_conform") # fmt:off workflow.connect([ (inputnode, anat_ref_dimensions, [('anat_files', 't1w_list')]), - (anat_ref_dimensions, denoise, [('t1w_valid_list', 'input_image')]), (anat_ref_dimensions, anat_conform, [ + ('t1w_valid_list', 'in_file'), ('target_zooms', 'target_zooms'), - ('target_shape', 'target_shape'), - ]), - (denoise, anat_conform, [('output_image', 'in_file')]), - (anat_ref_dimensions, outputnode, [ - ('out_report', 'out_report'), - ('t1w_valid_list', 'anat_valid_list'), - ]), + ('target_shape', 'target_shape')]), + (anat_ref_dimensions, outputnode, [('out_report', 'out_report'), + ('t1w_valid_list', 'anat_valid_list')]), ]) # fmt:on diff --git a/smriprep/workflows/base.py b/smriprep/workflows/base.py index 6e6f921be5..d3720398a7 100644 --- a/smriprep/workflows/base.py +++ b/smriprep/workflows/base.py @@ -50,6 +50,7 @@ def init_smriprep_wf( hires, layout, longitudinal, + fastsurfer, low_mem, omp_nthreads, output_dir, @@ -87,6 +88,7 @@ def init_smriprep_wf( hires=True, layout=BIDSLayout('.'), longitudinal=False, + fastsurfer=False, low_mem=False, omp_nthreads=1, output_dir='.', @@ -117,6 +119,8 @@ def init_smriprep_wf( longitudinal : :obj:`bool` Treat multiple sessions as longitudinal (may increase runtime) See sub-workflows for specific differences + fastsurfer : :obj:`bool` + Enable FastSurfer segmentation and surface reconstruction low_mem : :obj:`bool` Write uncompressed .nii files in some cases to reduce memory usage omp_nthreads : :obj:`int` @@ -162,6 +166,20 @@ def init_smriprep_wf( if fs_subjects_dir is not None: fsdir.inputs.subjects_dir = str(fs_subjects_dir.absolute()) + if fastsurfer: + freesurfer = False + fastsurfdir = pe.Node( + BIDSFreeSurferDir( + derivatives=output_dir, + freesurfer_home=os.getenv("FREESURFER_HOME"), + spaces=spaces.get_fs_spaces(), + ), + name="fastsurfdir_run_%s" % run_uuid.replace("-", "_"), + run_without_submitting=True + ) + if fs_subjects_dir is not None: + fastsurfdir.inputs.subjects_dir = str(fs_subjects_dir.absolute()) + for subject_id in subject_list: single_subject_wf = init_single_subject_wf( debug=debug, @@ -170,6 +188,7 @@ def init_smriprep_wf( hires=hires, layout=layout, longitudinal=longitudinal, + fastsurfer=fastsurfer, low_mem=low_mem, name="single_subject_%s_wf" % subject_id, omp_nthreads=omp_nthreads, @@ -181,7 +200,7 @@ def init_smriprep_wf( subject_id=subject_id, bids_filters=bids_filters, ) - + single_subject_wf.config["execution"]["crashdump_dir"] = os.path.join( output_dir, "smriprep", "sub-" + subject_id, "log", run_uuid ) @@ -191,6 +210,10 @@ def init_smriprep_wf( smriprep_wf.connect( fsdir, "subjects_dir", single_subject_wf, "inputnode.subjects_dir" ) + elif fastsurfer: + smriprep_wf.connect( + fastsurfdir, "subjects_dir", single_subject_wf, "inputnode.subjects_dir" + ) else: smriprep_wf.add_nodes([single_subject_wf]) @@ -205,6 +228,7 @@ def init_single_subject_wf( hires, layout, longitudinal, + fastsurfer, low_mem, name, omp_nthreads, @@ -244,6 +268,7 @@ def init_single_subject_wf( hires=True, layout=BIDSLayout('.'), longitudinal=False, + fastsurfer=False, low_mem=False, name='single_subject_wf', omp_nthreads=1, @@ -271,6 +296,8 @@ def init_single_subject_wf( longitudinal : :obj:`bool` Treat multiple sessions as longitudinal (may increase runtime) See sub-workflows for specific differences + fastsurfer : :obj:`bool` + Enable FastSurfer segmentation and surface reconstruction low_mem : :obj:`bool` Write uncompressed .nii files in some cases to reduce memory usage name : :obj:`str` @@ -351,6 +378,9 @@ def init_single_subject_wf( Path(output_dir) / "smriprep", subject_id, std_spaces, freesurfer ) + if fastsurfer: + freesurfer = False + inputnode = pe.Node( niu.IdentityInterface(fields=["subjects_dir"]), name="inputnode" ) @@ -405,6 +435,7 @@ def init_single_subject_wf( freesurfer=freesurfer, hires=hires, longitudinal=longitudinal, + fastsurfer=fastsurfer, name="anat_preproc_wf", t1w=subject_data["t1w"], t2w=subject_data["t2w"], diff --git a/smriprep/workflows/outputs.py b/smriprep/workflows/outputs.py index 80d907d648..e5f4d6357b 100644 --- a/smriprep/workflows/outputs.py +++ b/smriprep/workflows/outputs.py @@ -31,7 +31,7 @@ BIDS_TISSUE_ORDER = ("GM", "WM", "CSF") -def init_anat_reports_wf(*, freesurfer, output_dir, name="anat_reports_wf"): +def init_anat_reports_wf(*, freesurfer, fastsurfer, output_dir, name="anat_reports_wf"): """ Set up a battery of datasinks to store reports in the right location. @@ -39,6 +39,8 @@ def init_anat_reports_wf(*, freesurfer, output_dir, name="anat_reports_wf"): ---------- freesurfer : :obj:`bool` FreeSurfer was enabled + fastsurfer : :obj:`bool` + FastSurfer was enabled output_dir : :obj:`str` Directory in which to save derivatives name : :obj:`str` @@ -171,7 +173,7 @@ def init_anat_reports_wf(*, freesurfer, output_dir, name="anat_reports_wf"): ]) # fmt:on - if freesurfer: + if freesurfer or fastsurfer: from ..interfaces.reports import FSSurfaceReport recon_report = pe.Node(FSSurfaceReport(), name="recon_report") @@ -200,6 +202,7 @@ def init_anat_derivatives_wf( *, bids_root, freesurfer, + fastsurfer, num_t1w, t2w, output_dir, @@ -217,6 +220,8 @@ def init_anat_derivatives_wf( Root path of BIDS dataset freesurfer : :obj:`bool` FreeSurfer was enabled + fastsurfer : :obj:`bool` + FastSurfer was enabled num_t1w : :obj:`int` Number of T1w images output_dir : :obj:`str` @@ -574,7 +579,7 @@ def init_anat_derivatives_wf( ) # fmt:on - if not freesurfer: + if (not freesurfer) and (not fastsurfer): return workflow # T2w coregistration requires FreeSurfer surfaces, so only try to save if freesurfer diff --git a/smriprep/workflows/surfaces.py b/smriprep/workflows/surfaces.py index 2027d3cf97..cd2dbf76f1 100644 --- a/smriprep/workflows/surfaces.py +++ b/smriprep/workflows/surfaces.py @@ -24,7 +24,8 @@ Surface preprocessing workflows. **sMRIPrep** uses FreeSurfer to reconstruct surfaces from T1w/T2w -structural images. +structural images or FastSurfer for a machine learning-accelerated +alternative using only T1w structural images. """ import typing as ty @@ -39,6 +40,7 @@ ) from ..interfaces.freesurfer import ReconAll, MakeMidthickness +from ..interfaces import fastsurfer as fastsurf from niworkflows.engine.workflows import LiterateWorkflow as Workflow from niworkflows.interfaces.freesurfer import ( @@ -48,6 +50,283 @@ PatchedRobustRegister as RobustRegister, RefineBrainMask, ) + +from nipype import logging +# import os + +LOGGER = logging.getLogger("nipype.workflow") + + +def init_fastsurf_recon_wf(*, omp_nthreads, hires, name="fastsurf_recon_wf"): + r""" + Reconstruct anatomical surfaces using FastSurfer VINN and ``recon_surf``, + an alternative to FreeSurfer's ``recon-all``. + + First, FastSurfer VINN creates a segmentation using the T1w structural image. + This is followed by FastSurfer ``recon_surf`` processing of the surface. + + For example, a subject with only one session with a T1w image + would be processed by the following command:: + + $ /opt/FastSurfer/run_fastsurfer.sh \ + --t1 /sub-/anat/sub-_T1w.nii.gz \ + --sid sub- --sd /fastsurfer + + Similar to the Freesurfer workflow second phase, we then + import an externally computed skull-stripping mask. + This workflow refines the external brainmask using the internal mask + implicit the FastSurfer's ``aseg.mgz`` segmentation, + to reconcile ANTs' and FreeSurfer's brain masks. + + First, the ``aseg.mgz`` mask from FastSurfer is refined in two + steps, using binary morphological operations: + + 1. With a binary closing operation the sulci are included + into the mask. This results in a smoother brain mask + that does not exclude deep, wide sulci. + + 2. Fill any holes (typically, there could be a hole next to + the pineal gland and the corpora quadrigemina if the great + cerebral brain is segmented out). + + Second, the brain mask is grown, including pixels that have a high likelihood + to the GM tissue distribution: + + 3. Dilate and substract the brain mask, defining the region to search for candidate + pixels that likely belong to cortical GM. + + 4. Pixels found in the search region that are labeled as GM by ANTs + (during ``antsBrainExtraction.sh``) are directly added to the new mask. + + 5. Otherwise, estimate GM tissue parameters locally in patches of ``ww`` size, + and test the likelihood of the pixel to belong in the GM distribution. + + This procedure is inspired on mindboggle's solution to the problem: + https://github.com/nipy/mindboggle/blob/7f91faaa7664d820fe12ccc52ebaf21d679795e2/mindboggle/guts/segment.py#L1660 + + Memory annotations for FastSurfer are based off `their documentation + `_. + They recommend 8GB CPU RAM and 8GB NVIDIA GPU RAM to run a single batch. + + Workflow Graph + .. workflow:: + :graph2use: orig + :simple_form: yes + + from smriprep.workflows.surfaces import init_fastsurf_recon_wf + wf = init_fastsurf_recon_wf(omp_nthreads=1) + + Required arguments + ================== + sd + Output directory + sid + Subject ID for directory inside ``sd`` to be created + t1 + T1 full head input (not bias corrected, global path). + The 'network was trained with conformed images + (UCHAR, 256x256x256, 1 mm voxels and standard slice orientation). + These specifications are checked in the ``eval.py`` script and the image + is automatically conformed if it does not comply. + + Optional arguments + ================== + + Network specific arguments + -------------------------- + seg + Global path with filename of segmentation + (where and under which name to store it). + Default location + ``$SUBJECTS_DIR/$sid/mri/aparc.DKTatlas+aseg.deep.mgz`` + weights_sag + Pretrained weights of sagittal network. + Default + ``../checkpoints/Sagittal_Weights_FastSurferVINN/ckpts/Epoch_30_training_state.pkl`` + weights_ax + Pretrained weights of axial network. + Default + ``../checkpoints/Axial_Weights_FastSurferVINN/ckpts/Epoch_30_training_state.pkl`` + weights_cor + Pretrained weights of coronal network. + Default + ``../checkpoints/Coronal_Weights_FastSurferVINN/ckpts/Epoch_30_training_state.pkl`` + seg_log + Name and location for the log-file for the segmentation (FastSurferVINN). + Default '$SUBJECTS_DIR/$sid/scripts/deep-seg.log' + clean_seg + Flag to clean up FastSurferVINN segmentation + run_viewagg_on + Define where the view aggregation should be run on. + By default, the program checks if you have enough memory to run + the view aggregation on the gpu. The total memory is considered for this decision. + If this fails, or you actively overwrote the check with setting + ``run_viewagg_on cpu``, view agg is run on the cpu. + Equivalently, if you define ``--run_viewagg_on gpu``, view agg will be run on the gpu + (no memory check will be done). + no_cuda + Flag to disable CUDA usage in FastSurferVINN (no GPU usage, inference on CPU) + batch + Batch size for inference. Default 16. Lower this to reduce memory requirement + order + Order of interpolation for mri_convert T1 before segmentation + ``(0=nearest, 1=linear(default), 2=quadratic, 3=cubic)`` + + Surface pipeline arguments + -------------------------- + fstess + Use ``mri_tesselate`` instead of marching cube (default) for surface creation + fsqsphere + Use FreeSurfer default instead of + novel spectral spherical projection for qsphere + fsaparc + Use FS aparc segmentations in addition to DL prediction + (slower in this case and usually the mapped ones from the DL prediction are fine) + no_surfreg + Skip creating Surface-Atlas ``sphere.reg`` registration with FreeSurfer + (for cross-subject correspondence or other mappings) + parallel + Run both hemispheres in parallel + threads + Set openMP and ITK threads + + Other + ---- + py + which python version to use. + Default ``python3.6`` + seg_only + only run FastSurferVINN + (generate segmentation, do not run the surface pipeline) + surf_only + only run the surface pipeline ``recon_surf``. + The segmentation created by FastSurferVINN must already exist in this case. + + + """ + workflow = Workflow(name=name) + workflow.__desc__ = """\ +Brain surfaces were reconstructed using `recon-surf` [FastSurfer {fastsurfer_version}, +@fastsurfer], and the brain mask estimated +previously was refined with a custom variation of the method to reconcile +ANTs-derived and FastSurfer-derived segmentations of the cortical +gray-matter of Mindboggle [RRID:SCR_002438, @mindboggle]. +""".format( + fastsurfer_version="2.0.4" + ) + + inputnode = pe.Node( + niu.IdentityInterface( + fields=[ + "subjects_dir", + "subject_id", + "t1w", + "skullstripped_t1", + "corrected_t1", + "ants_segs", + ] + ), + name="inputnode", + ) + outputnode = pe.Node( + niu.IdentityInterface( + fields=[ + "subjects_dir", + "subject_id", + "t1w2fsnative_xfm", + "fsnative2t1w_xfm", + "surfaces", + "out_brainmask", + "out_aseg", + "out_aparc", + "morphometrics", + ] + ), + name="outputnode", + ) + + # outputnode = pe.Node(niu.IdentityInterface(["surfaces", "morphometrics"]), name="outputnode") + + fsnative2t1w_xfm = pe.Node( + RobustRegister(auto_sens=True, est_int_scale=True), name="fsnative2t1w_xfm") + t1w2fsnative_xfm = pe.Node( + LTAConvert(out_lta=True, invert=True), name="t1w2fsnative_xfm") + + # inputnode.in_file = inputnode.t1w + gifti_surface_wf = init_gifti_surface_wf(fastsurfer=True) + aseg_to_native_wf = init_segs_to_native_wf(fastsurfer=True) + aparc_to_native_wf = init_segs_to_native_wf(fastsurfer=True, segmentation="aparc_aseg") + refine = pe.Node(RefineBrainMask(), name="refine") + + recon_conf = pe.Node( + fastsurf.FastSurferSource(), name="recon_conf") + + fastsurf_recon = pe.Node( + fastsurf.FastSurfer(), + name="fastsurf_recon", + n_procs=omp_nthreads, + mem_gb=12, + ) + fastsurf_recon.interface._can_resume = False + fastsurf_recon.interface._always_run = True + + skull_strip_extern = pe.Node(FSInjectBrainExtracted(), name="skull_strip_extern") + + fsnative2t1w_xfm = pe.Node( + RobustRegister(auto_sens=True, est_int_scale=True), name="fsnative2t1w_xfm" + ) + t1w2fsnative_xfm = pe.Node( + LTAConvert(out_lta=True, invert=True), name="t1w2fsnative_xfm" + ) + + # fmt:off + workflow.connect([ + (inputnode, fastsurf_recon, [('subjects_dir', 'subjects_dir'), + ('subject_id', 'subject_id'), # ]), + ('t1w', 'T1_files')]), + # (inputnode, skull_strip_extern, [('subjects_dir', 'subjects_dir'), + # ('subject_id', 'subject_id')]), + # (inputnode, skull_strip_extern, [('skullstripped_t1', 'in_brain')]), + (fastsurf_recon, gifti_surface_wf, [ + ('subjects_dir', 'inputnode.subjects_dir'), + ('subject_id', 'inputnode.subject_id')]), + # (inputnode, skull_strip_extern, [('skullstripped_t1', 'in_brain')]), + # Construct transform from FreeSurfer conformed image to sMRIPrep + # reoriented image + (inputnode, fsnative2t1w_xfm, [('t1w', 'target_file')]), + (fastsurf_recon, fsnative2t1w_xfm, [('T1', 'source_file')]), + (fsnative2t1w_xfm, gifti_surface_wf, [ + ('out_reg_file', 'inputnode.fsnative2t1w_xfm')]), + (fsnative2t1w_xfm, t1w2fsnative_xfm, [('out_reg_file', 'in_lta')]), + # Refine ANTs mask, deriving new mask from FS' aseg + (inputnode, refine, [('corrected_t1', 'in_anat'), + ('ants_segs', 'in_ants')]), + (inputnode, aseg_to_native_wf, [('corrected_t1', 'inputnode.in_file')]), + (fastsurf_recon, aseg_to_native_wf, [ + ('subjects_dir', 'inputnode.subjects_dir'), + ('subject_id', 'inputnode.subject_id')]), + (fsnative2t1w_xfm, aseg_to_native_wf, [ + ('out_reg_file', 'inputnode.fsnative2t1w_xfm')]), + (inputnode, aparc_to_native_wf, [('corrected_t1', 'inputnode.in_file')]), + (fsnative2t1w_xfm, aparc_to_native_wf, [ + ('out_reg_file', 'inputnode.fsnative2t1w_xfm')]), + (aseg_to_native_wf, refine, [('outputnode.out_file', 'in_aseg')]), + + # Output + (inputnode, outputnode, [('subjects_dir', 'subjects_dir'), + ('subject_id', 'subject_id')]), + (gifti_surface_wf, outputnode, [('outputnode.surfaces', 'surfaces'), + ('outputnode.morphometrics', 'morphometrics')]), + (t1w2fsnative_xfm, outputnode, [('out_lta', 't1w2fsnative_xfm')]), + (fsnative2t1w_xfm, outputnode, [('out_reg_file', 'fsnative2t1w_xfm')]), + (refine, outputnode, [('out_file', 'out_brainmask')]), + (aseg_to_native_wf, outputnode, [('outputnode.out_file', 'out_aseg')]), + (aparc_to_native_wf, outputnode, [('outputnode.out_file', 'out_aparc')]), + ]) + # fmt:on + + return workflow + from ..interfaces.workbench import CreateSignedDistanceVolume @@ -237,10 +516,9 @@ def init_surface_recon_wf(*, omp_nthreads, hires, name="surface_recon_wf"): ) autorecon_resume_wf = init_autorecon_resume_wf(omp_nthreads=omp_nthreads) - gifti_surface_wf = init_gifti_surface_wf() - - aseg_to_native_wf = init_segs_to_native_wf() - aparc_to_native_wf = init_segs_to_native_wf(segmentation="aparc_aseg") + gifti_surface_wf = init_gifti_surface_wf(fastsurfer=False) + aseg_to_native_wf = init_segs_to_native_wf(fastsurfer=False) + aparc_to_native_wf = init_segs_to_native_wf(fastsurfer=False, segmentation="aparc_aseg") refine = pe.Node(RefineBrainMask(), name="refine") # fmt:off @@ -486,7 +764,7 @@ def _dedup(in_list): return workflow -def init_gifti_surface_wf(*, name="gifti_surface_wf"): +def init_gifti_surface_wf(*, fastsurfer, name="gifti_surface_wf"): r""" Prepare GIFTI surfaces from a FreeSurfer subjects directory. @@ -514,6 +792,8 @@ def init_gifti_surface_wf(*, name="gifti_surface_wf"): FreeSurfer subject ID fsnative2t1w_xfm LTA formatted affine transform file (inverse) + fastsurfer + Boolean to indicate FastSurfer surface processing Outputs ------- @@ -595,7 +875,7 @@ def init_gifti_surface_wf(*, name="gifti_surface_wf"): return workflow -def init_segs_to_native_wf(*, name="segs_to_native", segmentation="aseg"): +def init_segs_to_native_wf(*, fastsurfer, name="segs_to_native", segmentation="aseg"): """ Get a segmentation from FreeSurfer conformed space into native T1w space. @@ -622,6 +902,8 @@ def init_segs_to_native_wf(*, name="segs_to_native", segmentation="aseg"): FreeSurfer subject ID fsnative2t1w_xfm LTA-style affine matrix translating from FreeSurfer-conformed subject space to T1w + fastsurfer + Boolean to indicate FastSurfer processing Outputs ------- @@ -639,6 +921,7 @@ def init_segs_to_native_wf(*, name="segs_to_native", segmentation="aseg"): outputnode = pe.Node(niu.IdentityInterface(["out_file"]), name="outputnode") # Extract the aseg and aparc+aseg outputs fssource = pe.Node(nio.FreeSurferSource(), name="fs_datasource") + # Resample from T1.mgz to T1w.nii.gz, applying any offset in fsnative2t1w_xfm, # and convert to NIfTI while we're at it resample = pe.Node( @@ -838,7 +1121,7 @@ def init_anat_ribbon_wf(name="anat_ribbon_wf"): def init_morph_grayords_wf( grayord_density: ty.Literal['91k', '170k'], - name: str = "morph_grayords_wf", + name: str = "bold_grayords_wf", ): """ Sample Grayordinates files onto the fsLR atlas. @@ -858,7 +1141,7 @@ def init_morph_grayords_wf( grayord_density : :obj:`str` Either `91k` or `170k`, representing the total of vertices or *grayordinates*. name : :obj:`str` - Unique name for the subworkflow (default: ``"morph_grayords_wf"``) + Unique name for the subworkflow (default: ``"bold_grayords_wf"``) Inputs ------