feat(TOPUP): an initial implementation of SD estimation.

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: #92

Author: oesteban

.docker/fsl-6.0/bin/topup

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

.docker/fsl-6.0/lib/libgfortran.so.3

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

.docker/fsl-6.0/lib/libopenblas.so.0

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

Dockerfile

@@ -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/config/reports-spec.yml

@@ -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.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_1.cnf

@@ -12,13 +12,14 @@
     BIDSInfo, BIDSFreeSurferDir
 )
 from niworkflows.utils.misc import fix_multi_T1w_source_name
+from niworkflows.utils.spaces import Reference
 from smriprep.workflows.anatomical import init_anat_preproc_wf
 
-from niworkflows.utils.spaces import Reference
 from ..interfaces import DerivativesDataSink, BIDSDataGrabber
 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():
@@ -253,21 +254,27 @@ def init_single_subject_wf(subject_id):
     anat_preproc_wf.__postdesc__ = (anat_preproc_wf.__postdesc__ or "") + f"""
 Diffusion data preprocessing
 
-: For each of the {len(subject_data["dwi"])} dwi scans found per subject
- (across all sessions), the following preprocessing was performed."""
+: For each of the {len(subject_data["dwi"])} DWI scans found per subject
+ (across all sessions), the gradient table was vetted and converted into the *RASb*
+format (i.e., given in RAS+ scanner coordinates, normalized b-vectors and scaled b-values),
+and a *b=0* average for reference to the subsequent steps of preprocessing was calculated.
+"""
 
     layout = config.execution.layout
+    dwi_data = tuple([
+        (dwi, layout.get_metadata(dwi), layout.get_bvec(dwi), layout.get_bval(dwi))
+        for dwi in subject_data["dwi"]
+    ])
     inputnode = pe.Node(niu.IdentityInterface(fields=["dwi_data"]),
                         name="inputnode")
-    inputnode.iterables = [(
-        "dwi_data", tuple([
-            (dwi, layout.get_bvec(dwi), layout.get_bval(dwi),
-             layout.get_metadata(dwi)["PhaseEncodingDirection"])
-            for dwi in subject_data["dwi"]
-        ])
-    )]
+    inputnode.iterables = [("dwi_data", dwi_data)]
+
+    referencenode = pe.JoinNode(niu.IdentityInterface(
+        fields=["dwi_file", "metadata", "dwi_reference", "dwi_mask", "gradients_rasb"]),
+        name="referencenode", joinsource="inputnode", run_without_submitting=True)
+
     split_info = pe.Node(niu.Function(
-        function=_unpack, output_names=["dwi_file", "bvec", "bval", "pedir"]),
+        function=_unpack, output_names=["dwi_file", "metadata", "bvec", "bval"]),
         name="split_info", run_without_submitting=True)
 
     early_b0ref_wf = init_early_b0ref_wf()
@@ -276,8 +283,21 @@ def init_single_subject_wf(subject_id):
         (split_info, early_b0ref_wf, [("dwi_file", "inputnode.dwi_file"),
                                       ("bvec", "inputnode.in_bvec"),
                                       ("bval", "inputnode.in_bval")]),
+        (split_info, referencenode, [("dwi_file", "dwi_file"),
+                                     ("metadata", "metadata")]),
+        (early_b0ref_wf, referencenode, [
+            ("outputnode.dwi_reference", "dwi_reference"),
+            ("outputnode.dwi_mask", "dwi_mask"),
+            ("outputnode.gradients_rasb", "gradients_rasb"),
+        ])
     ])
 
+    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