This file is indexed.

/usr/bin/pynifti_pst is in python-nifti 0.20100607.1-4.

This file is owned by root:root, with mode 0o755.

The actual contents of the file can be viewed below.

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
#!/usr/bin/python2.7
### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ##
#
#    Copyright (C) 2007-2008 by
#    Michael Hanke <michael.hanke@gmail.com>
#
#    This is free software; you can redistribute it and/or
#    modify it under the terms of the MIT License.
#
#    This package is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the COPYING
#    file that comes with this package for more details.
#
### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ##


import os, sys
from optparse import OptionParser
import nifti.utils
import numpy as N

def loadEVFile( filename ):
    try:
        file = open( filename )
    except IOError:
        raise ValueError, "Cannot open EV file '%s'" % filename

    onsets = []

    # split each line and ignore empty lines
    for line in file:
        e = line.split()
        if len(e):
            onsets.append( float( e[0] ) )

    return onsets


def checkCmdLine(args, opts, min_args):
    if len(args) < min_args:
        raise RuntimeError, \
            "%s: error: Incorrect number of arguments " \
            "(supplied %i, needed %i)." % (sys.argv[0], len(args), min_args)


def main():
    parser = OptionParser( \
        usage="%prog [options] <4dimage> <outfile> <vol_id | filename> [...]", \
        version="%prog " + nifti.__version__, description="""\
%prog computes the peristimulus signal timecourse for all voxels in a volume at
once. Stimulus onsets can be specified as volume numbers or times (will be
converted into volume numbers using a supplied repetition time). Onsets can be
specified directly on the command line, but can also be read from (multiple)
files. Such file are assumed to list one onset per line (as the first value).
Empty lines are ignored. This enables %prog to use e.g. FSL's custom EV files.
If several files are specified, the read onsets are combined to a single onset
vector.
%prog writes a 4d timeseries image as output. This image can e.g. be loaded into
FSLView to look at each voxels signal timecourse in a certain condition by
simply clicking on it.
""" )

    # define options
    parser.add_option('--verbose', action='store_true', dest='verbose',
                      default=False, help='print status messages')
    parser.add_option('-n', '--nvols', action='store', dest='nvols',
                      default=10, type="int", help="""\
Set the length of the computed peristimulus signal timecourse (in volumes).
Default: 10
""" )
    parser.add_option('-t', '--times', action='store_true', dest='times',
                      default=False, help="""\
If supplied, the read values are treated as onset times and will be converted
to volume numbers. For each onset the volume that is closest in time will be
selected. Volumes are assumed to be recorded exactly (and completely) after
tr/2, e.g. if 'tr' is 2 secs the first volume is recorded at exactly one
second. Please see the --tr and --offset options to learn how to adjust the
conversion. """ )
    parser.add_option('--tr', action='store', dest='tr', default=None,
                      type="float", help="""\
Repetition time of the 4d image (temporal difference of two successive
volumes). This can be used to override the setting in the 4d image. The
repetition time is necessary to convert onset times into volume numbers.
If the '--times' option is not set this value has no effect. Please note
that repetitions time and the supplied onsets have to be in the same unit.

Please note, that if --times is given the tr has to be specified in the
same unit as the read onset times.
""" )
    parser.add_option('--offset', dest='offset', action='store', default=0.0,
                      type='float', help="""\
Constant offset applied to the onsets times when converting them to volume
numbers. Without setting '--times' this option has no effect'.
""" )
    parser.add_option('-p', '--percsigchg', action='store_true',
                      dest='perc_sig_chg', default=False, help="""\
Convert the computed timecourse to percent signal change relative to the first
(onset) volume. This might not be meaningful when --operation is set to
something different than 'mean'. Please note, that the shape of the computes
timeseries heavily depends on the first average volume. It might be more
meaningful to use a real baseline condition as origin. However, this is not
supported yet. Default: False
""" )
    parser.add_option('--printvoxel', action='store', dest='printvoxel',
                      default=None, help="""\
Print the peristimulus timeseries of a single voxel for all onsets separately.
This will print a matrix (onsets x time), where the number of columns is
identical to the value of --nvols and the number of rows corresponds to the
number of specified onsets. (e.g. 'z,y,x')
""" )
    parser.add_option('--operation', action='store', dest='operation',
                      default='mean', help="""\
Choose the math operation that is performed to compute the peristimulus
timeseries. By default this is the mean across all stimulations ('mean').
Other possibilities are the standard deviation ('std') and standard error
('sde').
""" )
    # parse options
    (opts, args) = parser.parse_args() # read sys.argv[1:] by default

    try:
        checkCmdLine( args, opts, 3 )
    except:
        print sys.exc_info()[1]
        sys.exit( 1 )

    if opts.verbose:
        print "Read 4d image '%s'" % args[0]
    nimg = nifti.NiftiImage( args[0] )

    if not nimg.timepoints > 1:
        print "%s: error: Need 4d image as input. " \
              "Supplied image only has one volume." \
                % sys.argv[0]
        sys.exit( 1 )

    # determine the requested function for timeseries calculation
    if opts.operation == 'mean':
        pst_fx = N.mean
    elif opts.operation == 'std' or opts.operation == 'sde':
        pst_fx = N.std
    else:
        print "'%s' is not a supported operation." % opts.operation
        sys.exit(1)
 
    outfilename = args[1]
    if opts.verbose:
        print "Output will be written to '%s'" % outfilename

    if opts.tr:
        tr = opts.tr
        if opts.verbose:
            print "Using provide repetition time: %f" % tr
    else:
        tr = nimg.rtime
        if opts.verbose:
            print "Using repetition from 4d image: %f" % tr

    onsets = []

    for src in args[2:]:
        if os.path.isfile( src ):
            if opts.verbose:
                print "Reading values from '%s'" % src
            onsets += loadEVFile( src )
        else:
            try:
                onsets.append( float( src ) )
            except ValueError:
                print "%s: error: '%s' cannot be converted into a " \
                      "floating-point value" % (sys.argv[0], src)
                sys.exit( 1 )

    if opts.times:
        if opts.verbose:
            print "Convert supplied onset times into volumes " \
                  "(tr: %f, offset: %f)" % (tr, opts.offset)

        onsetvols = nifti.utils.time2vol( onsets,
                                          tr,
                                          opts.offset,
                                          decimals = 0 ).astype('int')
    else:
        if opts.verbose:
            print "Verify onset volume numbers"
        onsetvols = [ int(i) for i in onsets ]

    if opts.verbose:
        print "Selected volumes (index starts at zero!):\n" + ', '.join([ str(o) for o in onsetvols])

    if opts.verbose:
        print "Compute mean peristimulus signal timecourse " \
              "(length: %i volumes)" % opts.nvols

    if opts.printvoxel:
        # get timeseries for each onset
        pst = nifti.utils.getPeristimulusTimeseries( nimg,
                                                     onsetvols,
                                                     opts.nvols,
                                                     tuple ).copy()
        # make matrix ( onset x time )
        voxel_pst = eval('pst[:,:,' + opts.printvoxel +'].T')

        for onset in voxel_pst:
            for t in onset:
                print t,
            print ''


    if opts.verbose:
        print "Compute peristimulus timeseries"

    pst = nifti.utils.getPeristimulusTimeseries( nimg,
                                                 onsetvols,
                                                 opts.nvols,
                                                 pst_fx ).copy()

    # divide by srqt of number of stimulations if standard error is requested
    if opts.operation == 'sde':
        pst /= N.sqrt( len(onsetvols) )

    if opts.perc_sig_chg:
        if opts.verbose:
            print "Convert to percent signal change"
        baseline = nifti.utils.getPeristimulusTimeseries( nimg,
                                                          onsetvols,
                                                          1,
                                                          N.mean )[0].copy()

        if opts.operation == 'mean':
            pst -= baseline

        # only divide if baseline is different from zero
        pst[:,baseline != 0] /= baseline[baseline != 0]
        pst *= 100.0

    if opts.verbose:
        print "Write output"
    onimg = nifti.NiftiImage(pst,nimg.header)
    onimg.save( outfilename )

    if opts.verbose:
        print "Done"


if __name__ == "__main__":
    main()