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