Thread Subject: Maximum inscribed quadrilateral inside the convex hull

Hi all ,
please i want to ask about the way to get the area of the maximum inscribed quadrilateral inside the convex hull .

Thnx in advance,
Ahmed

On Nov 23, 5:09 pm, "Ahmed Hegazy" <lord_of_the_rings_3...@yahoo.com>
wrote:
> Hi all ,
> please i want to ask about the way to get the area of the maximum inscribed quadrilateral inside the convex hull .
>
> Thnx in advance,
> Ahmed

-------------------------------------------------------------------
There may be more than one with the same area, especially if the shape
is symmetric. I don't know of any super efficient way. I know that
you can just use (clever) brute force to check all possibilities, like
this (take about 11 seconds for this 18 vertex point convex hull: (Be
sure to fix any lines broken into two by the newsreader).

clc;
close all;
clear all;
workspace;

% Run MATLAB Example 1 from the documentation
xx = -1:.05:1; yy = abs(sqrt(xx));
[x,y] = pol2cart(xx,yy);
k = convhull(x,y);
plot(x(k),y(k),'r-',x,y,'b+')

% Now let's take over coding.....
set(gcf, 'Position', get(0, 'ScreenSize')); % Maximize figure.
drawnow;
uiwait(msgbox('The convex hull is in red. Click OK to find the
largest quad (in about 11 seconds).'));
drawnow; % Force message box to go away.

% Now find the max area
tic;
cvX = x(k);
cvY = y(k);
numberOfPoints = length(cvX)
maxArea = 0;
% Be sure to skip any overlapping vertices to make it more efficient!
for k1 = 1 : numberOfPoints
  for k2 = 1 : numberOfPoints
    if k2 == k1, continue, end;
    for k3 = 1 : numberOfPoints
      if k3 == k2 || k3 == k1, continue, end;
      for k4 = 1 : numberOfPoints
        if k4 == k3 || k4 == k2 || k4 == k1, continue, end;
verticesX = [cvX(k1) cvX(k2) cvX(k3) cvX(k4)];
verticesY = [cvY(k1) cvY(k2) cvY(k3) cvY(k4)];
thisArea = polyarea(verticesX, verticesY);
if thisArea > maxArea
maxArea = thisArea;
indexesAtmax = [k1 k2 k3 k4];
end
end
end
  end
end

% Plot the quad inside the convex hull.
boxX = [cvX(indexesAtmax(1)), cvX(indexesAtmax(2)), cvX(indexesAtmax
(3)), cvX(indexesAtmax(4)), cvX(indexesAtmax(1))]
boxY = [cvY(indexesAtmax(1)), cvY(indexesAtmax(2)), cvY(indexesAtmax
(3)), cvY(indexesAtmax(4)), cvY(indexesAtmax(1))]
hold on;
plot(boxX, boxY, 'g');
message = sprintf('The largest quad is shown in green.\nThe max area =
%f, and occurs with indexes [%d %d %d %d]', ...
maxArea, indexesAtmax(1), indexesAtmax(2), indexesAtmax(3),
indexesAtmax(4));

toc;
msgbox(message);

"Ahmed Hegazy" <lord_of_the_rings_3000@yahoo.com> wrote in message <hef16f$rnt$1@fred.mathworks.com>...
> Hi all ,
> please i want to ask about the way to get the area of the maximum inscribed quadrilateral inside the convex hull .
>

I have wrote the code for largest triangle here. I believe you can modify the code for quadrilateral by adding the two triangles from other sides (Tmin + Tmax) to form a quadrilateral.

http://www.mathworks.com/matlabcentral/newsreader/view_thread/250533

Bruno

Here we go, I make a change for quadrilateral:

function result = hullmaxquadri(caseid, graphic, verbose)
% function result = hullmaxquadri(X)
%
% Find the largest quadrilateral inside the convex hull defined by X
% X is (2 x N) or (N x 2) array, coordinates of N-points in R^2,
% not necessary sorted.
%
% hullmaxquadri(caseid) to run test case (caseif 1 to 3)
%

try % turn off something in Bruno's computer
    bsxops(false);
end

if nargin<1 || isempty(caseid)
    caseid = 1;
end

if nargin<2 || isempty(graphic)
    graphic = 1;
end

if nargin<3 || isempty(verbose)
    verbose = 1;
end

% Data
if isscalar(caseid)
    switch caseid
        case 1,
            m = 1e5;
            x=randn(1,m)*3;
            y=randn(1,m);
        case 2,
            % LONG, 12 seconds
            theta=linspace(0,2*pi,1e3);
            theta(end)=[];
            x=cos(theta);
            y=sin(theta);
        case 3,
            % Medium, 250 ms
            theta=linspace(0,2*pi,2e2);
            theta(end)=[];
            x=cos(theta);
            y=2*sin(theta);
        otherwise
            error('unknown caseid');
    end
else
    X = caseid;
    if size(X,2)==2
        X = X.';
    end
    x = X(1,:);
    y = X(2,:);
end
% Make sure points are sorted and convex
k=convhull(x,y);
X=[x(k); y(k)];
clear x y;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Engine
% size
n=size(X,2);

% All two-vertices combo
[J I] = find(tril(ones(n),-1));

% Normal vector to all couple of points
normal = [X(2,J)-X(2,I); ...
          X(1,I)-X(1,J)].';

% The following block is used later to discard points too close to each
% other. The idea is take the triangle with the largest width enclosed
% inside the convex hull, then take its height as threshold. All the pairs
% of points having distance smaller than the threshold would not be good
% candidate since the height must be larger than the largest side of the
% convex hull.

% Find the distance threshold under which a pair of points are surely
% not solution (two corners of the triangle that maximizes the area)
lgt2 = sum(normal.^2,2);
[maxlgt imax] = max(lgt2); % locate where is the largest side?
% DX = bsxfun(@minus, X, X(:,I(imax)));
DX = X - repmat(X(:,I(imax)), [1 n]);
maxT = max(normal(imax,:)*DX); % how it expands on the orthogonal?
minT = min(normal(imax,:)*DX); % how it expands on the orthogonal?
threshold = 0.5*abs(maxT-minT) / sqrt(maxlgt);
tooclose = lgt2 < (threshold^2); % boolean flag, false for pair that
                                 % will be discarded
clear lgt2 % clean up

% Split with the same j
isplit = mat2cell(1:size(I,1),1,n-1:-1:1);
AREAMAX = -Inf;
innerlooptime=zeros(1,length(isplit));
% array used to move to the next point by indexing it
nextwrap = [2:n 1];
for k=1:length(isplit)
    tic
    idx = isplit{k};
    DX = bsxfun(@minus, X, X(:,I(idx(1))));
    %DX = X - repmat(X(:,I(1)), [1 n]);
    firstj = true;
    for j=idx
        if tooclose(j) % skip it
            continue
        end
        nj = normal(j,:);
        if firstj
            p = nj*DX;
            [Tmin imin] = min(p);
            [Tmax imax] = max(p);
            firstj=false;
        else
            % Compute area
            Tmin = nj*DX(:,imin);
            % This loop is O(1) complexity in average.
            while true
                next = nextwrap(imin);
                Tnext = nj*DX(:,next);
                if Tnext < Tmin
                    Tmin = Tnext;
                    imin = next;
                else
                    break
                end
            end
            Tmax = nj*DX(:,imax);
            while true
                next = nextwrap(imax);
                Tnext = nj*DX(:,next);
                if Tnext > Tmax
                    Tmax = Tnext;
                    imax = next;
                else
                    break
                end
            end
        end % if ~firstj

        Qmax = Tmax-Tmin;
        % Keep track the largest
        if abs(Qmax)>AREAMAX
           AREAMAX = abs(Qmax);
           IdxQuadMax = [I(j) imax J(j) imin];
        end
    end % for loop on j
    innerlooptime(k)=toc;
end % for loop on split

AREAMAX = AREAMAX/2;
CoordinatesQuad = X(:,IdxQuadMax);
clear normal I J

Perimeter = sum(sqrt(sum(diff(CoordinatesQuad(:,[1:end 1]),1,2).^2,1)));

result = struct('Hull',X,...
                'AREAMAX',AREAMAX, ...
                'Perimeter', Perimeter, ...
                'IdxQuadMax', IdxQuadMax, ...
                'CoordinatesQuad', CoordinatesQuad, ...
                'innerlooptime', innerlooptime);

%%% Verbose
if verbose
    fprintf('Dimensions = %d\n', n)
    fprintf('AREAMAX = %g\n', AREAMAX)
    fprintf('IdxQuadMax @ [%d %d %d %d]\n', IdxQuadMax)
    fprintf('XCoords = [%g %g %g %g]\n', CoordinatesQuad(1,:))
    fprintf('YCoords = [%g %g %g %g]\n', CoordinatesQuad(2,:))
end

%%%%%%%% Graphic
if graphic
    fig=figure(1);
    clf(fig);
    ax=subplot(1,2,1,'Parent',fig);
    plot(ax,n-1:-1:1,innerlooptime);
    xlabel(ax,'size');
    ylabel(ax,'Inner loop time [s]');
    ax=subplot(1,2,2,'Parent',fig);
    plot(ax,X(1,[1:end 1]),X(2,[1:end 1]),'k-',...
        X(1,IdxQuadMax([1:end 1])),X(2,IdxQuadMax([1:end 1])),'ro-');
    axis(ax,'equal')
end

end

% Bruno

ImageAnalyst <imageanalyst@mailinator.com> wrote in message <48a55ef2-d5dd-401f-a3d2-765184198699@m35g2000vbi.googlegroups.com>...

> -------------------------------------------------------------------
> There may be more than one with the same area, especially if the shape
> is symmetric. I don't know of any super efficient way. I know that
> you can just use (clever) brute force to check all possibilities, like
> this (take about 11 seconds for this 18 vertex point convex hull: (Be
> sure to fix any lines broken into two by the newsreader).

brute force method is indeed long... my code takes about 0.11 milisecond for the same.

Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <hef5mr$85i$1@fred.mathworks.com>...
> ImageAnalyst <imageanalyst@mailinator.com> wrote in message <48a55ef2-d5dd-401f-a3d2-765184198699@m35g2000vbi.googlegroups.com>...
>
> > -------------------------------------------------------------------
> > There may be more than one with the same area, especially if the shape
> > is symmetric. I don't know of any super efficient way. I know that
> > you can just use (clever) brute force to check all possibilities, like
> > this (take about 11 seconds for this 18 vertex point convex hull: (Be
> > sure to fix any lines broken into two by the newsreader).
>
> brute force method is indeed long... my code takes about 0.11 milisecond for the same.

Sorry 2 milisecond

Bruno

I didn't dig in to what you did but it must be very clever. I ran it
and it did indeed give the same result with the 18 vertex MATLAB demo
convex hull as the brute force method. Impressive. You also have an
impressively fast computer - like 80 times as fast as mine!
My brute force method = 11 seconds
Bruno's method on my computer = 0.167 seconds = 167 ms
Bruno's method on Bruno's computer = 0.002 seconds = 2 ms (Wow!)

I'd be interested in hearing from Ahmed what use he will make of the
quad - it's a need I've never heard before, and I was just wondering.

ImageAnalyst <imageanalyst@mailinator.com> wrote in message <c24f9022-8bad-4c0d-8597-8210d453688f@r31g2000vbi.googlegroups.com>...
> I didn't dig in to what you did but it must be very clever. I ran it
> and it did indeed give the same result with the 18 vertex MATLAB demo
> convex hull as the brute force method. Impressive. You also have an
> impressively fast computer - like 80 times as fast as mine!
> My brute force method = 11 seconds
> Bruno's method on my computer = 0.167 seconds = 167 ms
> Bruno's method on Bruno's computer = 0.002 seconds = 2 ms (Wow!)

ImageAnalyst, note that I discard the preprocessing step when calling convexhull before timing. My home computer is not fast one, a 2 year old Core 2 duo 2GHz laptop.

The "trick" is split the quadrilateral in two triangles and parametrize it by their common edge which comprises two vertexes on the hull. The two other opposite vertexes maximize the respective triangle areas and can be compute by constant CPU time if we track them carefully. So the algorithm essentially reduces the search to a common edge (N^2 complexity where N is the number of hull vertexes) which is done by listing all possible combination.

The awfully difficult case is when the hull has much more points - for example when running the test case hullmaxquadri(2). The hull comprises 1000 points distributed on a circle. On my computer it takes 5 seconds to find the maximum quadrilateral (which converge toward a square).

Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <hefvqj$qvn$1@fred.mathworks.com>...
> ImageAnalyst <imageanalyst@mailinator.com> wrote in message <c24f9022-8bad-4c0d-8597-8210d453688f@r31g2000vbi.googlegroups.com>...
> > I didn't dig in to what you did but it must be very clever. I ran it
> > and it did indeed give the same result with the 18 vertex MATLAB demo
> > convex hull as the brute force method. Impressive. You also have an
> > impressively fast computer - like 80 times as fast as mine!
> > My brute force method = 11 seconds
> > Bruno's method on my computer = 0.167 seconds = 167 ms
> > Bruno's method on Bruno's computer = 0.002 seconds = 2 ms (Wow!)


Thnx ImageAnalyst and bruno for your help realy i appreciate it alot .
i used Bruno's method....it is wonderful and really fast i didnot count the milliseconds :D but the problem was that i want to use it as a feature to detect the original shape of a given convex hull and i want to run this function for 420 training pattern ... the brute force approach was really slow and i waited alot and i couldnot see the results :D so i used Bruno's method. Thnx again for your help and sorry for the late reply .

I'm not sure how finding the largest inscribed quadrilateral helps you "detect the original shape of a given convex hull" - you said you already had the original convex hull to begin with.

"Image Analyst" <imageanalyst@mailinator.com> wrote in message <heum58$rrm$1@fred.mathworks.com>...
> I'm not sure how finding the largest inscribed quadrilateral helps you "detect the original shape of a given convex hull" - you said you already had the original convex hull to begin with.

It is possibly way to compress any convex shape into to 4 points for the purpose of training something, like making a machine to recognize a bad fortune cookie?

Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <heuper$hk5$1@fred.mathworks.com>...
> "Image Analyst" <imageanalyst@mailinator.com> wrote in message <heum58$rrm$1@fred.mathworks.com>...
> > I'm not sure how finding the largest inscribed quadrilateral helps you "detect the original shape of a given convex hull" - you said you already had the original convex hull to begin with.
>
> It is possibly way to compress any convex shape into to 4 points for the purpose of training something, like making a machine to recognize a bad fortune cookie?
>
> Bruno

I am working in a pattern recognition project , so after i got the convex hull of the image i want to check if it is a circle or a triangle or a rectangle , but unfortuantly the results for the triangle and rectangles was so close when testing it using some test cases , Do you know any features that can help in this classification :D:D:D:D?

If you compute the "circularity" you might be able to detect the
difference, alhtough I can still think of mathematically possible quad
that are pretty darn close to a triangle. But in practice, with your
actual images, it may be that they're all different enough to work.
You'd just have to measure a bunch and see if it works. The
circularity is gotten by assuming it's a circle
A = pi * Diameter^2 / 4
Now MATLAB returns the perimeter and area, so for a circle
P = pi * D
A = pi * (P/pi)^2 / 4 = P^2 / (4*pi)
If you divide these
1 = P^2/(4 * pi * A) = circularity
For a circle, circularity = 1 but the more tortuous the perimeter
gets, P increases while A doesn't (or increases less) and so the
circularity rises above 1. It could be that you calculate the
circularity and find, through checking test cases, that
for more or less circular objects, circ = 1 - 3
for rectangle-like objects, circ = 3 - 8
for triangular objects, circ > 8
or something like that. You'd have to play with the numbers that
define the ranges. If all you need is to detect whether something is
closer to a circle, rectangle, or a triangle, then you should try this
first.

"Ahmed Hegazy" <lord_of_the_rings_3000@yahoo.com> wrote in message <heuv0r$oar$1@fred.mathworks.com>...

> I am working in a pattern recognition project , so after i got the convex hull of the image i want to check if it is a circle or a triangle or a rectangle , but unfortuantly the results for the triangle and rectangles was so close when testing it using some test cases , Do you know any features that can help in this classification :D:D:D:D?

I would think the angles of the quadrilateral (some sort of discrete curvature) is a good test to distinguish triangle to rectangle. You can train a learning machine (ANN/SVM) I'm pretty sure it will perform well.

Bruno

How would one modify your code to find the largest rectangle inscribed within a drawn loop? The difference between this and the convex hull is that the convex hull finds the outermost points to create its loop while I need to consider the innermost points.

I.e. one draws a shape with the 'getPosition(imfreehand)' function and a rectangle completely contained within those points with the maximum area is outputted.

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <hef59m$d08$1@fred.mathworks.com>...
> Here we go, I make a change for quadrilateral:
>
> function result = hullmaxquadri(caseid, graphic, verbose)
> % function result = hullmaxquadri(X)
> %
> % Find the largest quadrilateral inside the convex hull defined by X
> % X is (2 x N) or (N x 2) array, coordinates of N-points in R^2,
> % not necessary sorted.
> %
> % hullmaxquadri(caseid) to run test case (caseif 1 to 3)
> %
>
> try % turn off something in Bruno's computer
> bsxops(false);
> end
>
> if nargin<1 || isempty(caseid)
> caseid = 1;
> end
>
> if nargin<2 || isempty(graphic)
> graphic = 1;
> end
>
> if nargin<3 || isempty(verbose)
> verbose = 1;
> end
>
> % Data
> if isscalar(caseid)
> switch caseid
> case 1,
> m = 1e5;
> x=randn(1,m)*3;
> y=randn(1,m);
> case 2,
> % LONG, 12 seconds
> theta=linspace(0,2*pi,1e3);
> theta(end)=[];
> x=cos(theta);
> y=sin(theta);
> case 3,
> % Medium, 250 ms
> theta=linspace(0,2*pi,2e2);
> theta(end)=[];
> x=cos(theta);
> y=2*sin(theta);
> otherwise
> error('unknown caseid');
> end
> else
> X = caseid;
> if size(X,2)==2
> X = X.';
> end
> x = X(1,:);
> y = X(2,:);
> end
> % Make sure points are sorted and convex
> k=convhull(x,y);
> X=[x(k); y(k)];
> clear x y;
>
> %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
> % Engine
> % size
> n=size(X,2);
>
> % All two-vertices combo
> [J I] = find(tril(ones(n),-1));
>
> % Normal vector to all couple of points
> normal = [X(2,J)-X(2,I); ...
> X(1,I)-X(1,J)].';
>
> % The following block is used later to discard points too close to each
> % other. The idea is take the triangle with the largest width enclosed
> % inside the convex hull, then take its height as threshold. All the pairs
> % of points having distance smaller than the threshold would not be good
> % candidate since the height must be larger than the largest side of the
> % convex hull.
>
> % Find the distance threshold under which a pair of points are surely
> % not solution (two corners of the triangle that maximizes the area)
> lgt2 = sum(normal.^2,2);
> [maxlgt imax] = max(lgt2); % locate where is the largest side?
> % DX = bsxfun(@minus, X, X(:,I(imax)));
> DX = X - repmat(X(:,I(imax)), [1 n]);
> maxT = max(normal(imax,:)*DX); % how it expands on the orthogonal?
> minT = min(normal(imax,:)*DX); % how it expands on the orthogonal?
> threshold = 0.5*abs(maxT-minT) / sqrt(maxlgt);
> tooclose = lgt2 < (threshold^2); % boolean flag, false for pair that
> % will be discarded
> clear lgt2 % clean up
>
> % Split with the same j
> isplit = mat2cell(1:size(I,1),1,n-1:-1:1);
> AREAMAX = -Inf;
> innerlooptime=zeros(1,length(isplit));
> % array used to move to the next point by indexing it
> nextwrap = [2:n 1];
> for k=1:length(isplit)
> tic
> idx = isplit{k};
> DX = bsxfun(@minus, X, X(:,I(idx(1))));
> %DX = X - repmat(X(:,I(1)), [1 n]);
> firstj = true;
> for j=idx
> if tooclose(j) % skip it
> continue
> end
> nj = normal(j,:);
> if firstj
> p = nj*DX;
> [Tmin imin] = min(p);
> [Tmax imax] = max(p);
> firstj=false;
> else
> % Compute area
> Tmin = nj*DX(:,imin);
> % This loop is O(1) complexity in average.
> while true
> next = nextwrap(imin);
> Tnext = nj*DX(:,next);
> if Tnext < Tmin
> Tmin = Tnext;
> imin = next;
> else
> break
> end
> end
> Tmax = nj*DX(:,imax);
> while true
> next = nextwrap(imax);
> Tnext = nj*DX(:,next);
> if Tnext > Tmax
> Tmax = Tnext;
> imax = next;
> else
> break
> end
> end
> end % if ~firstj
>
> Qmax = Tmax-Tmin;
> % Keep track the largest
> if abs(Qmax)>AREAMAX
> AREAMAX = abs(Qmax);
> IdxQuadMax = [I(j) imax J(j) imin];
> end
> end % for loop on j
> innerlooptime(k)=toc;
> end % for loop on split
>
> AREAMAX = AREAMAX/2;
> CoordinatesQuad = X(:,IdxQuadMax);
> clear normal I J
>
> Perimeter = sum(sqrt(sum(diff(CoordinatesQuad(:,[1:end 1]),1,2).^2,1)));
>
> result = struct('Hull',X,...
> 'AREAMAX',AREAMAX, ...
> 'Perimeter', Perimeter, ...
> 'IdxQuadMax', IdxQuadMax, ...
> 'CoordinatesQuad', CoordinatesQuad, ...
> 'innerlooptime', innerlooptime);
>
> %%% Verbose
> if verbose
> fprintf('Dimensions = %d\n', n)
> fprintf('AREAMAX = %g\n', AREAMAX)
> fprintf('IdxQuadMax @ [%d %d %d %d]\n', IdxQuadMax)
> fprintf('XCoords = [%g %g %g %g]\n', CoordinatesQuad(1,:))
> fprintf('YCoords = [%g %g %g %g]\n', CoordinatesQuad(2,:))
> end
>
> %%%%%%%% Graphic
> if graphic
> fig=figure(1);
> clf(fig);
> ax=subplot(1,2,1,'Parent',fig);
> plot(ax,n-1:-1:1,innerlooptime);
> xlabel(ax,'size');
> ylabel(ax,'Inner loop time [s]');
> ax=subplot(1,2,2,'Parent',fig);
> plot(ax,X(1,[1:end 1]),X(2,[1:end 1]),'k-',...
> X(1,IdxQuadMax([1:end 1])),X(2,IdxQuadMax([1:end 1])),'ro-');
> axis(ax,'equal')
> end
>
> end
>
> % Bruno

"Adam Geboff" wrote in message <ii5bqo$4rt$1@fred.mathworks.com>...
> How would one modify your code to find the largest rectangle inscribed within a drawn loop? The difference between this and the convex hull is that the convex hull finds the outermost points to create its loop while I need to consider the innermost points.
>
> I.e. one draws a shape with the 'getPosition(imfreehand)' function and a rectangle completely contained within those points with the maximum area is outputted.

It is useless attempt to modify my code, because it's specifically designed for the quadrilateral in a convex hull and I can't see how it can adapt to your problem. You better start from scratch.

Bruno

Perhaps it might work if you took the Euclidean Distance Transform
with bwdist() and then locate the highest point. This will be the
interior point farthest away from any boundary point. Then assume a
rectangle is going to be centered at that point with a half-diagonal
that goes from that point out to the closest boundary point. Once you
know the boundary point of that, it would be easy to find the other
point diagonally opposite from it and the two corners at
perpendiculars to that largest diagonal. This would give you the
inscribed rectangle with the largest diagonal. Now the rectangle with
the largest possible diagonal doesn't always have the max area but in
this case I believe it will because of the way the Euclidean Distance
Transform works - it finds the largest *smallest* diagonal that will
fit. Off the top of my head I can't think of any counterexamples
where it wouldn't work.