How can I detect the maximum-sized rectangle that I can draw onto the mask?

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:

enter image description here

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

enter image description here

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.

This will be an inefficient solution. Suppose your image corners are represented by combinations of `xmin,ymin,xmax,ymax`. Now, start from `xmin` at row `1`, traverse right, if encountered white pixel, break. Do this for each row and remember the maximum value you had to traverse over all the rows. That number (coordinate) will be your `xminRect` for the rectangle. Similarly traverse up from `ymax`, traverse down from `ymin` and traverse left from `xmax` to get `ymaxRect,yminRect,xmaxRect` respectively. — Autonomous, May 09 '15 at 00:10
@ParagS.Chandakkar That's probably the best solution you can hope for, provided the boundaries don't get too crazy. Make that into an answer? — eigenchris, May 09 '15 at 00:16
@eigenchris **I did not make that as an answer because it won't work** (just keeping that comment for reference and for this conversation to make sense). I devised a new strategy which I am writing it down now. By using earlier strategy, on the first iteration of finding `xminRect`, the code will reach *almost* the right top corner. I didn't think of this simple bug when I wrote that comment. — Autonomous, May 09 '15 at 02:11
@ParagS.Chandakkar Good point. I didn't take that area in the corner into account either. This seems to be a surprisingly hard problem for such a simple task. — eigenchris, May 09 '15 at 03:17
In the general case, the answer is not unique if you don't tell how you define the size. It could be the area, the perimeter, the diagonal or whatever else. — , May 09 '15 at 07:41
There's a really neat animation at the end of this method... http://d3plus.org/blog/behind-the-scenes/2014/07/08/largest-rect/ — Mark Setchell, May 09 '15 at 09:11

eigenchris · Accepted Answer · 2015-05-09T14:07:36.077

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);

enter image description here

Another example image:

enter image description here

Nice job dude. Love it. +1. This in a sense is kind of a messed up active contours.... very cool. — rayryeng, May 09 '15 at 05:59
@rayryeng Thanks. I'm a little ashamed at how long it took me to come to this simple solution. I was trying all sorts of unnecessary stuff! — eigenchris, May 09 '15 at 06:19
Hahaha no worries. I've done that too.... but when it comes to finding a solution (such as your answer), always remember [Occam's Razor](http://en.wikipedia.org/wiki/Occam%27s_razor) - The fewer assumptions you make on the problem when solving it, the better... which is what I see with your solution. — rayryeng, May 09 '15 at 06:25
This algorithm is wrong. It might not produce correct results — DanielHsH, May 09 '15 at 08:02
I've tested it with different inputs. This approach has worked with my case awesomely. — Erdi İzgi, May 09 '15 at 11:32
@DanielHsH I should probably have added the caveat that this isn't meant to be a universal solution. For example, if the area is not centered this algorithm will fail. But for examples similar to the one given in the question, this algorithm seems to work well and relatively and fast, as requested. — eigenchris, May 09 '15 at 12:49
Yes but for shapes like hourglass (with extreme distortion) it will produce very poor results. This approach works only for very easy masks. It will find some rectangle but not the maximal one (like you asked). Moreover, you cannot control which aspect ratio of a rectangle you want. — DanielHsH, May 09 '15 at 15:05

DanielHsH · Answer 2 · 2015-05-09T08:41:04.753

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

I'm still trying to understand your code but it seems to work well. Definitely a nice solution if is guaranteed to work in all cases. Resizing the image to a square is a clever trick too! — eigenchris, May 09 '15 at 13:54
Please read in wikipedia about dynamic programing. The trick is to solve the problem to a smaller image first. Say if you have an image NxN and someone gives you a solution for sub image (N-1)x(N-1) You can easily add the last line and row to the solution, and get the full solution. So I start with a tiny image (1xN and Nx1) and solve in the initialization step. The next step in the for loop assumes existing solution for MxM image and extends it to (M+1)xM, then to Mx(M+1) and finally to MxM. — DanielHsH, May 09 '15 at 14:58
Please note: My algorithm will not work if your mask has holes (like white are of a donuts shape) The mask must have 1 big blob without holes. It basically checks all possible solutions but does it in N^2 steps instead of Naive N^8 steps implementation (checking all possible combinations of 4 corners). That is why this algorithm guarantees the largest square for a convex shape of the blob. It will work pretty well for non convex shapes as well, but can give sub-optimal results in extremely difficult cases. — DanielHsH, May 09 '15 at 15:02
eigenchris, BTW - I upvoted your solution because it is very fast, intuitive and quite good for simple masks. — DanielHsH, May 09 '15 at 15:33
@DanielHsH As you said, "1 big convex white blob of any shape", so you need the white blob to be convex? If so, then the image in the question itself is not a convex white blob. — Autonomous, May 09 '15 at 18:47
For convext it is guaranteed to give a 100% correct result. For a bit non convex (like in the example) it will work fine, but for non convex with holes it will give errors. My algorithm checks that the diagonals of the rectangle are completely white (inside mask). One can easily change it to test if the sides are white instead of diagonals (or in addition). This will ensure also 100% correctness for non convex shapes but without holes. With holes it is much more difficult (can be solved in N^3 complexity instead of N^2) — DanielHsH, May 10 '15 at 06:37

score 1 · Answer 3 · edited Jun 20 '20 at 09:12

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.

How can I detect the maximum-sized rectangle that I can draw onto the mask?

3 Answers3

Discussion:

Linked

Related