Skip to content

Commit fb83460

Browse files
josephmjeoesteban
josephmje
authored and
oesteban
committed
enh: initial implementation of b0 reference and skull-stripping workflow
initial port from niworkflows add initial pre-mask step fix doc add pre_mask address review comments add period extract b0s first fix function name typo initial port from niworkflows add initial pre-mask step fix doc add pre_mask address review comments add period extract b0s first fix function name typo remove MatchHeader
1 parent 5282950 commit fb83460

File tree

2 files changed

+312
-0
lines changed

2 files changed

+312
-0
lines changed

‎dmriprep/interfaces/images.py

+62
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,62 @@
1+
#!/usr/bin/env python
2+
# -*- coding: utf-8 -*-
3+
# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
4+
# vi: set ft=python sts=4 ts=4 sw=4 et:
5+
"""
6+
Image tools interfaces
7+
~~~~~~~~~~~~~~~~~~~~~~
8+
9+
"""
10+
11+
from nipype import logging
12+
from nipype.interfaces.base import (
13+
TraitedSpec, BaseInterfaceInputSpec, SimpleInterface, File
14+
)
15+
16+
LOGGER = logging.getLogger('nipype.interface')
17+
18+
19+
class ExtractB0InputSpec(BaseInterfaceInputSpec):
20+
in_file = File(exists=True, mandatory=True, desc='dwi file')
21+
in_rasb = File(exists=True, mandatory=True,
22+
desc='File containing the gradient directions in RAS+ scanner coordinates')
23+
24+
25+
class ExtractB0OutputSpec(TraitedSpec):
26+
out_file = File(exists=True, desc='mean b0 file')
27+
28+
29+
class ExtractB0(SimpleInterface):
30+
input_spec = ExtractB0InputSpec
31+
output_spec = ExtractB0OutputSpec
32+
33+
def _run_interface(self, runtime):
34+
self._results['out_file'] = extract_b0(
35+
self.inputs.in_file,
36+
self.inputs.in_rasb,
37+
newpath=runtime.cwd)
38+
return runtime
39+
40+
41+
def extract_b0(in_file, in_rasb, b0_thres=50, newpath=None):
42+
"""Extract the *b0* volumes from a DWI dataset."""
43+
import numpy as np
44+
import nibabel as nib
45+
from nipype.utils.filemanip import fname_presuffix
46+
47+
out_file = fname_presuffix(
48+
in_file, suffix='_b0', newpath=newpath)
49+
50+
img = nib.load(in_file)
51+
data = img.get_fdata()
52+
53+
bvals = np.loadtxt(in_rasb, usecols=-1, skiprows=1)
54+
55+
b0 = data[..., bvals < b0_thres]
56+
57+
hdr = img.header.copy()
58+
hdr.set_data_shape(b0.shape)
59+
hdr.set_xyzt_units('mm')
60+
hdr.set_data_dtype(np.float32)
61+
nib.Nifti1Image(b0, img.affine, hdr).to_filename(out_file)
62+
return out_file

‎dmriprep/workflows/dwi/util.py

+250
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,250 @@
1+
# -*- coding: utf-8 -*-
2+
3+
"""
4+
Utility workflows
5+
^^^^^^^^^^^^^^^^^
6+
7+
.. autofunction:: init_dwi_reference_wf
8+
.. autofunction:: init_enhance_and_skullstrip_dwi_wf
9+
10+
"""
11+
12+
from nipype.pipeline import engine as pe
13+
from nipype.interfaces import utility as niu, fsl, afni
14+
15+
from niworkflows.engine.workflows import LiterateWorkflow as Workflow
16+
from niworkflows.interfaces.images import ValidateImage
17+
from niworkflows.interfaces.fixes import FixN4BiasFieldCorrection as N4BiasFieldCorrection
18+
from niworkflows.interfaces.masks import SimpleShowMaskRPT
19+
from niworkflows.interfaces.utils import CopyXForm
20+
21+
from ...interfaces.images import MeanB0
22+
23+
DEFAULT_MEMORY_MIN_GB = 0.01
24+
25+
26+
def init_dwi_reference_wf(omp_nthreads, dwi_file=None,
27+
name='dwi_reference_wf', gen_report=False):
28+
"""
29+
Build a workflow that generates a reference b0 image from a dwi series.
30+
31+
.. workflow::
32+
:graph2use: orig
33+
:simple_form: yes
34+
35+
from dmriprep.workflows.dwi.util import init_dwi_reference_wf
36+
wf = init_dwi_reference_wf(omp_nthreads=1)
37+
38+
**Parameters**
39+
40+
dwi_file : str
41+
dwi NIfTI file
42+
omp_nthreads : int
43+
Maximum number of threads an individual process may use
44+
name : str
45+
Name of workflow (default: ``dwi_reference_wf``)
46+
gen_report : bool
47+
Whether a mask report node should be appended in the end
48+
49+
**Inputs**
50+
51+
dwi_file
52+
dwi NIfTI file
53+
54+
**Outputs**
55+
56+
dwi_file
57+
Validated dwi NIfTI file
58+
raw_ref_image
59+
Reference image to which dwi series is motion corrected
60+
ref_image
61+
Contrast-enhanced reference image
62+
ref_image_brain
63+
Skull-stripped reference image
64+
dwi_mask
65+
Skull-stripping mask of reference image
66+
validation_report
67+
HTML reportlet indicating whether ``dwi_file`` had a valid affine
68+
69+
70+
**Subworkflows**
71+
72+
* :py:func:`~dmriprep.workflows.dwi.util.init_enhance_and_skullstrip_wf`
73+
74+
"""
75+
workflow = Workflow(name=name)
76+
workflow.__desc__ = """\
77+
First, a reference volume and its skull-stripped version were generated
78+
using a custom methodology taken from *fMRIPrep*.
79+
"""
80+
inputnode = pe.Node(niu.IdentityInterface(fields=['dwi_file', 'bvec_file', 'bval_file']),
81+
name='inputnode')
82+
outputnode = pe.Node(
83+
niu.IdentityInterface(fields=['dwi_file', 'raw_ref_image',
84+
'ref_image', 'ref_image_brain',
85+
'dwi_mask', 'validation_report']),
86+
name='outputnode')
87+
88+
# Simplify manually setting input image
89+
if dwi_file is not None:
90+
inputnode.inputs.dwi_file = dwi_file
91+
92+
validate = pe.Node(ValidateImage(), name='validate', mem_gb=DEFAULT_MEMORY_MIN_GB)
93+
94+
gen_ref = pe.Node(MeanB0(), name="gen_ref")
95+
96+
pre_mask = pe.Node(afni.Automask(outputtype="NIFTI_GZ"), name="pre_mask")
97+
98+
enhance_and_skullstrip_dwi_wf = init_enhance_and_skullstrip_dwi_wf(
99+
omp_nthreads=omp_nthreads)
100+
101+
workflow.connect([
102+
(inputnode, validate, [('dwi_file', 'in_file')]),
103+
(validate, gen_ref, [('out_file', 'in_file')]),
104+
(inputnode, gen_ref, [('bvec_file', 'in_bvec'),
105+
('bval_file', 'in_bval')]),
106+
(gen_ref, pre_mask, [('out_file', 'in_file')]),
107+
(gen_ref, enhance_and_skullstrip_dwi_wf, [('out_file', 'inputnode.in_file')]),
108+
(pre_mask, enhance_and_skullstrip_dwi_wf, [('out_file', 'inputnode.pre_mask')]),
109+
(validate, outputnode, [('out_file', 'dwi_file'),
110+
('out_report', 'validation_report')]),
111+
(gen_ref, outputnode, [('out_file', 'raw_ref_image')]),
112+
(enhance_and_skullstrip_dwi_wf, outputnode, [
113+
('outputnode.bias_corrected_file', 'ref_image'),
114+
('outputnode.mask_file', 'dwi_mask'),
115+
('outputnode.skull_stripped_file', 'ref_image_brain')]),
116+
])
117+
118+
if gen_report:
119+
mask_reportlet = pe.Node(SimpleShowMaskRPT(), name='mask_reportlet')
120+
workflow.connect([
121+
(enhance_and_skullstrip_dwi_wf, mask_reportlet, [
122+
('outputnode.bias_corrected_file', 'background_file'),
123+
('outputnode.mask_file', 'mask_file'),
124+
]),
125+
])
126+
127+
return workflow
128+
129+
130+
def init_enhance_and_skullstrip_dwi_wf(
131+
name='enhance_and_skullstrip_dwi_wf',
132+
omp_nthreads=1):
133+
"""
134+
This workflow takes in a dwi reference iamge and sharpens the histogram
135+
with the application of the N4 algorithm for removing the
136+
:abbr:`INU (intensity non-uniformity)` bias field and calculates a signal
137+
mask.
138+
139+
Steps of this workflow are:
140+
141+
1. Run ANTs' ``N4BiasFieldCorrection`` on the input
142+
dwi reference image and mask.
143+
2. Calculate a loose mask using FSL's ``bet``, with one mathematical morphology
144+
dilation of one iteration and a sphere of 6mm as structuring element.
145+
3. Mask the :abbr:`INU (intensity non-uniformity)`-corrected image
146+
with the latest mask calculated in 3), then use AFNI's ``3dUnifize``
147+
to *standardize* the T2* contrast distribution.
148+
4. Calculate a mask using AFNI's ``3dAutomask`` after the contrast
149+
enhancement of 4).
150+
5. Calculate a final mask as the intersection of 2) and 4).
151+
6. Apply final mask on the enhanced reference.
152+
153+
.. workflow ::
154+
:graph2use: orig
155+
:simple_form: yes
156+
157+
from dmriprep.workflows.dwi.util import init_enhance_and_skullstrip_dwi_wf
158+
wf = init_enhance_and_skullstrip_dwi_wf(omp_nthreads=1)
159+
160+
**Parameters**
161+
name : str
162+
Name of workflow (default: ``enhance_and_skullstrip_dwi_wf``)
163+
omp_nthreads : int
164+
number of threads available to parallel nodes
165+
166+
**Inputs**
167+
168+
in_file
169+
b0 image (single volume)
170+
pre_mask
171+
initial mask
172+
173+
**Outputs**
174+
175+
bias_corrected_file
176+
the ``in_file`` after `N4BiasFieldCorrection`_
177+
skull_stripped_file
178+
the ``bias_corrected_file`` after skull-stripping
179+
mask_file
180+
mask of the skull-stripped input file
181+
out_report
182+
reportlet for the skull-stripping
183+
184+
.. _N4BiasFieldCorrection: https://hdl.handle.net/10380/3053
185+
"""
186+
workflow = Workflow(name=name)
187+
inputnode = pe.Node(niu.IdentityInterface(fields=['in_file', 'pre_mask']),
188+
name='inputnode')
189+
outputnode = pe.Node(niu.IdentityInterface(fields=[
190+
'mask_file', 'skull_stripped_file', 'bias_corrected_file']), name='outputnode')
191+
192+
# Run N4 normally, force num_threads=1 for stability (images are small, no need for >1)
193+
n4_correct = pe.Node(N4BiasFieldCorrection(
194+
dimension=3, copy_header=True, bspline_fitting_distance=200), shrink_factor=2,
195+
name='n4_correct', n_procs=1)
196+
n4_correct.inputs.rescale_intensities = True
197+
198+
# Create a generous BET mask out of the bias-corrected EPI
199+
skullstrip_first_pass = pe.Node(fsl.BET(frac=0.2, mask=True),
200+
name='skullstrip_first_pass')
201+
bet_dilate = pe.Node(fsl.DilateImage(
202+
operation='max', kernel_shape='sphere', kernel_size=6.0,
203+
internal_datatype='char'), name='skullstrip_first_dilate')
204+
bet_mask = pe.Node(fsl.ApplyMask(), name='skullstrip_first_mask')
205+
206+
# Use AFNI's unifize for T2 constrast & fix header
207+
unifize = pe.Node(afni.Unifize(
208+
t2=True, outputtype='NIFTI_GZ',
209+
# Default -clfrac is 0.1, 0.4 was too conservative
210+
# -rbt because I'm a Jedi AFNI Master (see 3dUnifize's documentation)
211+
args='-clfrac 0.2 -rbt 18.3 65.0 90.0',
212+
out_file="uni.nii.gz"), name='unifize')
213+
fixhdr_unifize = pe.Node(CopyXForm(), name='fixhdr_unifize', mem_gb=0.1)
214+
215+
# Run ANFI's 3dAutomask to extract a refined brain mask
216+
skullstrip_second_pass = pe.Node(afni.Automask(dilate=1,
217+
outputtype='NIFTI_GZ'),
218+
name='skullstrip_second_pass')
219+
fixhdr_skullstrip2 = pe.Node(CopyXForm(), name='fixhdr_skullstrip2', mem_gb=0.1)
220+
221+
# Take intersection of both masks
222+
combine_masks = pe.Node(fsl.BinaryMaths(operation='mul'),
223+
name='combine_masks')
224+
225+
# Compute masked brain
226+
apply_mask = pe.Node(fsl.ApplyMask(), name='apply_mask')
227+
228+
workflow.connect([
229+
(inputnode, n4_correct, [('in_file', 'input_image'),
230+
('pre_mask', 'mask_image')]),
231+
(inputnode, fixhdr_unifize, [('in_file', 'hdr_file')]),
232+
(inputnode, fixhdr_skullstrip2, [('in_file', 'hdr_file')]),
233+
(n4_correct, skullstrip_first_pass, [('output_image', 'in_file')]),
234+
(skullstrip_first_pass, bet_dilate, [('mask_file', 'in_file')]),
235+
(bet_dilate, bet_mask, [('out_file', 'mask_file')]),
236+
(skullstrip_first_pass, bet_mask, [('out_file', 'in_file')]),
237+
(bet_mask, unifize, [('out_file', 'in_file')]),
238+
(unifize, fixhdr_unifize, [('out_file', 'in_file')]),
239+
(fixhdr_unifize, skullstrip_second_pass, [('out_file', 'in_file')]),
240+
(skullstrip_first_pass, combine_masks, [('mask_file', 'in_file')]),
241+
(skullstrip_second_pass, fixhdr_skullstrip2, [('out_file', 'in_file')]),
242+
(fixhdr_skullstrip2, combine_masks, [('out_file', 'operand_file')]),
243+
(fixhdr_unifize, apply_mask, [('out_file', 'in_file')]),
244+
(combine_masks, apply_mask, [('out_file', 'mask_file')]),
245+
(combine_masks, outputnode, [('out_file', 'mask_file')]),
246+
(apply_mask, outputnode, [('out_file', 'skull_stripped_file')]),
247+
(n4_correct, outputnode, [('output_image', 'bias_corrected_file')]),
248+
])
249+
250+
return workflow

0 commit comments

Comments
 (0)