Puncta Segmentation

This notebook demonstrates how to segment puncta from an image using the “Analyze Particles” ImageJ plugin, perform measurements and return those values as a pandas DataFrame.

First let’s initialize ImageJ in headless mode with the legacy layer enabled. See How to initialize PyImageJ for more information on PyImageJ’s initialization modes.

import imagej
import scyjava as sj

# initialize imagej
ij = imagej.init(mode='headless', add_legacy=True)
print(f"ImageJ version: {ij.getVersion()}")
ImageJ version: 2.9.0/1.53t

Now that ImageJ has been successfully initialized, let’s next import the Java classes (i.e. ImageJ resources) we need for the rest of the workflow.

# get additional resources
HyperSphereShape = sj.jimport('net.imglib2.algorithm.neighborhood.HyperSphereShape')
Overlay = sj.jimport('ij.gui.Overlay')
Table = sj.jimport('org.scijava.table.Table')
ParticleAnalyzer = sj.jimport('ij.plugin.filter.ParticleAnalyzer')

Next we load the data with ij.io().open() and then convert the 16-bit image to 32-bit with ImageJ ops (i.e. ij.op().convert().int32()). Alternatively you can load the data using other software packages such as scikit-image’s skimage.io.imread() which returns a NumPy array. To convert the NumPy array into an ImageJ Java image (e.g. a Dataset: net.imagej.Dataset) call ij.py.to_dataset().

The sample images used in this notebook are available on the PyImageJ GitHub repository here.

# load test data
ds_src = ij.io().open('sample-data/test_still.tif')
ds_src = ij.op().convert().int32(ds_src) # convert image to 32-bit
ij.py.show(ds_src, cmap='binary')
_images/3284dc0ccd5b9c1cdd3de715a38b25d83e027108daed8071afa694ee287f18e8.png

There are many ways to background subtract/supress noise in an image. In this example we take a duplicate of the input image, here ds_mean, apply a mean filter with a radius of 5 and then multiply it back to the original image. The trick to this operation improving the signal-to-noise ratio lies in the application of the mean filter. In the background area outside of the cell high value pixels are generally intermixed with low value ones. In contrast the high value pixels in the object of interest are surrounded by high value ones. When the mean filter is applied, the high value pixels in the background are averaged with low values and high values pixels in the object are averaged with other high value pixels. When this mean filtered image is multiplied with the original image, the background noise is multiplied by a “smaller” number than the objects/signal of interest. Note that these are destructive changes and the resulting background supressed image is not sutible for measurements. Instead, these processes are done to generate masks and labels which can then be applied over the original data to obtain measurements that are unmodified.

Here we are using ImageJ ops to apply the mean filter. You can also apply the filter using ij.py.run_plugin(imp, "Mean...", "radius=2"). Note that using the original ImageJ plugins/functions like “Mean…” require the ImagePlus image data type. In this case we can convert the ds_mean into an ImagePlus like so:

imp = ij.py.to_imageplus(ds_mean)
# supress background noise
mean_radius = HyperSphereShape(5)
ds_mean = ij.dataset().create(ds_src.copy())
ij.op().filter().mean(ds_mean, ds_src.copy(), mean_radius)
ds_mul = ds_src * ds_mean
ij.py.show(ds_mul, cmap='binary')
_images/49a88913fc45677d7d83f6f2f3305d2bbb70a1e2d3cd29498c0e49aab45e3732.png

Here we employ a gaussian blur subtraction to futher enhance the puncta for segmentation. For more information on gaussian blurs and how they can be used in image processing please checkout Peter Bankhead’s amazing online book.

# use gaussian subtraction to enhance puncta
img_blur = ij.op().filter().gauss(ds_mul.copy(), 1.2)
img_enhanced = ds_mul - img_blur
ij.py.show(img_enhanced, cmap='binary')
_images/9583d30fa8b51f785368f8f389e5dba8809d602cfbae693e9a8e4ccdaee9dcf6.png
# apply threshold
img_thres = ij.op().threshold().renyiEntropy(img_enhanced)
ij.py.show(img_thres)
_images/dc286e55f35b139021e7789cce1a1abd17519d0d12e50d98af0db87ce39672d9.png

Finally we convert the threshold image into an ImagePlus with ij.py.to_imageplus() and run the “Analyze Particles” plugin. The ResultsTable that is returned by the “Analyze Particles” plugin first needs to be converted to a SciJava table (i.e. org.scijava.table.Table) before it can be converted to a pandas DataFrame.

# convert ImgPlus to ImagePlus
imp_thres = ij.py.to_imageplus(img_thres)

# get ResultsTable and set ParticleAnalyzer
rt = ij.ResultsTable.getResultsTable()
ParticleAnalyzer.setResultsTable(rt)

# set measurements
ij.IJ.run("Set Measurements...", "area center shape")

# run the analyze particle plugin
ij.py.run_plugin(plugin="Analyze Particles...", args="clear", imp=imp_thres)

# convert results table -> scijava table -> pandas dataframe
sci_table = ij.convert().convert(rt, Table)
df = ij.py.from_java(sci_table)

# print dataframe
print(df)
Operating in headless mode - the original ImageJ will have limited functionality.
Operating in headless mode - the ResultsTable class will not be fully functional.
Operating in headless mode - the IJ class will not be fully functional.
      Area          XM          YM     Circ.        AR     Round  Solidity
0  89891.0  149.968879  149.990049  0.779948  1.000123  0.999877  0.998789