finding the axis for least moment of inertia of an object in 2D binary image

I have a 2D binarized image as shown below (with background black and a white object of interest).

Now I can easily find the centroid of the white object using 'regionprops'. I want to know how to find an axis passing through the centroid which gives least area moment of inertia of the white object ?

Answer 1

[DUY Nguyen]

Hi,

Firstly, you invert the image. Then you can calculate the principal axis passing through the centroid by this eq:

tan(2θ) = 2mxy / (mx^2 - my^2)

You can try this code below:

% Read the binary image and invert it
im = imread('binary_image.png');
im = imcomplement(im);

% Calculate the region properties of the object
stats = regionprops(im, 'centroid', 'area', 'majoraxislength', 'minoraxislength');

% Get the centroid coordinates
centroid = stats.Centroid;

% Get the major and minor axis lengths
major_axis = stats.MajorAxisLength / 2;
minor_axis = stats.MinorAxisLength / 2;

% Calculate the moments of the object
 Ixx = sum(sum((repmat((1:size(im,1))', 1, size(im,2)) - centroid(2)).^2 .* im));
Iyy  = sum(sum((repmat(1:size(im,2), size(im,1), 1) - centroid(1)).^2 .* im));
Ixy = sum(sum((repmat((1:size(im,1))', 1, size(im,2)) - centroid(2)) .* ...
    (repmat(1:size(im,2), size(im,1), 1) - centroid(1)) .* im));

% Calculate the angle of the principal axis
theta = atan2(2*Ixy, Iyy-Ixx)/2;

% Calculate the coordinates of the two points on the principal axis
x1 = centroid(1) + major_axis*cos(theta);
y1 = centroid(2) - major_axis*sin(theta);
x2 = centroid(1) - major_axis*cos(theta);
y2 = centroid(2) + major_axis*sin(theta);

% Plot the principal axis on top of the image
imshow(im);
hold on;
plot([x1, x2], [y1, y2], 'LineWidth', 2, 'Color', 'red');

[navin chandra]

Thanks for this answer. To check this code, I ran it through a simple object (perfect ellipse shape). The result obtained from the above code is shown in the figure below

With the red line being the axis identified by the code. But, I think this is not correct because, for an ellipse, the axis with least area moment of inertia should coincide with the major axis of the ellipse. But the above method calculates something else.

[navin chandra]

Okay, I figured out the problem. In the above code Ixx, Iyy and Ixy are being calculated about the axis passing through the origin (x=0,y=0). However, all these axes should pass from the centroid of the object. Anyways thanks for the approach. Here is the corrected code for anyone in need of similar thing-

% Read the binary image and invert it

im = imread('binary_image.tif');

im = imbinarize(im);

% im = imcomplement(im);

% Calculate the region properties of the object

stats = regionprops(im, 'centroid', 'area', 'majoraxislength', 'minoraxislength');

% Get the centroid coordinates

centroid = stats.Centroid;

% Get the major and minor axis lengths

major_axis = stats.MajorAxisLength / 2;

minor_axis = stats.MinorAxisLength / 2;

% Calculate the moments of the object

Ixx = sum(((1:size(im,1))'-centroid(2)).^2 * ones(1,size(im,2)) .* im, 'all');

Iyy = sum(ones(size(im,1),1) * ((1:size(im,2))-centroid(1)).^2 .* im, 'all');

Ixy = sum(repmat(((1:size(im,1))'-centroid(2)),1,size(im,2)) .* repmat((1:size(im,2))-centroid(1),size(im,1),1) .* im, 'all');

% Calculate the angle of the principal axis

theta = atan2(2*Ixy, Iyy-Ixx)/2;

% Calculate the coordinates of the two points on the principal axis

x1 = centroid(1) + major_axis*cos(theta);

y2 = centroid(2) - major_axis*sin(theta);

x2 = centroid(1) - major_axis*cos(theta);

y1 = centroid(2) + major_axis*sin(theta);

% Plot the principal axis on top of the image

imshow(im);

hold on;

plot([x1, x2], [y1, y2], 'LineWidth', 2, 'Color', 'red');

And the prediction for an ellipse is shown below-

[DUY Nguyen]

Thanks @navin chandra, I also updated the above code. You're correct!

[Selina]

Hey! I tried running the code using an ellipse (I created the shape in PowerPoint and added a black background). I then added a line in the code to change the png to a binary image and ran it, I added the code with the extra 2 lines at the end. However, it does not seem to pick up the axis correctly. Any suggestions? It seems to pick up the centroid of the white background rather than the ellipse which makes me think, the axis is for the white rather than the black?

% Read the binary image and invert it

I = imread('image_2.png');

im = imbinarize(I(:,:,3));

% im = imread('binary_image.tif');

% im = imbinarize(im);

im = imcomplement(im);

% Calculate the region properties of the object

stats = regionprops(im, 'centroid', 'area', 'majoraxislength', 'minoraxislength');

% Get the centroid coordinates

centroid = stats.Centroid;

% Get the major and minor axis lengths

major_axis = stats.MajorAxisLength / 2;

minor_axis = stats.MinorAxisLength / 2;

% Calculate the moments of the object

Ixx = sum(((1:size(im,1))'-centroid(2)).^2 * ones(1,size(im,2)) .* im, 'all');

Iyy = sum(ones(size(im,1),1) * ((1:size(im,2))-centroid(1)).^2 .* im, 'all');

Ixy = sum(repmat(((1:size(im,1))'-centroid(2)),1,size(im,2)) .* repmat((1:size(im,2))-centroid(1),size(im,1),1) .* im, 'all');

% Calculate the angle of the principal axis

theta = atan2(2*Ixy, Iyy-Ixx)/2;

% Calculate the coordinates of the two points on the principal axis

x1 = centroid(1) + major_axis*cos(theta);

y2 = centroid(2) - major_axis*sin(theta);

x2 = centroid(1) - major_axis*cos(theta);

y1 = centroid(2) + major_axis*sin(theta);

% Plot the principal axis on top of the image

imshow(im);

hold on;

plot([x1, x2], [y1, y2], 'LineWidth', 2, 'Color', 'red');

[navin chandra]

Yes @Selina, you are right. The code is written considering white pixels as the object of interest and black pixels as the background. The image you have posted is having the ellipse in black and the background in white, if you just invert that image then the code will work fine.

I can also see that you have used "im = imcomplement(im)", in your code. I think that if you skip this step then rest of the code will work fine.