Obtaining the best ellipse fit of an irregular shape

Question

Say that we have an irregular shape. How can we obtain the best ellipse fit of such shape? Does regionprops work for that? If so, can you kindly guide me on that?

Thanks.

Answer 1

Ok here is a bit of code to get you started. This is by no means perfect but it was fun doing it and I hope it will help you. The most important part of the code was actually written by @Amro in this great answer: check me!. I incorporated Amro's code to plot the ellipse (Thanks Amro!).

As I mentioned in my comment, you can use regionprops with a couple of parameters to estimate an ellipse that fit some shape. The names are self-explanatory and can be used to describe an ellipse. Here is my dummy shape which looks like a potato:

The parameters we will use in regionprops are:

• 'Centroid'
• 'MajorAxisLength'
• 'MinorAxisLength'
• 'Orientation' and
• 'PixelList'

Pixel list is used to plot the original shape of the image. Quick and dirty way to see how good is the fit.

So here is the code:

1) Read the image, convert it to black and white (i.e. binary) and apply regionprops:

clear all
clc

A = im2bw(imread('DummyEllipse.jpg'));

s = regionprops(A,'Centroid','MajorAxisLength','MinorAxisLength','Orientation','PixelList')

2) Recover the parameters from the s structure. Note that I access the 2nd element from the structure since the 1st one corresponds to the edges of the plotting area (I don't know how to call it sorry).:

PixList = s(2).PixelList;

x = s(2).Centroid(1);
y = s(2).Centroid(2);
a = s(2).MajorAxisLength/2;
b = s(2).MinorAxisLength/2;
angle = s(2).Orientation;
steps = 50;

3) Provide these parameters to Amro code to compute the ellipse:

%#  @param x     X coordinate
%#  @param y     Y coordinate
%#  @param a     Semimajor axis
%#  @param b     Semiminor axis
%#  @param angle Angle of the ellipse (in degrees)

beta = angle * (pi / 180);
sinbeta = sin(beta);
cosbeta = cos(beta);

alpha = linspace(0, 360, steps)' .* (pi / 180);
sinalpha = sin(alpha);
cosalpha = cos(alpha);

X = x + (a * cosalpha * cosbeta - b * sinalpha * sinbeta);
Y = y + (a * cosalpha * sinbeta + b * sinalpha * cosbeta);

4) Plot the results:

figure

plot(PixList(:,1),rot90(PixList(:,2),2),'-r') % flip the y coordinates because indexing in Matlab starts form the top of the image.
hold on
plot(X,Y,'k','LineWidth',2)

hold off

which gives me this:

So it's not perfect but hopefully it will help you get started.