Skip to content

Commit

Permalink
feat(TOPUP): an initial implementation of SD estimation.
Browse files Browse the repository at this point in the history
Adds a new subworkflow based on FSL TOPUP to integrate SD estimation
for the ds001771 dataset.

- [x] Pin niworkflows to current master (while I release 1.2.0rc5
  containing nipreps/niworkflows#503, nipreps/niworkflows#504, which
  are used here).
- [x] Create a new sdc estimation workflow, with the expectation of
  upstreaming it to SDCFlows.
- [x] Implement the barebones of how nipreps/sdcflows#101 could look
  like. Also to be upstreamed to SDCFlows when mature.
- [x] Stick TOPUP from FSL 6.0.3 in the Docker image, since topup from
  FSL 5.0.x is really unstable (for instance, it fails with a
  segmentation fault on the workflow of ds001771)

Resolves: nipreps#92
  • Loading branch information
oesteban committed Jun 2, 2020
1 parent 83f6ed1 commit 61a9125
Show file tree
Hide file tree
Showing 13 changed files with 274 additions and 10 deletions.
Binary file added .docker/fsl-6.0/bin/topup
Binary file not shown.
Binary file added .docker/fsl-6.0/lib/libgfortran.so.3
Binary file not shown.
Binary file added .docker/fsl-6.0/lib/libopenblas.so.0
Binary file not shown.
3 changes: 3 additions & 0 deletions Dockerfile
Original file line number Diff line number Diff line change
Expand Up @@ -92,6 +92,9 @@ ENV FSLDIR="/usr/share/fsl/5.0" \
AFNI_PLUGINPATH="/usr/lib/afni/plugins"
ENV PATH="/usr/lib/fsl/5.0:/usr/lib/afni/bin:$PATH"

COPY .docker/fsl-6.0/bin/topup /usr/share/fsl/5.0/bin/topup
COPY .docker/fsl-6.0/lib/* /usr/lib/fsl/5.0/

# Installing ANTs 2.2.0 (NeuroDocker build)
ENV ANTSPATH=/usr/lib/ants
RUN mkdir -p $ANTSPATH && \
Expand Down
13 changes: 13 additions & 0 deletions dmriprep/config/reports-spec.yml
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,19 @@ sections:
caption: Surfaces (white and pial) reconstructed with FreeSurfer (<code>recon-all</code>)
overlaid on the participant's T1w template.
subtitle: Surface reconstruction
- name: Fieldmaps
ordering: session,acquisition,run
reportlets:
- bids: {datatype: dwi, desc: pepolar, suffix: dwi}
caption: Inhomogeneities of the *B0* field introduce (oftentimes severe) spatial distortions
along the phase-encoding direction of the image. Utilizing two or more images with different
phase-encoding polarities (PEPolar) or directions, it is possible to estimate the inhomogeneity
of the field. The plot below shows a reference EPI (echo-planar imaging) volume generated
using two or more EPI images with varying phase-encoding blips.
description: Hover on the panels with the mouse pointer to also visualize the intensity of the
inhomogeneity of the field in Hertz.
static: false
subtitle: "Susceptibility-derived Distortion Correction (SDC): field inhomogeneity estimation"
- name: Diffusion
ordering: session,acquisition,run
reportlets:
Expand Down
26 changes: 26 additions & 0 deletions dmriprep/data/flirtsch/b02b0.cnf
Original file line number Diff line number Diff line change
@@ -0,0 +1,26 @@
# Resolution (knot-spacing) of warps in mm
--warpres=20,16,14,12,10,6,4,4,4
# Subsampling level (a value of 2 indicates that a 2x2x2 neighbourhood is collapsed to 1 voxel)
--subsamp=2,2,2,2,2,1,1,1,1
# FWHM of gaussian smoothing
--fwhm=8,6,4,3,3,2,1,0,0
# Maximum number of iterations
--miter=5,5,5,5,5,10,10,20,20
# Relative weight of regularisation
--lambda=0.005,0.001,0.0001,0.000015,0.000005,0.0000005,0.00000005,0.0000000005,0.00000000001
# If set to 1 lambda is multiplied by the current average squared difference
--ssqlambda=1
# Regularisation model
--regmod=bending_energy
# If set to 1 movements are estimated along with the field
--estmov=1,1,1,1,1,0,0,0,0
# 0=Levenberg-Marquardt, 1=Scaled Conjugate Gradient
--minmet=0,0,0,0,0,1,1,1,1
# Quadratic or cubic splines
--splineorder=3
# Precision for calculation and storage of Hessian
--numprec=double
# Linear or spline interpolation
--interp=spline
# If set to 1 the images are individually scaled to a common mean intensity
--scale=1
26 changes: 26 additions & 0 deletions dmriprep/data/flirtsch/b02b0_1.cnf
Original file line number Diff line number Diff line change
@@ -0,0 +1,26 @@
# Resolution (knot-spacing) of warps in mm
--warpres=20,16,14,12,10,6,4,4,4
# Subsampling level (a value of 2 indicates that a 2x2x2 neighbourhood is collapsed to 1 voxel)
--subsamp=1,1,1,1,1,1,1,1,1
# FWHM of gaussian smoothing
--fwhm=8,6,4,3,3,2,1,0,0
# Maximum number of iterations
--miter=5,5,5,5,5,10,10,20,20
# Relative weight of regularisation
--lambda=0.0005,0.0001,0.00001,0.0000015,0.0000005,0.0000005,0.00000005,0.0000000005,0.00000000001
# If set to 1 lambda is multiplied by the current average squared difference
--ssqlambda=1
# Regularisation model
--regmod=bending_energy
# If set to 1 movements are estimated along with the field
--estmov=1,1,1,1,1,0,0,0,0
# 0=Levenberg-Marquardt, 1=Scaled Conjugate Gradient
--minmet=0,0,0,0,0,1,1,1,1
# Quadratic or cubic splines
--splineorder=3
# Precision for calculation and storage of Hessian
--numprec=double
# Linear or spline interpolation
--interp=spline
# If set to 1 the images are individually scaled to a common mean intensity
--scale=1
26 changes: 26 additions & 0 deletions dmriprep/data/flirtsch/b02b0_2.cnf
Original file line number Diff line number Diff line change
@@ -0,0 +1,26 @@
# Resolution (knot-spacing) of warps in mm
--warpres=20,16,14,12,10,6,4,4,4
# Subsampling level (a value of 2 indicates that a 2x2x2 neighbourhood is collapsed to 1 voxel)
--subsamp=2,2,2,2,2,1,1,1,1
# FWHM of gaussian smoothing
--fwhm=8,6,4,3,3,2,1,0,0
# Maximum number of iterations
--miter=5,5,5,5,5,10,10,20,20
# Relative weight of regularisation
--lambda=0.005,0.001,0.0001,0.000015,0.000005,0.0000005,0.00000005,0.0000000005,0.00000000001
# If set to 1 lambda is multiplied by the current average squared difference
--ssqlambda=1
# Regularisation model
--regmod=bending_energy
# If set to 1 movements are estimated along with the field
--estmov=1,1,1,1,1,0,0,0,0
# 0=Levenberg-Marquardt, 1=Scaled Conjugate Gradient
--minmet=0,0,0,0,0,1,1,1,1
# Quadratic or cubic splines
--splineorder=3
# Precision for calculation and storage of Hessian
--numprec=double
# Linear or spline interpolation
--interp=spline
# If set to 1 the images are individually scaled to a common mean intensity
--scale=1
26 changes: 26 additions & 0 deletions dmriprep/data/flirtsch/b02b0_4.cnf
Original file line number Diff line number Diff line change
@@ -0,0 +1,26 @@
# Resolution (knot-spacing) of warps in mm
--warpres=20,16,14,12,10,6,4,4,4
# Subsampling level (a value of 2 indicates that a 2x2x2 neighbourhood is collapsed to 1 voxel)
--subsamp=4,4,2,2,2,1,1,1,1
# FWHM of gaussian smoothing
--fwhm=8,6,4,3,3,2,1,0,0
# Maximum number of iterations
--miter=5,5,5,5,5,10,10,20,20
# Relative weight of regularisation
--lambda=0.035,0.006,0.0001,0.000015,0.000005,0.0000005,0.00000005,0.0000000005,0.00000000001
# If set to 1 lambda is multiplied by the current average squared difference
--ssqlambda=1
# Regularisation model
--regmod=bending_energy
# If set to 1 movements are estimated along with the field
--estmov=1,1,1,1,1,0,0,0,0
# 0=Levenberg-Marquardt, 1=Scaled Conjugate Gradient
--minmet=0,0,0,0,0,1,1,1,1
# Quadratic or cubic splines
--splineorder=3
# Precision for calculation and storage of Hessian
--numprec=double
# Linear or spline interpolation
--interp=spline
# If set to 1 the images are individually scaled to a common mean intensity
--scale=1
7 changes: 7 additions & 0 deletions dmriprep/workflows/base.py
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,7 @@
from ..interfaces.reports import SubjectSummary, AboutSummary
from ..utils.bids import collect_data
from .dwi.base import init_early_b0ref_wf
from .fmap.base import init_fmap_estimation_wf


def init_dmriprep_wf():
Expand Down Expand Up @@ -292,6 +293,12 @@ def init_single_subject_wf(subject_id):
]),
])

fmap_estimation_wf = init_fmap_estimation_wf(subject_data["dwi"])
workflow.connect([
(referencenode, fmap_estimation_wf, [
("dwi_reference", "inputnode.dwi_reference"),
("dwi_mask", "inputnode.dwi_mask")]),
])
return workflow


Expand Down
Empty file.
136 changes: 136 additions & 0 deletions dmriprep/workflows/fmap/base.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,136 @@
"""Deploy a susceptibility-distortion estimation strategy."""
from ... import config
from pkg_resources import resource_filename as _pkg_fname
from nipype.pipeline import engine as pe
from nipype.interfaces import utility as niu

from niworkflows.engine.workflows import LiterateWorkflow as Workflow
from ...interfaces import DerivativesDataSink


def init_fmap_estimation_wf(
epi_targets,
generate_report=True,
name="fmap_estimation_wf",
):
"""
Setup a fieldmap estimation strategy and write results to derivatives folder.
Parameters
----------
participant_label : :obj:`str`
The particular subject for which the BIDS layout will be queried.
epi_targets : :obj:`list` of :obj:`os.pathlike`
A list of :abbr:`EPI (echo planar imaging)` scans that will be corrected for
susceptibility distortions with the estimated fieldmaps.
omp_nthreads : :obj:`int`
Number of CPUs available to individual processes for multithreaded execution.
debug : :obj:`bool`
Whether fast (and less accurate) execution parameters should be used whenever available.
name : :obj:`str`
A unique workflow name to build Nipype's workflow hierarchy.
"""
layout = config.execution.layout
wf = Workflow(name=name)

inputnode = pe.Node(niu.IdentityInterface(fields=["dwi_reference", "dwi_mask"]),
name="inputnode")
# Create one outputnode with a port for each potential EPI target
outputnode = pe.Node(niu.IdentityInterface(fields=[_fname2outname(p) for p in epi_targets]),
name="outputnode")

# Return identity transforms for all if fieldmaps are ignored
if "fieldmaps" in config.workflow.ignore:
raise NotImplementedError

# Set-up PEPOLAR estimators only with EPIs under fmap/
# fmap_epi = {f: layout.get_metadata(f)
# for f in layout.get(
# subject=participant_label, datatype="fmap",
# suffix="epi", extension=("nii", "nii.gz"))}

metadata = [layout.get_metadata(p) for p in epi_targets]
pedirs = [m.get("PhaseEncodingDirection", "unknown") for m in metadata]
if len(set(pedirs) - set(("unknown",))) > 1:
if "unknown" in pedirs or len(set(pe[0] for pe in set(pedirs))) > 1:
raise NotImplementedError

# Get EPI polarities and their metadata
sdc_estimate_wf = init_pepolar_estimate_wf()
sdc_estimate_wf.inputs.inputnode.metadata = metadata

wf.connect([
(inputnode, sdc_estimate_wf, [("dwi_reference", "inputnode.in_data")]),
])
if generate_report:
from sdcflows.interfaces.reportlets import FieldmapReportlet
pepolar_report = pe.Node(FieldmapReportlet(reference_label="SDC'd B0"),
name="pepolar_report")
ds_report_pepolar = pe.Node(DerivativesDataSink(
base_directory=str(config.execution.work_dir / "reportlets"), keep_dtype=True,
desc="pepolar"), name="ds_report_pepolar")
ds_report_pepolar.inputs.source_file = epi_targets[0]
wf.connect([
(sdc_estimate_wf, pepolar_report, [
("outputnode.fieldmap", "fieldmap"),
("outputnode.corrected", "reference"),
("outputnode.corrected_mask", "mask")]),
(pepolar_report, ds_report_pepolar, [("out_report", "in_file")]),
])

return wf


def init_pepolar_estimate_wf(generate_report=True, name="pepolar_estimate_wf"):
"""Initialize a barebones TOPUP implementation."""
from nipype.interfaces.afni import Automask
from nipype.interfaces.fsl.epi import TOPUP
from niworkflows.interfaces.nibabel import MergeSeries
from ...interfaces.images import RescaleB0
wf = Workflow(name=name)

inputnode = pe.Node(niu.IdentityInterface(fields=["metadata", "in_data"]),
name="inputnode")
outputnode = pe.Node(niu.IdentityInterface(fields=["fieldmap", "corrected", "corrected_mask"]),
name="outputnode")

concat_blips = pe.Node(MergeSeries(), name="concat_blips")

topup = pe.Node(TOPUP(config=_pkg_fname("dmriprep", "data/flirtsch/b02b0.cnf")),
name="topup")

pre_mask = pe.Node(Automask(dilate=1, outputtype="NIFTI_GZ"),
name="pre_mask")
rescale_corrected = pe.Node(RescaleB0(), name="rescale_corrected")
post_mask = pe.Node(Automask(outputtype="NIFTI_GZ"),
name="post_mask")
wf.connect([
(inputnode, concat_blips, [("in_data", "in_files")]),
(inputnode, topup, [(("metadata", _get_ro), "readout_times"),
(("metadata", _get_pedir), "encoding_direction")]),
(concat_blips, topup, [("out_file", "in_file")]),
(topup, pre_mask, [("out_corrected", "in_file")]),
(pre_mask, rescale_corrected, [("out_file", "mask_file")]),
(topup, rescale_corrected, [("out_corrected", "in_file")]),
(topup, outputnode, [("out_field", "fieldmap")]),
(rescale_corrected, post_mask, [("out_ref", "in_file")]),
(rescale_corrected, outputnode, [("out_ref", "corrected")]),
(post_mask, outputnode, [("out_file", "corrected_mask")]),
])

return wf


def _get_ro(metadata):
return [m["TotalReadoutTime"] for m in metadata]


def _get_pedir(metadata):
return [m["PhaseEncodingDirection"].replace("j", "y").replace("i", "x").replace("k", "z")
for m in metadata]


def _fname2outname(in_file):
from pathlib import Path
return Path(in_file).name.rstrip(".gz").rstrip(".nii").replace("-", "_")
21 changes: 11 additions & 10 deletions setup.cfg
Original file line number Diff line number Diff line change
Expand Up @@ -45,16 +45,6 @@ packages = find:
exclude =
*.tests

[options.package_data]
dmriprep =
VERSION
config/reports-spec.yml
data/boilerplate.bib
data/tests/config.toml
data/tests/THP/*
data/tests/THP/sub-THP0005/anat/*
data/tests/THP/sub-THP0005/dwi/*

[options.exclude_package_data]
* = tests

Expand Down Expand Up @@ -93,6 +83,17 @@ all =
%(popylar)s
%(tests)s

[options.package_data]
dmriprep =
VERSION
config/reports-spec.yml
data/boilerplate.bib
data/flirtsch/*.cnf
data/tests/config.toml
data/tests/THP/*
data/tests/THP/sub-THP0005/anat/*
data/tests/THP/sub-THP0005/dwi/*

[options.entry_points]
console_scripts =
dmriprep=dmriprep.cli.run:main
Expand Down

0 comments on commit 61a9125

Please sign in to comment.