Is the PrincipalAxisLength measurement equivalent to the object diameter in the regionprops3 output?

Question

I am using the regionprops3 command to compute properties of an object in MATLAB. I thought PrincipalAxisLength is the diameter of my object but it is not.

I created a binary image I containing an ellipsoid with radii (40, 10, 10). I get these values from regionprops3:

sum(I(:)) = 16741
stats.volume = 16741
stats.EquivDiameter = 31.739  % 4/3*pi*(31.739/2)^3=16741
stats.PrincipalAxisLength = [71.535, 17.936, 17.908]

If PrincipalAxisLength is the 3 diameters of an ellipsoid then it's volume is:

4/3*pi*71.535/2*17.936/2*17.908/2 = 12031

Which is not equivalent to above volume.

So, how can I compute the radii of the ellipsoid?

Answer 1

The documentation for the PrincipalAxisLength shape measurement describes it as (emphasis mine):

Length (in voxels) of the major axes of the ellipsoid that have the same normalized second central moments as the region, returned as 1-by-3 vector. regionprops3 sorts the values from highest to lowest.

In other words, it fits a multivariate probability distribution to the region and calculates the central moments for that distribution, returning the second central moment (i.e. variance) for the PrincipalAxisLength measurement.

The volume discrepancy you're seeing is due to the fact that these central moments don't define a bounding ellipsoid fit to the data in the region, just the variance of a probabilistic fit to the data distributed within it. For a case where your region contains an ellipsoid shape, it will underestimate the extent of the ellipsoid (and thus the total volume). Here's some code to visualize this:

% Create 3D binary data containing ellipsoid:
[X, Y, Z] = meshgrid(-50:50, -20:20, -20:20);
R = [40 10 10];
I = ((X./R(1)).^2 + (Y./R(2)).^2 + (Z./R(3)).^2 <= 1);

% Calculate statistics with regionprops3:
stats = regionprops3(I, 'all');
L = stats.PrincipalAxisLength./2;

% Create ellipsoid surface from second central moments and plot:
[x, y, z] = ellipsoid(0, 0, 0, L(1), L(2), L(3), 100);
surf(x, y, z, 'FaceAlpha', 0.5, 'FaceColor', 'r', 'EdgeColor', 'none');
axis equal;
hold on;

% Create image of middle slice through ellipsoid and plot:
imageSlice = 256.*uint8(I(:, :, 21));
imageSlice = cat(3, imageSlice, imageSlice, imageSlice);
image(-50:50, -20:20, imageSlice);
view(0, 90);

And here's the final plot you get:

If you're wanting to fit an ellipsoid surface to your binary object, you could probably use the PrincipalAxisLength measurement as a starting guess for the true axes lengths. There's this function on the File Exchange that does ellipsoid fitting, but that uses a set of surface points as opposed to a binary volume. Perhaps you can adapt that code to fit your needs.

Update

The link provided by Cris in a comment suggests that a clear mathematical relationship could be determined for how the second central moments and the ellipse radii are related. Although I haven't had a chance to work through the math yet, I noticed that simply scaling PrincipalAxisLength by sqrt(5/16) gives values tantalizingly close to R:

>> sqrt(5/16).*stats.PrincipalAxisLength

ans =

  40.008372204885049   9.970908565527971   9.970908565527971

Answer 2

I was facing the same problem some time ago for identifying ellipsoids corresponding to binary regions. By developing the maths it is possible retrieve a factor of sqrt(5) for converting the square root of eigen values to ellipsoid radiusses.

The regionprops3 function from Matlab uses a factor of 4, therefore the factor sqrt(5/16) fixes the difference.

The maths involve integration using spherical coordinates, and linearisation of trigonometric expression. As I could not find any other document or resource explaining it, I propose to write it here (hoping this is adequate...).

Development

Definitions

Let a>b>c be the three radiusses of the ellipsoid. We can assume the ellipsoid is centered and aligned with the three main axes, the longest radius along x, the smallest one along z. In that case, the eigen values correspond to the three second-order centered moments lambda_1 = mu200, lambda_2 = mu020, and lambda_2 = mu002.

The idea is then to compute explicitely these values from the parameters of the ellipsoid.

For mu200, the definition is:

mu200 = \int \int \int I(x,y,z) x^2 dx dy dz

with I(x,y,z) being 1 within the ellipsoid, and 0 otherwise.

Integration using spherical coordinates

For the integration, it is more convenient to use spherical coordinates (rho, theta, phi), corresponding to the radius, inclination with the vertical, and azimuth. They corresponds to:

x = rho * a * cos(phi) * sin(theta)
y = rho * b * sin(phi) * sin(theta)
z = rho * c * sin(theta)

The change of coordinates requires the use of the Jacobian of the transform. In that case, its determinant is

a * b * c * rho^2 * sin(theta)

The moment integral for mu200 is then

mu200 = a * b *c \int_0^{2*pi} \int_0^pi \int_0^1 (a * rho * cos(phi) * sin(theta) )^2 * rho^2 * sin(theta) drho dtheta dphi

Can be simplified into

mu200 = ( a^3 * b * c / 5 ) \int_0^{2*pi} \int_0^pi (sin(theta))^3 dtheta (cos(phi))^2 dphi

The integral over theta can be developed using the linearisation of sin^3 theta:

sin^3(theta) = (3/4) * sin(theta) - (1/4) * sin(theta)

Then the inner integral becomes:

\int_0^pi (sin(theta))^3 dtheta = 4/3

Incorporating into the full integral:

mu200 = ( a^3 * b * c / 5 ) * (4/3) \int_0^{2*pi} (cos(phi))^2 dphi
mu200 = ( a^3 * b * c / 5 ) * (4/3) * (1/2) 2*pi
mu200 = (4 * pi / 3) * (a^3 * b * c / 5)

Identification

Denoting by m000 = (4 * pi / 3) * a * b * c the volume of the ellipsoid, one gets:

mu200 = (a^2 / 5) * V

In regionprops3 function, the covariance matrix is obtained by normalizing by the volume.

We therefore have:

lambda_1 = a^2 / 5

and the relation:

a = sqrt(5) * sqrt(lambda1)