How can I detect the maximum-sized rectangle that

Question

I'm making an image processing project and I have stuck in one the project's steps. Here is the situation;

This is my mask:

and I want to detect the maximum-sized rectangle that can fit into this mask like this.

I'm using MATLAB for my project. Do you know any fast way to accomplish this aim. Any code sample, approach or technique would be great.

EDIT 1 : The two algorithms below are works with lot's of the cases. But both of them give wrong results in some difficult cases. I'am using both of them in my project.

Answer 1

Here is a correct answer. You must use dynamic programming! Other methods of direct calculation (like cutting 1 pixel from each edge) might produce sub-optimal results. My method guarantees that it selects the largest possible rectangle that fits in the mask. I assume that the mask has 1 big convex white blob of any shape with black background around it.

I wrote 2 methods. findRect() which finds the best possible square (starting on x,y with length l). The second method LargestInscribedImage() is an example of how to find any rectangle (of any aspect ratio). The trick is to resize the mask image, find a square and resize it back. In my example the method finds the larges rectangle that can be fit in the mask having the same aspect ration as the mask image. For example if the mask image is of size 100x200 pixels than the algorithm will find the largest rectangle having aspect ratio 1:2.

% ----------------------------------------------------------
function LargestInscribedImage()
% ----------------------------------------------------------
    close all
    im = double(imread('aa.bmp')); % Balck and white image of your mask
    im = im(:,:,1);  % If it is colored RGB take only one of the channels
    b = imresize(im,[size(im,1) size(im,1)]); Make the mask square by resizing it by its aspect ratio.
    SC = 1;  % Put 2..4 to scale down the image an speed up the algorithm

    [x1,y1,l1] = findRect(b,SC);          % Lunch the dyn prog algorithm
    [x2,y2,l2] = findRect(rot90(b),SC);   % rotate the image by 90deg and solve
    % Rotate back: x2,y2 according to rot90
    tmp = x2;
    x2 = size(im,1)/SC-y2-l2;
    y2 = tmp;

% Select the best solution of the above (for the original image and for the rotated by 90degrees
        if (l1>=l2)
           corn = sqCorn(x1,y1,l1);
        else
           corn = sqCorn(x2,y2,l2);
        end

        b = imresize(b,1/SC);
    figure;imshow(b>0); hold on;
    plot(corn(1,:),corn(2,:),'O')

    corn = corn*SC;
    corn(1,:) = corn(1,:)*size(im,2)/size(im,1);
    figure;imshow(im); hold on;
    plot(corn(1,:),corn(2,:),'O')
end

function corn = sqCorn(x,y,l)
     corn = [x,y;x,y+l;x+l,y;x+l,y+l]';
end
% ----------------------------------------------------------
function [x,y,l] = findRect(b,SC)
b = imresize(b,1/SC);
res = zeros(size(b,1),size(b,2),3);
% initialize first col
for i = 1:1:size(b,1)
    if (b(i,1) > 0)
       res(i,1,:) = [i,1,0];
    end
end
% initialize first row
for i = 1:1:size(b,2)
    if (b(1,i) > 0)
       res(1,i,:) = [1,i,0];
    end
end

% DynProg
for i = 2:1:size(b,1)
for j = 2:1:size(b,2)
    isWhite = b(i,j) > 0;
    if (~isWhite)
       res(i,j,:)=res(i-1,j-1,:); % copy
    else
        if (b(i-1,j-1)>0)  % continuous line
           lineBeg    = [res(i-1,j-1,1),res(i-1,j-1,2)]; 
           lineLenght = res(i-1,j-1,3); 
           if ((b(lineBeg(1),j)>0)&&(b(i,lineBeg(2))>0))  % if second diag is good
              res(i,j,:) = [lineBeg,lineLenght+1];
           else
              res(i,j,:)=res(i-1,j-1,:); % copy since line has ended
           end
        else
           res(i,j,:) = [i,j,0];         % Line start
        end
    end

end
end

% check last col
[maxValCol,WhereCol] = max(res(:,end,3));
% check last row
[maxValRow,WhereRow] = max(res(end,:,3));

% Find max
x= 0; y = 0; l = 0;
if (maxValCol>maxValRow)
   y = res(WhereCol,end,1);
   x = res(WhereCol,end,2);
   l = maxValCol;
else
   y = res(end,WhereRow,1);
   x = res(end,WhereRow,2);
   l = maxValRow;
end

    corn = [x,y;x,y+l;x+l,y;x+l,y+l]';
%    figure;imshow(b>0); hold on;
 %   plot(corn(1,:),corn(2,:),'O')
return;
end

Answer 2

The black boundaries in your image are curved and not closed. For example, in the top right corner, the black boundaries won't meet and form a closed contour. Therefore, a simple strategy in one of my comments will not work.

I am now providing you with a skeleton of a code which you can play with and add conditions as per your need. My idea is as follows:

To find left-side x-coordinate of the rectangle, first count the white pixels each column of the image contains:

%I assume that the image has already been converted to binary.

whitePixels=sum(img,1);

Then find the rate of change:

diffWhitePixels=diff(whitePixels);

If you see the bar plot of diffWhitePixels then you will observe various large entries (which indicate that the white region is still not in a straight line, and it is not a proper place to put the rectangles left vertical edge). Small entries (in your image, less than 5) indicate you can put the rectangle edge there.

You can do similar things to determine right, top and bottom edge positions of the rectangles.

Discussion:

First of all, the problem is ill-posed in my opinion. What do you mean by maximum-sized rectangle? Is it maximum area or length of side? In all possible cases, I don't think above method can get the correct answer. I can think of two or three cases right now where above method would fail, but it will at least give you the right answer on images similar to the given image, provided you adjust the values.

You can put some constraints once you know how your images are going to look. For example, if the black boundary curves inside, you can say that you don't want a column such as [0;0;...0;1;1;...0;0;...;0;1;1;...;1] i.e. zeros surrounded by ones. Another constraint could be how many black pixels do you want to allow? You can also crop the image till to remove extra black pixels. In your image, you can crop the image (programmatically) from the left and the bottom edge. Cropping an image is probably necessary, and definitely the better thing to do.

Answer 3

This approach starts with the entire image and shrinks each border in turn pixel-by-pixel until it finds an acceptable rectangle.

It takes ~0.02 seconds to run on the example image, so it's reasonably fast.

EDIT: I should clarify that this isn't meant to be a universal solution. This algorithm relies on the rectangle being centered and having roughly the same aspect ratio as the image itself. However, in the cases where it is appropriate, it is fast. @DanielHsH offered a solution which they claim works in all cases.

---

The code:

clear; clc;
tic;
%% // read image
imrgb= imread('box.png');
im = im2bw(rgb2gray(imrgb));    %// binarize image
im = 1-im;                      %// convert "empty" regions to 0 intensity
[rows,cols] = size(im);

%% // set up initial parameters
ULrow = 1;       %// upper-left row        (param #1)
ULcol = 1;       %// upper-left column     (param #2)
BRrow = rows;    %// bottom-right row      (param #3)
BRcol = cols;    %// bottom-right column   (param #4)

parameters = 1:4;   %// parameters left to be updated
pidx = 0;           %// index of parameter currently being updated

%% // shrink region until acceptable
while ~isempty(parameters); %// update until all parameters reach bounds

    %// 1. update parameter number
    pidx = pidx+1;
    pidx = mod( pidx-1, length(parameters) ) + 1;
    p = parameters(pidx);   %// current parameter number

    %// 2. update current parameter
    if p==1; ULrow = ULrow+1; end;
    if p==2; ULcol = ULcol+1; end;
    if p==3; BRrow = BRrow-1; end;
    if p==4; BRcol = BRcol-1; end;

    %// 3. grab newest part of region (row or column)
    if p==1; region = im(ULrow,ULcol:BRcol); end;
    if p==2; region = im(ULrow:BRrow,ULcol); end;
    if p==3; region = im(BRrow,ULcol:BRcol); end;
    if p==4; region = im(ULrow:BRrow,BRcol); end;

    %// 4. if the new region has only zeros, stop shrinking the current parameter
    if isempty(find(region,1))
        parameters(pidx) = [];
    end

end

toc;
params = [ULrow ULcol BRrow BRcol]
area = (BRrow-ULrow)*(BRcol-ULcol) 

---

The results for this image:

Elapsed time is 0.027032 seconds.

params =

    10    25   457   471


area =

      199362

Code to visualize results:

imrgb(params(1):params(3),params(2):params(4),1) = 0;
imrgb(params(1):params(3),params(2):params(4),2) = 255;
imrgb(params(1):params(3),params(2):params(4),3) = 255;
imshow(imrgb);

Another example image: