ENH: An initial implementation of SD estimation.

Dockerfile

@@ -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 && \

dmriprep/config/reports-spec.yml

@@ -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,run
+  reportlets:
+  - bids: {datatype: figures, desc: pepolar, suffix: fieldmap}
+    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:

dmriprep/data/flirtsch/b02b0.cnf

@@ -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

dmriprep/data/flirtsch/b02b0_1.cnf

@@ -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

dmriprep/data/flirtsch/b02b0_2.cnf

@@ -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

dmriprep/data/flirtsch/b02b0_4.cnf

@@ -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

dmriprep/data/flirtsch/b02b0_quick.cnf

@@ -0,0 +1,26 @@
+# Resolution (knot-spacing) of warps in mm
+--warpres=20,14,10
+# Subsampling level (a value of 2 indicates that a 2x2x2 neighbourhood is collapsed to 1 voxel)
+--subsamp=2,2,2
+# FWHM of gaussian smoothing
+--fwhm=8,2,1
+# Maximum number of iterations
+--miter=5,5,5
+# Relative weight of regularisation
+--lambda=0.005,0.000005,0.0000005
+# 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=0,0,0
+# 0=Levenberg-Marquardt, 1=Scaled Conjugate Gradient
+--minmet=0,0,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

dmriprep/workflows/base.py

@@ -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():
@@ -292,6 +293,13 @@ def init_single_subject_wf(subject_id):
         ]),
     ])
 
+    fmap_estimation_wf = init_fmap_estimation_wf(
+        subject_data["dwi"], debug=config.execution.debug)
+    workflow.connect([
+        (referencenode, fmap_estimation_wf, [
+            ("dwi_reference", "inputnode.dwi_reference"),
+            ("dwi_mask", "inputnode.dwi_mask")]),
+    ])
     return workflow
 
 

dmriprep/workflows/fmap/base.py

@@ -0,0 +1,145 @@
+"""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,
+    debug=False,
+    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")
+    wf.add_nodes([inputnode])  # TODO: remove when fully implemented
+    # 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:
+        return wf
+
+    # 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]
+    if any("TotalReadoutTime" not in m for m in metadata):
+        return wf
+
+    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(debug=debug)
+        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.output_dir), datatype="figures",
+                suffix="fieldmap", desc="pepolar", dismiss_entities=("acquisition", "dir")),
+                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(debug=False, 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 sdcflows.interfaces.fmap import get_trt
+    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")
+    readout_time = pe.MapNode(niu.Function(
+        input_names=["in_meta", "in_file"], function=get_trt), name="readout_time",
+        iterfield=["in_meta", "in_file"], run_without_submitting=True
+    )
+
+    topup = pe.Node(TOPUP(config=_pkg_fname(
+        "dmriprep", f"data/flirtsch/b02b0{'_quick' * debug}.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, readout_time, [("in_data", "in_file"),
+                                   ("metadata", "in_meta")]),
+        (inputnode, topup, [(("metadata", _get_pedir), "encoding_direction")]),
+        (readout_time, topup, [("out", "readout_times")]),
+        (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_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("-", "_")

setup.cfg

@@ -50,6 +50,7 @@ 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/*