A Fiji Scripting Tutorial
|
Most of what you want to do with an image exists in Fiji. This tutorial will provide you with the general idea of how Fiji works: how are its capabilities organized, and how can they be composed into a program. To learn about Fiji, we'll start the hard way: by programming. This tutorial will teach you both python and Fiji. |
|
|
Index
|
Tutorial created by Albert Cardona. Zurich, 2010-11-10. All source code is under the Public Domain. Remember: Fiji is just ImageJ (batteries included). See also:
|
|
1. Getting startedOpen the Script Editor by choosing "File - New - Script". |
|
|
|
Alternatively, use the Command finder: Push 'l' (letter L) and then start typing "scri". You will see a list of Fiji commands, getting shorter the more letters you type. When the "Script Editor" command is visible, push the up arrow to go to it, and then push return to launch it. (Or double-click on it.) The Command Finder is useful for invoking any Fiji command. |
|
|
2. Your first Fiji scriptWe'll need an image to work on: please open any image. |
|
|
|
This tutorial will use the programming language Python 2.5. We start by telling the opened Script Editor what language you want to write the script on: choose "Language - Python". |
|
|
|
Type in what you see on the image to the right into the Script Editor, and then push "Run", or choose "Run - Run", or control+R (command+R in MacOSX). Line by line:
So what is "imp"? "imp" is a commonly used name to refer to an instance of an ImagePlus. The ImagePlus is one of ImageJ's abstractions to represent an image. Every image that you open in ImageJ is an instance of ImagePlus. |
|
|
|
Saving an image with a file dialog The first action we'll do on our image is to save it. To do that, you could call "File - Save" from the menus. |
|
|
|
Saving an image directly to a file The point of running a script is to avoid human interaction. Notice that the '#' sign defines comments. In python, any text after a '#' is not executed. |
from ij import IJ
from ij.io import FileSaver
imp = IJ.getImage()
fs = FileSaver(imp)
# A known folder to store the image at:
folder = "/home/albert/Desktop/t2/fiji-tutorial"
filepath = folder + "/" + "boats.tif"
fs.saveAsTiff(filepath):
|
|
|
Saving an image ... checking first if it's a good idea. The FileSaver will overwrite whatever file exists at the file path that you give it. That is not always a good idea! Here, we write the same code but checking first:
And finally, if all expected preconditions hold, then we place the call to saveAsTiff. This script introduced three new programming items:
|
from ij import IJ
from ij.io import FileSaver
from os import path
imp = IJ.getImage()
fs = FileSaver(imp)
# A known folder to store the image at:
folder = "/home/albert/Desktop/t2/fiji-tutorial"
# Test if the folder exists before attempting to save the image:
if path.exists(folder) and path.isdir(folder):
print "folder exists:", folder
filepath = path.join(folder, "boats.tif") # Operating System-specific
if path.exists(filepath):
print "File exists! Not saving the image, would overwrite a file!"
elif fs.saveAsTiff(filepath):
print "File saved successfully at ", filepath
else:
print "Folder does not exist or it's not a folder!"
|
|
3. Inspecting properties and pixels of an imageAn image in ImageJ or Fiji is, internally, an instance of ImagePlus. In python, accessing fields of an instance is straightforward: just add a dot '.' between the variable "imp" and the field "title" to access. In the Fiji API documentation, if you don't see a specific field like width in a particular class, but there is a getWidth method, then from python they are one and the same. The image type Notice how we created a dictionary to hold key/value pairs: of the image type versus a text representation of that type. This dictionary (also called map or table in other programming languages) then lets us ask it for a specific image type (such as ImagePlus.GRAY8), and we get back the corresponding text, such as "8-bit". You may have realized by now that the ImagePlus.getType() (or what is the same in python: "imp.type") returns us any of the controled values of image type that an ImagePlus instance can take. These values are GRAY8, GRAY16, GRAY32, COLOR_RGB, and COLOR_256. What is the image type? It's the kind of pixel data that the image holds. It could be numbers from 0 to 255 (what fits in an 8-bit range), or from 0 to 65536 (values that fit in a 16-bit range), or could be three channels of 8-bit values (an RGB image), or floating-point values (32-bit). The COLOR_256 indicates an 8-bit image that has an associated look-up table: each pixel value does not represent an intensity, but rather it's associated with a color. The table of values versus colors is limited to 256, and hence these images may not look very well. For image processing, you should avoid COLOR_256 images (also known as "8-bit color" images). These images are meant for display in the web in ".gif" format, but have been superseeded by JPEG or PNG. The GRAY_8 ("8-bit"), GRAY_16 ("16-bit") and GRAY_32 ("32-bit") images may also be associated with a look-up table. For example, in a "green" look-up table on an 8-bit image, values of zero are black, values of 128 are darkish green, and the maximum value of 255 is fully pure green. |
from ij import IJ, ImagePlus
# Grab the last activated image
imp = IJ.getImage()
# Print image details
print "title:", imp.title
print "width:", imp.width
print "height:", imp.height
print "number of pixels:", imp.width * imp.height
print "number of slices:", imp.getNSlices()
print "number of channels:", imp.getNChannels()
print "number of time frames:", imp.getNFrames()
types = {ImagePlus.COLOR_RGB : "RGB",
ImagePlus.GRAY8 : "8-bit",
ImagePlus.GRAY16 : "16-bit",
ImagePlus.GRAY32 : "32-bit",
ImagePlus.COLOR_256 : "8-bit color"}
print "image type:", types[imp.type]
Started New_.py at Wed Nov 10 14:57:46 CET 2010
title: boats.gif
width: 720
height: 576
number of pixels: 414720
number of slices: 1
number of channels: 1
number of time frames: 1
image type: 8-bit
|
|
|
Obtaining pixel statistics of an image (and your first function) ImageJ / Fiji offers an ImageStatistics class that does all the work for us. Notice how we import the ImageStatistics namespace as "IS", i.e. we alias it--it's too long to type! The options variable is the bitwise-or combination of three different static fields of the ImageStatistics class. The final options is an integer that has specific bits set that indicate mean, median and min and max values. |
from ij import IJ
from ij.process import ImageStatistics as IS
# Grab the active image
imp = IJ.getImage()
# Get its ImageProcessor
ip = imp.getProcessor()
options = IS.MEAN | IS.MEDIAN | IS.MIN_MAX
stats = IS.getStatistics(ip, options, imp.getCalibration())
# print statistics on the image
print "Image statistics for", imp.title
print "Mean:", stats.mean
print "Median:", stats.median
print "Min and max:", stats.min, "-", stats.max
Started New_.py at Wed Nov 10 19:54:37 CET 2010
Image statistics for boats.gif
Mean: 120.026837384
Median: 138.0
Min and max: 3.0 - 220.0
|
|
|
Now, how about obtaining statistics for a lot of images?
So we define a folder that contains our images, and we loop the list of filenames that it has. For every filename that ends with ".tif", we load it as an ImagePlus, and handle it to the getStatistics function, which returns us the mean, median, and min and max values. (Note: if the images are stacks, use StackStatistics instead.) This script introduces a few new concepts:
See also the python documentation page on control flow, with explanations on the keywords if, else and elif, the for loop keyword and the break and continue keywords, defining a function with def, functions with variable number of arguments, anonymous functions (with the keyword lambda), and guidelines on coding style. |
from ij import IJ
from ij.process import ImageStatistics as IS
import os
options = IS.MEAN | IS.MEDIAN | IS.MIN_MAX
def getStatistics(imp):
""" Return statistics for the given ImagePlus """
global options
ip = imp.getProcessor()
stats = IS.getStatistics(ip, options, imp.getCalibration())
return stats.mean, stats.median, stats.min, stats.max
# Folder to read all images from:
folder = "/home/albert/Desktop/t2/fiji-tutorial"
# Get statistics for each image in the folder
# whose file extension is ".tif":
for filename in os.listdir(folder):
if filename.endswith(".tif"):
print "Processing", filename
imp = IJ.openImage(os.path.join(folder, filename))
if imp is None:
print "Could not open image from file:", filename
continue
mean, median, min, max = getStatistics(imp)
print "Image statistics for", imp.title
print "Mean:", mean
print "Median:", median
print "Min and max:", min, "-", max
else:
print "Ignoring", filename
|
|
|
Iterating pixels Iterating pixels is considered a low-level operation that you would seldom, if ever, have to do. But just so you can do it when you need to, here are various ways to iterate all pixels in an image. The three iteration methods:
The last should be your preferred method. There's the least opportunity for introducting an error, and it is very concise. Regarding the example given, keep in mind:
|
from ij import IJ
# Grab the active image
imp = IJ.getImage()
# Grab the image processor converted to float values
# to avoid problems with bytes
ip = imp.getProcessor().convertToFloat() # as a copy
# The pixels points to an array of floats
pixels = ip.getPixels()
print "Image is", imp.title, "of type", imp.type
# Obtain the minimum pixel value
# Method 1: the for loop, C style
minimum = Float.MAX_VALUE
for i in xrange(len(pixels)):
if pixels[i] < minimum:
minimum = pixels[i]
print "1. Minimum is:", minimum
# Method 2: iterate pixels as a list
minimum = Float.MAX_VALUE
for pix in pixels:
if pix < minimum:
minimum = pix
print "2. Minimum is:", minimum
# Method 3: apply the built-in min function
# to the first pair of pixels,
# and then to the result of that and the next pixel, etc.
minimum = reduce(min, pixels)
print "3. Minimum is:", minimum
Started New_.py at Wed Nov 10 20:49:31 CET 2010
Image is boats.gif of type 0
1. Minimum is: 3.0
2. Minimum is: 3.0
3. Minimum is: 3.0
|
|
|
On iterating or looping lists or collections of elements Ultimately all operations that involve iterating a list or a collection of elements can be done with the for looping construct. But in almost all occasions the for is not the best choice, neither regarding performance nor in clarity or conciseness. The latter is important to minimize the amount of errors that we may possibly introduce without noticing. There are three kinds of operations to perform on lists or collections: map, reduce, and filter. We show them here along with the equivalent for loop. |
||
|
A map operation takes a list of length N and returns another list also of length N, with the results of applying a function (that takes a single argument) to every element of the original list. For example, suppose you want to get a list of all images open in Fiji. With the for loop, we have to create first a list explictly and then append one by one every image. With list comprehension, the list is created directly and the logic of what goes in it is placed inside the square brackets--but it is still a loop. That is, it is still a sequential, unparallelizable operation. With the map, we obtain the list automatically by executing the function WM.getImage to every ID in the list of IDs. While this is a trivial example, suppose you were executing a complex operation on every element of a list or an array. If you were to redefine the map function to work in parallel, suddenly any map operation in your program will run faster, without you having to modify a single line of tested code! |
from ij import WindowManager as WM
# Method 1: with a 'for' loop
images = []
for id in WM.getIDList():
images.append(WM.getImage(id))
# Method 2: with list comprehension
images = [WM.getImage(id) for id in WM.getIDList()]
# Method 3: with a 'map' operation
images = map(WM.getImage, WM.getIDList())
|
|
|
A filter operation takes a list of length N and returns a shorter list, with anywhere from 0 to N elements. Only those elements of the original list that pass a test are placed in the new, returned list. For example, suppose you want to find the subset of opened images in Fiji whose title match a specific criterium. With the for loop, we have to create a new list first, and then append elements to that list as we iterate the list of images. The second variant of the for loop uses list comprehension. The code is reduced to a single short line, which is readable, but is still a python loop (with potentially lower performance). With the filter operation, we get the (potentially) shorter list automatically. The code is a single short line, instead of 4 lines! |
from ij import WindowManager as WM
# A list of all open images
imps = map(WM.getImage, WM.getIDList())
def match(imp):
""" Returns true if the image title contains the word 'cochlea'"""
return imp.title.find("cochlea") > -1
# Method 1: with a 'for' loop
# (We have to explicitly create a new list)
matching = []
for imp in imps:
if match(imp):
matching.append(imp)
# Method 2: with list comprehension
matching = [imp for imp in imps if match(imp)]
# Method 3: with a 'filter' operation
matching = filter(match, imps)
|
|
|
A reduce operation takes a list of length N and returns a single value. This value is composed from applying a function that takes two arguments to the first two elements of the list, then to the result of that and the next element, etc. Optionally an initial value may be provided, so that the cycle starts with that value and the first element of the list. For example, suppose you want to find the largest image, by area, from the list of all opened images in Fiji. With the for loop, we have to we have to keep track of which was the largest area in a pair of temporary variables. And even check whether the stored largest image is null! We could have initizalized the largestArea variable to the first element of the list, and then start looping at the second element by slicing the first element off the list (with "for imp in imps[1:]:"), but then we would have had to check if the list contains at least one element. With the reduce operation, we don't need any temporary variables. All we need is to define a helper function (which could have been an anonymous lambda function, but we defined it explicitly for extra clarity and reusability). |
from ij import IJ
from ij import WindowManager as WM
# A list of all open images
imps = map(WM.getImage, WM.getIDList())
def area(imp):
return imp.width * imp.height
# Method 1: with a 'for' loop
largest = None
largestArea = 0
for imp in imps:
if largest is None:
largest = imp
else:
a = area(imp)
if a > largestArea:
largest = imp
largestArea = a
# Method 2: with a 'reduce' operation
def largestImage(imp1, imp2):
return imp1 if area(imp1) > area(imp2) else imp2
largest = reduce(largestImage, imps)
|
|
|
Subtract the min value to every pixel First we obtain the minimum pixel value, using the reduce method explained just above. Then we subtract this minimum value to every pixel. We have two ways to do it:
With the first method, since the pixels array was already a copy (notice we called convertToFloat() on the ImageProcessor), we can use it to create a new ImagePlus with it without any unintended consequences. With the second method, the new list of pixels must be given to a new FloatProcessor instance, and with it, a new ImagePlus is created, of the same dimensions as the original. |
from ij import IJ
imp = IJ.getImage()
ip = imp.getProcessor().convertToFloat() # as a copy
pixels = ip.getPixels()
# Apply the built-in min function
# to the first pair of pixels,
# and then to the result of that and the next pixel, etc.
minimum = reduce(min, pixels)
# Method 1: subtract the minim from every pixel,
# in place, modifying the pixels array
for i in xrange(len(pixels)):
pixels[i] -= minimum
# ... and create a new image:
imp2 = ImagePlus(imp.title, ip)
# Method 2: subtract the minimum from every pixel
# and store the result in a new array
pixels3 = map(lambda x: x - minimum, pixels)
# ... and create a new image:
ip3 = FloatProcessor(ip.width, ip.height, pixels3, None)
imp3 = ImagePlus(imp.title, ip3)
# Show the images in an ImageWindow:
imp2.show()
imp3.show()
|
|
|
Reduce a pixel array to a single value: count pixels above a threshold Suppose you want to analyze a subset of pixels. For example, find out how many pixels have a value over a certain threshold. The reduce built-in function is made just for that. It takes a function with two arguments (the running count and the next pixel); the list or array of pixels; and an initial value (in this case, zero) for the first argument (the "count'), and will return a single value (the total count). In this example, we computed first the mean pixel intensity, and then filtered all pixels for those whose value is above the mean. Notice that we compute the mean by using the convenient built-in function sum, which is able to add all numbers contained in any kind of collection (be it a list, a native array, a set of unique elements, or the keys of a dictionary). We could imitate the built-in sum function with reduce(lambda s, x: s + x, pixels), but paying a price in performance. Notice we are using anonymous functions again (that is, functions that lack a name), declared in place with the lambda keyword. A function defined with def would do just fine as well. |
from ij import IJ
# Grab currently active image
imp = IJ.getImage()
ip = imp.getProcessor().convertToFloat()
pixels = ip.getPixels()
# Compute the mean value (sum of all divided by number of pixels)
mean = sum(pixels) / len(pixels)
# Count the number of pixels above the mean
n_pix_above = reduce(lambda count, a: count + 1 if a > mean else count, pixels, 0)
print "Mean value", mean
print "% pixels above mean:", n_pix_above / float(len(pixels)) * 100
Started New_.py at Thu Nov 11 01:50:49 CET 2010
Mean value 120.233899981
% pixels above mean: 66.4093846451
|
|
|
Another useful application of filtering pixels by their value: finding the coordinates of all pixels above a certain value (in this case, the mean), and then calculating their center of mass. The filter built-in function is made just for that. The indices of the pixels whose value is above the mean are collected in a list named "above", which is created by filtering the indices of all pixels (provided by the built-in function xrange). The filtering is done by an anonymous function declared with lambda, with a single argument: the index i of the pixel. Here, note that in ImageJ, the pixels of an image are stored in a linear array. The length of the array is width * height, and the pixels are stored as concatenated rows. Therefore, the modulus of dividing the index of a pixel by the width the image provides the X coordinate of a pixel. Similarly, the integer division of the index of a pixel by the width provides the Y coordinate. To compute the center of mass, there are two equivalent methods. The C-style method with a for loop, with every variable being declared prior to the loop, and then modified at each loop iteration and, after the loop, dividing the sum of coordinates by the number of coordinates (the length of the "above" list). For this example, this is the method with the best performance.The second method computes the X and Y coordinates of the center of mass with a single line of code for each. Notice that both lines are nearly identical, differing only in the body of the function mapped to the "above" list containing the indices of the pixels whose value is above the mean. While, in this case, the method is less performant due to repeated iteration of the list "above", the code is shorter, easier to read, and with far less opportunities for introducing errors. If the actual computation was far more expensive than the simple calculation of the coordinates of a pixel given its index in the array of pixels, this method would pay off for its clarity. |
from ij import IJ
# Grab the currently active image
imp = IJ.getImage()
ip = imp.getProcessor().convertToFloat()
pixels = ip.getPixels()
# Compute the mean value
mean = sum(pixels) / len(pixels)
# Obtain the list of indices of pixels whose value is above the mean
above = filter(lambda i: pixels[i] > mean, xrange(len(pixels)))
print "Number of pixels above mean value:", len(above)
# Measure the center of mass of all pixels above the mean
# The width of the image, necessary for computing the x,y coordinate of each pixel
width = imp.width
# Method 1: with a for loop
xc = 0
yc = 0
for i in above:
xc += i % width # the X coordinate of pixel at index i
yc += i / width # the Y coordinate of pixel at index i
xc = xc / len(above)
yc = yc / len(above)
print xc, yc
# Method 2: with sum and map
xc = sum(map(lambda i: i % width, above)) / len(above)
yc = sum(map(lambda i: i / width, above)) / len(above)
print xc, yc
|
|
|
The third method pushes the functional approach too far. While written in a single line, that doesn't mean it is clearer to read: it's intent is obfuscated by starting from the end: the list comprehension starts off by stating that each element (there are only two) of the list resulting from the reduce has to be divided by the length of the list of pixels "above", and only then we learn than the collection being iterated is the array of two coordinates, created at every iteration over the list "above", containing the sum of all coordinates for X and for Y. Notice that the reduce is invoked with three arguments, the third one being the list [0, 0] containing the initialization values of the sums. Confusing! Avoid writing code like this. Notice as well that, by creating a new list at every iteration step, this method is the least performant of all. The fourth method is a clean up of the third method. Notice that we import the partial function from the functools package. With it, we are able to create a version of the "accum" helper function that has a frozen "width" argument (also known as currying a function). In this way, the "accum" function is seen by the reduce as a two-argument function (which is what reduce needs here). While we regain the performance of the for loop, notice that now the code is just as long as with the for loop. The purpose of writing this example is to illustrate how one can write python code that doesn't use temporary variables, these generally being potential points of error in a computer program. It is always better to write lots of small functions that are easy to read, easy to test, free of side effects, documented, and that then can be used to assemble our program. |
# (Continues from above...)
# Method 3: iterating the list "above" just once
xc, yc = [d / len(above) for d in
reduce(lambda c, i: [c[0] + i % width, c[1] + i / width], above, [0, 0])]
print xc, yc
# Method 4: iterating the list "above" just once, more clearly and performant
from functools import partial
def accum(width, c, i):
""" Accumulate the sum of the X,Y coordinates of index i in the list c."""
c[0] += i % width
c[1] += i / width
return c
xy, yc = [d / len(above) for d in reduce(partial(accum, width), above, [0, 0])]
print xc, yc
|
|
4. Running ImageJ / Fiji plugins on an ImagePlusHere is an example plugin run programmatically: a median filter applied to the currently active image. The median filter, along with the mean, minimum, maximum, variance, remove outliers and despeckle menu commands, are implemented in the RankFilters class. |
from ij import IJ
from ij.plugin.filter import RankFilters
# Grab the active image
imp = IJ.getImage()
ip = imp.getProcessor().convertToFloat() # as a copy
# Remove noise by running a median filter
# with a radius of 2
radius = 2
RankFilters().rank(ip, radius, RankFilters.MEDIAN)
imp2 = ImagePlus(imp.title + " median filtered", ip)
imp2.show()
|
|
|
Finding the class that implements a specific ImageJ command When starting ImageJ/Fiji programming, the problem is not so much how to run a plugin on an image, as it is to find out which class implements which plugin. Here is a simple method to find out, via the Command Finder:
Notice that the plugin class comes with some text. For example:
FFT (in Process > FFT) [ij.plugin.FFT("fft")]
Inverse FFT (in Process > FFT) [ij.plugin.FFT("inverse")]
The above two commands are implemented by a single plugin (ij.plugin.FFT) whose run method accepts, like all PlugIn, a text string specifying the action: the fft, or the inverse. |
|
|
|
Finding the java documentation for any class Once you have found the PlugIn class that implements a specific command, you may want to use that class directly. The information is either in the online java documentation or in the source code. How to find these?
|
|
|
|
Figuring out what parameters a plugin requires To do that, we'll use the Macro Recorder. Make sure that an image is open. Then:
That is valid macro code, that ImageJ can execute. The first part is the command ("Median..."), the second part is the parameters that that command uses; in this case, just one ("radius=2"). If there were more parameters, they would be separated by spaces. |
|
|
|
Running a command on an image We can use these macro recordings to create jython code that executes a given plugin on a given image. Here is an example. Very simple! The IJ namespace has a function, run, that accepts an ImagePlus as first argument, then the name of the command to run, and then the macro-ready list of arguments that the command requires. |
from ij import IJ
# Grab the active image
imp = IJ.getImage()
# Run the median filter on it, with a radius of 2
IJ.run(imp, "Median...", "radius=2")
|
|
5. Creating images and regions of interest (ROIs)Create an image from scratch An ImageJ/Fiji image is composed of at least three objects:
In the example, we create an empty array of floats (see creating native arrays), and fill it in with random float values. Then we give it to a FloatProcessor instance, which is then wrapped by an ImagePlus instance. Voilà! |
from ij import ImagePlus
from ij.process import FloatProcessor
from array import zeros
from random import random
width = 1024
height = 1024
pixels = zeros('f', width * height)
for i in xrange(len(pixels)):
pixels[i] = random()
fp = FloatProcessor(width, height, pixels, None)
imp = ImagePlus("White noise", fp)
imp.show()
|
|
|
Fill a region of interest (ROI) with a given value To fill a region of interest in an image, we could iterate the pixels, find the pixels that lay within the bounds of interest, and set their values to a specified value. But that tedious and error prone. Much more effective is to create an instance of a Roi class or one of its subclasses (PolygonRoi, OvalRoi, ShapeRoi, etc.) and tell the ImageProcessor to fill that region. In this example, we create an image filled with white noise like before, and then define a rectangular region of interest in it, which is filled with a value of 2.0. The white noise is drawn from a random distribution whose values range from 0 to 1. When filling an area of the FloatProcessor with a value of 2.0, that is the new maximum value. The area with 2.0 pixel values will look white (look at the status bar): 
|
from ij import IJ, ImagePlus
from ij.process import FloatProcessor
from array import zeros
from random import random
from ij.gui import Roi, PolygonRoi
# Create a new ImagePlus filled with noise
width = 1024
height = 1024
pixels = zeros('f', width * height)
for i in xrange(len(pixels)):
pixels[i] = random()
fp = FloatProcessor(width, height, pixels, None)
imp = ImagePlus("Random", fp)
# Fill a rectangular region of interest
# with a value of 2:
roi = Roi(400, 200, 400, 300)
fp.setRoi(roi)
fp.setValue(2.0)
fp.fill()
# Fill a polygonal region of interest
# with a value of -3
xs = [234, 174, 162, 102, 120, 123, 153, 177, 171,
60, 0, 18, 63, 132, 84, 129, 69, 174, 150,
183, 207, 198, 303, 231, 258, 234, 276, 327,
378, 312, 228, 225, 246, 282, 261, 252]
ys = [48, 0, 60, 18, 78, 156, 201, 213, 270, 279,
336, 405, 345, 348, 483, 615, 654, 639, 495,
444, 480, 648, 651, 609, 456, 327, 330, 432,
408, 273, 273, 204, 189, 126, 57, 6]
proi = PolygonRoi(xs, ys, len(xs), Roi.POLYGON)
fp.setRoi(proi)
fp.setValue(-3)
fp.fill(proi.getMask()) # Attention!
imp.show()
|
|
6. Create and manipulate image stacks and hyperstacksLoad a color image stack and extract its green channel First we load the stack from the web--it's the "Fly Brain" sample image. Then we iterate its slices. Each slice is a ColorProcessor: wraps an integer array. Each integer is represented by 4 bytes, and the lower 3 bytes represent, respectively, the intensity values for red, green and blue. The upper most byte is usually reserved for alpha (the inverse of transparency), but ImageJ ignores it. Dealing with low-level details like that is not necessary. The ColorProcessor has a method, toFloat, that can give us a FloatProcessor for a specific color channel. Red is 0, green is 1, and blue is 2. Representing the color channel in floats is most convenient for further processing of the pixel values--won't overflow like a byte would. In this example, all we do is collect each slice into a list of slices we named greens. Then we add all the slices to a new ImageStack, and pass it to a new ImagePlus. Then we invoke the "Green" command on that ImagePlus instance, so that a linear green look-up table is assigned to it. And we show it. |
from ij import IJ, ImagePlus, ImageStack
# Load a stack of images: a fly brain, in RGB
imp = IJ.openImage("http://imagej.nih.gov/ij/images/flybrain.zip")
stack = imp.getImageStack()
print "number of slices:", imp.getNSlices()
# A list of green slices
greens = []
# Iterate each slice in the stack
for i in xrange(1, imp.getNSlices()+1):
# Get the ColorProcessor slice at index i
cp = stack.getProcessor(i)
# Get its green channel as a FloatProcessor
fp = cp.toFloat(1, None)
# ... and store it in a list
greens.append(fp)
# Create a new stack with only the green channel
stack2 = ImageStack(imp.width, imp.height)
for fp in greens:
stack2.addSlice(None, fp)
# Create a new image with the stack of green channel slices
imp2 = ImagePlus("Green channel", stack2)
# Set a green look-up table:
IJ.run(imp2, "Green", "")
imp2.show()
|
|
|
Convert an RGB stack to a 2-channel, 32-bit hyperstack We load an RGB stack--the "Fly brain" sample image, as before. Suppose we want to analyze each color channel independently: an RGB image doesn't really let us, without lots of low-level work to disentangle the different color values from each pixel value. So we convert the RGB stack to a hyperstack with two separate channels, where each channel slice is a 32-bit FloatProcessor. The first step is to create a new ImageStack instance, to hold all the slices that we'll need: one per color channel, times the number of slices. Realize that we could have 7 channels if we wanted, or 20, for each slice. As many as you want. The final step is to open the hyperstack. For that:
Open the "Image - Color - Channels Tool" and you'll see that the Composite image is set to show only the red channel--try checking the second channel as well. For a real-world example of a python script that uses hyperstacks, see the Correct_3D_drift.py script (available as the command "Plugins - Registration - Correct 3D drift"). |
from ij import IJ, ImagePlus, ImageStack, CompositeImage
# Load a stack of images: a fly brain, in RGB
imp = IJ.openImage("http://imagej.nih.gov/ij/images/flybrain.zip")
stack = imp.getImageStack()
# A new stack to hold the data of the hyperstack
stack2 = ImageStack(imp.width, imp.height)
# Convert each color slice in the stack
# to two 32-bit FloatProcessor slices
for i in xrange(1, imp.getNSlices()+1):
# Get the ColorProcessor slice at index i
cp = stack.getProcessor(i)
# Extract the red and green channels as FloatProcessor
red = cp.toFloat(0, None)
green = cp.toFloat(1, None)
# Add both to the new stack
stack2.addSlice(None, red)
stack2.addSlice(None, green)
# Create a new ImagePlus with the new stack
imp2 = ImagePlus("32-bit 2-channel composite", stack2)
imp2.setCalibration(imp.getCalibration().copy())
# Tell the ImagePlus to represent the slices in its stack
# in hyperstack form, and open it as a CompositeImage:
nChannels = 2 # two color channels
nSlices = stack.getSize() # the number of slices of the original stack
nFrames = 1 # only one time point
imp2.setDimensions(nChannels, nSlices, nFrames)
comp = CompositeImage(imp2, CompositeImage.COLOR)
comp.show()
|
|
7. Interacting with humans: file and option dialogs, messages, progress bars.Ask the user for a directory See DirectoryChooser. |
from ij.io import DirectoryChooser
dc = DirectoryChooser("Choose a folder")
folder = dc.getDirectory()
if folder is None:
print "User canceled the dialog!"
else
print "Selected folder:", folder
|
|
|
Ask the user for a file See OpenDialog and SaveDialog. |
from ij.io import OpenDialog
od = OpenDialog("Choose a file", None)
filename = od.getFileName()
if filename is None:
print "User canceled the dialog!"
else:
directory = od.getDirectory()
filepath = directory + filename
print "Selected file path:", filepath
|
|
|
Ask the user to enter a few parameters in a dialog There are more possibilities, but these are the basics. See GenericDialog. All plugins that use a GenericDialog are automatable. Remember how above we run a command on an image? When the names in the dialog fields match the names in the macro string, the dialog is fed in the values automatically. If a dialog field doesn't have a match, it takes the default value as defined in the dialog declaration. If a plugin was using a dialog like the one we built here, we would run it automatically like this:
args = "name='first' alpha=0.5 output='32-bit' scale=80"
IJ.run(imp, "Some PlugIn", args)
Above, leaving out the word 'optimize' means that it will use the default value (True) for it. |
from ij.gui import GenericDialog
def getOptions():
gd = GenericDialog("Options")
gd.addStringField("name", "Untitled")
gd.addNumericField("alpha", 0.25, 2) # show 2 decimals
gd.addCheckbox("optimize", True)
types = ["8-bit", "16-bit", "32-bit"]
gd.addChoice("output as", types, types[2])
gd.addSlider("scale", 1, 100, 100)
gd.showDialog()
#
if gd.wasCanceled():
print "User canceled dialog!"
return
# Read out the options
name = gd.getNextString()
alpha = gd.getNextNumber()
optimize = gd.getNextBoolean()
output = gd.getNextChoice()
scale = gd.getNextNumber()
return name, alpha, optimize, output, scale
options = getOptions()
if options is not None:
name, alpha, optimize, output, scale = options
print name, alpha, optimize, output, scale
|
|
|
Show a progress bar Will show a progress bar in the Fiji window. |
from ij import IJ
imp = IJ.getImage()
stack = imp.getImageStack()
for i in xrange(1, stack.getSize()+1):
# Report progress
IJ.showProgress(i, stack.getSize()+1)
# Do some processing
ip = stack.getProcessor(i)
# ...
# Signal completion
IJ.showProgress(1)
|
|
8. Turn your script into a pluginSave the script in Fiji's plugins folder or a subfolder, with:
Then run "Help - Update Menus", or restart Fiji. That's it! The script will appear as a regular menu command under "Plugins", and you'll be able to run it from the Command Finder. Where is the plugins folder?
See also the Fiji wiki on Jython Scripting. |
||
9. Lists, native arrays, and interacting with Java classesJython lists are passed as read-only arrays to Java classes Calling java classes and methods for jython is seamless: on the surface, there isn't any difference with calling jython classes and methods. But there is a subtle difference when calling java methods that expect native arrays. Jython will automatically present a jython list as a native array to the java method that expects it. But as read-only! In this example, we create an AffineTransform that specifies a translation. Then we give it:
The ability to pass jython lists as native arrays to java methods is extremely convenient, and we have used it in the example above to pass a list of strings to the GenericDialog addChoice method. |
from java.awt.geom import AffineTransform
from array import array
# A 2D point
x = 10
y = 40
# A transform that does a translation
# of dx=45, dy=56
aff = AffineTransform(1, 0, 0, 1, 45, 56)
# Create a point as a list of x,y
p = [x, y]
aff.transform(p, 0, p, 0, 1)
print p
# prints: [10, 40] -- the list p was not updated!
# Create a point as a native float array of x,y
q = array('f', [x, y])
aff.transform(q, 0, q, 0, 1)
print q
# prints: [55, 96] -- the native array q was updated properly
Started New_.py at Sat Nov 13 09:31:51 CET 2010
[10, 40]
array('f', [55.0, 96.0])
|
|
|
Creating native arrays: empty, or from a list The package array contains two functions:
The type of array is specified by the first argument. For primitive types (char, short, int, float, long, double), use a single character in quotes. See the list of all possible characters. Manipulating arrays is done in the same way that you would do in java. See lines 16--18 in the example. See also the documentation on how to create multidimensional native arrays with Jython. |
from array import array, zeros
from ij import ImagePlus
# An empty native float array of length 5
a = zeros('f', 5)
print a
# A native float array with values 0 to 9
b = array('f', [0, 1, 2, 3, 4])
print b
# An empty native ImagePlus array of length 5
imps = zeros(ImagePlus, 5)
print imps
# Assign the current image to the first element of the array
imps[0] = IJ.getImage()
print imps
# Length of an array
print "length:", len(imps)
Started New_.py at Sat Nov 13 09:40:00 CET 2010
array('f', [0.0, 0.0, 0.0, 0.0, 0.0])
array('f', [0.0, 1.0, 2.0, 3.0, 4.0])
array(ij.ImagePlus, [None, None, None, None, None])
array(ij.ImagePlus, [imp[boats.gif 720x576x1], None, None, None, None])
length: 5
|
|
10. Generic algorithms that work on images of any kind: using ImgLibImglib is a general-purpose software library for n-dimensional data processing, mostly oriented towards images. Scripting with Imglib greatly simplifies operations on images of different types (8-bit, 16-bit, color images, etc). Scripting in imglib is based around the Compute function, which composes images, functions and numbers into output images. |
||
|
Mathematical operations on images The script.imglib packages offers means to compute with images. There are three kinds of operations, each in its own package:
|
from script.imglib.math import Compute, Subtract
from script.imglib.color import Red, Green, Blue, RGBA
from script.imglib import ImgLib
# Open an RGB image stack
imp = IJ.openImage("http://imagej.nih.gov/ij/images/flybrain.zip")
# Wrap it as an Imglib image
img = ImgLib.wrap(imp)
# Example 1: subtract red from green channel
sub = Compute.inFloats(Subtract(Green(img), Red(img)))
ImgLib.wrap(sub).show()
# Example 2: subtract red from green channel, and compose a new RGBA image
rgb = RGBA(Red(img), Subtract(Green(img), Red(img)), Blue(img)).asImage()
ImgLib.wrap(rgb).show()
|
|
|
Using image math for flat-field correction In the example, we start by opening an image from the sample image collection of ImageJ.
With imglib, all the above operations happen in a pixel-by-pixel basis, and are computed as fast or faster than if you had manually hand-coded every operation. And multithreaded! |
from script.imglib.math import Compute, Divide, Multiply, Subtract
from script.imglib.algorithm import Gauss, Scale2D, Resample
from mpicbg.imglib import ImgLib
# 1. Open an image
img = ImgLib.wrap(IJ.openImage("http://imagej.nih.gov/ij/images/bridge.gif"))
# 2. Simulate a brighfield from a Gauss with a large radius
# (First scale down by 4x, then gauss of radius=20, then scale up)
brightfield = Resample(Gauss(Scale2D(img, 0.25), 20), img.getDimensions())
# 3. Simulate a perfect darkfield
darkfield = 0
# 4. Compute the mean pixel intensity value of the image
mean = reduce(lambda s, t: s + t.get(), img, 0) / img.size()
# 5. Correct the illumination
corrected = Compute.inFloats(Multiply(Divide(Subtract(img, brightfield),
Subtract(brightfield, darkfield)), mean))
# 6. ... and show it in ImageJ
ImgLib.wrap(corrected).show()
|
|
|
Extracting and manipulating image color channels: RGBA and HSB In the examples above we have already used the Red and Green functions. There's also Blue, Alpha, and a generic Channel that takes the channel index as argument--where red is 3, green is 2, blue is 1, and alpha is 4 (these numbers are related to the byte order in the 4-byte that makes up a 32-bit integer). The basic color operations have to do with extracting the color channel, for a particular color space (RGBA or HSB) The function RGBA takes from 1 to 4 arguments, and creates an RGBA image out of them. These arguments can be images, other functions, or numbers--for example, all pixels of a channel would have the value 255 (maximum intensity). |
from script.imglib.math import Compute, Subtract, Multiply
from script.imglib.color import Red, Blue, RGBA
from script.imglib.algorithm import Gauss, Dither
# Obtain a color image from the ImageJ samples
clown = ImgLib.wrap(IJ.openImage("http://imagej.nih.gov/ij/images/clown.jpg"))
# Example 1: compose a new image manipulating the color channels of the clown image:
img = RGBA(Gauss(Red(clown), 10), 40, Multiply(255, Dither(Blue(clown)))).asImage()
ImgLib.wrap(img).show()
| |
|
In the second example, we extract the HSB channels from the clown image. To the Hue channel (which is expressed in the range [0, 1]), we add 0.5. We've shifted the hue around a bit.
|
from script.imglib.math import Compute, Add, Subtract
from script.imglib.color import HSB, Hue, Saturation, Brightness
from script.imglib import ImgLib
# Obtain an image
img = ImgLib.wrap(IJ.openImage("http://imagej.nih.gov/ij/images/clown.jpg"))
# Obtain a new clown, whose hue has been shifted by half
# with the same saturation and brightness of the original
bluey = Compute.inRGBA(HSB(Add(Hue(img), 0.5), Saturation(img), Brightness(img)))
ImgLib.wrap(bluey).show()
|
|
|
In the third example, we apply a gamma correction to an RGB confocal stack. To correct the gamma, we must first extract each color channel from the image, and then apply the gamma to each channel independently. In this example we use a gamma of 0.5 for every channel. Of course you could apply different gamma values to each channel, or apply it only to specific channels. Notice how we use asImage() instead of Compute.inRGBA. The result is the same; the former is syntactic sugar of the latter.
|
# Correct gamma
from script.imglib.math import Min, Max, Exp, Multiply, Divide, Log
from script.imglib.color import RGBA, Red, Green, Blue
gamma = 0.5
img = ImgLib.wrap(IJ.getImage())
def g(channel, gamma):
""" Return a function that, when evaluated, computes the gamma of the given color channel.
If 'i' was the pixel value, then this function would do:
double v = Math.exp(Math.log(i/255.0) * gamma) * 255.0);
if (v < 0) v = 0;
if (v >255) v = 255;
"""
return Min(255, Max(0, Multiply(Exp(Multiply(gamma, Log(Divide(channel, 255)))), 255)))
corrected = RGBA(g(Red(img), gamma), g(Green(img), gamma), g(Blue(img), gamma)).asImage()
ImgLib.wrap(corrected).show()
|
|
|
Find cells in an 3D image stack by Difference of Gaussian, count them, and show them in 3D as spheres. First we define the cell diameter that we are looking for (5 microns; measure it with a line ROI over the image) and the minimum voxel intensity that will care about (in this case, anything under a value of 40 will be ignored). And we load the image of interest: a 3-color channel image of the first instar Drosophila larval brain. Then we scale down the image to make it isotropic: so that voxels have the same dimensions in all axes. We run the DoGPeaks ("Difference of Gaussian Peaks") with a pair of appropriate sigmas: the scaled diameter of the cell, and half that. The peaks are each a float[] array that specifies its coordinate. With these, we create Point3f instances, which we transport back to calibrated image coordinates. Finally, we show in the 3D Viewer the peaks as spheres, and the image as a 3D volume. 
|
# Load an image of the Drosophila larval fly brain and segment
# the 5-micron diameter cells present in the red channel.
from script.imglib.analysis import DoGPeaks
from script.imglib.color import Red
from script.imglib.algorithm import Scale2D
from script.imglib.math import Compute
from script.imglib import ImgLib
from ij3d import Image3DUniverse
from javax.vecmath import Color3f, Point3f
cell_diameter = 5 # in microns
minPeak = 40 # The minimum intensity for a peak to be considered so.
imp = IJ.openImage("http://pacific.mpi-cbg.de/samples/first-instar-brain.zip")
# Scale the X,Y axis down to isotropy with the Z axis
cal = imp.getCalibration()
scale2D = cal.pixelWidth / cal.pixelDepth
iso = Compute.inFloats(Scale2D(Red(ImgLib.wrap(imp)), scale2D))
# Find peaks by difference of Gaussian
sigma = (cell_diameter / cal.pixelWidth) * scale2D
peaks = DoGPeaks(iso, sigma, sigma * 0.5, minPeak, 1)
print "Found", len(peaks), "peaks"
# Convert the peaks into points in calibrated image space
ps = []
for peak in peaks:
p = Point3f(peak)
p.scale(cal.pixelWidth * 1/scale2D)
ps.append(p)
# Show the peaks as spheres in 3D, along with orthoslices:
univ = Image3DUniverse(512, 512)
univ.addIcospheres(ps, Color3f(1, 0, 0), 2, cell_diameter/2, "Cells").setLocked(True)
univ.addOrthoslice(imp).setLocked(True)
univ.show()
|
