minperpoly.m

%源程序
function [x, y] = minperpoly(B, cellsize) 
%MINPERPOLY Computes the minimum perimeter polygon. 
%   [X, Y] = MINPERPOLY(F, CELLSIZE) computes the vertices in [X, Y] 
%   of the minimum perimeter polygon of a single binary region or 
%   boundary in image B. The procedure is based on Slansky's 
%   shrinking rubber band approach. Parameter CELLSIZE determines the 
%   size of the square cells that enclose the boundary of the region 
%   in B. CELLSIZE must be a nonzero integer greater than 1. 
%    
%   The algorithm is applicable only to boundaries that are not 
%   self-intersecting and that do not have one-pixel-thick 
%   protrusions.  
 
%   Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins 
%   Digital Image Processing Using MATLAB, Prentice-Hall, 2004 
%   $Revision: 1.6 $  $Date: 2003/11/21 14:41:39 $ 
 
if cellsize <= 1 
   error('CELLSIZE must be an integer > 1.');  
end 
 
% Fill B in case the input was provided as a boundary. Later  
% the boundary will be extracted with 4-connectivity, which  
% is required by the algorithm. The use of bwperim assures  
% that 4-connectivity is preserved at this point.                
B = imfill(B, 'holes'); 
B = bwperim(B);                
[M, N] = size(B); 
 
% Increase image size so that the image is of size K-by-K 
% with (a) K >= max(M,N) and (b)  K/cellsize = a power of 2. 
K = nextpow2(max(M, N)/cellsize); 
K = (2^K)*cellsize; 
     
% Increase image size to nearest integer power of 2, by  
% appending zeros to the end of the image. This will allow  
% quadtree  decompositions as small as cells of size 2-by-2,  
% which is the smallest allowed value of cellsize. 
M = K - M; 
N = K - N; 
B = padarray(B, [M N], 'post'); % f is now of size K-by-K 
 
% Quadtree decomposition. 
Q = qtdecomp(B, 0, cellsize);  
 
% Get all the subimages of size cellsize-by-cellsize. 
[vals, r, c] = qtgetblk(B, Q, cellsize); 
 
% Get all the subimages that contain at least one black  
% pixel. These are the cells of the wall enclosing the boundary. 
I = find(sum(sum(vals(:, :, :)) >= 1)); 
x = r(I); 
y = c(I); 
 
% [x', y'] is a length(I)-by-2 array.  Each member of this array is 
% the left, top corner of a black cell of size cellsize-by-cellsize. 
% Fill the cells with black to form a closed border of black cells 
% around interior points. These cells are the cellular complex. 
for k = 1:length(I) 
   B(x(k):x(k) + cellsize-1, y(k):y(k) + cellsize-1) = 1; 
end 
     
BF = imfill(B, 'holes'); 
 
% Extract the points interior to the black border. This is the region 
% of interest around which the MPP will be found.  
B = BF & (~B); 
     
% Extract the 4-connected boundary. 
B = boundaries(B, 4, 'cw'); 
% Find the largest one in case of parasitic regions. 
J = cellfun('length', B); 
I = find(J == max(J)); 
B = B{I(1)}; 
 
% Function boundaries outputs the last coordinate pair equal to the 
% first.  Delete it.  
B = B(1:end-1,:); 
                   
% Obtain the xy coordinates of the boundary. 
x = B(:, 1); 
y = B(:, 2); 
 
% Find the smallest x-coordinate and corresponding 
% smallest y-coordinate.   
cx = find(x == min(x)); 
cy = find(y == min(y(cx))); 
 
% The cell with top leftmost corner at (x1, y1) below is the first 
% point considered by the algorithm.  The remaining points are 
% visited in the clockwise direction starting at (x1, y1).  
x1 = x(cx(1)); 
y1 = y(cy(1)); 
 
% Scroll data so that the first point is (x1, y1). 
I = find(x == x1 & y == y1); 
x = circshift(x, [-(I - 1), 0]); 
y = circshift(y, [-(I - 1), 0]); 
 
% The same shift applies to B. 
B = circshift(B, [-(I - 1), 0]); 
 
% Get the Freeman chain code.  The first row of B is the required 
% starting point. The first element of the code is the transition 
% between the 1st and 2nd element of B, the second element of  
% the code is the transition between the 2nd and 3rd elements of B, 
% and so on.  The last element of the code is the transition between 
% the last and 1st elements of B. The elements of B form a cw 
% sequence (see above), so we use 'same' for the direction in 
% function fchcode.  
code = fchcode(B, 4, 'same');  
code = code.fcc; 
 
% Follow the code sequence to extract the Black Dots, BD, (convex 
% corners) and White Dots, WD, (concave corners). The transitions are 
% as follows: 0-to-1=WD; 0-to-3=BD; 1-to-0=BD; 1-to-2=WD; 2-to-1=BD; 
% 2-to-3=WD; 3-to-0=WD; 3-to-2=dot.  The formula t = 2*first - second 
% gives the following unique values for these transitions: -1, -3, 2, 
% 0, 3, 1, 6, 4.  These are applicable to travel in the cw direction.  
% The WD's are displaced one-half a diagonal from the BD's to form 
% the half-cell expansion required in the algorithm.  
 
% Vertices will be computed as array "vertices" of dimension nv-by-3, 
% where nv is the number of vertices. The first two elements of any 
% row of array vertices are the (x,y) coordinates of the vertex 
% corresponding to that row, and the third element is 1 if the 
% vertex is convex (BD) or 2 if it is concave (WD). The first vertex 
% is known to be convex, so it is black.  
vertices = [x1, y1, 1]; 
n = 1; 
k = 1; 
for k = 2:length(code) 
   if code(k - 1) ~= code(k) 
      n = n + 1; 
      t = 2*code(k-1) - code(k); % t = value of formula. 
      if t == -3 | t == 2 | t == 3 | t == 4 % Convex: Black Dots. 
         vertices(n, 1:3) = [x(k), y(k), 1]; 
      elseif t == -1 | t == 0 | t == 1 | t == 6 % Concave: White Dots. 
         if t == -1 
            vertices(n, 1:3) = [x(k) - cellsize, y(k) - cellsize,2]; 
         elseif t==0 
            vertices(n, 1:3) = [x(k) + cellsize, y(k) - cellsize,2]; 
         elseif t==1 
            vertices(n, 1:3) = [x(k) + cellsize, y(k) + cellsize,2]; 
         else 
            vertices(n, 1:3) = [x(k) - cellsize, y(k) + cellsize,2]; 
         end 
      else 
         % Nothing to do here. 
      end 
   end 
end 
 
% The rest of minperpoly.m processes the vertices to 
% arrive at the MPP. 
 
flag = 1; 
while flag 
   % Determine which vertices lie on or inside the 
   % polygon whose vertices are the Black Dots. Delete all  
   % other points. 
   I = find(vertices(:, 3) == 1); 
   xv = vertices(I, 1); % Coordinates of the Black Dots. 
   yv = vertices(I, 2); 
   X = vertices(:, 1); % Coordinates of all vertices. 
   Y = vertices(:, 2); 
   IN = inpolygon(X, Y, xv, yv); 
   I = find(IN ~= 0);         
   vertices = vertices(I, :); 
    
   % Now check for any Black Dots that may have been turned into 
   % concave vertices after the previous deletion step. Delete 
   % any such Black Dots and recompute the polygon as in the 
   % previous section of code. When no more changes occur, set 
   % flag to 0, which causes the loop to terminate.  
   x = vertices(:, 1); 
   y = vertices(:, 2); 
   angles = polyangles(x, y); % Find all the interior angles. 
   I = find(angles > 180 & vertices(:, 3) == 1);   
   if isempty(I) 
      flag = 0; 
   else 
      J = 1:length(vertices); 
      for k = 1:length(I) 
         K = find(J ~= I(k)); 
         J = J(K); 
      end 
      vertices = vertices(J, :); 
   end 
end 
 
% Final pass to delete the vertices with angles of 180 degrees. 
x = vertices(:, 1); 
y = vertices(:, 2); 
angles = polyangles(x, y); 
I = find(angles ~= 180); 
 
% Vertices of the MPP: 
x = vertices(I, 1); 
y = vertices(I, 2);