Area of a particle in Imagej (actual particle area or area of the best fit ellipse?)

Dear all,
I am new to this forum and I would like to pose a question on *how to measure the area of particles. Hopefully you can help me and other readers to clarify the meaning of some particle’s parameters.

My goal is to calculate the area of each individual particle of a given image (attached is an example). After having obtained a binarized image with white as the background and black as the particles, I use the command “Analyze Particles” to count and measure the particles. In the result table, ImageJ reports the calculated parameters for each particle, including the area and the major and minor axis of the fit ellipse. Knowing the axes of the fit ellipse I calculated the area of the fit ellipse for each particle, to then find out that the area of the fit ellipse calculated by Imagej is equal to the area of the particle.

My questions are: 1) What area is ImageJ calculating? The actual area of the particle or the area of its best fit ellipse? 2) Does ImageJ draw the best fit ellipse of an object, so that the area of the fit ellipse is equal to the area of the object?

For my research it is interesting to know both the area of a particle and the area of its best fit ellipse to then calculate the difference between the two area values (which is something similar to the solidity parameter where the area of the particle is compared to the area of the particle’s convex hull).

Thanks in advance for your time
Oscar

Hi @Oscar_Amicis,

Welcome to the forum!

  1. The area that ImageJ calculates is the area of the particle itself. If the image is spatially calibrated, it will provide the area in physical units. Otherwise, the area is just the count of pixels comprising the object.
  2. In fact, that’s the definition it is using. It finds an ellipse that has the same area of the object while maintaining the same second order central moments (see this post and follow the link to the source code for more info).

In this case, the definition of “best fit” that ImageJ uses precludes the use of that metric, since both areas are always the same.

What definition of “best fit ellipse” did you have in mind?

Cheers,
Nico

1 Like

Hi @NicoDF
Thank you very much for your quick answer! That helped.
ImageJ applies then a fitting method so that the fit ellipse is scaled to have the same area of the original particle. There can be different ways to fit an ellipse to an irregular shaped object. For my purposes, the smallest ellipse that circumscribes the object would do. I found some useful information in the link you recommended about this topic.
Cheers

As discussed in that topic, that problem is harder than it seems. Out of curiosity, I’ve been doing some exploration for that particular question, and the problem boils down to two alternatives, which revolve around the convex hull of the object: either 3 points of the hull provide a minimum ellipse which encloses all of the other points, or no point triplet can provide such an ellipse, in which case you need a 4th point to define it. I have a macro that solves the first case, but I haven’t tackled the second case yet. I can post that partial solution if it is of help to you.

Cheers,
Nico

That’d be great! Thanks for exploring the problem.

it is very complicated indeed. Do you think it would be possible to use the minimum bounding box to draw an ellipse having as axes the length and the width of the rectangle? For my investigation that would also be a good way of having a consistent way for creating the fit ellipse to which compare the irregular shaped objects. In this case, is it possible to obtain the actual smallest bounding box of an object? Having a look here Bounding smallest tilted rectangle, Imagej creates a bounding box aligned to the x and y axes, which precludes to have the smallest bounding box due to orientation problems

Sorry for the delay, @Oscar_Amicis.

That’s what I think. I believe the fit rectangle with minimum area is the one aligned with the minimum Feret diameter (any reference that backs up my intuition would be helpful here).

Update: scratch that! I posted a corrected version below in the thread.



I wrote a function that creates such a rectangle. Here’s a example implementing it in the image you posted:

run("Select None");
if (roiManager("count")>0){
	roiManager("deselect");
	roiManager("delete");	
	}

setAutoThreshold("Default");
run("Analyze Particles...", "size=100-Infinity add");
resetThreshold();
roiManager("Show All without labels");

nrois = roiManager("count");
for (i = 0; i < nrois; i++) {
	roiManager("select", i);
	fitMinRectangle();
	run("Properties... ", "  stroke=orange");
	roiManager("add");
}

run("Select None");

// ******************* functions ***************************

function fitMinRectangle(){
	// fits a minimum area rectangle into a ROI, based on the minimum Feret diameter
	run("Convex Hull");
	getSelectionCoordinates(xp, yp);
	run("Undo"); // until run("Restore Selection"); for convex hull gets into production (ver >= 1.52u18)
	
	np = xp.length;
	
	minFD = getValue("Width") + getValue("Height");
	imin = -1;
	i2min = -1;
	jmin = -1;
	for (i = 0; i < np; i++) {
		maxLD=0;
		imax = -1;
		i2max = -1;
		jmax = -1;
		if(i<np-1) i2 = i + 1; else i2 = 0;
		for (j = 0; j < np; j++) {
			d = abs(perpDist(xp[i], yp[i], xp[i2], yp[i2], xp[j], yp[j]));
			if (maxLD < d) {
				maxLD = d;
				imax = i;
				jmax = j;
				i2max = i2;
				}
			}
		if (minFD > maxLD){
			minFD = maxLD;
			imin = imax;
			i2min = i2max;
			jmin = jmax;
			}
		}
	
	pd = perpDist(xp[imin], yp[imin], xp[i2min], yp[i2min], xp[jmin], yp[jmin]); //signed minFeret
	pairAngle = atan2( yp[i2min]- yp[imin], xp[i2min]- xp[imin]);
	minAngle = pairAngle + PI/2;
	
	hmin = 0;
	hmax = 0;
	
	for (i = 0; i < np; i++) {
		hd = parDist(xp[imin], yp[imin], xp[i2min], yp[i2min], xp[i], yp[i]);
		hmin = minOf(hmin, hd);
		hmax = maxOf(hmax, hd);
		}
	
	nxp=newArray(4);
	nyp=newArray(4);
	
	nxp[0] = xp[imin] + cos(pairAngle) * hmax;
	nyp[0] = yp[imin] + sin(pairAngle) * hmax;
	
	nxp[1] = nxp[0] + cos(minAngle) * pd; 
	nyp[1] = nyp[0] + sin(minAngle) * pd;
	
	nxp[2] = nxp[1] + cos(pairAngle) * (hmin-hmax); 
	nyp[2] = nyp[1] + sin(pairAngle) * (hmin-hmax);
	
	nxp[3] = nxp[2] + cos(minAngle) * - pd; 
	nyp[3] = nyp[2] + sin(minAngle) * - pd;
	
	makeSelection("polygon", nxp, nyp);
	}

function dist2(x1,y1,x2,y2){
	return pow(x1-x2, 2) + pow(y1-y2, 2);
	}

function perpDist(p1x, p1y, p2x, p2y, x, y){
	// signed distance from a point (x,y) to a line passing through p1 and p2
	return ((p2x - p1x)*(y - p1y) - (x - p1x)*(p2y - p1y))/sqrt(dist2(p1x, p1y, p2x, p2y));
	}
	
function parDist(p1x, p1y, p2x, p2y, x, y){
	// signed projection of vector (x,y)-p1 into a line passing through p1 and p2
	return ((p2x - p1x)*(x - p1x) + (y - p1y)*(p2y - p1y))/sqrt(dist2(p1x, p1y, p2x, p2y));
	}

Let me know if it works for your goals.

Cheers,
Nico

2 Likes

Hello Nico -

As a geometrical aside, the minimum-area rectangle need not be
aligned with the minimum Feret diameter. I’ve created an annoying
shape (“pesky_shape”) as a counter-example.

Below is an illustrative script. In it I create the pesky_shape
Roi. Its conventional bounding box (the bounding rectangle
aligned with the x and y axes) is the Feret-aligned rectangle.
I then create an approximation to the actual minimum-area
rectangle.

For reasons I don’t really understand, when I run your
fitMinRectangle() IJM function on the pesky_shape Roi,
I don’t get the x-y-aligned, Feret-aligned rectangle, but one
that is rotated clockwise a bit. (I suspect this has to do with
some artifacts introduced by the ShapeRoi conversion in the
Roi “Combine” operation I use.)

I measure the “Area” of the three rectangles, and, indeed, the
approximate minimum-area rectangle has the smallest area of
the three. These Areas are:

   naive bounding box:       8000
   minimum-area rectangle:   7966
   fitMinRectangle():        8020

Here is the jython script:

from ij import IJ
from ij.gui import EllipseRoi
from ij.gui import PolygonRoi
from ij.gui import Roi
from ij.plugin.frame import RoiManager

imp = IJ.createImage ('pesky_shape', '8-bit ramp', 512, 512, 1)
imp.getProcessor().multiply (0.5)
imp.getProcessor().add (63)
imp.show()

rm = RoiManager()

x0 = 256.0
y0 = 256.0
scale = 10.0

w = 0.2
l = 10.0
o = 1.0
er = 2.0
ar = 1.0 / 1.1

x1 = 0.0
y1 = -er / ar
x2 = 0.0
y2 = -y1

x1 = x0 + scale * x1
y1 = y0 + scale * y1
x2 = x0 + scale * x2
y2 = y0 + scale * y2

ell = EllipseRoi (x1, y1, x2, y2, ar)

xr1 = o - w
yr1 = -l
xr2 = -xr1
yr2 = 0.0
width = w
height = l

xr1 = x0 + scale * xr1
yr1 = y0 + scale * yr1
xr2 = x0 + scale * xr2
yr2 = y0 + scale * yr2
width *= scale
height *= scale

r1 = Roi (xr1, yr1, width, height)
r2 = Roi (xr2, yr2, width, height)

rm.addRoi (ell)
rm.addRoi (r1)
rm.addRoi (r2)

# create "pesky_shape" roi by combining ell, r1, and r2
rm.runCommand ('Select All')
rm.runCommand ('Combine')
rm.runCommand ('Delete')
rm.runCommand ('Add')

# create pesky_shape's naive bounding box
xrb = x0 - scale * er
yrb = y0 - scale * l
wrb = scale * 2.0 * er
hrb = scale * 2.0 * l
rb = Roi (xrb, yrb, wrb, hrb)
rm.addRoi (rb)

# create approximate minimum-area rectangle for pesky_shape
l2 = 9.91
er2 = 2.01
xrb2 = x0 - scale * er2
yrb2 = y0 - scale * l2
wrb2 = scale * 2.0 * er2
hrb2 = scale * 2.0 * l2
rb2 = Roi (xrb2, yrb2, wrb2, hrb2)
imp.setRoi (rb2)
IJ.run ('Rotate...', '  angle=-5.50')
rm.runCommand ('Add')

# save as overlay
IJ.run ('Add Selection...')

# PolygonRoi from running Nico's fitMinRectangle() on pesky_shape
xp = [232.05567932128906, 271.9927673339844, 279.9474792480469, 240.01039123535156]
yp = [355.3670959472656, 356.952392578125, 156.5536346435547, 154.96835327148438]
poly = PolygonRoi (xp, yp, Roi.POLYGON)
rm.addRoi (poly)

# measure naive bounding box
rm.select (1)
rm.runCommand ('Measure')   # Area = 8000
# measure approximate minimum-area rectangle
rm.select (2)
rm.runCommand ('Measure')   # Area = 7966
# measure fitMinRectangle() result
rm.select (3)
rm.runCommand ('Measure')   # Area = 8020

# (duplicate and) save as png with pesky_shape and minimum-area rectangle
rm.runCommand ('Deselect')
impdup = imp.duplicate()
impdup.show()
rm.select (0)
IJ.run ('Add Selection...')
IJ.saveAs (impdup, 'PNG', 'pesky_shape')
impdup.close()

# show pesky_shape roi inside of minimum-area rectangle overlay
rm.select (0)

And here is an illustrative image produced by the script:

Thanks, mm

2 Likes

Hey @mountain_man!

Thank you so much for looking into this! You’re totally right! (and I was totally wrong, which I love!). There was certainly room from improvement :smiley:.

I could not find a counter-example to my logic, but it bugged me not to be able to provide a clear proof. I was planning to brute-force some examples during this weekend to see if it actually held, but there’s no need to.

I did my due dilligence (sorry for the initial lazyness), and indeed, the right solution does not need to align to the minimum Feret diameter:

H. Freeman and R. Shapira. 1975. Determining the minimum-area encasing rectangle for an arbitrary closed curve. Commun. ACM 18, 7 (July 1975), 409–413. DOI:https://doi.org/10.1145/360881.360919

The minimum area rectangle has a side that lies on one of the edges of the convex hull.

Fortunately, the right solution is closely based on the algorithm I used to find the minFeret in the first place, so all I had to do is rearrange my code to look for the actual minimum area. It adds a little overhead, since it has perform the rotating calipers routine for each side of convex hull (instead of a single time for minFeret), but it works ok:

Min area: 7910

So, here’s the link to the updated function: https://github.com/ndefrancesco/macro-frenzy/blob/master/geometry/fitting/fitMinRectangle.ijm

And here’s the corrected code for @Oscar_Amicis:

run("Select None");
if (roiManager("count")>0){
	roiManager("deselect");
	roiManager("delete");	
	}

setAutoThreshold("Default");
run("Analyze Particles...", "size=100-Infinity add");
resetThreshold();
roiManager("Show All without labels");

nrois = roiManager("count");
for (i = 0; i < nrois; i++) {
	roiManager("select", i);
	fitMinRectangle();
	run("Properties... ", "  stroke=orange");
	roiManager("add");
}

run("Select None");

// ******************* functions ***************************

function fitMinRectangle(){
	/*	Fits a minimum area rectangle into a ROI.
	 * 	
	 * 	It searches the for the minimum area bounding rectangles among the ones that have a side 
	 * 	which is collinear with an edge of the convex hull.
	 *	
	 * 	Concept based on:
	 * 	H. Freeman and R. Shapira. 1975. Determining the minimum-area encasing rectangle for an arbitrary 
	 * 	closed curve. Commun. ACM 18, 7 (July 1975), 409–413. DOI:https://doi.org/10.1145/360881.360919
	 * 	
	 * 	Many thanks to @mountain_man (https://forum.image.sc/u/mountain_man) for helping me correct
	 * 	my initial (wrong) intuition about the right way to tackle this problem!
	 * 	
	 */

	run("Convex Hull");
	getSelectionCoordinates(xp, yp);
	run("Undo"); // until run("Restore Selection"); for convex hull gets into production (ver >= 1.52u18)
	
	np = xp.length;

	minArea = 2 * getValue("Width") * getValue("Height");
	minFD = getValue("Width") + getValue("Height"); // FD now stands for first diameter :)
	imin = -1;
	i2min = -1;
	jmin = -1;
	for (i = 0; i < np; i++) {
		maxLD = 0;
		imax = -1;
		i2max = -1;
		jmax = -1;
		if(i<np-1) i2 = i + 1; else i2 = 0;
		
		for (j = 0; j < np; j++) {
			d = abs(perpDist(xp[i], yp[i], xp[i2], yp[i2], xp[j], yp[j]));
			if (maxLD < d) {
				maxLD = d;
				imax = i;
				jmax = j;
				i2max = i2;
				}
			}
			
		hmin = 0;
		hmax = 0;
		
		for (k = 0; k < np; k++) { // rotating calipers
			hd = parDist(xp[imax], yp[imax], xp[i2max], yp[i2max], xp[k], yp[k]);
			hmin = minOf(hmin, hd);
			hmax = maxOf(hmax, hd);
			}
		
		area = maxLD * (hmax - hmin);
		
		if (minArea > area){
			
			minArea = area;
			minFD = maxLD;
			min_hmin = hmin;
			min_hmax = hmax;

			imin = imax;
			i2min = i2max;
			jmin = jmax;
			}
		}
	
	pd = perpDist(xp[imin], yp[imin], xp[i2min], yp[i2min], xp[jmin], yp[jmin]); // signed feret diameter
	pairAngle = atan2( yp[i2min]- yp[imin], xp[i2min]- xp[imin]);
	minAngle = pairAngle + PI/2;

	
	
	nxp=newArray(4);
	nyp=newArray(4);
	
	nxp[0] = xp[imin] + cos(pairAngle) * min_hmax;
	nyp[0] = yp[imin] + sin(pairAngle) * min_hmax;
	
	nxp[1] = nxp[0] + cos(minAngle) * pd; 
	nyp[1] = nyp[0] + sin(minAngle) * pd;
	
	nxp[2] = nxp[1] + cos(pairAngle) * (min_hmin - min_hmax); 
	nyp[2] = nyp[1] + sin(pairAngle) * (min_hmin - min_hmax);
	
	nxp[3] = nxp[2] + cos(minAngle) * - pd; 
	nyp[3] = nyp[2] + sin(minAngle) * - pd;
	
	makeSelection("polygon", nxp, nyp);
	}

function dist2(x1,y1,x2,y2){
	return pow(x1-x2, 2) + pow(y1-y2, 2);
	}

function perpDist(p1x, p1y, p2x, p2y, x, y){
	// signed distance from a point (x,y) to a line passing through p1 and p2
	return ((p2x - p1x)*(y - p1y) - (x - p1x)*(p2y - p1y))/sqrt(dist2(p1x, p1y, p2x, p2y));
	}
	
function parDist(p1x, p1y, p2x, p2y, x, y){
	// signed projection of vector (x,y)-p1 into a line passing through p1 and p2
	return ((p2x - p1x)*(x - p1x) + (y - p1y)*(p2y - p1y))/sqrt(dist2(p1x, p1y, p2x, p2y));
	}

which now yields:

I’ve updated the post to point to this correction.

Thanks again for your interest and time invested! I really love this kind of exchanges!

Cheers!
Nico

2 Likes

Hello Nico -

Thanks for the old-school reference. Like the song goes, everything
old is new again!

Thanks, mm

1 Like

Hi @NicoDF,
I really appreciate your help and interest in this, many thanks again! It’s a great solution for my study. I will work on it and let you know if everything works fine.

Thank you @mountain_man for the geometrical insight!

cheers
OA

2 Likes

Well, this escalated quickly:

"Fit Rectangle" is now a command on ImageJ/Fiji :blush:

You’ll have to update ImageJ to the latest daily build (1.52u24) to give it a try.

As always, huge thanks to @Wayne for the support.

Cheers!
Nico

5 Likes

Very good @NicoDF and @mountain_man
Mathieu

1 Like

Hello Nico (and Wayne) -

Well, looky that … Thanks Wayne!

Thanks, mm

1 Like

I am not sure if I am doing something wrong, but the fit rectangle does not tilt the box to fit the object I have selected.

Also, is there a way to print out the lengths of the box created in NicoDF’s macro? The box it creates is exactly what I was looking for but I need to know the dimensions of the box, not just the area.

This thread has been very helpful in my project, thank you.

Hi @JD_Welch,

That’s not the way it works: this command fits a (possibly tilted) rectangle on a previouly defined ROI.
So, for example, you create a ROI by using the magic wand on this image, and then run fit rectangle to obtain the tilted (minimum area) rectangle. If you run “Measure” on that rectangle, you’ll get the parameters RRLength and RRwidth, which caracterize its shape:

Furthermore, if you run “Duplicate” with the Roi active, you’ll get a rotated duplicate, which can be useful for neatly tiling detected objects:

imagen

Cheers!
Nico

3 Likes

Hi Nico,

I just came across the Fit Rectangle function and find it very helpful!
However, the duplicate is not rotated when I use this function so it always contains image areas outside of the selected ROI. Do you know what could be the reason and how to fix this?

Edit: It is working now so no more need to answer!

Cheers,
Anne

1 Like