-
Notifications
You must be signed in to change notification settings - Fork 1
/
sliceMotion4d_py2
executable file
·193 lines (164 loc) · 10 KB
/
sliceMotion4d_py2
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
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
#!/usr/bin/env python
#import IPython
import matplotlib
matplotlib.use('Agg') #non-interactive backend for plotting
import numpy as np
import os
import argparse
from nipy import load_image, save_image
from nipy.algorithms.registration import SpaceTimeRealign
from sys import exit
import matplotlib.pyplot as plt
parser = argparse.ArgumentParser(prog='sliceMotion4d')
parser.add_argument('--inputs', '-i', nargs='+', required=True, help='One or more 4d fMRI inputs to realign and interpolate (multiple runs are co-registered to each other)')
parser.add_argument('--tr', '-t', required=True, type=float, help='The repetition time (TR) of the acquisition sequence in seconds.')
parser.add_argument('--loops', '-l', required=False, type=str, default='5,1', help='Number of within-run realignment loops to run. Default is 5,1 indicating 5 loops at coarse subsampling, 1 at fine (see --speedup).')
parser.add_argument('--speedup', required=False, type=str, default='5,2', help='Comma-separated list of spatial subsampling factor for computing affine transforms. Default is 5,2, indicating 5x subsampling (coarse), then 2x subsampling (fine)')
parser.add_argument('--bw_loops', required=False, type=str, default='5,1', help='Number of between-run realignment loops to run (only relevant with multiple --inputs). Default is 5,1 indicating 5 between-run loops at coarse subsampling, 1 at fine (see --speedup).')
parser.add_argument('--slice_times', '-s', nargs='+', required=True, help='Time of acquisition for each slice (starting with bottom slice). Common options: ascending, descending, interleaved. Can also be multiple arguments specifying time of each slice (in seconds): 0 0.25 0.5 0.75 ... Can also be a comma-delimited file containing the exact timings (in seconds): sliceTimes.1D')
parser.add_argument('--slice_direction', '-d', required=False, type=str, default='z', help='Direction of slice acquisition: x, y, or z. Default is z')
parser.add_argument('--prefix', '-p', required=False, type=str, default='mt_', help='File prefix for resampled data. Default is mt_')
parser.add_argument('--mats', '-m', action="store_true", required=False, default=False, help='Whether to output affine transformation matrices.')
parser.add_argument('--dtype', required=False, type=str, default="int16", help='Data type for saved nifti images (int16, float32, float64). Default is int16')
parser.add_argument('--plots', action="store_true", required=False, default=False, help='Create a plot of motion parameters as a .png file.')
parser.add_argument('--siemens', action="store_true", required=False, default=False, help='Whether data were acquired on a Siemens scanner, which has a strange slice order for interleaved sequences depending on number of slices.')
parser.add_argument('--refscan', required=False, type=int, default=None, help='Index of reference scan (first volume = 0) to which all volumes are aligned.')
# for testing
# cmdInput = parser.parse_args('--inputs 1.nii 2.nii 3.nii --tr 1.5 --slice_times 0 0.5 1.0 1.5 2 2.5 3'.split())
cmdInput = parser.parse_args()
slice_times=None
if len(cmdInput.slice_times) == 1:
if os.path.isfile(cmdInput.slice_times[0]):
print "Using custom slice timings from csv file: " + cmdInput.slice_times[0]
slice_times = np.loadtxt(cmdInput.slice_times[0], comments="#", delimiter=",", unpack=False)
print slice_times
elif not cmdInput.slice_times[0].lower() in ['ascending', 'seqasc', 'descending', 'seqdesc', 'interleaved']:
print "--slice_times must be ascending, descending, interleaved, or a numeric array of slice ordering"
exit(1)
else:
if cmdInput.slice_times[0].lower() == 'interleaved':
if cmdInput.siemens:
print "Using Siemens interleaved ascending: 0,2,4,1,3,5 for odd n; 1,3,5,0,2,4 for even n"
slice_times="asc_alt_siemens" #0,2,4,1,3,5 for odd n, 1,3,5,0,2,4 for even; default behavior on Siemens
else:
print "Assuming default interleaved ascending: 0,2,4,...,1,3,5,..."
slice_times="asc_alt_2"
elif cmdInput.slice_times[0].lower() == 'seqasc':
slice_times="ascending" #convert slicetimer nomenclature
elif cmdInput.slice_times[0].lower() == 'seqdesc':
slice_times="descending"
else:
slice_times=cmdInput.slice_times[0].lower()
else:
print "Processing custom slice timings based on command line input."
slice_times=[float(i) for i in cmdInput.slice_times] #convert to list of floats (in seconds)
if not cmdInput.slice_direction.lower() in ['x', 'y', 'z']:
print "--slice_direction must be x, y, or z"
exit(1)
else:
if cmdInput.slice_direction.lower() == "x":
slice_direction=0
elif cmdInput.slice_direction.lower() == "y":
slice_direction=1
elif cmdInput.slice_direction.lower() == "z":
slice_direction=2
wiLoops=map(int, cmdInput.loops.split(',')) #convert input to list
speedup=map(int, cmdInput.speedup.split(','))
bwLoops=map(int, cmdInput.bw_loops.split(','))
#verify comparable lengths for loops and speedup
l1 = len(wiLoops)
if any(len(lst) != l1 for lst in [speedup, bwLoops]):
print "--loops, --speedup, and --bw_loops must all be of the same length"
exit(1)
# process data type
if not cmdInput.dtype.lower() in ('int16', 'float32', 'float64'):
print "--dtype must be one of int16, float32, float64."
exit(1)
else:
if cmdInput.dtype.lower() == "int16":
dtype=np.int16
elif cmdInput.dtype.lower() == "float32":
dtype=np.float32
elif cmdInput.dtype.lower() == "float64":
dtype=np.float64
runs = [ load_image(f) for f in cmdInput.inputs ]
R = SpaceTimeRealign(runs, tr=cmdInput.tr, slice_times=slice_times, slice_info=slice_direction)
# estimate motion and slice timing realignment
# hard code multi-run alignment to use 5 loops
# Jul2014: For now, estimating transforms in two steps: at 5x subsampling (spatially coarse), then at 2x subsampling (spatially fine)
# This is slower, but should improve accuracy of image coregistration. At some point, need to allow for tuple syntax
R.estimate(loops=wiLoops, between_loops=bwLoops, refscan=cmdInput.refscan, speedup=speedup)
# resample data
ra_runs = R.resample()
#save resampled images
for i, corrImage in enumerate(ra_runs):
#save motion estimates, stored in realign object _transforms
motion = R._transforms[i]
#pull apart path and filename
dname = os.path.dirname(cmdInput.inputs[i])
fname = os.path.basename(cmdInput.inputs[i])
#trim off .nii or .nii.gz extension
iname = fname.split('.')
if len(iname) > 2 and (iname[-1].lower() == "gz" and iname[-2].lower() == "nii"):
basename = cmdInput.prefix + '.'.join(iname[:-2]) #drop last two dotted pieces: .nii.gz
elif len(iname) > 1 and iname[-1].lower() in ("nii", "hdr", "img"):
basename = cmdInput.prefix + '.'.join(iname[:-1]) #drop last dotted piece: .nii
else:
print "Can't determine file name of input properly."
exit(1)
#save realigned image
#set data type, e.g., float32 to reduce file size (default from algorithm is float64)
corrImage._data = corrImage._data.astype(dtype)
save_image(corrImage, os.path.join(dname, '%s%s' % (basename, ".nii.gz")))
mname=basename + '.par'
mfile = open(os.path.join(dname, mname), 'w')
#Make mats dir if needed
if cmdInput.mats:
tdir = os.path.join(dname, 'mc_mats')
if not os.path.exists(tdir):
os.makedirs(tdir)
#N.B. nipy rotation parameters are not stored in standard radians with euler angles
#http://mail.scipy.org/pipermail/nipy-devel/2012-October/008396.html
#for consistency with rapidArt, output as is [ translation, rotation ]??
motionparams=np.zeros((len(motion), 6), dtype=np.float32)
for j, mo in enumerate(motion):
params = ['%.10f' % item for item in np.hstack((mo.rotation, mo.translation))]
motionparams[j,:] = params
string = ' '.join(params) + '\n'
mfile.write(string)
if cmdInput.mats:
if len(cmdInput.inputs) > 1:
np.savetxt(os.path.join(tdir, '%s_mot%.4d.mat' % (basename, j)), mo.as_affine(), fmt='%.8f')
else:
np.savetxt(os.path.join(tdir, 'MAT_%.4d' % (j)), mo.as_affine(), fmt='%.8f') #MAT_0000 to match mcflirt
mfile.close()
#IPython.embed()
if cmdInput.plots:
#motion parameter plots
fig01 = plt.figure()
ax01 = fig01.add_subplot(1, 2, 1)
ax01.plot(range(0, len(motion)), motionparams[:,3], label="Tx", linewidth=1.6)
ax01.plot(range(0, len(motion)), motionparams[:,4], label="Ty", linewidth=1.6)
ax01.plot(range(0, len(motion)), motionparams[:,5], label="Tz", linewidth=1.6)
ax01.set_xlabel("Volume number")
ax01.set_ylabel("Translation (mm)")
#ax01.set_title('Translation')
#lgd01 = ax01.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
lgd01=ax01.legend(loc='upper center', bbox_to_anchor=(0.5, 1.145),
ncol=3, fancybox=True, shadow=True)
#rotation
ax02 = fig01.add_subplot(1, 2, 2)
ax02.plot(range(0, len(motion)), motionparams[:,0], label="Rx", linewidth=1.6)
ax02.plot(range(0, len(motion)), motionparams[:,1], label="Ry", linewidth=1.6)
ax02.plot(range(0, len(motion)), motionparams[:,2], label="Rz", linewidth=1.6)
ax02.set_xlabel("Volume number")
ax02.set_ylabel("Rotation (radians)")
#ax02.set_title('Rotation')
#lgd02=ax02.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
lgd02=ax02.legend(loc='upper center', bbox_to_anchor=(0.5, 1.145),
ncol=3, fancybox=True, shadow=True)
fig01.suptitle('Motion estimates for: ' + fname, fontsize=21) #overall title
fig01.set_size_inches(9, 5) #set figure size
fig01.tight_layout() #minimize overlap between subplots
plt.subplots_adjust(top=0.83) #move down panels to make space for overall title
fig01.savefig(os.path.join(dname, basename + '_motion.png'), dpi=300) #, bbox_extra_artists=(lgd01, lgd02,), bbox_inches="tight")