/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()
|