目录
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格式返回该图像作为分割结果。