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:
- 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?
- 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