1import sys as sys
2import numpy as np
3
4#np.set_printoptions(linewidth=200)
5
6def read_fsl_eddy_parameters(my_file):
7    '''Read in an FSL (eddy-)produced *.eddy_parameters file and output a
8    simple text file of 6 columns resembling the output of 3dvolreg: 3
9    translations (mm) and 3 rotations (deg).
10
11    Tested/used on FSL v5.0.11; may work on eddy_parameters files from
12    different versions of the same software.
13
14    Note that reordering occurs.  The first 6 columns of FSL/eddy
15    *eddy_parameters files are in the following order:
16
17    rot x
18    rot y
19    rot z
20    trans x
21    trans y
22    trans z
23
24    3dvolreg order is:
25
26    rot z
27    rot x
28    rot y
29    trans z
30    trans x
31    trans y
32
33    '''
34
35    fff = open( my_file ,'r')
36    raw_data = fff.readlines()
37    fff.close()
38
39    data_list1 = []
40    for line in raw_data:
41        aa = np.zeros(6)
42        x = w.split()
43        Nx = len(x)
44        if Nx < 6:
45            sys.exit("** Problem in this line:\n %s\n\n" % line)
46        for i in range (6):
47            aa[i] = float(x[i])
48
49        # convert rads to degs
50        aa[3:]*= 180./np.pi
51
52        # convert to 3dvolreg ordering
53        bb = Conv_TRxyz_to_RTzyx(aa)
54
55        data_list1.append(bb)
56
57    return np.array(data_list1)
58
59def read_tortoise_transformations(my_file):
60    '''Read in a TORTOISE (DIFFPREP-)produced *_transformations.txt file
61    and output a simple text file of 6 columns resembling the output
62    of 3dvolreg: 3 rotations (deg) and 3 translations (mm).
63
64    Tested/used on TORTOISE v3.1; may work on transformations.txt
65    files from different versions of the same software.
66
67    Note that reordering occurs.  The first 6 columns of
68    TORTOISE/DIFFPREP transformation.txt files are in the following
69    order:
70
71    Column 0 : x (for axial data = RL, in mm)
72    Column 1 : y (for axial data = AP, in mm)
73    Column 2 : z (for axial data = IS, in mm)
74    Column 3 : x (Euler angles, in rads)
75    Column 4 : y (Euler angles, in rads)
76    Column 5 : z (Euler angles, in rads)
77
78    3dvolreg order is:
79
80    rot z
81    rot x
82    rot y
83    trans z
84    trans x
85    trans y
86
87    '''
88
89    fff = open( my_file ,'r')
90    raw_data = fff.readlines()
91    fff.close()
92
93    data_list1 = []
94    for line in raw_data:
95        aa = np.zeros(6)
96        w = line.translate(None,"[]")
97        x = w.split(',')
98        Nx = len(x)
99        if Nx < 6:
100            sys.exit("** Problem in this line:\n %s\n\n" % line)
101        for i in range (6):
102            aa[i] = float(x[i])
103
104        # convert rads to degs
105        aa[3:]*= 180./np.pi
106
107        # convert to 3dvolreg ordering
108        bb = Conv_TRxyz_to_RTzyx(aa)
109
110        data_list1.append(bb)
111
112    return np.array(data_list1)
113
114# ------------ convert between volreg ordering and trans_xyz rot_xyz -----
115
116def Conv_TRxyz_to_RTzyx(x):
117    '''trans (xyz) rot (xyz) ordering to volreg; inverse of
118onv_RTzyx_to_TRxyz
119
120    '''
121
122    if len(x) - 6:
123        sys.exit("** Problem converting matrices! Wrong number of components")
124
125    y = np.zeros(6)
126
127    y[0] = x[5]
128    y[1] = x[3]
129    y[2] = x[4]
130    y[3] = x[2]
131    y[4] = x[0]
132    y[5] = x[1]
133
134    return y
135
136def Conv_RTzyx_to_TRxyz(x):
137    '''volreg ordering to trans (xyz) rot (xyz); inverse of
138 Conv_TRxyz_to_RTzyx
139
140    '''
141
142    if len(x) - 6:
143        sys.exit("** Problem converting matrices! Wrong number of components")
144
145    y = np.zeros(6)
146
147    y[0] = x[4]
148    y[1] = x[5]
149    y[2] = x[3]
150    y[3] = x[1]
151    y[4] = x[2]
152    y[5] = x[0]
153
154    return y
155
156# --------------------- for RMS calcs ----------------------
157
158def Rx(x):
159
160    out = np.zeros((3,3))
161
162    out[0,0] = 1
163    out[1,1] = np.cos(x)
164    out[2,2] = out[1,1]
165    out[2,1] = np.sin(x)
166    out[1,2] = -out[2,1]
167
168    return out
169
170def Ry(x):
171
172    out = np.zeros((3,3))
173
174    out[1,1] = 1
175    out[0,0] = np.cos(x)
176    out[2,2] = out[0,0]
177    out[0,2] = np.sin(x)
178    out[2,0] = -out[0,2]
179
180    return out
181
182def Rz(x):
183
184    out = np.zeros((3,3))
185
186    out[2,2] = 1
187    out[0,0] = np.cos(x)
188    out[1,1] = out[0,0]
189    out[1,0] = np.sin(x)
190    out[0,1] = -out[1,0]
191
192    return out
193
194# input angles in order of z-, y- and x- rots.
195# calculation would be done as if L-multiplying in order of
196# Rx, then Ry, then Rz
197# can use optional 4th argument if input angles are in 'rad'
198# and not (default) 'deg'
199def Rzyx(c, b, a, atype='deg'):
200
201    if atype=='deg':
202        fac = np.pi/180.
203        a*=fac
204        b*=fac
205        c*=fac
206    elif not(atype =='rad'):
207        print "Error! unknown angle specification '%s'!", atype
208        sys.exit(5)
209
210    YX = np.dot(Ry(b), Rx(a))
211    ZYX = np.dot(Rz(c), YX)
212
213    return ZYX
214
215# Following Reuter et al. 2014.  Need 't' and 'R' to have the same
216# units (mm). For humans, 'R' here should prob be between 60-75 mm.
217def calc_RMS_from_MtR(M, t, R):
218
219    Mid = M - np.identity(3)
220    MidTMid = np.dot(np.transpose(Mid), Mid)
221    TrMids = np.trace(MidTMid)
222
223    out = R * R * TrMids / 5.
224    out+= np.dot(np.transpose(t), t)
225
226    return np.sqrt(out)
227
228def calc_RMS_from_6mot(x, brrad, atype='deg'):
229    '''Main function to calculate RMS from six solid body parameters:
230three translation (in mm), and 3 rotation (in deg, def; can also put
231in 'rad').
232
233INPUT
234
235     x : [req] An array of 6 numbers, which should be 3 translation in mm
236         and 3 rotation in either deg (def) or rad (signalled with 3rd
237         arg).  Order of components should be in same order as outputs
238         of 3dvolreg:
239
240             roll pitch yaw dS  dL  dP
241
242         where "roll" is rotation about I-S axis (shaking head "no"),
243         "pitch" is rotation around L-R axis (nodding, which is
244         common), and yaw is rotation about A-P axis.
245
246 brrad : [req] Single number, the length scale brain; e.g., approx radius of
247         the brain, like r~(3.*V/(4*np.pi))**(1./3). Required.
248
249 atype : [opt] Specify whether rotations are either 'deg' (def) or
250         'rad'.
251
252OUTPUT
253
254    Single floating point number, the RMS.
255
256    '''
257
258    Nx = len(x)
259
260    if not( Nx == 6 ):
261        sys.exit("** Error: wrong length %d in input array" % Nx)
262
263    # convert volreg ordering to Tx, Ty, Tz, Rx, Ry, Rz
264    y = Conv_RTzyx_to_TRxyz(x)
265
266    # NOTE: while rotations and translations can't be mixed, the order
267    # of components within rotations isn't so important, since the
268    # trace is taken, and likewise for rotations since it is a
269    # dot-product.  Also note, annoying in Taylor et al. (2016), the
270    # RMS formula is written incorrectly (it should be
271    # (r**2)/5*(M-I)... etc.; there is an extra "plus" incorrectly),
272    # but the actual calculation was correct.
273    MM = LSF.Rzyx(y[5], y[4], y[3], atype=atype)
274    tt = np.array(y[:3])
275    RR = brrad
276    y  = LSF.calc_RMS(MM,tt,RR)
277
278    return y
279
280
281
282
283# ===============================================================
284
285#if __name__ == "__main__":
286
287
288
289