MINC/Tutorials/PythonSliceLoop

A slice loop example using pyminc

When you open a volume in pyminc, none of the data is read until you access the .data variable. This means that you can open lots of volumes at the same time, and then selectively get a hyperslab out in order to conserve memory. The example below reads a bunch of volumes to find the mode (most common voxel) at every point - a simplistic voxel voting scheme for segmentation. A few comments are in the text, more explanation follows:

#!/usr/bin/env python

# all the imports go here
from pyminc.volumes.factory import *
from numpy import *
from scipy.stats import *
from optparse import OptionParser

if __name__ == "__main__":

    # some argument handling
    usage = "Usage text"
    description = "Description text"
    
    parser = OptionParser(usage=usage, description=description)
    parser.add_option("--clobber", dest="clobber",
                      help="clobber output file",
                      type="string")

    (options, args) = parser.parse_args()

    if len(args) < 4:
        parser.error("Incorrect number of arguments")

    outfilename = args[-1]
    propfilename = args[-2]
    # clobber check should go here
    
    # create an array of volume handles
    volhandles = []

    nfiles = len(args)-2
    for i in range( nfiles ):
        volhandles.append(volumeFromFile(args[i], dtype='ubyte'))

    outfile = volumeFromInstance(volhandles[0], outfilename)
    proportion = volumeFromInstance(volhandles[0], propfilename)

    # hold one slice from all volumes
    sliceArray = zeros( (nfiles,
                         volhandles[0].sizes[1],
                         volhandles[0].sizes[2]))
    
    # now loop over the slices to get on slice at a time from all volumes - then take the mode                
    for i in range(volhandles[0].sizes[0]):
        print "SLICE: %i" % i
        for j in range(nfiles):
            t = volhandles[j].getHyperslab((i,0,0),
                                           (1,volhandles[0].sizes[1],
                                            volhandles[0].sizes[2]))
            t.shape = (volhandles[0].sizes[1], volhandles[0].sizes[2])
            sliceArray[j::] = t
            
        # take the mode
        m = mode(sliceArray)
        # and write the mode and its proportion to the output files
        outfile.data[i::] = m[0]
        proportion.data[i::] = m[1] / nfiles

    proportion.writeFile()
    proportion.closeVolume()
    outfile.writeFile()
    outfile.closeVolume()

The key steps in this slice loop is to create an array of all the MINC volume handles, an an array which holds one slice for all volumes, and then to loop over all slices. At every slice the array is refilled with the current slice from each volume, the stat (mode in this case) is computed, and the results inserted into the output volumes.