Skip to content

Commit eba2228

Browse files
oestebanoesteban
committed
ENH: Deduplicating magnitude handling and fieldmap postprocessing workflows
This PR splits the magnitude processing and the fieldmap postprocessing that is shared between phasediff (and phase1/2) and direct B0 mapping approaches. It also addresses some issues on the direct B0 mapping approach, where the fieldmap was converted between units and subject to PRELUDE (which is not correct as these fieldmaps are not phase-wrapped). Adding a test case of SE fieldmaps would be great for a future PR. Close #17
1 parent be7058e commit eba2228

File tree

3 files changed

+253
-108
lines changed

3 files changed

+253
-108
lines changed

‎sdcflows/workflows/fmap.py

Lines changed: 19 additions & 66 deletions
Original file line numberDiff line numberDiff line change
@@ -19,16 +19,11 @@
1919
"""
2020

2121
from nipype.pipeline import engine as pe
22-
from nipype.interfaces import utility as niu, fsl, ants
23-
from niflow.nipype1.workflows.dmri.fsl.utils import demean_image, cleanup_edge_pipeline
22+
from nipype.interfaces import utility as niu, fsl
2423
from niworkflows.engine.workflows import LiterateWorkflow as Workflow
25-
from niworkflows.interfaces.bids import DerivativesDataSink
2624
from niworkflows.interfaces.images import IntraModalMerge
27-
from niworkflows.interfaces.masks import BETRPT
2825

29-
from ..interfaces.fmap import (
30-
FieldEnhance, FieldToRadS, FieldToHz
31-
)
26+
from .gre import init_fmap_postproc_wf, init_magnitude_wf
3227

3328

3429
def init_fmap_wf(omp_nthreads, fmap_bspline, name='fmap_wf'):
@@ -81,70 +76,28 @@ def init_fmap_wf(omp_nthreads, fmap_bspline, name='fmap_wf'):
8176
outputnode = pe.Node(niu.IdentityInterface(fields=['fmap', 'fmap_ref', 'fmap_mask']),
8277
name='outputnode')
8378

84-
# Merge input magnitude images
85-
magmrg = pe.Node(IntraModalMerge(), name='magmrg')
79+
magnitude_wf = init_magnitude_wf(omp_nthreads=omp_nthreads)
80+
workflow.connect([
81+
(inputnode, magnitude_wf, [('magnitude', 'inputnode.magnitude')]),
82+
(magnitude_wf, outputnode, [('outputnode.fmap_mask', 'fmap_mask'),
83+
('outputnode.fmap_ref', 'fmap_ref')]),
84+
])
85+
8686
# Merge input fieldmap images
8787
fmapmrg = pe.Node(IntraModalMerge(zero_based_avg=False, hmc=False),
8888
name='fmapmrg')
89-
90-
# de-gradient the fields ("bias/illumination artifact")
91-
n4_correct = pe.Node(ants.N4BiasFieldCorrection(dimension=3, copy_header=True),
92-
name='n4_correct', n_procs=omp_nthreads)
93-
bet = pe.Node(BETRPT(generate_report=True, frac=0.6, mask=True),
94-
name='bet')
95-
ds_report_fmap_mask = pe.Node(DerivativesDataSink(
96-
desc='brain', suffix='mask'), name='ds_report_fmap_mask',
97-
run_without_submitting=True)
89+
applymsk = pe.Node(fsl.ApplyMask(), name='applymsk')
90+
fmap_postproc_wf = init_fmap_postproc_wf(omp_nthreads=omp_nthreads,
91+
fmap_bspline=fmap_bspline)
9892

9993
workflow.connect([
100-
(inputnode, magmrg, [('magnitude', 'in_files')]),
10194
(inputnode, fmapmrg, [('fieldmap', 'in_files')]),
102-
(magmrg, n4_correct, [('out_file', 'input_image')]),
103-
(n4_correct, bet, [('output_image', 'in_file')]),
104-
(bet, outputnode, [('mask_file', 'fmap_mask'),
105-
('out_file', 'fmap_ref')]),
106-
(inputnode, ds_report_fmap_mask, [('fieldmap', 'source_file')]),
107-
(bet, ds_report_fmap_mask, [('out_report', 'in_file')]),
95+
(fmapmrg, applymsk, [('out_file', 'in_file')]),
96+
(magnitude_wf, applymsk, [('outputnode.fmap_mask', 'mask_file')]),
97+
(applymsk, fmap_postproc_wf, [('out_file', 'inputnode.fmap')]),
98+
(magnitude_wf, fmap_postproc_wf, [
99+
('outputnode.fmap_mask', 'inputnode.fmap_mask'),
100+
('outputnode.fmap_ref', 'inputnode.fmap_ref')]),
101+
(fmap_postproc_wf, outputnode, [('outputnode.out_fmap', 'fmap')]),
108102
])
109-
110-
if fmap_bspline:
111-
# despike_threshold=1.0, mask_erode=1),
112-
fmapenh = pe.Node(FieldEnhance(unwrap=False, despike=False),
113-
name='fmapenh', mem_gb=4, n_procs=omp_nthreads)
114-
115-
workflow.connect([
116-
(bet, fmapenh, [('mask_file', 'in_mask'),
117-
('out_file', 'in_magnitude')]),
118-
(fmapmrg, fmapenh, [('out_file', 'in_file')]),
119-
(fmapenh, outputnode, [('out_file', 'fmap')]),
120-
])
121-
122-
else:
123-
torads = pe.Node(FieldToRadS(), name='torads')
124-
prelude = pe.Node(fsl.PRELUDE(), name='prelude')
125-
tohz = pe.Node(FieldToHz(), name='tohz')
126-
127-
denoise = pe.Node(fsl.SpatialFilter(operation='median', kernel_shape='sphere',
128-
kernel_size=3), name='denoise')
129-
demean = pe.Node(niu.Function(function=demean_image), name='demean')
130-
cleanup_wf = cleanup_edge_pipeline(name='cleanup_wf')
131-
132-
applymsk = pe.Node(fsl.ApplyMask(), name='applymsk')
133-
134-
workflow.connect([
135-
(bet, prelude, [('mask_file', 'mask_file'),
136-
('out_file', 'magnitude_file')]),
137-
(fmapmrg, torads, [('out_file', 'in_file')]),
138-
(torads, tohz, [('fmap_range', 'range_hz')]),
139-
(torads, prelude, [('out_file', 'phase_file')]),
140-
(prelude, tohz, [('unwrapped_phase_file', 'in_file')]),
141-
(tohz, denoise, [('out_file', 'in_file')]),
142-
(denoise, demean, [('out_file', 'in_file')]),
143-
(demean, cleanup_wf, [('out', 'inputnode.in_file')]),
144-
(bet, cleanup_wf, [('mask_file', 'inputnode.in_mask')]),
145-
(cleanup_wf, applymsk, [('outputnode.out_file', 'in_file')]),
146-
(bet, applymsk, [('mask_file', 'mask_file')]),
147-
(applymsk, outputnode, [('out_file', 'fmap')]),
148-
])
149-
150103
return workflow

‎sdcflows/workflows/gre.py

Lines changed: 219 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,219 @@
1+
# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
2+
# vi: set ft=python sts=4 ts=4 sw=4 et:
3+
"""Processing phase-difference (aka :abbr:`GRE (gradient-recalled echo)`) fieldmaps.
4+
5+
.. _gre-fieldmaps:
6+
7+
Workflows for processing :abbr:`GRE (gradient recalled echo)` fieldmaps
8+
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
9+
10+
Workflows for preparing the magnitude part of :abbr:`GRE (gradient-recalled echo)` fieldmap
11+
images and cleaning up the fieldmaps created from the phases or phasediff.
12+
13+
"""
14+
15+
from nipype.pipeline import engine as pe
16+
from nipype.interfaces import utility as niu, fsl, ants
17+
from niflow.nipype1.workflows.dmri.fsl.utils import cleanup_edge_pipeline
18+
from niworkflows.engine.workflows import LiterateWorkflow as Workflow
19+
from niworkflows.interfaces.images import IntraModalMerge
20+
from niworkflows.interfaces.masks import BETRPT
21+
22+
23+
def init_magnitude_wf(omp_nthreads, name='magnitude_wf'):
24+
"""
25+
Prepare the magnitude part of :abbr:`GRE (gradient-recalled echo)` fieldmaps.
26+
27+
Average (if not done already) the magnitude part of the
28+
:abbr:`GRE (gradient recalled echo)` images, run N4 to
29+
correct for B1 field nonuniformity, and skull-strip the
30+
preprocessed magnitude.
31+
32+
Workflow Graph
33+
.. workflow ::
34+
:graph2use: orig
35+
:simple_form: yes
36+
37+
from sdcflows.workflows.fmap import init_magnitude_wf
38+
wf = init_magnitude_wf(omp_nthreads=6)
39+
40+
Parameters
41+
----------
42+
omp_nthreads : int
43+
Maximum number of threads an individual process may use
44+
name : str
45+
Name of workflow (default: ``prepare_magnitude_w``)
46+
47+
Inputs
48+
------
49+
magnitude : pathlike
50+
Path to the corresponding magnitude path(s).
51+
52+
Outputs
53+
-------
54+
fmap_ref : pathlike
55+
Path to the fieldmap reference calculated in this workflow.
56+
fmap_mask : pathlike
57+
Path to a binary brain mask corresponding to the reference above.
58+
59+
"""
60+
workflow = Workflow(name=name)
61+
inputnode = pe.Node(
62+
niu.IdentityInterface(fields=['magnitude']), name='inputnode')
63+
outputnode = pe.Node(
64+
niu.IdentityInterface(fields=['fmap_ref', 'fmap_mask', 'mask_report']),
65+
name='outputnode')
66+
67+
# Merge input magnitude images
68+
magmrg = pe.Node(IntraModalMerge(), name='magmrg')
69+
70+
# de-gradient the fields ("bias/illumination artifact")
71+
n4_correct = pe.Node(ants.N4BiasFieldCorrection(dimension=3, copy_header=True),
72+
name='n4_correct', n_procs=omp_nthreads)
73+
bet = pe.Node(BETRPT(generate_report=True, frac=0.6, mask=True),
74+
name='bet')
75+
76+
workflow.connect([
77+
(inputnode, magmrg, [('magnitude', 'in_files')]),
78+
(magmrg, n4_correct, [('out_file', 'input_image')]),
79+
(n4_correct, bet, [('output_image', 'in_file')]),
80+
(bet, outputnode, [('mask_file', 'fmap_mask'),
81+
('out_file', 'fmap_ref'),
82+
('out_report', 'mask_report')]),
83+
])
84+
return workflow
85+
86+
87+
def init_fmap_postproc_wf(omp_nthreads, fmap_bspline, median_kernel_size=3,
88+
name='fmap_postproc_wf'):
89+
"""
90+
Postprocess a B0 map estimated elsewhere.
91+
92+
This workflow denoises (mostly via smoothing) a B0 fieldmap.
93+
94+
Workflow Graph
95+
.. workflow ::
96+
:graph2use: orig
97+
:simple_form: yes
98+
99+
from sdcflows.workflows.fmap import init_fmap_postproc_wf
100+
wf = init_fmap_postproc_wf(omp_nthreads=6, fmap_bspline=False)
101+
102+
Parameters
103+
----------
104+
omp_nthreads : int
105+
Maximum number of threads an individual process may use
106+
fmap_bspline : bool
107+
Whether the fieldmap should be smoothed and extrapolated to off-brain regions
108+
using B-Spline basis.
109+
median_kernel_size : int
110+
Size of the kernel when smoothing is done with a median filter.
111+
name : str
112+
Name of workflow (default: ``fmap_postproc_wf``)
113+
114+
Inputs
115+
------
116+
fmap_mask : pathlike
117+
A brain binary mask corresponding to this fieldmap.
118+
fmap_ref : pathlike
119+
A preprocessed magnitude/reference image for the fieldmap.
120+
fmap : pathlike
121+
A B0-field nonuniformity map (aka fieldmap) estimated elsewhere.
122+
123+
Outputs
124+
-------
125+
out_fmap : pathlike
126+
Postprocessed fieldmap.
127+
128+
"""
129+
workflow = Workflow(name=name)
130+
inputnode = pe.Node(niu.IdentityInterface(
131+
fields=['fmap_mask', 'fmap_ref', 'fmap']), name='inputnode')
132+
outputnode = pe.Node(niu.IdentityInterface(fields=['out_fmap']),
133+
name='outputnode')
134+
if fmap_bspline:
135+
from ..interfaces.fmap import FieldEnhance
136+
# despike_threshold=1.0, mask_erode=1),
137+
fmapenh = pe.Node(
138+
FieldEnhance(unwrap=False, despike=False),
139+
name='fmapenh', mem_gb=4, n_procs=omp_nthreads)
140+
141+
workflow.connect([
142+
(inputnode, fmapenh, [('fmap_mask', 'in_mask'),
143+
('fmap_ref', 'in_magnitude'),
144+
('fmap_hz', 'in_file')]),
145+
(fmapenh, outputnode, [('out_file', 'out_fmap')]),
146+
])
147+
148+
else:
149+
recenter = pe.Node(niu.Function(function=_recenter),
150+
name='recenter', run_without_submitting=True)
151+
denoise = pe.Node(fsl.SpatialFilter(
152+
operation='median', kernel_shape='sphere',
153+
kernel_size=median_kernel_size), name='denoise')
154+
demean = pe.Node(niu.Function(function=_demean), name='demean')
155+
cleanup_wf = cleanup_edge_pipeline(name="cleanup_wf")
156+
157+
workflow.connect([
158+
(inputnode, cleanup_wf, [('fmap_mask', 'inputnode.in_mask')]),
159+
(inputnode, recenter, [('fmap', 'in_file')]),
160+
(recenter, denoise, [('out', 'in_file')]),
161+
(denoise, demean, [('out_file', 'in_file')]),
162+
(demean, cleanup_wf, [('out', 'inputnode.in_file')]),
163+
(cleanup_wf, outputnode, [('outputnode.out_file', 'out_fmap')]),
164+
])
165+
166+
return workflow
167+
168+
169+
def _recenter(in_file):
170+
"""Recenter the phase-map distribution to the -pi..pi range."""
171+
from os import getcwd
172+
import numpy as np
173+
import nibabel as nb
174+
from nipype.utils.filemanip import fname_presuffix
175+
176+
nii = nb.load(in_file)
177+
data = nii.get_fdata(dtype='float32')
178+
msk = data != 0
179+
data[msk] -= np.median(data[msk])
180+
181+
out_file = fname_presuffix(in_file, suffix='_recentered',
182+
newpath=getcwd())
183+
nb.Nifti1Image(data, nii.affine, nii.header).to_filename(out_file)
184+
return out_file
185+
186+
187+
def _demean(in_file, in_mask=None, usemode=True):
188+
"""
189+
Subtract the median (since it is robuster than the mean) from a map.
190+
191+
Parameters
192+
----------
193+
usemode : bool
194+
Use the mode instead of the median (should be even more robust
195+
against outliers).
196+
197+
"""
198+
from os import getcwd
199+
import numpy as np
200+
import nibabel as nb
201+
from nipype.utils.filemanip import fname_presuffix
202+
203+
nii = nb.load(in_file)
204+
data = nii.get_fdata(dtype='float32')
205+
206+
msk = np.ones_like(data, dtype=bool)
207+
if in_mask is not None:
208+
msk[nb.load(in_mask).get_fdata(dtype='float32') < 1e-4] = False
209+
210+
if usemode:
211+
from scipy.stats import mode
212+
data[msk] -= mode(data[msk], axis=None)[0][0]
213+
else:
214+
data[msk] -= np.median(data[msk], axis=None)
215+
216+
out_file = fname_presuffix(in_file, suffix='_demean',
217+
newpath=getcwd())
218+
nb.Nifti1Image(data, nii.affine, nii.header).to_filename(out_file)
219+
return out_file

0 commit comments

Comments
 (0)