www.pudn.com > lianma.rar > bsubsamp.m, change:2004-12-16,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.'); end if gridsep/2 ~= round(gridsep/2) error('GRIDSEP must be an even integer.') end % 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, :); end % 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)) <= ... DIST); % 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; end % 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;