Can anyone suggest a fast method of estimating an ellipse bounded by a region?

I have a 2D binary image containing a single object, and I'm trying to determine the largest ellipse (by area) that fits entirely within that object. I've included an example below:

There are many techniques that will quickly fit an ellipse in a best-fit sense (like least-squares) to the object, but this does not ensure that the ellipse will be entirely within the bounds of the object - in fact, unless the object is a perfect ellipse, they usually guarantee that it won't.

The data set this is to work on requires hundreds, or potentially thousands, of "noisy" ellipses to be processed. Operators are likely to expect results in a minute or two.

I have tried using a cost function and search algorithms to refine a starting estimate, but this is too slow. I've also been using a home-made iterative technique. I fit an ellipse to the edge pixels, then discard those outside and re-fit using only the edges within the ellipse until there are no more edges outside. This, however, is quite slow and I'm not sure it is guaranteed to get the optimal answer.

Any help would be appreciated. I'm thinking that I may be re-inventing the wheel here.

Thanks in advance!

[Tim]

Thanks to those who answered. I ended up using a different strategy which worked surprisingly quickly given how inefficient it was.

 - Find all edge points
 - Estimate ellipse (can use Least Squares or one of many implementations)
 - Find focii of estimated ellipse
 - Calculate distances of each edge pixel from focii to determine
 - If every edge is outside of / on edge of ellipse (with some tolerance),
   stop
 - Remove edge furthest from focii and repeat

Answer 1

[Image Analyst]

You can use morphology if you have the Image Processing Toolbox. Here's how I'd proceed

1. Find the centroid using bwulterode(). You can't use regionprops like you'd think because that won't get you the centroid of the ellipse - it will get you the centroid of the entire tortuous shape. That's why we use morphology. Alternatively you can calculate the Euclidean Distance Transform with bwdist() and look for the max.
2. Now that you have the center of the ellipse you need to calculate and examine the Euclidean Distance Transform (EDT) to find out the minor axis length of the ellipse. This will be the max of your EDT.
3. Now call bwboundaries() to get all the perimeter coordinates. Run around the perimeter and find out which pixel has the same distance to the centroid as the minor axis length. Now you know which pixel is on the minor axis and so you can calculate the angle of the minor axis - call it alpha. OK, we're getting there...
4. Next you need to run around the perimeter again and find the major axis. Look at this web page: http://www.maa.org/joma/volume8/kalman/general.html. Look at equation 1 for a tilted ellipse. You know everything in that equation except b, the major axis length. You know x, y, a, and alpha. So now you have a quadratic equation for b so you can solve for b at every location.
5. There will be one b for which all perimeter points are outside the ellipse so you can check which b that is. Then you have your equation.

[Liang]

Great.

[Matt J]

I can't really see how bwulterode leads to the centroid. When I do

 bwulterode(Image)

on the OP's posted image, it produces an almost completely eroded result.

[Image Analyst]

It's supposed to. It eats away layer upon layer of binary blobs until only a single pixel remains at the center of each blob. If you don't see the dot, then maybe it's being lost upon subsampling for display and you'd have to check the image variable in the variable editor. But I haven't actually tried it with his image. I think that the "tail" should go away and not influence the centroid of the rest of the blob, but I'm not absolutely 100% sure.

[Matt J]

I think that the "tail" should go away and not influence the centroid of the rest of the blob, but I'm not absolutely 100% sure.

No, it doesn't. It does produce dots (several of them !?!) near the center of the large central blob, but also in the tail as well. That part makes sense, since wherever there is a narrow neck in a blob, the neck will lead to a local region max in the EDT.

Answer 2

[Matt J]

Here's another method based on FMINSEARCH. I'm sure it can be further speed-optimized by replacing the conv2 operations with fft-based convolutions.

 function maxellipse(Image)

    ImageArea=nnz(Image);

    [m,n]=size(Image);
    cxy=[m,n]/2+.5;

    [x,y]=ndgrid((1:m)-cxy(1),(1:n)-cxy(2));
     xygrid=[x(:),y(:)].';

     ImgAdjusted=single(Image*100-99);
     fun=@(p) objective(p,ImgAdjusted,xygrid, ImageArea); %tuning params

               %%Initial Guess
              S=regionprops(Image,'Orientation',...
                           'MajorAxisLength','MinorAxisLength' );

              p0(1)=S.MajorAxisLength/3;
              p0(2)=S.MinorAxisLength/3;
              p0(3)=S.Orientation;     

     p=fminsearch(fun,p0(:));

     [~,report]=fun(p);

     imagesc(report.ellipsemask+Image);

function  [cost,report]=objective(p,Image,xygrid, ImageArea)

   a=abs(p(1)); %ellipse length a 
   b=abs(p(2)); %ellipse length b
   theta=p(3);

   R=@(t)[cosd(t), -sind(t); sind(t) cosd(t) ];
    R=R(-theta);

   D=diag(1./([a,b]/2).^2);

   strel =  sum( ((R*D*R.')*xygrid).*xygrid, 1 ) <=1; %mask
    strel=reshape(single(strel), size(Image));

    map=conv2(Image,strel,'same');
    overlap=max(map(:));

    area=min(pi*a*b, ImageArea );

    cost=-(area+overlap);

    if nargout>1
        Map=conv2(double(Image>0),strel,'same')>=nnz(strel);
        idx=find(Map);
        Map(idx(2:end))=0;
        ImageEroded=Image & ~(conv2(double(Map),strel,'same')>0);
        BW=Image& ~ImageEroded;
        report=regionprops(BW,'Centroid','Orientation',...
                              'MajorAxisLength','MinorAxisLength' );
        report.ellipsemask=BW;
    end