Using Image Processing to calculate Moment of Inertia

I have made code that would (hopefully) have calculated the moment of inertia of different infill density 3D printed beams. This code is not working however. I am not quite sure where to go from here. I've turned the image into a binary matrix and attempted to apply inorder to solve for the moment of intertia with the trapz function.

How could my code be corrected to fix this?

close all; clear all;
set(0,'defaultTextFontName','Times')
set(0,'defaultTextFontSize',10)
set(0,'defaultAxesFontName','Times')
set(0,'defaultAxesFontSize',12)
set(0,'defaultLineLineWidth',1.5)
% Load the image and convert to binary 
Im_50 = imread('50percent.JPG');
B_50 = imbinarize(Im_50);
B1_50 = B_50(:,:,1);
b_50 = 0.0256286; % Width
x_50 = 0:(b_50/length(B_50(1,:,1))):b_50;
h_50 = 0.009652; % Thickness M
y_50 = 0:(h_50/length(B_50(:,1,1))):h_50;
rows = length(B_50(:,1,1));
columns = length(B_50(1,:,1));
deltax = zeros(rows, columns);
deltay = zeros(rows, columns);
midRow = mean([1, rows]);
midCol = mean([1, columns]);
for col = 1 : columns
    for row = 1 : rows
        deltax(row, col) = col - midCol;
        deltay(row, col) = row - midRow;
    end
end
deltay = (deltay/deltay(1,1)).*h_50;
deltax = (deltax/deltax(1,1)).*b_50;
% Matix with zeros where there is no fill and distance from origin where there is fill
hy_50 = deltay.*B1_50(:,:,1);
bx_50 = deltax.*B1_50(:,:,1);
I_50 = trapz(hy_50.^2,bx_50);

Answer 1

[Ayush Modi]

Hi,

As per my understanding, you would like to calculate the moment of inertia of different infill density 3D printed beams. You can achieve this by summing up the contributions of all the pixels that have material and multiplying by the square of their distance from the neutral axis (centroid).

Here is an example showing how you can calculate the moment of inertia:

% Number of rows and columns in the binary image
rows = size(B1_50, 1);
columns = size(B1_50, 2);

% Calculate pixel size in each dimension
dx = b_50 / columns;
dy = h_50 / rows;

% Create matrices for the x and y coordinates of each pixel
[x, y] = meshgrid(dx/2:dx:b_50-dx/2, dy/2:dy:h_50-dy/2);

% Calculate the centroid of the binary image
area = sum(B1_50, 'all');
x_centroid = sum(sum(x .* B1_50)) / area;
y_centroid = sum(sum(y .* B1_50)) / area;

% Calculate the distance of each pixel to the centroid
x_dist = x - x_centroid;
y_dist = y - y_centroid;

% Calculate the second moment of area (moment of inertia) for each row/column
I_xx_per_row = sum((y_dist.^2) .* B1_50, 2) * dx; % Moment of inertia contribution per row
I_yy_per_col = sum((x_dist.^2) .* B1_50, 1) * dy; % Moment of inertia contribution per column

% Integrate the moment of inertia contributions using trapz
I_xx = trapz(I_xx_per_row) * dy;
I_yy = trapz(I_yy_per_col) * dx;

% Display the calculated moments of inertia
disp(['I_xx = ', num2str(I_xx), ' m^4']);
disp(['I_yy = ', num2str(I_yy), ' m^4']);

This script assumes that the beam is symmetric about the x-axis and that the neutral axis is at the centre of the beam's height. If the neutral axis is not at the centre, you will need to find the centroid of the shape first and then calculate the moment of inertia about that axis.

I hope this resolves the issue you were facing.