a heavily annotated python script to breakout and combine channels of LSM confocal images

convert-lsm.py

import sys
import os
import scipy
from numpy import array, dstack
from pylsm import lsmreader

# all of the above imports pull into our program features we use of other
# modules. Ideally anything not defined in our program can be seen above

# sys.argv contains a list of all the files passed to our script when it is called

if len(sys.argv) == 1:
    # by default, it always contains the name of the current script file as the
    # first item.  any additional items are files passed to our script as
    # arguments. If there is only one item, no files were given, so print some
    # simple help
    print "no files given, use as convert-lsm.py file [file2 filen...]"
    # sys.exit will stop the script immediately
    sys.exit(1)

# --------------- creating a blank layer
# We need to make a blank, all black, layer for the RGB channels we want blank
#
# this is done right at the beginning, since all blank layers we need later
# are the same, so there is no reason to keep recreating it inside any of our loops
# below - you can skip reading this now, until it is needed later, and it will
# be more clear.

# create a 1 dimensional array of the right length - this is done by taking
# a single value as a simple list '[0]' and multiplying it by the total number
# pixels in your image
blank_layer = array([0] * 1024 * 1024)

# if we were to look at this interactively, we would see it is a single long list
# In [3]: blank_layer.shape
# Out[3]: (1048576,)

# the reshape command turns it into a 2 dimensional array
# note that the reshape command does not reshape "in place" but returns a new
# array that is reshaped - so we reassign it to "blank_layer" since we don't care
# about keeping the 1D layer around at all.
#
blank_layer = blank_layer.reshape(1024, 1024)
# --------------- done creating a blank layer

# the [1:] syntax is called slicing, and gives us everything from the second item on
# which is all the files passed to us. This for-loop is repeated for each file
# given - you may choose to call this with *.lsm in a folder
for filename in sys.argv[1:]:

    # gets a special image object
    lsmimage = lsmreader.Lsmimage(filename)

    # open the image to read the data
    lsmimage.open()

    # here we create a list to hold each channel's data
    channel_data = []

    # grab each channels data as a 2 dimmensional array (1024, 1024)
    for current_channel in [0, 1, 2, 3]:  # This can also be done as for x in range(4)
        image_data = lsmimage.get_image(channel=current_channel)
        # add the data to our list
        channel_data.append(image_data)

    # after the 'for loop', channel_data is a list of data arrays:
    #   [<channel 0 data>, <channel 1 data>, ... etc]

    # since we are done with the lsm file image, we close it to clean up some
    # memory:
    lsmimage.close()

    # The task now is to combine these in various ways
    # we need to find a technique of saving data into image files - here is how I 
    # figured out how to do it:
    # googled: "scipy save new image"
    # http://stackoverflow.com/questions/902761/saving-a-numpy-array-as-an-image
    #
    # then looking at docs for this imsave function:
    # http://docs.scipy.org/doc/scipy/reference/generated/scipy.misc.imsave.html
    #
    # we see that if the data is 2 dimensional, it will save as greyscale, but
    # if it is 3 dimensional, it will save as RGB
    #
    # so we need to take our channels and combine them into a 3D array
    #
    # I needed to learn more about working with numpy arrays so...
    # googled: working with numpy arrays
    # http://wiki.scipy.org/Tentative_NumPy_Tutorial#head-6a1bc005bd80e1b19f812e1e64e0d25d50f99fe2
    # this tought me about how to create the blank layer above, and stack
    # three layers to make the RGB 3D array

    # now we need to go through each of our combinations
    # to get one of our channels, we can access it like:
    # channel_data[0] will get us channel 0 (which is the first)
    # channels are 0, 1, 2, 3 since python is 0-based index

    # we want 0+1, 0+2, 0+3
    # we only have RGB to work with
    # so in this case I'm always going to use channel 0 for red, and the other as green

    for other in [1, 2, 3]:
        # dsatck is a numpy command to depth stack arrays on a new axis
        # it is what we use to take our 2d layers, and create a 3d array
        # we reuse the blank layer created above as needed
        rgb_data = dstack((channel_data[0], channel_data[other], blank_layer))
        # we want to create a new filename
        fname_without_extension, extension = os.path.splitext(filename)
        # we use a special builtin string formatting template like thing in python
        # to replace {} with other values
        # we add 1 to other, to convert back to the 1 based index system which I'm
        # guessing is what you normally think of for lsm channels
        new_name = "{}-1-{}.tiff".format(fname_without_extension, other + 1)
        # now that we have the new name, and the RGB data, we save out the file
        # scipy automatically creates the right format based on the extension used
        scipy.misc.imsave(new_name, rgb_data)