python-nifti 0.20100607.1-4 / usr / bin / pynifti_pst

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