Fit ellipse to objects in image

Hello, I am trying to fit ellipses to many different objects within a binary images. The ultimate goal is calculate the average eccentricity of the objects, but I can currently only plot a boundary around the objects. Any suggestions would be much appreciated. Image below (metallic sample with visible porosity).

clear
clc
close all
[f,p]=uigetfile('*.tif','select image file');
file_name=strcat(p,f)
BW = imbinarize(rgb2gray(imread(file_name)),0.68);
[B,L,N,A] = bwboundaries(BW);
figure; imshow(BW); hold on;
% Loop through object boundaries

for k = 1:N
    % Boundary k is the parent of a hole if the k-th column
    % of the adjacency matrix A contains a non-zero element
    if (nnz(A(:,k)) > 0)
        boundary = B{k};
        % Loop through the children of boundary k
        area1=polyarea(boundary(:,2), boundary(:,1));
        for l = find(A(:,k))'
                boundary = B{l};
                plot(boundary(:,2),...
                boundary(:,1),'g','LineWidth',2);
            areaholes(l)=polyarea(boundary(:,2), boundary(:,1));
            
        end
    end
end

Answer 1

[Kevin Holly]

I would suggest using regionprops. Steve Eddins walks through the process on his blog.

https://blogs.mathworks.com/steve/2015/08/17/ellipse-visualization-and-regionprops/

[Kevin Holly]

regionprops also can provide the eccentricity value.

[Sean Dobson]

Thank you for your help!

Answer 2

[Matt J]

T=regionprops('table',~BW,'MajorAxisLength','MinorAxisLength','Eccentricity');

meanEccentricity=mean(T.Eccentricity);

Areas=pi*T.MajorAxis.Length*T.MinorAxisLength/4;

[Sean Dobson]

Thank you for your help!

Answer 3

[Matt J]

You cna use ellipticalFit(), distributed here:

https://www.mathworks.com/matlabcentral/fileexchange/32108-fitting-a-conic-to-a-given-set-of-points-using-levenberg-marquardt-method

BW0 = imbinarize(rgb2gray(imread('https://www.mathworks.com/matlabcentral/answers/uploaded_files/856545/image.jpeg')),0.68);

BW=bwareaopen(~BW0,15);
B=bwboundaries(BW,'noholes');
figure; imshow(BW0); hold on;
% Loop through object boundaries

N=numel(B);

for k = 1:N

    efit=ellipticalFit(flipud(B{k}.'));

    hold on
    fimplicit(efit,'Color','r','LineWidth',3);
    hold off

    areaholes(k)=pi*efit.a*efit.b/4;

end

axis([0.0070    0.7170    0.7402    1.2728]*1000)