diff --git a/.circleci/config.yml b/.circleci/config.yml index 702cbdd901..be7a222e18 100644 --- a/.circleci/config.yml +++ b/.circleci/config.yml @@ -12,6 +12,9 @@ _machine_defaults: &machine_defaults _python_defaults: &python_defaults docker: - image: cimg/python:3.10.9 + auth: + username: $DOCKER_USER + password: $DOCKER_PAT working_directory: /tmp/src/smriprep _docker_auth: &docker_auth @@ -43,6 +46,9 @@ _pull_from_registry: &pull_from_registry docs_deploy: &docs docker: - image: node:8.10.0 + auth: + username: $DOCKER_USER + password: $DOCKER_PAT working_directory: /tmp/gh-pages steps: - run: @@ -282,8 +288,10 @@ jobs: paths: - /home/circleci/.local - test_wrapper: + test: <<: *machine_defaults + environment: + - FS_LICENSE: /tmp/fslicense/license.txt steps: - attach_workspace: at: /tmp @@ -323,6 +331,24 @@ jobs: which smriprep-docker smriprep-docker -i nipreps/smriprep:latest --help smriprep-docker -i nipreps/smriprep:latest --version + - restore_cache: + keys: + - testdata-v2-{{ .Branch }}-{{ epoch }} + - testdata-v2-{{ .Branch }} + - testdata-v2- + - restore_cache: + keys: + - templateflow-v1-{{ .Branch }}-{{ epoch }} + - templateflow-v1-{{ .Branch }} + - templateflow-v1- + - run: + name: Run Pytest + no_output_timeout: 2h + command: | + bash /tmp/src/smriprep/.circleci/pytest.sh + - codecov/upload: + file: /tmp/data/pytest_cov.xml + flags: pytest test_deploy_pypi: <<: *python_defaults @@ -414,9 +440,10 @@ jobs: - run: *pull_from_registry - restore_cache: keys: - - ds005-anat-v9-{{ .Branch }}-{{ epoch }} - - ds005-anat-v9-{{ .Branch }} - - ds005-anat-v9-master + - ds005-anat-v0-{{ .Branch }}-{{ epoch }} + - ds005-anat-v0-{{ .Branch }} + - ds005-anat-v0-master + - ds005-anat-v0-next - restore_cache: keys: - testdata-v2-{{ .Branch }}-{{ epoch }} @@ -470,7 +497,7 @@ jobs: rm -rf /tmp/ds005/work/reportlets rm -rf /tmp/ds005/work/smriprep_wf/fsdir_run_*/ - save_cache: - key: ds005-anat-v9-{{ .Branch }}-{{ epoch }} + key: ds005-anat-v0-{{ .Branch }}-{{ epoch }} paths: - /tmp/ds005/work @@ -772,8 +799,6 @@ workflows: only: /.*/ - get_data: - context: - - fs-license filters: branches: ignore: @@ -793,10 +818,11 @@ workflows: tags: only: /.*/ - - test_wrapper: + - test: context: - nipreps-common requires: + - get_data - build filters: branches: @@ -852,7 +878,7 @@ workflows: - build_docs - ds005 - ds054 - - test_wrapper + - test filters: branches: only: /master/ @@ -865,7 +891,7 @@ workflows: - ds005 - ds054 - test_deploy_pypi - - test_wrapper + - test filters: branches: ignore: /.*/ diff --git a/.circleci/ds005_outputs.txt b/.circleci/ds005_outputs.txt index 41d92ff92a..779c44abc6 100644 --- a/.circleci/ds005_outputs.txt +++ b/.circleci/ds005_outputs.txt @@ -16,6 +16,7 @@ smriprep/sub-01/anat/sub-01_desc-brain_mask.json smriprep/sub-01/anat/sub-01_desc-brain_mask.nii.gz smriprep/sub-01/anat/sub-01_desc-preproc_T1w.json smriprep/sub-01/anat/sub-01_desc-preproc_T1w.nii.gz +smriprep/sub-01/anat/sub-01_desc-ribbon_mask.json smriprep/sub-01/anat/sub-01_desc-ribbon_mask.nii.gz smriprep/sub-01/anat/sub-01_dseg.nii.gz smriprep/sub-01/anat/sub-01_from-fsnative_to-T1w_mode-image_xfm.txt @@ -25,7 +26,9 @@ smriprep/sub-01/anat/sub-01_hemi-L_desc-reg_sphere.surf.gii smriprep/sub-01/anat/sub-01_hemi-L_inflated.surf.gii smriprep/sub-01/anat/sub-01_hemi-L_midthickness.surf.gii smriprep/sub-01/anat/sub-01_hemi-L_pial.surf.gii +smriprep/sub-01/anat/sub-01_hemi-L_space-fsLR_desc-msmsulc_sphere.surf.gii smriprep/sub-01/anat/sub-01_hemi-L_space-fsLR_desc-reg_sphere.surf.gii +smriprep/sub-01/anat/sub-01_hemi-L_sphere.surf.gii smriprep/sub-01/anat/sub-01_hemi-L_sulc.shape.gii smriprep/sub-01/anat/sub-01_hemi-L_thickness.shape.gii smriprep/sub-01/anat/sub-01_hemi-L_white.surf.gii @@ -34,12 +37,20 @@ smriprep/sub-01/anat/sub-01_hemi-R_desc-reg_sphere.surf.gii smriprep/sub-01/anat/sub-01_hemi-R_inflated.surf.gii smriprep/sub-01/anat/sub-01_hemi-R_midthickness.surf.gii smriprep/sub-01/anat/sub-01_hemi-R_pial.surf.gii +smriprep/sub-01/anat/sub-01_hemi-R_space-fsLR_desc-msmsulc_sphere.surf.gii smriprep/sub-01/anat/sub-01_hemi-R_space-fsLR_desc-reg_sphere.surf.gii +smriprep/sub-01/anat/sub-01_hemi-R_sphere.surf.gii smriprep/sub-01/anat/sub-01_hemi-R_sulc.shape.gii smriprep/sub-01/anat/sub-01_hemi-R_thickness.shape.gii smriprep/sub-01/anat/sub-01_hemi-R_white.surf.gii smriprep/sub-01/anat/sub-01_label-CSF_probseg.nii.gz smriprep/sub-01/anat/sub-01_label-GM_probseg.nii.gz smriprep/sub-01/anat/sub-01_label-WM_probseg.nii.gz +smriprep/sub-01/anat/sub-01_space-fsLR_den-91k_curv.dscalar.nii +smriprep/sub-01/anat/sub-01_space-fsLR_den-91k_curv.json +smriprep/sub-01/anat/sub-01_space-fsLR_den-91k_sulc.dscalar.nii +smriprep/sub-01/anat/sub-01_space-fsLR_den-91k_sulc.json +smriprep/sub-01/anat/sub-01_space-fsLR_den-91k_thickness.dscalar.nii +smriprep/sub-01/anat/sub-01_space-fsLR_den-91k_thickness.json smriprep/sub-01.html /tmp/ds005/derivatives diff --git a/.circleci/ds005_run.sh b/.circleci/ds005_run.sh index 864ea87f52..a90682b518 100644 --- a/.circleci/ds005_run.sh +++ b/.circleci/ds005_run.sh @@ -19,4 +19,5 @@ docker run -it -e FMRIPREP_DEV=1 -u $(id -u) \ --ncpus 2 --omp-nthreads 2 -vv \ --fs-license-file /tmp/fslicense/license.txt \ --fs-subjects-dir /tmp/ds005/freesurfer \ + --cifti-output \ ${@:1} diff --git a/.circleci/nipype.cfg b/.circleci/nipype.cfg index 10451477e5..ed07e98c3c 100644 --- a/.circleci/nipype.cfg +++ b/.circleci/nipype.cfg @@ -1,4 +1,4 @@ [execution] stop_on_first_crash = true poll_sleep_duration = 0.01 -hash_method = content \ No newline at end of file +hash_method = content diff --git a/.circleci/pytest.sh b/.circleci/pytest.sh new file mode 100644 index 0000000000..ff1114a540 --- /dev/null +++ b/.circleci/pytest.sh @@ -0,0 +1,13 @@ +#!/bin/bash +docker run --rm -it \ + -v /tmp/data:/tmp/data:rw \ + -v /tmp/fslicense:/tmp/fslicense:ro \ + -v /tmp/templateflow:/home/smriprep/.cache/templateflow \ + -v /tmp/src/smriprep/docker/multiproc.coveragerc:/tmp/multiproc.coveragerc:ro \ + -v /tmp/src/smriprep/.circleci/nipype.cfg:/home/smriprep/.nipype/nipype.cfg \ + -e FS_LICENSE=/tmp/fslicense/license.txt \ + --entrypoint=pytest \ + nipreps/smriprep:latest \ + -v --doctest-modules --pyargs smriprep \ + --cov smriprep --cov-report=xml:/tmp/data/pytest_cov.xml \ + ${@:1} diff --git a/.github/workflows/pythonpackage.yml b/.github/workflows/pythonpackage.yml index b1904f5cc3..353477f5d1 100644 --- a/.github/workflows/pythonpackage.yml +++ b/.github/workflows/pythonpackage.yml @@ -5,7 +5,7 @@ on: branches: [ '*' ] tags: [ '*' ] pull_request: - branches: [ master, 'maint/*' ] + branches: [ master, 'maint/*', 'next' ] schedule: - cron: '0 0 * * *' @@ -65,7 +65,7 @@ jobs: runs-on: ubuntu-latest strategy: matrix: - python-version: [3.8, 3.9, "3.10", "3.11"] + python-version: ["3.10", "3.11"] install: [repo] include: - python-version: "3.11" diff --git a/.gitignore b/.gitignore index 083f736b3a..44286ed081 100644 --- a/.gitignore +++ b/.gitignore @@ -1,4 +1,4 @@ -# Auto-generated by .maint/make_gitignore on Fri May 19 03:22:33 PM EDT 2023 +# Auto-generated by .maint/make_gitignore on Wed Nov 8 04:05:04 PM EST 2023 # setuptools_scm _version.py diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index b803d4cc23..bed9fef756 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -1,11 +1,11 @@ -exclude: ".*/data/.*" repos: - repo: https://github.com/pre-commit/pre-commit-hooks rev: v4.4.0 hooks: - id: trailing-whitespace - exclude: .gitignore + exclude: '.gitignore|.*\.gii$' - id: end-of-file-fixer + exclude: 'data/itkIdentityTransform.txt|\.svg$' - id: check-yaml - id: check-added-large-files - id: check-toml diff --git a/Dockerfile b/Dockerfile index a14ec26dfc..47d63718f6 100644 --- a/Dockerfile +++ b/Dockerfile @@ -56,27 +56,6 @@ COPY docker/files/freesurfer7.3.2-exclude.txt /usr/local/etc/freesurfer7.3.2-exc 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 -# AFNI -FROM downloader as afni -# Bump the date to current to update AFNI -RUN echo "2023.04.04" -RUN mkdir -p /opt/afni-latest \ - && curl -fsSL --retry 5 https://afni.nimh.nih.gov/pub/dist/tgz/linux_openmp_64.tgz \ - | tar -xz -C /opt/afni-latest --strip-components 1 \ - --exclude "linux_openmp_64/*.gz" \ - --exclude "linux_openmp_64/funstuff" \ - --exclude "linux_openmp_64/shiny" \ - --exclude "linux_openmp_64/afnipy" \ - --exclude "linux_openmp_64/lib/RetroTS" \ - --exclude "linux_openmp_64/lib_RetroTS" \ - --exclude "linux_openmp_64/meica.libs" \ - # Keep only what we use - && find /opt/afni-latest -type f -not \( \ - -name "3dTshift" -or \ - -name "3dUnifize" -or \ - -name "3dAutomask" -or \ - -name "3dvolreg" \) -delete - # Connectome Workbench 1.5.0 FROM downloader as workbench RUN mkdir /opt/workbench && \ @@ -115,7 +94,9 @@ ENV DEBIAN_FRONTEND="noninteractive" \ LANG="en_US.UTF-8" \ LC_ALL="en_US.UTF-8" -# Some baseline tools; bc is needed for FreeSurfer, so don't drop it +# Some baseline tools +# bc is needed for FreeSurfer +# libglu1-mesa is needed for Connectome Workbench RUN apt-get update && \ apt-get install -y --no-install-recommends \ bc \ @@ -123,48 +104,14 @@ RUN apt-get update && \ curl \ git \ gnupg \ + libglu1-mesa \ lsb-release \ netbase \ xvfb && \ apt-get clean && rm -rf /var/lib/apt/lists/* /tmp/* /var/tmp/* -# Configure PPAs for libpng12 and libxp6 -RUN GNUPGHOME=/tmp gpg --keyserver hkps://keyserver.ubuntu.com --no-default-keyring --keyring /usr/share/keyrings/linuxuprising.gpg --recv 0xEA8CACC073C3DB2A \ - && GNUPGHOME=/tmp gpg --keyserver hkps://keyserver.ubuntu.com --no-default-keyring --keyring /usr/share/keyrings/zeehio.gpg --recv 0xA1301338A3A48C4A \ - && echo "deb [signed-by=/usr/share/keyrings/linuxuprising.gpg] https://ppa.launchpadcontent.net/linuxuprising/libpng12/ubuntu jammy main" > /etc/apt/sources.list.d/linuxuprising.list \ - && echo "deb [signed-by=/usr/share/keyrings/zeehio.gpg] https://ppa.launchpadcontent.net/zeehio/libxp/ubuntu jammy main" > /etc/apt/sources.list.d/zeehio.list - -# Dependencies for AFNI; requires a discontinued multiarch-support package from bionic (18.04) -RUN apt-get update -qq \ - && apt-get install -y -q --no-install-recommends \ - ed \ - gsl-bin \ - libglib2.0-0 \ - libglu1-mesa-dev \ - libglw1-mesa \ - libgomp1 \ - libjpeg62 \ - libpng12-0 \ - libxm4 \ - libxp6 \ - netpbm \ - tcsh \ - xfonts-base \ - xvfb \ - && 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 \ - && dpkg -i /tmp/multiarch.deb \ - && rm /tmp/multiarch.deb \ - && apt-get install -f \ - && apt-get clean && rm -rf /var/lib/apt/lists/* /tmp/* /var/tmp/* \ - && gsl2_path="$(find / -name 'libgsl.so.19' || printf '')" \ - && if [ -n "$gsl2_path" ]; then \ - ln -sfv "$gsl2_path" "$(dirname $gsl2_path)/libgsl.so.0"; \ - fi \ - && ldconfig - # Install files from stages COPY --from=freesurfer /opt/freesurfer /opt/freesurfer -COPY --from=afni /opt/afni-latest /opt/afni-latest COPY --from=workbench /opt/workbench /opt/workbench # Simulate SetUpFreeSurfer.sh @@ -184,11 +131,6 @@ ENV PERL5LIB="$MINC_LIB_DIR/perl5/5.8.5" \ MNI_PERL5LIB="$MINC_LIB_DIR/perl5/5.8.5" \ PATH="$FREESURFER_HOME/bin:$FREESURFER_HOME/tktools:$MINC_BIN_DIR:$PATH" -# AFNI config -ENV PATH="/opt/afni-latest:$PATH" \ - AFNI_IMSAVE_WARNINGS="NO" \ - AFNI_PLUGINPATH="/opt/afni-latest" - # Workbench config ENV PATH="/opt/workbench/bin_linux64:$PATH" \ LD_LIBRARY_PATH="/opt/workbench/lib_linux64:$LD_LIBRARY_PATH" diff --git a/Makefile b/Makefile index 501071c083..89e6439f24 100644 --- a/Makefile +++ b/Makefile @@ -15,4 +15,3 @@ docker: --build-arg BUILD_DATE=`date -u +"%Y-%m-%dT%H:%M:%SZ"` \ --build-arg VCS_REF=`git rev-parse --short HEAD` \ --build-arg VERSION=`python setup.py --version` . - diff --git a/README.rst b/README.rst index c2aac2eb0b..0d7f2f12f4 100644 --- a/README.rst +++ b/README.rst @@ -7,7 +7,7 @@ sMRIPrep: Structural MRI PREProcessing pipeline .. image:: https://circleci.com/gh/nipreps/smriprep/tree/master.svg?style=shield :target: https://circleci.com/gh/nipreps/smriprep/tree/master - + .. image:: https://codecov.io/gh/nipreps/smriprep/branch/master/graph/badge.svg :target: https://codecov.io/gh/nipreps/smriprep :alt: Coverage report @@ -15,7 +15,7 @@ sMRIPrep: Structural MRI PREProcessing pipeline .. image:: https://img.shields.io/pypi/v/smriprep.svg :target: https://pypi.python.org/pypi/smriprep/ :alt: Latest Version - + .. image:: https://img.shields.io/badge/doi-10.1038%2Fs41592--018--0235--4-blue.svg :target: https://doi.org/10.1038/s41592-018-0235-4 :alt: Published in Nature Methods diff --git a/docs/_static/alabaster.css b/docs/_static/alabaster.css index 0eddaeb07d..969ce31601 100644 --- a/docs/_static/alabaster.css +++ b/docs/_static/alabaster.css @@ -698,4 +698,4 @@ nav#breadcrumbs li+li:before { div.related { display: none; } -} \ No newline at end of file +} diff --git a/docs/_static/basic.css b/docs/_static/basic.css index 0807176ec0..b8dbd02b95 100644 --- a/docs/_static/basic.css +++ b/docs/_static/basic.css @@ -673,4 +673,4 @@ div.math:hover a.headerlink { #top-link { display: none; } -} \ No newline at end of file +} diff --git a/docs/_static/documentation_options.js b/docs/_static/documentation_options.js index d47d15f79f..3ed4a153d9 100644 --- a/docs/_static/documentation_options.js +++ b/docs/_static/documentation_options.js @@ -7,4 +7,4 @@ var DOCUMENTATION_OPTIONS = { HAS_SOURCE: true, SOURCELINK_SUFFIX: '.txt', NAVIGATION_WITH_KEYS: false, -}; \ No newline at end of file +}; diff --git a/docs/_static/language_data.js b/docs/_static/language_data.js index 5266fb19ec..f891ac1928 100644 --- a/docs/_static/language_data.js +++ b/docs/_static/language_data.js @@ -13,7 +13,7 @@ var stopwords = ["a","and","are","as","at","be","but","by","for","if","in","into","is","it","near","no","not","of","on","or","such","that","the","their","then","there","these","they","this","to","was","will","with"]; -/* Non-minified version JS is _stemmer.js if file is provided */ +/* Non-minified version JS is _stemmer.js if file is provided */ /** * Porter Stemmer */ @@ -293,5 +293,3 @@ function splitQuery(query) { } return result; } - - diff --git a/docs/_static/pygments.css b/docs/_static/pygments.css index dd6621d884..19bf02a9e8 100644 --- a/docs/_static/pygments.css +++ b/docs/_static/pygments.css @@ -74,4 +74,4 @@ .highlight .vg { color: #000000 } /* Name.Variable.Global */ .highlight .vi { color: #000000 } /* Name.Variable.Instance */ .highlight .vm { color: #000000 } /* Name.Variable.Magic */ -.highlight .il { color: #990000 } /* Literal.Number.Integer.Long */ \ No newline at end of file +.highlight .il { color: #990000 } /* Literal.Number.Integer.Long */ diff --git a/docs/_templates/apidoc/module.rst_t b/docs/_templates/apidoc/module.rst_t index 2490278551..49cca84bb8 100644 --- a/docs/_templates/apidoc/module.rst_t +++ b/docs/_templates/apidoc/module.rst_t @@ -6,4 +6,3 @@ {%- for option in automodule_options %} :{{ option }}: {%- endfor %} - diff --git a/docs/_templates/apidoc/toc.rst_t b/docs/_templates/apidoc/toc.rst_t index f0877eeb2f..878540cefe 100644 --- a/docs/_templates/apidoc/toc.rst_t +++ b/docs/_templates/apidoc/toc.rst_t @@ -5,4 +5,3 @@ {% for docname in docnames %} {{ docname }} {%- endfor %} - diff --git a/docs/api.rst b/docs/api.rst index b0c8205218..db8c6c94d5 100644 --- a/docs/api.rst +++ b/docs/api.rst @@ -8,4 +8,4 @@ Information on specific functions, classes, and methods. api/smriprep.cli api/smriprep.interfaces api/smriprep.utils - api/smriprep.workflows \ No newline at end of file + api/smriprep.workflows diff --git a/docs/conf.py b/docs/conf.py index 8405e5a2e3..0e0fefeea2 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -64,7 +64,6 @@ "seaborn", "skimage", "svgutils", - "templateflow", "transforms3d", ] diff --git a/docs/installation.rst b/docs/installation.rst index a4c6d204b1..13431762be 100644 --- a/docs/installation.rst +++ b/docs/installation.rst @@ -56,10 +56,10 @@ You can check your `Docker Engine`_ installation running their ``hello-world`` i $ docker run --rm hello-world If you have a functional installation, then you should obtain the following output. :: - + Hello from Docker! This message shows that your installation appears to be working correctly. - + To generate this message, Docker took the following steps: 1. The Docker client contacted the Docker daemon. 1. The Docker daemon pulled the "hello-world" image from the Docker Hub. @@ -68,20 +68,20 @@ If you have a functional installation, then you should obtain the following outp executable that produces the output you are currently reading. 1. The Docker daemon streamed that output to the Docker client, which sent it to your terminal. - + To try something more ambitious, you can run an Ubuntu container with: $ docker run -it ubuntu bash - + Share images, automate workflows, and more with a free Docker ID: https://hub.docker.com/ - + For more examples and ideas, visit: https://docs.docker.com/get-started/ After checking your Docker Engine is capable of running Docker images, then go ahead and `check out our documentation `__ to run the *sMRIPrep* image. -The list of Docker images ready to use is found at the `Docker Hub`_, +The list of Docker images ready to use is found at the `Docker Hub`_, under the ``nipreps/smriprep`` identifier. The ``smriprep-docker`` wrapper diff --git a/docs/usage.rst b/docs/usage.rst index 462502c542..f191307152 100644 --- a/docs/usage.rst +++ b/docs/usage.rst @@ -94,4 +94,3 @@ would be equivalent to the latest example: :: -v /path/to_output/dir:/out nipreps/smriprep:1.0.0 \ /data /out participant ... - diff --git a/pyproject.toml b/pyproject.toml index aa6f260430..f844c0cdd6 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -17,7 +17,7 @@ classifiers = [ "Programming Language :: Python :: 3.10", ] license = {file = "LICENSE"} -requires-python = ">=3.8" +requires-python = ">=3.10" dependencies = [ "importlib_resources >= 1.3; python_version < '3.9'", "indexed_gzip >= 0.8.8", diff --git a/smriprep/cli/run.py b/smriprep/cli/run.py index ffff602d21..10ce182814 100644 --- a/smriprep/cli/run.py +++ b/smriprep/cli/run.py @@ -72,7 +72,7 @@ def get_parser(): "output_dir", action="store", type=Path, - help="the output path for the outcomes of preprocessing and visual " "reports", + help="the output path for the outcomes of preprocessing and visual reports", ) parser.add_argument( "analysis_level", @@ -93,6 +93,15 @@ def get_parser(): help="a space delimited list of participant identifiers or a single " "identifier (the sub- prefix can be removed)", ) + g_bids.add_argument( + "-d", + "--derivatives", + action="store", + metavar="PATH", + type=Path, + nargs="*", + help="Search PATH(s) for pre-computed derivatives.", + ) g_bids.add_argument( "--bids-filter-file", action="store", @@ -132,7 +141,7 @@ def get_parser(): g_perfm.add_argument( "--low-mem", action="store_true", - help="attempt to reduce memory usage (will increase disk usage " "in working directory)", + help="attempt to reduce memory usage (will increase disk usage in working directory)", ) g_perfm.add_argument( "--use-plugin", @@ -205,6 +214,17 @@ def get_parser(): help="Path to existing FreeSurfer subjects directory to reuse. " "(default: OUTPUT_DIR/freesurfer)", ) + g_fs.add_argument( + "--cifti-output", + nargs="?", + const="91k", + default=False, + choices=("91k", "170k"), + type=str, + help="Output morphometry as CIFTI dense scalars. " + "Optionally, the number of grayordinate can be specified " + "(default is 91k, which equates to 2mm resolution)", + ) # Surface generation xor g_surfs = parser.add_argument_group("Surface preprocessing options") @@ -242,7 +262,8 @@ def get_parser(): "--fast-track", action="store_true", default=False, - help="fast-track the workflow by searching for existing derivatives.", + help="fast-track the workflow by searching for existing derivatives. " + "(DEPRECATED for --derivatives).", ) g_other.add_argument( "--resource-monitor", @@ -274,7 +295,7 @@ def get_parser(): "--stop-on-first-crash", action="store_true", default=False, - help="Force stopping on first crash, even if a work directory" " was specified.", + help="Force stopping on first crash, even if a work directory was specified.", ) g_other.add_argument( "--notrack", @@ -421,6 +442,7 @@ def build_workflow(opts, retval): from shutil import copyfile from os import cpu_count import uuid + import warnings from time import strftime from subprocess import check_call, CalledProcessError, TimeoutExpired from pkg_resources import resource_filename as pkgrf @@ -562,11 +584,20 @@ def build_workflow(opts, retval): ), ) + derivatives = opts.derivatives or [] + if opts.fast_track: + # XXX Makes strong assumption of legacy layout + smriprep_dir = str(output_dir / "smriprep") + warnings.warn( + f"Received DEPRECATED --fast-track flag. Adding {smriprep_dir} to --derivatives list." + ) + derivatives.append(smriprep_dir) + # Build main workflow retval["workflow"] = init_smriprep_wf( sloppy=opts.sloppy, debug=False, - fast_track=opts.fast_track, + derivatives=derivatives, freesurfer=opts.run_reconall, fs_subjects_dir=opts.fs_subjects_dir, hires=opts.hires, @@ -584,6 +615,7 @@ def build_workflow(opts, retval): subject_list=subject_list, work_dir=str(work_dir), bids_filters=bids_filters, + cifti_output=opts.cifti_output, ) retval["return_code"] = 0 diff --git a/smriprep/data/io_spec.json b/smriprep/data/io_spec.json index 1185a54db7..6b8fcd50fc 100644 --- a/smriprep/data/io_spec.json +++ b/smriprep/data/io_spec.json @@ -3,94 +3,158 @@ "baseline": { "preproc": { "datatype": "anat", + "space": null, "desc": "preproc", - "suffix": "T1w" + "suffix": "T1w", + "extension": [ + ".nii.gz", + ".nii" + ] }, "mask": { "datatype": "anat", + "space": null, "desc": "brain", - "suffix": "mask" + "suffix": "mask", + "extension": [ + ".nii.gz", + ".nii" + ] }, "dseg": { "datatype": "anat", - "suffix": "dseg" + "space": null, + "desc": null, + "suffix": "dseg", + "extension": [ + ".nii.gz", + ".nii" + ] }, "tpms": { "datatype": "anat", - "label": ["GM", "WM", "CSF"], - "suffix": "probseg" + "space": null, + "label": [ + "GM", + "WM", + "CSF" + ], + "suffix": "probseg", + "extension": [ + ".nii.gz", + ".nii" + ] } }, - "std_xfms": { - "anat2std_xfm": { + "transforms": { + "forward": { "datatype": "anat", - "extension": "h5", + "extension": [ + "h5", + "txt" + ], "from": "T1w", - "to": [], + "to": null, "suffix": "xfm", "mode": "image" }, - "std2anat_xfm": { + "reverse": { "datatype": "anat", - "extension": "h5", - "from": [], + "extension": [ + "h5", + "txt" + ], + "from": null, "to": "T1w", "suffix": "xfm", "mode": "image" } }, "surfaces": { - "t1w_aseg": { + "white": { "datatype": "anat", - "desc": "aseg", - "suffix": "dseg" + "hemi": ["L", "R"], + "space": null, + "suffix": "white", + "extension": ".surf.gii" }, - "t1w_aparc": { + "pial": { "datatype": "anat", - "desc": "aparcaseg", - "suffix": "dseg" + "hemi": ["L", "R"], + "space": null, + "suffix": "pial", + "extension": ".surf.gii" }, - "t1w2fsnative_xfm": { + "midthickness": { "datatype": "anat", - "from": "T1w", - "to": "fsnative", - "extension": "txt", - "suffix": "xfm", - "mode": "image" + "hemi": ["L", "R"], + "space": null, + "suffix": "midthickness", + "extension": ".surf.gii" }, - "fsnative2t1w_xfm": { + "sphere": { "datatype": "anat", - "from": "fsnative", - "to": "T1w", - "extension": "txt", - "suffix": "xfm", - "mode": "image" + "hemi": ["L", "R"], + "space": null, + "desc": null, + "suffix": "sphere", + "extension": ".surf.gii" }, - "surfaces": { + "thickness": { "datatype": "anat", "hemi": ["L", "R"], - "extension": "surf.gii", - "suffix": [ "inflated", "midthickness", "pial", "white"] + "space": null, + "suffix": "thickness", + "extension": ".shape.gii" }, - "morphometrics": { + "sulc": { "datatype": "anat", "hemi": ["L", "R"], - "extension": "shape.gii", - "suffix": ["thickness", "sulc", "curv"] + "space": null, + "suffix": "sulc", + "extension": ".shape.gii" }, - "anat_ribbon": { + "curv": { "datatype": "anat", - "desc": "ribbon", - "suffix": "mask", - "extension": "nii.gz" + "hemi": ["L", "R"], + "space": null, + "suffix": "curv", + "extension": ".shape.gii" + }, + "sphere_reg": { + "datatype": "anat", + "hemi": ["L", "R"], + "space": null, + "desc": "reg", + "suffix": "sphere", + "extension": ".surf.gii" }, "sphere_reg_fsLR": { - "datatype" : "anat", + "datatype": "anat", "hemi": ["L", "R"], "space": "fsLR", "desc": "reg", "suffix": "sphere", - "extension": "surf.gii" + "extension": ".surf.gii" + }, + "sphere_reg_msm": { + "datatype": "anat", + "hemi": ["L", "R"], + "space": "fsLR", + "desc": "msmsulc", + "suffix": "sphere", + "extension": ".surf.gii" + } + }, + "masks": { + "anat_ribbon": { + "datatype": "anat", + "desc": "ribbon", + "suffix": "mask", + "extension": [ + ".nii.gz", + ".nii" + ] } } }, diff --git a/smriprep/interfaces/freesurfer.py b/smriprep/interfaces/freesurfer.py index b363c92008..9d00063350 100644 --- a/smriprep/interfaces/freesurfer.py +++ b/smriprep/interfaces/freesurfer.py @@ -293,13 +293,21 @@ def _gen_filename(self, name): return super()._gen_filename(name) -class MakeMidthickness(nwfs.MakeMidthickness): +class MakeMidthicknessInputSpec(nwfs._MakeMidthicknessInputSpec, fs.base.FSTraitedSpecOpenMP): + pass + + +class MakeMidthickness(nwfs.MakeMidthickness, fs.base.FSCommandOpenMP): """Patched MakeMidthickness interface Ensures output filenames are specified with hemisphere labels, when appropriate. This may not cover all use-cases in MRIsExpand, but we're just making midthickness files. + This interface also sets the OMP_NUM_THREADS environment variable to floor(1.5x) the + number of threads requested by the user, as tests indicate that cores are underutilized + by a factor of 2/3. + >>> from smriprep.interfaces.freesurfer import MakeMidthickness >>> mris_expand = MakeMidthickness(thickness=True, distance=0.5) >>> mris_expand.inputs.in_file = 'lh.white' @@ -316,9 +324,17 @@ class MakeMidthickness(nwfs.MakeMidthickness): 'mris_expand -thickness lh.white 0.5 rh.graymid' """ + input_spec = MakeMidthicknessInputSpec + def _format_arg(self, name, trait_spec, value): # FreeSurfer at some point changed whether it would add the hemi label onto the # surface. Therefore we'll do it ourselves. if name == "out_name": value = self._associated_file(self.inputs.in_file, value) return super()._format_arg(name, trait_spec, value) + + def _num_threads_update(self): + if self.inputs.num_threads: + self.inputs.environ.update( + {"OMP_NUM_THREADS": str(self.inputs.num_threads * 3 // 2)} + ) diff --git a/smriprep/interfaces/gifti.py b/smriprep/interfaces/gifti.py index 22f9925201..3d9d33c008 100644 --- a/smriprep/interfaces/gifti.py +++ b/smriprep/interfaces/gifti.py @@ -5,7 +5,7 @@ from nipype.interfaces.base import File, SimpleInterface, TraitedSpec, isdefined, traits -class InvertShapeInputSpec(TraitedSpec): +class MetricMathInputSpec(TraitedSpec): subject_id = traits.Str(desc='subject ID') hemisphere = traits.Enum( "L", @@ -13,50 +13,68 @@ class InvertShapeInputSpec(TraitedSpec): mandatory=True, desc='hemisphere', ) - shape = traits.Str(desc='name of shape to invert') - shape_file = File(exists=True, mandatory=True, desc='input GIFTI file') + metric = traits.Str(desc='name of metric to invert') + metric_file = File(exists=True, mandatory=True, desc='input GIFTI file') + operation = traits.Enum( + "invert", + "abs", + "bin", + mandatory=True, + desc='operation to perform', + ) -class InvertShapeOutputSpec(TraitedSpec): - shape_file = File(desc='output GIFTI file') +class MetricMathOutputSpec(TraitedSpec): + metric_file = File(desc='output GIFTI file') -class InvertShape(SimpleInterface): - """Prepare GIFTI shape file for use in MSMSulc +class MetricMath(SimpleInterface): + """Prepare GIFTI metric file for use in MSMSulc This interface mirrors the action of the following portion of FreeSurfer2CaretConvertAndRegisterNonlinear.sh:: - wb_command -set-structure ${shape_file} CORTEX_[LEFT|RIGHT] - wb_command -metric-math "var * -1" ${shape_file} -var var ${shape_file} - wb_command -set-map-names ${shape_file} -map 1 ${subject}_[L|R]_${shape} + wb_command -set-structure ${metric_file} CORTEX_[LEFT|RIGHT] + wb_command -metric-math "var * -1" ${metric_file} -var var ${metric_file} + wb_command -set-map-names ${metric_file} -map 1 ${subject}_[L|R]_${metric} + # If abs: + wb_command -metric-math "abs(var)" ${metric_file} -var var ${metric_file} We do not add palette information to the output file. """ - input_spec = InvertShapeInputSpec - output_spec = InvertShapeOutputSpec + input_spec = MetricMathInputSpec + output_spec = MetricMathOutputSpec def _run_interface(self, runtime): - subject, hemi, shape = self.inputs.subject_id, self.inputs.hemisphere, self.inputs.shape + subject, hemi, metric = self.inputs.subject_id, self.inputs.hemisphere, self.inputs.metric if not isdefined(subject): subject = 'sub-XYZ' - img = nb.GiftiImage.from_filename(self.inputs.shape_file) + img = nb.GiftiImage.from_filename(self.inputs.metric_file) # wb_command -set-structure img.meta["AnatomicalStructurePrimary"] = {'L': 'CortexLeft', 'R': 'CortexRight'}[hemi] darray = img.darrays[0] # wb_command -set-map-names meta = darray.meta - meta['Name'] = f"{subject}_{hemi}_{shape}" - - # wb_command -metric-math "var * -1" - inv = -darray.data + meta['Name'] = f"{subject}_{hemi}_{metric}" + + datatype = darray.datatype + if self.inputs.operation == "abs": + # wb_command -metric-math "abs(var)" + data = abs(darray.data) + elif self.inputs.operation == "invert": + # wb_command -metric-math "var * -1" + data = -darray.data + elif self.inputs.operation == "bin": + # wb_command -metric-math "var > 0" + data = darray.data > 0 + datatype = 'uint8' darray = nb.gifti.GiftiDataArray( - inv, + data, intent=darray.intent, - datatype=darray.datatype, + datatype=datatype, encoding=darray.encoding, endian=darray.endian, coordsys=darray.coordsys, @@ -65,7 +83,7 @@ def _run_interface(self, runtime): ) img.darrays[0] = darray - out_filename = os.path.join(runtime.cwd, f"{subject}.{hemi}.{shape}.native.shape.gii") + out_filename = os.path.join(runtime.cwd, f"{subject}.{hemi}.{metric}.native.shape.gii") img.to_filename(out_filename) - self._results["shape_file"] = out_filename + self._results["metric_file"] = out_filename return runtime diff --git a/smriprep/interfaces/surf.py b/smriprep/interfaces/surf.py index bd91db0966..369dd638c0 100644 --- a/smriprep/interfaces/surf.py +++ b/smriprep/interfaces/surf.py @@ -22,10 +22,12 @@ # """Handling surfaces.""" import os +from pathlib import Path from typing import Optional import numpy as np import nibabel as nb +import nitransforms as nt from nipype.interfaces.base import ( BaseInterfaceInputSpec, @@ -40,7 +42,7 @@ class _NormalizeSurfInputSpec(BaseInterfaceInputSpec): in_file = File(mandatory=True, exists=True, desc="Freesurfer-generated GIFTI file") - transform_file = File(exists=True, desc="FSL or LTA affine transform file") + transform_file = File(exists=True, desc="FSL, LTA or ITK affine transform file") class _NormalizeSurfOutputSpec(TraitedSpec): @@ -155,7 +157,35 @@ def _run_interface(self, runtime): return runtime -def normalize_surfs(in_file: str, transform_file: str, newpath: Optional[str] = None) -> str: +class MakeRibbonInputSpec(TraitedSpec): + white_distvols = traits.List( + File(exists=True), minlen=2, maxlen=2, desc="White matter distance volumes" + ) + pial_distvols = traits.List( + File(exists=True), minlen=2, maxlen=2, desc="Pial matter distance volumes" + ) + + +class MakeRibbonOutputSpec(TraitedSpec): + ribbon = File(desc="Binary ribbon mask") + + +class MakeRibbon(SimpleInterface): + """Create a binary ribbon mask from white and pial distance volumes.""" + + input_spec = MakeRibbonInputSpec + output_spec = MakeRibbonOutputSpec + + def _run_interface(self, runtime): + self._results["ribbon"] = make_ribbon( + self.inputs.white_distvols, self.inputs.pial_distvols, newpath=runtime.cwd + ) + return runtime + + +def normalize_surfs( + in_file: str, transform_file: str | None, newpath: Optional[str] = None +) -> str: """ Update GIFTI metadata and apply rigid coordinate correction. @@ -168,29 +198,44 @@ def normalize_surfs(in_file: str, transform_file: str, newpath: Optional[str] = """ img = nb.load(in_file) - transform = load_transform(transform_file) + if transform_file is None: + transform = nt.linear.Affine() + else: + xfm_fmt = { + ".txt": "itk", + ".mat": "fsl", + ".lta": "fs", + }[Path(transform_file).suffix] + transform = nt.linear.load(transform_file, fmt=xfm_fmt) pointset = img.get_arrays_from_intent("NIFTI_INTENT_POINTSET")[0] - if not np.allclose(transform, np.eye(4)): - pointset.data = nb.affines.apply_affine(transform, pointset.data) - - # mris_convert --to-scanner removes VolGeom transform from coordinates, - # but not metadata. - # We could set to default LIA affine, but there seems little advantage - for XYZC in "XYZC": - for RAS in "RAS": - pointset.meta.pop(f"VolGeom{XYZC}_{RAS}", None) + if not np.allclose(transform.matrix, np.eye(4)): + pointset.data = transform.map(pointset.data, inverse=True) fname = os.path.basename(in_file) - if "midthickness" in fname.lower() or "graymid" in fname.lower(): + if "graymid" in fname.lower(): + # Rename graymid to midthickness + fname = fname.replace("graymid", "midthickness") + if "midthickness" in fname.lower(): pointset.meta.setdefault("AnatomicalStructureSecondary", "MidThickness") pointset.meta.setdefault("GeometricType", "Anatomical") # FreeSurfer incorrectly uses "Sphere" for spherical surfaces if pointset.meta.get("GeometricType") == "Sphere": pointset.meta["GeometricType"] = "Spherical" - - if newpath is not None: + else: + # mris_convert --to-scanner removes VolGeom transform from coordinates, + # but not metadata. + # We could set to default LIA affine, but there seems little advantage. + # + # Following the lead of HCP pipelines, we only adjust the coordinates + # for anatomical surfaces. To ensure consistent treatment by FreeSurfer, + # we leave the metadata for spherical surfaces intact. + for XYZC in "XYZC": + for RAS in "RAS": + pointset.meta.pop(f"VolGeom{XYZC}_{RAS}", None) + + if newpath is None: newpath = os.getcwd() out_file = os.path.join(newpath, fname) img.to_filename(out_file) @@ -214,37 +259,32 @@ def fix_gifti_metadata(in_file: str, newpath: Optional[str] = None) -> str: if pointset.meta.get("GeometricType") == "Sphere": pointset.meta["GeometricType"] = "Spherical" - if newpath is not None: + if newpath is None: newpath = os.getcwd() out_file = os.path.join(newpath, os.path.basename(in_file)) img.to_filename(out_file) return out_file -def load_transform(fname: str) -> np.ndarray: - """Load affine transform from file +def make_ribbon( + white_distvols: list[str], + pial_distvols: list[str], + newpath: Optional[str] = None, +) -> str: + base_img = nb.load(white_distvols[0]) + header = base_img.header + header.set_data_dtype("uint8") - Parameters - ---------- - fname : str or None - Filename of an LTA or FSL-style MAT transform file. - If ``None``, return an identity transform + ribbons = [ + (np.array(nb.load(white).dataobj) > 0) & (np.array(nb.load(pial).dataobj) < 0) + for white, pial in zip(white_distvols, pial_distvols) + ] - Returns - ------- - affine : (4, 4) numpy.ndarray - """ - if fname is None: - return np.eye(4) - - if fname.endswith(".mat"): - return np.loadtxt(fname) - elif fname.endswith(".lta"): - with open(fname, "rb") as fobj: - for line in fobj: - if line.startswith(b"1 4 4"): - break - lines = fobj.readlines()[:4] - return np.genfromtxt(lines) - - raise ValueError("Unknown transform type; pass FSL (.mat) or LTA (.lta)") + if newpath is None: + newpath = os.getcwd() + out_file = os.path.join(newpath, "ribbon.nii.gz") + + ribbon_data = ribbons[0] | ribbons[1] + ribbon = base_img.__class__(ribbon_data, base_img.affine, header) + ribbon.to_filename(out_file) + return out_file diff --git a/smriprep/interfaces/tests/__init__.py b/smriprep/interfaces/tests/__init__.py new file mode 100644 index 0000000000..e69de29bb2 diff --git a/smriprep/interfaces/tests/data/sub-fsaverage_res-4_desc-cropped_ribbon.nii.gz b/smriprep/interfaces/tests/data/sub-fsaverage_res-4_desc-cropped_ribbon.nii.gz new file mode 100644 index 0000000000..7f8f3fefda Binary files /dev/null and b/smriprep/interfaces/tests/data/sub-fsaverage_res-4_desc-cropped_ribbon.nii.gz differ diff --git a/smriprep/interfaces/tests/data/sub-fsaverage_res-4_hemi-L_desc-cropped_pialdist.nii.gz b/smriprep/interfaces/tests/data/sub-fsaverage_res-4_hemi-L_desc-cropped_pialdist.nii.gz new file mode 100644 index 0000000000..b5904717c7 Binary files /dev/null and b/smriprep/interfaces/tests/data/sub-fsaverage_res-4_hemi-L_desc-cropped_pialdist.nii.gz differ diff --git a/smriprep/interfaces/tests/data/sub-fsaverage_res-4_hemi-L_desc-cropped_wmdist.nii.gz b/smriprep/interfaces/tests/data/sub-fsaverage_res-4_hemi-L_desc-cropped_wmdist.nii.gz new file mode 100644 index 0000000000..0d77f1a0e1 Binary files /dev/null and b/smriprep/interfaces/tests/data/sub-fsaverage_res-4_hemi-L_desc-cropped_wmdist.nii.gz differ diff --git a/smriprep/interfaces/tests/data/sub-fsaverage_res-4_hemi-R_desc-cropped_pialdist.nii.gz b/smriprep/interfaces/tests/data/sub-fsaverage_res-4_hemi-R_desc-cropped_pialdist.nii.gz new file mode 100644 index 0000000000..217c530c18 Binary files /dev/null and b/smriprep/interfaces/tests/data/sub-fsaverage_res-4_hemi-R_desc-cropped_pialdist.nii.gz differ diff --git a/smriprep/interfaces/tests/data/sub-fsaverage_res-4_hemi-R_desc-cropped_wmdist.nii.gz b/smriprep/interfaces/tests/data/sub-fsaverage_res-4_hemi-R_desc-cropped_wmdist.nii.gz new file mode 100644 index 0000000000..a645d827dd Binary files /dev/null and b/smriprep/interfaces/tests/data/sub-fsaverage_res-4_hemi-R_desc-cropped_wmdist.nii.gz differ diff --git a/smriprep/interfaces/tests/test_surf.py b/smriprep/interfaces/tests/test_surf.py new file mode 100644 index 0000000000..d8fedca5ce --- /dev/null +++ b/smriprep/interfaces/tests/test_surf.py @@ -0,0 +1,36 @@ +import nibabel as nb +from nipype.pipeline import engine as pe +import numpy as np + +from ...data import load_resource +from ..surf import MakeRibbon + + +def test_MakeRibbon(tmp_path): + res_template = "{path}/sub-fsaverage_res-4_hemi-{hemi}_desc-cropped_{surf}dist.nii.gz" + white, pial = [ + [ + load_resource( + res_template.format(path="../interfaces/tests/data", hemi=hemi, surf=surf) + ) + for hemi in "LR" + ] + for surf in ("wm", "pial") + ] + + make_ribbon = pe.Node( + MakeRibbon(white_distvols=white, pial_distvols=pial), + name="make_ribbon", + base_dir=tmp_path, + ) + + result = make_ribbon.run() + + ribbon = nb.load(result.outputs.ribbon) + expected = nb.load( + load_resource("../interfaces/tests/data/sub-fsaverage_res-4_desc-cropped_ribbon.nii.gz") + ) + + assert ribbon.shape == expected.shape + assert np.allclose(ribbon.affine, expected.affine) + assert np.array_equal(ribbon.dataobj, expected.dataobj) diff --git a/smriprep/utils/bids.py b/smriprep/utils/bids.py index 99cfd8ef53..5b9174d6dd 100644 --- a/smriprep/utils/bids.py +++ b/smriprep/utils/bids.py @@ -21,111 +21,13 @@ # https://www.nipreps.org/community/licensing/ # """Utilities to handle BIDS inputs.""" -from collections import defaultdict from pathlib import Path from json import loads from pkg_resources import resource_filename as pkgrf -from bids.layout.writing import build_path +from bids.layout import BIDSLayout -def get_outputnode_spec(): - """ - Generate outputnode's fields from I/O spec file. - - Examples - -------- - >>> get_outputnode_spec() # doctest: +NORMALIZE_WHITESPACE - ['t1w_preproc', 't1w_mask', 't1w_dseg', 't1w_tpms', - 'std_preproc', 'std_mask', 'std_dseg', 'std_tpms', - 'anat2std_xfm', 'std2anat_xfm', - 't1w_aseg', 't1w_aparc', - 't1w2fsnative_xfm', 'fsnative2t1w_xfm', - 'surfaces', 'morphometrics', 'anat_ribbon', 'sphere_reg_fsLR'] - - """ - spec = loads(Path(pkgrf("smriprep", "data/io_spec.json")).read_text())["queries"] - fields = ["_".join((m, s)) for m in ("t1w", "std") for s in spec["baseline"].keys()] - fields += [s for s in spec["std_xfms"].keys()] - fields += [s for s in spec["surfaces"].keys()] - return fields - - -def predict_derivatives(subject_id, output_spaces, freesurfer): - """ - Generate a list of the files that should be found in the output folder. - - The prediction of outputs serves two purposes: - - * Anticipates the set of outputs that will be generated. - * Informs the decision of whether the workflow should be staged or some - sections can be skip. - - Parameters - ---------- - subject_id : :obj:`str` - A subject id - output_spaces : :obj:`list` - TemplateFlow identifiers of the requested output spaces - freesurfer : :obj:`bool` - Whether the ``--fs-no-reconall`` flag was used. - - Examples - -------- - >>> predict_derivatives('01', ['MNI152NLin2009cAsym'], False) # doctest: +NORMALIZE_WHITESPACE - ['sub-01/anat/sub-01_desc-brain_mask.nii.gz', - 'sub-01/anat/sub-01_desc-preproc_T1w.nii.gz', - 'sub-01/anat/sub-01_dseg.nii.gz', - 'sub-01/anat/sub-01_from-MNI152NLin2009cAsym_to-T1w_mode-image_xfm.h5', - 'sub-01/anat/sub-01_from-T1w_to-MNI152NLin2009cAsym_mode-image_xfm.h5', - 'sub-01/anat/sub-01_label-CSF_probseg.nii.gz', - 'sub-01/anat/sub-01_label-GM_probseg.nii.gz', - 'sub-01/anat/sub-01_label-WM_probseg.nii.gz', - 'sub-01/anat/sub-01_space-MNI152NLin2009cAsym_desc-brain_mask.nii.gz', - 'sub-01/anat/sub-01_space-MNI152NLin2009cAsym_desc-preproc_T1w.nii.gz', - 'sub-01/anat/sub-01_space-MNI152NLin2009cAsym_dseg.nii.gz', - 'sub-01/anat/sub-01_space-MNI152NLin2009cAsym_label-CSF_probseg.nii.gz', - 'sub-01/anat/sub-01_space-MNI152NLin2009cAsym_label-GM_probseg.nii.gz', - 'sub-01/anat/sub-01_space-MNI152NLin2009cAsym_label-WM_probseg.nii.gz'] - - """ - spec = loads(Path(pkgrf("smriprep", "data/io_spec.json")).read_text()) - - def _normalize_q(query, space=None): - query = query.copy() - query["subject"] = subject_id - if space is not None: - query["space"] = space - if "from" in query and not query["from"]: - query["from"] = output_spaces - if "to" in query and not query["to"]: - query["to"] = output_spaces - return query - - queries = [_normalize_q(q, space=None) for q in spec["queries"]["baseline"].values()] - - queries += [ - _normalize_q(q, space=s) - for s in output_spaces - for q in spec["queries"]["baseline"].values() - ] - queries += [_normalize_q(q) for q in spec["queries"]["std_xfms"].values()] - if freesurfer: - queries += [_normalize_q(q) for q in spec["queries"]["surfaces"].values()] - - output = [] - for q in queries: - paths = build_path(q, spec["patterns"]) - if isinstance(paths, str): - output.append(paths) - elif paths: - output += paths - - return sorted(output) - - -def collect_derivatives( - derivatives_dir, subject_id, std_spaces, freesurfer, spec=None, patterns=None -): +def collect_derivatives(derivatives_dir, subject_id, std_spaces, spec=None, patterns=None): """Gather existing derivatives and compose a cache.""" if spec is None or patterns is None: _spec, _patterns = tuple( @@ -137,64 +39,37 @@ def collect_derivatives( if patterns is None: patterns = _patterns - derivs_cache = defaultdict(list, {}) - derivatives_dir = Path(derivatives_dir) + layout = BIDSLayout(derivatives_dir, config=["bids", "derivatives"], validate=False) - def _check_item(item): + derivs_cache = {} + for k, q in spec["baseline"].items(): + q["subject"] = subject_id + item = layout.get(return_type='filename', **q) if not item: - return None - - if isinstance(item, str): - item = [item] - - result = [] - for i in item: - if not (derivatives_dir / i).exists(): - i = i.rstrip(".gz") - if not (derivatives_dir / i).exists(): - return None - result.append(str(derivatives_dir / i)) + continue - return result - - for space in [None] + std_spaces: - for k, q in spec["baseline"].items(): - q["subject"] = subject_id - if space is not None: - q["space"] = space - item = _check_item(build_path(q, patterns, strict=True)) - if not item: - return None - - if space: - derivs_cache["std_%s" % k] += item if len(item) == 1 else [item] - else: - derivs_cache["t1w_%s" % k] = item[0] if len(item) == 1 else item + derivs_cache["t1w_%s" % k] = item[0] if len(item) == 1 else item + transforms = derivs_cache.setdefault('transforms', {}) for space in std_spaces: - for k, q in spec["std_xfms"].items(): + for k, q in spec["transforms"].items(): + q = q.copy() q["subject"] = subject_id q["from"] = q["from"] or space q["to"] = q["to"] or space - item = _check_item(build_path(q, patterns)) + item = layout.get(return_type='filename', **q) if not item: - return None - derivs_cache[k] += item + continue + transforms.setdefault(space, {})[k] = item[0] if len(item) == 1 else item - derivs_cache = dict(derivs_cache) # Back to a standard dictionary - - if freesurfer: - for k, q in spec["surfaces"].items(): - q["subject"] = subject_id - item = _check_item(build_path(q, patterns)) - if not item: - return None + for k, q in spec["surfaces"].items(): + q["subject"] = subject_id + item = layout.get(return_type='filename', **q) + if not item or len(item) != 2: + continue - if len(item) == 1: - item = item[0] - derivs_cache[k] = item + derivs_cache[k] = sorted(item) - derivs_cache["template"] = std_spaces return derivs_cache diff --git a/smriprep/workflows/anatomical.py b/smriprep/workflows/anatomical.py index 738ec03941..a6b90f6d68 100644 --- a/smriprep/workflows/anatomical.py +++ b/smriprep/workflows/anatomical.py @@ -21,6 +21,8 @@ # https://www.nipreps.org/community/licensing/ # """Anatomical reference preprocessing workflows.""" +import typing as ty + from pkg_resources import resource_filename as pkgr from nipype import logging @@ -43,19 +45,41 @@ ) from niworkflows.interfaces.header import ValidateImage from niworkflows.interfaces.images import TemplateDimensions, Conform +from niworkflows.interfaces.nibabel import ApplyMask, Binarize from niworkflows.interfaces.nitransforms import ConcatenateXFMs -from niworkflows.interfaces.utility import KeySelect -from niworkflows.utils.misc import fix_multi_T1w_source_name, add_suffix +from niworkflows.utils.spaces import SpatialReferences, Reference +from niworkflows.utils.misc import add_suffix from niworkflows.anat.ants import init_brain_extraction_wf, init_n4_only_wf -from ..utils.bids import get_outputnode_spec +from ..interfaces import DerivativesDataSink from ..utils.misc import apply_lut as _apply_bids_lut, fs_isRunning as _fs_isRunning -from .norm import init_anat_norm_wf -from .outputs import init_anat_reports_wf, init_anat_derivatives_wf +from .fit.registration import init_register_template_wf +from .outputs import ( + init_anat_reports_wf, + init_ds_anat_volumes_wf, + init_ds_grayord_metrics_wf, + init_ds_surface_metrics_wf, + init_ds_surfaces_wf, + init_ds_template_wf, + init_ds_mask_wf, + init_ds_dseg_wf, + init_ds_tpms_wf, + init_ds_template_registration_wf, + init_ds_fs_registration_wf, + init_anat_second_derivatives_wf, + init_template_iterator_wf, +) from .surfaces import ( init_anat_ribbon_wf, - init_sphere_reg_wf, - init_surface_recon_wf, + init_fsLR_reg_wf, + init_gifti_morphometrics_wf, + init_hcp_morphometrics_wf, + init_gifti_surfaces_wf, init_morph_grayords_wf, + init_msm_sulc_wf, + init_surface_derivatives_wf, + init_surface_recon_wf, + init_refinement_wf, + init_resample_midthickness_wf, ) LOGGER = logging.getLogger("nipype.workflow") @@ -63,24 +87,369 @@ def init_anat_preproc_wf( *, - bids_root, - freesurfer, - hires, - longitudinal, - t1w, - t2w, - omp_nthreads, - output_dir, - skull_strip_mode, - skull_strip_template, - spaces, - cifti_output=False, - debug=False, - sloppy=False, - msm_sulc=False, - existing_derivatives=None, - name="anat_preproc_wf", - skull_strip_fixed_seed=False, + bids_root: str, + output_dir: str, + freesurfer: bool, + hires: bool, + longitudinal: bool, + msm_sulc: bool, + t1w: list, + t2w: list, + skull_strip_mode: str, + skull_strip_template: Reference, + spaces: SpatialReferences, + precomputed: dict, + omp_nthreads: int, + debug: bool = False, + sloppy: bool = False, + cifti_output: ty.Literal["91k", "170k", False] = False, + name: str = "anat_preproc_wf", + skull_strip_fixed_seed: bool = False, +): + """ + Stage the anatomical preprocessing steps of *sMRIPrep*. + + This workflow is a compatibility wrapper around :py:func:`init_anat_fit_wf` + that emits all derivatives that were present in sMRIPrep 0.9.x and before. + + Workflow Graph + .. workflow:: + :graph2use: orig + :simple_form: yes + + from niworkflows.utils.spaces import SpatialReferences, Reference + from smriprep.workflows.anatomical import init_anat_preproc_wf + spaces = SpatialReferences(spaces=['MNI152NLin2009cAsym', 'fsaverage5']) + spaces.checkpoint() + wf = init_anat_preproc_wf( + bids_root='.', + output_dir='.', + freesurfer=True, + hires=True, + longitudinal=False, + msm_sulc=False, + t1w=['t1w.nii.gz'], + t2w=[], + skull_strip_mode='force', + skull_strip_template=Reference('OASIS30ANTs'), + spaces=spaces, + precomputed={}, + omp_nthreads=1, + ) + + + Parameters + ---------- + bids_root : :obj:`str` + Path of the input BIDS dataset root + output_dir : :obj:`str` + Directory in which to save derivatives + freesurfer : :obj:`bool` + Enable FreeSurfer surface reconstruction (increases runtime by 6h, + at the very least) + hires : :obj:`bool` + Enable sub-millimeter preprocessing in FreeSurfer + longitudinal : :obj:`bool` + Create unbiased structural template, regardless of number of inputs + (may increase runtime) + t1w : :obj:`list` + List of T1-weighted structural images. + skull_strip_mode : :obj:`str` + Determiner for T1-weighted skull stripping (`force` ensures skull stripping, + `skip` ignores skull stripping, and `auto` automatically ignores skull stripping + if pre-stripped brains are detected). + skull_strip_template : :py:class:`~niworkflows.utils.spaces.Reference` + Spatial reference to use in atlas-based brain extraction. + spaces : :py:class:`~niworkflows.utils.spaces.SpatialReferences` + Object containing standard and nonstandard space specifications. + precomputed : :obj:`dict` + Dictionary mapping output specification attribute names and + paths to precomputed derivatives. + omp_nthreads : :obj:`int` + Maximum number of threads an individual process may use + debug : :obj:`bool` + Enable debugging outputs + sloppy: :obj:`bool` + Quick, impercise operations. Used to decrease workflow duration. + name : :obj:`str`, optional + Workflow name (default: anat_fit_wf) + skull_strip_fixed_seed : :obj:`bool` + Do not use a random seed for skull-stripping - will ensure + run-to-run replicability when used with --omp-nthreads 1 + (default: ``False``). + + Inputs + ------ + t1w + List of T1-weighted structural images + t2w + List of T2-weighted structural images + roi + A mask to exclude regions during standardization + flair + List of FLAIR images + subjects_dir + FreeSurfer SUBJECTS_DIR + subject_id + FreeSurfer subject ID + + Outputs + ------- + t1w_preproc + The T1w reference map, which is calculated as the average of bias-corrected + and preprocessed T1w images, defining the anatomical space. + t1w_mask + Brain (binary) mask estimated by brain extraction. + t1w_dseg + Brain tissue segmentation of the preprocessed structural image, including + gray-matter (GM), white-matter (WM) and cerebrospinal fluid (CSF). + t1w_tpms + List of tissue probability maps corresponding to ``t1w_dseg``. + template + List of template names to which the structural image has been registered + anat2std_xfm + List of nonlinear spatial transforms to resample data from subject + anatomical space into standard template spaces. Collated with template. + std2anat_xfm + List of nonlinear spatial transforms to resample data from standard + template spaces into subject anatomical space. Collated with template. + subjects_dir + FreeSurfer SUBJECTS_DIR; use as input to a node to ensure that it is run after + FreeSurfer reconstruction is completed. + subject_id + FreeSurfer subject ID; use as input to a node to ensure that it is run after + FreeSurfer reconstruction is completed. + fsnative2t1w_xfm + ITK-style affine matrix translating from FreeSurfer-conformed subject space to T1w + + """ + workflow = Workflow(name=name) + + inputnode = pe.Node( + niu.IdentityInterface(fields=["t1w", "t2w", "roi", "flair", "subjects_dir", "subject_id"]), + name="inputnode", + ) + outputnode = pe.Node( + niu.IdentityInterface( + fields=[ + "template", + "subjects_dir", + "subject_id", + "t1w_preproc", + "t1w_mask", + "t1w_dseg", + "t1w_tpms", + "anat2std_xfm", + "std2anat_xfm", + "fsnative2t1w_xfm", + "t1w_aparc", + "t1w_aseg", + "sphere_reg", + "sphere_reg_fsLR", + ] + ), + name="outputnode", + ) + + anat_fit_wf = init_anat_fit_wf( + bids_root=bids_root, + output_dir=output_dir, + freesurfer=freesurfer, + hires=hires, + longitudinal=longitudinal, + msm_sulc=msm_sulc, + skull_strip_mode=skull_strip_mode, + skull_strip_template=skull_strip_template, + spaces=spaces, + t1w=t1w, + t2w=t2w, + precomputed=precomputed, + debug=debug, + sloppy=sloppy, + omp_nthreads=omp_nthreads, + skull_strip_fixed_seed=skull_strip_fixed_seed, + ) + template_iterator_wf = init_template_iterator_wf(spaces=spaces) + ds_std_volumes_wf = init_ds_anat_volumes_wf( + bids_root=bids_root, + output_dir=output_dir, + name="ds_std_volumes_wf", + ) + + workflow.connect([ + (inputnode, anat_fit_wf, [ + ("t1w", "inputnode.t1w"), + ("t2w", "inputnode.t2w"), + ("roi", "inputnode.roi"), + ("flair", "inputnode.flair"), + ("subjects_dir", "inputnode.subjects_dir"), + ("subject_id", "inputnode.subject_id"), + ]), + (anat_fit_wf, outputnode, [ + ("outputnode.template", "template"), + ("outputnode.subjects_dir", "subjects_dir"), + ("outputnode.subject_id", "subject_id"), + ("outputnode.t1w_preproc", "t1w_preproc"), + ("outputnode.t1w_mask", "t1w_mask"), + ("outputnode.t1w_dseg", "t1w_dseg"), + ("outputnode.t1w_tpms", "t1w_tpms"), + ("outputnode.anat2std_xfm", "anat2std_xfm"), + ("outputnode.std2anat_xfm", "std2anat_xfm"), + ("outputnode.fsnative2t1w_xfm", "fsnative2t1w_xfm"), + ("outputnode.sphere_reg", "sphere_reg"), + (f"outputnode.sphere_reg_{'msm' if msm_sulc else 'fsLR'}", "sphere_reg_fsLR"), + ("outputnode.anat_ribbon", "anat_ribbon"), + ]), + (anat_fit_wf, template_iterator_wf, [ + ("outputnode.template", "inputnode.template"), + ("outputnode.anat2std_xfm", "inputnode.anat2std_xfm"), + ]), + (anat_fit_wf, ds_std_volumes_wf, [ + ('outputnode.t1w_valid_list', 'inputnode.source_files'), + ("outputnode.t1w_preproc", "inputnode.t1w_preproc"), + ("outputnode.t1w_mask", "inputnode.t1w_mask"), + ("outputnode.t1w_dseg", "inputnode.t1w_dseg"), + ("outputnode.t1w_tpms", "inputnode.t1w_tpms"), + ]), + (template_iterator_wf, ds_std_volumes_wf, [ + ("outputnode.std_t1w", "inputnode.ref_file"), + ("outputnode.anat2std_xfm", "inputnode.anat2std_xfm"), + ("outputnode.space", "inputnode.space"), + ("outputnode.cohort", "inputnode.cohort"), + ("outputnode.resolution", "inputnode.resolution"), + ]), + ]) # fmt:skip + + if freesurfer: + anat_second_derivatives_wf = init_anat_second_derivatives_wf( + bids_root=bids_root, + output_dir=output_dir, + cifti_output=cifti_output, + ) + surface_derivatives_wf = init_surface_derivatives_wf( + cifti_output=cifti_output, + ) + ds_surfaces_wf = init_ds_surfaces_wf( + bids_root=bids_root, output_dir=output_dir, surfaces=["inflated"] + ) + ds_curv_wf = init_ds_surface_metrics_wf( + bids_root=bids_root, output_dir=output_dir, metrics=["curv"], name="ds_curv_wf" + ) + + workflow.connect([ + (anat_fit_wf, surface_derivatives_wf, [ + ('outputnode.t1w_preproc', 'inputnode.reference'), + ('outputnode.subjects_dir', 'inputnode.subjects_dir'), + ('outputnode.subject_id', 'inputnode.subject_id'), + ('outputnode.fsnative2t1w_xfm', 'inputnode.fsnative2t1w_xfm'), + ]), + (anat_fit_wf, ds_surfaces_wf, [ + ('outputnode.t1w_valid_list', 'inputnode.source_files'), + ]), + (surface_derivatives_wf, ds_surfaces_wf, [ + ('outputnode.inflated', 'inputnode.inflated'), + ]), + (anat_fit_wf, ds_curv_wf, [ + ('outputnode.t1w_valid_list', 'inputnode.source_files'), + ]), + (surface_derivatives_wf, ds_curv_wf, [ + ('outputnode.curv', 'inputnode.curv'), + ]), + (anat_fit_wf, anat_second_derivatives_wf, [ + ('outputnode.t1w_valid_list', 'inputnode.source_files'), + ]), + (surface_derivatives_wf, anat_second_derivatives_wf, [ + ('outputnode.out_aseg', 'inputnode.t1w_fs_aseg'), + ('outputnode.out_aparc', 'inputnode.t1w_fs_aparc'), + ]), + (surface_derivatives_wf, outputnode, [ + ('outputnode.out_aseg', 't1w_aseg'), + ('outputnode.out_aparc', 't1w_aparc'), + ]), + ]) # fmt:skip + + if cifti_output: + hcp_morphometrics_wf = init_hcp_morphometrics_wf(omp_nthreads=omp_nthreads) + resample_midthickness_wf = init_resample_midthickness_wf(grayord_density=cifti_output) + morph_grayords_wf = init_morph_grayords_wf( + grayord_density=cifti_output, omp_nthreads=omp_nthreads + ) + + ds_grayord_metrics_wf = init_ds_grayord_metrics_wf( + bids_root=bids_root, + output_dir=output_dir, + metrics=['curv', 'thickness', 'sulc'], + cifti_output=cifti_output, + ) + + workflow.connect([ + (anat_fit_wf, hcp_morphometrics_wf, [ + ('outputnode.subject_id', 'inputnode.subject_id'), + ('outputnode.sulc', 'inputnode.sulc'), + ('outputnode.thickness', 'inputnode.thickness'), + ('outputnode.midthickness', 'inputnode.midthickness'), + ]), + (surface_derivatives_wf, hcp_morphometrics_wf, [ + ('outputnode.curv', 'inputnode.curv'), + ]), + (anat_fit_wf, resample_midthickness_wf, [ + ('outputnode.midthickness', 'inputnode.midthickness'), + ( + f"outputnode.sphere_reg_{'msm' if msm_sulc else 'fsLR'}", + "inputnode.sphere_reg_fsLR", + ), + ]), + (anat_fit_wf, morph_grayords_wf, [ + ('outputnode.midthickness', 'inputnode.midthickness'), + ( + f"outputnode.sphere_reg_{'msm' if msm_sulc else 'fsLR'}", + "inputnode.sphere_reg_fsLR", + ), + ]), + (hcp_morphometrics_wf, morph_grayords_wf, [ + ('outputnode.curv', 'inputnode.curv'), + ('outputnode.sulc', 'inputnode.sulc'), + ('outputnode.thickness', 'inputnode.thickness'), + ('outputnode.roi', 'inputnode.roi'), + ]), + (resample_midthickness_wf, morph_grayords_wf, [ + ('outputnode.midthickness_fsLR', 'inputnode.midthickness_fsLR'), + ]), + (anat_fit_wf, ds_grayord_metrics_wf, [ + ('outputnode.t1w_valid_list', 'inputnode.source_files'), + ]), + (morph_grayords_wf, ds_grayord_metrics_wf, [ + ('outputnode.curv_fsLR', 'inputnode.curv'), + ('outputnode.curv_metadata', 'inputnode.curv_metadata'), + ('outputnode.thickness_fsLR', 'inputnode.thickness'), + ('outputnode.thickness_metadata', 'inputnode.thickness_metadata'), + ('outputnode.sulc_fsLR', 'inputnode.sulc'), + ('outputnode.sulc_metadata', 'inputnode.sulc_metadata'), + ]), + ]) # fmt:skip + + return workflow + + +def init_anat_fit_wf( + *, + bids_root: str, + output_dir: str, + freesurfer: bool, + hires: bool, + longitudinal: bool, + msm_sulc: bool, + t1w: list, + t2w: list, + skull_strip_mode: str, + skull_strip_template: Reference, + spaces: SpatialReferences, + precomputed: dict, + omp_nthreads: int, + debug: bool = False, + sloppy: bool = False, + name="anat_fit_wf", + skull_strip_fixed_seed: bool = False, ): """ Stage the anatomical preprocessing steps of *sMRIPrep*. @@ -101,19 +470,23 @@ def init_anat_preproc_wf( :simple_form: yes from niworkflows.utils.spaces import SpatialReferences, Reference - from smriprep.workflows.anatomical import init_anat_preproc_wf - wf = init_anat_preproc_wf( + from smriprep.workflows.anatomical import init_anat_fit_wf + wf = init_anat_fit_wf( bids_root='.', + output_dir='.', freesurfer=True, hires=True, longitudinal=False, + msm_sulc=True, t1w=['t1w.nii.gz'], - t2w=[], - omp_nthreads=1, - output_dir='.', + t2w=['t2w.nii.gz'], skull_strip_mode='force', skull_strip_template=Reference('OASIS30ANTs'), spaces=SpatialReferences(spaces=['MNI152NLin2009cAsym', 'fsaverage5']), + precomputed={}, + debug=False, + sloppy=False, + omp_nthreads=1, ) @@ -121,9 +494,8 @@ def init_anat_preproc_wf( ---------- bids_root : :obj:`str` Path of the input BIDS dataset root - existing_derivatives : :obj:`dict` or None - Dictionary mapping output specification attribute names and - paths to corresponding derivatives. + output_dir : :obj:`str` + Directory in which to save derivatives freesurfer : :obj:`bool` Enable FreeSurfer surface reconstruction (increases runtime by 6h, at the very least) @@ -134,24 +506,25 @@ def init_anat_preproc_wf( (may increase runtime) t1w : :obj:`list` List of T1-weighted structural images. - omp_nthreads : :obj:`int` - Maximum number of threads an individual process may use - output_dir : :obj:`str` - Directory in which to save derivatives + skull_strip_mode : :obj:`str` + Determiner for T1-weighted skull stripping (`force` ensures skull stripping, + `skip` ignores skull stripping, and `auto` automatically ignores skull stripping + if pre-stripped brains are detected). skull_strip_template : :py:class:`~niworkflows.utils.spaces.Reference` Spatial reference to use in atlas-based brain extraction. spaces : :py:class:`~niworkflows.utils.spaces.SpatialReferences` Object containing standard and nonstandard space specifications. + precomputed : :obj:`dict` + Dictionary mapping output specification attribute names and + paths to precomputed derivatives. + omp_nthreads : :obj:`int` + Maximum number of threads an individual process may use debug : :obj:`bool` Enable debugging outputs sloppy: :obj:`bool` Quick, impercise operations. Used to decrease workflow duration. name : :obj:`str`, optional - Workflow name (default: anat_preproc_wf) - skull_strip_mode : :obj:`str` - Determiner for T1-weighted skull stripping (`force` ensures skull stripping, - `skip` ignores skull stripping, and `auto` automatically ignores skull stripping - if pre-stripped brains are detected). + Workflow name (default: anat_fit_wf) skull_strip_fixed_seed : :obj:`bool` Do not use a random seed for skull-stripping - will ensure run-to-run replicability when used with --omp-nthreads 1 @@ -177,8 +550,6 @@ def init_anat_preproc_wf( t1w_preproc The T1w reference map, which is calculated as the average of bias-corrected and preprocessed T1w images, defining the anatomical space. - t1w_brain - Skull-stripped ``t1w_preproc`` t1w_mask Brain (binary) mask estimated by brain extraction. t1w_dseg @@ -186,33 +557,25 @@ def init_anat_preproc_wf( gray-matter (GM), white-matter (WM) and cerebrospinal fluid (CSF). t1w_tpms List of tissue probability maps corresponding to ``t1w_dseg``. - std_preproc - T1w reference resampled in one or more standard spaces. - std_mask - Mask of skull-stripped template, in MNI space - std_dseg - Segmentation, resampled into MNI space - std_tpms - List of tissue probability maps in MNI space - subjects_dir - FreeSurfer SUBJECTS_DIR + t1w_valid_list + List of input T1w images accepted for preprocessing. If t1w_preproc is + precomputed, this is always a list containing that image. + template + List of template names to which the structural image has been registered anat2std_xfm - Nonlinear spatial transform to resample imaging data given in anatomical space - into standard space. + List of nonlinear spatial transforms to resample data from subject + anatomical space into standard template spaces. Collated with template. std2anat_xfm - Inverse transform of the above. + List of nonlinear spatial transforms to resample data from standard + template spaces into subject anatomical space. Collated with template. + subjects_dir + FreeSurfer SUBJECTS_DIR; use as input to a node to ensure that it is run after + FreeSurfer reconstruction is completed. subject_id - FreeSurfer subject ID - t1w2fsnative_xfm - LTA-style affine matrix translating from T1w to - FreeSurfer-conformed subject space + FreeSurfer subject ID; use as input to a node to ensure that it is run after + FreeSurfer reconstruction is completed. fsnative2t1w_xfm - LTA-style affine matrix translating from FreeSurfer-conformed - subject space to T1w - surfaces - GIFTI surfaces (gray/white boundary, midthickness, pial, inflated) - morphometrics - GIFTIs of cortical thickness, curvature, and sulcal depth + ITK-style affine matrix translating from FreeSurfer-conformed subject space to T1w See Also -------- @@ -222,317 +585,528 @@ def init_anat_preproc_wf( """ workflow = Workflow(name=name) num_t1w = len(t1w) - desc = """ + desc = f""" Anatomical data preprocessing -: """ - desc += """\ -A total of {num_t1w} T1-weighted (T1w) images were found within the input -BIDS dataset.""".format( - num_t1w=num_t1w - ) - +: A total of {num_t1w} T1-weighted (T1w) images were found within the input +BIDS dataset.""" + + have_t1w = "t1w_preproc" in precomputed + have_t2w = "t2w_preproc" in precomputed + have_mask = "t1w_mask" in precomputed + have_dseg = "t1w_dseg" in precomputed + have_tpms = "t1w_tpms" in precomputed + + # Organization + # ------------ + # This workflow takes the usual (inputnode -> graph -> outputnode) format + # The graph consists of (input -> compute -> datasink -> buffer) units, + # and all inputs to outputnode are buffer. + # If precomputed inputs are found, then these units are replaced with (buffer) + # At the time of writing, t1w_mask is an exception, which takes the form + # (t1w_buffer -> refined_buffer -> datasink -> outputnode) + # All outputnode components should therefore point to files in the input or + # output directories. inputnode = pe.Node( niu.IdentityInterface(fields=["t1w", "t2w", "roi", "flair", "subjects_dir", "subject_id"]), name="inputnode", ) - outputnode = pe.Node( niu.IdentityInterface( - fields=["template", "subjects_dir", "subject_id", "t2w_preproc", "sphere_reg_fsLR"] - + get_outputnode_spec() + fields=[ + # Primary derivatives + "t1w_preproc", + "t2w_preproc", + "t1w_mask", + "t1w_dseg", + "t1w_tpms", + "anat2std_xfm", + "fsnative2t1w_xfm", + # Surface and metric derivatives for fsLR resampling + "white", + "pial", + "midthickness", + "sphere", + "thickness", + "sulc", + "sphere_reg", + "sphere_reg_fsLR", + "sphere_reg_msm", + "anat_ribbon", + # Reverse transform; not computable from forward transform + "std2anat_xfm", + # Metadata + "template", + "subjects_dir", + "subject_id", + "t1w_valid_list", + ] ), name="outputnode", ) + # If all derivatives exist, inputnode could go unconnected, so add explicitly + workflow.add_nodes([inputnode]) + + # Stage 1 inputs (filtered) + sourcefile_buffer = pe.Node( + niu.IdentityInterface(fields=["source_files"]), + name="sourcefile_buffer", + ) + + # Stage 2 results + t1w_buffer = pe.Node( + niu.IdentityInterface(fields=["t1w_preproc", "t1w_mask", "t1w_brain", "ants_seg"]), + name="t1w_buffer", + ) + # Stage 3 results + seg_buffer = pe.Node( + niu.IdentityInterface(fields=["t1w_dseg", "t1w_tpms"]), + name="seg_buffer", + ) + # Stage 4 results: collated template names, forward and reverse transforms + template_buffer = pe.Node(niu.Merge(2), name="template_buffer") + anat2std_buffer = pe.Node(niu.Merge(2), name="anat2std_buffer") + std2anat_buffer = pe.Node(niu.Merge(2), name="std2anat_buffer") + + # Stage 6 results: Refined stage 2 results; may be direct copy if no refinement + refined_buffer = pe.Node( + niu.IdentityInterface(fields=["t1w_mask", "t1w_brain"]), + name="refined_buffer", + ) + + # Stage 8 results: GIFTI surfaces + surfaces_buffer = pe.Node( + niu.IdentityInterface( + fields=["white", "pial", "midthickness", "sphere", "sphere_reg", "thickness", "sulc"] + ), + name="surfaces_buffer", + ) + + # Stage 9 and 10 results: fsLR sphere registration + fsLR_buffer = pe.Node(niu.IdentityInterface(fields=["sphere_reg_fsLR"]), name="fsLR_buffer") + msm_buffer = pe.Node(niu.IdentityInterface(fields=["sphere_reg_msm"]), name="msm_buffer") - # Connect reportlets workflows + # fmt:off + workflow.connect([ + (seg_buffer, outputnode, [ + ("t1w_dseg", "t1w_dseg"), + ("t1w_tpms", "t1w_tpms"), + ]), + (anat2std_buffer, outputnode, [("out", "anat2std_xfm")]), + (std2anat_buffer, outputnode, [("out", "std2anat_xfm")]), + (template_buffer, outputnode, [("out", "template")]), + (sourcefile_buffer, outputnode, [("source_files", "t1w_valid_list")]), + (surfaces_buffer, outputnode, [ + ("white", "white"), + ("pial", "pial"), + ("midthickness", "midthickness"), + ("sphere", "sphere"), + ("sphere_reg", "sphere_reg"), + ("thickness", "thickness"), + ("sulc", "sulc"), + ]), + (fsLR_buffer, outputnode, [("sphere_reg_fsLR", "sphere_reg_fsLR")]), + (msm_buffer, outputnode, [("sphere_reg_msm", "sphere_reg_msm")]), + ]) + # fmt:on + + # Reporting anat_reports_wf = init_anat_reports_wf( + spaces=spaces, freesurfer=freesurfer, output_dir=output_dir, ) # fmt:off workflow.connect([ (outputnode, anat_reports_wf, [ - ('t1w_preproc', 'inputnode.t1w_preproc'), - ('t1w_mask', 'inputnode.t1w_mask'), - ('t1w_dseg', 'inputnode.t1w_dseg')]), + ("t1w_valid_list", "inputnode.source_file"), + ("t1w_preproc", "inputnode.t1w_preproc"), + ("t1w_mask", "inputnode.t1w_mask"), + ("t1w_dseg", "inputnode.t1w_dseg"), + ("template", "inputnode.template"), + ("anat2std_xfm", "inputnode.anat2std_xfm"), + ("subjects_dir", "inputnode.subjects_dir"), + ("subject_id", "inputnode.subject_id"), + ]), ]) # fmt:on - if existing_derivatives is not None: - LOGGER.log( - 25, - "Anatomical workflow will reuse prior derivatives found in the output folder (%s).", - output_dir, - ) - desc += """ -Anatomical preprocessing was reused from previously existing derivative objects.\n""" - workflow.__desc__ = desc - - templates = existing_derivatives.pop("template") - templatesource = pe.Node(niu.IdentityInterface(fields=["template"]), name="templatesource") - templatesource.iterables = [("template", templates)] - outputnode.inputs.template = templates - - for field, value in existing_derivatives.items(): - setattr(outputnode.inputs, field, value) - - anat_reports_wf.inputs.inputnode.source_file = [existing_derivatives["t1w_preproc"]] + # Stage 1: Conform images and validate + # If desc-preproc_T1w.nii.gz is provided, just validate it + anat_validate = pe.Node(ValidateImage(), name="anat_validate", run_without_submitting=True) + if not have_t1w: + LOGGER.info("ANAT Stage 1: Adding template workflow") + ants_ver = ANTsInfo.version() or "(version unknown)" + desc += f"""\ + {"Each" if num_t1w > 1 else "The"} T1w image was corrected for intensity +non-uniformity (INU) with `N4BiasFieldCorrection` [@n4], distributed with ANTs {ants_ver} +[@ants, RRID:SCR_004757]""" + desc += ".\n" if num_t1w > 1 else ", and used as T1w-reference throughout the workflow.\n" - stdselect = pe.Node( - KeySelect(fields=["std_preproc", "std_mask"], keys=templates), - name="stdselect", - run_without_submitting=True, + anat_template_wf = init_anat_template_wf( + longitudinal=longitudinal, + omp_nthreads=omp_nthreads, + num_files=num_t1w, + contrast="T1w", + name="anat_template_wf", ) + ds_template_wf = init_ds_template_wf(output_dir=output_dir, num_t1w=num_t1w) + # fmt:off workflow.connect([ - (inputnode, outputnode, [('subjects_dir', 'subjects_dir'), - ('subject_id', 'subject_id')]), - (inputnode, anat_reports_wf, [ - ('subjects_dir', 'inputnode.subjects_dir'), - ('subject_id', 'inputnode.subject_id')]), - (templatesource, stdselect, [('template', 'key')]), - (outputnode, stdselect, [('std_preproc', 'std_preproc'), - ('std_mask', 'std_mask')]), - (stdselect, anat_reports_wf, [ - ('key', 'inputnode.template'), - ('std_preproc', 'inputnode.std_t1w'), - ('std_mask', 'inputnode.std_mask'), + (inputnode, anat_template_wf, [("t1w", "inputnode.anat_files")]), + (anat_template_wf, anat_validate, [("outputnode.anat_ref", "in_file")]), + (anat_template_wf, sourcefile_buffer, [ + ("outputnode.anat_valid_list", "source_files"), + ]), + (anat_template_wf, anat_reports_wf, [ + ("outputnode.out_report", "inputnode.t1w_conform_report"), ]), + (anat_template_wf, ds_template_wf, [ + ("outputnode.anat_realign_xfm", "inputnode.t1w_ref_xfms"), + ]), + (sourcefile_buffer, ds_template_wf, [("source_files", "inputnode.source_files")]), + (t1w_buffer, ds_template_wf, [("t1w_preproc", "inputnode.t1w_preproc")]), + (ds_template_wf, outputnode, [("outputnode.t1w_preproc", "t1w_preproc")]), ]) # fmt:on - return workflow - - # The workflow is not cached. - desc += ( - """ -All of them were corrected for intensity non-uniformity (INU) -""" - if num_t1w > 1 - else """\ -The T1-weighted (T1w) image was corrected for intensity non-uniformity (INU) + else: + LOGGER.info("ANAT Found preprocessed T1w - skipping Stage 1") + desc += """ A preprocessed T1w image was provided as a precomputed input +and used as T1w-reference throughout the workflow. """ - ) - desc += """\ -with `N4BiasFieldCorrection` [@n4], distributed with ANTs {ants_ver} \ -[@ants, RRID:SCR_004757]""" - desc += ".\n" if num_t1w > 1 else ", and used as T1w-reference throughout the workflow.\n" - desc += """\ + anat_validate.inputs.in_file = precomputed["t1w_preproc"] + sourcefile_buffer.inputs.source_files = [precomputed["t1w_preproc"]] + + # fmt:off + workflow.connect([ + (anat_validate, t1w_buffer, [("out_file", "t1w_preproc")]), + (t1w_buffer, outputnode, [("t1w_preproc", "t1w_preproc")]), + ]) + # fmt:on + + # Stage 2: INU correction and masking + # We always need to generate t1w_brain; how to do that depends on whether we have + # a pre-corrected T1w or precomputed mask, or are given an already masked image + if not have_mask: + LOGGER.info("ANAT Stage 2: Preparing brain extraction workflow") + if skull_strip_mode == "auto": + run_skull_strip = all(_is_skull_stripped(img) for img in t1w) + else: + run_skull_strip = {"force": True, "skip": False} + + # Brain extraction + if run_skull_strip: + desc += f"""\ The T1w-reference was then skull-stripped with a *Nipype* implementation of -the `antsBrainExtraction.sh` workflow (from ANTs), using {skullstrip_tpl} +the `antsBrainExtraction.sh` workflow (from ANTs), using {skull_strip_template.fullname} as target template. +""" + brain_extraction_wf = init_brain_extraction_wf( + in_template=skull_strip_template.space, + template_spec=skull_strip_template.spec, + atropos_use_random_seed=not skull_strip_fixed_seed, + omp_nthreads=omp_nthreads, + normalization_quality="precise" if not sloppy else "testing", + ) + # fmt:off + workflow.connect([ + (anat_validate, brain_extraction_wf, [("out_file", "inputnode.in_files")]), + (brain_extraction_wf, t1w_buffer, [ + ("outputnode.out_mask", "t1w_mask"), + (("outputnode.out_file", _pop), "t1w_brain"), + ("outputnode.out_segm", "ants_seg"), + ]), + ]) + if not have_t1w: + workflow.connect([ + (brain_extraction_wf, t1w_buffer, [ + (("outputnode.bias_corrected", _pop), "t1w_preproc"), + ]), + ]) + # fmt:on + # Determine mask from T1w and uniformize + elif not have_t1w: + LOGGER.info("ANAT Stage 2: Skipping skull-strip, INU-correction only") + desc += """\ +The provided T1w image was previously skull-stripped; a brain mask was +derived from the input image. +""" + n4_only_wf = init_n4_only_wf( + omp_nthreads=omp_nthreads, + atropos_use_random_seed=not skull_strip_fixed_seed, + ) + # fmt:off + workflow.connect([ + (anat_validate, n4_only_wf, [("out_file", "inputnode.in_files")]), + (n4_only_wf, t1w_buffer, [ + (("outputnode.bias_corrected", _pop), "t1w_preproc"), + ("outputnode.out_mask", "t1w_mask"), + (("outputnode.out_file", _pop), "t1w_brain"), + ("outputnode.out_segm", "ants_seg"), + ]), + ]) + # fmt:on + # Binarize the already uniformized image + else: + LOGGER.info("ANAT Stage 2: Skipping skull-strip, generating mask from input") + desc += """\ +The provided T1w image was previously skull-stripped; a brain mask was +derived from the input image. +""" + binarize = pe.Node(Binarize(thresh_low=2), name="binarize") + # fmt:off + workflow.connect([ + (anat_validate, binarize, [("out_file", "in_file")]), + (anat_validate, t1w_buffer, [("out_file", "t1w_brain")]), + (binarize, t1w_buffer, [("out_file", "t1w_mask")]), + ]) + # fmt:on + + ds_t1w_mask_wf = init_ds_mask_wf( + bids_root=bids_root, + output_dir=output_dir, + mask_type="brain", + name="ds_t1w_mask_wf", + ) + # fmt:off + workflow.connect([ + (sourcefile_buffer, ds_t1w_mask_wf, [("source_files", "inputnode.source_files")]), + (refined_buffer, ds_t1w_mask_wf, [("t1w_mask", "inputnode.mask_file")]), + (ds_t1w_mask_wf, outputnode, [("outputnode.mask_file", "t1w_mask")]), + ]) + # fmt:on + else: + LOGGER.info("ANAT Found brain mask") + desc += """\ +A pre-computed brain mask was provided as input and used throughout the workflow. +""" + t1w_buffer.inputs.t1w_mask = precomputed["t1w_mask"] + # If we have a mask, always apply it + apply_mask = pe.Node(ApplyMask(in_mask=precomputed["t1w_mask"]), name="apply_mask") + workflow.connect([(anat_validate, apply_mask, [("out_file", "in_file")])]) + # Run N4 if it hasn't been pre-run + if not have_t1w: + LOGGER.info("ANAT Skipping skull-strip, INU-correction only") + n4_only_wf = init_n4_only_wf( + omp_nthreads=omp_nthreads, + atropos_use_random_seed=not skull_strip_fixed_seed, + ) + # fmt:off + workflow.connect([ + (apply_mask, n4_only_wf, [("out_file", "inputnode.in_files")]), + (n4_only_wf, t1w_buffer, [ + (("outputnode.bias_corrected", _pop), "t1w_preproc"), + (("outputnode.out_file", _pop), "t1w_brain"), + ]), + ]) + # fmt:on + else: + LOGGER.info("ANAT Skipping Stage 2") + workflow.connect([(apply_mask, t1w_buffer, [("out_file", "t1w_brain")])]) + workflow.connect([(refined_buffer, outputnode, [("t1w_mask", "t1w_mask")])]) + + # Stage 3: Segmentation + if not (have_dseg and have_tpms): + LOGGER.info("ANAT Stage 3: Preparing segmentation workflow") + fsl_ver = fsl.FAST().version or "(version unknown)" + desc += f"""\ Brain tissue segmentation of cerebrospinal fluid (CSF), white-matter (WM) and gray-matter (GM) was performed on -the brain-extracted T1w using `fast` [FSL {fsl_ver}, RRID:SCR_002823, -@fsl_fast]. +the brain-extracted T1w using `fast` [FSL {fsl_ver}, RRID:SCR_002823, @fsl_fast]. """ - - workflow.__desc__ = desc.format( - ants_ver=ANTsInfo.version() or "(version unknown)", - fsl_ver=fsl.FAST().version or "(version unknown)", - num_t1w=num_t1w, - skullstrip_tpl=skull_strip_template.fullname, - ) - - buffernode = pe.Node( - niu.IdentityInterface(fields=["t1w_brain", "t1w_mask"]), name="buffernode" - ) - - # 1. Anatomical reference generation - average input T1w images. - anat_template_wf = init_anat_template_wf( - longitudinal=longitudinal, - omp_nthreads=omp_nthreads, - num_files=num_t1w, - contrast="T1w", - ) - - anat_validate = pe.Node(ValidateImage(), name="anat_validate", run_without_submitting=True) - - # 2. Brain-extraction and INU (bias field) correction. - if skull_strip_mode == "auto": - import numpy as np - import nibabel as nb - - def _is_skull_stripped(imgs): - """Check if T1w images are skull-stripped.""" - - def _check_img(img): - data = np.abs(nb.load(img).get_fdata(dtype=np.float32)) - sidevals = ( - data[0, :, :].sum() - + data[-1, :, :].sum() - + data[:, 0, :].sum() - + data[:, -1, :].sum() - + data[:, :, 0].sum() - + data[:, :, -1].sum() - ) - return sidevals < 10 - - return all(_check_img(img) for img in imgs) - - skull_strip_mode = _is_skull_stripped(t1w) - - if skull_strip_mode in (True, "skip"): - brain_extraction_wf = init_n4_only_wf( - omp_nthreads=omp_nthreads, - atropos_use_random_seed=not skull_strip_fixed_seed, + fast = pe.Node( + fsl.FAST(segments=True, no_bias=True, probability_maps=True), + name="fast", + mem_gb=3, + ) + lut_t1w_dseg = pe.Node(niu.Function(function=_apply_bids_lut), name="lut_t1w_dseg") + lut_t1w_dseg.inputs.lut = (0, 3, 1, 2) # Maps: 0 -> 0, 3 -> 1, 1 -> 2, 2 -> 3. + fast2bids = pe.Node( + niu.Function(function=_probseg_fast2bids), + name="fast2bids", + run_without_submitting=True, ) + workflow.connect([(refined_buffer, fast, [("t1w_brain", "in_files")])]) + + # fmt:off + if not have_dseg: + ds_dseg_wf = init_ds_dseg_wf(output_dir=output_dir) + workflow.connect([ + (fast, lut_t1w_dseg, [("partial_volume_map", "in_dseg")]), + (sourcefile_buffer, ds_dseg_wf, [("source_files", "inputnode.source_files")]), + (lut_t1w_dseg, ds_dseg_wf, [("out", "inputnode.t1w_dseg")]), + (ds_dseg_wf, seg_buffer, [("outputnode.t1w_dseg", "t1w_dseg")]), + ]) + if not have_tpms: + ds_tpms_wf = init_ds_tpms_wf(output_dir=output_dir) + workflow.connect([ + (fast, fast2bids, [("partial_volume_files", "inlist")]), + (sourcefile_buffer, ds_tpms_wf, [("source_files", "inputnode.source_files")]), + (fast2bids, ds_tpms_wf, [("out", "inputnode.t1w_tpms")]), + (ds_tpms_wf, seg_buffer, [("outputnode.t1w_tpms", "t1w_tpms")]), + ]) + # fmt:on else: - brain_extraction_wf = init_brain_extraction_wf( - in_template=skull_strip_template.space, - template_spec=skull_strip_template.spec, - atropos_use_random_seed=not skull_strip_fixed_seed, + LOGGER.info("ANAT Skipping Stage 3") + if have_dseg: + LOGGER.info("ANAT Found discrete segmentation") + desc += "Precomputed discrete tissue segmentations were provided as inputs.\n" + seg_buffer.inputs.t1w_dseg = precomputed["t1w_dseg"] + if have_tpms: + LOGGER.info("ANAT Found tissue probability maps") + desc += "Precomputed tissue probabiilty maps were provided as inputs.\n" + seg_buffer.inputs.t1w_tpms = precomputed["t1w_tpms"] + + # Stage 4: Normalization + templates = [] + found_xfms = {} + for template in spaces.get_spaces(nonstandard=False, dim=(3,)): + xfms = precomputed.get("transforms", {}).get(template, {}) + if set(xfms) != {"forward", "reverse"}: + templates.append(template) + else: + found_xfms[template] = xfms + + template_buffer.inputs.in1 = list(found_xfms) + anat2std_buffer.inputs.in1 = [xfm["forward"] for xfm in found_xfms.values()] + std2anat_buffer.inputs.in1 = [xfm["reverse"] for xfm in found_xfms.values()] + + if templates: + LOGGER.info(f"ANAT Stage 4: Preparing normalization workflow for {templates}") + register_template_wf = init_register_template_wf( + sloppy=sloppy, omp_nthreads=omp_nthreads, - normalization_quality="precise" if not sloppy else "testing", + templates=templates, ) + ds_template_registration_wf = init_ds_template_registration_wf(output_dir=output_dir) - # 4. Spatial normalization - anat_norm_wf = init_anat_norm_wf( - sloppy=sloppy, - omp_nthreads=omp_nthreads, - templates=spaces.get_spaces(nonstandard=False, dim=(3,)), - ) + # fmt:off + workflow.connect([ + (inputnode, register_template_wf, [('roi', 'inputnode.lesion_mask')]), + (t1w_buffer, register_template_wf, [('t1w_preproc', 'inputnode.moving_image')]), + (refined_buffer, register_template_wf, [('t1w_mask', 'inputnode.moving_mask')]), + (sourcefile_buffer, ds_template_registration_wf, [ + ("source_files", "inputnode.source_files") + ]), + (register_template_wf, ds_template_registration_wf, [ + ('outputnode.template', 'inputnode.template'), + ('outputnode.anat2std_xfm', 'inputnode.anat2std_xfm'), + ('outputnode.std2anat_xfm', 'inputnode.std2anat_xfm'), + ]), + (register_template_wf, template_buffer, [('outputnode.template', 'in2')]), + (ds_template_registration_wf, std2anat_buffer, [('outputnode.std2anat_xfm', 'in2')]), + (ds_template_registration_wf, anat2std_buffer, [('outputnode.anat2std_xfm', 'in2')]), + ]) + # fmt:on + if found_xfms: + LOGGER.info(f"ANAT Stage 4: Found pre-computed registrations for {found_xfms}") - # fmt:off - workflow.connect([ - # Step 1. - (inputnode, anat_template_wf, [('t1w', 'inputnode.anat_files')]), - (anat_template_wf, anat_validate, [ - ('outputnode.anat_ref', 'in_file')]), - (anat_validate, brain_extraction_wf, [ - ('out_file', 'inputnode.in_files')]), - (brain_extraction_wf, outputnode, [ - (('outputnode.bias_corrected', _pop), 't1w_preproc')]), - (anat_template_wf, outputnode, [ - ('outputnode.anat_realign_xfm', 't1w_ref_xfms')]), - (buffernode, outputnode, [('t1w_brain', 't1w_brain'), - ('t1w_mask', 't1w_mask')]), - # Steps 2, 3 and 4 - (inputnode, anat_norm_wf, [ - (('t1w', fix_multi_T1w_source_name), 'inputnode.orig_t1w'), - ('roi', 'inputnode.lesion_mask')]), - (brain_extraction_wf, anat_norm_wf, [ - (('outputnode.bias_corrected', _pop), 'inputnode.moving_image')]), - (buffernode, anat_norm_wf, [('t1w_mask', 'inputnode.moving_mask')]), - (anat_norm_wf, outputnode, [ - ('poutputnode.standardized', 'std_preproc'), - ('poutputnode.std_mask', 'std_mask'), - ('poutputnode.std_dseg', 'std_dseg'), - ('poutputnode.std_tpms', 'std_tpms'), - ('outputnode.template', 'template'), - ('outputnode.anat2std_xfm', 'anat2std_xfm'), - ('outputnode.std2anat_xfm', 'std2anat_xfm'), - ]), - ]) - # fmt:on + # Do not attempt refinement (Stage 6, below) + if have_mask or not freesurfer: + # fmt:off + workflow.connect([ + (t1w_buffer, refined_buffer, [ + ("t1w_mask", "t1w_mask"), + ("t1w_brain", "t1w_brain"), + ]), + ]) + # fmt:on - # Change LookUp Table - BIDS wants: 0 (bg), 1 (gm), 2 (wm), 3 (csf) - lut_t1w_dseg = pe.Node(niu.Function(function=_apply_bids_lut), name="lut_t1w_dseg") + workflow.__desc__ = desc - # fmt:off - workflow.connect([ - (lut_t1w_dseg, anat_norm_wf, [ - ('out', 'inputnode.moving_segmentation')]), - (lut_t1w_dseg, outputnode, [('out', 't1w_dseg')]), - ]) - # fmt:on + if not freesurfer: + LOGGER.info("ANAT Skipping Stages 5+") + return workflow - # Connect reportlets - # fmt:off - workflow.connect([ - (inputnode, anat_reports_wf, [('t1w', 'inputnode.source_file')]), - (outputnode, anat_reports_wf, [ - ('std_preproc', 'inputnode.std_t1w'), - ('std_mask', 'inputnode.std_mask'), - ]), - (anat_template_wf, anat_reports_wf, [ - ('outputnode.out_report', 'inputnode.t1w_conform_report')]), - (anat_norm_wf, anat_reports_wf, [ - ('poutputnode.template', 'inputnode.template')]), - ]) - # fmt:on + fs_isrunning = pe.Node( + niu.Function(function=_fs_isRunning), overwrite=True, name="fs_isrunning" + ) + fs_isrunning.inputs.logger = LOGGER - # Write outputs ############################################3 - anat_derivatives_wf = init_anat_derivatives_wf( - bids_root=bids_root, - freesurfer=freesurfer, - num_t1w=num_t1w, - t2w=t2w, - output_dir=output_dir, - spaces=spaces, - cifti_output=cifti_output, + # Stage 5: Surface reconstruction (--fs-no-reconall not set) + LOGGER.info("ANAT Stage 5: Preparing surface reconstruction workflow") + surface_recon_wf = init_surface_recon_wf( + name="surface_recon_wf", + omp_nthreads=omp_nthreads, + hires=hires, + precomputed=precomputed, ) # fmt:off workflow.connect([ - # Connect derivatives - (anat_template_wf, anat_derivatives_wf, [ - ('outputnode.anat_valid_list', 'inputnode.source_files')]), - (anat_norm_wf, anat_derivatives_wf, [ - ('outputnode.template', 'inputnode.template'), - ('outputnode.anat2std_xfm', 'inputnode.anat2std_xfm'), - ('outputnode.std2anat_xfm', 'inputnode.std2anat_xfm') + (inputnode, fs_isrunning, [ + ('subjects_dir', 'subjects_dir'), + ('subject_id', 'subject_id'), ]), - (outputnode, anat_derivatives_wf, [ - ('t1w_ref_xfms', 'inputnode.t1w_ref_xfms'), - ('t1w_preproc', 'inputnode.t1w_preproc'), - ('t1w_mask', 'inputnode.t1w_mask'), - ('t1w_dseg', 'inputnode.t1w_dseg'), - ('t1w_tpms', 'inputnode.t1w_tpms'), - ('t2w_preproc', 'inputnode.t2w_preproc'), + (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')]), + (t1w_buffer, surface_recon_wf, [('t1w_brain', 'inputnode.skullstripped_t1')]), + (surface_recon_wf, outputnode, [ + ('outputnode.subjects_dir', 'subjects_dir'), + ('outputnode.subject_id', 'subject_id'), ]), ]) # fmt:on - # XXX Keeping FAST separate so that it's easier to swap in ANTs or FreeSurfer - - # Brain tissue segmentation - FAST produces: 0 (bg), 1 (wm), 2 (csf), 3 (gm) - t1w_dseg = pe.Node( - fsl.FAST(segments=True, no_bias=True, probability_maps=True), - name="t1w_dseg", - mem_gb=3, - ) - lut_t1w_dseg.inputs.lut = (0, 3, 1, 2) # Maps: 0 -> 0, 3 -> 1, 1 -> 2, 2 -> 3. - fast2bids = pe.Node( - niu.Function(function=_probseg_fast2bids), - name="fast2bids", - run_without_submitting=True, - ) - - # fmt:off - workflow.connect([ - (buffernode, t1w_dseg, [('t1w_brain', 'in_files')]), - (t1w_dseg, lut_t1w_dseg, [('partial_volume_map', 'in_dseg')]), - (t1w_dseg, fast2bids, [('partial_volume_files', 'inlist')]), - (fast2bids, anat_norm_wf, [('out', 'inputnode.moving_tpms')]), - (fast2bids, outputnode, [('out', 't1w_tpms')]), - ]) - # fmt:on - if not freesurfer: # Flag --fs-no-reconall is set - return + fsnative_xfms = precomputed.get("transforms", {}).get("fsnative") + if not fsnative_xfms: + ds_fs_registration_wf = init_ds_fs_registration_wf(output_dir=output_dir) # fmt:off workflow.connect([ - (brain_extraction_wf, buffernode, [ - (('outputnode.out_file', _pop), 't1w_brain'), - ('outputnode.out_mask', 't1w_mask')]), + (sourcefile_buffer, ds_fs_registration_wf, [ + ("source_files", "inputnode.source_files"), + ]), + (surface_recon_wf, ds_fs_registration_wf, [ + ('outputnode.fsnative2t1w_xfm', 'inputnode.fsnative2t1w_xfm'), + ]), + (ds_fs_registration_wf, outputnode, [ + ('outputnode.fsnative2t1w_xfm', 'fsnative2t1w_xfm'), + ]), ]) # fmt:on - return workflow + elif "reverse" in fsnative_xfms: + LOGGER.info("ANAT Found fsnative-T1w transform - skipping registration") + outputnode.inputs.fsnative2t1w_xfm = fsnative_xfms["reverse"] + else: + raise RuntimeError( + "Found a T1w-to-fsnative transform without the reverse. Time to handle this." + ) - # check for older IsRunning files and remove accordingly - fs_isrunning = pe.Node( - niu.Function(function=_fs_isRunning), overwrite=True, name="fs_isrunning" - ) - fs_isrunning.inputs.logger = LOGGER + if not have_mask: + LOGGER.info("ANAT Stage 6: Preparing mask refinement workflow") + # Stage 6: Refine ANTs mask with FreeSurfer segmentation + refinement_wf = init_refinement_wf() + applyrefined = pe.Node(fsl.ApplyMask(), name="applyrefined") - # 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 - ) - applyrefined = pe.Node(fsl.ApplyMask(), name="applyrefined") - sphere_reg_wf = init_sphere_reg_wf(msm_sulc=msm_sulc, sloppy=sloppy, name="sphere_reg_wf") + # fmt:off + workflow.connect([ + (surface_recon_wf, refinement_wf, [ + ('outputnode.subjects_dir', 'inputnode.subjects_dir'), + ('outputnode.subject_id', 'inputnode.subject_id'), + ('outputnode.fsnative2t1w_xfm', 'inputnode.fsnative2t1w_xfm'), + ]), + (t1w_buffer, refinement_wf, [ + ('t1w_preproc', 'inputnode.reference_image'), + ('ants_seg', 'inputnode.ants_segs'), + ]), + (t1w_buffer, applyrefined, [('t1w_preproc', 'in_file')]), + (refinement_wf, applyrefined, [('outputnode.out_brainmask', 'mask_file')]), + (refinement_wf, refined_buffer, [('outputnode.out_brainmask', 't1w_mask')]), + (applyrefined, refined_buffer, [('out_file', 't1w_brain')]), + ]) + # fmt:on + else: + LOGGER.info("ANAT Found brain mask - skipping Stage 6") - if t2w: + if t2w and not have_t2w: + LOGGER.info("ANAT Stage 7: Creating T2w template") t2w_template_wf = init_anat_template_wf( longitudinal=longitudinal, omp_nthreads=omp_nthreads, @@ -561,6 +1135,14 @@ def _check_img(img): ), name="t2w_resample", ) + + ds_t2w_preproc = pe.Node( + DerivativesDataSink(base_directory=output_dir, desc="preproc", compress=True), + name="ds_t2w_preproc", + run_without_submitting=True, + ) + ds_t2w_preproc.inputs.SkullStripped = False + # fmt:off workflow.connect([ (inputnode, t2w_template_wf, [('t2w', 'inputnode.anat_files')]), @@ -577,93 +1159,193 @@ def _check_img(img): (('outputnode.bias_corrected', _pop), 'reference_image'), ]), (t2wtot1w_xfm, t2w_resample, [('out_xfm', 'transforms')]), - (t2w_resample, outputnode, [('output_image', 't2w_preproc')]), + (inputnode, ds_t2w_preproc, [('t2w', 'source_file')]), + (t2w_resample, ds_t2w_preproc, [('output_image', 'in_file')]), + (ds_t2w_preproc, outputnode, [('out_file', 't2w_preproc')]), ]) # fmt:on + elif not t2w: + LOGGER.info("ANAT No T2w images provided - skipping Stage 7") + else: + LOGGER.info("ANAT Found preprocessed T2w - skipping Stage 7") + + # Stages 8-10: Surface conversion and registration + # sphere_reg is needed to generate sphere_reg_fsLR + # sphere and sulc are needed to generate sphere_reg_msm + # white, pial, midthickness and thickness are needed to resample in the cortical ribbon + # TODO: Consider paring down or splitting into a subworkflow that can be called on-demand + # A subworkflow would still need to check for precomputed outputs + needed_anat_surfs = ["white", "pial", "midthickness"] + needed_metrics = ["thickness", "sulc"] + needed_spheres = ["sphere_reg", "sphere"] + + # Detect pre-computed surfaces + found_surfs = { + surf: sorted(precomputed[surf]) + for surf in needed_anat_surfs + needed_metrics + needed_spheres + if len(precomputed.get(surf, [])) == 2 + } + if found_surfs: + LOGGER.info(f"ANAT Stage 8: Found pre-converted surfaces for {list(found_surfs)}") + surfaces_buffer.inputs.trait_set(**found_surfs) + + # Stage 8: Surface conversion + surfs = [surf for surf in needed_anat_surfs if surf not in found_surfs] + spheres = [sphere for sphere in needed_spheres if sphere not in found_surfs] + if surfs or spheres: + LOGGER.info(f"ANAT Stage 8: Creating GIFTI surfaces for {surfs + spheres}") + if surfs: + gifti_surfaces_wf = init_gifti_surfaces_wf(surfaces=surfs) + ds_surfaces_wf = init_ds_surfaces_wf( + bids_root=bids_root, output_dir=output_dir, surfaces=surfs + ) + # fmt:off + workflow.connect([ + (surface_recon_wf, gifti_surfaces_wf, [ + ('outputnode.subject_id', 'inputnode.subject_id'), + ('outputnode.subjects_dir', 'inputnode.subjects_dir'), + ('outputnode.fsnative2t1w_xfm', 'inputnode.fsnative2t1w_xfm'), + ]), + (gifti_surfaces_wf, surfaces_buffer, [ + (f'outputnode.{surf}', surf) for surf in surfs + ]), + (sourcefile_buffer, ds_surfaces_wf, [("source_files", "inputnode.source_files")]), + (gifti_surfaces_wf, ds_surfaces_wf, [ + (f'outputnode.{surf}', f'inputnode.{surf}') for surf in surfs + ]), + ]) + # fmt:on + if spheres: + gifti_spheres_wf = init_gifti_surfaces_wf( + surfaces=spheres, to_scanner=False, name='gifti_spheres_wf' + ) + ds_spheres_wf = init_ds_surfaces_wf( + bids_root=bids_root, output_dir=output_dir, surfaces=spheres, name='ds_spheres_wf' + ) + # fmt:off + workflow.connect([ + (surface_recon_wf, gifti_spheres_wf, [ + ('outputnode.subject_id', 'inputnode.subject_id'), + ('outputnode.subjects_dir', 'inputnode.subjects_dir'), + # No transform for spheres, following HCP pipelines' lead + ]), + (gifti_spheres_wf, surfaces_buffer, [ + (f'outputnode.{sphere}', sphere) for sphere in spheres + ]), + (sourcefile_buffer, ds_spheres_wf, [("source_files", "inputnode.source_files")]), + (gifti_spheres_wf, ds_spheres_wf, [ + (f'outputnode.{sphere}', f'inputnode.{sphere}') for sphere in spheres + ]), + ]) + # fmt:on + metrics = [metric for metric in needed_metrics if metric not in found_surfs] + if metrics: + LOGGER.info(f"ANAT Stage 8: Creating GIFTI metrics for {metrics}") + gifti_morph_wf = init_gifti_morphometrics_wf(morphometrics=metrics) + ds_morph_wf = init_ds_surface_metrics_wf( + bids_root=bids_root, output_dir=output_dir, metrics=metrics, name='ds_morph_wf' + ) - # 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'), - ('outputnode.out_segm', 'inputnode.ants_segs'), - (('outputnode.bias_corrected', _pop), 'inputnode.corrected_t1')]), - (brain_extraction_wf, applyrefined, [ - (('outputnode.bias_corrected', _pop), 'in_file')]), - (surface_recon_wf, applyrefined, [ - ('outputnode.out_brainmask', 'mask_file')]), - (surface_recon_wf, outputnode, [ - ('outputnode.subjects_dir', 'subjects_dir'), - ('outputnode.subject_id', 'subject_id'), - ('outputnode.t1w2fsnative_xfm', 't1w2fsnative_xfm'), - ('outputnode.fsnative2t1w_xfm', 'fsnative2t1w_xfm'), - ('outputnode.surfaces', 'surfaces'), - ('outputnode.morphometrics', 'morphometrics'), - ('outputnode.out_aseg', 't1w_aseg'), - ('outputnode.out_aparc', 't1w_aparc')]), - (surface_recon_wf, anat_ribbon_wf, [ - ('outputnode.surfaces', 'inputnode.surfaces'), - ('outputnode.out_brainmask', 'inputnode.t1w_mask')]), - (anat_ribbon_wf, outputnode, [ - ("outputnode.anat_ribbon", "anat_ribbon")]), - (applyrefined, buffernode, [('out_file', 't1w_brain')]), - (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'), - ('outputnode.sphere_reg_fsLR', 'sphere_reg_fsLR'), - ]), - (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, anat_derivatives_wf, [ - ('t1w2fsnative_xfm', 'inputnode.t1w2fsnative_xfm'), - ('fsnative2t1w_xfm', 'inputnode.fsnative2t1w_xfm'), - ('surfaces', 'inputnode.surfaces'), - ('morphometrics', 'inputnode.morphometrics'), - ('sphere_reg', 'inputnode.sphere_reg'), - ('sphere_reg_fsLR', 'inputnode.sphere_reg_fsLR'), - ]), - (anat_ribbon_wf, anat_derivatives_wf, [ - ("outputnode.anat_ribbon", "inputnode.anat_ribbon"), - ]), - ]) - # fmt:on - - if cifti_output: - morph_grayords_wf = init_morph_grayords_wf(grayord_density=cifti_output) - anat_derivatives_wf.get_node('inputnode').inputs.cifti_density = cifti_output # fmt:off workflow.connect([ - (surface_recon_wf, morph_grayords_wf, [ + (surface_recon_wf, gifti_morph_wf, [ ('outputnode.subject_id', 'inputnode.subject_id'), ('outputnode.subjects_dir', 'inputnode.subjects_dir'), ]), - (morph_grayords_wf, anat_derivatives_wf, [ - ("outputnode.cifti_morph", "inputnode.cifti_morph"), - ("outputnode.cifti_metadata", "inputnode.cifti_metadata"), + (gifti_morph_wf, surfaces_buffer, [ + (f'outputnode.{metric}', metric) for metric in metrics + ]), + (sourcefile_buffer, ds_morph_wf, [("source_files", "inputnode.source_files")]), + (gifti_morph_wf, ds_morph_wf, [ + (f'outputnode.{metric}', f'inputnode.{metric}') for metric in metrics + ]), + ]) + # fmt:on + + if "anat_ribbon" not in precomputed: + LOGGER.info("ANAT Stage 8a: Creating cortical ribbon mask") + anat_ribbon_wf = init_anat_ribbon_wf() + ds_ribbon_mask_wf = init_ds_mask_wf( + bids_root=bids_root, + output_dir=output_dir, + mask_type="ribbon", + name="ds_ribbon_mask_wf", + ) + # fmt:off + workflow.connect([ + (t1w_buffer, anat_ribbon_wf, [ + ("t1w_preproc", "inputnode.ref_file"), + ]), + (surfaces_buffer, anat_ribbon_wf, [ + ("white", "inputnode.white"), + ("pial", "inputnode.pial"), + ]), + (sourcefile_buffer, ds_ribbon_mask_wf, [("source_files", "inputnode.source_files")]), + (anat_ribbon_wf, ds_ribbon_mask_wf, [ + ("outputnode.anat_ribbon", "inputnode.mask_file"), + ]), + (ds_ribbon_mask_wf, outputnode, [("outputnode.mask_file", "anat_ribbon")]), + ]) + # fmt:on + else: + LOGGER.info("ANAT Stage 8a: Found pre-computed cortical ribbon mask") + outputnode.inputs.anat_ribbon = precomputed["anat_ribbon"] + + # Stage 9: Baseline fsLR registration + if len(precomputed.get("sphere_reg_fsLR", [])) < 2: + LOGGER.info("ANAT Stage 9: Creating fsLR registration sphere") + fsLR_reg_wf = init_fsLR_reg_wf() + ds_fsLR_reg_wf = init_ds_surfaces_wf( + bids_root=bids_root, + output_dir=output_dir, + surfaces=["sphere_reg_fsLR"], + name='ds_fsLR_reg_wf', + ) + + # fmt:off + workflow.connect([ + (surfaces_buffer, fsLR_reg_wf, [('sphere_reg', 'inputnode.sphere_reg')]), + (sourcefile_buffer, ds_fsLR_reg_wf, [("source_files", "inputnode.source_files")]), + (fsLR_reg_wf, ds_fsLR_reg_wf, [ + ('outputnode.sphere_reg_fsLR', 'inputnode.sphere_reg_fsLR') ]), + (ds_fsLR_reg_wf, fsLR_buffer, [('outputnode.sphere_reg_fsLR', 'sphere_reg_fsLR')]), ]) # fmt:on + else: + LOGGER.info("ANAT Stage 9: Found pre-computed fsLR registration sphere") + fsLR_buffer.inputs.sphere_reg_fsLR = sorted(precomputed["sphere_reg_fsLR"]) + + # Stage 10: MSMSulc + if msm_sulc and len(precomputed.get("sphere_reg_msm", [])) < 2: + LOGGER.info("ANAT Stage 10: Creating MSM-Sulc registration sphere") + msm_sulc_wf = init_msm_sulc_wf(sloppy=sloppy) + ds_msmsulc_wf = init_ds_surfaces_wf( + bids_root=bids_root, + output_dir=output_dir, + surfaces=["sphere_reg_msm"], + name='ds_msmsulc_wf', + ) + + # fmt:off + workflow.connect([ + (surfaces_buffer, msm_sulc_wf, [ + ('sulc', 'inputnode.sulc'), + ('sphere', 'inputnode.sphere'), + ]), + (fsLR_buffer, msm_sulc_wf, [('sphere_reg_fsLR', 'inputnode.sphere_reg_fsLR')]), + (sourcefile_buffer, ds_msmsulc_wf, [("source_files", "inputnode.source_files")]), + (msm_sulc_wf, ds_msmsulc_wf, [ + ('outputnode.sphere_reg_fsLR', 'inputnode.sphere_reg_msm') + ]), + (ds_msmsulc_wf, msm_buffer, [('outputnode.sphere_reg_msm', 'sphere_reg_msm')]), + ]) + # fmt:on + elif msm_sulc: + LOGGER.info("ANAT Stage 10: Found pre-computed MSM-Sulc registration sphere") + msm_buffer.inputs.sphere_reg_msm = sorted(precomputed["sphere_reg_msm"]) + else: + LOGGER.info("ANAT Stage 10: MSM-Sulc disabled") return workflow @@ -723,7 +1405,7 @@ def init_anat_template_wf( workflow = Workflow(name=name) if num_files > 1: - fs_ver = fs.Info().looseversion() or "" + fs_ver = fs.Info().looseversion() or "(version unknown)" workflow.__desc__ = f"""\ An anatomical {contrast}-reference map was computed after registration of {num_files} {contrast} images (after INU-correction) using @@ -921,3 +1603,20 @@ def _split_segments(in_file): def _probseg_fast2bids(inlist): """Reorder a list of probseg maps from FAST (CSF, WM, GM) to BIDS (GM, WM, CSF).""" return (inlist[1], inlist[2], inlist[0]) + + +def _is_skull_stripped(img): + """Check if T1w images are skull-stripped.""" + import nibabel as nb + import numpy as np + + data = nb.load(img).dataobj + sidevals = ( + np.abs(data[0, :, :]).sum() + + np.abs(data[-1, :, :]).sum() + + np.abs(data[:, 0, :]).sum() + + np.abs(data[:, -1, :]).sum() + + np.abs(data[:, :, 0]).sum() + + np.abs(data[:, :, -1]).sum() + ) + return sidevals < 10 diff --git a/smriprep/workflows/base.py b/smriprep/workflows/base.py index 9ef3884b2f..0c2dc70262 100644 --- a/smriprep/workflows/base.py +++ b/smriprep/workflows/base.py @@ -23,7 +23,6 @@ """*sMRIPrep* base processing workflows.""" import sys import os -from pathlib import Path from copy import deepcopy from nipype import __version__ as nipype_ver @@ -45,7 +44,7 @@ def init_smriprep_wf( *, sloppy, debug, - fast_track, + derivatives, freesurfer, fs_subjects_dir, hires, @@ -63,6 +62,7 @@ def init_smriprep_wf( subject_list, work_dir, bids_filters, + cifti_output, ): """ Create the execution graph of *sMRIPrep*, with a sub-workflow for each subject. @@ -81,10 +81,12 @@ def init_smriprep_wf( os.environ['FREESURFER_HOME'] = os.getcwd() from smriprep.workflows.base import init_smriprep_wf from niworkflows.utils.spaces import SpatialReferences, Reference + spaces = SpatialReferences(spaces=['MNI152NLin2009cAsym', 'fsaverage5']) + spaces.checkpoint() wf = init_smriprep_wf( sloppy=False, debug=False, - fast_track=False, + derivatives=[], freesurfer=True, fs_subjects_dir=None, hires=True, @@ -98,10 +100,11 @@ def init_smriprep_wf( skull_strip_fixed_seed=False, skull_strip_mode='force', skull_strip_template=Reference('OASIS30ANTs'), - spaces=SpatialReferences(spaces=['MNI152NLin2009cAsym', 'fsaverage5']), + spaces=spaces, subject_list=['smripreptest'], work_dir='.', bids_filters=None, + cifti_output=None, ) Parameters @@ -110,7 +113,7 @@ def init_smriprep_wf( Quick, impercise operations. Used to decrease workflow duration. debug : :obj:`bool` Enable debugging outputs - fast_track : :obj:`bool` + derivatives : :obj:`list` of directories Fast-track the workflow by searching for existing derivatives. freesurfer : :obj:`bool` Enable FreeSurfer surface reconstruction (may increase runtime) @@ -175,7 +178,7 @@ def init_smriprep_wf( sloppy=sloppy, debug=debug, freesurfer=freesurfer, - fast_track=fast_track, + derivatives=derivatives, hires=hires, layout=layout, longitudinal=longitudinal, @@ -190,6 +193,7 @@ def init_smriprep_wf( spaces=spaces, subject_id=subject_id, bids_filters=bids_filters, + cifti_output=cifti_output, ) single_subject_wf.config["execution"]["crashdump_dir"] = os.path.join( @@ -209,7 +213,7 @@ def init_single_subject_wf( *, sloppy, debug, - fast_track, + derivatives, freesurfer, hires, layout, @@ -225,6 +229,7 @@ def init_single_subject_wf( spaces, subject_id, bids_filters, + cifti_output, ): """ Create a single subject workflow. @@ -247,11 +252,13 @@ def init_single_subject_wf( from niworkflows.utils.spaces import SpatialReferences, Reference from smriprep.workflows.base import init_single_subject_wf BIDSLayout = namedtuple('BIDSLayout', ['root']) + spaces = SpatialReferences(spaces=['MNI152NLin2009cAsym', 'fsaverage5']) + spaces.checkpoint() wf = init_single_subject_wf( sloppy=False, debug=False, freesurfer=True, - fast_track=False, + derivatives=[], hires=True, layout=BIDSLayout('.'), longitudinal=False, @@ -263,9 +270,10 @@ def init_single_subject_wf( skull_strip_fixed_seed=False, skull_strip_mode='force', skull_strip_template=Reference('OASIS30ANTs'), - spaces=SpatialReferences(spaces=['MNI152NLin2009cAsym', 'fsaverage5']), + spaces=spaces, subject_id='test', bids_filters=None, + cifti_output=None, ) Parameters @@ -274,8 +282,8 @@ def init_single_subject_wf( Quick, impercise operations. Used to decrease workflow duration. debug : :obj:`bool` Enable debugging outputs - fast_track : :obj:`bool` - If ``True``, attempt to collect previously run derivatives. + derivatives : :obj:`list` of directories + Fast-track the workflow by searching for existing derivatives. freesurfer : :obj:`bool` Enable FreeSurfer surface reconstruction (may increase runtime) hires : :obj:`bool` @@ -356,14 +364,13 @@ def init_single_subject_wf( """ - deriv_cache = None - if fast_track: - from ..utils.bids import collect_derivatives + from ..utils.bids import collect_derivatives - std_spaces = spaces.get_spaces(nonstandard=False, dim=(3,)) - deriv_cache = collect_derivatives( - Path(output_dir) / "smriprep", subject_id, std_spaces, freesurfer - ) + deriv_cache = {} + std_spaces = spaces.get_spaces(nonstandard=False, dim=(3,)) + std_spaces.append("fsnative") + for deriv_dir in derivatives: + deriv_cache.update(collect_derivatives(deriv_dir, subject_id, std_spaces)) inputnode = pe.Node(niu.IdentityInterface(fields=["subjects_dir"]), name="inputnode") @@ -412,7 +419,7 @@ def init_single_subject_wf( bids_root=layout.root, sloppy=sloppy, debug=debug, - existing_derivatives=deriv_cache, + precomputed=deriv_cache, freesurfer=freesurfer, hires=hires, longitudinal=longitudinal, @@ -426,7 +433,7 @@ def init_single_subject_wf( skull_strip_mode=skull_strip_mode, skull_strip_template=skull_strip_template, spaces=spaces, - cifti_output=False, # Enabling this needs a CLI flag + cifti_output=cifti_output, ) # fmt:off diff --git a/smriprep/workflows/fit/__init__.py b/smriprep/workflows/fit/__init__.py new file mode 100644 index 0000000000..e69de29bb2 diff --git a/smriprep/workflows/norm.py b/smriprep/workflows/fit/registration.py similarity index 63% rename from smriprep/workflows/norm.py rename to smriprep/workflows/fit/registration.py index 46f89d5bf4..a8ed1926ee 100644 --- a/smriprep/workflows/norm.py +++ b/smriprep/workflows/fit/registration.py @@ -32,16 +32,15 @@ from templateflow.api import get_metadata from niworkflows.engine.workflows import LiterateWorkflow as Workflow from niworkflows.interfaces.norm import SpatialNormalization -from niworkflows.interfaces.fixes import FixHeaderApplyTransforms as ApplyTransforms -from ..interfaces.templateflow import TemplateFlowSelect, TemplateDesc +from ...interfaces.templateflow import TemplateFlowSelect, TemplateDesc -def init_anat_norm_wf( +def init_register_template_wf( *, sloppy, omp_nthreads, templates, - name="anat_norm_wf", + name="register_template_wf", ): """ Build an individual spatial normalization workflow using ``antsRegistration``. @@ -51,8 +50,8 @@ def init_anat_norm_wf( :graph2use: orig :simple_form: yes - from smriprep.workflows.norm import init_anat_norm_wf - wf = init_anat_norm_wf( + from smriprep.workflows.fit.registration import init_register_template_wf + wf = init_register_template_wf( sloppy=False, omp_nthreads=1, templates=['MNI152NLin2009cAsym', 'MNI152NLin6Asym'], @@ -84,37 +83,18 @@ def init_anat_norm_wf( ------ moving_image The input image that will be normalized to standard space. - moving_mask - A precise brain mask separating skull/skin/fat from brain - structures. - moving_segmentation - A brain tissue segmentation of the ``moving_image``. - moving_tpms - tissue probability maps (TPMs) corresponding to the - ``moving_segmentation``. lesion_mask (optional) A mask to exclude regions from the cost-function input domain to enable standardization of lesioned brains. - orig_t1w - The original T1w image from the BIDS structure. template Template name and specification Outputs ------- - standardized - The T1w after spatial normalization, in template space. anat2std_xfm The T1w-to-template transform. std2anat_xfm The template-to-T1w transform. - std_mask - The ``moving_mask`` in template space (matches ``standardized`` output). - std_dseg - The ``moving_segmentation`` in template space (matches ``standardized`` - output). - std_tpms - The ``moving_tpms`` in template space (matches ``standardized`` output). template Template name extracted from the input parameter ``template``, for further use in downstream nodes. @@ -167,9 +147,6 @@ def init_anat_norm_wf( "lesion_mask", "moving_image", "moving_mask", - "moving_segmentation", - "moving_tpms", - "orig_t1w", "template", ] ), @@ -179,15 +156,11 @@ def init_anat_norm_wf( out_fields = [ "anat2std_xfm", - "standardized", "std2anat_xfm", - "std_dseg", - "std_mask", - "std_tpms", "template", "template_spec", ] - poutputnode = pe.Node(niu.IdentityInterface(fields=out_fields), name="poutputnode") + outputnode = _make_outputnode(workflow, out_fields, joinsource='inputnode') split_desc = pe.Node(TemplateDesc(), run_without_submitting=True, name="split_desc") @@ -213,73 +186,43 @@ def init_anat_norm_wf( mem_gb=2, ) - # Resample T1w-space inputs - tpl_moving = pe.Node( - ApplyTransforms( - dimension=3, - default_value=0, - float=True, - interpolation="LanczosWindowedSinc", - ), - name="tpl_moving", - ) - - std_mask = pe.Node(ApplyTransforms(interpolation="MultiLabel"), name="std_mask") - std_dseg = pe.Node(ApplyTransforms(interpolation="MultiLabel"), name="std_dseg") - - std_tpms = pe.MapNode( - ApplyTransforms(dimension=3, default_value=0, float=True, interpolation="Gaussian"), - iterfield=["input_image"], - name="std_tpms", - ) - # fmt:off workflow.connect([ (inputnode, split_desc, [('template', 'template')]), - (inputnode, poutputnode, [('template', 'template')]), (inputnode, trunc_mov, [('moving_image', 'op1')]), (inputnode, registration, [ ('moving_mask', 'moving_mask'), ('lesion_mask', 'lesion_mask')]), - (inputnode, tpl_moving, [('moving_image', 'input_image')]), - (inputnode, std_mask, [('moving_mask', 'input_image')]), - (split_desc, tf_select, [('name', 'template'), - ('spec', 'template_spec')]), - (split_desc, registration, [('name', 'template'), - ('spec', 'template_spec')]), - (tf_select, tpl_moving, [('t1w_file', 'reference_image')]), - (tf_select, std_mask, [('t1w_file', 'reference_image')]), - (tf_select, std_dseg, [('t1w_file', 'reference_image')]), - (tf_select, std_tpms, [('t1w_file', 'reference_image')]), + (split_desc, tf_select, [ + ('name', 'template'), + ('spec', 'template_spec'), + ]), + (split_desc, registration, [ + ('name', 'template'), + ('spec', 'template_spec'), + ]), (trunc_mov, registration, [ ('output_image', 'moving_image')]), - (registration, tpl_moving, [('composite_transform', 'transforms')]), - (registration, std_mask, [('composite_transform', 'transforms')]), - (inputnode, std_dseg, [('moving_segmentation', 'input_image')]), - (registration, std_dseg, [('composite_transform', 'transforms')]), - (inputnode, std_tpms, [('moving_tpms', 'input_image')]), - (registration, std_tpms, [('composite_transform', 'transforms')]), - (registration, poutputnode, [ + (split_desc, outputnode, [ + ('name', 'template'), + ('spec', 'template_spec'), + ]), + (registration, outputnode, [ ('composite_transform', 'anat2std_xfm'), - ('inverse_composite_transform', 'std2anat_xfm')]), - (tpl_moving, poutputnode, [('output_image', 'standardized')]), - (std_mask, poutputnode, [('output_image', 'std_mask')]), - (std_dseg, poutputnode, [('output_image', 'std_dseg')]), - (std_tpms, poutputnode, [('output_image', 'std_tpms')]), - (split_desc, poutputnode, [('spec', 'template_spec')]), - ]) - # fmt:on - - # Provide synchronized output - outputnode = pe.JoinNode( - niu.IdentityInterface(fields=out_fields), - name="outputnode", - joinsource="inputnode", - ) - # fmt:off - workflow.connect([ - (poutputnode, outputnode, [(f, f) for f in out_fields]), + ('inverse_composite_transform', 'std2anat_xfm'), + ]), ]) # fmt:on return workflow + + +def _make_outputnode(workflow, out_fields, joinsource): + if joinsource: + pout = pe.Node(niu.IdentityInterface(fields=out_fields), name="poutputnode") + out = pe.JoinNode( + niu.IdentityInterface(fields=out_fields), name="outputnode", joinsource=joinsource + ) + workflow.connect([(pout, out, [(f, f) for f in out_fields])]) + return pout + return pe.Node(niu.IdentityInterface(fields=out_fields), name="outputnode") diff --git a/smriprep/workflows/outputs.py b/smriprep/workflows/outputs.py index 975f2535fe..e6df726074 100644 --- a/smriprep/workflows/outputs.py +++ b/smriprep/workflows/outputs.py @@ -21,17 +21,23 @@ # https://www.nipreps.org/community/licensing/ # """Writing outputs.""" +import typing as ty + from nipype.pipeline import engine as pe from nipype.interfaces import utility as niu -from niworkflows.interfaces.nibabel import ApplyMask +from niworkflows.interfaces.fixes import FixHeaderApplyTransforms as ApplyTransforms +from niworkflows.interfaces.nibabel import ApplyMask, GenerateSamplingReference +from niworkflows.interfaces.space import SpaceDataSource +from niworkflows.interfaces.utility import KeySelect from niworkflows.engine.workflows import LiterateWorkflow as Workflow from ..interfaces import DerivativesDataSink +from ..interfaces.templateflow import TemplateFlowSelect BIDS_TISSUE_ORDER = ("GM", "WM", "CSF") -def init_anat_reports_wf(*, freesurfer, output_dir, name="anat_reports_wf"): +def init_anat_reports_wf(*, spaces, freesurfer, output_dir, name="anat_reports_wf"): """ Set up a battery of datasinks to store reports in the right location. @@ -73,19 +79,18 @@ def init_anat_reports_wf(*, freesurfer, output_dir, name="anat_reports_wf"): SimpleBeforeAfterRPT as SimpleBeforeAfter, ) from niworkflows.interfaces.reportlets.masks import ROIsPlot - from ..interfaces.templateflow import TemplateFlowSelect workflow = Workflow(name=name) inputfields = [ "source_file", - "t1w_conform_report", "t1w_preproc", - "t1w_dseg", "t1w_mask", + "t1w_dseg", "template", - "std_t1w", - "std_mask", + "anat2std_xfm", + # May be missing + "t1w_conform_report", "subject_id", "subjects_dir", ] @@ -124,44 +129,77 @@ def init_anat_reports_wf(*, freesurfer, output_dir, name="anat_reports_wf"): ]) # fmt:on - # Generate reportlets showing spatial normalization - tf_select = pe.Node( - TemplateFlowSelect(resolution=1), name="tf_select", run_without_submitting=True - ) - norm_msk = pe.Node( - niu.Function( - function=_rpt_masks, - output_names=["before", "after"], - input_names=["mask_file", "before", "after", "after_mask"], - ), - name="norm_msk", - ) - norm_rpt = pe.Node(SimpleBeforeAfter(), name="norm_rpt", mem_gb=0.1) - norm_rpt.inputs.after_label = "Participant" # after + if getattr(spaces, "_cached") is not None and spaces.cached.references: + template_iterator_wf = init_template_iterator_wf(spaces=spaces) + t1w_std = pe.Node( + ApplyTransforms( + dimension=3, + default_value=0, + float=True, + interpolation="LanczosWindowedSinc", + ), + name="t1w_std", + ) + mask_std = pe.Node( + ApplyTransforms( + dimension=3, + default_value=0, + float=True, + interpolation="MultiLabel", + ), + name="mask_std", + ) - ds_std_t1w_report = pe.Node( - DerivativesDataSink(base_directory=output_dir, suffix="T1w", datatype="figures"), - name="ds_std_t1w_report", - run_without_submitting=True, - ) + # Generate reportlets showing spatial normalization + norm_msk = pe.Node( + niu.Function( + function=_rpt_masks, + output_names=["before", "after"], + input_names=["mask_file", "before", "after", "after_mask"], + ), + name="norm_msk", + ) + norm_rpt = pe.Node(SimpleBeforeAfter(), name="norm_rpt", mem_gb=0.1) + norm_rpt.inputs.after_label = "Participant" # after - # fmt:off - workflow.connect([ - (inputnode, tf_select, [(('template', _drop_cohort), 'template'), - (('template', _pick_cohort), 'cohort')]), - (inputnode, norm_rpt, [('template', 'before_label')]), - (inputnode, norm_msk, [('std_t1w', 'after'), - ('std_mask', 'after_mask')]), - (tf_select, norm_msk, [('t1w_file', 'before'), - ('brain_mask', 'mask_file')]), - (norm_msk, norm_rpt, [('before', 'before'), - ('after', 'after')]), - (inputnode, ds_std_t1w_report, [ - (('template', _fmt), 'space'), - ('source_file', 'source_file')]), - (norm_rpt, ds_std_t1w_report, [('out_report', 'in_file')]), - ]) - # fmt:on + ds_std_t1w_report = pe.Node( + DerivativesDataSink(base_directory=output_dir, suffix="T1w", datatype="figures"), + name="ds_std_t1w_report", + run_without_submitting=True, + ) + + # fmt:off + workflow.connect([ + (inputnode, template_iterator_wf, [ + ("template", "inputnode.template"), + ("anat2std_xfm", "inputnode.anat2std_xfm"), + ]), + (inputnode, t1w_std, [("t1w_preproc", "input_image")]), + (inputnode, mask_std, [("t1w_mask", "input_image")]), + (template_iterator_wf, t1w_std, [ + ("outputnode.anat2std_xfm", "transforms"), + ("outputnode.std_t1w", "reference_image"), + ]), + (template_iterator_wf, mask_std, [ + ("outputnode.anat2std_xfm", "transforms"), + ("outputnode.std_t1w", "reference_image"), + ]), + (template_iterator_wf, norm_rpt, [('outputnode.space', 'before_label')]), + (t1w_std, norm_msk, [('output_image', 'after')]), + (mask_std, norm_msk, [('output_image', 'after_mask')]), + (template_iterator_wf, norm_msk, [ + ('outputnode.std_t1w', 'before'), + ('outputnode.std_mask', 'mask_file'), + ]), + (norm_msk, norm_rpt, [ + ('before', 'before'), + ('after', 'after'), + ]), + (inputnode, ds_std_t1w_report, [('source_file', 'source_file')]), + (template_iterator_wf, ds_std_t1w_report, [('outputnode.space', 'space')]), + (norm_rpt, ds_std_t1w_report, [('out_report', 'in_file')]), + ]) + # fmt:on if freesurfer: from ..interfaces.reports import FSSurfaceReport @@ -186,42 +224,26 @@ def init_anat_reports_wf(*, freesurfer, output_dir, name="anat_reports_wf"): return workflow -def init_anat_derivatives_wf( +def init_ds_template_wf( *, - bids_root, - freesurfer, num_t1w, - t2w, output_dir, - spaces, - cifti_output, - name="anat_derivatives_wf", - tpm_labels=BIDS_TISSUE_ORDER, + name="ds_template_wf", ): """ - Set up a battery of datasinks to store derivatives in the right location. + Save the subject-specific template Parameters ---------- - bids_root : :obj:`str` - Root path of BIDS dataset - freesurfer : :obj:`bool` - FreeSurfer was enabled num_t1w : :obj:`int` Number of T1w images output_dir : :obj:`str` Directory in which to save derivatives name : :obj:`str` - Workflow name (default: anat_derivatives_wf) - tpm_labels : :obj:`tuple` - Tissue probability maps in order - cifti_output : :obj:`bool` - Whether the ``--cifti-output`` flag was set. + Workflow name (default: ds_template_wf) Inputs ------ - template - Template space and specifications source_files List of input T1w images t1w_ref_xfms @@ -229,87 +251,26 @@ def init_anat_derivatives_wf( t1w_preproc The T1w reference map, which is calculated as the average of bias-corrected and preprocessed T1w images, defining the anatomical space. - t1w_mask - Mask of the ``t1w_preproc`` - t1w_dseg - Segmentation in T1w space - t1w_tpms - Tissue probability maps in T1w space - t2w_preproc - The preprocessed T2w image, bias-corrected and resampled into anatomical space. - anat2std_xfm - Nonlinear spatial transform to resample imaging data given in anatomical space - into standard space. - std2anat_xfm - Inverse transform of ``anat2std_xfm`` - std_t1w - T1w reference resampled in one or more standard spaces. - std_mask - Mask of skull-stripped template, in standard space - std_dseg - Segmentation, resampled into standard space - std_tpms - Tissue probability maps in standard space - t1w2fsnative_xfm - LTA-style affine matrix translating from T1w to - FreeSurfer-conformed subject space - fsnative2t1w_xfm - LTA-style affine matrix translating from FreeSurfer-conformed - subject space to T1w - surfaces - GIFTI surfaces (gray/white boundary, midthickness, pial, inflated) - morphometrics - GIFTIs of cortical thickness, curvature, and sulcal depth - anat_ribbon - Cortical ribbon volume in T1w space - t1w_fs_aseg - FreeSurfer's aseg segmentation, in native T1w space - t1w_fs_aparc - FreeSurfer's aparc+aseg segmentation, in native T1w space - cifti_morph - Morphometric CIFTI-2 dscalar files - cifti_density - Grayordinate density - cifti_metadata - JSON files containing metadata dictionaries - """ - from niworkflows.interfaces.utility import KeySelect + Outputs + ------- + t1w_preproc + The location in the output directory of the preprocessed T1w image + """ workflow = Workflow(name=name) inputnode = pe.Node( niu.IdentityInterface( fields=[ - "template", "source_files", "t1w_ref_xfms", "t1w_preproc", - "t1w_mask", - "t1w_dseg", - "t1w_tpms", - "t2w_preproc", - "anat2std_xfm", - "std2anat_xfm", - "t1w2fsnative_xfm", - "fsnative2t1w_xfm", - "surfaces", - "sphere_reg", - "sphere_reg_fsLR", - "morphometrics", - "anat_ribbon", - "t1w_fs_aseg", - "t1w_fs_aparc", - 'cifti_metadata', - 'cifti_density', - 'cifti_morph', ] ), name="inputnode", ) - - raw_sources = pe.Node(niu.Function(function=_bids_relative), name="raw_sources") - raw_sources.inputs.bids_root = bids_root + outputnode = pe.Node(niu.IdentityInterface(fields=["t1w_preproc"]), name="outputnode") ds_t1w_preproc = pe.Node( DerivativesDataSink(base_directory=output_dir, desc="preproc", compress=True), @@ -318,72 +279,14 @@ def init_anat_derivatives_wf( ) ds_t1w_preproc.inputs.SkullStripped = False - ds_t1w_mask = pe.Node( - DerivativesDataSink(base_directory=output_dir, desc="brain", suffix="mask", compress=True), - name="ds_t1w_mask", - run_without_submitting=True, - ) - ds_t1w_mask.inputs.Type = "Brain" - - ds_t1w_dseg = pe.Node( - DerivativesDataSink(base_directory=output_dir, suffix="dseg", compress=True), - name="ds_t1w_dseg", - run_without_submitting=True, - ) - - ds_t1w_tpms = pe.Node( - DerivativesDataSink(base_directory=output_dir, suffix="probseg", compress=True), - name="ds_t1w_tpms", - run_without_submitting=True, - ) - ds_t1w_tpms.inputs.label = tpm_labels - # fmt:off workflow.connect([ - (inputnode, raw_sources, [('source_files', 'in_files')]), (inputnode, ds_t1w_preproc, [('t1w_preproc', 'in_file'), ('source_files', 'source_file')]), - (inputnode, ds_t1w_mask, [('t1w_mask', 'in_file'), - ('source_files', 'source_file')]), - (inputnode, ds_t1w_tpms, [('t1w_tpms', 'in_file'), - ('source_files', 'source_file')]), - (inputnode, ds_t1w_dseg, [('t1w_dseg', 'in_file'), - ('source_files', 'source_file')]), - (raw_sources, ds_t1w_mask, [('out', 'RawSources')]), + (ds_t1w_preproc, outputnode, [("out_file", "t1w_preproc")]), ]) # fmt:on - # Transforms - if spaces.get_spaces(nonstandard=False, dim=(3,)): - ds_std2t1w_xfm = pe.MapNode( - DerivativesDataSink(base_directory=output_dir, to="T1w", mode="image", suffix="xfm"), - iterfield=("in_file", "from"), - name="ds_std2t1w_xfm", - run_without_submitting=True, - ) - - ds_t1w2std_xfm = pe.MapNode( - DerivativesDataSink( - base_directory=output_dir, mode="image", suffix="xfm", **{"from": "T1w"} - ), - iterfield=("in_file", "to"), - name="ds_t1w2std_xfm", - run_without_submitting=True, - ) - - # fmt:off - workflow.connect([ - (inputnode, ds_t1w2std_xfm, [ - ('anat2std_xfm', 'in_file'), - (('template', _combine_cohort), 'to'), - ('source_files', 'source_file')]), - (inputnode, ds_std2t1w_xfm, [ - ('std2anat_xfm', 'in_file'), - (('template', _combine_cohort), 'from'), - ('source_files', 'source_file')]), - ]) - # fmt:on - if num_t1w > 1: # Please note the dictionary unpacking to provide the from argument. # It is necessary because from is a protected keyword (not allowed as argument name). @@ -407,266 +310,783 @@ def init_anat_derivatives_wf( ]) # fmt:on - # Write derivatives in standard spaces specified by --output-spaces - if getattr(spaces, "_cached") is not None and spaces.cached.references: - from niworkflows.interfaces.space import SpaceDataSource - from niworkflows.interfaces.nibabel import GenerateSamplingReference - from niworkflows.interfaces.fixes import ( - FixHeaderApplyTransforms as ApplyTransforms, - ) - - from ..interfaces.templateflow import TemplateFlowSelect + return workflow - spacesource = pe.Node(SpaceDataSource(), name="spacesource", run_without_submitting=True) - spacesource.iterables = ( - "in_tuple", - [(s.fullname, s.spec) for s in spaces.cached.get_standard(dim=(3,))], - ) - gen_tplid = pe.Node( - niu.Function(function=_fmt_cohort), - name="gen_tplid", - run_without_submitting=True, - ) +def init_ds_mask_wf( + *, + bids_root: str, + output_dir: str, + mask_type: str, + name="ds_mask_wf", +): + """ + Save the subject brain mask - select_xfm = pe.Node( - KeySelect(fields=["anat2std_xfm"]), - name="select_xfm", - run_without_submitting=True, - ) - select_tpl = pe.Node(TemplateFlowSelect(), name="select_tpl", run_without_submitting=True) + Parameters + ---------- + bids_root : :obj:`str` + Root path of BIDS dataset + output_dir : :obj:`str` + Directory in which to save derivatives + name : :obj:`str` + Workflow name (default: ds_mask_wf) - gen_ref = pe.Node(GenerateSamplingReference(), name="gen_ref", mem_gb=0.01) + Inputs + ------ + source_files + List of input T1w images + mask_file + Mask to save - # Mask T1w preproc images - mask_t1w = pe.Node(ApplyMask(), name='mask_t1w') + Outputs + ------- + mask_file + The location in the output directory of the mask - # Resample T1w-space inputs - anat2std_t1w = pe.Node( - ApplyTransforms( - dimension=3, - default_value=0, - float=True, - interpolation="LanczosWindowedSinc", - ), - name="anat2std_t1w", - ) + """ + workflow = Workflow(name=name) - anat2std_mask = pe.Node(ApplyTransforms(interpolation="MultiLabel"), name="anat2std_mask") - anat2std_dseg = pe.Node(ApplyTransforms(interpolation="MultiLabel"), name="anat2std_dseg") - anat2std_tpms = pe.MapNode( - ApplyTransforms(dimension=3, default_value=0, float=True, interpolation="Gaussian"), - iterfield=["input_image"], - name="anat2std_tpms", - ) + inputnode = pe.Node( + niu.IdentityInterface(fields=["source_files", "mask_file"]), + name="inputnode", + ) + outputnode = pe.Node(niu.IdentityInterface(fields=["mask_file"]), name="outputnode") - ds_std_t1w = pe.Node( - DerivativesDataSink( - base_directory=output_dir, - desc="preproc", - compress=True, - ), - name="ds_std_t1w", - run_without_submitting=True, - ) - ds_std_t1w.inputs.SkullStripped = True + raw_sources = pe.Node(niu.Function(function=_bids_relative), name="raw_sources") + raw_sources.inputs.bids_root = bids_root - ds_std_mask = pe.Node( - DerivativesDataSink( - base_directory=output_dir, desc="brain", suffix="mask", compress=True - ), - name="ds_std_mask", - run_without_submitting=True, - ) - ds_std_mask.inputs.Type = "Brain" + ds_mask = pe.Node( + DerivativesDataSink( + base_directory=output_dir, + desc=mask_type, + suffix="mask", + compress=True, + ), + name="ds_t1w_mask", + run_without_submitting=True, + ) + if mask_type == "brain": + ds_mask.inputs.Type = "Brain" + else: + ds_mask.inputs.Type = "ROI" - ds_std_dseg = pe.Node( - DerivativesDataSink(base_directory=output_dir, suffix="dseg", compress=True), - name="ds_std_dseg", - run_without_submitting=True, - ) + # fmt:off + workflow.connect([ + (inputnode, raw_sources, [('source_files', 'in_files')]), + (inputnode, ds_mask, [('mask_file', 'in_file'), + ('source_files', 'source_file')]), + (raw_sources, ds_mask, [('out', 'RawSources')]), + (ds_mask, outputnode, [('out_file', 'mask_file')]), + ]) + # fmt:on - ds_std_tpms = pe.Node( - DerivativesDataSink(base_directory=output_dir, suffix="probseg", compress=True), - name="ds_std_tpms", - run_without_submitting=True, - ) + return workflow - # CRITICAL: the sequence of labels here (CSF-GM-WM) is that of the output of FSL-FAST - # (intensity mean, per tissue). This order HAS to be matched also by the ``tpms`` - # output in the data/io_spec.json file. - ds_std_tpms.inputs.label = tpm_labels - # fmt:off - workflow.connect([ - (inputnode, mask_t1w, [('t1w_preproc', 'in_file'), - ('t1w_mask', 'in_mask')]), - (mask_t1w, anat2std_t1w, [('out_file', 'input_image')]), - (inputnode, anat2std_mask, [('t1w_mask', 'input_image')]), - (inputnode, anat2std_dseg, [('t1w_dseg', 'input_image')]), - (inputnode, anat2std_tpms, [('t1w_tpms', 'input_image')]), - (inputnode, gen_ref, [('t1w_preproc', 'moving_image')]), - (inputnode, select_xfm, [ - ('anat2std_xfm', 'anat2std_xfm'), - ('template', 'keys')]), - (spacesource, gen_tplid, [('space', 'template'), - ('cohort', 'cohort')]), - (gen_tplid, select_xfm, [('out', 'key')]), - (spacesource, select_tpl, [('space', 'template'), - ('cohort', 'cohort'), - (('resolution', _no_native), 'resolution')]), - (spacesource, gen_ref, [(('resolution', _is_native), 'keep_native')]), - (select_tpl, gen_ref, [('t1w_file', 'fixed_image')]), - (anat2std_t1w, ds_std_t1w, [('output_image', 'in_file')]), - (anat2std_mask, ds_std_mask, [('output_image', 'in_file')]), - (anat2std_dseg, ds_std_dseg, [('output_image', 'in_file')]), - (anat2std_tpms, ds_std_tpms, [('output_image', 'in_file')]), - (select_tpl, ds_std_mask, [(('brain_mask', _drop_path), 'RawSources')]), - ]) - workflow.connect( - # Connect apply transforms nodes - [ - (gen_ref, n, [('out_file', 'reference_image')]) - for n in (anat2std_t1w, anat2std_mask, anat2std_dseg, anat2std_tpms) - ] - + [ - (select_xfm, n, [('anat2std_xfm', 'transforms')]) - for n in (anat2std_t1w, anat2std_mask, anat2std_dseg, anat2std_tpms) - ] - # Connect the source_file input of these datasinks - + [ - (inputnode, n, [('source_files', 'source_file')]) - for n in (ds_std_t1w, ds_std_mask, ds_std_dseg, ds_std_tpms) - ] - # Connect the space input of these datasinks - + [ - (spacesource, n, [ - ('space', 'space'), ('cohort', 'cohort'), ('resolution', 'resolution') - ]) - for n in (ds_std_t1w, ds_std_mask, ds_std_dseg, ds_std_tpms) - ] - ) - # fmt:on +def init_ds_dseg_wf(*, output_dir, name="ds_dseg_wf"): + """ + Save discrete segmentations - if not freesurfer: - return workflow + Parameters + ---------- + output_dir : :obj:`str` + Directory in which to save derivatives + name : :obj:`str` + Workflow name (default: ds_dseg_wf) - # T2w coregistration requires FreeSurfer surfaces, so only try to save if freesurfer - if t2w: - ds_t2w_preproc = pe.Node( - DerivativesDataSink( - base_directory=output_dir, - desc="preproc", - suffix="T2w", - compress=True, - ), - name="ds_t2w_preproc", - run_without_submitting=True, - ) - ds_t2w_preproc.inputs.SkullStripped = False - ds_t2w_preproc.inputs.source_file = t2w + Inputs + ------ + source_files + List of input T1w images + t1w_dseg + Segmentation in T1w space - # fmt:off - workflow.connect([ - (inputnode, ds_t2w_preproc, [('t2w_preproc', 'in_file')]), - ]) - # fmt:on + Outputs + ------- + t1w_dseg + The location in the output directory of the discrete segmentation - from niworkflows.interfaces.nitransforms import ConcatenateXFMs - from niworkflows.interfaces.surf import Path2BIDS + """ + workflow = Workflow(name=name) - # FS native space transforms - lta2itk_fwd = pe.Node(ConcatenateXFMs(), name="lta2itk_fwd", run_without_submitting=True) - lta2itk_inv = pe.Node(ConcatenateXFMs(), name="lta2itk_inv", run_without_submitting=True) - ds_t1w_fsnative = pe.Node( - DerivativesDataSink( - base_directory=output_dir, - mode="image", - to="fsnative", - suffix="xfm", - extension="txt", - **{"from": "T1w"}, - ), - name="ds_t1w_fsnative", - run_without_submitting=True, + inputnode = pe.Node( + niu.IdentityInterface(fields=["source_files", "t1w_dseg"]), + name="inputnode", ) - ds_fsnative_t1w = pe.Node( + outputnode = pe.Node(niu.IdentityInterface(fields=["t1w_dseg"]), name="outputnode") + + ds_t1w_dseg = pe.Node( DerivativesDataSink( base_directory=output_dir, - mode="image", - to="T1w", - suffix="xfm", - extension="txt", - **{"from": "fsnative"}, + suffix="dseg", + compress=True, + dismiss_entities=["desc"], ), - name="ds_fsnative_t1w", - run_without_submitting=True, - ) - # Surfaces - name_surfs = pe.MapNode( - Path2BIDS(), iterfield="in_file", name="name_surfs", run_without_submitting=True - ) - ds_surfs = pe.MapNode( - DerivativesDataSink(base_directory=output_dir, extension=".surf.gii"), - iterfield=["in_file", "hemi", "suffix"], - name="ds_surfs", + name="ds_t1w_dseg", run_without_submitting=True, ) - # Sphere registrations - name_regs = pe.MapNode( - Path2BIDS(), iterfield="in_file", name="name_regs", run_without_submitting=True - ) - ds_regs = pe.MapNode( - DerivativesDataSink( - base_directory=output_dir, - desc="reg", - suffix="sphere", - extension=".surf.gii", + + # fmt:off + workflow.connect([ + (inputnode, ds_t1w_dseg, [('t1w_dseg', 'in_file'), + ('source_files', 'source_file')]), + (ds_t1w_dseg, outputnode, [('out_file', 't1w_dseg')]), + ]) + # fmt:on + + return workflow + + +def init_ds_tpms_wf(*, output_dir, name="ds_tpms_wf", tpm_labels=BIDS_TISSUE_ORDER): + """ + Save tissue probability maps + + Parameters + ---------- + output_dir : :obj:`str` + Directory in which to save derivatives + name : :obj:`str` + Workflow name (default: anat_derivatives_wf) + tpm_labels : :obj:`tuple` + Tissue probability maps in order + + Inputs + ------ + source_files + List of input T1w images + t1w_tpms + Tissue probability maps in T1w space + + Outputs + ------- + t1w_tpms + The location in the output directory of the tissue probability maps + + """ + workflow = Workflow(name=name) + + inputnode = pe.Node( + niu.IdentityInterface(fields=["source_files", "t1w_tpms"]), + name="inputnode", + ) + outputnode = pe.Node(niu.IdentityInterface(fields=["t1w_tpms"]), name="outputnode") + + ds_t1w_tpms = pe.Node( + DerivativesDataSink( + base_directory=output_dir, + suffix="probseg", + compress=True, + dismiss_entities=["desc"], ), - iterfield=["in_file", "hemi"], - name="ds_regs", + name="ds_t1w_tpms", run_without_submitting=True, ) - name_reg_fsLR = pe.MapNode( - Path2BIDS(), iterfield="in_file", name="name_reg_fsLR", run_without_submitting=True + ds_t1w_tpms.inputs.label = tpm_labels + + # fmt:off + workflow.connect([ + (inputnode, ds_t1w_tpms, [('t1w_tpms', 'in_file'), + ('source_files', 'source_file')]), + (ds_t1w_tpms, outputnode, [("out_file", "t1w_tpms")]), + ]) + # fmt:on + + return workflow + + +def init_ds_template_registration_wf( + *, + output_dir, + name="ds_template_registration_wf", +): + """ + Save template registration transforms + + Parameters + ---------- + output_dir : :obj:`str` + Directory in which to save derivatives + name : :obj:`str` + Workflow name (default: anat_derivatives_wf) + + Inputs + ------ + template + Template space and specifications + source_files + List of input T1w images + anat2std_xfm + Nonlinear spatial transform to resample imaging data given in anatomical space + into standard space. + std2anat_xfm + Inverse transform of ``anat2std_xfm`` + + """ + workflow = Workflow(name=name) + + inputnode = pe.Node( + niu.IdentityInterface(fields=["template", "source_files", "anat2std_xfm", "std2anat_xfm"]), + name="inputnode", + ) + outputnode = pe.Node( + niu.IdentityInterface(fields=["anat2std_xfm", "std2anat_xfm"]), + name="outputnode", ) - ds_reg_fsLR = pe.MapNode( + + ds_std2t1w_xfm = pe.MapNode( DerivativesDataSink( base_directory=output_dir, - space="fsLR", - desc="reg", - suffix="sphere", - extension=".surf.gii", + to="T1w", + mode="image", + suffix="xfm", + dismiss_entities=["desc"], ), - iterfield=["in_file", "hemi"], - name="ds_reg_fsLR", + iterfield=("in_file", "from"), + name="ds_std2t1w_xfm", run_without_submitting=True, ) - # Morphometrics - name_morphs = pe.MapNode( - Path2BIDS(), - iterfield="in_file", - name="name_morphs", + + ds_t1w2std_xfm = pe.MapNode( + DerivativesDataSink( + base_directory=output_dir, + mode="image", + suffix="xfm", + dismiss_entities=["desc"], + **{"from": "T1w"}, + ), + iterfield=("in_file", "to"), + name="ds_t1w2std_xfm", run_without_submitting=True, ) - ds_morphs = pe.MapNode( - DerivativesDataSink(base_directory=output_dir, extension=".shape.gii"), - iterfield=["in_file", "hemi", "suffix"], - name="ds_morphs", + + # fmt:off + workflow.connect([ + (inputnode, ds_t1w2std_xfm, [ + ('anat2std_xfm', 'in_file'), + (('template', _combine_cohort), 'to'), + ('source_files', 'source_file')]), + (inputnode, ds_std2t1w_xfm, [ + ('std2anat_xfm', 'in_file'), + (('template', _combine_cohort), 'from'), + ('source_files', 'source_file')]), + (ds_t1w2std_xfm, outputnode, [("out_file", "anat2std_xfm")]), + (ds_std2t1w_xfm, outputnode, [("out_file", "std2anat_xfm")]), + ]) + # fmt:on + + return workflow + + +def init_ds_fs_registration_wf( + *, + output_dir, + name="ds_fs_registration_wf", +): + """ + Save rigid registration between subject anatomical template and FreeSurfer T1.mgz + + Parameters + ---------- + output_dir : :obj:`str` + Directory in which to save derivatives + name : :obj:`str` + Workflow name (default: ds_fs_registration_wf) + + Inputs + ------ + source_files + List of input T1w images + fsnative2t1w_xfm + LTA-style affine matrix translating from FreeSurfer-conformed + subject space to T1w + + Outputs + ------- + t1w2fsnative_xfm + LTA-style affine matrix translating from T1w to + FreeSurfer-conformed subject space + fsnative2t1w_xfm + LTA-style affine matrix translating from FreeSurfer-conformed + subject space to T1w + + """ + workflow = Workflow(name=name) + + inputnode = pe.Node( + niu.IdentityInterface(fields=["source_files", "fsnative2t1w_xfm"]), + name="inputnode", + ) + outputnode = pe.Node( + niu.IdentityInterface(fields=["fsnative2t1w_xfm", "t1w2fsnative_xfm"]), + name="outputnode", + ) + + from niworkflows.interfaces.nitransforms import ConcatenateXFMs + + # FS native space transforms + lta2itk = pe.Node(ConcatenateXFMs(inverse=True), name="lta2itk", run_without_submitting=True) + ds_t1w_fsnative = pe.Node( + DerivativesDataSink( + base_directory=output_dir, + mode="image", + to="fsnative", + suffix="xfm", + extension="txt", + **{"from": "T1w"}, + ), + name="ds_t1w_fsnative", run_without_submitting=True, ) - # Ribbon volume - ds_anat_ribbon = pe.Node( + ds_fsnative_t1w = pe.Node( DerivativesDataSink( base_directory=output_dir, - desc="ribbon", - suffix="mask", - extension=".nii.gz", + mode="image", + to="T1w", + suffix="xfm", + extension="txt", + **{"from": "fsnative"}, + ), + name="ds_fsnative_t1w", + run_without_submitting=True, + ) + + # fmt:off + workflow.connect([ + (inputnode, lta2itk, [('fsnative2t1w_xfm', 'in_xfms')]), + (inputnode, ds_t1w_fsnative, [('source_files', 'source_file')]), + (lta2itk, ds_t1w_fsnative, [('out_inv', 'in_file')]), + (inputnode, ds_fsnative_t1w, [('source_files', 'source_file')]), + (lta2itk, ds_fsnative_t1w, [('out_xfm', 'in_file')]), + (ds_fsnative_t1w, outputnode, [('out_file', 'fsnative2t1w_xfm')]), + (ds_t1w_fsnative, outputnode, [('out_file', 't1w2fsnative_xfm')]), + ]) + # fmt:on + return workflow + + +def init_ds_surfaces_wf( + *, + bids_root: str, + output_dir: str, + surfaces: list[str], + name="ds_surfaces_wf", +) -> Workflow: + """ + Save GIFTI surfaces + + Parameters + ---------- + bids_root : :class:`str` + Root path of BIDS dataset + output_dir : :class:`str` + Directory in which to save derivatives + surfaces : :class:`str` + List of surfaces to generate DataSinks for + name : :class:`str` + Workflow name (default: ds_surfaces_wf) + + Inputs + ------ + source_files + List of input T1w images + ```` + Left and right GIFTIs for each surface passed to ``surfaces`` + + Outputs + ------- + ```` + Left and right GIFTIs in ``output_dir`` for each surface passed to ``surfaces`` + + """ + workflow = Workflow(name=name) + + inputnode = pe.Node( + niu.IdentityInterface(fields=["source_files"] + surfaces), + name="inputnode", + ) + outputnode = pe.Node(niu.IdentityInterface(fields=surfaces), name="outputnode") + + for surf in surfaces: + ds_surf = pe.MapNode( + DerivativesDataSink( + base_directory=output_dir, + hemi=["L", "R"], + suffix=surf.split("_")[0], # Split for sphere_reg and sphere_reg_fsLR + extension=".surf.gii", + ), + iterfield=("in_file", "hemi"), + name=f"ds_{surf}", + run_without_submitting=True, + ) + if surf.startswith("sphere_reg"): + if surf == "sphere_reg_msm": + ds_surf.inputs.desc = "msmsulc" + ds_surf.inputs.space = "fsLR" + else: + ds_surf.inputs.desc = "reg" + if surf == "sphere_reg_fsLR": + ds_surf.inputs.space = "fsLR" + + # fmt:off + workflow.connect([ + (inputnode, ds_surf, [(surf, 'in_file'), ('source_files', 'source_file')]), + (ds_surf, outputnode, [('out_file', surf)]), + ]) + # fmt:on + + return workflow + + +def init_ds_surface_metrics_wf( + *, + bids_root: str, + output_dir: str, + metrics: list[str], + name="ds_surface_metrics_wf", +) -> Workflow: + """ + Save GIFTI surface metrics + + Parameters + ---------- + bids_root : :class:`str` + Root path of BIDS dataset + output_dir : :class:`str` + Directory in which to save derivatives + metrics : :class:`str` + List of metrics to generate DataSinks for + name : :class:`str` + Workflow name (default: ds_surface_metrics_wf) + + Inputs + ------ + source_files + List of input T1w images + ```` + Left and right GIFTIs for each metric passed to ``metrics`` + + Outputs + ------- + ```` + Left and right GIFTIs in ``output_dir`` for each metric passed to ``metrics`` + + """ + workflow = Workflow(name=name) + + inputnode = pe.Node( + niu.IdentityInterface(fields=["source_files"] + metrics), + name="inputnode", + ) + outputnode = pe.Node(niu.IdentityInterface(fields=metrics), name="outputnode") + + for metric in metrics: + ds_surf = pe.MapNode( + DerivativesDataSink( + base_directory=output_dir, + hemi=["L", "R"], + suffix=metric, + extension=".shape.gii", + ), + iterfield=("in_file", "hemi"), + name=f"ds_{metric}", + run_without_submitting=True, + ) + + # fmt:off + workflow.connect([ + (inputnode, ds_surf, [(metric, 'in_file'), ('source_files', 'source_file')]), + (ds_surf, outputnode, [('out_file', metric)]), + ]) + # fmt:on + + return workflow + + +def init_ds_grayord_metrics_wf( + *, + bids_root: str, + output_dir: str, + metrics: list[str], + cifti_output: ty.Literal["91k", "170k"], + name="ds_grayord_metrics_wf", +) -> Workflow: + """ + Save CIFTI-2 surface metrics + + Parameters + ---------- + bids_root : :class:`str` + Root path of BIDS dataset + output_dir : :class:`str` + Directory in which to save derivatives + metrics : :class:`str` + List of metrics to generate DataSinks for + cifti_output : :class:`str` + Density of CIFTI-2 files to save + name : :class:`str` + Workflow name (default: ds_surface_metrics_wf) + + Inputs + ------ + source_files + List of input T1w images + ```` + CIFTI-2 scalar file for each metric passed to ``metrics`` + ``_metadata`` + JSON file containing metadata for each metric passed to ``metrics`` + + Outputs + ------- + ```` + CIFTI-2 scalar file in ``output_dir`` for each metric passed to ``metrics`` + + """ + workflow = Workflow(name=name) + + inputnode = pe.Node( + niu.IdentityInterface( + fields=["source_files"] + metrics + [f"{m}_metadata" for m in metrics] + ), + name="inputnode", + ) + outputnode = pe.Node(niu.IdentityInterface(fields=metrics), name="outputnode") + + for metric in metrics: + ds_metric = pe.Node( + DerivativesDataSink( + base_directory=output_dir, + space='fsLR', + density=cifti_output, + suffix=metric, + compress=False, + extension=".dscalar.nii", + ), + name=f"ds_{metric}", + run_without_submitting=True, + ) + + workflow.connect([ + (inputnode, ds_metric, [ + ('source_files', 'source_file'), + (metric, 'in_file'), + ((f'{metric}_metadata', _read_json), 'meta_dict'), + ]), + (ds_metric, outputnode, [('out_file', metric)]), + ]) # fmt:skip + + return workflow + + +def init_ds_anat_volumes_wf( + *, + bids_root: str, + output_dir: str, + name="ds_anat_volumes_wf", + tpm_labels=BIDS_TISSUE_ORDER, +) -> pe.Workflow: + + workflow = pe.Workflow(name=name) + inputnode = pe.Node( + niu.IdentityInterface( + fields=[ + # Original T1w image + 'source_files', + # T1w-space images + 't1w_preproc', + 't1w_mask', + 't1w_dseg', + 't1w_tpms', + # Template + 'ref_file', + 'anat2std_xfm', + # Entities + 'space', + 'cohort', + 'resolution', + ] + ), + name='inputnode', + ) + + raw_sources = pe.Node(niu.Function(function=_bids_relative), name='raw_sources') + raw_sources.inputs.bids_root = bids_root + + gen_ref = pe.Node(GenerateSamplingReference(), name="gen_ref", mem_gb=0.01) + + # Mask T1w preproc images + mask_t1w = pe.Node(ApplyMask(), name='mask_t1w') + + # Resample T1w-space inputs + anat2std_t1w = pe.Node( + ApplyTransforms( + dimension=3, + default_value=0, + float=True, + interpolation="LanczosWindowedSinc", + ), + name="anat2std_t1w", + ) + + anat2std_mask = pe.Node(ApplyTransforms(interpolation="MultiLabel"), name="anat2std_mask") + anat2std_dseg = pe.Node(ApplyTransforms(interpolation="MultiLabel"), name="anat2std_dseg") + anat2std_tpms = pe.MapNode( + ApplyTransforms(dimension=3, default_value=0, float=True, interpolation="Gaussian"), + iterfield=["input_image"], + name="anat2std_tpms", + ) + + ds_std_t1w = pe.Node( + DerivativesDataSink( + base_directory=output_dir, + desc="preproc", compress=True, ), - name="ds_anat_ribbon", + name="ds_std_t1w", + run_without_submitting=True, + ) + ds_std_t1w.inputs.SkullStripped = True + + ds_std_mask = pe.Node( + DerivativesDataSink( + base_directory=output_dir, desc="brain", suffix="mask", compress=True + ), + name="ds_std_mask", run_without_submitting=True, ) + ds_std_mask.inputs.Type = "Brain" + + ds_std_dseg = pe.Node( + DerivativesDataSink(base_directory=output_dir, suffix="dseg", compress=True), + name="ds_std_dseg", + run_without_submitting=True, + ) + + ds_std_tpms = pe.Node( + DerivativesDataSink(base_directory=output_dir, suffix="probseg", compress=True), + name="ds_std_tpms", + run_without_submitting=True, + ) + + # CRITICAL: the sequence of labels here (CSF-GM-WM) is that of the output of FSL-FAST + # (intensity mean, per tissue). This order HAS to be matched also by the ``tpms`` + # output in the data/io_spec.json file. + ds_std_tpms.inputs.label = tpm_labels + + workflow.connect([ + (inputnode, gen_ref, [ + ("ref_file", "fixed_image"), + (("resolution", _is_native), "keep_native"), + ]), + (inputnode, mask_t1w, [ + ("t1w_preproc", "in_file"), + ("t1w_mask", "in_mask"), + ]), + (mask_t1w, anat2std_t1w, [("out_file", "input_image")]), + (inputnode, anat2std_mask, [("t1w_mask", "input_image")]), + (inputnode, anat2std_dseg, [("t1w_dseg", "input_image")]), + (inputnode, anat2std_tpms, [("t1w_tpms", "input_image")]), + (inputnode, gen_ref, [("t1w_preproc", "moving_image")]), + (anat2std_t1w, ds_std_t1w, [("output_image", "in_file")]), + (anat2std_mask, ds_std_mask, [("output_image", "in_file")]), + (anat2std_dseg, ds_std_dseg, [("output_image", "in_file")]), + (anat2std_tpms, ds_std_tpms, [("output_image", "in_file")]), + ]) # fmt:skip + + workflow.connect( + # Connect apply transforms nodes + [ + (gen_ref, n, [('out_file', 'reference_image')]) + for n in (anat2std_t1w, anat2std_mask, anat2std_dseg, anat2std_tpms) + ] + + [ + (inputnode, n, [('anat2std_xfm', 'transforms')]) + for n in (anat2std_t1w, anat2std_mask, anat2std_dseg, anat2std_tpms) + ] + + [ + (inputnode, n, [ + ('source_files', 'source_file'), + ('space', 'space'), + ('cohort', 'cohort'), + ('resolution', 'resolution'), + ]) + for n in (ds_std_t1w, ds_std_mask, ds_std_dseg, ds_std_tpms) + ] + ) # fmt:skip + + return workflow + + +def init_anat_second_derivatives_wf( + *, + bids_root: str, + output_dir: str, + cifti_output: ty.Literal["91k", "170k", False], + name="anat_second_derivatives_wf", + tpm_labels=BIDS_TISSUE_ORDER, +): + """ + Set up a battery of datasinks to store derivatives in the right location. + + Parameters + ---------- + bids_root : :obj:`str` + Root path of BIDS dataset + output_dir : :obj:`str` + Directory in which to save derivatives + name : :obj:`str` + Workflow name (default: anat_derivatives_wf) + tpm_labels : :obj:`tuple` + Tissue probability maps in order + + Inputs + ------ + template + Template space and specifications + source_files + List of input T1w images + t1w_preproc + The T1w reference map, which is calculated as the average of bias-corrected + and preprocessed T1w images, defining the anatomical space. + t1w_mask + Mask of the ``t1w_preproc`` + t1w_dseg + Segmentation in T1w space + t1w_tpms + Tissue probability maps in T1w space + anat2std_xfm + Nonlinear spatial transform to resample imaging data given in anatomical space + into standard space. + surfaces + GIFTI surfaces (gray/white boundary, midthickness, pial, inflated) + morphometrics + GIFTIs of cortical thickness, curvature, and sulcal depth + t1w_fs_aseg + FreeSurfer's aseg segmentation, in native T1w space + t1w_fs_aparc + FreeSurfer's aparc+aseg segmentation, in native T1w space + cifti_morph + Morphometric CIFTI-2 dscalar files + cifti_metadata + JSON files containing metadata dictionaries + + """ + workflow = Workflow(name=name) + + inputnode = pe.Node( + niu.IdentityInterface( + fields=[ + "template", + "source_files", + "t1w_fs_aseg", + "t1w_fs_aparc", + ] + ), + name="inputnode", + ) + + raw_sources = pe.Node(niu.Function(function=_bids_relative), name="raw_sources") + raw_sources.inputs.bids_root = bids_root # Parcellations ds_t1w_fsaseg = pe.Node( @@ -684,60 +1104,98 @@ def init_anat_derivatives_wf( # fmt:off workflow.connect([ - (inputnode, lta2itk_fwd, [('t1w2fsnative_xfm', 'in_xfms')]), - (inputnode, lta2itk_inv, [('fsnative2t1w_xfm', 'in_xfms')]), - (inputnode, ds_t1w_fsnative, [('source_files', 'source_file')]), - (lta2itk_fwd, ds_t1w_fsnative, [('out_xfm', 'in_file')]), - (inputnode, ds_fsnative_t1w, [('source_files', 'source_file')]), - (lta2itk_inv, ds_fsnative_t1w, [('out_xfm', 'in_file')]), - (inputnode, name_surfs, [('surfaces', 'in_file')]), - (inputnode, ds_surfs, [('surfaces', 'in_file'), - ('source_files', 'source_file')]), - (name_surfs, ds_surfs, [('hemi', 'hemi'), - ('suffix', 'suffix')]), - (inputnode, name_regs, [('sphere_reg', 'in_file')]), - (inputnode, ds_regs, [('sphere_reg', 'in_file'), - ('source_files', 'source_file')]), - (name_regs, ds_regs, [('hemi', 'hemi')]), - (inputnode, name_reg_fsLR, [('sphere_reg_fsLR', 'in_file')]), - (inputnode, ds_reg_fsLR, [('sphere_reg_fsLR', 'in_file'), - ('source_files', 'source_file')]), - (name_reg_fsLR, ds_reg_fsLR, [('hemi', 'hemi')]), - (inputnode, name_morphs, [('morphometrics', 'in_file')]), - (inputnode, ds_morphs, [('morphometrics', 'in_file'), - ('source_files', 'source_file')]), - (name_morphs, ds_morphs, [('hemi', 'hemi'), - ('suffix', 'suffix')]), (inputnode, ds_t1w_fsaseg, [('t1w_fs_aseg', 'in_file'), ('source_files', 'source_file')]), (inputnode, ds_t1w_fsparc, [('t1w_fs_aparc', 'in_file'), ('source_files', 'source_file')]), - (inputnode, ds_anat_ribbon, [('anat_ribbon', 'in_file'), - ('source_files', 'source_file')]), + ]) + # fmt:on + return workflow + + +def init_template_iterator_wf(*, spaces, name="template_iterator_wf"): + """Prepare the necessary components to resample an image to a template space + + This produces a workflow with an unjoined iterable named "spacesource". + + It takes as input a collated list of template specifiers and transforms to that + space. + + The fields in `outputnode` can be used as if they come from a single template. + """ + workflow = pe.Workflow(name=name) + + inputnode = pe.Node( + niu.IdentityInterface(fields=["template", "anat2std_xfm"]), + name="inputnode", + ) + outputnode = pe.Node( + niu.IdentityInterface( + fields=[ + "space", + "resolution", + "cohort", + "anat2std_xfm", + "std_t1w", + "std_mask", + ], + ), + name="outputnode", + ) + + spacesource = pe.Node(SpaceDataSource(), name="spacesource", run_without_submitting=True) + spacesource.iterables = ( + "in_tuple", + [(s.fullname, s.spec) for s in spaces.cached.get_standard(dim=(3,))], + ) + + gen_tplid = pe.Node( + niu.Function(function=_fmt_cohort), + name="gen_tplid", + run_without_submitting=True, + ) + + select_xfm = pe.Node( + KeySelect(fields=["anat2std_xfm"]), + name="select_xfm", + run_without_submitting=True, + ) + select_tpl = pe.Node( + TemplateFlowSelect(resolution=1), name="select_tpl", run_without_submitting=True + ) + + # fmt:off + workflow.connect([ + (inputnode, select_xfm, [ + ('anat2std_xfm', 'anat2std_xfm'), + ('template', 'keys'), + ]), + (spacesource, gen_tplid, [ + ('space', 'template'), + ('cohort', 'cohort'), + ]), + (gen_tplid, select_xfm, [('out', 'key')]), + (spacesource, select_tpl, [ + ('space', 'template'), + ('cohort', 'cohort'), + (('resolution', _no_native), 'resolution'), + ]), + (spacesource, outputnode, [ + ("space", "space"), + ("resolution", "resolution"), + ("cohort", "cohort"), + ]), + (select_xfm, outputnode, [ + ("anat2std_xfm", "anat2std_xfm"), + ]), + (select_tpl, outputnode, [ + ("t1w_file", "std_t1w"), + ("brain_mask", "std_mask"), + ]), ]) # fmt:on - if cifti_output: - ds_cifti_morph = pe.MapNode( - DerivativesDataSink( - base_directory=output_dir, - suffix=['curv', 'sulc', 'thickness'], - compress=False, - space='fsLR', - ), - name='ds_cifti_morph', - run_without_submitting=True, - iterfield=["in_file", "meta_dict", "suffix"], - ) - # fmt:off - workflow.connect([ - (inputnode, ds_cifti_morph, [('cifti_morph', 'in_file'), - ('source_files', 'source_file'), - ('cifti_density', 'density'), - (('cifti_metadata', _read_jsons), 'meta_dict')]) - ]) - # fmt:on return workflow @@ -746,7 +1204,12 @@ def _bids_relative(in_files, bids_root): if not isinstance(in_files, (list, tuple)): in_files = [in_files] - in_files = [str(Path(p).relative_to(bids_root)) for p in in_files] + ret = [] + for file in in_files: + try: + ret.append(str(Path(file).relative_to(bids_root))) + except ValueError: + ret.append(file) return in_files @@ -836,8 +1299,8 @@ def _combine_cohort(in_template): return [_combine_cohort(v) for v in in_template] -def _read_jsons(in_file): +def _read_json(in_file): from json import loads from pathlib import Path - return [loads(Path(f).read_text()) for f in in_file] + return loads(Path(in_file).read_text()) diff --git a/smriprep/workflows/surfaces.py b/smriprep/workflows/surfaces.py index 450ecae6e8..71fe9d3557 100644 --- a/smriprep/workflows/surfaces.py +++ b/smriprep/workflows/surfaces.py @@ -28,31 +28,42 @@ """ import typing as ty -from nipype.pipeline import engine as pe + +from nipype.interfaces import freesurfer as fs +from nipype.interfaces import io as nio +from nipype.interfaces import utility as niu from nipype.interfaces.base import Undefined -from nipype.interfaces import ( - fsl, - io as nio, - utility as niu, - freesurfer as fs, - workbench as wb, +from nipype.pipeline import engine as pe +from niworkflows.engine.workflows import LiterateWorkflow as Workflow +from niworkflows.interfaces.freesurfer import FSDetectInputs, FSInjectBrainExtracted +from niworkflows.interfaces.freesurfer import PatchedRobustRegister as RobustRegister +from niworkflows.interfaces.freesurfer import RefineBrainMask +from niworkflows.interfaces.nitransforms import ConcatenateXFMs +from niworkflows.interfaces.utility import KeySelect +from niworkflows.interfaces.workbench import ( + MetricDilate, + MetricFillHoles, + MetricMask, + MetricRemoveIslands, + MetricResample, ) -from ..data import load_resource -from ..interfaces.freesurfer import ReconAll, MakeMidthickness +from smriprep.interfaces.surf import MakeRibbon +from smriprep.interfaces.workbench import SurfaceResample -from niworkflows.engine.workflows import LiterateWorkflow as Workflow -from niworkflows.interfaces.freesurfer import ( - FSDetectInputs, - FSInjectBrainExtracted, - PatchedLTAConvert as LTAConvert, - PatchedRobustRegister as RobustRegister, - RefineBrainMask, -) +from ..data import load_resource +from ..interfaces.freesurfer import MakeMidthickness, ReconAll +from ..interfaces.gifti import MetricMath from ..interfaces.workbench import CreateSignedDistanceVolume -def init_surface_recon_wf(*, omp_nthreads, hires, name="surface_recon_wf"): +def init_surface_recon_wf( + *, + omp_nthreads: int, + hires: bool, + precomputed: dict, + name="surface_recon_wf", +): r""" Reconstruct anatomical surfaces using FreeSurfer's ``recon-all``. @@ -115,7 +126,7 @@ def init_surface_recon_wf(*, omp_nthreads, hires, name="surface_recon_wf"): :simple_form: yes from smriprep.workflows.surfaces import init_surface_recon_wf - wf = init_surface_recon_wf(omp_nthreads=1, hires=True) + wf = init_surface_recon_wf(omp_nthreads=1, hires=True, precomputed={}) Parameters ---------- @@ -134,10 +145,6 @@ def init_surface_recon_wf(*, omp_nthreads, hires, name="surface_recon_wf"): List of FLAIR images skullstripped_t1 Skull-stripped T1-weighted image (or mask of image) - ants_segs - Brain tissue segmentation from ANTS ``antsBrainExtraction.sh`` - corrected_t1 - INU-corrected, merged T1-weighted image subjects_dir FreeSurfer SUBJECTS_DIR subject_id @@ -149,26 +156,12 @@ def init_surface_recon_wf(*, omp_nthreads, hires, name="surface_recon_wf"): FreeSurfer SUBJECTS_DIR subject_id FreeSurfer subject ID - t1w2fsnative_xfm - LTA-style affine matrix translating from T1w to FreeSurfer-conformed subject space fsnative2t1w_xfm LTA-style affine matrix translating from FreeSurfer-conformed subject space to T1w - surfaces - GIFTI surfaces for gray/white matter boundary, pial surface, - midthickness (or graymid) surface, and inflated surfaces - out_brainmask - Refined brainmask, derived from FreeSurfer's ``aseg`` volume - out_aseg - FreeSurfer's aseg segmentation, in native T1w space - out_aparc - FreeSurfer's aparc+aseg segmentation, in native T1w space - morphometrics - GIFTIs of cortical thickness, curvature, and sulcal depth See also -------- * :py:func:`~smriprep.workflows.surfaces.init_autorecon_resume_wf` - * :py:func:`~smriprep.workflows.surfaces.init_gifti_surface_wf` """ workflow = Workflow(name=name) @@ -189,8 +182,6 @@ def init_surface_recon_wf(*, omp_nthreads, hires, name="surface_recon_wf"): "t2w", "flair", "skullstripped_t1", - "corrected_t1", - "ants_segs", "subjects_dir", "subject_id", ] @@ -204,18 +195,6 @@ def init_surface_recon_wf(*, omp_nthreads, hires, name="surface_recon_wf"): "subject_id", "t1w2fsnative_xfm", "fsnative2t1w_xfm", - "surfaces", - "out_brainmask", - "out_aseg", - "out_aparc", - "morphometrics", - "midthickness", - "pial", - "white", - "inflated", - "thickness", - "sulc", - "curv", ] ), name="outputnode", @@ -237,19 +216,27 @@ def init_surface_recon_wf(*, omp_nthreads, hires, name="surface_recon_wf"): 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" + autorecon_resume_wf = init_autorecon_resume_wf(omp_nthreads=omp_nthreads) + + get_surfaces = pe.Node(nio.FreeSurferSource(), name="get_surfaces") + + midthickness = pe.MapNode( + MakeMidthickness(thickness=True, distance=0.5, out_name="midthickness"), + iterfield="in_file", + name="midthickness", + n_procs=min(omp_nthreads, 12), ) - t1w2fsnative_xfm = pe.Node(LTAConvert(out_lta=True, invert=True), name="t1w2fsnative_xfm") - autorecon_resume_wf = init_autorecon_resume_wf(omp_nthreads=omp_nthreads) - gifti_surface_wf = init_gifti_surface_wf() + save_midthickness = pe.Node(nio.DataSink(parameterization=False), name="save_midthickness") - aseg_to_native_wf = init_segs_to_native_wf() - aparc_to_native_wf = init_segs_to_native_wf(segmentation="aparc_aseg") - refine = pe.Node(RefineBrainMask(), name="refine") + sync = pe.Node( + niu.Function( + function=_extract_fs_fields, + output_names=['subjects_dir', 'subject_id'], + ), + name="sync", + ) - # fmt:off workflow.connect([ # Configuration (inputnode, recon_config, [('t1w', 't1w_list'), @@ -262,9 +249,6 @@ def init_surface_recon_wf(*, omp_nthreads, hires, name="surface_recon_wf"): ('subject_id', 'subject_id')]), (skull_strip_extern, autorecon_resume_wf, [('subjects_dir', 'inputnode.subjects_dir'), ('subject_id', 'inputnode.subject_id')]), - (autorecon_resume_wf, gifti_surface_wf, [ - ('outputnode.subjects_dir', 'inputnode.subjects_dir'), - ('outputnode.subject_id', 'inputnode.subject_id')]), # Reconstruction phases (inputnode, autorecon1, [('t1w', 'T1_files')]), (inputnode, fov_check, [('t1w', 'in_files')]), @@ -277,43 +261,143 @@ def init_surface_recon_wf(*, omp_nthreads, hires, name="surface_recon_wf"): (inputnode, skull_strip_extern, [('skullstripped_t1', 'in_brain')]), (recon_config, autorecon_resume_wf, [('use_t2w', 'inputnode.use_T2'), ('use_flair', 'inputnode.use_FLAIR')]), - # Construct transform from FreeSurfer conformed image to sMRIPrep - # reoriented image - (inputnode, fsnative2t1w_xfm, [('t1w', 'target_file')]), - (autorecon1, fsnative2t1w_xfm, [('T1', 'source_file')]), - (fsnative2t1w_xfm, t1w2fsnative_xfm, [('out_reg_file', 'in_lta')]), + # Generate mid-thickness surfaces + (autorecon_resume_wf, get_surfaces, [ + ('outputnode.subjects_dir', 'subjects_dir'), + ('outputnode.subject_id', 'subject_id'), + ]), + (autorecon_resume_wf, save_midthickness, [ + ('outputnode.subjects_dir', 'base_directory'), + ('outputnode.subject_id', 'container'), + ]), + (get_surfaces, midthickness, [ + ('white', 'in_file'), + ('graymid', 'graymid'), + ]), + (midthickness, save_midthickness, [('out_file', 'surf.@graymid')]), + # Output + (save_midthickness, sync, [('out_file', 'filenames')]), + (sync, outputnode, [('subjects_dir', 'subjects_dir'), + ('subject_id', 'subject_id')]), + ]) # fmt:skip + + if "fsnative" not in precomputed.get("transforms", {}): + fsnative2t1w_xfm = pe.Node( + RobustRegister(auto_sens=True, est_int_scale=True), name="fsnative2t1w_xfm" + ) + + workflow.connect([ + (inputnode, fsnative2t1w_xfm, [('t1w', 'target_file')]), + (autorecon1, fsnative2t1w_xfm, [('T1', 'source_file')]), + (fsnative2t1w_xfm, outputnode, [('out_reg_file', 'fsnative2t1w_xfm')]), + ]) # fmt:skip + + return workflow + + +def init_refinement_wf(*, name="refinement_wf"): + r""" + Refine ANTs brain extraction with FreeSurfer segmentation + + Workflow Graph + .. workflow:: + :graph2use: orig + :simple_form: yes + + from smriprep.workflows.surfaces import init_refinement_wf + wf = init_refinement_wf() + + Inputs + ------ + subjects_dir + FreeSurfer SUBJECTS_DIR + subject_id + FreeSurfer subject ID + fsnative2t1w_xfm + LTA-style affine matrix translating from FreeSurfer-conformed subject space to T1w + reference_image + Input + t2w + List of T2-weighted structural images (only first used) + flair + List of FLAIR images + skullstripped_t1 + Skull-stripped T1-weighted image (or mask of image) + ants_segs + Brain tissue segmentation from ANTS ``antsBrainExtraction.sh`` + corrected_t1 + INU-corrected, merged T1-weighted image + subjects_dir + FreeSurfer SUBJECTS_DIR + subject_id + FreeSurfer subject ID + + Outputs + ------- + subjects_dir + FreeSurfer SUBJECTS_DIR + subject_id + FreeSurfer subject ID + t1w2fsnative_xfm + LTA-style affine matrix translating from T1w to FreeSurfer-conformed subject space + fsnative2t1w_xfm + LTA-style affine matrix translating from FreeSurfer-conformed subject space to T1w + out_brainmask + Refined brainmask, derived from FreeSurfer's ``aseg`` volume + + See also + -------- + * :py:func:`~smriprep.workflows.surfaces.init_autorecon_resume_wf` + + """ + workflow = Workflow(name=name) + workflow.__desc__ = """\ +Brain surfaces were reconstructed using `recon-all` [FreeSurfer {fs_ver}, +RRID:SCR_001847, @fs_reconall], and the brain mask estimated +previously was refined with a custom variation of the method to reconcile +ANTs-derived and FreeSurfer-derived segmentations of the cortical +gray-matter of Mindboggle [RRID:SCR_002438, @mindboggle]. +""".format( + fs_ver=fs.Info().looseversion() or "" + ) + + inputnode = pe.Node( + niu.IdentityInterface( + fields=[ + "reference_image", + "ants_segs", + "fsnative2t1w_xfm", + "subjects_dir", + "subject_id", + ] + ), + name="inputnode", + ) + outputnode = pe.Node( + niu.IdentityInterface( + fields=[ + "out_brainmask", + ] + ), + name="outputnode", + ) + + aseg_to_native_wf = init_segs_to_native_wf() + refine = pe.Node(RefineBrainMask(), name="refine") + + # fmt:off + workflow.connect([ # Refine ANTs mask, deriving new mask from FS' aseg - (inputnode, refine, [('corrected_t1', 'in_anat'), + (inputnode, aseg_to_native_wf, [ + ('subjects_dir', 'inputnode.subjects_dir'), + ('subject_id', 'inputnode.subject_id'), + ('reference_image', 'inputnode.in_file'), + ('fsnative2t1w_xfm', 'inputnode.fsnative2t1w_xfm'), + ]), + (inputnode, refine, [('reference_image', 'in_anat'), ('ants_segs', 'in_ants')]), - (inputnode, aseg_to_native_wf, [('corrected_t1', 'inputnode.in_file')]), - (autorecon_resume_wf, aseg_to_native_wf, [ - ('outputnode.subjects_dir', 'inputnode.subjects_dir'), - ('outputnode.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')]), - (autorecon_resume_wf, aparc_to_native_wf, [ - ('outputnode.subjects_dir', 'inputnode.subjects_dir'), - ('outputnode.subject_id', 'inputnode.subject_id')]), - (fsnative2t1w_xfm, aparc_to_native_wf, [('out_reg_file', 'inputnode.fsnative2t1w_xfm')]), (aseg_to_native_wf, refine, [('outputnode.out_file', 'in_aseg')]), - - # Output - (autorecon_resume_wf, outputnode, [('outputnode.subjects_dir', 'subjects_dir'), - ('outputnode.subject_id', 'subject_id')]), - (gifti_surface_wf, outputnode, [('outputnode.surfaces', 'surfaces'), - ('outputnode.morphometrics', 'morphometrics'), - ('outputnode.midthickness', 'midthickness'), - ('outputnode.pial', 'pial'), - ('outputnode.white', 'white'), - ('outputnode.inflated', 'inflated'), - ('outputnode.thickness', 'thickness'), - ('outputnode.sulc', 'sulc'), - ('outputnode.curv', 'curv')]), - (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 @@ -493,36 +577,124 @@ def _dedup(in_list): return workflow -def init_sphere_reg_wf( +def init_surface_derivatives_wf( *, - msm_sulc: bool = False, - sloppy: bool = False, - name: str = "sphere_reg_wf", + cifti_output: ty.Literal["91k", "170k", False] = False, + name="surface_derivatives_wf", ): - """Generate GIFTI registration files to fsLR space""" - from ..interfaces.surf import FixGiftiMetadata - from ..interfaces.workbench import SurfaceSphereProjectUnproject + r""" + Generate sMRIPrep derivatives from FreeSurfer derivatives + Workflow Graph + .. workflow:: + :graph2use: orig + :simple_form: yes + + from smriprep.workflows.surfaces import init_surface_derivatives_wf + wf = init_surface_derivatives_wf() + + Inputs + ------ + reference + Reference image in native T1w space, for defining a resampling grid + fsnative2t1w_xfm + LTA-style affine matrix translating from FreeSurfer-conformed subject space to T1w + subjects_dir + FreeSurfer SUBJECTS_DIR + subject_id + FreeSurfer subject ID + + Outputs + ------- + surfaces + GIFTI surfaces for gray/white matter boundary, pial surface, + midthickness (or graymid) surface, and inflated surfaces + morphometrics + GIFTIs of cortical thickness, curvature, and sulcal depth + out_aseg + FreeSurfer's aseg segmentation, in native T1w space + out_aparc + FreeSurfer's aparc+aseg segmentation, in native T1w space + + See also + -------- + * :py:func:`~smriprep.workflows.surfaces.init_gifti_surface_wf` + + """ workflow = Workflow(name=name) inputnode = pe.Node( - niu.IdentityInterface(["subjects_dir", "subject_id", "sulc"]), + niu.IdentityInterface( + fields=[ + "subjects_dir", + "subject_id", + "fsnative2t1w_xfm", + "reference", + ] + ), name="inputnode", ) outputnode = pe.Node( - niu.IdentityInterface(["sphere_reg", "sphere_reg_fsLR"]), name="outputnode" + niu.IdentityInterface( + fields=[ + "inflated", + "curv", + "out_aseg", + "out_aparc", + ] + ), + name="outputnode", ) - get_surfaces = pe.Node(nio.FreeSurferSource(), name="get_surfaces") + gifti_surfaces_wf = init_gifti_surfaces_wf(surfaces=["inflated"]) + gifti_morph_wf = init_gifti_morphometrics_wf(morphometrics=["curv"]) + aseg_to_native_wf = init_segs_to_native_wf() + aparc_to_native_wf = init_segs_to_native_wf(segmentation="aparc_aseg") - # Via FreeSurfer2CaretConvertAndRegisterNonlinear.sh#L270-L273 - # - # See https://github.com/DCAN-Labs/DCAN-HCP/tree/9291324 - sphere_reg_gii = pe.MapNode( - fs.MRIsConvert(out_datatype="gii"), iterfield="in_file", name="sphere_reg_gii" - ) + # fmt:off + workflow.connect([ + # Configuration + (inputnode, gifti_surfaces_wf, [ + ('subjects_dir', 'inputnode.subjects_dir'), + ('subject_id', 'inputnode.subject_id'), + ('fsnative2t1w_xfm', 'inputnode.fsnative2t1w_xfm'), + ]), + (inputnode, gifti_morph_wf, [ + ('subjects_dir', 'inputnode.subjects_dir'), + ('subject_id', 'inputnode.subject_id'), + ]), + (inputnode, aseg_to_native_wf, [ + ('subjects_dir', 'inputnode.subjects_dir'), + ('subject_id', 'inputnode.subject_id'), + ('reference', 'inputnode.in_file'), + ('fsnative2t1w_xfm', 'inputnode.fsnative2t1w_xfm'), + ]), + (inputnode, aparc_to_native_wf, [ + ('subjects_dir', 'inputnode.subjects_dir'), + ('subject_id', 'inputnode.subject_id'), + ('reference', 'inputnode.in_file'), + ('fsnative2t1w_xfm', 'inputnode.fsnative2t1w_xfm'), + ]), + + # Output + (gifti_surfaces_wf, outputnode, [('outputnode.inflated', 'inflated')]), + (aseg_to_native_wf, outputnode, [('outputnode.out_file', 'out_aseg')]), + (aparc_to_native_wf, outputnode, [('outputnode.out_file', 'out_aparc')]), + (gifti_morph_wf, outputnode, [('outputnode.curv', 'curv')]), + ]) + # fmt:on + + return workflow + + +def init_fsLR_reg_wf(*, name="fsLR_reg_wf"): + """Generate GIFTI registration files to fsLR space""" + from ..interfaces.workbench import SurfaceSphereProjectUnproject + + workflow = Workflow(name=name) - fix_reg_meta = pe.MapNode(FixGiftiMetadata(), iterfield="in_file", name="fix_reg_meta") + inputnode = pe.Node(niu.IdentityInterface(["sphere_reg", "sulc"]), name="inputnode") + outputnode = pe.Node(niu.IdentityInterface(["sphere_reg_fsLR"]), name="outputnode") # Via # ${CARET7DIR}/wb_command -surface-sphere-project-unproject @@ -547,44 +719,16 @@ def init_sphere_reg_wf( # fmt:off workflow.connect([ - (inputnode, get_surfaces, [ - ('subjects_dir', 'subjects_dir'), - ('subject_id', 'subject_id'), - ]), - (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')]), + (inputnode, project_unproject, [('sphere_reg', 'sphere_in')]), + (project_unproject, outputnode, [('sphere_out', 'sphere_reg_fsLR')]), ]) # 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(sloppy=sloppy) - # 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(*, sloppy: bool = False, name: str = 'msm_sulc_wf'): """Run MSMSulc registration to fsLR surfaces, per hemisphere.""" - from ..interfaces.gifti import InvertShape from ..interfaces.msm import MSM from ..interfaces.workbench import ( SurfaceAffineRegression, @@ -631,8 +775,8 @@ def init_msm_sulc_wf(*, sloppy: bool = False, name: str = 'msm_sulc_wf'): # 2) Invert sulc # wb_command -metric-math "-1 * var" ... invert_sulc = pe.MapNode( - InvertShape(shape='sulc'), - iterfield=['shape_file', 'hemisphere'], + MetricMath(metric='sulc', operation='invert'), + iterfield=['metric_file', 'hemisphere'], name='invert_sulc', ) invert_sulc.inputs.hemisphere = ['L', 'R'] @@ -665,35 +809,39 @@ def init_msm_sulc_wf(*, sloppy: bool = False, name: str = 'msm_sulc_wf'): ('sphere_reg_fsLR', 'target_surface'), ]), (inputnode, apply_surface_affine, [('sphere', 'in_surface')]), - (inputnode, invert_sulc, [('sulc', 'shape_file')]), + (inputnode, invert_sulc, [('sulc', 'metric_file')]), (regress_affine, apply_surface_affine, [('out_affine', 'in_affine')]), (apply_surface_affine, modify_sphere, [('out_surface', 'in_surface')]), (modify_sphere, msmsulc, [('out_surface', 'in_mesh')]), - (invert_sulc, msmsulc, [('shape_file', 'in_data')]), + (invert_sulc, msmsulc, [('metric_file', 'in_data')]), (msmsulc, outputnode, [('warped_mesh', 'sphere_reg_fsLR')]), ]) # fmt:skip return workflow -def init_gifti_surface_wf(*, name="gifti_surface_wf"): +def init_gifti_surfaces_wf( + *, + surfaces: ty.List[str] = ["pial", "midthickness", "inflated", "white"], + to_scanner: bool = True, + name: str = "gifti_surface_wf", +): r""" Prepare GIFTI surfaces from a FreeSurfer subjects directory. - If midthickness (or graymid) surfaces do not exist, they are generated and - saved to the subject directory as ``lh/rh.midthickness``. - These, along with the gray/white matter boundary (``lh/rh.white``), pial - surfaces (``lh/rh.pial``) and inflated surfaces (``lh/rh.inflated``) are - converted to GIFTI files. - Additionally, the vertex coordinates are :py:class:`recentered - ` to align with native T1w space. + The default surfaces are ``lh/rh.pial``, ``lh/rh.midthickness``, + ``lh/rh.inflated``, and ``lh/rh.white``. + + Vertex coordinates are :py:class:`transformed + ` to align with native T1w space + when ``fsnative2t1w_xfm`` is provided. Workflow Graph .. workflow:: :graph2use: orig :simple_form: yes - from smriprep.workflows.surfaces import init_gifti_surface_wf - wf = init_gifti_surface_wf() + from smriprep.workflows.surfaces import init_gifti_surfaces_wf + wf = init_gifti_surfaces_wf() Inputs ------ @@ -702,33 +850,17 @@ def init_gifti_surface_wf(*, name="gifti_surface_wf"): subject_id FreeSurfer subject ID fsnative2t1w_xfm - LTA formatted affine transform file (inverse) + LTA formatted affine transform file Outputs ------- - midthickness - Left and right midthickness (or graymid) surface GIFTIs - pial - Left and right pial surface GIFTIs - white - Left and right white surface GIFTIs - inflated - Left and right inflated surface GIFTIs surfaces - GIFTI surfaces for gray/white matter boundary, pial surface, - midthickness (or graymid) surface, and inflated surfaces - thickness - Left and right cortical thickness GIFTIs - sulc - Left and right sulcal depth map GIFTIs - curv - Left and right curvature map GIFTIs - morphometrics - GIFTIs of cortical thickness, curvature, and sulcal depth + GIFTI surfaces for all requested surfaces + ```` + Left and right GIFTIs for each surface passed to ``surfaces`` """ - from ..interfaces.freesurfer import MRIsConvertData - from ..interfaces.surf import NormalizeSurf, AggregateSurfaces + from ..interfaces.surf import NormalizeSurf workflow = Workflow(name=name) @@ -736,45 +868,113 @@ def init_gifti_surface_wf(*, name="gifti_surface_wf"): niu.IdentityInterface(["subjects_dir", "subject_id", "fsnative2t1w_xfm"]), name="inputnode", ) - outputnode = pe.Node( - niu.IdentityInterface([ - "pial", - "white", - "inflated", - "midthickness", - "thickness", - "sulc", - "curv", - # Preserve grouping - "surfaces", - "morphometrics", - ]), - name="outputnode" - ) + outputnode = pe.Node(niu.IdentityInterface(["surfaces", *surfaces]), name="outputnode") - get_surfaces = pe.Node(nio.FreeSurferSource(), name="get_surfaces") - - midthickness = pe.MapNode( - MakeMidthickness(thickness=True, distance=0.5, out_name="midthickness"), - iterfield="in_file", - name="midthickness", + get_surfaces = pe.Node( + niu.Function(function=_get_surfaces, output_names=surfaces), + name="get_surfaces", ) - - save_midthickness = pe.Node(nio.DataSink(parameterization=False), name="save_midthickness") + get_surfaces.inputs.surfaces = surfaces surface_list = pe.Node( - niu.Merge(4, ravel_inputs=True), + niu.Merge(len(surfaces), ravel_inputs=True), name="surface_list", run_without_submitting=True, ) fs2gii = pe.MapNode( - fs.MRIsConvert(out_datatype="gii", to_scanner=True), + fs.MRIsConvert(out_datatype="gii", to_scanner=to_scanner), iterfield="in_file", name="fs2gii", ) fix_surfs = pe.MapNode(NormalizeSurf(), iterfield="in_file", name="fix_surfs") - surfmorph_list = pe.Node( - niu.Merge(3, ravel_inputs=True), + + surface_groups = pe.Node( + niu.Split(splits=[2] * len(surfaces)), + name="surface_groups", + run_without_submitting=True, + ) + + workflow.connect([ + (inputnode, get_surfaces, [ + ('subjects_dir', 'subjects_dir'), + ('subject_id', 'subject_id'), + ]), + (inputnode, fix_surfs, [ + ('fsnative2t1w_xfm', 'transform_file'), + ]), + (get_surfaces, surface_list, [ + (surf, f'in{i}') for i, surf in enumerate(surfaces, start=1) + ]), + (surface_list, fs2gii, [('out', 'in_file')]), + (fs2gii, fix_surfs, [('converted', 'in_file')]), + (fix_surfs, outputnode, [('out_file', 'surfaces')]), + (fix_surfs, surface_groups, [('out_file', 'inlist')]), + (surface_groups, outputnode, [ + (f'out{i}', surf) for i, surf in enumerate(surfaces, start=1) + ]), + ]) # fmt:skip + return workflow + + +def init_gifti_morphometrics_wf( + *, + morphometrics: ty.List[str] = ["thickness", "curv", "sulc"], + name: str = "gifti_morphometrics_wf", +): + r""" + Prepare GIFTI shape files from morphometrics found in a FreeSurfer subjects + directory. + + The default morphometrics are ``lh/rh.thickness``, ``lh/rh.curv``, and + ``lh/rh.sulc``. + + Workflow Graph + .. workflow:: + :graph2use: orig + :simple_form: yes + + from smriprep.workflows.surfaces import init_gifti_morphometrics_wf + wf = init_gifti_morphometrics_wf() + + Inputs + ------ + subjects_dir + FreeSurfer SUBJECTS_DIR + subject_id + FreeSurfer subject ID + fsnative2t1w_xfm + LTA formatted affine transform file (inverse) + + Outputs + ------- + morphometrics + GIFTI shape files for all requested morphometrics + ```` + Left and right GIFTIs for each morphometry type passed to ``morphometrics`` + + """ + from ..interfaces.freesurfer import MRIsConvertData + + workflow = Workflow(name=name) + + inputnode = pe.Node( + niu.IdentityInterface(["subjects_dir", "subject_id"]), + name="inputnode", + ) + outputnode = pe.Node( + niu.IdentityInterface( + [ + "morphometrics", + *morphometrics, + ] + ), + name="outputnode", + ) + + get_subject = pe.Node(nio.FreeSurferSource(), name="get_surfaces") + + morphometry_list = pe.Node( + niu.Merge(len(morphometrics), ravel_inputs=True), name="surfmorph_list", run_without_submitting=True, ) @@ -784,47 +984,175 @@ def init_gifti_surface_wf(*, name="gifti_surface_wf"): name="morphs2gii", ) - agg_surfaces = pe.Node(AggregateSurfaces(), name="agg_surfaces") + morph_groups = pe.Node( + niu.Split(splits=[2] * len(morphometrics)), + name="morph_groups", + run_without_submitting=True, + ) # fmt:off workflow.connect([ - (inputnode, get_surfaces, [('subjects_dir', 'subjects_dir'), - ('subject_id', 'subject_id')]), - (inputnode, save_midthickness, [('subjects_dir', 'base_directory'), - ('subject_id', 'container')]), - # Generate midthickness surfaces and save to FreeSurfer derivatives - (get_surfaces, midthickness, [('white', 'in_file'), - ('graymid', 'graymid')]), - (midthickness, save_midthickness, [('out_file', 'surf.@graymid')]), - # Produce valid GIFTI surface files (dense mesh) - (get_surfaces, surface_list, [('white', 'in1'), - ('pial', 'in2'), - ('inflated', 'in3')]), - (save_midthickness, surface_list, [('out_file', 'in4')]), - (surface_list, fs2gii, [('out', 'in_file')]), - (fs2gii, fix_surfs, [('converted', 'in_file')]), - (inputnode, fix_surfs, [('fsnative2t1w_xfm', 'transform_file')]), - (fix_surfs, outputnode, [('out_file', 'surfaces')]), - (get_surfaces, surfmorph_list, [('thickness', 'in1'), - ('sulc', 'in2'), - ('curv', 'in3')]), - (surfmorph_list, morphs2gii, [('out', 'scalarcurv_file')]), + (inputnode, get_subject, [('subjects_dir', 'subjects_dir'), + ('subject_id', 'subject_id')]), + (get_subject, morphometry_list, [ + ((morph, _sorted_by_basename), f'in{i}') + for i, morph in enumerate(morphometrics, start=1) + ]), + (morphometry_list, morphs2gii, [('out', 'scalarcurv_file')]), (morphs2gii, outputnode, [('converted', 'morphometrics')]), # Output individual surfaces as well - (fix_surfs, agg_surfaces, [('out_file', 'surfaces')]), - (morphs2gii, agg_surfaces, [('converted', 'morphometrics')]), - (agg_surfaces, outputnode, [('pial', 'pial'), - ('white', 'white'), - ('inflated', 'inflated'), - ('midthickness', 'midthickness'), - ('thickness', 'thickness'), - ('sulc', 'sulc'), - ('curv', 'curv')]), + (morphs2gii, morph_groups, [('converted', 'inlist')]), + (morph_groups, outputnode, [ + (f'out{i}', surf) for i, surf in enumerate(morphometrics, start=1) + ]), ]) # fmt:on return workflow +def init_hcp_morphometrics_wf( + *, + omp_nthreads: int, + name: str = "hcp_morphometrics_wf", +) -> Workflow: + """Convert FreeSurfer-style morphometrics to HCP-style. + + The HCP uses different conventions for curvature and sulcal depth from FreeSurfer, + and performs some geometric preprocessing steps to get final metrics and a + subject-specific cortical ROI. + + Workflow Graph + + .. workflow:: + + from smriprep.workflows.surfaces import init_hcp_morphometrics_wf + wf = init_hcp_morphometrics_wf(omp_nthreads=1) + + Parameters + ---------- + omp_nthreads + Maximum number of threads an individual process may use + + Inputs + ------ + subject_id + FreeSurfer subject ID + thickness + FreeSurfer thickness file in GIFTI format + curv + FreeSurfer curvature file in GIFTI format + sulc + FreeSurfer sulcal depth file in GIFTI format + + Outputs + ------- + thickness + HCP-style thickness file in GIFTI format + curv + HCP-style curvature file in GIFTI format + sulc + HCP-style sulcal depth file in GIFTI format + roi + HCP-style cortical ROI file in GIFTI format + """ + DEFAULT_MEMORY_MIN_GB = 0.01 + + workflow = Workflow(name=name) + + inputnode = pe.Node( + niu.IdentityInterface(fields=['subject_id', 'thickness', 'curv', 'sulc', 'midthickness']), + name='inputnode', + ) + + itersource = pe.Node( + niu.IdentityInterface(fields=['hemi']), + name='itersource', + iterables=[('hemi', ['L', 'R'])], + ) + + outputnode = pe.JoinNode( + niu.IdentityInterface(fields=["thickness", "curv", "sulc", "roi"]), + name="outputnode", + joinsource="itersource", + ) + + select_surfaces = pe.Node( + KeySelect( + fields=[ + "thickness", + "curv", + "sulc", + "midthickness", + ], + keys=["L", "R"], + ), + name="select_surfaces", + run_without_submitting=True, + ) + + # HCP inverts curvature and sulcal depth relative to FreeSurfer + invert_curv = pe.Node(MetricMath(metric='curv', operation='invert'), name='invert_curv') + invert_sulc = pe.Node(MetricMath(metric='sulc', operation='invert'), name='invert_sulc') + # Thickness is presumably already positive, but HCP uses abs(-thickness) + abs_thickness = pe.Node(MetricMath(metric='thickness', operation='abs'), name='abs_thickness') + + # Native ROI is thickness > 0, with holes and islands filled + initial_roi = pe.Node(MetricMath(metric='roi', operation='bin'), name='initial_roi') + fill_holes = pe.Node(MetricFillHoles(), name="fill_holes", mem_gb=DEFAULT_MEMORY_MIN_GB) + native_roi = pe.Node(MetricRemoveIslands(), name="native_roi", mem_gb=DEFAULT_MEMORY_MIN_GB) + + # Dilation happens separately from ROI creation + dilate_curv = pe.Node( + MetricDilate(distance=10, nearest=True), + name="dilate_curv", + n_procs=omp_nthreads, + mem_gb=DEFAULT_MEMORY_MIN_GB, + ) + dilate_thickness = pe.Node( + MetricDilate(distance=10, nearest=True), + name="dilate_thickness", + n_procs=omp_nthreads, + mem_gb=DEFAULT_MEMORY_MIN_GB, + ) + + workflow.connect([ + (inputnode, select_surfaces, [ + ('thickness', 'thickness'), + ('curv', 'curv'), + ('sulc', 'sulc'), + ('midthickness', 'midthickness'), + ]), + (inputnode, invert_curv, [('subject_id', 'subject_id')]), + (inputnode, invert_sulc, [('subject_id', 'subject_id')]), + (inputnode, abs_thickness, [('subject_id', 'subject_id')]), + (itersource, select_surfaces, [('hemi', 'key')]), + (itersource, invert_curv, [('hemi', 'hemisphere')]), + (itersource, invert_sulc, [('hemi', 'hemisphere')]), + (itersource, abs_thickness, [('hemi', 'hemisphere')]), + (select_surfaces, invert_curv, [('curv', 'metric_file')]), + (select_surfaces, invert_sulc, [('sulc', 'metric_file')]), + (select_surfaces, abs_thickness, [('thickness', 'metric_file')]), + (select_surfaces, dilate_curv, [('midthickness', 'surf_file')]), + (select_surfaces, dilate_thickness, [('midthickness', 'surf_file')]), + (invert_curv, dilate_curv, [('metric_file', 'in_file')]), + (abs_thickness, dilate_thickness, [('metric_file', 'in_file')]), + (dilate_curv, outputnode, [('out_file', 'curv')]), + (dilate_thickness, outputnode, [('out_file', 'thickness')]), + (invert_sulc, outputnode, [('metric_file', 'sulc')]), + # Native ROI file from thickness + (inputnode, initial_roi, [('subject_id', 'subject_id')]), + (itersource, initial_roi, [('hemi', 'hemisphere')]), + (abs_thickness, initial_roi, [('metric_file', 'metric_file')]), + (select_surfaces, fill_holes, [('midthickness', 'surface_file')]), + (select_surfaces, native_roi, [('midthickness', 'surface_file')]), + (initial_roi, fill_holes, [('metric_file', 'metric_file')]), + (fill_holes, native_roi, [('out_file', 'metric_file')]), + (native_roi, outputnode, [('out_file', 'roi')]), + ]) # fmt:skip + + return workflow + + def init_segs_to_native_wf(*, name="segs_to_native", segmentation="aseg"): """ Get a segmentation from FreeSurfer conformed space into native T1w space. @@ -867,6 +1195,9 @@ 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") + + lta = pe.Node(ConcatenateXFMs(out_fmt="fs"), name="lta", run_without_submitting=True) + # 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( @@ -897,9 +1228,12 @@ def _sel(x): (inputnode, fssource, [ ('subjects_dir', 'subjects_dir'), ('subject_id', 'subject_id')]), - (inputnode, resample, [('in_file', 'target_file'), - ('fsnative2t1w_xfm', 'lta_file')]), + (inputnode, lta, [('in_file', 'reference'), + ('fsnative2t1w_xfm', 'in_xfms')]), + (fssource, lta, [('T1', 'moving')]), + (inputnode, resample, [('in_file', 'target_file')]), (fssource, resample, [(segmentation, 'source_file')]), + (lta, resample, [('out_xfm', 'lta_file')]), (resample, outputnode, [('transformed_file', 'out_file')]), ]) # fmt:on @@ -907,51 +1241,39 @@ def _sel(x): def init_anat_ribbon_wf(name="anat_ribbon_wf"): + """Create anatomical ribbon mask + + Workflow Graph + + .. workflow:: + :graph2use: orig + :simple_form: yes + + from smriprep.workflows.surfaces import init_anat_ribbon_wf + wf = init_anat_ribbon_wf() + + Inputs + ------ + white + Left and right gray/white surfaces (as GIFTI files) + pial + Left and right pial surfaces (as GIFTI files) + ref_file + Reference image (one 3D volume) to define the target space + + Outputs + ------- + anat_ribbon + Cortical gray matter mask, sampled into ``ref_file`` space + """ DEFAULT_MEMORY_MIN_GB = 0.01 workflow = pe.Workflow(name=name) inputnode = pe.Node( - niu.IdentityInterface( - fields=[ - "surfaces", - "t1w_mask", - ] - ), + niu.IdentityInterface(fields=["white", "pial", "ref_file"]), name="inputnode", ) - outputnode = pe.Node( - niu.IdentityInterface( - fields=[ - "anat_ribbon", - ] - ), - name="outputnode", - ) - - # 0, 1 = wm; 2, 3 = pial; 6, 7 = mid - # note that order of lh / rh within each surf type is not guaranteed due to use - # of unsorted glob by FreeSurferSource prior, but we can do a sort - # to ensure consistent ordering - select_wm = pe.Node( - niu.Select(index=[0, 1]), - name="select_wm", - mem_gb=DEFAULT_MEMORY_MIN_GB, - run_without_submitting=True, - ) - - select_pial = pe.Node( - niu.Select(index=[2, 3]), - name="select_pial", - mem_gb=DEFAULT_MEMORY_MIN_GB, - run_without_submitting=True, - ) - - select_midthick = pe.Node( - niu.Select(index=[6, 7]), - name="select_midthick", - mem_gb=DEFAULT_MEMORY_MIN_GB, - run_without_submitting=True, - ) + outputnode = pe.Node(niu.IdentityInterface(fields=["anat_ribbon"]), name="outputnode") create_wm_distvol = pe.MapNode( CreateSignedDistanceVolume(), @@ -965,107 +1287,109 @@ def init_anat_ribbon_wf(name="anat_ribbon_wf"): name="create_pial_distvol", ) - thresh_wm_distvol = pe.MapNode( - fsl.maths.MathsCommand(args="-thr 0 -bin -mul 255"), - iterfield=["in_file"], - name="thresh_wm_distvol", - mem_gb=DEFAULT_MEMORY_MIN_GB, - ) + make_ribbon = pe.Node(MakeRibbon(), name="make_ribbon", mem_gb=DEFAULT_MEMORY_MIN_GB) - uthresh_pial_distvol = pe.MapNode( - fsl.maths.MathsCommand(args="-uthr 0 -abs -bin -mul 255"), - iterfield=["in_file"], - name="uthresh_pial_distvol", - mem_gb=DEFAULT_MEMORY_MIN_GB, + # fmt: off + workflow.connect( + [ + (inputnode, create_wm_distvol, [ + ("white", "surf_file"), + ("ref_file", "ref_file"), + ]), + (inputnode, create_pial_distvol, [ + ("pial", "surf_file"), + ("ref_file", "ref_file"), + ]), + (create_wm_distvol, make_ribbon, [("out_file", "white_distvols")]), + (create_pial_distvol, make_ribbon, [("out_file", "pial_distvols")]), + (make_ribbon, outputnode, [("ribbon", "anat_ribbon")]), + ] ) + # fmt: on + return workflow - bin_wm_distvol = pe.MapNode( - fsl.maths.UnaryMaths(operation="bin"), - iterfield=["in_file"], - name="bin_wm_distvol", - mem_gb=DEFAULT_MEMORY_MIN_GB, - ) - bin_pial_distvol = pe.MapNode( - fsl.maths.UnaryMaths(operation="bin"), - iterfield=["in_file"], - name="bin_pial_distvol", - mem_gb=DEFAULT_MEMORY_MIN_GB, - ) +def init_resample_midthickness_wf( + grayord_density: ty.Literal['91k', '170k'], + name: str = "resample_midthickness_wf", +): + """ + Resample subject midthickness surface to specified density. - split_wm_distvol = pe.Node( - niu.Split(splits=[1, 1]), - name="split_wm_distvol", - mem_gb=DEFAULT_MEMORY_MIN_GB, - run_without_submitting=True, - ) + Workflow Graph + .. workflow:: + :graph2use: colored + :simple_form: yes - merge_wm_distvol_no_flatten = pe.Node( - niu.Merge(2), - no_flatten=True, - name="merge_wm_distvol_no_flatten", - mem_gb=DEFAULT_MEMORY_MIN_GB, - run_without_submitting=True, - ) + from smriprep.workflows.surfaces import init_resample_midthickness_wf + wf = init_resample_midthickness_wf(grayord_density="91k") - make_ribbon_vol = pe.MapNode( - fsl.maths.MultiImageMaths(op_string="-mas %s -mul 255"), - iterfield=["in_file", "operand_files"], - name="make_ribbon_vol", - mem_gb=DEFAULT_MEMORY_MIN_GB, - ) + Parameters + ---------- + grayord_density : :obj:`str` + Either `91k` or `170k`, representing the total of vertices or *grayordinates*. + name : :obj:`str` + Unique name for the subworkflow (default: ``"resample_midthickness_wf"``) - bin_ribbon_vol = pe.MapNode( - fsl.maths.UnaryMaths(operation="bin"), - iterfield=["in_file"], - name="bin_ribbon_vol", - mem_gb=DEFAULT_MEMORY_MIN_GB, - ) + Inputs + ------ + midthickness + GIFTI surface mesh corresponding to the midthickness surface + sphere_reg_fsLR + GIFTI surface mesh corresponding to the subject's fsLR registration sphere - split_squeeze_ribbon_vol = pe.Node( - niu.Split(splits=[1, 1], squeeze=True), - name="split_squeeze_ribbon", - mem_gb=DEFAULT_MEMORY_MIN_GB, - run_without_submitting=True, - ) + Outputs + ------- + midthickness + GIFTI surface mesh corresponding to the midthickness surface, resampled to fsLR + """ + import templateflow.api as tf + from niworkflows.engine.workflows import LiterateWorkflow as Workflow - combine_ribbon_vol_hemis = pe.Node( - fsl.maths.BinaryMaths(operation="add"), - name="combine_ribbon_vol_hemis", - mem_gb=DEFAULT_MEMORY_MIN_GB, + workflow = Workflow(name=name) + + fslr_density = "32k" if grayord_density == "91k" else "59k" + + inputnode = pe.Node( + niu.IdentityInterface(fields=["midthickness", "sphere_reg_fsLR"]), + name="inputnode", ) - # make HCP-style ribbon volume in T1w space - workflow.connect( - [ - (inputnode, select_wm, [("surfaces", "inlist")]), - (inputnode, select_pial, [("surfaces", "inlist")]), - (inputnode, select_midthick, [("surfaces", "inlist")]), - (select_wm, create_wm_distvol, [(("out", _sorted_by_basename), "surf_file")]), - (inputnode, create_wm_distvol, [("t1w_mask", "ref_file")]), - (select_pial, create_pial_distvol, [(("out", _sorted_by_basename), "surf_file")]), - (inputnode, create_pial_distvol, [("t1w_mask", "ref_file")]), - (create_wm_distvol, thresh_wm_distvol, [("out_file", "in_file")]), - (create_pial_distvol, uthresh_pial_distvol, [("out_file", "in_file")]), - (thresh_wm_distvol, bin_wm_distvol, [("out_file", "in_file")]), - (uthresh_pial_distvol, bin_pial_distvol, [("out_file", "in_file")]), - (bin_wm_distvol, split_wm_distvol, [("out_file", "inlist")]), - (split_wm_distvol, merge_wm_distvol_no_flatten, [("out1", "in1")]), - (split_wm_distvol, merge_wm_distvol_no_flatten, [("out2", "in2")]), - (bin_pial_distvol, make_ribbon_vol, [("out_file", "in_file")]), - (merge_wm_distvol_no_flatten, make_ribbon_vol, [("out", "operand_files")]), - (make_ribbon_vol, bin_ribbon_vol, [("out_file", "in_file")]), - (bin_ribbon_vol, split_squeeze_ribbon_vol, [("out_file", "inlist")]), - (split_squeeze_ribbon_vol, combine_ribbon_vol_hemis, [("out1", "in_file")]), - (split_squeeze_ribbon_vol, combine_ribbon_vol_hemis, [("out2", "operand_file")]), - (combine_ribbon_vol_hemis, outputnode, [("out_file", "anat_ribbon")]), - ] + outputnode = pe.Node(niu.IdentityInterface(fields=["midthickness_fsLR"]), name="outputnode") + + resampler = pe.MapNode( + SurfaceResample(method='BARYCENTRIC'), + iterfield=['surface_in', 'current_sphere', 'new_sphere'], + name="resampler", ) + resampler.inputs.new_sphere = [ + str( + tf.get( + template='fsLR', + density=fslr_density, + suffix='sphere', + hemi=hemi, + space=None, + extension='.surf.gii', + ) + ) + for hemi in ['L', 'R'] + ] + + workflow.connect([ + (inputnode, resampler, [ + ('midthickness', 'surface_in'), + ('sphere_reg_fsLR', 'current_sphere'), + ]), + (resampler, outputnode, [('surface_out', 'midthickness_fsLR')]), + ]) # fmt:skip + return workflow def init_morph_grayords_wf( grayord_density: ty.Literal['91k', '170k'], + omp_nthreads: int, name: str = "morph_grayords_wf", ): """ @@ -1079,7 +1403,7 @@ def init_morph_grayords_wf( :simple_form: yes from smriprep.workflows.surfaces import init_morph_grayords_wf - wf = init_morph_grayords_wf(grayord_density="91k") + wf = init_morph_grayords_wf(grayord_density="91k", omp_nthreads=1) Parameters ---------- @@ -1090,163 +1414,168 @@ def init_morph_grayords_wf( Inputs ------ - subject_id : :obj:`str` - FreeSurfer subject ID - subjects_dir : :obj:`str` - FreeSurfer SUBJECTS_DIR + curv + HCP-style curvature file in GIFTI format + sulc + HCP-style sulcal depth file in GIFTI format + thickness + HCP-style thickness file in GIFTI format + roi + HCP-style cortical ROI file in GIFTI format + midthickness + GIFTI surface mesh corresponding to the midthickness surface + sphere_reg_fsLR + GIFTI surface mesh corresponding to the subject's fsLR registration sphere Outputs ------- - cifti_morph : :obj:`list` of :obj:`str` - Paths of CIFTI dscalar files - cifti_metadata : :obj:`list` of :obj:`str` - Paths to JSON files containing metadata corresponding to ``cifti_morph`` - + curv_fsLR + HCP-style curvature file in CIFTI-2 format, resampled to fsLR + curv_metadata + Path to JSON file containing metadata corresponding to ``curv_fsLR`` + sulc_fsLR + HCP-style sulcal depth file in CIFTI-2 format, resampled to fsLR + sulc_metadata + Path to JSON file containing metadata corresponding to ``sulc_fsLR`` + thickness_fsLR + HCP-style thickness file in CIFTI-2 format, resampled to fsLR """ + 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"""\ -*Grayordinate* "dscalar" files [@hcppipelines] containing {grayord_density} samples were -also generated using the highest-resolution ``fsaverage`` as an intermediate standardized -surface space. +*Grayordinate* "dscalar" files containing {grayord_density} samples were +resampled onto fsLR using the Connectome Workbench [@hcppipelines]. """ fslr_density = "32k" if grayord_density == "91k" else "59k" inputnode = pe.Node( - niu.IdentityInterface(fields=["subject_id", "subjects_dir"]), + niu.IdentityInterface( + fields=[ + "curv", + "sulc", + "thickness", + "roi", + "midthickness", + "midthickness_fsLR", + "sphere_reg_fsLR", + ] + ), name="inputnode", ) + # Note we have to use a different name than itersource to + # avoid duplicates in the workflow graph, even though the + # previous was already Joined + hemisource = pe.Node( + niu.IdentityInterface(fields=['hemi']), + name='hemisource', + iterables=[('hemi', ['L', 'R'])], + ) + outputnode = pe.Node( - niu.IdentityInterface(fields=["cifti_morph", "cifti_metadata"]), + niu.IdentityInterface( + fields=[ + "curv_fsLR", + "curv_metadata", + "sulc_fsLR", + "sulc_metadata", + "thickness_fsLR", + "thickness_metadata", + ] + ), name="outputnode", ) - get_surfaces = pe.Node(nio.FreeSurferSource(), name="get_surfaces") - - surfmorph_list = pe.Node( - niu.Merge(3, ravel_inputs=True), - name="surfmorph_list", + atlases = load_resource('atlases') + select_surfaces = pe.Node( + KeySelect( + fields=[ + "curv", + "sulc", + "thickness", + "roi", + "midthickness", + "midthickness_fsLR", + "sphere_reg_fsLR", + "template_sphere", + "template_roi", + ], + keys=["L", "R"], + ), + name="select_surfaces", run_without_submitting=True, ) - - surf2surf = pe.MapNode( - fs.SurfaceTransform(target_subject="fsaverage", target_type="gii"), - iterfield=["source_file", "hemi"], - name="surf2surf", - mem_gb=0.01, - ) - surf2surf.inputs.hemi = ["lh", "rh"] * 3 - - # Setup Workbench command. LR ordering for hemi can be assumed, as it is imposed - # by the iterfield of the MapNode in the surface sampling workflow above. - resample = pe.MapNode( - wb.MetricResample(method="ADAP_BARY_AREA", area_metrics=True), - name="resample", - iterfield=[ - "in_file", - "out_file", - "new_sphere", - "new_area", - "current_sphere", - "current_area", - ], - ) - resample.inputs.current_sphere = [ - str( - tf.get( - "fsaverage", - hemi=hemi, - density="164k", - desc="std", - suffix="sphere", - extension=".surf.gii", - ) - ) - for hemi in "LR" - ] * 3 - resample.inputs.current_area = [ - str( - tf.get( - "fsaverage", - hemi=hemi, - density="164k", - desc="vaavg", - suffix="midthickness", - extension=".shape.gii", - ) - ) - for hemi in "LR" - ] * 3 - resample.inputs.new_sphere = [ + select_surfaces.inputs.template_sphere = [ str( tf.get( - "fsLR", - space="fsaverage", - hemi=hemi, + template='fsLR', density=fslr_density, - suffix="sphere", - extension=".surf.gii", - ) - ) - for hemi in "LR" - ] * 3 - resample.inputs.new_area = [ - str( - tf.get( - "fsLR", + suffix='sphere', hemi=hemi, - density=fslr_density, - desc="vaavg", - suffix="midthickness", - extension=".shape.gii", + space=None, + extension='.surf.gii', ) ) - for hemi in "LR" - ] * 3 - resample.inputs.out_file = [ - f"space-fsLR_hemi-{h}_den-{grayord_density}_{morph}.shape.gii" - # Order: curv-L, curv-R, sulc-L, sulc-R, thickness-L, thickness-R - for morph in ('curv', 'sulc', 'thickness') - for h in "LR" + for hemi in ['L', 'R'] + ] + select_surfaces.inputs.template_roi = [ + str(atlases / f'L.atlasroi.{fslr_density}_fs_LR.shape.gii'), + str(atlases / f'R.atlasroi.{fslr_density}_fs_LR.shape.gii'), ] - gen_cifti = pe.MapNode( - GenerateDScalar( - grayordinates=grayord_density, - ), - iterfield=['scalar_name', 'scalar_surfs'], - name="gen_cifti", - ) - gen_cifti.inputs.scalar_name = ['curv', 'sulc', 'thickness'] - - # fmt: off workflow.connect([ - (inputnode, get_surfaces, [ - ('subject_id', 'subject_id'), - ('subjects_dir', 'subjects_dir'), - ]), - (inputnode, surf2surf, [ - ('subject_id', 'source_subject'), - ('subjects_dir', 'subjects_dir'), + (inputnode, select_surfaces, [ + ('curv', 'curv'), + ('sulc', 'sulc'), + ('thickness', 'thickness'), + ('roi', 'roi'), + ('midthickness', 'midthickness'), + ('midthickness_fsLR', 'midthickness_fsLR'), + ('sphere_reg_fsLR', 'sphere_reg_fsLR'), ]), - (get_surfaces, surfmorph_list, [ - (('curv', _sorted_by_basename), 'in1'), - (('sulc', _sorted_by_basename), 'in2'), - (('thickness', _sorted_by_basename), 'in3'), - ]), - (surfmorph_list, surf2surf, [('out', 'source_file')]), - (surf2surf, resample, [('out_file', 'in_file')]), - (resample, gen_cifti, [ - (("out_file", _collate), "scalar_surfs")]), - (gen_cifti, outputnode, [("out_file", "cifti_morph"), - ("out_metadata", "cifti_metadata")]), - ]) - # fmt: on + (hemisource, select_surfaces, [('hemi', 'key')]), + ]) # fmt:skip + + for metric in ('curv', 'sulc', 'thickness'): + resampler = pe.Node( + MetricResample(method='ADAP_BARY_AREA', area_surfs=True), + name=f"resample_{metric}", + n_procs=omp_nthreads, + ) + mask_fsLR = pe.Node( + MetricMask(), + name=f"mask_{metric}", + n_procs=omp_nthreads, + ) + cifti_metric = pe.JoinNode( + GenerateDScalar(grayordinates=grayord_density, scalar_name=metric), + name=f"cifti_{metric}", + joinfield=['scalar_surfs'], + joinsource="hemisource", + ) + + workflow.connect([ + (select_surfaces, resampler, [ + (metric, 'in_file'), + ('sphere_reg_fsLR', 'current_sphere'), + ('template_sphere', 'new_sphere'), + ('midthickness', 'current_area'), + ('midthickness_fsLR', 'new_area'), + ('roi', 'roi_metric'), + ]), + (select_surfaces, mask_fsLR, [('template_roi', 'mask')]), + (resampler, mask_fsLR, [('out_file', 'in_file')]), + (mask_fsLR, cifti_metric, [('out_file', 'scalar_surfs')]), + (cifti_metric, outputnode, [ + ('out_file', f'{metric}_fsLR'), + ('out_metadata', f'{metric}_metadata'), + ]), + ]) # fmt:skip return workflow @@ -1273,3 +1602,60 @@ def _sorted_by_basename(inlist): def _collate(files): return [files[i : i + 2] for i in range(0, len(files), 2)] + + +def _extract_fs_fields(filenames: str | list[str]) -> tuple[str, str]: + from pathlib import Path + + if isinstance(filenames, str): + filenames = [filenames] + paths = [Path(fn) for fn in filenames] + sub_dir = paths[0].parent.parent + subjects_dir, subject_id = sub_dir.parent, sub_dir.name + assert all(path == subjects_dir / subject_id / 'surf' / path.name for path in paths) + return str(subjects_dir), subject_id + + +def _get_surfaces(subjects_dir: str, subject_id: str, surfaces: list[str]) -> tuple[list[str]]: + """ + Get a list of FreeSurfer surface files for a given subject. + + If ``midthickness`` is requested but not present in the directory, + ``graymid`` will be returned instead. For surfaces with dots (``.``) in + their names, pass the name with underscores (``_``). + + Parameters + ---------- + subjects_dir + FreeSurfer SUBJECTS_DIR + subject_id + FreeSurfer subject ID + surfaces + List of surfaces to fetch + + Returns + ------- + tuple + A list of surfaces for each requested surface, sorted + + """ + from pathlib import Path + + expanded_surfaces = surfaces.copy() + if "midthickness" in surfaces: + expanded_surfaces.append("graymid") + + surf_dir = Path(subjects_dir) / subject_id / "surf" + all_surfs = { + surface: sorted( + str(fn) + for fn in surf_dir.glob(f"[lr]h.{surface.replace('_', '.')}") + ) + for surface in expanded_surfaces + } + + if all_surfs.get("graymid") and not all_surfs.get("midthickness"): + all_surfs["midthickness"] = all_surfs.pop("graymid") + + ret = tuple(all_surfs[surface] for surface in surfaces) + return ret if len(ret) > 1 else ret[0] diff --git a/smriprep/workflows/tests/__init__.py b/smriprep/workflows/tests/__init__.py new file mode 100644 index 0000000000..e69de29bb2 diff --git a/smriprep/workflows/tests/test_surfaces.py b/smriprep/workflows/tests/test_surfaces.py new file mode 100644 index 0000000000..3589624c7e --- /dev/null +++ b/smriprep/workflows/tests/test_surfaces.py @@ -0,0 +1,56 @@ +import os +from pathlib import Path +from shutil import which + +import nibabel as nb +import numpy as np +import pytest +from nipype.pipeline import engine as pe + +from ..surfaces import init_anat_ribbon_wf, init_gifti_surfaces_wf +from ...data import load_resource + + +def test_ribbon_workflow(tmp_path: Path): + """Create ribbon mask for fsaverage subject""" + + for command in ("mris_convert", "wb_command"): + if not which(command): + pytest.skip(f"Could not find {command} in PATH") + + if not os.path.exists(os.getenv('SUBJECTS_DIR')): + pytest.skip("Could not find $SUBJECTS_DIR") + + # Low-res file that includes the fsaverage surfaces in its bounding box + # We will use it both as a template and a comparison. + test_ribbon = load_resource( + "../interfaces/tests/data/sub-fsaverage_res-4_desc-cropped_ribbon.nii.gz" + ) + + gifti_surfaces_wf = init_gifti_surfaces_wf(surfaces=['white', 'pial']) + anat_ribbon_wf = init_anat_ribbon_wf() + anat_ribbon_wf.inputs.inputnode.ref_file = test_ribbon + + gifti_surfaces_wf.inputs.inputnode.subjects_dir = os.getenv('SUBJECTS_DIR') + gifti_surfaces_wf.inputs.inputnode.subject_id = 'fsaverage' + + wf = pe.Workflow(name='test_ribbon_wf', base_dir=tmp_path) + # fmt: off + wf.connect([ + (gifti_surfaces_wf, anat_ribbon_wf, [ + ('outputnode.white', 'inputnode.white'), + ('outputnode.pial', 'inputnode.pial'), + ]), + ]) + # fmt: on + result = wf.run() + + make_ribbon = next(node for node in result.nodes() if node.name == 'make_ribbon') + + expected = nb.load(test_ribbon) + ribbon = nb.load(make_ribbon.result.outputs.ribbon) + + assert ribbon.shape == expected.shape + assert np.allclose(ribbon.affine, expected.affine) + # Mask data is binary, so we can use np.array_equal + assert np.array_equal(ribbon.dataobj, expected.dataobj) diff --git a/wrapper/LICENSE b/wrapper/LICENSE index c4140b1eb0..26d2ab656b 100644 --- a/wrapper/LICENSE +++ b/wrapper/LICENSE @@ -24,4 +24,4 @@ DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE -OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. \ No newline at end of file +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. diff --git a/wrapper/smriprep_docker.py b/wrapper/smriprep_docker.py index bbb71c280a..7775316e04 100755 --- a/wrapper/smriprep_docker.py +++ b/wrapper/smriprep_docker.py @@ -165,7 +165,7 @@ def _get_posargs(usage): # Make sure we're not clobbering options we don't mean to overlap = set(w_flags).intersection(t_flags) expected_overlap = { - "h", "version", "w", "fs-license-file", "use-plugin", "fs-subjects-dir" + "h", "version", "w", "d", "fs-license-file", "use-plugin", "fs-subjects-dir" } assert overlap == expected_overlap, "Clobbering options: {}".format( @@ -280,6 +280,14 @@ def _is_file(path, parser): type=os.path.abspath, help="path where intermediate results should be stored", ) + g_wrap.add_argument( + "-d", + "--derivatives", + nargs="+", + metavar="PATH", + action=ToDict, + help="Search PATH(s) for pre-computed derivatives.", + ) g_wrap.add_argument( "--fs-license-file", metavar="PATH", @@ -453,6 +461,13 @@ def main(): command.extend(['-v', '{}:/opt/subjects'.format(opts.fs_subjects_dir)]) unknown_args.extend(['--fs-subjects-dir', '/opt/subjects']) + # Patch derivatives for searching + if opts.derivatives: + unknown_args.append("--derivatives") + for deriv, deriv_path in opts.derivatives.items(): + command.extend(['-v', '{}:/deriv/{}:ro'.format(deriv_path, deriv)]) + unknown_args.append("/deriv/{}".format(deriv)) + if opts.config: command.extend(["-v", ":".join((opts.config, "/home/smriprep/.nipype/nipype.cfg", "ro"))])