Matlab最小面积包围四边形

对于存在透视变换的物体，提取时最小面积包围矩形不能满足要求，google到一个求最小面积包围四边形的算法，虽然速度较慢。以提取书本为例，实验结果和代码如下。

booktest.m

I = imread('book1.jpg');
figure(1);
imshow(I);
Ir = I(:,:,1);
Ismall = imresize(Ir,0.25);
bw =~imbinarize(Ismall);
imshow(bw);
bw2 = bwareaopen(bw,1000);
imshow(bw2);

stats = regionprops(bw2,'BoundingBox','ConvexHull');
hold on;
for i=1:numel(stats)
    %凸包
    h1=impoly(gca,stats(i).ConvexHull);
    api = iptgetapi(h1);
    api.setColor('red');
    %包围矩形
    h2=imrect(gca,stats(i).BoundingBox);
    api = iptgetapi(h2);
    api.setColor('blue');
    %包围盒（最小面积包围矩形）
    point = stats(i).ConvexHull;
    c = minBoundingBox(point');
    plot(c(1,[1:end 1]),c(2,[1:end 1]),'g')
    %最小面积包围四边形
    [qx,qy] = minboundquad(point(:,1),point(:,2));
    h3=impoly(gca,[qx',qy']);
    api = iptgetapi(h3);
    api.setColor('yellow');
end
hold off;

%整理四个点位置
%记最左点为1，同最左的话，靠下为1
%       2
%1
%            3
%    4

quadVetices = [qx(1:4)',qy(1:4)'];
quadVetices  = quadVetices .* 4;
[quadVeticesSort,index] = sortrows(quadVetices);
if quadVeticesSort(1,1) == quadVeticesSort(2,1)
    if quadVeticesSort(1,2) < quadVeticesSort(2,2)
        k = index(1);
    else
        k = index(2);
    end
else
    k = index(1);
end
i = mod((k-1):(k+2),4)+1;

rRect = quadVetices(i,:);
%计算边长
len = zeros(4,1);
for i = 1:3
    len(i) = sqrt(power(rRect(i,1)-rRect(i+1,1),2)+power(rRect(i,2)-rRect(i+1,2),2));
end
len(4) = sqrt(power(rRect(4,1)-rRect(1,1),2)+power(rRect(4,2)-rRect(1,2),2));
dsth = (len(1)+len(3))/2;
dstw = (len(2)+len(4))/2;
dstRect  = [0,dsth;0,0;dstw,0;dstw,dsth;];
%纠正透视变换
tform = estimateGeometricTransform(rRect, dstRect,...
    'projective');
imgBp = imwarp(I, tform,...
    'OutputView',imref2d([floor(dsth),floor(dstw)]));

imshow(imgBp);

2D minimal bounding box-minBoundingBox()
http://cn.mathworks.com/matlabcentral/fileexchange/31126-2d-minimal-bounding-box?focused=5188215&tab=function
A suite of minimal bounding objects-minboundquad()
http://cn.mathworks.com/matlabcentral/fileexchange/34767-a-suite-of-minimal-bounding-objects?focused=3820662&tab=function
这个代码包中还包含了最小面积包围三角形，包围圆，椭圆等多种形状。

现在总是报警告，暂时没时间细看了

In minboundquad (line 181)
In booktest (line 27)
警告: 矩阵为奇异工作精度。

minboundquad.m
我稍微改了下，
在第125行增加nedges = size(edges,1);

function [qx,qy,quadarea] = minboundquad(x,y)
% minboundquad: Compute the minimum area bounding quadrilateral of points in the plane
% usage: [qx,qy] = minboundquad(x,y)
%
% arguments: (input)
%  x,y - vectors of points, describing points in the plane as
%        (x,y) pairs. x and y must be the same size.
%
% arguments: (output)
%  qx,qy - 5x1 vectors of points that define the minimum
%        area bounding quadrilateral.
%
% WARNING: Most 2-d convex hulls are fairly small.
% However, this code will be O(N^4), where N is the
% number of distinct edges in the convex hull of
% your data. So finding the bounding quadrilateral
% of 1000 points around the perimeter of a circle 
% will take much time.
%
% Example usage:
%  x = randn(50,1);
%  y = randn(50,1);
%  [qx,qy] = minboundquad(x,y);
%  plot(x,y,'ro',qx,qy,'b-')
%
%
% See also: minboundcircle, minboundrect, minboundtri
%
%
% Author: John D'Errico
% E-mail: woodchips@rochester.rr.com
% Release: 1.0
% Release date: 3/14/07

% preprocess data
x=x(:);
y=y(:);

% not many error checks to worry about
n = length(x);
if n~=length(y)
  error('MINBOUNDQUAD:size','x and y must be the same sizes')
end

% start out with the convex hull of the points to
% reduce the problem dramatically. Note that any
% points in the interior of the convex hull are
% never needed, so we drop them.
xy = [x,y];
if n>3
  edges = convhull(x,y);
  % convhull returns a list of points around
  % the perimeter. I prefer the convhulln form,
  % where I have an explicit list of edges.
  % Make it so.
  edges = [edges(1:(end-1)),edges(2:end)];
elseif n==3
  % it is a triangle. I don't care how we
  % traverse it. Replicate a vertex into
  % a quadrilateral.
  qx = xy([1 2 3 3 1],1);
  qy = xy([1 2 3 3 1],2);
  quadarea = polyarea(qx,qy);
  return
elseif n==2
  % a single edge
  qx = xy([1 2 2 2 1],1);
  qy = xy([1 2 2 2 1],2);
  % no area to be found
  quadarea = 0;
  return
elseif n==1
  % a single point
  qx = xy([1 1 1 1 1],1);
  qy = xy([1 1 1 1 1],2);
  quadarea = 0;
  return
else
  % empty begets empty
  qx = [];
  qy = [];
  quadarea = 0;
  return
end
nedges = size(edges,1);

% now we must find the bounding quadrilateral of those
% that remain.

% special case small numbers of points. If we trip any
% of these cases, then we are done, so return.
if nedges == 3
  qx = xy([1 2 3 1],1);
  qy = xy([1 2 3 1],2);
  return
elseif nedges == 4
  qx = xy([1 2 3 4 1],1);
  qy = xy([1 2 3 4 1],2);
  return
end
% more than 4 points.

% get the angle of each edge of the hull polygon
edgeangles = atan2(xy(edges(:,2),2) - xy(edges(:,1),2), ...
  xy(edges(:,2),1) - xy(edges(:,1),1));

% (These next two steps are probably superfluous.)
% work with all positive angles
k = edgeangles < 0;
edgeangles(k) = edgeangles(k) + 2*pi;
% sort the edges into increasing order of angle
[edgeangles,tags] = sort(edgeangles);
edges = edges(tags,:);

% Look for consecutive edges that have the same
% angles. This test will generally only trip if
% the data set has multiple collinear points
% around the perimeter.
angletol = eps*100;
k = diff(edgeangles) < angletol;
edges(k,:) = [];
edgeangles(k) = [];

% there are nchoosek(nedges,4) sets of edges to
% worry about
nedges = size(edges,1);
edgelist = nchoosek(1:nedges,4);

% The edges are now sorted in counter-clockwise
% order around the convex hull. We can toss out any
% combination of edges where the last edge angle
% minus the first is less than 180 degrees
% (i.e., pi radians.)
k = (edgeangles(edgelist(:,4)) - edgeangles(edgelist(:,1)) <= pi);
edgelist(k,:) = [];

% how many edges remain that can form a valid
% quadrilateral?
nquads = size(edgelist,1);

% test each set of 4 edges
quadarea = inf;
qxi = zeros(1,5);
qyi = zeros(1,5);
qx = qxi;
qy = qyi;
for i = 1:nquads
  % find the intersections of each consecutive
  % pair of edges.
  edgeind = edgelist(i,:);
  edgesi = edges(edgeind([1 2 3 4 1]),:);

  if any(diff(edgeangles(edgeind)) > pi)
    % if one of the consecutive angles is too
    % large, then this set of edges will be a
    % failed quadrilateral.
    continue
  end

  for j = 1:4
    % Does this pair of edges share a node from
    % the convex hull?
    jplus1 = j + 1;
    shared = intersect(edgesi(j,:),edgesi(jplus1,:));
    if ~isempty(shared)
      % there was a shared node between these edges
      qxi(j) = xy(shared,1);
      qyi(j) = xy(shared,2);
    else
      % no shared node, so we must find the
      % intersection of the edges by extrapolation
      % of the lines that contain these edges to
      % see where they intersect.
      A = xy(edgesi(j,1),:);
      B = xy(edgesi(j,2),:);
      C = xy(edgesi(jplus1,1),:);
      D = xy(edgesi(jplus1,2),:);

      % solve for the line parameters that correspond
      % to the intersection point
      ts = [(A-B)',(D-C)']\(A-C)';
      % recover the intersection point
      Q = A + (B-A)*ts(1);

      qxi(j) = Q(1);
      qyi(j) = Q(2);
    end
  end
  % wrap the polygon
  qxi(5) = qxi(1);
  qyi(5) = qyi(1);

  % compute the area. Simplest is to use polyarea.
  % I might want to test to see if it is faster
  % to compute the area using other methods though.
  A_i = polyarea(qxi,qyi);

  if (A_i < quadarea)
    % keep this one
    quadarea = A_i;

    qx = qxi;
    qy = qyi;
  end

  % plot the points, the current quad, and the best quad
%  plot(xy(:,1),xy(:,2),'ko')
%  hold on
%  plot([xy(edgesi(1:4,1),1),xy(edgesi(1:4,2),1)]', ...
%    [xy(edgesi(1:4,1),2),xy(edgesi(1:4,2),2)]','g*-','linewidth',8)
%  plot(qx,qy,'b-',qxi,qyi,'r:')
%  hold off
end

% plot the points and the quad
plot(xy(:,1),xy(:,2),'r.',qx,qy,'b-')