How can I display multiple Images and plot the centroid ?

I have this code which finds the centroid and plots it, how can i do the same but for multiple images and display them?

I have also posted my attempt which is a mess, it shows no errors, but it only displays the last image due to the for loop I'm guesssing.

clear;
close all;
clc;


I = imread('C:\DilatedRedBeanv0.png');
Ibw = im2bw(I);
Ibw = imfill(Ibw,'holes');
Ilabel = bwlabel(Ibw);
stat = regionprops(Ilabel,'centroid');
imshow(I); hold on;
for x = 1: numel(stat)
   plot(stat(x).Centroid(1),stat(x).Centroid(2),'ro');
end


X = (stat(x).Centroid(1));
Y = (stat(x).Centroid(2));

[rows, columns, numberOfColorChannels] = size(Ibw);  %finding center of image
y_center = rows / 2;
x_center = columns / 2;

Pos= [X Y;x_center y_center];
EuclideanDistance=pdist(Pos, 'euclidean');  %Calculates euclidean distance from centroid to the center of image

CityblockDistance=pdist(Pos, 'cityblock'); %Calculates cityblock distance from centroid to the center of image
CityblockFormula = (X-x_center)+(Y-y_center);

Distance= sqrt((x_center-X)^2+(y_center-Y)^2);



%MY ATTEMPT AT PLOTTING THE CENTROID IN  MULTIPLE IMAGES

clear;
close all;
clc;


I = imread('C:\DilatedRedBeanv0.png');%Red Bean
I2 = imread('C:\ClosedGreenBeanv0.png'); % Green Bean
I3 = imread('C:\DilatedWhiteBeanv0.png'); %White Bean
I4 = imread('C:\DilatedYellowBeanv0.png'); %Yellow Bean

%Making the images Binary

Ibw = im2bw(I);
Ibw2 = im2bw(I2);
Ibw3 = im2bw(I3);
Ibw4 = im2bw(I4);


Ibw = imfill(Ibw,'holes');
Ibw2 = imfill(Ibw2,'holes');
Ibw3 = imfill(Ibw3,'holes');
Ibw4 = imfill(Ibw4,'holes');


Ilabel = bwlabel(Ibw);
Ilabel2 = bwlabel(Ibw2);
Ilabel3 = bwlabel(Ibw3);
Ilabel4 = bwlabel(Ibw4);


stat = regionprops(Ilabel,'centroid');
imshow(I); hold on;
for x = 1: numel(stat)
   plot(stat(x).Centroid(1),stat(x).Centroid(2),'ro');
end
stat2 = regionprops(Ilabel2,'centroid');
imshow(I2); hold on;
for x2 = 1: numel(stat2)
   plot(stat2(x2).Centroid(1),stat2(x2).Centroid(2),'ro');
end
stat3 = regionprops(Ilabel3,'centroid');
imshow(I3); hold on;
for x3 = 1: numel(stat3)
   plot(stat3(x3).Centroid(1),stat3(x3).Centroid(2),'ro');
end
stat4 = regionprops(Ilabel4,'centroid');
imshow(I4); hold on;
for x4 = 1: numel(stat4)
   plot(stat4(x4).Centroid(1),stat4(x4).Centroid(2),'ro');
end


X = (stat(x).Centroid(1));
Y = (stat(x).Centroid(2));

X2 = (stat2(x).Centroid(1));
Y2 = (stat2(x).Centroid(2));

X3 = (stat3(x).Centroid(1));
Y3 = (stat3(x).Centroid(2));

X4 = (stat4(x).Centroid(1));
Y4 = (stat4(x).Centroid(2));

[rows, columns, numberOfColorChannels] = size(Ibw); %finding center of image
y_center = rows / 2;
x_center = columns / 2;

[rows, columns, numberOfColorChannels2] = size(Ibw2); 
y2_center = rows / 2;
x2_center = columns / 2;

[rows, columns, numberOfColorChannels3] = size(Ibw3);
y3_center = rows / 2;
x3_center = columns / 2;

[rows, columns, numberOfColorChannels4] = size(Ibw4);
y4_center = rows / 2;
x4_center = columns / 2;

[Rik]

[rescued from spam filter]

Answer 1

[Meg Noah]

I played with the first part to make some more plot options (I had to do the poor man's distances because I can't afford the stats toolbox required to do basic image processing that is not included in the image processing toolbox):

clear all;
close all;
clc;

I = imread('DilatedRedBeanv0.png');
Ibw = im2bw(I);
Ibw = imfill(Ibw,'holes');
Ilabel = bwlabel(Ibw);
stat = regionprops(Ilabel,'centroid');
fig1 = figure('color','white','position',[33 348 680 570]);
imagesc(1:size(I,2),1:size(I,1),I); hold on;
set(gca,'ydir','normal');
axis equal; axis tight;
for x = 1: numel(stat)
   plot(stat(x).Centroid(1),stat(x).Centroid(2),'ro','MarkerFaceColor','r', ...
       'DisplayName',['Dilated Red Bean v0 Centroid = (' ...
       num2str(stat(x).Centroid(1)) ',' ...
       num2str(stat(x).Centroid(2)) ')']);
end
X = (stat(x).Centroid(1));
Y = (stat(x).Centroid(2));
[rows, columns, numberOfColorChannels] = size(Ibw);  %finding center of image
y_center = rows / 2;
x_center = columns / 2;
Pos= [X Y;x_center y_center];
EuclideanDistance = sqrt((X-x_center).^2+(Y-y_center).^2);
CityblockDistance = sum(abs(X-x_center) + abs(Y-y_center));
% EuclideanDistance=pdist(Pos, 'euclidean');  %Calculates euclidean distance from centroid to the center of image
% CityblockDistance=pdist(Pos, 'cityblock'); %Calculates cityblock distance from centroid to the center of image
CityblockFormula = (X-x_center)+(Y-y_center);
Distance= sqrt((x_center-X)^2+(y_center-Y)^2);

% distance to center
plot([x_center X],[y_center Y],'c','DisplayName',['Euclidean Distance = ' num2str(EuclideanDistance)]);
plot([x_center x_center x_center X],[y_center Y Y Y],'g','DisplayName',['Cityblock Distance = ' num2str(CityblockDistance)]);
legend('location','best');

So this function simplifies things:

function [I,Ilabel,stat] = processBeanImage(filename)
I = imread(filename);
%Making the images Binary
Ibw = im2bw(I);
Ibw = imfill(Ibw,'holes');
Ilabel = bwlabel(Ibw);
stat = regionprops(Ilabel,'centroid');
end

Then your processing multiple beans is:

fig2 = figure('color','white','position',[33 0 900 900]);

fileList = {'DilatedRedBeanv0.png';'ClosedGreenBeanv0.png'; ...
    'DilatedWhiteBeanv0.png';'DilatedYellowBeanv0.png';};
for ifile = 1:length(fileList)
    [I,Ilabel,stat] = processBeanImage(fileList{ifile});
    subplot(2,2,ifile)
    imagesc(Ilabel); hold on;
    title([fileList{ifile} ' has ' num2str(numel(stat)) ' centroids']);
    for j = 1: numel(stat)
        plot(stat(j).Centroid(1),stat(j).Centroid(2),'ro');
    end
end