This file is indexed.

/usr/bin/mvpa-prep-fmri is in python-mvpa 0.4.8-1.

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
#!/usr/bin/python2.7
# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
# vi: set ft=python sts=4 ts=4 sw=4 et:
### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ##
#
#   See COPYING file distributed along with the PyMVPA package for the
#   copyright and license terms.
#
### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ##
"""Tiny tool to prepare a directory for a typical analysis of fMRI data with
PyMVPA. Tools from the FSL suite will be used for preprocessing. It takes a 4D
fMRI timeseries as input and performs the following steps:

- extract an example volume

- perform motion correction using the example volume as reference

- conservative skull-stripping and brain mask generation

- masking of the motion-corrected timeseries with the brain mask

All results will be stored either in the current directory, or in a
subdirectory with the subject ID (if specified)."""

import sys
import os
from subprocess import call

import numpy as N

from mvpa.misc.cmdline import parser, opt
from mvpa.base import verbose, externals, error
import mvpa

if __debug__:
    from mvpa.base import debug


_EXFUNC_CONV_DICT = {'last' :  lambda x: x-1,
                     'first':  lambda x: 0,
                     'middle': lambda x: int(x/2)}
"""Dictionary to get exemplar volume given a literal string"""

def prepParser(parser):
    # use module docstring for help output
    parser.usage = "%s [OPTIONS] <fmri-data>\n\n" % sys.argv[0] + __doc__
    parser.version = "%prog " + mvpa.__version__

    parser.add_option(opt.verbose)
    parser.add_option(opt.help)

    parser.add_option("-s", "--subject-id",
                      action="store", type="string", dest="subj",
                      default=None,
                      help="Subject ID used as output path")

    parser.add_option("-e", "--example-func-vol",
                      action="store", type="string", dest="exfunc",
                      default='middle',
                      help="Volume (numeric ID or 'last', 'first', 'middle') "
                           "to be used as an example functional image. "
                           "Default: 10")

    parser.add_option("-m", "--mcflirt-options",
                      action="store", type="string", dest="mcflirt_opts",
                      default='',
                      help="Options for MCFLIRT. '-plots' is auto-added ")

    parser.add_option("-p", "--mcflirt-plots",
                      action="store_true", dest="mcflirt_plots",
                      help="Create a .pdf with plots of motion parameters")

    parser.add_option("-b", "--bet-options",
                      action="store", type="string", dest="bet_opts",
                      default='-f 0.3',
                      help="Options for BET. '-m' is auto-added. "
                           "Default: '-f 0.3' for a safe guess of the brain "
                           "outline")



def main():
    """
    """
    prepParser(parser)
    (options, infiles) = parser.parse_args()

    # late import of pynifti to be able to get help output without a big
    # external dep.
    externals.exists('nifti', raiseException=True)
    from nifti import NiftiImage

    if len(infiles) > 1 or not len(infiles):
        error("%s needs exactly one input fMRI image as argument. "
                "Got %s" % (sys.argv[0], str(infiles)))

    func_fname = infiles[0]

    # compressed or uncompressed? decide by input image
    # XXX maybe add override option
    if func_fname.lower().endswith('nii.gz'):
        nii_ext = '.nii.gz'
        verbose(2, "Output files will be compressed NIfTI images")
    else:
        nii_ext = '.nii'
        verbose(2, "Output files will be uncompressed NIfTI images")

    # determine output path
    if not options.subj is None:
        opath = options.subj
    else:
        opath = os.path.curdir

    if not os.path.exists(opath):
        verbose(1, "Create output directory '%s'" % opath)
        os.makedirs(opath)
    else:
        verbose(2, "Using output path '%s'" % opath)

    verbose(2, "Load image file from '%s'" % func_fname)
    func_nim = NiftiImage(func_fname, load=True)

    # process exfunc option
    exfunc = options.exfunc.lower()
    timepoints = func_nim.timepoints
    if exfunc in _EXFUNC_CONV_DICT.keys():
        exfuncid = _EXFUNC_CONV_DICT[exfunc](timepoints)
    else:
        try:
            exfuncid = int(exfunc)
        except ValueError, e:
            error("Failed to convert '%s' into numerical id of "
                  "volume." % (exfunc))

    if exfuncid >= timepoints or exfuncid < 0:
        error("Example functional volume id must be in the "
              "range 0 .. %d. Got %d." % (timepoints-1, exfuncid))

    verbose(2, "Extract volume %i as example volume" % exfuncid)
    ef_nim = NiftiImage(func_nim.data[exfuncid], func_nim.header)
    ef_nim.save(os.path.join(opath, 'example_func' + nii_ext))

    # close input file -- will operate on motion-corrected one later on
    del func_nim

    mcflirt_call = \
        ' '.join(
           ['mcflirt',
            '-in ' + func_fname,
            '-out ' + os.path.join(opath, 'func_mc'),
            '-reffile ' + os.path.join(opath, 'example_func'),
            '-verbose 0',
            '-plots',
            options.mcflirt_opts]).strip()

    verbose(2, "Perform motion correction ('%s')" % mcflirt_call)

    # run MCFLIRT (silence stderr; 5 being some random file descriptor)
    if call(mcflirt_call, shell=True, stderr=None):
        error("MCFLIRT failed to perform the motion correction.")

    if options.mcflirt_plots:
        verbose(2, "Plot motion parameters estimates")
        externals.exists('pylab', raiseException=True)
        mc = McFlirtParams(os.path.join(opath, 'func_mc.par'))

        for k, (title, fields, ylabel) in enumerate(
            (('Translation', ('x', 'y', 'z'), 'mm'),
             ('Rotation', ('rot1', 'rot2', 'rot3'), 'radians'))):
            P.subplot(211+k)
            P.title(title)
            P.plot([0, timepoints], [0, 0], '0.6')
            for dim in fields:
                P.plot(mc[dim], label=dim)
            P.legend()
            P.axis('tight')
            P.ylabel(ylabel)

        P.gcf().savefig(os.path.join(opath, 'func_mc.pdf'))

    bet_call = \
        ' '.join(
            ['bet',
             os.path.join(opath, 'example_func'),
             os.path.join(opath, 'example_func_brain'),
             '-m',
             options.bet_opts]).strip()

    verbose(2, "Determine brain mask in functional space ('%s')" % bet_call)

    # run BET (silence stderr; 5 being some random file descriptor)
    if call(bet_call, shell=True, stderr=None):
        error("BET failed to perform the skull stripping.")


    verbose(2, "Threshold image background using brain mask")
    mask_nim = NiftiImage(os.path.join(opath, 'example_func_brain_mask'))
    func_nim = NiftiImage(os.path.join(opath, 'func_mc'))
    # special case: single slice mask
    if len(mask_nim.extent) < 3:
        func_nim.data[:, N.asarray([mask_nim.data]) == 0] = 0
    else:
        func_nim.data[:, mask_nim.data == 0] = 0
    func_nim.save()


if __name__ == '__main__':
    main()