Fiji-ops regions and channels

fiji
imglib2
analyze-particles
imagej
roi
python
imagej-ops
scripting

#1

Hi everyone,
I’m trying to die into the Ops of imageJ/Fiji and I’m struggling with a script where I’m trying to merge a few of the examples to basically find objects in a 4D (xycz) image and measure a few things on those. A short version of the code here:

#@ OpService ops
#@ DatasetService ds
#@ net.imagej.Dataset image_in

#@ OUTPUT imp_thresholded
#@ OUTPUT labelingIndex

from net.imagej.axis import Axes
from net.imagej.ops import Ops
from net.imglib2.algorithm.labeling.ConnectedComponents import StructuringElement
from net.imglib2.type.logic import BitType
from net.imglib2.roi import Regions
from net.imglib2.roi.labeling import LabelRegions

# Get dimension indexes
x_dim_index = image_in.dimensionIndex(Axes.X)
y_dim_index = image_in.dimensionIndex(Axes.Y)
z_dim_index = image_in.dimensionIndex(Axes.Z)
c_dim_index = image_in.dimensionIndex(Axes.CHANNEL)

# Get dimension lengths
x_length = image_in.dimension(x_dim_index)
y_length = image_in.dimension(y_dim_index)
z_length = image_in.dimension(z_dim_index)
c_length = image_in.dimension(c_dim_index)
    
# create the otsu op
otsu = ops.op(Ops.Threshold.Otsu, image_in)

# create memory for the thresholded image
imp_thresholded = ops.create().img(image_in, BitType())

# call threshold op slice-wise for the defined axes
ops.slice(imp_thresholded, image_in, otsu, [x_dim_index, y_dim_index, z_dim_index])

# call connected components to label each connected region
labeling = ops.labeling().cca(imp_thresholded, StructuringElement.FOUR_CONNECTED)

# get the index image (each object will have a unique gray level)
labelingIndex = labeling.getIndexImg()

# get the collection of regions and loop through them
regions = LabelRegions(labeling)

for region in regions:
    size = region.size()
    if size > 10:
	    print '----------'
	    print 'Label: ', region.getLabel()
	    print 'Size:', size
	
	    # Calculate the center of mass using region.getCenterOfMass()
	    center_of_mass = region.getCenterOfMass() 
	    print 'CoM x:', center_of_mass.getDoublePosition(x_dim_index)
	    print 'CoM y:', center_of_mass.getDoublePosition(y_dim_index)
	    print 'CoM z:', center_of_mass.getDoublePosition(z_dim_index)
	
	    # Get the mean intensity the region by "sampling" the intensity of the input image at the region pixels
	    intensity = ops.stats().mean(Regions.sample(region, image_in)).getRealFloat()
	    print "intensity:", intensity
	    
	    center_of_gravity_op = ops.op(Ops.Geometric.CenterOfGravity, Regions.sample(region, image_in))
	    
	    print 'CoG:', center_of_gravity_op.calculate()

imp_thresholded = ds.create(imp_thresholded)
labelingIndex = ds.create(labelingIndex)

I have a few questions:

  1. When I run this on my images or on the sample ‘Confocal_series’ I find that the Center of Mass calculated with region.getCenterOfMass() is slightly different than the Center of Gravity calculated through the Ops.Geometric.CenterOfGravity. I’m not sure I’m doing something wrong, sampling multiple channels,… I guess I could crop the channels and loop through them to be sure about this but, in the principle, how does this work?
  2. Is there an elegant way to extract the channel to which a region belongs to? I can only think about sampling intensity on the labelingIndex image and compare that with the label and see if those match, but it seems very unelegant to me…

Can someone give me some clues here.

Thanks very much, Julio


#2

I think the issue with your code is that you’re doing the connected-component analysis (cca) across all dimensions, thereby assigning the same label to all objects that “touch” in the channel dimension.

You could try to apply the same strategy of slicing that you used for thresholding, but I’m not sure if the labels would be generated starting from 1 for both slices then, which could lead to different problems (i.e. the same label for disparate objects in the two channels). So you might have to do the looping over regions also separately for both channels then.

If all my assumptions are correct, both of your points 1) and 2) should be addressed with these changes.


You’re using different ways of accessing ops here. When working with scripts, you usually don’t need type safety, and could also use the less type-safe but maybe more convenient notation to call ops by their name:

otsu = ops.op("threshold.otsu", image_in)
imp_thresholded = ops.run("create.img", image_in, BitType())

… which wouldn’t require importing Ops at least.

You shouldn’t need to create Datasets here, as the output processing should be smart enough to find a way to display your outputs. If that’s not the case, maybe there’s a bug that can be fixed…


#3

Hi,

Thanks for the advice. I modified the code to slice through the thresholded image to get the regions but I run into issues when trying to do so.

java.lang.ClassCastException: net.imglib2.view.IntervalView cannot be cast to net.imglib2.roi.labeling.ImgLabeling

Something doesn’t seem to be the right type but I cannot get around this… What am I doing wrong?
I tried to pass an intervalView to the ops.slice in the input and/or output without success

This is the modified part and below the whole code in case needed:

[...]
# create memory for the labeling image. Different options to test
imp_labels = ops.run("create.img", image_in)
intervals = [x_length, y_length, c_length, z_length]
imp_labels_interval = ops.transform().intervalView(imp_labels, [0,0,0,0], intervals)
imp_thresholded_interval = ops.transform().intervalView(imp_thresholded, [0,0,0,0], intervals)

# create the labeling op
labeling_op = ops.op("labeling.cca", imp_thresholded, StructuringElement.FOUR_CONNECTED)

# call the labeling op slice-wise per channel
# NOTE: here I tried either imp_thresholded_interval and imp_labels_interval as input and output
ops.slice(imp_labels, imp_thresholded, labeling_op, [x_dim_index, y_dim_index, z_dim_index])
[...]

and the full script:

#@ OpService ops
#@ DatasetService ds
#@ net.imagej.Dataset image_in

#@ OUTPUT imp_thresholded
#@ OUTPUT labelingIndex

from net.imagej.axis import Axes
from net.imglib2.algorithm.labeling.ConnectedComponents import StructuringElement
from net.imglib2.type.logic import BitType
from net.imglib2.roi import Regions
from net.imglib2.roi.labeling import LabelRegions, ImgLabeling

# Get dimension indexes
x_dim_index = image_in.dimensionIndex(Axes.X)
y_dim_index = image_in.dimensionIndex(Axes.Y)
z_dim_index = image_in.dimensionIndex(Axes.Z)
c_dim_index = image_in.dimensionIndex(Axes.CHANNEL)

# Get dimension lengths
x_length = image_in.dimension(x_dim_index)
y_length = image_in.dimension(y_dim_index)
z_length = image_in.dimension(z_dim_index)
c_length = image_in.dimension(c_dim_index)

# create memory for the thresholded image
imp_thresholded = ops.run("create.img", image_in, BitType())

# create the otsu op
otsu_op = ops.op("threshold.otsu", image_in)

# call threshold op slice-wise per channel
ops.slice(imp_thresholded, image_in, otsu_op, [x_dim_index, y_dim_index, z_dim_index])

# create memory for the labeling image
imp_labels = ops.run("create.img", image_in)
intervals = [x_length, y_length, c_length, z_length]
imp_labels_interval = ops.transform().intervalView(imp_labels, [0,0,0,0], intervals)
imp_thresholded_interval = ops.transform().intervalView(imp_thresholded, [0,0,0,0], intervals)

print 'imp_labels', type(imp_labels)
print 'imp_labels_interval', type(imp_labels_interval)
print 'imp_thresholded', type(imp_thresholded)
print 'imp_thresholded_interval', type(imp_thresholded_interval)

# create the labeling op
labeling_op = ops.op("labeling.cca", imp_thresholded, StructuringElement.FOUR_CONNECTED)

# call the labeling op slice-wise per channel
# NOTE: here I tried either imp_thresholded_interval and imp_labels_interval as input and output
ops.slice(imp_labels, imp_thresholded, labeling_op, [x_dim_index, y_dim_index, z_dim_index])

# get the index image (each object will have a unique gray level)
labelingIndex = imp_labels.getIndexImg()

# get the collection of regions and loop through them
regions = LabelRegions(imp_labels)

for region in regions:
    size = region.size()
    if size > 100:
	    print '----------'
	    print 'Label: ', region.getLabel()
	    print 'Size:', size
	
	    # Calculate the center of mass using region.getCenterOfMass()
	    center_of_mass = region.getCenterOfMass() 
	    print 'CoM x:', center_of_mass.getDoublePosition(x_dim_index)
	    print 'CoM y:', center_of_mass.getDoublePosition(y_dim_index)
	    print 'CoM z:', center_of_mass.getDoublePosition(z_dim_index)
	
	    # Get the mean intensity the region by "sampling" the intensity of the input image at the region pixels
	    intensity = ops.stats().mean(Regions.sample(region, image_in)).getRealFloat()
	    print "intensity:", intensity
	    
	    center_of_gravity_op = ops.op('geom.centerOfGravity', Regions.sample(region, image_in))
	    
	    print 'CoG:', center_of_gravity_op.calculate()

Cheers, Julio