Skip to content

Instantly share code, notes, and snippets.

@ptone
Created June 12, 2013 01:16
Show Gist options
  • Save ptone/5762215 to your computer and use it in GitHub Desktop.
Save ptone/5762215 to your computer and use it in GitHub Desktop.
a heavily annotated python script to breakout and combine channels of LSM confocal images
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)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment