mvpa.misc.fsl.base

1 # emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- 2 # vi: set ft=python sts=4 ts=4 sw=4 et: 3 ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ## 4 # 5 # See COPYING file distributed along with the PyMVPA package for the 6 # copyright and license terms. 7 # 8 ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ## 9 """Tiny snippets to interface with FSL easily.""" 10 11 __docformat__ = 'restructuredtext' 12 13 import numpy as N 14 15 from mvpa.misc.io import ColumnData 16 from mvpa.misc.support import Event 17 18 if __debug__: 19 from mvpa.base import debug 20 21

22 -class FslEV3(ColumnData):

23 """IO helper to read FSL's EV3 files. 24 25 This is a three-column textfile format that is used to specify stimulation 26 protocols for fMRI data analysis in FSL's FEAT module. 27 28 Data is always read as `float`. 29 """

30 - def __init__(self, source):

31 """Read and write FSL EV3 files. 32 33 :Parameter: 34 source: filename of an EV3 file 35 """ 36 # init data from known format 37 ColumnData.__init__(self, source, 38 header=['onsets', 'durations', 'intensities'], 39 sep=None, dtype=float)

40 41

42 - def getNEVs(self):

43 """Returns the number of EVs in the file. 44 """ 45 return self.nrows

46 47

48 - def getEV(self, evid):

49 """Returns a tuple of (onset time, simulus duration, intensity) for a 50 certain EV. 51 """ 52 return (self['onsets'][evid], 53 self['durations'][evid], 54 self['intensities'][evid])

55 56

57 - def tofile(self, filename):

58 """Write data to a FSL EV3 file. 59 """ 60 ColumnData.tofile(self, filename, 61 header=False, 62 header_order=['onsets', 'durations', 'intensities'], 63 sep=' ')

64 65

66 - def toEvents(self, **kwargs):

67 """Convert into a list of `Event` instances. 68 69 :Parameters: 70 kwargs 71 Any keyword arugment provided would be replicated, through all 72 the entries. Useful to specify label or even a chunk 73 """ 74 return \ 75 [Event(onset=self['onsets'][i], 76 duration=self['durations'][i], 77 features=[self['intensities'][i]], 78 **kwargs) 79 for i in xrange(self.nevs)]

80 81 82 onsets = property(fget=lambda self: self['onsets']) 83 durations = property(fget=lambda self: self['durations']) 84 intensities = property(fget=lambda self: self['intensities']) 85 nevs = property(fget=getNEVs)

86 87 88

89 -class McFlirtParams(ColumnData):

90 """Read and write McFlirt's motion estimation parameters from and to text 91 files. 92 """ 93 header_def = ['rot1', 'rot2', 'rot3', 'x', 'y', 'z'] 94

95 - def __init__(self, source):

96 """Initialize McFlirtParams 97 98 :Parameter: 99 source: str 100 Filename of a parameter file. 101 """ 102 ColumnData.__init__(self, source, 103 header=McFlirtParams.header_def, 104 sep=None, dtype=float)

105 106

107 - def tofile(self, filename):

108 """Write motion parameters to file. 109 """ 110 ColumnData.tofile(self, filename, 111 header=False, 112 header_order=McFlirtParams.header_def, 113 sep=' ')

114 115

116 - def plot(self):

117 """Produce a simple plot of the estimated translation and rotation 118 parameters using. 119 120 You still need to can pylab.show() or pylab.savefig() if you want to 121 see/get anything. 122 """ 123 # import internally as it takes some time and might not be needed most 124 # of the time 125 import pylab as P 126 127 # translations subplot 128 P.subplot(211) 129 P.plot(self.x) 130 P.plot(self.y) 131 P.plot(self.z) 132 P.ylabel('Translations in mm') 133 P.legend(('x', 'y', 'z'), loc=0) 134 135 # rotations subplot 136 P.subplot(212) 137 P.plot(self.rot1) 138 P.plot(self.rot2) 139 P.plot(self.rot3) 140 P.ylabel('Rotations in rad') 141 P.legend(('rot1', 'rot2', 'rot3'), loc=0)

142 143

144 - def toarray(self):

145 """Returns the data as an array with six columns (same order as in file). 146 """ 147 import numpy as N 148 149 # return as array with time axis first 150 return N.array([self[i] for i in McFlirtParams.header_def], 151 dtype='float').T

152 153

154 -class FslGLMDesign(object):

155 """Load FSL GLM design matrices from file. 156 157 Be aware that such a desig matrix has its regressors in columns and the 158 samples in its rows. 159 """

160 - def __init__(self, source):

161 """ 162 :Parameter: 163 source: filename 164 Compressed files will be read as well, if their filename ends with 165 '.gz'. 166 """ 167 # XXX maybe load from array as well 168 self._loadFile(source)

169 170

171 - def _loadFile(self, fname):

172 """Helper function to load GLM definition from a file. 173 """ 174 # header info 175 nwaves = 0 176 ntimepoints = 0 177 matrix_offset = 0 178 179 # open the file compressed or not 180 if fname.endswith('.gz'): 181 import gzip 182 fh = gzip.open(fname, 'r') 183 else: 184 fh = open(fname, 'r') 185 186 # read header 187 for i, line in enumerate(fh): 188 if line.startswith('/NumWaves'): 189 nwaves = int(line.split()[1]) 190 if line.startswith('/NumPoints'): 191 ntimepoints = int(line.split()[1]) 192 if line.startswith('/PPheights'): 193 self.ppheights = [float(i) for i in line.split()[1:]] 194 if line.startswith('/Matrix'): 195 matrix_offset = i + 1 196 197 # done with the header, now revert to NumPy's loadtxt for convenience 198 fh.close() 199 self.mat = N.loadtxt(fname, skiprows=matrix_offset) 200 201 # checks 202 if not self.mat.shape == (ntimepoints, nwaves): 203 raise IOError, "Design matrix file '%s' did not contain expected " \ 204 "matrix size (expected %s, got %s)" \ 205 % (fname, str((ntimepoints, nwaves)), self.mat.shape)

206 207

208 - def plot(self, style='lines', **kwargs):

209 """Visualize the design matrix. 210 211 :Parameters: 212 style: 'lines', 'matrix' 213 **kwargs: 214 Additional arguments will be passed to the corresponding matplotlib 215 plotting functions 'plot()' and 'pcolor()' for 'lines' and 'matrix' 216 plots respectively. 217 """ 218 # import internally as it takes some time and might not be needed most 219 # of the time 220 import pylab as P 221 222 if style == 'lines': 223 # common y-axis 224 yax = N.arange(0, self.mat.shape[0]) 225 axcenters = [] 226 col_offset = max(self.ppheights) 227 228 # for all columns 229 for i in xrange(self.mat.shape[1]): 230 axcenter = i * col_offset 231 P.plot(self.mat[:, i] + axcenter, yax, **kwargs) 232 axcenters.append(axcenter) 233 234 P.xticks(N.array(axcenters), range(self.mat.shape[1])) 235 elif style == 'matrix': 236 P.pcolor(self.mat, **kwargs) 237 ticks = N.arange(1, self.mat.shape[1]+1) 238 P.xticks(ticks - 0.5, ticks) 239 else: 240 raise ValueError, "Unknown plotting style '%s'" % style 241 242 # labels and turn y-axis upside down 243 P.ylabel('Samples (top to bottom)') 244 P.xlabel('Regressors') 245 P.ylim(self.mat.shape[0],0)

246 247

248 -def read_fsl_design(fsf_file):

249 """Reads an FSL FEAT design.fsf file and return the content as a dictionary. 250 251 :Parameters: 252 fsf_file : filename, file-like 253 """ 254 # This function was originally contributed by Russell Poldrack 255 256 if isinstance(fsf_file, basestring): 257 infile = open(fsf_file, 'r') 258 else: 259 infile = fsf_file 260 261 # target dict 262 fsl = {} 263 264 # loop over all lines 265 for line in infile: 266 line = line.strip() 267 # if there is nothing on the line, do nothing 268 if not line or line[0] == '#': 269 continue 270 271 # strip leading TCL 'set' 272 key, value = line.split(None, 2)[1:] 273 274 # fixup the 'y-' thing 275 if value == 'y-': 276 value = "y" 277 278 # special case of variable keyword 279 if line.count('_files('): 280 # e.g. feat_files(1) -> feat_files 281 key = key.split('(')[0] 282 283 # decide which type we have for the value 284 # int? 285 if value.isdigit(): 286 fsl[key] = int(value) 287 else: 288 # float? 289 try: 290 fsl[key] = float(value) 291 except ValueError: 292 # must be string then, but... 293 # sometimes there are quotes, sometimes not, but if the value 294 # should be a string we remove them, since the value is already 295 # of this type 296 fsl[key] = value.strip('"') 297 298 return fsl

299

Source Code for Module mvpa.misc.fsl.base