← Back to branch summary

~ubuntu-branches/ubuntu/trusty/pynifti/trusty

« back to all changes in this revision

Viewing changes to bin/pynifti_pst

  • Committer: Bazaar Package Importer
  • Author(s): Michael Hanke, Michael Hanke, Yaroslav Halchenko
  • Date: 2007-09-17 08:32:59 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: james.westby@ubuntu.com-20070917083259-8eg0pofquqdn7w00
Tags: 0.20070917.1-1
http://bugs.debian.org/442175
[ Michael Hanke ]
* Bugfix: Can now update NIfTI header data when no filename is set
  (Closes: #442175).
* Unloading of image data without a filename set is no checked and prevented
  as it would damage data integrity and the image data could not be
  recovered.

[ Yaroslav Halchenko ]
* Added 'pixdim' property.

Show diffs side-by-side

added added

removed removed

Lines of Context:
bin/pynifti_pst
1
1
#!/usr/bin/python
2
 
### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ###
 
2
### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ##
3
3
#
4
4
#    Copyright (C) 2007 by
5
5
#    Michael Hanke <michael.hanke@gmail.com>
13
13
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
14
14
#    Lesser General Public License for more details.
15
15
#
16
 
### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ###
17
 
# SVN version control block - do not edit manually
18
 
# $Id$
19
 
### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ###
 
16
### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ##
20
17
 
21
18
 
22
19
import os, sys
23
20
from optparse import OptionParser
24
21
import nifti.utils
25
 
 
 
22
import numpy as np
26
23
 
27
24
def loadEVFile( filename ):
28
25
    try:
44
41
def checkCmdLine(args, opts, min_args):
45
42
    if len(args) < min_args:
46
43
        raise RuntimeError, \
47
 
            "%s: error: Incorrect number of arguments (supplied %i, needed min. %i)." \
48
 
                % (sys.argv[0], len(args), min_args)
 
44
            "%s: error: Incorrect number of arguments " \
 
45
            "(supplied %i, needed %i)." % (sys.argv[0], len(args), min_args)
49
46
 
50
47
 
51
48
def main():
52
49
    parser = OptionParser( \
53
 
               usage="%prog [options] <4d image> <outfile> <vol_id | filename> [ ... ]", \
54
 
               version="%prog 0.20070424", description="""\
 
50
        usage="%prog [options] <4dimage> <outfile> <vol_id | filename> [...]", \
 
51
        version="%prog 0.20070424", description="""\
55
52
%prog computes the peristimulus signal timecourse for all voxels in a volume at
56
53
once. Stimulus onsets can be specified as volume numbers or times (will be
57
54
converted into volume numbers using a supplied repetition time). Onsets can be
60
57
Empty lines are ignored. This enables %prog to use e.g. FSL's custom EV files.
61
58
If several files are specified, the read onsets are combined to a single onset
62
59
vector.
63
 
%prog write a 4d timeserie image as output. This image can e.g. be loaded into
 
60
%prog writes a 4d timeseries image as output. This image can e.g. be loaded into
64
61
FSLView to look at each voxels signal timecourse in a certain condition by
65
62
simply clicking on it.
66
63
""" )
67
64
 
68
65
    # define options
69
 
    parser.add_option('--verbose', action='store_true', dest='verbose', default=False,
70
 
                      help='print status messages')
71
 
    parser.add_option('-n', '--nvols', action='store', dest='nvols', default=10,
72
 
                      type="int", help="""\
73
 
set the length of the computed peristimulus signal timecourse (in volumes).
74
 
Default: 10 """ )
75
 
    parser.add_option('-t', '--times', action='store_true', dest='times', default=False,
76
 
                      help="""\
77
 
if supplied, the read values are treated as onset times and will be converted
 
66
    parser.add_option('--verbose', action='store_true', dest='verbose',
 
67
                      default=False, help='print status messages')
 
68
    parser.add_option('-n', '--nvols', action='store', dest='nvols',
 
69
                      default=10, type="int", help="""\
 
70
Set the length of the computed peristimulus signal timecourse (in volumes).
 
71
Default: 10
 
72
""" )
 
73
    parser.add_option('-t', '--times', action='store_true', dest='times',
 
74
                      default=False, help="""\
 
75
If supplied, the read values are treated as onset times and will be converted
78
76
to volume numbers. For each onset the volume that is closest in time will be
79
77
selected. Volumes are assumed to be recorded exactly (and completely) after
80
78
tr/2, e.g. if 'tr' is 2 secs the first volume is recorded at exactly one
82
80
conversion. """ )
83
81
    parser.add_option('--tr', action='store', dest='tr', default=None,
84
82
                      type="float", help="""\
85
 
repetition time of the 4d image (temporal difference of two successive
 
83
Repetition time of the 4d image (temporal difference of two successive
86
84
volumes). This can be used to override the setting in the 4d image. The
87
85
repetition time is necessary to convert onset times into volume numbers.
88
86
If the '--times' option is not set this value has no effect. Please note
89
87
that repetitions time and the supplied onsets have to be in the same unit.
 
88
 
 
89
Please note, that if --times is given the tr has to be specified in the
 
90
same unit as the read onset times.
90
91
""" )
91
92
    parser.add_option('--offset', dest='offset', action='store', default=0.0,
92
93
                      type='float', help="""\
93
 
constant offset applied to the onsets times when converting them to volume
 
94
Constant offset applied to the onsets times when converting them to volume
94
95
numbers. Without setting '--times' this option has no effect'.
95
96
""" )
96
 
    parser.add_option('-p', '--percsigchg', action='store_true', dest='perc_sig_chg', default=False,
97
 
                      help="""\
98
 
convert the computed timecourse to percent signal change relative to the first
99
 
(onset) volume. Default: False
 
97
    parser.add_option('-p', '--percsigchg', action='store_true',
 
98
                      dest='perc_sig_chg', default=False, help="""\
 
99
Convert the computed timecourse to percent signal change relative to the first
 
100
(onset) volume. This might not be meaningful when --operation is set to
 
101
something different than 'mean'. Please note, that the shape of the computes
 
102
timeseries heavily depends on the first average volume. It might be more
 
103
meaningful to use a real baseline condition as origin. However, this is not
 
104
supported yet. Default: False
 
105
""" )
 
106
    parser.add_option('--printvoxel', action='store', dest='printvoxel',
 
107
                      default=None, help="""\
 
108
Print the peristimulus timeseries of a single voxel for all onsets separately.
 
109
This will print a matrix (onsets x time), where the number of columns is
 
110
identical to the value of --nvols and the number of rows corresponds to the
 
111
number of specified onsets. (e.g. 'z,y,x')
 
112
""" )
 
113
    parser.add_option('--operation', action='store', dest='operation',
 
114
                      default='mean', help="""\
 
115
Choose the math operation that is performed to compute the peristimulus
 
116
timeseries. By default this is the mean across all stimulations ('mean').
 
117
Other possibilities are the standard deviation ('std') and standard error
 
118
('sde').
100
119
""" )
101
120
    # parse options
102
121
    (opts, args) = parser.parse_args() # read sys.argv[1:] by default
112
131
    nimg = nifti.NiftiImage( args[0] )
113
132
 
114
133
    if not nimg.timepoints > 1:
115
 
        print "%s: error: Need 4d image as input. Supplied image only has one volume." \
 
134
        print "%s: error: Need 4d image as input. " \
 
135
              "Supplied image only has one volume." \
116
136
                % sys.argv[0]
117
137
        sys.exit( 1 )
118
138
 
 
139
    # determine the requested function for timeseries calculation
 
140
    if opts.operation == 'mean':
 
141
        pst_fx = np.mean
 
142
    elif opts.operation == 'std' or opts.operation == 'sde':
 
143
        pst_fx = np.std
 
144
    else:
 
145
        print "'%s' is not a supported operation." % opts.operation
 
146
        sys.exit(1)
 
147
 
119
148
    outfilename = args[1]
120
149
    if opts.verbose:
121
150
        print "Output will be written to '%s'" % outfilename
140
169
            try:
141
170
                onsets.append( float( src ) )
142
171
            except ValueError:
143
 
                print "%s: error: '%s' cannot be converted into a floating-point value" \
144
 
                        % (sys.argv[0], src)
 
172
                print "%s: error: '%s' cannot be converted into a " \
 
173
                      "floating-point value" % (sys.argv[0], src)
145
174
                sys.exit( 1 )
146
175
 
147
176
    if opts.times:
148
177
        if opts.verbose:
149
 
            print "Convert supplied onset times into volumes (tr: %f, offset: %f)" \
150
 
                    % (tr, opts.offset)
 
178
            print "Convert supplied onset times into volumes " \
 
179
                  "(tr: %f, offset: %f)" % (tr, opts.offset)
151
180
 
152
 
        onsetvols = nifti.utils.time2vol( onsets, tr, opts.offset, decimals = 0 ).astype('int')
 
181
        onsetvols = nifti.utils.time2vol( onsets,
 
182
                                          tr,
 
183
                                          opts.offset,
 
184
                                          decimals = 0 ).astype('int')
153
185
    else:
154
186
        if opts.verbose:
155
187
            print "Verify onset volume numbers"
156
188
        onsetvols = [ int(i) for i in onsets ]
157
189
 
158
190
    if opts.verbose:
159
 
        print "Selected volumes:\n" + ', '.join([ str(o) for o in onsetvols])
160
 
 
161
 
    print "Compute mean peristimulus signal timecourse (length: %i volumes)" \
162
 
            % opts.nvols
 
191
        print "Selected volumes (index starts at zero!):\n" + ', '.join([ str(o) for o in onsetvols])
 
192
 
 
193
    if opts.verbose:
 
194
        print "Compute mean peristimulus signal timecourse " \
 
195
              "(length: %i volumes)" % opts.nvols
 
196
 
 
197
    if opts.printvoxel:
 
198
        # get timeseries for each onset
 
199
        pst = nifti.utils.getPeristimulusTimeseries( nimg,
 
200
                                                     onsetvols,
 
201
                                                     opts.nvols,
 
202
                                                     tuple ).copy()
 
203
        # make matrix ( onset x time )
 
204
        voxel_pst = eval('pst[:,:,' + opts.printvoxel +'].T')
 
205
 
 
206
        for onset in voxel_pst:
 
207
            for t in onset:
 
208
                print t,
 
209
            print ''
163
210
 
164
211
 
165
212
    if opts.verbose:
166
213
        print "Compute peristimulus timeseries"
167
 
    pst = nifti.utils.getPeristimulusTimeseries(nimg, onsetvols, opts.nvols).copy()
 
214
 
 
215
    pst = nifti.utils.getPeristimulusTimeseries( nimg,
 
216
                                                 onsetvols,
 
217
                                                 opts.nvols,
 
218
                                                 pst_fx ).copy()
 
219
 
 
220
    # divide by srqt of number of stimulations if standard error is requested
 
221
    if opts.operation == 'sde':
 
222
        pst /= np.sqrt( len(onsetvols) )
168
223
 
169
224
    if opts.perc_sig_chg:
170
225
        if opts.verbose:
171
226
            print "Convert to percent signal change"
172
 
        baseline = pst[0].copy()
173
 
        pst -= baseline
174
 
        pst /= baseline
 
227
        baseline = nifti.utils.getPeristimulusTimeseries( nimg,
 
228
                                                          onsetvols,
 
229
                                                          1,
 
230
                                                          np.mean )[0].copy()
 
231
 
 
232
        if opts.operation == 'mean':
 
233
            pst -= baseline
 
234
 
 
235
        # only divide if baseline is different from zero
 
236
        pst[:,baseline != 0] /= baseline[baseline != 0]
175
237
        pst *= 100.0
176
238
 
177
239
    if opts.verbose: