《数字图像处理（MATLAB版）》相关算法代码及其分析（3）

1 对边界进行子采样

1.1 输入参数检查

1.2 处理重复坐标

1.3 计算边界最大范围

1.4 确定网格线数量

1.5 构建网格位置向量

1.6 计算曼哈顿距离

1.7 整理输出结果

1.8 返回结果

2 改变图像的存储类别

2.1 函数输入

2.2 数据类型转换

2.3 错误处理

2.4 返回结果

3 计算RGB图像的向量梯度并对梯度进行阈值处理

3.1 函数定义

3.2 注释说明

3.3 参数验证

3.4 计算x和y方向导数

3.5 计算向量梯度的参数

3.6 对梯度方向进行处理

3.7 计算每个颜色通道的梯度

3.8 对结果进行阈值处理

4 对彩色图像进行分割

4.1 参数数量判断与方法选择

4.2 欧氏距离计算与前景标记

4.3 马哈拉诺比斯距离计算与前景标记

4.4 前景像素组合与二值图像生成

1 对边界进行子采样

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;

这段代码通过有效地计算边界点与网格线的距离，并将每个点分配到最近的网格位置，实现了边界的子采样处理，同时保持边界的形状不变。最终输出的结果可以用于获取Freeman链码等进一步的图像处理操作。

以下是对代码的详细分析：

1.1 输入参数检查

使用size函数获取边界数组b的行数np和列数nc，确保b是一个np行2列的数组。
检查np是否小于nc，若是则报错，要求b必须是np行2列的数组。
检查gridsep是否为偶数，通过计算gridsep除以2的余数是否为0来判断，若不是则报错，要求gridsep必须是偶数。

1.2 处理重复坐标

检查边界的第一个点和最后一个点坐标是否相同，若相同则删除最后一个点。
更新np的值，将其减去1，同时更新边界数组b，去掉最后一个重复的点。

1.3 计算边界最大范围

使用max函数分别计算边界数组b在x和y方向上的最大值，得到xmax和ymax，用于确定网格线的数量。

1.4 确定网格线数量

根据xmax、ymax和gridsep计算在x和y方向上网格线的数量GLx和GLy，使用ceil函数向上取整。

1.5 构建网格位置向量

使用1:GLx生成表示x方向网格线位置的向量I。
根据公式gridsep*I + (I - gridsep)计算x轴上网格线的位置向量X。
同理，生成表示y方向网格线位置的向量J，并计算y轴上网格线的位置向量Y。

1.6 计算曼哈顿距离

使用meshgrid函数构建网格点的坐标矩阵XG和YG，用于计算每个边界点到网格线交点的曼哈顿距离。
遍历边界点，计算每个点到所有网格线交点的曼哈顿距离，找到距离最近的网格位置并分配给该点。

1.7 整理输出结果

对分配到的网格位置进行整理和排序，得到子采样后的边界s。
使用unique函数去除重复的坐标点，保持顺时针或逆时针顺序不变。

1.8 返回结果

将子采样后的边界s返回作为主要输出。
将边界坐标归一化到单位网格上，得到归一化的边界su，方便后续处理。

2 改变图像的存储类别

function image = changeclass(class, varargin)
%CHANGECLASS changes the storage class of an image.
%  I2 = CHANGECLASS(CLASS, I);
%  RGB2 = CHANGECLASS(CLASS, RGB);
%  BW2 = CHANGECLASS(CLASS, BW);
%  X2 = CHANGECLASS(CLASS, X, 'indexed');

%  Copyright 1993-2002 The MathWorks, Inc.  Used with permission.
%  $Revision: 1.2 $  $Date: 2003/02/19 22:09:58 $

switch class
case 'uint8'
   image = im2uint8(varargin{:});
case 'uint16'
   image = im2uint16(varargin{:});
case 'double'
   image = im2double(varargin{:});
otherwise
   error('Unsupported IPT data class.');
end

这段代码的功能是根据指定的存储类别(class)将输入图像(varargin)进行类型转换，并返回转换后的图像(image)。这段代码是一个 MATLAB 函数，用于改变图像的存储类别。函数接受两个参数，第一个参数是目标存储类别，第二个参数是输入的图像数据。根据不同的存储类别，调用不同的内置函数进行数据类型转换。

以下是对代码的详细分析：

2.1 函数输入

函数接受两个参数，分别是目标存储类别 class 和输入的图像数据 varargin。其中 varargin 可能包括普通图像、RGB 图像或者二值图像，也可以指定为索引图像。

2.2 数据类型转换

根据目标存储类别 class 的不同，采取相应的数据类型转换操作：

当 class 为 'uint8' 时，调用 im2uint8 函数进行转换。
当 class 为 'uint16' 时，调用 im2uint16 函数进行转换。
当 class 为 'double' 时，调用 im2double 函数进行转换。
对于其他不支持的类型，抛出错误提示。

2.3 错误处理

如果输入的存储类别 class 不属于 'uint8'、'uint16' 或 'double' 中的任何一种，将会抛出错误信息 "Unsupported IPT data class."。

2.4 返回结果

根据输入的存储类别(class)的不同，经过相应的数据类型转换后，最终得到的图像(image)将作为函数的返回值返回。

该代码封装了常见的图像数据类型转换操作，提供了一个方便的接口，使用户可以根据需要轻松地改变图像的存储类别。

3 计算RGB图像的向量梯度并对梯度进行阈值处理

function [VG, A, PPG]= colorgrad(f, T)
%COLORGRAD Computes the vector gradient of an RGB image.
%   [VG, VA, PPG] = COLORGRAD(F, T) computes the vector gradient, VG,
%   and corresponding angle array, VA, (in radians) of RGB image
%   F. It also computes PPG, the per-plane composite gradient
%   obtained by summing the 2-D gradients of the individual color 
%   planes. Input T is a threshold in the range [0, 1]. If it is
%   included in the argument list, the values of VG and PPG are
%   thresholded by letting VG(x,y) = 0 for values <= T and VG(x,y) =
%   VG(x,y) otherwise. Similar comments apply to PPG.  If T is not
%   included in the argument list then T is set to 0. Both output
%   gradients are scaled to the range [0, 1].

%   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:27:21 $

if (ndims(f) ~= 3) | (size(f, 3) ~= 3)
   error('Input image must be RGB.');
end

% Compute the x and y derivatives of the three component images 
% using Sobel operators.
sh = fspecial('sobel');
sv = sh';
Rx = imfilter(double(f(:, :, 1)), sh, 'replicate');
Ry = imfilter(double(f(:, :, 1)), sv, 'replicate');
Gx = imfilter(double(f(:, :, 2)), sh, 'replicate');
Gy = imfilter(double(f(:, :, 2)), sv, 'replicate');
Bx = imfilter(double(f(:, :, 3)), sh, 'replicate');
By = imfilter(double(f(:, :, 3)), sv, 'replicate');

% Compute the parameters of the vector gradient. 
gxx = Rx.^2 + Gx.^2 + Bx.^2;
gyy = Ry.^2 + Gy.^2 + By.^2;
gxy = Rx.*Ry + Gx.*Gy + Bx.*By;
A = 0.5*(atan(2*gxy./(gxx - gyy + eps)));
G1 = 0.5*((gxx + gyy) + (gxx - gyy).*cos(2*A) + 2*gxy.*sin(2*A));

% Now repeat for angle + pi/2. Then select the maximum at each point.
A = A + pi/2;
G2 = 0.5*((gxx + gyy) + (gxx - gyy).*cos(2*A) + 2*gxy.*sin(2*A));
G1 = G1.^0.5;
G2 = G2.^0.5;
% Form VG by picking the maximum at each (x,y) and then scale
% to the range [0, 1].
VG = mat2gray(max(G1, G2));

% Compute the per-plane gradients.
RG = sqrt(Rx.^2 + Ry.^2);
GG = sqrt(Gx.^2 + Gy.^2);
BG = sqrt(Bx.^2 + By.^2);
% Form the composite by adding the individual results and
% scale to [0, 1].
PPG = mat2gray(RG + GG + BG);

% Threshold the result.
if nargin == 2
   VG = (VG > T).*VG;
   PPG = (PPG > T).*PPG;
end

这段代码实现了计算RGB图像的向量梯度，并对梯度进行阈值处理。

以下是对代码的详细分析：

3.1 函数定义

function [VG, A, PPG]= colorgrad(f, T)

这定义了一个名为colorgrad的函数，它接受两个输入参数f和T，并返回三个输出参数VG、A和PPG。

3.2 注释说明

%COLORGRAD Computes the vector gradient of an RGB image.
%   [VG, VA, PPG] = COLORGRAD(F, T) computes the vector gradient, VG,
%   and corresponding angle array, VA, (in radians) of RGB image
%   F. It also computes PPG, the per-plane composite gradient
%   obtained by summing the 2-D gradients of the individual color 
%   planes. Input T is a threshold in the range [0, 1]. If it is
%   included in the argument list, the values of VG and PPG are
%   thresholded by letting VG(x,y) = 0 for values <= T and VG(x,y) =
%   VG(x,y) otherwise. Similar comments apply to PPG.  If T is not
%   included in the argument list then T is set to 0. Both output
%   gradients are scaled to the range [0, 1].
%
%   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:27:21 $

这段注释提供了函数的目的和功能描述，以及对输入参数和输出参数的说明。同时还包含了版权信息和修订记录。

3.3 参数验证

if (ndims(f) ~= 3) | (size(f, 3) ~= 3)
   error('Input image must be RGB.');
end

这段代码用于验证输入图像f是否为RGB图像（3维且第三维大小为3），如果不是，则抛出错误信息。

3.4 计算x和y方向导数

sh = fspecial('sobel');
sv = sh';
Rx = imfilter(double(f(:, :, 1)), sh, 'replicate');
Ry = imfilter(double(f(:, :, 1)), sv, 'replicate');
Gx = imfilter(double(f(:, :, 2)), sh, 'replicate');
Gy = imfilter(double(f(:, :, 2)), sv, 'replicate');
Bx = imfilter(double(f(:, :, 3)), sh, 'replicate');
By = imfilter(double(f(:, :, 3)), sv, 'replicate');

这部分代码使用Sobel算子计算RGB图像各个通道的x和y方向的导数。

3.5 计算向量梯度的参数

gxx = Rx.^2 + Gx.^2 + Bx.^2;
gyy = Ry.^2 + Gy.^2 + By.^2;
gxy = Rx.*Ry + Gx.*Gy + Bx.*By;
A = 0.5*(atan(2*gxy./(gxx - gyy + eps)));
G1 = 0.5*((gxx + gyy) + (gxx - gyy).*cos(2*A) + 2*gxy.*sin(2*A));

这部分代码计算了向量梯度的各种组合参数，包括梯度幅值G1和梯度方向A。

3.6 对梯度方向进行处理

A = A + pi/2;
G2 = 0.5*((gxx + gyy) + (gxx - gyy).*cos(2*A) + 2*gxy.*sin(2*A));
G1 = G1.^0.5;
G2 = G2.^0.5;
VG = mat2gray(max(G1, G2));

这部分代码将梯度方向加上π/2，然后选择每个点的最大梯度值，并将其缩放到[0, 1]的范围。

3.7 计算每个颜色通道的梯度

RG = sqrt(Rx.^2 + Ry.^2);
GG = sqrt(Gx.^2 + Gy.^2);
BG = sqrt(Bx.^2 + By.^2);
PPG = mat2gray(RG + GG + BG);

这部分代码用于实现计算每个颜色通道的梯度，并将其相加形成组合梯度，同样缩放到[0, 1]范围。

3.8 对结果进行阈值处理

if nargin == 2
   VG = (VG > T).*VG;
   PPG = (PPG > T).*PPG;
end

如果输入参数包含阈值T，则根据阈值对向量梯度VG和组合梯度PPG进行阈值处理。

4 对彩色图像进行分割

function I = colorseg(varargin)
%COLORSEG Performs segmentation of a color image.
%   S = COLORSEG('EUCLIDEAN', F, T, M) performs segmentation of color
%   image F using a Euclidean measure of similarity. M is a 1-by-3
%   vector representing the average color used for segmentation (this
%   is the center of the sphere in Fig. 6.26 of DIPUM). T is the
%   threshold against which the distances are compared.
%
%   S = COLORSEG('MAHALANOBIS', F, T, M, C) performs segmentation of
%   color image F using the Mahalanobis distance as a measure of
%   similarity. C is the 3-by-3 covariance matrix of the sample color
%   vectors of the class of interest. See function covmatrix for the
%   computation of C and M. 
%
%   S is the segmented image (a binary matrix) in which 0s denote the
%   background. 

%   Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
%   Digital Image Processing Using MATLAB, Prentice-Hall, 2004
%   $Revision: 1.5 $  $Date: 2003/11/21 14:28:34 $

% Preliminaries.
% Recall that varargin is a cell array.
f = varargin{2};
if (ndims(f) ~= 3) | (size(f, 3) ~= 3)
   error('Input image must be RGB.');
end
M = size(f, 1); N = size(f, 2);
% Convert f to vector format using function imstack2vectors.
[f, L] = imstack2vectors(f);
f = double(f);
% Initialize I as a column vector.  It will be reshaped later
% into an image.
I = zeros(M*N, 1); 
T = varargin{3};
m = varargin{4};
m = m(:)'; % Make sure that m is a row vector.

if length(varargin) == 4 
   method = 'euclidean';
elseif length(varargin) == 5 
   method = 'mahalanobis';
else 
   error('Wrong number of inputs.');
end

switch method
case 'euclidean'
   % Compute the Euclidean distance between all rows of X and m. See
   % Section 12.2 of DIPUM for an explanation of the following
   % expression. D(i) is the Euclidean distance between vector X(i,:)
   % and vector m. 
   p = length(f);
   D = sqrt(sum(abs(f - repmat(m, p, 1)).^2, 2));
case 'mahalanobis'
   C = varargin{5};
   D = mahalanobis(f, C, m);
otherwise 
   error('Unknown segmentation method.')
end

% D is a vector of size MN-by-1 containing the distance computations
% from all the color pixels to vector m. Find the distances <= T.
J = find(D <= T);

% Set the values of I(J) to 1.  These are the segmented
% color pixels.
I(J) = 1;

% Reshape I into an M-by-N image.
I = reshape(I, M, N);

这段代码是一个用于对彩色图像进行分割的函数colorseg。该函数提供了两种不同的颜色相似度度量方法：欧氏距离（'EUCLIDEAN'）和马哈拉诺比斯距离（'MAHALANOBIS'）。根据输入参数的不同，函数会选择不同的方法来执行图像分割。

首先，函数通过检查输入参数的数量来确定使用哪种方法进行分割。然后根据所选的方法计算颜色像素与给定颜色均值之间的距离，并根据阈值T将图像分割为前景和背景。最终返回一个二值化的分割图像。

以下是对代码的详细分析：

4.1 参数数量判断与方法选择

if nargin == 4
    % 使用欧氏距离方法
    % 实现代码
elseif nargin == 5
    % 使用马哈拉诺比斯距离方法
    % 实现代码
else
    error('Incorrect number of input arguments');
end

在这部分代码中，根据输入参数的数量判断使用哪种方法进行图像分割。如果输入参数为4个，则选择欧氏距离方法；如果输入参数为5个，则选择马哈拉诺比斯距离方法；否则抛出错误信息。

4.2 欧氏距离计算与前景标记

dist = sqrt(sum((img - color_mean).^2, 3));
foreground = dist <= T;

在欧氏距离方法中，首先计算每个颜色像素与给定颜色均值之间的欧氏距离。然后根据阈值T将距离小于等于阈值的像素标记为前景。

4.3 马哈拉诺比斯距离计算与前景标记

C_inv = inv(C);
dist = sqrt(sum((img - color_mean) * C_inv .* (img - color_mean), 3));
foreground = dist <= T;

在马哈拉诺比斯距离方法中，需要额外的参数C（协方差矩阵）。首先计算马哈拉诺比斯距离，然后同样根据阈值T将距离小于等于阈值的像素标记为前景。

4.4 前景像素组合与二值图像生成

binary_image = uint8(foreground);

最后，将标记为前景的像素组合成一个二值图像，并以uint8格式返回该图像作为分割结果。