Skip to content

Commit

Permalink
Merge pull request #50 from josephmje/enh/skullstrip
Browse files Browse the repository at this point in the history
ENH: b0 reference and skullstrip workflow
  • Loading branch information
oesteban committed Jan 15, 2020
2 parents 5bd7fcf + 51da317 commit c695f4c
Show file tree
Hide file tree
Showing 11 changed files with 687 additions and 107 deletions.
18 changes: 18 additions & 0 deletions .circleci/config.yml
Expand Up @@ -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:
Expand Down
1 change: 1 addition & 0 deletions dmriprep/config/__init__.py
Expand Up @@ -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']
10 changes: 10 additions & 0 deletions dmriprep/config/reports-spec.yml
Expand Up @@ -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}
Empty file.
145 changes: 145 additions & 0 deletions dmriprep/interfaces/images.py
@@ -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
43 changes: 0 additions & 43 deletions dmriprep/interfaces/reports.py
Expand Up @@ -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>
Expand Down Expand Up @@ -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)
102 changes: 38 additions & 64 deletions dmriprep/workflows/base.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
Expand All @@ -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(
Expand Down Expand Up @@ -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

Expand Down
9 changes: 9 additions & 0 deletions dmriprep/workflows/dwi/__init__.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',
]

0 comments on commit c695f4c

Please sign in to comment.