Clear Filters
Clear Filters

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

24 views (last 30 days)
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 ?

Accepted Answer

DUY Nguyen
DUY Nguyen on 3 Mar 2023
Edited: DUY Nguyen on 3 Mar 2023
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');
  5 Comments
Selina
Selina on 1 Dec 2023
Edited: Selina on 1 Dec 2023
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
navin chandra on 3 Dec 2023
Edited: navin chandra on 3 Dec 2023
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.

Sign in to comment.

More Answers (0)

Products


Release

R2020b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!