-
Notifications
You must be signed in to change notification settings - Fork 15
/
phdiff.py
142 lines (114 loc) · 5.97 KB
/
phdiff.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
# vi: set ft=python sts=4 ts=4 sw=4 et:
"""
Phase-difference B0 estimation.
.. _sdc_phasediff :
Phase-difference B0 estimation
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
The field inhomogeneity inside the scanner (fieldmap) is proportional to the
phase drift between two subsequent :abbr:`GRE (gradient recall echo)`
sequence.
This corresponds to `this section of the BIDS specification
<https://bids-specification.readthedocs.io/en/stable/04-modality-specific-files/01-magnetic-resonance-imaging-data.html#two-phase-images-and-two-magnitude-images>`__.
"""
from nipype.interfaces import fsl, utility as niu
from nipype.pipeline import engine as pe
from niworkflows.engine.workflows import LiterateWorkflow as Workflow
from ..interfaces.fmap import Phasediff2Fieldmap, PhaseMap2rads, SubtractPhases
from .gre import init_fmap_postproc_wf, init_magnitude_wf
def init_phdiff_wf(omp_nthreads, name='phdiff_wf'):
r"""
Distortion correction of EPI sequences using phase-difference maps.
Estimates the fieldmap using a phase-difference image and one or more
magnitude images corresponding to two or more :abbr:`GRE (Gradient Echo sequence)`
acquisitions.
The most delicate bit of this workflow is the phase-unwrapping process: phase maps
are clipped in the range :math:`[0 \dotsb 2\pi )`.
To find the integer number of offsets that make a region continously smooth with
its neighbour, FSL PRELUDE is run [Jenkinson2003]_.
FSL PRELUDE takes wrapped maps in the range 0 to 6.28, `as per the user guide
<https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/FUGUE/Guide#Step_2_-_Getting_.28wrapped.29_phase_in_radians>`__.
For the phase-difference maps, recentering back to :math:`[-\pi \dotsb \pi )` is necessary.
After some massaging and the application of the effective echo spacing parameter,
the phase-difference maps can be converted into a *B0 field map* in Hz units.
The `original code was taken from nipype
<https://github.com/nipy/nipype/blob/0.12.1/nipype/workflows/dmri/fsl/artifacts.py#L514>`_.
Workflow Graph
.. workflow ::
:graph2use: orig
:simple_form: yes
from aslprep.sdcflows.workflows.phdiff import init_phdiff_wf
wf = init_phdiff_wf(omp_nthreads=1)
Parameters
----------
omp_nthreads : int
Maximum number of threads an individual process may use
Inputs
------
magnitude : list of os.pathlike
List of path(s) the GRE magnitude maps.
phasediff : list of tuple(os.pathlike, dict)
List containing one GRE phase-difference map with its corresponding metadata
(requires ``EchoTime1`` and ``EchoTime2``), or the phase maps for the two
subsequent echoes, with their metadata (requires ``EchoTime``).
Outputs
-------
fmap_ref : pathlike
The average magnitude image, skull-stripped
fmap_mask : pathlike
The brain mask applied to the fieldmap
fmap : pathlike
The estimated fieldmap in Hz
References
----------
.. [Jenkinson2003] Jenkinson, M. (2003) Fast, automated, N-dimensional phase-unwrapping
algorithm. MRM 49(1):193-197. doi:`10.1002/mrm.10354 <10.1002/mrm.10354>`__.
"""
workflow = Workflow(name=name)
workflow.__desc__ = """\
A B0-nonuniformity map (or *fieldmap*) was estimated based on a phase-difference map
calculated with a dual-echo GRE (gradient-recall echo) sequence, processed with a
custom workflow of *SDCFlows* inspired by the
[`epidewarp.fsl` script](http://www.nmr.mgh.harvard.edu/~greve/fbirn/b0/epidewarp.fsl)
and further improvements in HCP Pipelines [@hcppipelines].
"""
inputnode = pe.Node(niu.IdentityInterface(fields=['magnitude', 'phasediff']),
name='inputnode')
outputnode = pe.Node(niu.IdentityInterface(
fields=['fmap', 'fmap_ref', 'fmap_mask']), name='outputnode')
split = pe.MapNode(niu.Function(function=_split, output_names=['map_file', 'meta']),
iterfield=['phasediff'], run_without_submitting=True, name='split')
magnitude_wf = init_magnitude_wf(omp_nthreads=omp_nthreads)
# phase diff -> radians
phmap2rads = pe.MapNode(PhaseMap2rads(), name='phmap2rads',
iterfield=['in_file'], run_without_submitting=True)
# FSL PRELUDE will perform phase-unwrapping
prelude = pe.Node(fsl.PRELUDE(), name='prelude')
calc_phdiff = pe.Node(SubtractPhases(), name='calc_phdiff',
run_without_submitting=True)
fmap_postproc_wf = init_fmap_postproc_wf(omp_nthreads=omp_nthreads,
fmap_bspline=False)
compfmap = pe.Node(Phasediff2Fieldmap(), name='compfmap')
workflow.connect([
(inputnode, split, [('phasediff', 'phasediff')]),
(inputnode, magnitude_wf, [('magnitude', 'inputnode.magnitude')]),
(magnitude_wf, prelude, [('outputnode.fmap_ref', 'magnitude_file'),
('outputnode.fmap_mask', 'mask_file')]),
(split, phmap2rads, [('map_file', 'in_file')]),
(phmap2rads, calc_phdiff, [('out_file', 'in_phases')]),
(split, calc_phdiff, [('meta', 'in_meta')]),
(calc_phdiff, prelude, [('phase_diff', 'phase_file')]),
(prelude, fmap_postproc_wf, [('unwrapped_phase_file', 'inputnode.fmap')]),
(calc_phdiff, fmap_postproc_wf, [('metadata', 'inputnode.metadata')]),
(magnitude_wf, fmap_postproc_wf, [
('outputnode.fmap_mask', 'inputnode.fmap_mask'),
('outputnode.fmap_ref', 'inputnode.fmap_ref')]),
(fmap_postproc_wf, compfmap, [('outputnode.out_fmap', 'in_file'),
('outputnode.metadata', 'metadata')]),
(compfmap, outputnode, [('out_file', 'fmap')]),
(magnitude_wf, outputnode, [('outputnode.fmap_ref', 'fmap_ref'),
('outputnode.fmap_mask', 'fmap_mask')]),
])
return workflow
def _split(phasediff):
return phasediff