Merge pull request #50 from josephmje/enh/skullstrip

ENH: b0 reference and skullstrip workflow

Author: oesteban

.circleci/config.yml

@@ -251,6 +251,24 @@ jobs:
          key: THP002-anat-v00-{{ .Branch }}-{{ .Revision }}-{{ epoch }}
          paths:
             - /tmp/THP002/work
+      - run:
+          name: Run full diffusion workflow on THP002
+          no_output_timeout: 2h
+          command: |
+            mkdir -p /tmp/THP002/work /tmp/THP002/derivatives
+            sudo setfacl -d -m group:$(id -gn):rwx /tmp/THP002/derivatives && \
+                sudo setfacl -m group:$(id -gn):rwx /tmp/THP002/derivatives
+            sudo setfacl -d -m group:$(id -gn):rwx /tmp/THP002/work && \
+                sudo setfacl -m group:$(id -gn):rwx /tmp/THP002/work
+            docker run -e FS_LICENSE=$FS_LICENSE --rm \
+                -v /tmp/data/THP002:/data \
+                -v /tmp/THP002/derivatives:/out \
+                -v /tmp/fslicense/license.txt:/tmp/fslicense/license.txt:ro \
+                -v /tmp/config/nipype.cfg:/home/dmriprep/.nipype/nipype.cfg \
+                -v /tmp/THP002/work:/work \
+                --user $(id -u):$(id -g) \
+                nipreps/dmriprep:latest /data /out participant -vv --sloppy \
+                --notrack --skip-bids-validation -w /work --omp-nthreads 2 --nprocs 2
       - store_artifacts:
           path: /tmp/THP002/derivatives/dmriprep
       - run:

dmriprep/config/__init__.py

@@ -2,4 +2,5 @@
 # vi: set ft=python sts=4 ts=4 sw=4 et:
 """Settings."""
 
+DEFAULT_MEMORY_MIN_GB = 0.01
 NONSTANDARD_REFERENCES = ['anat', 'T1w', 'dwi', 'fsnative']

dmriprep/config/reports-spec.yml

@@ -26,6 +26,16 @@ sections:
     caption: Surfaces (white and pial) reconstructed with FreeSurfer (<code>recon-all</code>)
       overlaid on the participant's T1w template.
     subtitle: Surface reconstruction
+- name: Diffusion
+  ordering: session,acquisition,run
+  reportlets:
+  - bids: {datatype: dwi, desc: validation, suffix: dwi}
+  - bids: {datatype: dwi, desc: brain, suffix: mask}
+    caption: Average b=0 that serves for reference in early preprocessing steps.
+    descriptions: The reference b=0 is obtained as the voxel-wise median across
+      all b=0 found in the dataset, after accounting for signal drift.
+      The red contour shows the brain mask calculated using this reference b=0.
+    subtitle: Reference b=0 and brain mask
 - name: About
   reportlets:
   - bids: {datatype: anat, desc: about, suffix: T1w}

dmriprep/data/tests/dwi_mask.nii.gz

@@ -0,0 +1,145 @@
+"""Image tools interfaces."""
+import numpy as np
+import nibabel as nb
+from nipype.utils.filemanip import fname_presuffix
+from nipype import logging
+from nipype.interfaces.base import (
+    traits, TraitedSpec, BaseInterfaceInputSpec, SimpleInterface, File
+)
+
+LOGGER = logging.getLogger('nipype.interface')
+
+
+class _ExtractB0InputSpec(BaseInterfaceInputSpec):
+    in_file = File(exists=True, mandatory=True, desc='dwi file')
+    b0_ixs = traits.List(traits.Int, mandatory=True,
+                         desc='Index of b0s')
+
+
+class _ExtractB0OutputSpec(TraitedSpec):
+    out_file = File(exists=True, desc='b0 file')
+
+
+class ExtractB0(SimpleInterface):
+    """
+    Extract all b=0 volumes from a dwi series.
+
+    Example
+    -------
+    >>> os.chdir(tmpdir)
+    >>> extract_b0 = ExtractB0()
+    >>> extract_b0.inputs.in_file = str(data_dir / 'dwi.nii.gz')
+    >>> extract_b0.inputs.b0_ixs = [0, 1, 2]
+    >>> res = extract_b0.run()  # doctest: +SKIP
+
+    """
+
+    input_spec = _ExtractB0InputSpec
+    output_spec = _ExtractB0OutputSpec
+
+    def _run_interface(self, runtime):
+        self._results['out_file'] = extract_b0(
+            self.inputs.in_file,
+            self.inputs.b0_ixs,
+            newpath=runtime.cwd)
+        return runtime
+
+
+def extract_b0(in_file, b0_ixs, newpath=None):
+    """Extract the *b0* volumes from a DWI dataset."""
+    out_file = fname_presuffix(
+        in_file, suffix='_b0', newpath=newpath)
+
+    img = nb.load(in_file)
+    data = img.get_fdata(dtype='float32')
+
+    b0 = data[..., b0_ixs]
+
+    hdr = img.header.copy()
+    hdr.set_data_shape(b0.shape)
+    hdr.set_xyzt_units('mm')
+    hdr.set_data_dtype(np.float32)
+    nb.Nifti1Image(b0, img.affine, hdr).to_filename(out_file)
+    return out_file
+
+
+class _RescaleB0InputSpec(BaseInterfaceInputSpec):
+    in_file = File(exists=True, mandatory=True, desc='b0s file')
+    mask_file = File(exists=True, mandatory=True, desc='mask file')
+
+
+class _RescaleB0OutputSpec(TraitedSpec):
+    out_ref = File(exists=True, desc='One average b0 file')
+    out_b0s = File(exists=True, desc='series of rescaled b0 volumes')
+
+
+class RescaleB0(SimpleInterface):
+    """
+    Rescale the b0 volumes to deal with average signal decay over time.
+
+    Example
+    -------
+    >>> os.chdir(tmpdir)
+    >>> rescale_b0 = RescaleB0()
+    >>> rescale_b0.inputs.in_file = str(data_dir / 'dwi.nii.gz')
+    >>> rescale_b0.inputs.mask_file = str(data_dir / 'dwi_mask.nii.gz')
+    >>> res = rescale_b0.run()  # doctest: +SKIP
+
+    """
+
+    input_spec = _RescaleB0InputSpec
+    output_spec = _RescaleB0OutputSpec
+
+    def _run_interface(self, runtime):
+        self._results['out_b0s'] = rescale_b0(
+            self.inputs.in_file,
+            self.inputs.mask_file,
+            newpath=runtime.cwd
+        )
+        self._results['out_ref'] = median(
+            self._results['out_b0s'],
+            newpath=runtime.cwd
+        )
+        return runtime
+
+
+def rescale_b0(in_file, mask_file, newpath=None):
+    """Rescale the input volumes using the median signal intensity."""
+    out_file = fname_presuffix(
+        in_file, suffix='_rescaled_b0', newpath=newpath)
+
+    img = nb.load(in_file)
+    if img.dataobj.ndim == 3:
+        return in_file
+
+    data = img.get_fdata(dtype='float32')
+    mask_img = nb.load(mask_file)
+    mask_data = mask_img.get_fdata(dtype='float32')
+
+    median_signal = np.median(data[mask_data > 0, ...], axis=0)
+    rescaled_data = 1000 * data / median_signal
+    hdr = img.header.copy()
+    nb.Nifti1Image(rescaled_data, img.affine, hdr).to_filename(out_file)
+    return out_file
+
+
+def median(in_file, newpath=None):
+    """Average a 4D dataset across the last dimension using median."""
+    out_file = fname_presuffix(
+        in_file, suffix='_b0ref', newpath=newpath)
+
+    img = nb.load(in_file)
+    if img.dataobj.ndim == 3:
+        return in_file
+    if img.shape[-1] == 1:
+        nb.squeeze_image(img).to_filename(out_file)
+        return out_file
+
+    median_data = np.median(img.get_fdata(dtype='float32'),
+                            axis=-1)
+
+    hdr = img.header.copy()
+    hdr.set_xyzt_units('mm')
+    hdr.set_data_dtype(np.float32)
+    nb.Nifti1Image(median_data, img.affine, hdr).to_filename(out_file)
+    return out_file

dmriprep/interfaces/images.py

@@ -23,13 +23,6 @@
 \t</ul>
 """
 
-DWI_TEMPLATE = """\t\t<h3 class="elem-title">Summary</h3>
-\t\t<ul class="elem-desc">
-\t\t\t<li>Susceptibility distortion correction: {sdc}</li>
-\t\t\t<li>Registration: {registration}</li>
-\t\t</ul>
-"""
-
 ABOUT_TEMPLATE = """\t<ul>
 \t\t<li>dMRIPrep version: {version}</li>
 \t\t<li>dMRIPrep command: <code>{command}</code></li>
@@ -128,39 +121,3 @@ def _generate_segment(self):
         return ABOUT_TEMPLATE.format(version=self.inputs.version,
                                      command=self.inputs.command,
                                      date=time.strftime("%Y-%m-%d %H:%M:%S %z"))
-
-
-class DiffusionSummaryInputSpec(BaseInterfaceInputSpec):
-    distortion_correction = traits.Str(desc='Susceptibility distortion correction method',
-                                       mandatory=True)
-    pe_direction = traits.Enum(None, 'i', 'i-', 'j', 'j-', mandatory=True,
-                               desc='Phase-encoding direction detected')
-    registration = traits.Enum('FSL', 'FreeSurfer', mandatory=True,
-                               desc='Diffusion/anatomical registration method')
-    fallback = traits.Bool(desc='Boundary-based registration rejected')
-    registration_dof = traits.Enum(6, 9, 12, desc='Registration degrees of freedom',
-                                   mandatory=True)
-
-
-class DiffusionSummary(SummaryInterface):
-    input_spec = DiffusionSummaryInputSpec
-
-    def _generate_segment(self):
-        dof = self.inputs.registration_dof
-        reg = {
-            'FSL': [
-                'FSL <code>flirt</code> with boundary-based registration'
-                ' (BBR) metric - %d dof' % dof,
-                'FSL <code>flirt</code> rigid registration - 6 dof'],
-            'FreeSurfer': [
-                'FreeSurfer <code>bbregister</code> '
-                '(boundary-based registration, BBR) - %d dof' % dof,
-                'FreeSurfer <code>mri_coreg</code> - %d dof' % dof],
-        }[self.inputs.registration][self.inputs.fallback]
-        if self.inputs.pe_direction is None:
-            pedir = 'MISSING - Assuming Anterior-Posterior'
-        else:
-            pedir = {'i': 'Left-Right', 'j': 'Anterior-Posterior'}[self.inputs.pe_direction[0]]
-
-        return DWI_TEMPLATE.format(
-            pedir=pedir, sdc=self.inputs.distortion_correction, registration=reg)

dmriprep/interfaces/reports.py

@@ -1,13 +1,4 @@
-# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
-# vi: set ft=python sts=4 ts=4 sw=4 et:
-"""
-dMRIPrep base processing workflows.
-
-.. autofunction:: init_dmriprep_wf
-.. autofunction:: init_single_subject_wf
-
-"""
-
+"""dMRIPrep base processing workflows."""
 import sys
 import os
 from packaging.version import Version
@@ -29,7 +20,7 @@
 from ..interfaces.reports import SubjectSummary, AboutSummary
 from ..utils.bids import collect_data
 from ..__about__ import __version__
-# from .dwi import init_dwi_preproc_wf
+from .dwi import init_dwi_preproc_wf
 
 
 def init_dmriprep_wf(
@@ -435,59 +426,42 @@ def init_single_subject_wf(
     if anat_only:
         return workflow
 
-    # for dwi_file in subject_data['dwi']:
-    #     dwi_preproc_wf = init_dwi_preproc_wf(
-    #         aroma_melodic_dim=aroma_melodic_dim,
-    #         bold2t1w_dof=bold2t1w_dof,
-    #         bold_file=bold_file,
-    #         cifti_output=cifti_output,
-    #         debug=debug,
-    #         dummy_scans=dummy_scans,
-    #         err_on_aroma_warn=err_on_aroma_warn,
-    #         fmap_bspline=fmap_bspline,
-    #         fmap_demean=fmap_demean,
-    #         force_syn=force_syn,
-    #         freesurfer=freesurfer,
-    #         ignore=ignore,
-    #         layout=layout,
-    #         low_mem=low_mem,
-    #         medial_surface_nan=medial_surface_nan,
-    #         num_bold=len(subject_data['bold']),
-    #         omp_nthreads=omp_nthreads,
-    #         output_dir=output_dir,
-    #         output_spaces=output_spaces,
-    #         reportlets_dir=reportlets_dir,
-    #         regressors_all_comps=regressors_all_comps,
-    #         regressors_fd_th=regressors_fd_th,
-    #         regressors_dvars_th=regressors_dvars_th,
-    #         t2s_coreg=t2s_coreg,
-    #         use_aroma=use_aroma,
-    #         use_syn=use_syn,
-    #     )
-
-    #     workflow.connect([
-    #         (anat_preproc_wf, dwi_preproc_wf,
-    #          [(('outputnode.t1_preproc', _pop), 'inputnode.t1_preproc'),
-    #           ('outputnode.t1_brain', 'inputnode.t1_brain'),
-    #           ('outputnode.t1_mask', 'inputnode.t1_mask'),
-    #           ('outputnode.t1_seg', 'inputnode.t1_seg'),
-    #           ('outputnode.t1_aseg', 'inputnode.t1_aseg'),
-    #           ('outputnode.t1_aparc', 'inputnode.t1_aparc'),
-    #           ('outputnode.t1_tpms', 'inputnode.t1_tpms'),
-    #           ('outputnode.template', 'inputnode.template'),
-    #           ('outputnode.forward_transform', 'inputnode.anat2std_xfm'),
-    #           ('outputnode.reverse_transform', 'inputnode.std2anat_xfm'),
-    #           ('outputnode.joint_template', 'inputnode.joint_template'),
-    #           ('outputnode.joint_forward_transform', 'inputnode.joint_anat2std_xfm'),
-    #           ('outputnode.joint_reverse_transform', 'inputnode.joint_std2anat_xfm'),
-    #           # Undefined if --no-freesurfer, but this is safe
-    #           ('outputnode.subjects_dir', 'inputnode.subjects_dir'),
-    #           ('outputnode.subject_id', 'inputnode.subject_id'),
-    #           ('outputnode.t1_2_fsnative_forward_transform',
-    #            'inputnode.t1_2_fsnative_forward_transform'),
-    #           ('outputnode.t1_2_fsnative_reverse_transform',
-    #            'inputnode.t1_2_fsnative_reverse_transform')]),
-    #     ])
+    for dwi_file in subject_data['dwi']:
+        dwi_preproc_wf = init_dwi_preproc_wf(
+            dwi_file=dwi_file,
+            debug=debug,
+            force_syn=force_syn,
+            ignore=ignore,
+            layout=layout,
+            low_mem=low_mem,
+            num_dwi=len(subject_data['dwi']),
+            omp_nthreads=omp_nthreads,
+            output_dir=output_dir,
+            reportlets_dir=reportlets_dir,
+            use_syn=use_syn,
+        )
+
+        workflow.connect([
+            (anat_preproc_wf, dwi_preproc_wf,
+             [(('outputnode.t1w_preproc', _pop), 'inputnode.t1w_preproc'),
+              ('outputnode.t1w_brain', 'inputnode.t1w_brain'),
+              ('outputnode.t1w_mask', 'inputnode.t1w_mask'),
+              ('outputnode.t1w_dseg', 'inputnode.t1w_dseg'),
+              ('outputnode.t1w_aseg', 'inputnode.t1w_aseg'),
+              ('outputnode.t1w_aparc', 'inputnode.t1w_aparc'),
+              ('outputnode.t1w_tpms', 'inputnode.t1w_tpms'),
+              ('outputnode.template', 'inputnode.template'),
+              ('outputnode.anat2std_xfm', 'inputnode.anat2std_xfm'),
+              ('outputnode.std2anat_xfm', 'inputnode.std2anat_xfm'),
+              ('outputnode.joint_template', 'inputnode.joint_template'),
+              ('outputnode.joint_anat2std_xfm', 'inputnode.joint_anat2std_xfm'),
+              ('outputnode.joint_std2anat_xfm', 'inputnode.joint_std2anat_xfm'),
+              # Undefined if --fs-no-reconall, but this is safe
+              ('outputnode.subjects_dir', 'inputnode.subjects_dir'),
+              ('outputnode.subject_id', 'inputnode.subject_id'),
+              ('outputnode.t1w2fsnative_xfm', 'inputnode.t1w2fsnative_xfm'),
+              ('outputnode.fsnative2t1w_xfm', 'inputnode.fsnative2t1w_xfm')]),
+        ])
 
     return workflow
 

dmriprep/workflows/base.py

@@ -0,0 +1,9 @@
+"""Pre-processing dMRI workflows."""
+
+from .base import init_dwi_preproc_wf
+from .util import init_dwi_reference_wf
+
+__all__ = [
+    'init_dwi_preproc_wf',
+    'init_dwi_reference_wf',
+]