Maximum inscribed quadrilateral inside the convex hull
Ahmed Hegazy
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
ImageAnalyst
wrote:
-------------------------------------------------------------------
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);
Bruno Luong
> 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
Bruno Luong
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
Bruno Luong
> -------------------------------------------------------------------
> 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
Sorry 2 milisecond
Bruno
ImageAnalyst
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.
Bruno Luong
> 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
Ahmed Hegazy
> ImageAnalyst <imagea...@mailinator.com> wrote in message <c24f9022-8bad-4c0d...@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 .
Image Analyst
Bruno Luong
> 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
Ahmed Hegazy
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?
ImageAnalyst
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.
Bruno Luong
> 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
Adam Geboff
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.l...@fogale.findmycountry> wrote in message <hef59m$d08$1...@fred.mathworks.com>...
Bruno Luong
> 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
ImageAnalyst
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.