Skip to content

Commit 0ad0736

Browse files
oestebanoesteban
committed
enh:new all-corrections workflows and example
When merged the new dipy interfaces I've written, this example should be updated adding a denoising node.
1 parent b9f1c25 commit 0ad0736

File tree

4 files changed

+364
-17
lines changed

4 files changed

+364
-17
lines changed

‎examples/dmri_preprocessing.py

Lines changed: 147 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,147 @@
1+
#!/usr/bin/env python
2+
# -*- coding: utf-8 -*-
3+
# @Author: oesteban
4+
# @Date: 2014-08-31 20:32:22
5+
# @Last Modified by: oesteban
6+
# @Last Modified time: 2014-09-01 18:12:14
7+
"""
8+
===================
9+
dMRI: Preprocessing
10+
===================
11+
12+
Introduction
13+
============
14+
15+
This script, dmri_preprocessing.py, demonstrates how to prepare dMRI data
16+
for tractography and connectivity analysis with nipype.
17+
18+
We perform this analysis using the FSL course data, which can be acquired from here:
19+
http://www.fmrib.ox.ac.uk/fslcourse/fsl_course_data2.tar.gz
20+
21+
Can be executed in command line using ``python dmri_preprocessing.py``
22+
23+
24+
Import necessary modules from nipype.
25+
"""
26+
27+
import nipype.interfaces.io as nio # Data i/o
28+
import nipype.interfaces.utility as util # utility
29+
import nipype.pipeline.engine as pe # pypeline engine
30+
import nipype.interfaces.fsl as fsl
31+
import nipype.algorithms.misc as misc
32+
import os # system functions
33+
34+
35+
"""
36+
Load specific nipype's workflows for preprocessing of dMRI data:
37+
:class:`nipype.workflows.dmri.preprocess.epi.all_peb_pipeline`,
38+
as data include a *b0* volume with reverse encoding direction
39+
(*P>>>A*, or *y*), in contrast with the general acquisition encoding
40+
that is *A>>>P* or *-y* (in RAS systems).
41+
"""
42+
43+
from nipype.workflows.dmri.preprocess.epi import all_peb_pipeline
44+
45+
46+
"""
47+
Map field names into individual subject runs
48+
"""
49+
50+
info = dict(dwi=[['subject_id', 'dwidata']],
51+
bvecs=[['subject_id','bvecs']],
52+
bvals=[['subject_id','bvals']],
53+
dwi_rev=[['subject_id','nodif_PA']])
54+
55+
infosource = pe.Node(interface=util.IdentityInterface(fields=['subject_id']),
56+
name="infosource")
57+
58+
# Set the subject 1 identifier in subject_list,
59+
# we choose the preproc dataset as it contains uncorrected files.
60+
subject_list = ['subj1_preproc']
61+
62+
63+
"""Here we set up iteration over all the subjects. The following line
64+
is a particular example of the flexibility of the system. The
65+
``datasource`` attribute ``iterables`` tells the pipeline engine that
66+
it should repeat the analysis on each of the items in the
67+
``subject_list``. In the current example, the entire first level
68+
preprocessing and estimation will be repeated for each subject
69+
contained in subject_list.
70+
"""
71+
72+
infosource.iterables = ('subject_id', subject_list)
73+
74+
75+
"""
76+
Now we create a :class:`nipype.interfaces.io.DataGrabber` object and
77+
fill in the information from above about the layout of our data. The
78+
:class:`~nipype.pipeline.engine.Node` module wraps the interface object
79+
and provides additional housekeeping and pipeline specific
80+
functionality.
81+
"""
82+
83+
datasource = pe.Node(nio.DataGrabber(infields=['subject_id'],
84+
outfields=info.keys()), name = 'datasource')
85+
86+
datasource.inputs.template = "%s/%s"
87+
88+
# This needs to point to the fdt folder you can find after extracting
89+
# http://www.fmrib.ox.ac.uk/fslcourse/fsl_course_data2.tar.gz
90+
datasource.inputs.base_directory = os.path.abspath('fdt1')
91+
datasource.inputs.field_template = dict(dwi='%s/%s.nii.gz', dwi_rev='%s/%s.nii.gz')
92+
datasource.inputs.template_args = info
93+
datasource.inputs.sort_filelist = True
94+
95+
96+
"""
97+
An inputnode is used to pass the data obtained by the data grabber to the actual processing functions
98+
"""
99+
100+
inputnode = pe.Node(util.IdentityInterface(fields=["dwi", "bvecs", "bvals",
101+
"dwi_rev"]), name="inputnode")
102+
103+
104+
"""
105+
106+
Setup for dMRI preprocessing
107+
============================
108+
109+
In this section we initialize the appropriate workflow for preprocessing of diffusion images.
110+
Particularly, we look into the ``acqparams.txt`` file of the selected subject to gather the
111+
encoding direction, acceleration factor (in parallel sequences it is > 1), and readout time or
112+
echospacing.
113+
114+
"""
115+
116+
epi_AP = {'echospacing': 66.5e-3, 'acc_factor': 1, 'enc_dir': 'y-'}
117+
epi_PA = {'echospacing': 66.5e-3, 'acc_factor': 1, 'enc_dir': 'y'}
118+
prep = all_peb_pipeline(epi_params=epi_AP, altepi_params=epi_PA)
119+
120+
121+
"""
122+
Connect nodes in workflow
123+
=========================
124+
125+
We create a higher level workflow to connect the nodes. Please excuse the author for
126+
writing the arguments of the ``connect`` function in a not-standard fashion with readability
127+
aims.
128+
"""
129+
130+
wf = pe.Workflow(name="dMRI_Preprocessing")
131+
wf.base_dir = os.path.abspath('preprocessing_dmri_tutorial')
132+
wf.connect([
133+
(infosource, datasource, [('subject_id', 'subject_id')])
134+
,(datasource, prep, [('dwi', 'inputnode.in_file'),
135+
('dwi_rev', 'inputnode.alt_file'),
136+
('bvals', 'inputnode.in_bval'),
137+
('bvecs', 'inputnode.in_bvec')])
138+
])
139+
140+
141+
"""
142+
Run the workflow as command line executable
143+
"""
144+
145+
if __name__ == '__main__':
146+
wf.run()
147+
wf.write_graph()

‎nipype/workflows/dmri/preprocess/__init__.py

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1 +1,2 @@
1-
from epi import hmc_pipeline, ecc_pipeline, sdc_fmb
1+
from epi import (all_fmb_pipeline, all_peb_pipeline,
2+
hmc_pipeline, ecc_pipeline, sdc_fmb, sdc_peb)

‎nipype/workflows/dmri/preprocess/epi.py

Lines changed: 200 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -5,11 +5,126 @@
55
import nipype.pipeline.engine as pe
66
import nipype.interfaces.utility as niu
77
import nipype.interfaces.fsl as fsl
8+
import nipype.interfaces.freesurfer as fs
89
import nipype.interfaces.ants as ants
910
import os
1011

1112
from .utils import *
1213

14+
15+
def all_fmb_pipeline(name='hmc_sdc_ecc'):
16+
"""
17+
Builds a pipeline including three artifact corrections: head-motion correction (HMC),
18+
susceptibility-derived distortion correction (SDC), and Eddy currents-derived distortion
19+
correction (ECC).
20+
"""
21+
inputnode = pe.Node(niu.IdentityInterface(fields=['in_file', 'in_bvec', 'in_bval',
22+
'bmap_pha', 'bmap_mag']), name='inputnode')
23+
24+
outputnode = pe.Node(niu.IdentityInterface(fields=['out_file', 'out_mask', 'out_bvec']),
25+
name='outputnode')
26+
27+
28+
avg_b0_0 = pe.Node(niu.Function(input_names=['in_dwi', 'in_bval'],
29+
output_names=['out_file'], function=b0_average), name='b0_avg_pre')
30+
avg_b0_1 = pe.Node(niu.Function(input_names=['in_dwi', 'in_bval'],
31+
output_names=['out_file'], function=b0_average), name='b0_avg_post')
32+
bet_dwi0 = pe.Node(fsl.BET(frac=0.3, mask=True, robust=True), name='bet_dwi_pre')
33+
bet_dwi1 = pe.Node(fsl.BET(frac=0.3, mask=True, robust=True), name='bet_dwi_post')
34+
35+
hmc = hmc_pipeline()
36+
sdc = sdc_fmb()
37+
ecc = ecc_pipeline()
38+
39+
regrid = pe.Node(fs.MRIConvert(vox_size=(2.0, 2.0, 2.0), out_orientation='RAS'),
40+
name='Reslice')
41+
42+
wf = pe.Workflow('dMRI_Artifacts')
43+
wf.connect([
44+
(inputnode, hmc, [('in_file', 'inputnode.in_file'),
45+
('in_bvec', 'inputnode.in_bvec')])
46+
,(inputnode, avg_b0_0, [('in_file', 'in_dwi'),
47+
('in_bval', 'in_bval')])
48+
,(avg_b0_0, bet_dwi0, [('out_file','in_file')])
49+
,(bet_dwi0, hmc, [('mask_file', 'inputnode.in_mask')])
50+
51+
,(hmc, sdc, [('outputnode.out_file', 'inputnode.in_file')])
52+
,(bet_dwi0, sdc, [('mask_file', 'inputnode.in_mask')])
53+
,(inputnode, sdc, [('in_bval', 'inputnode.in_bval'),
54+
('bmap_pha', 'inputnode.bmap_pha'),
55+
('bmap_mag', 'inputnode.bmap_mag')])
56+
,(inputnode, ecc, [('in_bval', 'inputnode.in_bval')])
57+
,(bet_dwi0, ecc, [('mask_file', 'inputnode.in_mask')])
58+
,(sdc, ecc, [('outputnode.out_file', 'inputnode.in_file')])
59+
,(hmc, outputnode, [('outputnode.out_bvec', 'out_bvec')])
60+
,(ecc, regrid, [('outputnode.out_file', 'in_file')])
61+
,(regrid, outputnode, [('out_file', 'out_file')])
62+
,(regrid, avg_b0_1, [('out_file', 'in_dwi')])
63+
,(inputnode, avg_b0_1, [('in_bval', 'in_bval')])
64+
,(avg_b0_1, bet_dwi1, [('out_file','in_file')])
65+
,(bet_dwi1, outputnode, [('mask_file', 'out_mask')])
66+
])
67+
return wf
68+
69+
70+
def all_peb_pipeline(name='hmc_sdc_ecc',
71+
epi_params=dict(echospacing=0.77e-3,
72+
acc_factor=3,
73+
enc_dir='y-'),
74+
altepi_params=dict(echospacing=0.77e-3,
75+
acc_factor=3,
76+
enc_dir='y')):
77+
"""
78+
Builds a pipeline including three artifact corrections: head-motion correction (HMC),
79+
susceptibility-derived distortion correction (SDC), and Eddy currents-derived distortion
80+
correction (ECC).
81+
"""
82+
inputnode = pe.Node(niu.IdentityInterface(fields=['in_file', 'in_bvec', 'in_bval',
83+
'alt_file']), name='inputnode')
84+
85+
outputnode = pe.Node(niu.IdentityInterface(fields=['out_file', 'out_mask', 'out_bvec']),
86+
name='outputnode')
87+
88+
avg_b0_0 = pe.Node(niu.Function(input_names=['in_dwi', 'in_bval'],
89+
output_names=['out_file'], function=b0_average), name='b0_avg_pre')
90+
avg_b0_1 = pe.Node(niu.Function(input_names=['in_dwi', 'in_bval'],
91+
output_names=['out_file'], function=b0_average), name='b0_avg_post')
92+
bet_dwi0 = pe.Node(fsl.BET(frac=0.3, mask=True, robust=True), name='bet_dwi_pre')
93+
bet_dwi1 = pe.Node(fsl.BET(frac=0.3, mask=True, robust=True), name='bet_dwi_post')
94+
95+
hmc = hmc_pipeline()
96+
sdc = sdc_peb(epi_params=epi_params, altepi_params=altepi_params)
97+
ecc = ecc_pipeline()
98+
99+
regrid = pe.Node(fs.MRIConvert(vox_size=(2.0, 2.0, 2.0), out_orientation='RAS'),
100+
name='Reslice')
101+
102+
wf = pe.Workflow('dMRI_Artifacts')
103+
wf.connect([
104+
(inputnode, hmc, [('in_file', 'inputnode.in_file'),
105+
('in_bvec', 'inputnode.in_bvec')])
106+
,(inputnode, avg_b0_0, [('in_file', 'in_dwi'),
107+
('in_bval', 'in_bval')])
108+
,(avg_b0_0, bet_dwi0, [('out_file','in_file')])
109+
,(bet_dwi0, hmc, [('mask_file', 'inputnode.in_mask')])
110+
,(hmc, sdc, [('outputnode.out_file', 'inputnode.in_file')])
111+
,(bet_dwi0, sdc, [('mask_file', 'inputnode.in_mask')])
112+
,(inputnode, sdc, [('in_bval', 'inputnode.in_bval'),
113+
('alt_file', 'inputnode.alt_file')])
114+
,(inputnode, ecc, [('in_bval', 'inputnode.in_bval')])
115+
,(bet_dwi0, ecc, [('mask_file', 'inputnode.in_mask')])
116+
,(sdc, ecc, [('outputnode.out_file', 'inputnode.in_file')])
117+
,(hmc, outputnode, [('outputnode.out_bvec', 'out_bvec')])
118+
,(ecc, regrid, [('outputnode.out_file', 'in_file')])
119+
,(regrid, outputnode, [('out_file', 'out_file')])
120+
,(regrid, avg_b0_1, [('out_file', 'in_dwi')])
121+
,(inputnode, avg_b0_1, [('in_bval', 'in_bval')])
122+
,(avg_b0_1, bet_dwi1, [('out_file','in_file')])
123+
,(bet_dwi1, outputnode, [('mask_file', 'out_mask')])
124+
])
125+
return wf
126+
127+
13128
def hmc_pipeline(name='motion_correct'):
14129
"""
15130
HMC stands for head-motion correction.
@@ -116,7 +231,7 @@ def ecc_pipeline(name='eddy_correct'):
116231
[Rohde04]_.
117232
118233
A list of rigid transformation matrices can be provided, sourcing from a
119-
:func:`.motion_correct` workflow, to initialize registrations in a *motion free*
234+
:func:`.hmc_pipeline` workflow, to initialize registrations in a *motion free*
120235
framework.
121236
122237
A list of affine transformation matrices is available as output, so that transforms
@@ -344,6 +459,90 @@ def sdc_fmb(name='fmb_correction',
344459
return wf
345460

346461

462+
def sdc_peb(name='peb_correction',
463+
epi_params=dict(echospacing=0.77e-3,
464+
acc_factor=3,
465+
enc_dir='y-'),
466+
altepi_params=dict(echospacing=0.77e-3,
467+
acc_factor=3,
468+
enc_dir='y')):
469+
"""
470+
SDC stands for susceptibility distortion correction. PEB stands for phase-encoding-based.
471+
472+
The phase-encoding-based (PEB) method implements SDC by acquiring diffusion images with
473+
two different enconding directions [Andersson2003]_. The most typical case is acquiring
474+
with opposed phase-gradient blips (e.g. *A>>>P* and *P>>>A*, or equivalently, *-y* and *y*)
475+
as in [Chiou2000]_, but it is also possible to use orthogonal configurations [Cordes2000]_
476+
(e.g. *A>>>P* and *L>>>R*, or equivalently *-y* and *x*).
477+
This workflow uses the implementation of FSL
478+
(`TOPUP <http://fsl.fmrib.ox.ac.uk/fsl/fslwiki/TOPUP>`_).
479+
480+
Example
481+
-------
482+
483+
>>> from nipype.workflows.dmri.preprocess.epi import sdc_peb
484+
>>> peb = sdc_peb()
485+
>>> peb.inputs.inputnode.in_file = 'diffusion.nii'
486+
>>> peb.inputs.inputnode.in_bval = 'diffusion.bval'
487+
>>> peb.inputs.inputnode.in_mask = 'mask.nii'
488+
>>> peb.run() # doctest: +SKIP
489+
490+
.. admonition:: References
491+
492+
.. [Andersson2003] Andersson JL et al., `How to correct susceptibility distortions in
493+
spin-echo echo-planar images: application to diffusion tensor imaging
494+
<http://dx.doi.org/10.1016/S1053-8119(03)00336-7>`_.
495+
Neuroimage. 2003 Oct;20(2):870-88. doi: 10.1016/S1053-8119(03)00336-7
496+
497+
.. [Cordes2000] Cordes D et al., Geometric distortion correction in EPI using two
498+
images with orthogonal phase-encoding directions, in Proc. ISMRM (8), p.1712,
499+
Denver, US, 2000.
500+
501+
.. [Chiou2000] Chiou JY, and Nalcioglu O, A simple method to correct off-resonance
502+
related distortion in echo planar imaging, in Proc. ISMRM (8), p.1712, Denver, US,
503+
2000.
504+
505+
"""
506+
507+
inputnode = pe.Node(niu.IdentityInterface(fields=['in_file', 'in_bval', 'in_mask',
508+
'alt_file', 'ref_num']),
509+
name='inputnode')
510+
outputnode = pe.Node(niu.IdentityInterface(fields=['out_file', 'out_vsm']),
511+
name='outputnode')
512+
513+
b0_ref = pe.Node(fsl.ExtractROI(t_size=1), name='b0_ref')
514+
b0_alt = pe.Node(fsl.ExtractROI(t_size=1), name='b0_alt')
515+
b0_comb = pe.Node(niu.Merge(2), name='b0_list')
516+
b0_merge = pe.Node(fsl.Merge(dimension='t'), name='b0_merged')
517+
518+
topup = pe.Node(fsl.TOPUP(), name='topup')
519+
topup.inputs.encoding_direction = [epi_params['enc_dir'], altepi_params['enc_dir']]
520+
topup.inputs.readout_times = [epi_params['echospacing']/(1.0*epi_params['acc_factor']),
521+
altepi_params['echospacing']/(1.0*altepi_params['acc_factor'])]
522+
unwarp = pe.Node(fsl.ApplyTOPUP(in_index=[1,2]), name='unwarp')
523+
dwi_comb = pe.Node(niu.Merge(2), name='DWIcomb')
524+
525+
wf = pe.Workflow(name=name)
526+
wf.connect([
527+
(inputnode, b0_ref, [('in_file','in_file'),
528+
(('ref_num', _checkrnum),'t_min')])
529+
,(inputnode, b0_alt, [('alt_file','in_file'),
530+
(('ref_num', _checkrnum),'t_min')])
531+
,(inputnode, dwi_comb, [('in_file','in1'),
532+
('alt_file','in2') ] )
533+
,(b0_ref, b0_comb, [('roi_file','in1')] )
534+
,(b0_alt, b0_comb, [('roi_file','in2')] )
535+
,(b0_comb, b0_merge, [('out', 'in_files')] )
536+
,(b0_merge, topup, [('merged_file','in_file')])
537+
,(topup, unwarp, [('out_fieldcoef','in_topup_fieldcoef'),
538+
('out_movpar','in_topup_movpar'),
539+
('out_enc_file','encoding_file')])
540+
,(dwi_comb, unwarp, [('out','in_files')])
541+
,(unwarp, outputnode, [('out_corrected','out_file')])
542+
])
543+
return wf
544+
545+
347546
def _checkrnum(ref_num):
348547
from nipype.interfaces.base import isdefined
349548
if (ref_num is None) or not isdefined(ref_num):

0 commit comments

Comments
 (0)