python-mvpa 0.4.7-2ubuntu1 / usr / bin / atlaslabeler

#!/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.
#
### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ##
"""Script to provide anatomical labels for the voxels, or their statistics """

import re, sys

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

if externals.exists('nifti', raiseException=True):
    from nifti import NiftiImage

if __debug__:
    from mvpa.base import debug

# from rumba.tools.misc import *
from mvpa.atlases.transformation import *

from mvpa.atlases import Atlas, ReferencesAtlas, FSLProbabilisticAtlas, \
     KNOWN_ATLASES, KNOWN_ATLAS_FAMILIES, XMLAtlasException

from optparse import OptionParser, Option

import numpy as N
#import numpy.linalg as la
# to read in transformation matrix

try:
    import psyco
    psyco.profile()
except:
    pass


def selectVoxelsFromVolumeIteratorNumPY(volFileName, lt=None, ut=None):
    """
    Generator which returns value + coordinates with values of non-0 entries
    from the `volFileName`

    :Returns:
     tuple with 0th entry value, the others are voxel coordinates
    More effective than previous loopy iteration since uses numpy's where
    function,    but for now is limited only to non-0 voxels selection
    """
    try:
        volFile = NiftiImage(volFileName)
    except:
        raise IOError("Cannot open image file %s" % volFileName)

    volData = volFile.data
    voxdim = volFile.voxdim
    if lt is None and ut is None:
        mask = volData != 0.0
    elif lt is None and ut is not None:
        mask = volData <= ut
    elif lt is not None and ut is None:
        mask = volData >= lt
    else:
        mask = N.logical_and(volData >= lt, volData <= ut)

    matchingVoxels = N.where(mask)
    # qf = volFile.qform
    # qfi = la.inv(qf)
    for ivoxel in xrange(len(matchingVoxels[0])):
        # in reverse order since numpy struct has order t, z, y, x
        voxel = ( matchingVoxels[2][ivoxel],
                  matchingVoxels[1][ivoxel],
                  matchingVoxels[0][ivoxel] )
        value = volData[ voxel[-1], voxel[-2], voxel[-3] ]
        yield (value, voxel[0], voxel[1], voxel[2], 0) #voxel[-4])


def parsedCoordinatesIterator(
    parseString="^\s*(?P<x>\S+)[ \t,](?P<y>\S+)[ \t,](?P<z>\S+)\s*$",
    inputStream=sys.stdin):
    """Iterator to provide coordinates/values parsed from the string stream,
    most often from the stdin
    """
    parser = re.compile(parseString)
    for line in inputStream.readlines():
        line = line.strip()
        match = parser.match(line)
        if not match:
            if __debug__:
                debug('ATL', "Line '%s' did not match '%s'"
                      % (line, parseString))
        else:
            r = match.groupdict()
            if r.has_key('v'): v = float(r['v'])
            else:              v = 0.0
            if r.has_key('t'): t = float(r['t'])
            else:              t = 0.0
            yield (v, float(r['x']), float(r['y']), float(r['z']), t)


# XXX helper to process labels... move me
def presentLabels(labels):
    if isinstance(labels, list):
        res = []
        for label in labels:
            # XXX warning -- some inconsistencies in atlas.py
            #     need refactoring
            s = label['label'] #.text
            if label.has_key('prob') and not options.createSummary:
                s += "(%d%%%%)" % label['prob']
            res += [s]
        if res == []:
            res = ['None']
        return '/'.join(res)
    else:
        if options.abbreviatedLabels:
            return labels['label'].abbr
        else:
            return labels['label'].text



#def processCmdLine():
parser.usage = "%s [OPTIONS] [input_file.nii.gz]" % sys.argv[0] + """
Examples:

%s -s -A talairach-dist -d 10 -R Closest\ Gray -l Structure,Brodmann\ area  -cC mask.nii.gz

produces a summary per each structure and brodmann area, for each voxel
looking within 10mm radius for the closest gray matter voxel.
""" % (sys.argv[0], )
parser.version = "%prog " + mvpa.__version__

# can't use due to conflict with -d (debug and distance)
#parser.option_groups = [opts.common]

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

parser.add_option("-a", "--atlas-file",
                  action="store", type="string", dest="atlasFile",
                  default=None,
                  help="Atlas file to use. Overrides --atlas-path and --atlas")

parser.add_option("--atlas-path",
                  action="store", type="string", dest="atlasPath",
                  default=None,
                  help=r"Path to the atlas files. '%(name)s' will be replaced"
                       " with the atlas name. See -A. Defaults depend on the"
                       " atlas family.")

parser.add_option("-A", "--atlas",
                  action="store", type="choice", dest="atlasName",
                  default="talairach", choices=KNOWN_ATLASES.keys(),
                  help="Atlas to use. Choices: %s"
                       % ', '.join(KNOWN_ATLASES.keys()))

parser.add_option("-i", "--input-coordinates-file",
                  action="store", type="string", dest="inputCoordFile",
                  default=None,
                  help="Fetch coordinates from ASCII file")

parser.add_option("-o", "--output-file",
                  action="store", type="string", dest="outputFile",
                  default=None,
                  help="Output file. Otherwise standard output")

parser.add_option("-d", "--max-distance",
                  action="store", type="float", dest="maxDistance",
                  default=0,
                  help="When working with reference/distance atlases, what"
                  " maximal distance to use to look for the voxel of interest")

parser.add_option("-T", "--transformation-file",
                  type="string", dest="transformationFile",
                  help="First transformation to apply to the data. Usually"+
                  " should be subject -> standard(MNI) transformation")

parser.add_option("-s", "--summary",
                  action="count", dest="createSummary", default=0,
                  help="Either to create a summary instead of dumping voxels."
                  " Use multiple -s for greater verbose summary")


parser.add_option("--ss", "--sort-summary-by",
                  type="choice", dest="sortSummaryBy", default="name",
                  choices=['name', 'count', 'a-p'],
                  help="How to sort summary entries. "
                  " a-p sorts anterior-posterior order")

parser.add_option("--dumpmap-file",
                  action="store", dest="dumpmapFile", default=None,
                  help="If original data is given as image file, dump indexes"
                  " per each treholded voxels into provided here output file")

parser.add_option("-l", "--levels",
                  type="string", dest="levels", default=None,
                  help="Indexes of levels which to print, or based on which "
                  "to create a summary (for a summary levels=4 is default). "
                  "To get listing of known for the atlas levels, use '-l list'")

parser.add_option("--mni2tal",
                  type="choice",
                  choices=["matthewbrett", "lancaster07fsl",
                           "lancaster07pooled", "meyerlindenberg98"],
                  dest="MNI2TalTransformation", default="matthewbrett",
                  help="Choose between available transformations from mni "
                  "2 talairach space")

parser.add_option("--thr", "--lthr", "--lower-threshold",
                  action="store", type="float", dest="lowerThreshold",
                  default=None,
                  help="Lower threshold for voxels to output")

parser.add_option("--uthr", "--upper-threshold",
                  action="store", type="float", dest="upperThreshold",
                  default=None,
                  help="Upper threshold for voxels to output")

parser.add_option("--abbr", "--abbreviated-labels",
                  action="store_true", dest="abbreviatedLabels",
                  help="Manipulate with abbreviations for labels instead of"
                  " full names, if the atlas has such")

# Parameters to be inline with older talairachlabel

parser.add_option("-c", "--tc", "--show-target-coord",
                  action="store_true", dest="showTargetCoordinates",
                  help="Show target coordinates")


parser.add_option("--tv", "--show-target-voxel",
                  action="store_true", dest="showTargetVoxel",
                  help="Show target coordinates")

parser.add_option("--rc", "--show-referenced-coord",
                  action="store_true", dest="showReferencedCoordinates",
                  help="Show referenced coordinates/distance in case if we are"
                  " working with reference atlas")

parser.add_option("-C", "--oc", "--show-orig-coord",
                  action="store_true", dest="showOriginalCoordinates",
                  help="Show original coordinates")

parser.add_option("-V", "--show-values",
                  action="store_true", dest="showValues",
                  help="Show values")


parser.add_option("-I", "--input-space",
                  action="store", type="string", dest="inputSpace",
                  default="MNI",
                  help="Space in which input volume/coordinates provided in. For instance Talairach/MNI")

parser.add_option("-F", "--forbid-direct-mapping",
                  action="store_true", dest="forbidDirectMapping",
                  default=False,
                  help="If volume is provided it first tries to do direct "
                  "mapping voxel-2-voxel if there is no transformation file "
                  "given. This option forbids such behavior and does "
                  "coordinates mapping anyway.")

parser.add_option("-t", "--talairach",
                  action="store_true", dest="coordInTalairachSpace",
                  default=False,
                  help="Coordinates are in talairach space (1x1x1mm)," +
                  " otherwise assumes in mni space (2x2x2mm)."
                  " Shortcut for '-I Talairach'")

parser.add_option("-H", "--half-voxel-correction",
                  action="store_false", dest="halfVoxelCorrection",
                  default=True,
                  help="Adjust coord by 0.5mm after transformation to " + \
                  "Tal space. Please use -H to turn such adjustment off")

parser.add_option("-r", "--relative-to-origin",
                  action="store_true", dest="coordRelativeToOrigin",
                  help="Coords are relative to the origin standard form" +
                  " ie in spatial units (mm), otherwise the default assumes" +
                  " raw voxel dimensions")

parser.add_option("--input-line-format",
                  action="store", type="string", dest="inputLineFormat",
                  default=r"^\s*(?P<x>\S+)[ \t,]+(?P<y>\S+)[ \t,]+(?P<z>\S+)\s*$",
                  help="Format of the input lines (if ASCII input is provided)")

# Specific atlas options
# TODO : group into options groups

# Reference atlas
parser.add_option("-R", "--reference",
                  action="store", type="string", dest="referenceLevel",
                  default="Closest Gray",
                  help="Which level to reference in the case of reference"
                  " atlas")

# Probabilistic atlases
parser.add_option("--prob-thr",
                  action="store", type="float", dest="probThr",
                  default=25.0,
                  help="At what probability (in %) to threshold in "
                  "probabilistic atlases (e.g. FSL)")

parser.add_option("--prob-strategy",
                  action="store", type="choice", dest="probStrategy",
                  choices=['all', 'max'], default='max',
                  help="What strategy to use for reporting. 'max' would report"
                  " single area (above threshold) with maximal probabilitity")


(options, infiles) = parser.parse_args()
#atlas.relativeToOrigin = options.coordRelativeToOrigin

if len(infiles)>1:
    print "We cannot handle multiple input files at once"
    sys.exit(1)

fileIn = None
coordT = None
niftiInput = None
# Setup coordinates read-in
#
# compatibility with older talairachlabel
if options.inputCoordFile:
    fileIn = file(options.inputCoordFile)
    coordsIterator = parsedCoordinatesIterator(options.inputLineFormat,
                                               fileIn)
# input is stdin
elif len(infiles)==0:
    coordsIterator = parsedCoordinatesIterator(options.inputLineFormat)
else:
    if len(infiles)>1:
        print "Just a single file should be provided at the command line"
        sys.exit(1)
    infile = infiles[0]
# got a volume/file to process
    try:
        if __debug__:
            debug('ATL', "Testing if 0th element in the list a volume")
        niftiInput = NiftiImage(infile)
        if __debug__:
            debug('ATL', "Yes it is")
        # if we got here -- it is a proper volume
        # XXX ask Michael to remove nasty warning message
        coordsIterator = selectVoxelsFromVolumeIteratorNumPY(
            infile, options.lowerThreshold, options.upperThreshold)
        assert(coordT is None)
        coordT = Linear(niftiInput.qform)
        # previous iterator returns space coordinates
        options.coordRelativeToOrigin = True
    except:
        if __debug__:
            debug('ATL', "No it is not")
        fileIn = file(infile)
        coordsIterator = parsedCoordinatesIterator(
            options.inputLineFormat, fileIn)


# Open and initialize atlas lookup
if options.atlasFile is None:
    if options.atlasPath is None:
        options.atlasPath = KNOWN_ATLASES[options.atlasName]
    options.atlasFile = options.atlasPath % ( {'name': options.atlasName} )

if not options.forbidDirectMapping \
       and niftiInput is not None and not options.transformationFile:
    akwargs = {
        'resolution': niftiInput.pixdim[0],
        'query_voxel': True }
    verbose(1, "Will attempt direct mapping from input voxels into atlas "
               "voxels at resolution %.2f" % akwargs['resolution'])

    atlas = Atlas(options.atlasFile, **akwargs)

    # verify that we got the same qforms in atlas and in the data file
    if atlas.space != options.inputSpace:
        verbose(0,
            "Cannot do direct mapping between input image in %s space and"
            " atlas in %s space. Use -I switch to override input space if"
            " it misspecified, or use -T to provide transformation. Trying"
            " to proceed" %(options.inputSpace, atlas.space), 1)
        atlas.query_voxel = False
    elif not (niftiInput.qform == atlas._image.qform).all():
        warning(
            "Cannot do direct mapping between files with different qforms."
            " Please provide original transformation (-T)."
            "\n Input qform:\n%s\n Atlas qform: \n%s"
            %(niftiInput.qform, atlas._image.qform), 1)
        # reset variables
        atlas.query_voxel = False
    else:
        coordT = None
else:
    atlas = Atlas(options.atlasFile)


if isinstance(atlas, ReferencesAtlas):
    options.referenceLevel = options.referenceLevel.replace('/', ' ')
    atlas.setReferenceLevel(options.referenceLevel)
    atlas.distance = options.maxDistance
else:
    options.showReferencedCoordinates = False

if isinstance(atlas, FSLProbabilisticAtlas):
    atlas.strategy = options.probStrategy
    atlas.thr = options.probThr

## If not in Talairach -- in MNI with voxel size 2x2x2
# Original talairachlabel assumed that if respective to origin -- voxels were
# scaled already.
#if options.coordInTalairachSpace:
#   voxelSizeOriginal = N.array([1, 1, 1])
#else:
#   voxelSizeOriginal = N.array([2, 2, 2])

if options.coordInTalairachSpace:
        options.inputSpace = "Talairach"

if not (options.inputSpace == atlas.space or
        (options.inputSpace in ["MNI", "Talairach"] and
         atlas.space == "Talairach")):
    raise XMLAtlasException("Unknown space '%s' which is not the same as atlas"
                            "space '%s' either" % ( inputSpace, atlas.space ))

if atlas.query_voxel:
    # we do direct mapping
    coordT = None
else:
    verbose(2, "Chaining needed transformations")
    # by default -- no transformation
    if options.transformationFile:
        externals.exists('scipy', raiseException=True)
        from scipy.io import read_array

        transfMatrix = read_array(options.transformationFile)
        coordT = Linear(transfMatrix, previous=coordT)
        verbose(2, "coordT got linear transformation from file %s" %
                   options.transformationFile)

    voxelOriginOriginal = None
    voxelSizeOriginal = None

    if not options.coordRelativeToOrigin:
        if options.inputSpace == "Talairach":
            # assume that atlas is in Talairach space already
            voxelOriginOriginal = atlas.origin
            voxelSizeOriginal = N.array([1, 1, 1])
        elif options.inputSpace == "MNI":
            # need to adjust for MNI origin as it was thought to be at
            # in terms of voxels
            #voxelOriginOriginal = N.array([46, 64, 37])
            voxelOriginOriginal = N.array([45, 63, 36])
            voxelSizeOriginal = N.array([2.0, 2.0, 2.0])
            warning("Assuming elderly sizes for MNI volumes with"
                       " origin %s and sizes %s" %\
                       ( `voxelOriginOriginal`, `voxelSizeOriginal`))


    if not (voxelOriginOriginal is None and voxelSizeOriginal is None):
        verbose(2, "Assigning origin adjusting transformation with"+\
                " origin=%s and voxelSize=%s" %\
                ( `voxelOriginOriginal`, `voxelSizeOriginal`))

        coordT = SpaceTransformation(origin=voxelOriginOriginal,
                                     voxelSize=voxelSizeOriginal,
                                     toRealSpace=True, previous=coordT)

    # besides adjusting for different origin we need to transform into
    # Talairach space
    if options.inputSpace == "MNI" and atlas.space == "Talairach":
        verbose(2, "Assigning transformation %s" %
                   options.MNI2TalTransformation)
        # What transformation to use
        coordT = {"matthewbrett": MNI2Tal_MatthewBrett,
                  "lancaster07fsl":  MNI2Tal_Lancaster07FSL,
                  "lancaster07pooled":  MNI2Tal_Lancaster07pooled,
                  "meyerlindenberg98":  MNI2Tal_MeyerLindenberg98,
                  "yohflirt": MNI2Tal_YOHflirt
                  }\
                  [options.MNI2TalTransformation](previous=coordT)

    if options.inputSpace == "MNI" and options.halfVoxelCorrection:
        originCorrection = N.array([0.5, 0.5, 0.5])
    else:
        # perform transformation any way to convert to voxel space (integers)
        originCorrection = None

    # To be closer to what original talairachlabel did -- add 0.5 to each coord
    coordT = SpaceTransformation(origin=originCorrection, voxelSize=None,
                                     toRealSpace=False, previous = coordT)

if options.createSummary:
    summary = {}
    if options.levels is None:
        options.levels = str(min(4, atlas.Nlevels-1))
if options.levels is None:
    options.levels = range(atlas.Nlevels)
elif isinstance(options.levels, basestring):
    if options.levels == 'list':
        print "Known levels and their indicies:\n" + atlas.levelsListing()
        sys.exit(0)
    slevels = options.levels.split(',')
    options.levels = []
    for level in slevels:
        try:
            int_level = int(level)
        except:
            if atlas.levels_dict.has_key(level):
                int_level = atlas.levels_dict[level].index
            else:
                raise RuntimeError(
                    "Unknown level '%s'. " % level +
                    "Known levels and their indicies:\n"
                    + atlas.levelsListing())
        options.levels += [int_level]
else:
    raise ValueError("Don't know how to handle list of levels %s."
                     "Example is '1,2,3'" % (options.levels,))

verbose(3, "Operating on following levels: %s" % options.levels)
# assign levels to the atlas
atlas.levels = options.levels

if options.outputFile:
    output = open(options.outputFile, 'w')
else:
    output = sys.stdout

# validity check
if options.dumpmapFile:
    if niftiInput is None:
        raise RuntimeError, "You asked to dump indexes into the volume, " \
              "but input wasn't a volume"
        sys.exit(1)
    ni_dump = NiftiImage(infile)
    ni_dump_data = N.zeros((len(options.levels),) + ni_dump.data.shape[:3])


# Read coordinates
numVoxels = 0
for c in coordsIterator:

    value, coord_orig, t = c[0], c[1:4], c[4]
    if __debug__:
        debug('ATL', "Obtained coord_orig=%s with value %s"
              % (repr(coord_orig), value))

    lt, ut = options.lowerThreshold, options.upperThreshold
    if lt is not None and value < lt:
        verbose(5, "Value %s is less than lower threshold %s, thus voxel "
                "is skipped" % (value, options.lowerThreshold))
        continue
    if ut is not None and value > ut:
        verbose(5, "Value %s is greater than upper threshold %s, thus voxel "
                "is skipped" % (value, options.upperThreshold))
        continue

    numVoxels += 1

    # Apply necessary transformations
    coord = coord_orig = N.array(coord_orig)

    if coordT:
        coord = coordT[ coord_orig ]

    # Query label
    voxel = atlas[ coord ]
    voxel['coord_orig'] = coord_orig
    voxel['value'] = value
    voxel['t'] = t
    if options.createSummary:
        summaryIndex = ""
        voxel_labels = voxel["labels"]
        for i,ind in enumerate(options.levels):
            voxel_label = voxel_labels[i]
            text = presentLabels(voxel_label)
            #if len(voxel_label):
            #   assert(voxel_label['index'] == ind)
            summaryIndex += text + " / "
        if not summary.has_key(summaryIndex):
            summary[summaryIndex] = {'values':[], 'max':value,
                                     'maxcoord':coord_orig}
            if voxel.has_key('voxel_referenced'):
                summary[summaryIndex]['distances'] = []
        summary_ = summary[summaryIndex]
        summary_['values'].append(value)
        if summary_['max'] < value:
            summary_['max'] = value
            summary_['maxcoord'] = coord_orig
        if voxel.has_key('voxel_referenced'):
            if voxel['voxel_referenced'] and voxel['distance']>=1e-3:
                verbose(5, 'Appending distance %e for voxel at %s'
                        % (voxel['distance'], voxel['coord_orig']))
                summary_['distances'].append(voxel['distance'])
    else:
        # Display while reading/processing
        first, out = True, ""

        if options.showValues:
            out += "%(value)5.2f "
        if options.showOriginalCoordinates:
            out += "%(coord_orig)s ->"
        if options.showReferencedCoordinates:
            out += " %(voxel_referenced)s=>%(distance).2f=>%(voxel_queried)s ->"
        if options.showTargetCoordinates:
            out += " %(coord_queried)s: "
            #out += "(%d,%d,%d): " % tuple(map(lambda x:int(round(x)),coord))
        if options.showTargetVoxel:
            out += " %(voxel_queried)s ->"

        if options.levels is None:
            options.levels = range(len(voxel['labels']))

        labels = [presentLabels(voxel['labels'][i]) for i in options.levels]
        out += ','.join(labels)
        #if options.abbreviatedLabels:
        #   out += ','.join([l.abbr for l in labels])
        #else:
        #   out += ','.join([l.text for l in labels])
        #try:
        output.write(out % voxel + "\n")
        #except:
        #    import pydb
        #    pydb.debugger()

    if options.dumpmapFile:
        try:
            ni_dump_data[:, coord_orig[-1], coord_orig[-2], coord_orig[-3]] = \
              [voxel['labels'][i]['label'].index
               for i,ind in enumerate(options.levels)]
        except Exception, e:
            import pydb
            pydb.debugger()

# if we opened any file -- close it
if fileIn:
    fileIn.close()

if options.dumpmapFile:
    ni_dump = NiftiImage(ni_dump_data, ni_dump.header)
    ni_dump.save(options.dumpmapFile)


def statistics(values):
    N_ = len(values)
    if N_==0:
        return 0, None, None, None, None, ""
    mean = N.mean(values)
    std = N.std(values)
    minv = N.min(values)
    maxv = N.max(values)
    ssummary = "[%3.2f : %3.2f] %3.2f+-%3.2f" % (minv, maxv, mean, std)
    return N_, mean, std, minv, maxv, ssummary


def generateSummary(summary, output):
    """Output the summary
    """
    # Sort either by the name (then ascending) or by the number of
    # elements (then descending)
    sort_keys = [(k, len(v['values']), v['maxcoord'][1])
                 for k,v in summary.iteritems()]
    sort_index, sort_reverse = {
        'name' : (0, False),
        'count': (1, True),
        'a-p': (2, True)}[options.sortSummaryBy]
    sort_keys.sort(cmp=lambda x,y: cmp(x[sort_index], y[sort_index]),
                   reverse=sort_reverse)
    # and here are the keys
    keys = [x[0] for x in sort_keys]
    maxkeylength = max (map(len, keys))

    # may be I should have simply made a counter ;-)
    total = sum(map(lambda x:len(x['values']), summary.values()))
    for index in keys:
        summary_ = summary[index]
        values = summary_['values']
        N, mean, std, minv, maxv, ssummary = statistics(values)
        # print "N=", N
        msg = "%%%ds:" % maxkeylength
        output.write(msg % index)
        output.write("%4d/%4.1f%% items" \
                     % (N, 100.0*N/total ))

        if options.createSummary>1:
            output.write(" %s" % ssummary)

        if options.createSummary>2:
            output.write(" max at %s" % summary_['maxcoord'])

        if options.createSummary>3 and summary_.has_key('distances'):
            # if we got statistics over referenced voxels
            Nr, mean, std, minv, maxv, ssummary = \
                statistics(summary_['distances'])
            Nr = len(summary_['distances'])
            # print "N=", N, " Nr=", Nr
            output.write(" Referenced: %d/%d%% Distances: %s" \
                         % (Nr, int(Nr*100.0 / N), \
                            ssummary))
        output.write("\n")
        # output might fail to flush, like in the case with broken pipe
        # -- imho that is not a big deal, ie not worth scaring the user
        try:
            output.flush()
        except IOError:
            pass
    output.write("-----\n")
    output.write("TOTAL: %d items\n" % total)

if options.createSummary:
    if numVoxels == 0:
        verbose(1, "No matching voxels were found.")
    else:
        generateSummary(summary, output)

if options.outputFile:
    output.close()