用于查找不连续的最常用的方法是对图像运行掩膜计算。就是将掩膜区与所包含的图像的灰度级乘积的和。
在该模板下，给定一个非负阈值，设|R|≥T则可以说该点被检测出来：
对于掩膜最重要的要求是，在一个孤立点的时候，掩膜的响应必须最强，而在亮度不变的区域中，响应为零。
可以通过该方法实现g = abs(imfilter(double(f), w)) >= T
为了便于观察，图片被放大：
f = imread('./point.tif');
w = [-1 -1 -1; -1 8 -1; -1 -1 -1];
subplot(1, 2, 1);
imshow(f);
g = abs(imfilter(double(f), w));
T = max(g(:));
g = g >= T;
subplot(1, 2, 2);
imshow(g);
当然还有的点检测方法是在大小为m*n的所有领域中寻找一些点，他们的最大值和最小值之差超过了T，该方法可以通过g = imsubtract(ordfilt2(f， m*n， ones(m, n)), ordfilt2(f, 1, ones(m*n))); g = g >= T实现。
也是通过掩膜计算，其掩膜为：
这个效果就不展示了。
图像中的边缘检测最好的方法还是检测亮度值的不连续性，不连续性意味着一阶二阶导数变化的比较大。那么二维函数的梯度可以写为：
∇
f
=
[
G
x
G
y
]
=
[
∂
f
∂
x
∂
f
∂
x
]
∇f = \left[ \frac{G_x}{Gy} \right]=\left[ \frac{\frac{∂f}{∂x}}{\frac{∂f}{∂x}} \right]
∇f=[GyGx]=[∂x∂f∂x∂f]
其幅值为：
∇
f
=
m
a
g
(
∇
f
)
=
[
G
x
2
+
G
y
2
]
1
2
∇f = mag(∇f) = [G_x^2 + G_y^2]^{\frac{1}{2}}
∇f=mag(∇f)=[Gx2+Gy2]21
为了简化计算可以直接近似于其绝对值或平方和。
二阶导数一般用于拉普拉斯算子计算。但拉普拉斯算子自身很少被直接用做边缘检测,因为二阶导数对噪声具有无法接受的敏感性,它的幅度会产生双边缘，而且它不能检测边缘的方向。然而，正像本节稍后讨论的那样，当与其他边缘检测技术组合使用时,拉普拉斯算子是一种有效的补充方法。例如,虽然它的双边缘使得它不适合直接用于边缘检测，但该性质可用于边缘定位。
所以 边缘检测的方法为：
所以edge函数就基于该方法，其中语法为:[g, t] = edge(f, 'method', parameters)
其中一些边缘检测器为：
用Sobel算子为例, 其调用语法为：[g, t] = edge(f, 'sobel', T, dir)
其中，f是输入图像，T是一个指定的阈值，dir指定检测边缘的首选方向: ‘horizontal’、‘vertical’或’ both’（默认值)。如先前说明的那样，g是在被检测到边缘的位置处为1而在其他位置为0的逻辑类图像。输出参数t是可选的，它是函数edge所用的阈值。若指定了T的值，则t=T。否则，若T未被赋值(或为空，[])，则函数edge会令t等于它自动确定的一个阈值，然后用于边缘检测。
在输出参量中要包括t的主要原因之一是为了得到一个阈值的初始值。若使用语法g = edge(f)或[g,t] = egde (f)，则函数edge会默认使用Sobel检测器。
大一的时候，我自己给出的一个理解，在这里：PythonCV学习记录9——霍夫变换
然后hough变换得通过自己的代码实现：
function [h, theta, rho] = hough(f, dtheta, drho)
%HOUGH Hough transform.
% [H, THETA, RHO] = HOUGH(F, DTHETA, DRHO) computes the Hough
% transform of the image F. DTHETA specifies the spacing (in
% degrees) of the Hough transform bins along the theta axis. DRHO
% specifies the spacing of the Hough transform bins along the rho
% axis. H is the Hough transform matrix. It is NRHO-by-NTHETA,
% where NRHO = 2*ceil(norm(size(F))/DRHO) - 1, and NTHETA =
% 2*ceil(90/DTHETA). Note that if 90/DTHETA is not an integer, the
% actual angle spacing will be 90 / ceil(90/DTHETA).
%
% THETA is an NTHETA-element vector containing the angle (in
% degrees) corresponding to each column of H. RHO is an
% NRHO-element vector containing the value of rho corresponding to
% each row of H.
%
% [H, THETA, RHO] = HOUGH(F) computes the Hough transform using
% DTHETA = 1 and DRHO = 1.
% Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
% Digital Image Processing Using MATLAB, Prentice-Hall, 2004
% $Revision: 1.4 $ $Date: 2003/10/26 22:33:44 $
if nargin < 3
drho = 1;
end
if nargin < 2
dtheta = 1;
end
f = double(f);
[M,N] = size(f);
theta = linspace(-90, 0, ceil(90/dtheta) + 1);
theta = [theta -fliplr(theta(2:end - 1))];
ntheta = length(theta);
D = sqrt((M - 1)^2 + (N - 1)^2);
q = ceil(D/drho);
nrho = 2*q - 1;
rho = linspace(-q*drho, q*drho, nrho);
[x, y, val] = find(f);
x = x - 1; y = y - 1;
% Initialize output.
h = zeros(nrho, length(theta));
% To avoid excessive memory usage, process 1000 nonzero pixel
% values at a time.
for k = 1:ceil(length(val)/1000)
first = (k - 1)*1000 + 1;
last = min(first+999, length(x));
x_matrix = repmat(x(first:last), 1, ntheta);
y_matrix = repmat(y(first:last), 1, ntheta);
val_matrix = repmat(val(first:last), 1, ntheta);
theta_matrix = repmat(theta, size(x_matrix, 1), 1)*pi/180;
rho_matrix = x_matrix.*cos(theta_matrix) + ...
y_matrix.*sin(theta_matrix);
slope = (nrho - 1)/(rho(end) - rho(1));
rho_bin_index = round(slope*(rho_matrix - rho(1)) + 1);
theta_bin_index = repmat(1:ntheta, size(x_matrix, 1), 1);
% Take advantage of the fact that the SPARSE function, which
% constructs a sparse matrix, accumulates values when input
% indices are repeated. That抯 the behavior we want for the
% Hough transform. We want the output to be a full (nonsparse)
% matrix, however, so we call function FULL on the output of
% SPARSE.
h = h + full(sparse(rho_bin_index(:), theta_bin_index(:), ...
val_matrix(:), nrho, ntheta));
end
画出每一个点的关于ρ和θ的函数线：
f = zeros(101, 101);
f(1, 1) = 1;
f(101, 1) = 1;
f(1, 101) = 1;
f(101, 101) = 1;
f(51, 51) = 1;
[H, theta, rho] = hough(f);
imshow(H, [ ]);
最终的效果是这样：
其中曲线的交点，就是图像中的点。
其原理为：
自定义函数：houghpeaks
function [r, c, hnew] = houghpeaks(h, numpeaks, threshold, nhood)
%HOUGHPEAKS Detect peaks in Hough transform.
% [R, C, HNEW] = HOUGHPEAKS(H, NUMPEAKS, THRESHOLD, NHOOD) detects
% peaks in the Hough transform matrix H. NUMPEAKS specifies the
% maximum number of peak locations to look for. Values of H below
% THRESHOLD will not be considered to be peaks. NHOOD is a
% two-element vector specifying the size of the suppression
% neighborhood. This is the neighborhood around each peak that is
% set to zero after the peak is identified. The elements of NHOOD
% must be positive, odd integers. R and C are the row and column
% coordinates of the identified peaks. HNEW is the Hough transform
% with peak neighborhood suppressed.
%
% If NHOOD is omitted, it defaults to the smallest odd values >=
% size(H)/50. If THRESHOLD is omitted, it defaults to
% 0.5*max(H(:)). If NUMPEAKS is omitted, it defaults to 1.
% 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 13:34:50 $
if nargin < 4
nhood = size(h)/50;
% Make sure the neighborhood size is odd.
nhood = max(2*ceil(nhood/2) + 1, 1);
end
if nargin < 3
threshold = 0.5 * max(h(:));
end
if nargin < 2
numpeaks = 1;
end
done = false;
hnew = h; r = []; c = [];
while ~done
[p, q] = find(hnew == max(hnew(:)));
p = p(1); q = q(1);
if hnew(p, q) >= threshold
r(end + 1) = p; c(end + 1) = q;
% Suppress this maximum and its close neighbors.
p1 = p - (nhood(1) - 1)/2; p2 = p + (nhood(1) - 1)/2;
q1 = q - (nhood(2) - 1)/2; q2 = q + (nhood(2) - 1)/2;
[pp, qq] = ndgrid(p1:p2,q1:q2);
pp = pp(:); qq = qq(:);
% Throw away neighbor coordinates that are out of bounds in
% the rho direction.
badrho = find((pp < 1) | (pp > size(h, 1)));
pp(badrho) = []; qq(badrho) = [];
% For coordinates that are out of bounds in the theta
% direction, we want to consider that H is antisymmetric
% along the rho axis for theta = +/- 90 degrees.
theta_too_low = find(qq < 1);
qq(theta_too_low) = size(h, 2) + qq(theta_too_low);
pp(theta_too_low) = size(h, 1) - pp(theta_too_low) + 1;
theta_too_high = find(qq > size(h, 2));
qq(theta_too_high) = qq(theta_too_high) - size(h, 2);
pp(theta_too_high) = size(h, 1) - pp(theta_too_high) + 1;
% Convert to linear indices to zero out all the values.
hnew(sub2ind(size(hnew), pp, qq)) = 0;
done = length(r) == numpeaks;
else
done = true;
end
end
一旦在Hough变换中识别出了一组候选的峰被，则它们还要留待确定是否存在与这些峰值相关的线段以及它们的起始和结束位置。对每一个峰值来说,第一步是找到图像中影响到峰值的每一个非零值点的位置。为达到这一目的，我们编写了函数houghpixels。
function [r, c] = houghpixels(f, theta, rho, rbin, cbin)
%HOUGHPIXELS Compute image pixels belonging to Hough transform bin.
% [R, C] = HOUGHPIXELS(F, THETA, RHO, RBIN, CBIN) computes the
% row-column indices (R, C) for nonzero pixels in image F that map
% to a particular Hough transform bin, (RBIN, CBIN). RBIN and CBIN
% are scalars indicating the row-column bin location in the Hough
% transform matrix returned by function HOUGH. THETA and RHO are
% the second and third output arguments from the HOUGH function.
% Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
% Digital Image Processing Using MATLAB, Prentice-Hall, 2004
% $Revision: 1.4 $ $Date: 2003/10/26 22:35:03 $
[x, y, val] = find(f);
x = x - 1; y = y - 1;
theta_c = theta(cbin) * pi / 180;
rho_xy = x*cos(theta_c) + y*sin(theta_c);
nrho = length(rho);
slope = (nrho - 1)/(rho(end) - rho(1));
rho_bin_index = round(slope*(rho_xy - rho(1)) + 1);
idx = find(rho_bin_index == rbin);
r = x(idx) + 1; c = y(idx) + 1;
然后通过这段代码找到的相关像素组合成线段，所以得再编写一个houghlines，其过程为：
function lines = houghlines(f,theta,rho,rr,cc,fillgap,minlength)
%HOUGHLINES Extract line segments based on the Hough transform.
% LINES = HOUGHLINES(F, THETA, RHO, RR, CC, FILLGAP, MINLENGTH)
% extracts line segments in the image F associated with particular
% bins in a Hough transform. THETA and RHO are vectors returned by
% function HOUGH. Vectors RR and CC specify the rows and columns
% of the Hough transform bins to use in searching for line
% segments. If HOUGHLINES finds two line segments associated with
% the same Hough transform bin that are separated by less than
% FILLGAP pixels, HOUGHLINES merges them into a single line
% segment. FILLGAP defaults to 20 if omitted. Merged line
% segments less than MINLENGTH pixels long are discarded.
% MINLENGTH defaults to 40 if omitted.
%
% LINES is a structure array whose length equals the number of
% merged line segments found. Each element of the structure array
% has these fields:
%
% point1 End-point of the line segment; two-element vector
% point2 End-point of the line segment; two-element vector
% length Distance between point1 and point2
% theta Angle (in degrees) of the Hough transform bin
% rho Rho-axis position of the Hough transform bin
% Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
% Digital Image Processing Using MATLAB, Prentice-Hall, 2004
% $Revision: 1.4 $ $Date: 2003/10/26 22:34:10 $
if nargin < 6
fillgap = 20;
end
if nargin < 7
minlength = 40;
end
numlines = 0; lines = struct;
for k = 1:length(rr)
rbin = rr(k); cbin = cc(k);
% Get all pixels associated with Hough transform cell.
[r, c] = houghpixels(f, theta, rho, rbin, cbin);
if isempty(r)
continue
end
% Rotate the pixel locations about (1,1) so that they lie
% approximately along a vertical line.
omega = (90 - theta(cbin)) * pi / 180;
T = [cos(omega) sin(omega); -sin(omega) cos(omega)];
xy = [r - 1 c - 1] * T;
x = sort(xy(:,1));
% Find the gaps larger than the threshold.
diff_x = [diff(x); Inf];
idx = [0; find(diff_x > fillgap)];
for p = 1:length(idx) - 1
x1 = x(idx(p) + 1); x2 = x(idx(p + 1));
linelength = x2 - x1;
if linelength >= minlength
point1 = [x1 rho(rbin)]; point2 = [x2 rho(rbin)];
% Rotate the end-point locations back to the original
% angle.
Tinv = inv(T);
point1 = point1 * Tinv; point2 = point2 * Tinv;
numlines = numlines + 1;
lines(numlines).point1 = point1 + 1;
lines(numlines).point2 = point2 + 1;
lines(numlines).length = linelength;
lines(numlines).theta = theta(cbin);
lines(numlines).rho = rho(rbin);
end
end
end
举个例子：
f = imread('buliding.tif');
[H, theta, rho] = hough(f);
[r, c, hnew] = houghpeaks(H, 5);
lines = houghlines(f, theta, rho, r, c);
imshow(f), hold on
for k = 1:length(lines)
xy = [lines(k).point1; lines(k).point2];
plot(xy(:, 2), xy(:, 1), 'Linewidth', 4, 'Color', [.6 .6 .6]);
end
除了我们手动指定阈值，还可以将阈值算法写入程序中，让他自定义计算阈值
选取阈值的一种方法是目视检查直方图。所示的直方图有两个不同的模式;正如结果所示，我们可很容易地选取一个阈值T来分开它们。另一种选择T的方法是反复试验，挑选不同的阈值，直到观测者觉得产生了较好的结果时为止。这在交互式环境下特别有效，例如，这种方法允许用户使用一个图形控件（如滑动条）来改变阈值并可立即看到结果。比如图示这样的时候：
当然算法可以这样迭代：
graythresh函数，通过Otsu方法计算阈值。具体过程可以百度，是计算选择最大类间方差σB2的阈值k实现的，调用方法：T = graythresh(f)比如这个图，它的直方图为：
我们通过全局阈值算法或者Otsu算法来计算阈值：
f = imread('calc.tif');
T = 0.5*(double(min(f(:))) + double(max(f(:))));
done = false;
while ~done
g = f >= T;
Tnext = 0.5*(mean(f(g)) + mean(f(~g)));
done = abs(T - Tnext) < 0.5;
T = Tnext;
end
T2 = graythresh(f) * 255;
display([T, T2]);
% 140.5727 140.0000
有些情况，全局阈值处理方法在背景照明不均匀时有可能无效。在这种情况下,一种常用的处理方法是针对照明问题做预处理以补偿图像,然后再对预处理后的图像采用全局阈值处理。通过应用一个形态学顶帽算子并对得到的结果使用函数graythresh来计算的。我们可以证明这种处理等同于使用局部变化的阈值函数T(x,y)对f(x, y)进行阈值处理:
g
(
x
,
y
)
=
x
=
{
1
若
f
(
x
,
y
)
≥
T
(
x
,
y
)
0
若
f
(
x
,
y
)
<
T
(
x
,
y
)
g(x,y) = x = \begin{cases} 1 &\text{若} f(x,y)≥T(x,y) \\ 0 &\text{若} f(x,y)<T(x,y) \end{cases}
g(x,y)=x={10若f(x,y)≥T(x,y)若f(x,y)<T(x,y)
其中：
T
(
x
,
y
)
=
f
o
(
x
,
y
)
+
T
o
T(x,y)=f_o(x,y)+T_o
T(x,y)=fo(x,y)+To
图像fo(x,y)是形态学开运算，常数To是对fo应用函数garythresh后的结果。
前面的区域分割是基于灰度级的不连续性来查找区域之间的边界。而该节则是使用像素的属性分布(亮度值等)的阈值来完成的。
若用R表示整个图像区域。则可以将分割看成R的n个子区域，记为R1，R2……Rn，且他们满足：
这些条件指明：
条件1指出分割必须是完全的，即每个点都必须在一个区域中。条件2要求区域中的点应该是按预先定义好的方式(如4连接或8连接）连接的。条件3说明区域之间不相交。条件4说明分割区域中的像素点必须满足的性质——例如，若所有Ri中的像素有相同的灰度级，则P(Ri)=TRUE。最后，条件5指出，邻近区域Ri和Rj在谓词P上的意义是不-样的
顾名思义，区域生长是根据预先定义的生长准则来把像素或子区域集合成较大区域的处理方法。基本处理方法是以–组“种子”点开始来形成生长区域，即将那些预定义属性类似于种子的邻域像素附加到每个种子上（如指定的灰度级或颜色)。那么就类似于前面的形态学的重构。那么我们就需要知道开始成长的条件，以及结束成长的条件和成长过程中的迭代算法。
由一个或多个开始点组成的集合的选择通常可基于问题的性质。当没有先验的信息可用时，一种处理方法是：在每一个像素上计算网一组属性，在生长过程期间，这些属性最终将用于把像素分配到区域中。若这些计算的结果显示了一簇值,则应把具有这些特性的像系放在可以作为种子的这些簇的质心。
若在区域生长处理中未使用连通性(邻接)信息，则描述符会严生销误的错误。例知，仅使用三个不同的亮度值来显示像素的一个随机排列。若不考虑连通性而组合有着相同灰度级的像素以形成一个“区域”，则会产生无意义的分割结果。
停止规则的公式表达。一般来说,当不再有像素满足该区域所包含的准则时，生长区域的过程就会停止。准则（如亮度值、纹理、色彩)实际上是局部的，在区域生长的“历史”中不予考虑。增加区域生长算法能力的附加准则利用了大小的概念，如被生长像素和候选像素之间的相似性(如候选像素的亮度与生长区域的平均亮度的比较),正生长区域的形状等。这些类型的描述符的应用基于这样一个假设，即期望结果的模型至少是部分可用的。
所以书上写个一个regiongrow函数：
[g, NR, SI, TI] = regiongrow(f, S, T)
其中f是被分割的图像，S是一个与f大小相同的数组或是一个标量。数组中的1则是种子。若是标量，则表示一个亮度值，所有等于该亮度值的像素点都会成为种子。T也是一个与f相同大小的数组或者标量。若T是数组，则是为每一个像素点确定一个阈值，若是标量，则是全局阈值。
在输出中, g是分割后的图像，每个区域的成员都用整数标出。参数NR是不同区域的数目。参数SI是一幅包含有种子点的图像,参数TI是一幅图像，该图像中包含在经过连通性处理前通过阈值测试的像素。SI和TI的大小均与f相同。
function [g, NR, SI, TI] = regiongrow(f, S, T)
%REGIONGROW Perform segmentation by region growing.
% [G, NR, SI, TI] = REGIONGROW(F, SR, T). S can be an array (the
% same size as F) with a 1 at the coordinates of every seed point
% and 0s elsewhere. S can also be a single seed value. Similarly,
% T can be an array (the same size as F) containing a threshold
% value for each pixel in F. T can also be a scalar, in which
% case it becomes a global threshold.
%
% On the output, G is the result of region growing, with each
% region labeled by a different integer, NR is the number of
% regions, SI is the final seed image used by the algorithm, and TI
% is the image consisting of the pixels in F that satisfied the
% threshold test.
% Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
% Digital Image Processing Using MATLAB, Prentice-Hall, 2004
% $Revision: 1.4 $ $Date: 2003/10/26 22:35:37 $
f = double(f);
% If S is a scalar, obtain the seed image.
if numel(S) == 1
SI = f == S;
S1 = S;
else
% S is an array. Eliminate duplicate, connected seed locations
% to reduce the number of loop executions in the following
% sections of code.
SI = bwmorph(S, 'shrink', Inf);
J = find(SI);
S1 = f(J); % Array of seed values.
end
TI = false(size(f));
for K = 1:length(S1)
seedvalue = S1(K);
S = abs(f - seedvalue) <= T;
TI = TI | S;
end
% Use function imreconstruct with SI as the marker image to
% obtain the regions corresponding to each seed in S. Function
% bwlabel assigns a different integer to each connected region.
[g, NR] = bwlabel(imreconstruct(SI, TI));
该例子是通过全局阈值法来设置种子点，然后在进行生长的结果：
我们刚刚讨论过从一组种子点来生长区域的过程。另一种方法是再把图像细分为一组任意且互不相连的区域,然后在试图满足基础公式节规定的条件下合并或分离这些区域。
令R代表整个图像区域并选择一个谓词P。分割R的一种方法是把它再细分为更小的象限区域，以便对任何区域Ri都有P(Ri) = TRUE。我们从整个区域开始。若P(Ri)=FALSE，则把图像分成象限。若对任何一个象限来说Р都是FALSE，则再把象限细分为子象限，依次类推。这种特别的分离技术有一种方便的表示方法——四叉树，即一棵树，树中的每一个节点都恰好有四个后代，如图10.16所示(对应于四叉树的节点有时称为四叉区域或四叉图像)。注意,树的根对应于整个图像，每个节点对应于一个再分为四个后代节点的节点分支。在这种情况下，只有R4被进一步细分了。
若仅使用了分离,则最终的部分通常包含具有相同性质的邻近区域。这种缺点可通过合并以及分离来弥补。若要满足10.4.1节中的约束条件，则要求只合并邻近区域，组合的像素满足谓词Р。换言之，仅当P(Rj∪Rk)= TRUE时，两个邻近区域Rj和Rk才会合并。
即过程为：
IPT中实现四叉树分解的函数是qtdecomp。其语法为：s = qtdecomp (f, split_test, parameters)
其中，f是输人图像，s是包含四叉树结构的稀疏矩阵。若s ( k , m)非零，则(k, m)是分解中的一个左上角块，且该块的大小是s (k, m)。函数split_test(见下例中的函数splitmerge)用来确定一个区域是否被分离，parameters(用逗号分开)是函数split_test所要求的附加参数这种机理和前面所讨论的函数coltfilt相似。
为了在四叉树分解中得到实际的四叉区域像素值，我们使用函数qtgetblk，其语法为[ vals, r, c] - qtgetblk (f, s, m)
其中，vals是一个数组，它包含f的四叉树分解中大小为m x m的块的值，s是由qtdecomp返回的稀疏矩阵。参数r和c是包含左上角块的行坐标和列坐标的向量。
我们通过编写一个基本的分离和合并M函数来说明函数qtdecomp的用法,该函数使用前面讨论的简化,若其中的两个区域均满足谓词,则这两个区域将被合并。我们所调用的函数splitmerge的语法为：g splitmerge (f, mindim, predicate)。
其中,f是输入图像,g是输出图像，其中的每一个连接区域都用不同的整数来标注。参数mindim定义分解中所允许的最小块;该参数必须是2的正整数次幂。
函数predicate是一个用户定义的函数，它必须包括在MATLAB路径中。其语法为flag =spredicate (region)
若region中的像素满足函数中由代码所定义的谓词，则函数就必须被写成返回true（逻辑1);否则，flag 的值就必须为false（逻辑0)。
function g = splitmerge(f, mindim, fun)
%SPLITMERGE Segment an image using a split-and-merge algorithm.
% G = SPLITMERGE(F, MINDIM, @PREDICATE) segments image F by using a
% split-and-merge approach based on quadtree decomposition. MINDIM
% (a positive integer power of 2) specifies the minimum dimension
% of the quadtree regions (subimages) allowed. If necessary, the
% program pads the input image with zeros to the nearest square
% size that is an integer power of 2. This guarantees that the
% algorithm used in the quadtree decomposition will be able to
% split the image down to blocks of size 1-by-1. The result is
% cropped back to the original size of the input image. In the
% output, G, each connected region is labeled with a different
% integer.
%
% Note that in the function call we use @PREDICATE for the value of
% fun. PREDICATE is a function in the MATLAB path, provided by the
% user. Its syntax is
%
% FLAG = PREDICATE(REGION) which must return TRUE if the pixels
% in REGION satisfy the predicate defined by the code in the
% function; otherwise, the value of FLAG must be FALSE.
%
% The following simple example of function PREDICATE is used in
% Example 10.9 of the book. It sets FLAG to TRUE if the
% intensities of the pixels in REGION have a standard deviation
% that exceeds 10, and their mean intensity is between 0 and 125.
% Otherwise FLAG is set to false.
%
% function flag = predicate(region)
% sd = std2(region);
% m = mean2(region);
% flag = (sd > 10) & (m > 0) & (m < 125);
% 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/10/26 22:36:01 $
% Pad image with zeros to guarantee that function qtdecomp will
% split regions down to size 1-by-1.
Q = 2^nextpow2(max(size(f)));
[M, N] = size(f);
f = padarray(f, [Q - M, Q - N], 'post');
%Perform splitting first.
S = qtdecomp(f, @split_test, mindim, fun);
% Now merge by looking at each quadregion and setting all its
% elements to 1 if the block satisfies the predicate.
% Get the size of the largest block. Use full because S is sparse.
Lmax = full(max(S(:)));
% Set the output image initially to all zeros. The MARKER array is
% used later to establish connectivity.
g = zeros(size(f));
MARKER = zeros(size(f));
% Begin the merging stage.
for K = 1:Lmax
[vals, r, c] = qtgetblk(f, S, K);
if ~isempty(vals)
% Check the predicate for each of the regions
% of size K-by-K with coordinates given by vectors
% r and c.
for I = 1:length(r)
xlow = r(I); ylow = c(I);
xhigh = xlow + K - 1; yhigh = ylow + K - 1;
region = f(xlow:xhigh, ylow:yhigh);
flag = feval(fun, region);
if flag
g(xlow:xhigh, ylow:yhigh) = 1;
MARKER(xlow, ylow) = 1;
end
end
end
end
% Finally, obtain each connected region and label it with a
% different integer value using function bwlabel.
g = bwlabel(imreconstruct(MARKER, g));
% Crop and exit
g = g(1:M, 1:N);
%-------------------------------------------------------------------%
function v = split_test(B, mindim, fun)
% THIS FUNCTION IS PART OF FUNCTION SPLIT-MERGE. IT DETERMINES
% WHETHER QUADREGIONS ARE SPLIT. The function returns in v
% logical 1s (TRUE) for the blocks that should be split and
% logical 0s (FALSE) for those that should not.
% Quadregion B, passed by qtdecomp, is the current decomposition of
% the image into k blocks of size m-by-m.
% k is the number of regions in B at this point in the procedure.
k = size(B, 3);
% Perform the split test on each block. If the predicate function
% (fun) returns TRUE, the region is split, so we set the appropriate
% element of v to TRUE. Else, the appropriate element of v is set to
% FALSE.
v(1:k) = false;
for I = 1:k
quadregion = B(:, :, I);
if size(quadregion, 1) <= mindim
v(I) = false;
continue
end
flag = feval(fun, quadregion);
if flag
v(I) = true;
end
end
用该函数实现的效果如下：
理解分水岭变换要求我们把灰度级图像看做是一个拓扑表面,其中f (x,y)的值是被解释为高度。例如,我们可以把图10.18(a)所示的简单图像看成是图10.18(b)所示的三维表面。若雨水降落在该表面上,则雨水将流向标注为汇水盆地的两个区域中。若雨水恰好降落在标注的分水岭脊线上,则雨水流向两个汇水盆地的概率是相同的。分水岭变换会在灰度级图像中找到汇水盆地和脊线。在求解图像问题方面，关键概念是将初始图像变成另一幅图像，在变换后的图像中,汇水盆地就是我们想要识别的对象或区域。
在分割中,与分水岭变换配合使用的常用工具是距离变换。二值图像的距离变换是一个相对简单的概念:它是每一个像素到最近非零值像素的距离。图10.19示例了距离变换。图10.19(a)显示了一个小的二值图像矩阵。图10.19(b)显示了相应的距离变换。注意，值为1的像素的距离变换为0。距离变换可以使用IPT 函数bwdist加以计算，该函数的调用语法为D=bwdist(f)，效果为这样：
这里举出一个分水岭分割二值图像的例子：
f = imread('disk.tif');
se = strel('disk', 4);
f = imopen(f, se);
subplot(2,3,1);
imshow(f);
g = im2bw(f, graythresh(f));
subplot(2,3,2);
imshow(g);
subplot(2,3,3);
gc = ~g;
imshow(gc);
subplot(2,3,4);
D = bwdist(gc);
imshow(D);
subplot(2,3,5);
L = watershed(-D);
imshow(L);
subplot(2,3,6);
w = L == 0;
g2 = g & ~w;
imshow(g2);
可能是图的原因，并没有书中的效果：
在为分割使用分水岭变换之间,通常要使用梯度幅度来预处理图像。梯度幅度图像在沿对象的边缘处有较高的像素值，而在其他地方则有较低的像素值。理想情况下，分水岭变换会在沿对象边缘处产生分水岭脊线。下例说明了这一概念。
但并没有例图，所以这个例子并不好实现，思路即为先通过形态学梯度的方法求出梯度，然后再和上面的例子一样，求出分水岭山脊线即可。
分水岭变换直接用于梯度图像时,噪声和梯度的其他局部不规则性常常会导致过分割。其导致的问题可能会非常严重,以至于产生不可用的结果。按照现在的思路,这将意味着具有大量的分割区域。解决该问题的一种实际方法是把加人一个预处理阶段,以将其他知识带到分割过程中,从而限制允许的区域数目。
用于控制过分割的一种方法基于标记符的概念。标记符是一个属于一幅图像的连接分量。我们希望有一个内部标记符集合(处在每一个感兴趣对象的内部)和一个外部标记符集合(包含在背景中)。这些标记符然后使用下例中描述的过程来修改梯度图像。用于计算内部和外部标记符的方法有许多，其中包括前面描述的线性滤波、非线性滤波及形态学处理。对于特定的应用,我们选择的方法高度依赖于与应用相关的图像的特性。
效果：
代码：
f = imread('bubble.tif');
fd = double(f);
h = fspecial('sobel');
t = imfilter(fd, h, 'replicate');
g = sqrt(t .^2 + t .^2);
L = watershed(g);
wr = L == 0;
rm = imregionalmin(g);
im = imextendedmin(f, 2);
Lim = watershed(bwdist(im));
fim = f;
fim(im) = 175;
em = Lim == 0;
g2 = imimposemin(g, im | em);
L2 = watershed(g2);
f2 = f;
f2(L2 == 0) = 255;
subplot(2,4,1);
imshow(f);
title('f');
subplot(2,4,2);
imshow(wr);
title('wr');
subplot(2,4,3);
imshow(rm);
title('rm');
subplot(2,4,4);
imshow(fim);
title('fim');
subplot(2,4,5);
imshow(em);
title('em');
subplot(2,4,6);
imshow(L2);
title('L2');
subplot(2,4,7);
imshow(f2);
title('f2');
解释：
rm = imregionalmin(g);该函数计算图像中大量局部最小区域的位置，从而了解为什么函数watershed会产生如此多的小汇水盆地。
图rm中显示的多数局部最小区域位置非常浅，并且表示了与我们的分割问题不相关的细节。为消除这些无关的小区域，我们使用IPT函数imextendedmin，该函数可计算图像中的“低点”集合，即比周围点更深的点的集合(通过某个高度阈值)。
其中，f是一幅灰度级图像,h是高度阙值，im是一幅二值图像，该二值图像的前景像素标记了深局部最小区域的位置。fim = f; fim(im) = 175; em = Lim == 0;
最后两行将在原图像上以灰色气泡的形式叠加扩展的局部最小区域位置，如图fm所示。我们看到，得到的气泡确实很合理地标记了我们想要分割的对象。
图em显显示了二值图像em中的分水岭脊线。因为这些脊线位于由im标记的暗气泡间的中间位置，所以它们应该是很好的外部标记符。
给出内部和外部标记符后，然后就可使用它们来修改梯度图像，方法是使用称为强制最小( minima imposition ）的过程。强制最小技术修改了灰度级图像，以便局部最小区域仅出现在标记的位置。其他像素值将按需要“上推”，以便删除其他的局部最小区域。IPT函数imimposemin可实现这种技术。
g2 = imimposemin(f, mask);其中，f是一幅灰度级图像，mask是一幅二值图像，该二值图像的前景像素标记出了输出图像mp中局部最小区域的期望位置。通过在内部和外部标记符的位置覆盖局部最小区域，我们可修改梯度图像:g2 = imimposemin(g, im | em);。
然后再对g2 进行分水岭操作，即可得到最后的f2.