www.pudn.com > bsubsamp.rar > bsubsamp.m, change:2006-06-01,size:3966b

function [s, su] = bsubsamp(b, gridsep) 
%BSUBSAMP Subsample a boundary. 
%   [S, SU] = BSUBSAMP(B, GRIDSEP) subsamples the boundary B by 
%   assigning each of its points to the grid node to which it is 
%   closest.  The grid is specified by GRIDSEP, which is the 
%   separation in pixels between the grid lines. For example, if 
%   GRIDSEP = 2, there are two pixels in between grid lines. So, for 
%   instance, the grid points in the first row would be at (1,1), 
%   (1,4), (1,6), ..., and similarly in the y direction. The value 
%   of GRIDSEP must be an even integer. The boundary is specified by 
%   a set of coordinates in the form of an np-by-2 array.  It is 
%   assumed that the boundary is one pixel thick.  
%   Output S is the subsampled boundary. Output SU is normalized so 
%   that the grid separation is unity.  This is useful for obtaining 
%   the Freeman chain code of the subsampled boundary. 
%   Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins 
%   Digital Image Processing Using MATLAB, Prentice-Hall, 2004 
%   $Revision: 1.8 $  $Date: 2004/11/04 20:17:59 $ 
% Check input. 
[np, nc] = size(b); 
if np < nc  
   error('B must be of size np-by-2.');  
if gridsep/2 ~= round(gridsep/2)  
   error('GRIDSEP must be an even integer.') 
% Some boundary tracing programs, such as boundaries.m, end with  
% the beginning, resulting in a sequence in which the coordinates 
% of the first and last points are the same. If this is the case  
% in b, eliminate the last point. 
if isequal(b(1, :), b(np, :)) 
   np = np - 1; 
   b = b(1:np, :); 
% Find the max x and y spanned by the boundary. 
xmax = max(b(:, 1)); 
ymax = max(b(:, 2)); 
% Determine the integral number of grid lines with gridsep points in 
% between them that encompass the intervals [1,xmax], [1,ymax]. 
GLx = ceil((xmax + gridsep)/(gridsep + 1)); 
GLy = ceil((ymax + gridsep)/(gridsep + 1)); 
% Form vectors of x and y grid locations. 
I = 1:GLx; 
% Vector of grid line locations intersecting x-axis. 
X(I) = gridsep*I + (I - gridsep);  
J = 1:GLy; 
% Vector of grid line locations intersecting y-axis. 
Y(J) = gridsep*J + (J - gridsep);  
% Compute both components of the cityblock distance between each 
% element of b and all the grid-line intersections.  Assign each 
% point to the grid location for which each comp of the cityblock 
% distance was <= gridsep/2. Note the use of meshgrid to 
% optimize the code. Keep in mind that meshgrid requires that columns 
% be listed first (see Chapter 2). 
DIST = gridsep/2; 
[YG, XG] = meshgrid(Y, X); 
Q = 1; 
for k=1:np 
   [I,J] = find(abs(XG - b(k, 1)) <= DIST & abs(YG - b(k, 2)) <= ... 
   % If point b(k,:) is equidistant from two or more grid intersections, 
   % assign the point arbitrarily to the first one: 
   I = I(1); 
   J = J(1); 
   ord = k; % To keep track of order of input coordinates. 
   d1(Q, :) = cat(2, Y(J), ord); 
   d2(Q, :) = cat(2, X(I), ord); 
   Q = Q + 1; 
% d is the set of points assigned to the new grid with line 
% separation of gridsep. Note that it is formed as d=(d2,d1) to 
% compensate for the coordinate transposition inherent in using 
% meshgrid (see Chapter 2).  
d = cat(2, d2(:, 1), d1); % The second column of d1 & d2 is ord. 
% Sort the points using the values in ord, which is the last col in 
% d. 
d = fliplr(d); % So the last column becomes first. 
d = sortrows(d); 
d = fliplr(d); % Flip back. 
% Eliminate duplicate rows in the first two components of  
% d to create the output. The cw or ccw order MUST be preserved. 
s = d(:, 1:2); 
[s, m, n] = unique(s, 'rows'); 
% Function unique sorts the data--Restore to original order 
% by using the contents of m. 
s = [s, m]; 
s = fliplr(s); 
s = sortrows(s); 
s = fliplr(s); 
s = s(:, 1:2); 
% Scale to unit grid so that can use directly to obtain Freeman 
% chain code.  The shape does not change. 
su = round(s./gridsep) + 1;