~jtv/corpusfiltergraph/cross-python

« back to all changes in this revision

Viewing changes to trunk/lib/corpusfg/io.py

  • Committer: tahoar
  • Date: 2012-05-02 15:46:23 UTC
  • Revision ID: svn-v4:bc069b21-dff4-4e29-a776-06a4e04bad4e::266
new layout. need to update code to use the new layout

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
1
#! /usr/bin/env python
 
2
# -*- coding, utf8 -*-
 
3
 
 
4
# Copyright (c) 2000-2007 Gary Strangman.  All rights reserved
 
5
#
 
6
# Disclaimer
 
7
#
 
8
# This software is provided "as-is".  There are no expressed or implied
 
9
# warranties of any kind, including, but not limited to, the warranties
 
10
# of merchantability and fittness for a given application.  In no event
 
11
# shall Gary Strangman be liable for any direct, indirect, incidental,
 
12
# special, exemplary or consequential damages (including, but not limited
 
13
# to, loss of use, data or profits, or business interruption) however
 
14
# caused and on any theory of liability, whether in contract, strict
 
15
# liability or tort (including negligence or otherwise) arising in any way
 
16
# out of the use of this software, even if advised of the possibility of
 
17
# such damage.
 
18
#
 
19
# Comments and/or additions are welcome (send e-mail to:
 
20
# strang@nmr.mgh.harvard.edu).
 
21
#
 
22
'''
 
23
Defines a number of functions for pseudo-command-line OS functionality.
 
24
 
 
25
        cd(directory)
 
26
        pwd                 <-- can be used WITHOUT parens
 
27
        ls(d='.')
 
28
        rename(from,to)
 
29
        get(namepatterns,verbose=1)
 
30
        getstrings(namepatterns,verbose=1)
 
31
        put(outlist,filename,writetype='w')
 
32
        aget(namepatterns,verbose=1)
 
33
        aput(outarray,filename,writetype='w')
 
34
        bget(filename,numslices=1,xsize=64,ysize=64)
 
35
        braw(filename,btype)
 
36
        bput(outarray,filename,writeheader=0,packstring='h',writetype='wb')
 
37
        mrget(filename)
 
38
        find_dirs(sourcedir)
 
39
'''
 
40
 
 
41
## CHANGES:
 
42
## =======
 
43
## 07-11-26 ... more numpy conversion work
 
44
## 06-08-07 ... converted to numpy, changed version to 0.6
 
45
## 06-02-03 ... added add2afnihistory() to load modify afni HISTORY_NOTEs,
 
46
##              and added that option to array2afni output
 
47
## 04-06-14 ... added getafniparam() to load in afni values from HEAD files
 
48
## 03-04-09 ... fixed brikget to load based on datatype, changed array2afni
 
49
##              so that the default sliceorder is altplus, not seqplus
 
50
## 02-11-20 ... added binget(), binput(), array2afni(), version 0.5
 
51
## 02-10-20 ... added find_dirs() function, changed version to 0.4
 
52
## 01-11-15 ... changed aput() and put() to accept a delimiter
 
53
## 01-04-19 ... added oneperline option to put() function
 
54
## 99-11-07 ... added DAs quick flat-text-file loaders, load() and fload()
 
55
## 99-11-01 ... added version number (0.1) for distribution
 
56
## 99-08-30 ... Put quickload in here
 
57
## 99-06-27 ... Changed bget thing back ... confused ...
 
58
## 99-06-24 ... exchanged xsize and ysize in bget for non-square images (NT??)
 
59
##              modified bget to raise an IOError when file not found
 
60
## 99-06-12 ... added load() and save() aliases for aget() and aput() (resp.)
 
61
## 99-04-13 ... changed aget() to ignore (!!!!) lines beginning with # or %
 
62
## 99-01-17 ... changed get() so ints come in as ints (not floats)
 
63
##
 
64
 
 
65
 
 
66
 
 
67
try:
 
68
        import mmapfile
 
69
except:
 
70
        pass
 
71
 
 
72
import pstat
 
73
import glob, re, string, types, os, struct, copy, time, tempfile, sys
 
74
from types import *
 
75
import numpy as N
 
76
 
 
77
__version__ = 0.6
 
78
 
 
79
def wrap(f):
 
80
        '''
 
81
Wraps a function so that if it's entered *by itself*
 
82
in the interpreter without ()'s, it gets called anyway
 
83
'''
 
84
        class W:
 
85
                def __init__(self, f):
 
86
                        self.f = f
 
87
                def __repr__(self):
 
88
                        x =apply(self.f)
 
89
                        if x:
 
90
                                return repr(x)
 
91
                        else:
 
92
                                return ''
 
93
        return W(f)
 
94
 
 
95
 
 
96
def cd(directory):
 
97
        '''
 
98
Changes the working python directory for the interpreter.
 
99
 
 
100
Usage:   cd(directory)      where 'directory' is a string
 
101
'''
 
102
        os.chdir(directory)
 
103
        return
 
104
 
 
105
 
 
106
def pwd():
 
107
        '''
 
108
Changes the working python directory for the interpreter.
 
109
 
 
110
Usage:   pwd       (no parens needed)
 
111
'''
 
112
        return os.getcwd()
 
113
pwd = wrap(pwd)
 
114
 
 
115
 
 
116
def ls(d='.'):
 
117
        '''
 
118
Produces a directory listing.  Default is the current directory.
 
119
 
 
120
Usage:   ls(d='.')
 
121
'''
 
122
        os.system('ls '+d)
 
123
        return None
 
124
 
 
125
 
 
126
def rename(source, dest):
 
127
        '''
 
128
Renames files specified by UNIX inpattern to those specified by UNIX
 
129
outpattern.  Can only handle a single '*' in the two patterns!!!
 
130
 
 
131
Usage:   rename(source, dest)     e.g., rename('*.txt', '*.c')
 
132
'''
 
133
        infiles = glob.glob(source)
 
134
        outfiles = []
 
135
        incutindex = string.index(source,'*')
 
136
        outcutindex = string.index(source,'*')
 
137
        findpattern1 = source[0:incutindex]
 
138
        findpattern2 = source[incutindex+1:]
 
139
        replpattern1 = dest[0:incutindex]
 
140
        replpattern2 = dest[incutindex+1:]
 
141
        for fname in infiles:
 
142
                if incutindex > 0:
 
143
                        newname = re.sub(findpattern1,replpattern1,fname,1)
 
144
                if outcutindex < len(dest)-1:
 
145
                        if incutindex > 0:
 
146
                                lastone = string.rfind(newname,replpattern2)
 
147
                                newname = newname[0:lastone] + re.sub(findpattern2,replpattern2,fname[lastone:],1)
 
148
                        else:
 
149
                                lastone = string.rfind(fname,findpattern2)
 
150
                                if lastone <> -1:
 
151
                                        newname = fname[0:lastone]
 
152
                                        newname = newname + re.sub(findpattern2,replpattern2,fname[lastone:],1)
 
153
                print fname, newname
 
154
                os.rename(fname,newname)
 
155
        return
 
156
 
 
157
 
 
158
def get(namepatterns,verbose=1):
 
159
        '''
 
160
Loads a list of lists from text files (specified by a UNIX-style
 
161
wildcard filename pattern) and converts all numeric values to floats.
 
162
Uses the glob module for filename pattern conversion.  Loaded filename
 
163
is printed if verbose=1.
 
164
 
 
165
Usage:   get(namepatterns,verbose=1)
 
166
Returns: a 1D or 2D list of lists from whitespace delimited text files
 
167
                 specified by namepatterns; numbers that can be converted to floats
 
168
                 are so converted
 
169
'''
 
170
        fnames = []
 
171
        if type(namepatterns) in [ListType,TupleType]:
 
172
                for item in namepatterns:
 
173
                        fnames = fnames + glob.glob(item)
 
174
        else:
 
175
                fnames = glob.glob(namepatterns)
 
176
        
 
177
        if len(fnames) == 0:
 
178
                if verbose:
 
179
                        print 'NO FILENAMES MATCH ('+namepatterns+') !!'
 
180
                return None
 
181
 
 
182
        if verbose:
 
183
                print fnames             # so user knows what has been loaded
 
184
        elements = []
 
185
        for i in range(len(fnames)):
 
186
                file = open(fnames[i])
 
187
                newelements = map(string.split,file.readlines())
 
188
                for i in range(len(newelements)):
 
189
                        for j in range(len(newelements[i])):
 
190
                                try:
 
191
                                        newelements[i][j] = string.atoi(newelements[i][j])
 
192
                                except ValueError:
 
193
                                        try:
 
194
                                                newelements[i][j] = string.atof(newelements[i][j])
 
195
                                        except:
 
196
                                                pass
 
197
                elements = elements + newelements
 
198
        if len(elements)==1:  elements = elements[0]
 
199
        return elements
 
200
 
 
201
 
 
202
def getstrings(namepattern,verbose=1):
 
203
        '''
 
204
Loads a (set of) text file(s), with all elements left as string type.
 
205
Uses UNIX-style wildcards (i.e., function uses glob).  Loaded filename
 
206
is printed if verbose=1.
 
207
 
 
208
Usage:   getstrings(namepattern, verbose=1)
 
209
Returns: a list of strings, one per line in each text file specified by
 
210
                 namepattern
 
211
'''
 
212
        fnames = glob.glob(namepattern)
 
213
        if len(fnames) == 0:
 
214
                if verbose:
 
215
                        print 'NO FILENAMES MATCH ('+namepattern+') !!'
 
216
                return None
 
217
        if verbose:
 
218
                print fnames
 
219
        elements = []
 
220
        for filename in fnames:
 
221
                file = open(filename)
 
222
                newelements = map(string.split,file.readlines())
 
223
                elements = elements + newelements
 
224
        return elements
 
225
 
 
226
 
 
227
def put(outlist,fname,writetype='w',oneperline=0,delimit=' '):
 
228
        '''
 
229
Writes a passed mixed-type list (str and/or numbers) to an output
 
230
file, and then closes the file.  Default is overwrite the destination
 
231
file.
 
232
 
 
233
Usage:   put(outlist,fname,writetype='w',oneperline=0,delimit=' ')
 
234
Returns: None
 
235
'''
 
236
        if type(outlist) in [N.ndarray]:
 
237
                aput(outlist,fname,writetype)
 
238
                return
 
239
        if type(outlist[0]) not in [ListType,TupleType]:  # 1D list
 
240
                outfile = open(fname,writetype)
 
241
                if not oneperline:
 
242
                        outlist = pstat.list2string(outlist,delimit)
 
243
                        outfile.write(outlist)
 
244
                        outfile.write('\n')
 
245
                else:  # they want one element from the list on each file line
 
246
                        for item in outlist:
 
247
                                outfile.write(str(item)+'\n')
 
248
                outfile.close()
 
249
        else:                                             # 2D list (list-of-lists)
 
250
                outfile = open(fname,writetype)
 
251
                for row in outlist:
 
252
                        outfile.write(pstat.list2string(row,delimit))
 
253
                        outfile.write('\n')
 
254
                outfile.close()
 
255
        return None
 
256
 
 
257
 
 
258
def isstring(x):
 
259
        if type(x)==StringType:
 
260
                return 1
 
261
        else:
 
262
                return 0
 
263
 
 
264
 
 
265
 
 
266
def aget(namepattern,verbose=1):
 
267
        '''
 
268
Loads an array from 2D text files (specified by a UNIX-style wildcard
 
269
filename pattern).  ONLY 'GET' FILES WITH EQUAL NUMBERS OF COLUMNS
 
270
ON EVERY ROW (otherwise returned array will be zero-dimensional).
 
271
 
 
272
Usage:   aget(namepattern)
 
273
Returns: an array of integers, floats or objects (type='O'), depending on the
 
274
                 contents of the files specified by namepattern
 
275
'''
 
276
        fnames = glob.glob(namepattern)
 
277
        if len(fnames) == 0:
 
278
                if verbose:
 
279
                        print 'NO FILENAMES MATCH ('+namepattern+') !!'
 
280
                        return None
 
281
        if verbose:
 
282
                print fnames
 
283
        elements = []
 
284
        for filename in fnames:
 
285
                file = open(filename)
 
286
                newelements = file.readlines()
 
287
                del_list = []
 
288
                for row in range(len(newelements)):
 
289
                        if (newelements[row][0]=='%' or newelements[row][0]=='#'
 
290
                                or len(newelements[row])==1 or newelements[row][0]=='\r'):
 
291
                                del_list.append(row)
 
292
                del_list.reverse()
 
293
                for i in del_list:
 
294
                        newelements.pop(i)
 
295
                newelements = map(string.split,newelements)
 
296
                for i in range(len(newelements)):
 
297
                        for j in range(len(newelements[i])):
 
298
                                try:
 
299
                                        newelements[i][j] = string.atof(newelements[i][j])
 
300
                                except:
 
301
                                        pass
 
302
                elements = elements + newelements
 
303
        for row in range(len(elements)):
 
304
                if N.add.reduce(N.array(map(isstring,elements[row])))==len(elements[row]):
 
305
                        print "A row of strings was found.  Returning a LIST."
 
306
                        return elements
 
307
        try:
 
308
                elements = N.array(elements)
 
309
        except TypeError:
 
310
                elements = N.array(elements,dtype='O')
 
311
        return elements
 
312
 
 
313
 
 
314
def aput(outarray,fname,writetype='w',delimit=' '):
 
315
        '''
 
316
Sends passed 1D or 2D array to an output file and closes the file.
 
317
 
 
318
Usage:   aput(outarray,fname,writetype='w',delimit=' ')
 
319
Returns: None
 
320
'''
 
321
        outfile = open(fname,writetype)
 
322
        if len(outarray.shape) == 1:
 
323
                outarray = outarray[N.newaxis,:]
 
324
        if len(outarray.shape) > 2:
 
325
                raise TypeError, "put() and aput() require 1D or 2D arrays.  Otherwise use some kind of pickling."
 
326
        else: # must be a 2D array
 
327
                for row in outarray:
 
328
                        outfile.write(string.join(map(str,row),delimit))
 
329
                        outfile.write('\n')
 
330
                outfile.close()
 
331
        return None
 
332
 
 
333
 
 
334
def bget(imfile,shp=None,unpackstr=N.int16,bytesperpixel=2.0,sliceinit=0):
 
335
        '''
 
336
Reads in a binary file, typically with a .bshort or .bfloat extension.
 
337
If so, the last 3 parameters are set appropriately.  If not, the last 3
 
338
parameters default to reading .bshort files (2-byte integers in big-endian
 
339
binary format).
 
340
 
 
341
Usage:   bget(imfile,shp=None,unpackstr=N.int16,bytesperpixel=2.0,sliceinit=0)
 
342
'''
 
343
        if imfile[:3] == 'COR':
 
344
                return CORget(imfile)
 
345
        if imfile[-2:] == 'MR':
 
346
                return mrget(imfile,unpackstr)
 
347
        if imfile[-4:] == 'BRIK':
 
348
                return brikget(imfile,unpackstr,shp)
 
349
        if imfile[-3:] in ['mnc','MNC','inc','INC']:
 
350
                return mincget(imfile,unpackstr,shp)
 
351
        if imfile[-3:] == 'img':
 
352
                return mghbget(imfile,unpackstr,shp)
 
353
        if imfile[-6:] == 'bshort' or imfile[-6:] == 'bfloat':
 
354
                if shp == None:
 
355
                        return mghbget(imfile,unpackstr=unpackstr,bytesperpixel=bytesperpixel,sliceinit=sliceinit)
 
356
                else:
 
357
                        return mghbget(imfile,shp[0],shp[1],shp[2],unpackstr,bytesperpixel,sliceinit)
 
358
 
 
359
 
 
360
def CORget(infile):
 
361
        '''
 
362
Reads a binary COR-nnn file (flattening file).
 
363
 
 
364
Usage:   CORget(imfile)
 
365
Returns: 2D array of 16-bit ints
 
366
'''
 
367
        d=braw(infile,N.int8)
 
368
        d.shape = (256,256)
 
369
        d = N.where(d>=0,d,256+d)
 
370
        return d
 
371
 
 
372
 
 
373
def mincget(imfile,unpackstr=N.int16,shp=None):
 
374
        '''
 
375
Loads in a .MNC file.
 
376
 
 
377
Usage:  mincget(imfile,unpackstr=N.int16,shp=None)  default shp = -1,20,64,64
 
378
'''
 
379
        if shp == None:
 
380
                shp = (-1,20,64,64)
 
381
        os.system('mincextract -short -range 0 4095 -image_range 0 4095 ' +
 
382
                          imfile+' > minctemp.bshort')
 
383
        try:
 
384
                d = braw('minctemp.bshort',unpackstr)
 
385
        except:
 
386
                print "Couldn't find file:  "+imfile
 
387
                raise IOError, "Couldn't find file in mincget()"
 
388
 
 
389
        print shp, d.shape
 
390
        d.shape = shp
 
391
        os.system('rm minctemp.bshort')
 
392
        return d
 
393
 
 
394
 
 
395
def brikget(imfile,unpackstr=N.int16,shp=None):
 
396
        '''
 
397
Gets an AFNI BRIK file.
 
398
 
 
399
Usage:  brikget(imfile,unpackstr=N.int16,shp=None)  default shp: (-1,48,61,51)
 
400
'''
 
401
        if shp == None:
 
402
                shp = (-1,48,61,51)
 
403
        try:
 
404
                file = open(imfile, "rb")
 
405
        except:
 
406
                print "Couldn't find file:  "+imfile
 
407
                raise IOError, "Couldn't find file in brikget()"
 
408
        try:
 
409
                header = imfile[0:-4]+'HEAD'
 
410
                lines = open(header).readlines()
 
411
                for i in range(len(lines)):
 
412
                        if string.find(lines[i],'DATASET_DIMENSIONS') <> -1:
 
413
                                dims = string.split(lines[i+2][0:string.find(lines[i+2],' 0')])
 
414
                                dims = map(string.atoi,dims)
 
415
                        if string.find(lines[i],'BRICK_FLOAT_FACS') <> -1:
 
416
                                count = string.atoi(string.split(lines[i+1])[2])
 
417
                                mults = []
 
418
                                for j in range(int(N.ceil(count/5.))):
 
419
                                        mults += map(string.atof,string.split(lines[i+2+j]))
 
420
                                mults = N.array(mults)
 
421
                        if string.find(lines[i],'BRICK_TYPES') <> -1:
 
422
                                first5 = lines[i+2]
 
423
                                first5 = map(string.atoi,string.split(first5))
 
424
                                if first5[0] == 0:
 
425
                                        unpackstr = N.uint8
 
426
                                elif first5[0] == 1:
 
427
                                        unpackstr = N.int16
 
428
                                elif first5[0] == 3:
 
429
                                        unpackstr = N.float32
 
430
                                elif first5[0] == 5:
 
431
                                        unpackstr = N.complex32
 
432
                dims.reverse()
 
433
                shp = [-1]+dims
 
434
        except IOError:
 
435
                print "No header file.  Continuing ..."
 
436
        lines = None
 
437
 
 
438
        print shp
 
439
        print 'Using unpackstr:',unpackstr  #,', bytesperpixel=',bytesperpixel
 
440
 
 
441
        file = open(imfile, "rb")
 
442
        bdata = file.read()
 
443
 
 
444
        # the > forces big-endian (for or from Sun/SGI)
 
445
        bdata = N.fromstring(bdata,unpackstr)
 
446
#    littleEndian = ( struct.pack('i',1)==struct.pack('<i',1) )
 
447
        if (max(bdata)>1e30):
 
448
                bdata = bdata.byteswap()
 
449
        try:
 
450
                bdata.shape = shp
 
451
        except:
 
452
                print 'Incorrect shape ...',shp,len(bdata)
 
453
                raise ValueError, 'Incorrect shape for file size'
 
454
        if len(bdata) == 1:
 
455
                bdata = bdata[0]
 
456
 
 
457
        if N.sum(mults) == 0:
 
458
                return bdata
 
459
        try:
 
460
                multshape = [1]*len(bdata.shape)
 
461
                for i in range(len(bdata.shape)):
 
462
                        if len(mults) == bdata.shape[i]:
 
463
                                multshape[i] = len(mults)
 
464
                                break
 
465
                mults.shape = multshape
 
466
                return bdata*mults
 
467
        except:
 
468
                return bdata
 
469
 
 
470
def mghbget(imfile,numslices=-1,xsize=64,ysize=64,
 
471
                   unpackstr=N.int16,bytesperpixel=2.0,sliceinit=0):
 
472
        '''
 
473
Reads in a binary file, typically with a .bshort or .bfloat extension.
 
474
If so, the last 3 parameters are set appropriately.  If not, the last 3
 
475
parameters default to reading .bshort files (2-byte integers in big-endian
 
476
binary format).
 
477
 
 
478
Usage:   mghbget(imfile, numslices=-1, xsize=64, ysize=64,
 
479
                                unpackstr=N.int16, bytesperpixel=2.0, sliceinit=0)
 
480
'''
 
481
        try:
 
482
                file = open(imfile, "rb")
 
483
        except:
 
484
                print "Couldn't find file:  "+imfile
 
485
                raise IOError, "Couldn't find file in bget()"
 
486
        try:
 
487
                header = imfile[0:-6]+'hdr'
 
488
                vals = get(header,0)  # '0' means no missing-file warning msg
 
489
                if type(vals[0]) == ListType:  # it's an extended header
 
490
                        xsize = int(vals[0][0])
 
491
                        ysize = int(vals[0][1])
 
492
                        numslices = int(vals[0][2])
 
493
                else:
 
494
                        xsize = int(vals[0])
 
495
                        ysize = int(vals[1])
 
496
                        numslices = int(vals[2])
 
497
        except:
 
498
                print "No header file.  Continuing ..."
 
499
 
 
500
        suffix = imfile[-6:]
 
501
        if suffix == 'bshort':
 
502
                pass
 
503
        elif suffix[-3:] == 'img':
 
504
                pass
 
505
        elif suffix == 'bfloat':
 
506
                unpackstr = N.float32
 
507
                bytesperpixel = 4.0
 
508
                sliceinit = 0.0
 
509
        else:
 
510
                print 'Not a bshort, bfloat or img file.'
 
511
                print 'Using unpackstr:',unpackstr,', bytesperpixel=',bytesperpixel
 
512
 
 
513
        imsize = xsize*ysize
 
514
        file = open(imfile, "rb")
 
515
        bdata = file.read()
 
516
 
 
517
        numpixels = len(bdata) / bytesperpixel
 
518
        if numpixels%1 != 0:
 
519
                raise ValueError, "Incorrect file size in fmri.bget()"
 
520
        else:  # the > forces big-endian (for or from Sun/SGI)
 
521
                bdata = N.fromstring(bdata,unpackstr)
 
522
#        littleEndian = ( struct.pack('i',1)==struct.pack('<i',1) )
 
523
#        if littleEndian:
 
524
#            bdata = bdata.byteswap()
 
525
                if (max(bdata)>1e30):
 
526
                        bdata = bdata.byteswap()
 
527
        if suffix[-3:] == 'img':
 
528
                if numslices == -1:
 
529
                        numslices = len(bdata)/8200  # 8200=(64*64*2)+8 bytes per image
 
530
                        xsize = 64
 
531
                        ysize = 128
 
532
                slices = N.zeros((numslices,xsize,ysize),N.int32)
 
533
                for i in range(numslices):
 
534
                        istart = i*8 + i*xsize*ysize
 
535
                        iend = i*8 + (i+1)*xsize*ysize
 
536
                        print i, istart,iend
 
537
                        slices[i] = N.reshape(N.array(bdata[istart:iend]),(xsize,ysize))
 
538
        else:
 
539
                if numslices == 1:
 
540
                        slices = N.reshape(N.array(bdata),[xsize,ysize])
 
541
                else:
 
542
                        slices = N.reshape(N.array(bdata),[numslices,xsize,ysize])
 
543
        if len(slices) == 1:
 
544
                slices = slices[0]
 
545
        return slices
 
546
 
 
547
 
 
548
def braw(fname,btype,shp=None):
 
549
        '''
 
550
Opens a binary file, unpacks it, and returns a flat array of the
 
551
type specified.  Use Numeric types ... N.float32, N.int64, etc.
 
552
 
 
553
Usage:   braw(fname,btype,shp=None)
 
554
Returns: flat array of floats, or ints (if btype=N.int16)
 
555
'''
 
556
        file = open(fname,'rb')
 
557
        bdata = file.read()
 
558
        bdata = N.fromstring(bdata,btype)
 
559
#    littleEndian = ( struct.pack('i',1)==struct.pack('<i',1) )
 
560
#    if littleEndian:
 
561
#        bdata = bdata.byteswap()  # didn't used to need this with '>' above
 
562
        if (max(bdata)>1e30):
 
563
                bdata = bdata.byteswap()
 
564
        if shp:
 
565
                try:
 
566
                        bdata.shape = shp
 
567
                        return bdata
 
568
                except:
 
569
                        pass
 
570
        return N.array(bdata)
 
571
 
 
572
 
 
573
def glget(fname,btype):
 
574
        '''
 
575
Load in a file containing pixels from glReadPixels dump.
 
576
 
 
577
Usage:   glget(fname,btype)
 
578
Returns: array of 'btype elements with shape 'shape', suitable for im.ashow()
 
579
'''
 
580
        d = braw(fname,btype)
 
581
        d = d[8:]
 
582
        f = open(fname,'rb')
 
583
        shp = f.read(8)
 
584
        f.close()
 
585
        shp = N.fromstring(shp,N.int32)
 
586
        shp[0],shp[1] = shp[1],shp[0]
 
587
        try:
 
588
                carray = N.reshape(d,shp)
 
589
                return
 
590
        except:
 
591
                pass
 
592
        try:
 
593
                r = d[0::3]+0
 
594
                g = d[1::3]+0
 
595
                b = d[2::3]+0
 
596
                r.shape = shp
 
597
                g.shape = shp
 
598
                b.shape = shp
 
599
                carray = N.array([r,g,b])
 
600
        except:
 
601
                outstr = "glget: shape not correct for data of length "+str(len(d))
 
602
                raise ValueError, outstr
 
603
        return carray
 
604
 
 
605
 
 
606
def mget(fname,btype):
 
607
        '''
 
608
Load in a file that was saved from matlab
 
609
 
 
610
Usage:   mget(fname,btype)
 
611
'''
 
612
        d = braw(fname,btype)
 
613
        try:
 
614
                header = fname[0:-6]+'hdr'
 
615
                vals = get(header,0)  # '0' means no missing-file warning msg
 
616
                if type(vals[0]) == ListType:  # it's an extended header
 
617
                        xsize = int(vals[0][0])
 
618
                        ysize = int(vals[0][1])
 
619
                        numslices = int(vals[0][2])
 
620
                else:
 
621
                        xsize = int(vals[0])
 
622
                        ysize = int(vals[1])
 
623
                        numslices = int(vals[2])
 
624
                print xsize,ysize,numslices, d.shape
 
625
        except:
 
626
                print "No header file.  Continuing ..."
 
627
        if numslices == 1:
 
628
                d.shape = [ysize,xsize]
 
629
                return N.transpose(d)*1
 
630
        else:
 
631
                d.shape = [numslices,ysize,xsize]
 
632
                return N.transpose(d)*1
 
633
 
 
634
 
 
635
def mput(outarray,fname,writeheader=0,btype=N.int16):
 
636
        '''
 
637
Save a file for use in matlab.
 
638
'''
 
639
        outarray = N.transpose(outarray)
 
640
        outdata = N.ravel(outarray).astype(btype)
 
641
        outdata = outdata.tostring()
 
642
        outfile = open(fname,'wb')
 
643
        outfile.write(outdata)
 
644
        outfile.close()
 
645
        if writeheader == 1:
 
646
                try:
 
647
                        suffixindex = string.rfind(fname,'.')
 
648
                        hdrname = fname[0:suffixindex]
 
649
                except ValueError:
 
650
                        hdrname = fname
 
651
                if len(outarray.shape) == 2:
 
652
                        hdr = [outarray.shape[1],outarray.shape[0], 1, 0]
 
653
                else:
 
654
                        hdr = [outarray.shape[2],outarray.shape[1],outarray.shape[0], 0,'\n']
 
655
                print hdrname+'.hdr'
 
656
                outfile = open(hdrname+'.hdr','w')
 
657
                outfile.write(pstat.list2string(hdr))
 
658
                outfile.close()
 
659
        return None
 
660
 
 
661
 
 
662
def bput(outarray,fname,writeheader=0,packtype=N.int16,writetype='wb'):
 
663
        '''
 
664
Writes the passed array to a binary output file, and then closes
 
665
the file.  Default is overwrite the destination file.
 
666
 
 
667
Usage:   bput (outarray,filename,writeheader=0,packtype=N.int16,writetype='wb')
 
668
'''
 
669
        suffix = fname[-6:]
 
670
        if suffix == 'bshort':
 
671
                packtype = N.int16
 
672
        elif suffix == 'bfloat':
 
673
                packtype = N.float32
 
674
        else:
 
675
                print 'Not a bshort or bfloat file.  Using packtype=',packtype
 
676
 
 
677
        outdata = N.ravel(outarray).astype(packtype)
 
678
#    littleEndian = ( struct.pack('i',1)==struct.pack('<i',1) )
 
679
#    if littleEndian:
 
680
#        outdata = outdata.byteswap()
 
681
        outdata = outdata.tostring()
 
682
        outfile = open(fname,writetype)
 
683
        outfile.write(outdata)
 
684
        outfile.close()
 
685
        if writeheader == 1:
 
686
                try:
 
687
                        suffixindex = string.rfind(fname,'.')
 
688
                        hdrname = fname[0:suffixindex]
 
689
                except ValueError:
 
690
                        hdrname = fname
 
691
                if len(outarray.shape) == 2:
 
692
                        hdr = [outarray.shape[0],outarray.shape[1], 1, 0]
 
693
                else:
 
694
                        hdr = [outarray.shape[1],outarray.shape[2],outarray.shape[0], 0,'\n']
 
695
                print hdrname+'.hdr'
 
696
                outfile = open(hdrname+'.hdr','w')
 
697
                outfile.write(pstat.list2string(hdr))
 
698
                outfile.close()
 
699
        return None
 
700
 
 
701
 
 
702
def mrget(fname,datatype=N.int16):
 
703
        '''
 
704
Opens a binary .MR file and clips off the tail data portion of it, returning
 
705
the result as an array.
 
706
 
 
707
Usage:   mrget(fname,datatype=N.int16)
 
708
'''
 
709
        d = braw(fname,datatype)
 
710
        if len(d) > 512*512:
 
711
                return N.reshape(d[-512*512:],(512,512))
 
712
        elif len(d) > 320*320:
 
713
                return N.reshape(d[-320*320:],(320,320))
 
714
        elif len(d) > 256*256:
 
715
                return N.reshape(d[-256*256:],(256,256))
 
716
        elif len(d) > 128*128:
 
717
                return N.reshape(d[-128*128:],(128,128))
 
718
        elif len(d) > 64*64:
 
719
                return N.reshape(d[-64*64:],(64,64))
 
720
        else:
 
721
                return N.reshape(d[-32*32:],(32,32))
 
722
 
 
723
 
 
724
def quickload(fname,linestocut=4):
 
725
        '''
 
726
Quickly loads in a long text file, chopping off first n 'linestocut'.
 
727
 
 
728
Usage:   quickload(fname,linestocut=4)
 
729
Returns: array filled with data in fname
 
730
'''
 
731
        f = open(fname,'r')
 
732
        d = f.readlines()
 
733
        f.close()
 
734
        print fname,'read in.'
 
735
        d = d[linestocut:]
 
736
        d = map(string.split,d)
 
737
        print 'Done with string.split on lines.'
 
738
        for i in range(len(d)):
 
739
                d[i] = map(string.atoi,d[i])
 
740
        print 'Conversion to ints done.'
 
741
        return N.array(d)
 
742
 
 
743
def writedelimited(listoflists, delimiter, file, writetype='w'):
 
744
        '''
 
745
Writes a list of lists in columns, separated by character(s) delimiter
 
746
to specified file.  File-overwrite is the default.
 
747
 
 
748
Usage:   writedelimited(listoflists,delimiter,filename,writetype='w')
 
749
Returns: None
 
750
'''
 
751
        if type(listoflists[0]) not in [ListType,TupleType]:
 
752
                listoflists = [listoflists]
 
753
        outfile = open(file,writetype)
 
754
        rowstokill = []
 
755
        list2print = copy.deepcopy(listoflists)
 
756
        for i in range(len(listoflists)):
 
757
                if listoflists[i] == ['\n'] or listoflists[i]=='\n' or listoflists[i]=='dashes':
 
758
                        rowstokill = rowstokill + [i]
 
759
        rowstokill.reverse()
 
760
        for row in rowstokill:
 
761
                del list2print[row]
 
762
        maxsize = [0]*len(list2print[0])
 
763
        for row in listoflists:
 
764
                if row == ['\n'] or row == '\n':
 
765
                        outfile.write('\n')
 
766
                elif row == ['dashes'] or row == 'dashes':
 
767
                        dashes = [0]*len(maxsize)
 
768
                        for j in range(len(maxsize)):
 
769
                                dashes[j] = '------'
 
770
                        outfile.write(pstat.linedelimited(dashes,delimiter))
 
771
                else:
 
772
                        outfile.write(pstat.linedelimited(row,delimiter))
 
773
                outfile.write('\n')
 
774
        outfile.close()
 
775
        return None
 
776
 
 
777
def writecc(listoflists,file,writetype='w',extra=2):
 
778
        '''
 
779
Writes a list of lists to a file in columns, customized by the max
 
780
size of items within the columns (max size of items in col, +2 characters)
 
781
to specified file.  File-overwrite is the default.
 
782
 
 
783
Usage:   writecc(listoflists,file,writetype='w',extra=2)
 
784
Returns: None
 
785
'''
 
786
        if type(listoflists[0]) not in [ListType,TupleType]:
 
787
                listoflists = [listoflists]
 
788
        outfile = open(file,writetype)
 
789
        rowstokill = []
 
790
        list2print = copy.deepcopy(listoflists)
 
791
        for i in range(len(listoflists)):
 
792
                if listoflists[i] == ['\n'] or listoflists[i]=='\n' or listoflists[i]=='dashes':
 
793
                        rowstokill = rowstokill + [i]
 
794
        rowstokill.reverse()
 
795
        for row in rowstokill:
 
796
                del list2print[row]
 
797
        maxsize = [0]*len(list2print[0])
 
798
        for col in range(len(list2print[0])):
 
799
                items = pstat.colex(list2print,col)
 
800
                items = map(pstat.makestr,items)
 
801
                maxsize[col] = max(map(len,items)) + extra
 
802
        for row in listoflists:
 
803
                if row == ['\n'] or row == '\n':
 
804
                        outfile.write('\n')
 
805
                elif row == ['dashes'] or row == 'dashes':
 
806
                        dashes = [0]*len(maxsize)
 
807
                        for j in range(len(maxsize)):
 
808
                                dashes[j] = '-'*(maxsize[j]-2)
 
809
                        outfile.write(pstat.lineincustcols(dashes,maxsize))
 
810
                else:
 
811
                        outfile.write(pstat.lineincustcols(row,maxsize))
 
812
                outfile.write('\n')
 
813
        outfile.close()
 
814
        return None
 
815
 
 
816
 
 
817
def writefc(listoflists,colsize,file,writetype='w'):
 
818
        '''
 
819
Writes a list of lists to a file in columns of fixed size.  File-overwrite
 
820
is the default.
 
821
 
 
822
Usage:   writefc(listoflists,colsize,file,writetype='w')
 
823
Returns: None
 
824
'''
 
825
        if type(listoflists) == N.ndarray:
 
826
                listoflists = listoflists.tolist()
 
827
        if type(listoflists[0]) not in [ListType,TupleType]:
 
828
                listoflists = [listoflists]
 
829
        outfile = open(file,writetype)
 
830
        rowstokill = []
 
831
        list2print = copy.deepcopy(listoflists)
 
832
        for i in range(len(listoflists)):
 
833
                if listoflists[i] == ['\n'] or listoflists[i]=='\n' or listoflists[i]=='dashes':
 
834
                        rowstokill = rowstokill + [i]
 
835
        rowstokill.reverse()
 
836
        for row in rowstokill:
 
837
                del list2print[row]
 
838
        n = [0]*len(list2print[0])
 
839
        for row in listoflists:
 
840
                if row == ['\n'] or row == '\n':
 
841
                        outfile.write('\n')
 
842
                elif row == ['dashes'] or row == 'dashes':
 
843
                        dashes = [0]*colsize
 
844
                        for j in range(len(n)):
 
845
                                dashes[j] = '-'*(colsize)
 
846
                        outfile.write(pstat.lineincols(dashes,colsize))
 
847
                else:
 
848
                        outfile.write(pstat.lineincols(row,colsize))
 
849
                outfile.write('\n')
 
850
        outfile.close()
 
851
        return None
 
852
 
 
853
 
 
854
def load(fname,lines_to_ignore=4,type='i'):
 
855
        '''
 
856
Load in huge, flat, 2D text files.  Can handle differing line-lengths AND
 
857
can strip #/% on UNIX (or with a better NT grep).  Requires wc, grep, and
 
858
mmapfile.lib/.pyd. Type can be 'i', 'f' or 'd', for ints, floats or doubles,
 
859
respectively.  Lines_to_ignore determines how many lines at the start of the
 
860
file to ignore (required for non-working grep).
 
861
 
 
862
Usage:   load(fname,lines_to_ignore=4,type='i')
 
863
Returns: numpy array of specified type
 
864
'''
 
865
        start = time.time()      ## START TIMER
 
866
        if type == 'i':
 
867
                intype = int
 
868
        elif type in ['f','d']:
 
869
                intype = float
 
870
        else:
 
871
                raise ValueError, "type can be 'i', 'f' or 'd' in load()"
 
872
 
 
873
        ## STRIP OUT % AND # LINES
 
874
        tmpname = tempfile.mktemp()
 
875
        if sys.platform == 'win32':
 
876
                # NT VERSION OF GREP DOESN'T DO THE STRIPPING ... SIGH
 
877
                cmd = "grep.exe -v \'%\' "+fname+" > "+tmpname
 
878
                print cmd
 
879
                os.system(cmd)
 
880
        else:
 
881
                # UNIX SIDE SHOULD WORK
 
882
                cmd = "cat "+fname+" | grep -v \'%\' |grep -v \'#\' > "+tmpname
 
883
                print cmd
 
884
                os.system(cmd)
 
885
 
 
886
        ## GET NUMBER OF ROWS, COLUMNS AND LINE-LENGTH, USING WC
 
887
        wc = string.split(os.popen("wc "+tmpname).read())
 
888
        numlines = int(wc[0]) - lines_to_ignore
 
889
        tfp = open(tmpname)
 
890
        if lines_to_ignore <> 0:
 
891
                for i in range(lines_to_ignore):
 
892
                        junk = tfp.readline()
 
893
        numcols = len(string.split(tfp.readline())) #int(float(wc[1])/numlines)
 
894
        tfp.close()
 
895
 
 
896
        ## PREPARE INPUT SPACE
 
897
        a = N.zeros((numlines*numcols), type)
 
898
        block = 65536  # chunk to read, in bytes
 
899
        data = mmapfile.mmapfile(tmpname, '', 0)
 
900
        if lines_to_ignore <> 0 and sys.platform == 'win32':
 
901
                for i in range(lines_to_ignore):
 
902
                        junk = data.readline()
 
903
        i = 0
 
904
        d = ' '
 
905
        carryover = ''
 
906
        while len(d) <> 0:
 
907
                d = carryover + data.read(block)
 
908
                cutindex = string.rfind(d,'\n')
 
909
                carryover = d[cutindex+1:]
 
910
                d = d[:cutindex+1]
 
911
                d = map(intype,string.split(d))
 
912
                a[i:i+len(d)] = d
 
913
                i = i + len(d)
 
914
        end = time.time()
 
915
        print "%d sec" % round(end-start,2)
 
916
        data.close()
 
917
        os.remove(tmpname)
 
918
        return N.reshape(a,[numlines,numcols])
 
919
 
 
920
 
 
921
def find_dirs(sourcedir):
 
922
        '''Finds and returns all directories in sourcedir
 
923
 
 
924
Usage:   find_dirs(sourcedir)
 
925
Returns: list of directory names (potentially empty)
 
926
'''
 
927
        files = os.listdir(sourcedir)
 
928
        dirs = []
 
929
        for fname in files:
 
930
                if os.path.isdir(os.path.join(sourcedir,fname)):
 
931
                        dirs.append(fname)
 
932
        return dirs
 
933
 
 
934
 
 
935
# ALIASES ...
 
936
save = aput
 
937
 
 
938
 
 
939
 
 
940
def binget(fname,btype=None):
 
941
        '''
 
942
Loads a binary file from disk. Assumes associated hdr file is in same
 
943
location. You can force an unpacking type, or else it tries to figure
 
944
it out from the filename (4th-to-last character). Hence, readable file
 
945
formats are ...
 
946
 
 
947
1bin=int8, sbin=int16, ibin=int32, fbin=float32, dbin=float64, etc.
 
948
 
 
949
Usage:   binget(fname,btype=None)
 
950
Returns: data in file fname of type btype
 
951
'''
 
952
        file = open(fname,'rb')
 
953
        bdata = file.read()
 
954
        file.close()
 
955
 
 
956
        # if none given, assume character preceeding 'bin' is the unpacktype
 
957
        if not btype:
 
958
                btype = fname[-4]
 
959
        try:
 
960
                bdata = N.fromstring(bdata,btype)
 
961
        except:
 
962
                raise ValueError, "Bad unpacking type."
 
963
 
 
964
        # force the data on disk to be LittleEndian (for more efficient PC/Linux use)
 
965
        if not N.little_endian:
 
966
                bdata = bdata.byteswap()
 
967
 
 
968
        try:
 
969
                header = fname[:-3]+'hdr'
 
970
                vals = get(header,0)  # '0' means no missing-file warning msg
 
971
                print vals
 
972
                if type(vals[0]) == ListType:  # it's an extended header
 
973
                        xsize = int(vals[0][0])
 
974
                        ysize = int(vals[0][1])
 
975
                        numslices = int(vals[0][2])
 
976
                else:
 
977
                        bdata.shape = vals
 
978
        except:
 
979
                print "No (or bad) header file. Returning unshaped array."
 
980
        return N.array(bdata)
 
981
 
 
982
 
 
983
 
 
984
def binput(outarray,fname,packtype=None,writetype='wb'):
 
985
        '''
 
986
Unravels outarray and writes the data to a file, always in LittleEndian
 
987
format, along with a header file containing the original data shape. Default
 
988
is overwrite the destination file. Tries to figure out packtype from
 
989
4th-to-last character in filename. Thus, the routine understands these
 
990
file formats ...
 
991
 
 
992
1bin=int8, sbin=int16, ibin=int32, fbin=float32, dbin=float64, etc.
 
993
 
 
994
Usage:  binput(outarray,filename,packtype=None,writetype='wb')
 
995
'''
 
996
        if not packtype:
 
997
                packtype = fname[-4]
 
998
 
 
999
        # a speck of error checking
 
1000
        if packtype == N.int16 and outarray.dtype.char == 'f':
 
1001
                # check to see if there's data loss
 
1002
                if max(N.ravel(outarray)) > 32767 or min(N.ravel(outarray))<-32768:
 
1003
                        print "*** WARNING: CONVERTING FLOAT DATA TO OUT-OF RANGE INT16 DATA"
 
1004
        outdata = N.ravel(outarray).astype(packtype)
 
1005
 
 
1006
        # force the data on disk to be little_endian (for more efficient PC/Linux use)
 
1007
        if not N.little_endian:
 
1008
                outdata = outdata.byteswap()
 
1009
        outdata = outdata.tostring()
 
1010
        outfile = open(fname,writetype)
 
1011
        outfile.write(outdata)
 
1012
        outfile.close()
 
1013
 
 
1014
        # Now, write the header file
 
1015
        try:
 
1016
                suffixindex = string.rfind(fname,'.')
 
1017
                hdrname = fname[0:suffixindex+2]+'hdr'  # include .s or .f or .1 or whatever
 
1018
        except ValueError:
 
1019
                hdrname = fname
 
1020
        hdr = outarray.shape
 
1021
        print hdrname
 
1022
        outfile = open(hdrname,'w')
 
1023
        outfile.write(pstat.list2string(hdr))
 
1024
        outfile.close()
 
1025
        return None
 
1026
 
 
1027
def getafniparam(headfilename,paramname):
 
1028
        '''
 
1029
Loads in an AFNI header file, and returns the values of 'paramname'.
 
1030
 
 
1031
Usage:   getafniparam(headfile,paramname)
 
1032
Returns: appropriate "type" for params, or None if fails
 
1033
'''
 
1034
        if headfilename[-4:] == 'BRIK':  # if asked for BRIK, change it to HEAD
 
1035
                headfilename = headfilename[:-4]+'HEAD'
 
1036
        d = get(headfilename)
 
1037
        lines = open(headfilename,'r').readlines()
 
1038
        for i in range(len(lines)):
 
1039
                if string.find(lines[i],paramname) <> -1:
 
1040
                        count = d[i+1][-1]
 
1041
                        gotten = 0
 
1042
                        result = []
 
1043
                        for j in range(i+2,len(lines)):
 
1044
                                for k in range(len(d[j])):
 
1045
                                        if type(d[j][k]) == StringType:
 
1046
                                                result = d[j][k][1:count]
 
1047
                                                return result
 
1048
                                        else:
 
1049
                                                result.append(d[j][k])
 
1050
                                                gotten += 1
 
1051
                                if gotten == count:
 
1052
                                        break
 
1053
                        return result
 
1054
        return None
 
1055
 
 
1056
 
 
1057
def add2afnihistory(headfilename,newtext):
 
1058
        '''
 
1059
Adds 'newtext' to HISTORY_NOTE in afni file specified in headfilename.
 
1060
 
 
1061
Usage:   add2afnihistory(headfile,newtext)
 
1062
Returns: None
 
1063
'''
 
1064
        if headfilename[-4:] == 'BRIK':  # if asked for BRIK, change it to HEAD
 
1065
                headfilename = headfilename[:-4]+'HEAD'
 
1066
        d = get(headfilename)
 
1067
        lines = open(headfilename,'r').readlines()
 
1068
        for i in range(len(lines)):
 
1069
                if string.find(lines[i],'HISTORY_NOTE') <> -1:
 
1070
                        bytecount = d[i+1][-1]
 
1071
                        oldstr = lines[i+2][:-2]
 
1072
                        date = '[python:***  %s] ' %time.asctime()
 
1073
                        lines[i+2] = oldstr +'\\n' +date +newtext +'~\n'
 
1074
                        lines[i+1] = '  count = %s\n' %str(len(lines[i+2]))
 
1075
        f = open(headfilename,'w')
 
1076
        f.writelines(lines)
 
1077
        f.close()
 
1078
        return
 
1079
 
 
1080
 
 
1081
def array2afni(d,brikprefix,voltype=None,TR=2000,sliceorder='seqplus',geomparent=None,view=None,corrlength=1,briklabels=None,historytext=None):
 
1082
        '''
 
1083
Converts an array 'd' to an AFNI BRIK/HEAD combo via putbin and to3d. Tries to
 
1084
guess the AFNI volume type
 
1085
 
 
1086
voltype = {'-anat','-epan','-fim'}
 
1087
geomparent = filename of the afni BRIK file with the same geometry
 
1088
view = {'tlrc', 'acpc' or 'orig'}
 
1089
corrlength = # of images used in the (single-timeseries) correlation (for fico)
 
1090
briklabels = list of names (strings) to use for brick labels
 
1091
historytext = string to be appended to the history file, if any
 
1092
 
 
1093
Usage:   array2afni(d,brikprefix,voltype=None,TR=2000,
 
1094
                                        sliceorder='seqplus',geomparent=None,view=None,
 
1095
                                        corrlength=1,briklabels=None,historytext=None)
 
1096
Returns: None
 
1097
'''
 
1098
        # converts numpy typecode()s into appropriate strings for to3d command line
 
1099
        typecodemapping = {'c':'b',  # character
 
1100
                                           'B':'b',  # UnsignedInt8
 
1101
                                           'f':'f',  # float0, float8, float16, float32
 
1102
                                           'd':'f',  # float64
 
1103
                                           'b':'b',  # int0, int8
 
1104
                                           'h':'',   # int16
 
1105
                                           'i':'i',  # int32
 
1106
                                           'l':'i'}  # int
 
1107
 
 
1108
        # Verify that the data is proper size (3- or 4-D)
 
1109
        if len(d.shape) not in [3,4]:
 
1110
                raise ValueError, "A 3D or 4D array is required for array2afni() ... %s" %d.shape
 
1111
 
 
1112
        # Save out the array to a binary file, homebrew style
 
1113
        if d.dtype.char == N.float64:
 
1114
                outcode = 'f'
 
1115
        else:
 
1116
                outcode = d.dtype.char
 
1117
        tmpoutname = 'afnitmp.%sbin' % outcode
 
1118
        binput(d.astype(outcode),tmpoutname)
 
1119
        if not voltype:
 
1120
                if len(d.shape) == 3:  # either anatomy or functional
 
1121
                        if d.dtype.char in ['s','i','l']:  # if floats, assume functional
 
1122
                                voltype = '-anat'
 
1123
                        else:
 
1124
                                voltype = '-fim'
 
1125
                else:  # 4D dataset, must be anatomical timeseries (epan)
 
1126
                        voltype = '-anat'
 
1127
        if voltype[0] != '-':
 
1128
                voltype = '-'+voltype
 
1129
        if len(d.shape) == 3:  # either anatomy or functional
 
1130
                timepts = 1
 
1131
                slices = d.shape[0]
 
1132
                timestr = ''
 
1133
        elif len(d.shape) == 4:
 
1134
                if voltype=='-fico':
 
1135
                        timepts = 1
 
1136
                        d = N.reshape(d,[d.shape[0]*d.shape[1],d.shape[2],d.shape[3]])
 
1137
                        slices = d.shape[0]
 
1138
                        timestr = '-statpar %s 1 1 ' % corrlength
 
1139
                else:
 
1140
                        timepts = d.shape[0]
 
1141
                        slices = d.shape[1]
 
1142
                        timestr = '-time:zt %d %d %0.3f %s ' % (slices,timepts,TR,sliceorder)
 
1143
 
 
1144
        cmd = 'to3d %s -prefix %s -session . ' % (voltype, brikprefix)
 
1145
        if not view:
 
1146
                view = 'orig'
 
1147
        cmd += '-view %s ' % view
 
1148
        if geomparent:
 
1149
                cmd += '-geomparent %s ' % geomparent
 
1150
        cmd += timestr
 
1151
        cmd += '3D%s:0:0:%d:%d:%d:%s' % (typecodemapping[d.dtype.char],d.shape[-1],d.shape[-2],slices*timepts,tmpoutname)
 
1152
        print cmd
 
1153
        os.system(cmd)
 
1154
        os.remove(tmpoutname)
 
1155
        os.remove(tmpoutname[:-3]+'hdr')
 
1156
 
 
1157
        if len(d.shape)==4 and briklabels:
 
1158
                names = ''
 
1159
                for label in briklabels:
 
1160
                        names += str(label)+'~'
 
1161
                count = len(names)
 
1162
                appendstr = '''\n\ntype = string-attribute
 
1163
name = BRICK_LABS
 
1164
count = %s
 
1165
'%s''' % (count, names)
 
1166
                f = open('%s+%s.HEAD' %(brikprefix,view), 'a')
 
1167
                f.write(appendstr)
 
1168
                f.close()
 
1169
 
 
1170
                if historytext:
 
1171
                        add2afnihistory('%s+%s.HEAD'%(brikprefix,view),historytext)