数字图像处理MATLAB学习笔记(八)
Morphological Image Processing形态学图像处理
1. Preliminaries
1.1 Some Basic Concepts from Set Theory
简单的、基本的不阐述。
满足一定特殊条件的像素坐标集合B记为 B = { w ∣ condition } B=\{w|\text{condition}\} B={w∣condition}
所有像素坐标的集合不属于集合A,即A的补集,记为 A C = { w ∣ w ∉ A } A^C=\{w|w\notin A\} AC={w∣w∈/A}
morphological operation通常还需要两个算子,它们特别针对元素是坐标的集合。
集合B的反射 B ^ \hat{B} B^定义为 B ^ = { w ∣ w = − b , b ∈ B } \hat{B}=\{w|w=-b,b\in B\} B^={w∣w=−b,b∈B}
点 z = ( z 1 , z 2 ) z=(z_1,z_2) z=(z1,z2)对集合A的平移表示为 ( A ) z (A)_z (A)z,定义为 ( A ) z = { c ∣ c = a + z , a ∈ A } (A)_z=\{c|c=a+z,a\in A\} (A)z={c∣c=a+z,a∈A}
1.2 Binary Images, Sets and Logical Operators 二值图像、集合、逻辑算子
形态学中,二值图像被看作是前景(单值)像素的集合,集合的元素属于 Z 2 Z^2 Z2。
在MATLAB中,逻辑运算对应关系如下表所示:
Set Operation | MATLAB Expression for Binary Images | Name |
---|---|---|
A ∩ B {A\cap B} A∩B | A & B | AND |
A ∪ B {A\cup B} A∪B | A | B | OR |
A C {A^C} AC | ~A | NOT |
A − B {A-B} A−B | A & ~B | DIFFERENCE |
2. Dilation and Erosion膨胀和腐蚀
2.1 Dilation膨胀
Dilation是使图像中的目标“生长”或“变粗”的操作。
图像A被structuring element结构元B膨胀,定义为 A ⊕ B = { z ∣ ( B ^ z ∩ A ≠ ϕ ) } A\oplus B=\{z|(\hat{B}_z\cap A \ne\phi)\} A⊕B={z∣(B^z∩A=ϕ)}。A被B膨胀是由所有结构元的原点位置组成的集合,反射并平移后的B至少与A的一个元素重叠。膨胀满足结合律和交换律。
也就是说,图像A根据结构元B及其膨胀原点位置进行膨胀操作。
在MATLAB中,使用function imdilate:
D = imdilate(A, B)
这里假定A和B都是二值的,那么结构元B中的原点会被自动计算,依据的公式为
floor((size(B) + 1)/2)
这一操作得到包含结构元中点坐标的二维向量。如果需要用原点不在中心的结构元来处理,方法是用0填充B,以使原始的中心移到希望的位置。
Example An application of dilation:
我们通过如下结构元对图像进行膨胀处理,结构元B为:
0
1
0
1
1
1
0
1
0
\begin{matrix} 0&1&0\\ 1&\color{red}\bold{1}&1\\ 0&1&0 \end{matrix}
010111010
原点为中心位置的1
>> A = imread('Fig0906(a).tif');
>> B = [0 1 0; 1 1 1; 0 1 0];
>> D = imdilate(A, B);
>> subplot(121);imshow(A);axis off
>> subplot(122);imshow(D);axis off

2.2 Structuring Element Decomposition结构元的分解
假定结构元B可以描述为结构元 B 1 B_1 B1和 B 2 B_2 B2的膨胀: B = B 1 ⊕ B 2 B=B_1\oplus B_2 B=B1⊕B2
因膨胀满足结合律,我们可以得到 A ⊕ B = A ⊕ ( B 1 ⊕ B 2 ) = ( A ⊕ B 1 ) ⊕ B 2 A\oplus B=A\oplus(B_1\oplus B_2)=(A\oplus B_1)\oplus B_2 A⊕B=A⊕(B1⊕B2)=(A⊕B1)⊕B2。结合律很重要,因为计算膨胀时需要的时间正比于结构元中非零像素的个数。
换句话说,一个如下的5x5数组的结构元可以分解为如下两个结构元的膨胀:
[
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
]
=
[
1
1
1
1
1
]
⊕
[
1
1
1
1
1
]
\begin{bmatrix} 1&1&1&1&1\\ 1&1&1&1&1\\ 1&1&\color{red}\bold{1}&1&1\\ 1&1&1&1&1\\ 1&1&1&1&1 \end{bmatrix}= \begin{bmatrix} 1&1&\color{red}\bold{1}&1&1 \end{bmatrix}\oplus \begin{bmatrix} 1\\ 1\\ \color{red}\bold{1}\\ 1\\ 1 \end{bmatrix}
⎣⎢⎢⎢⎢⎡1111111111111111111111111⎦⎥⎥⎥⎥⎤=[11111]⊕⎣⎢⎢⎢⎢⎡11111⎦⎥⎥⎥⎥⎤
这样进行行列分解,总元素数目由原来的25减少为10,即先用行结构元膨胀,再用列结构元膨胀。这样计算,时间和空间上的开销都会相对节省一些,具有重大意义。
2.3 Function strel
MATLAB中使用function strel构建不同的结构元:
Se = strel(shape, parameters)
其中,shape是期望返回的结构元名称,parameters是生成对应结构元需要传入的参数。
Syntax Form | Description |
---|---|
se = strel(‘diamond’, R) | 构造扁平的菱形结构元。规定从结构元原点到菱形最远点的距离为R |
se = strel(‘disk’, R) | 构造半径为R的扁平的圆形结构元。 |
se = strel(‘line’, LEN, DEG) | 构造扁平的直线型结构元。规定从水平轴开始以逆时针方向度量,LEN代表长度,DEG表示线倾斜的角度。 |
se = strel(‘octagon’, R) | 创建扁平的八边形结构元。规定沿水平轴和垂直轴度量,结构元原点到八边形边的距离为R。R必须为3的非负倍数。 |
se = strel(‘pair’, OFFSET) | 构造包括两个成员的扁平的结构元。一个成员在原始位置,第二个成员在向量OFFSET规定的位置。OFFSET必须是整形的二元向量。 |
se = strel(‘periodicline’, P, V) | 构造包含2*P+1个成员的扁平结构元。V是一个二元向量,包括行列值均为整数的位移。结构元的一个成员在原始位置,其他成员都分别在1*V、-1*V、2*V、-2*V、…、P*V和-P*V处。 |
se = strel(‘rectangle’, MN) | 构造扁平的矩形结构元,MN规定了矩形的大小。MN必须为二元的非负整形矢量,MN中第一个元素是结构元中的行数,第二个元素是结构元中的列数。 |
se = strel(‘square’, W) | 构造方形结构元,其中宽度为W。W必须为非负整形变量。 |
se = strel(‘arbitrary’, NHOOD);se = strel(NHOOD) | 构造人系形状的结构元。NHOOD是规定了形状的0/1矩阵。 |
Example Structuring element decomposition using function strel:
考虑用strel 函数构造菱形结构元:
>> se = strel('diamond', 5)
se =
strel is a diamond shaped structuring element with properties:
Neighborhood: [11×11 logical]
Dimensionality: 2
>> se.Neighborhood
ans =
11×11 logical array
0 0 0 0 0 1 0 0 0 0 0
0 0 0 0 1 1 1 0 0 0 0
0 0 0 1 1 1 1 1 0 0 0
0 0 1 1 1 1 1 1 1 0 0
0 1 1 1 1 1 1 1 1 1 0
1 1 1 1 1 1 1 1 1 1 1
0 1 1 1 1 1 1 1 1 1 0
0 0 1 1 1 1 1 1 1 0 0
0 0 0 1 1 1 1 1 0 0 0
0 0 0 0 1 1 1 0 0 0 0
0 0 0 0 0 1 0 0 0 0 0
2.4 Erosion 腐蚀
腐蚀“收缩”或“细化”二值图像中的物体。像膨胀一样,收缩的方法和程度由结构元控制。
图像A被结构元B腐蚀,定义为 A ⊖ B = { z ∣ ( B ) z ⊆ A } A\ominus B=\{z|(B)_z\subseteq A\} A⊖B={z∣(B)z⊆A}。在进行结构元检查时,需要检查是否完全匹配,完全匹配的前景图像才会在对应的结构元原点处赋值,其余区域均为0(针对二值图像)。
这个公式表明A被B腐蚀是包含在A中的B由z平移的所有点z的集合。因为B包含在A中,所以不会和A的背景元素有交集。所以腐蚀定义公式可以为: A ⊖ B = { z ∣ ( B ) z ∩ A = ϕ } A\ominus B=\{z|(B)_z\cap A=\phi\} A⊖B={z∣(B)z∩A=ϕ}。在这里,A被B腐蚀是所有结构元的原点位置不与A背景重叠的B的部分。
Example An illustration of erosion:
腐蚀用工具箱函数imerode执行,语法与imdilate语法相同。
>> A = imread('Fig0908(a).tif');
>> se = strel('disk', 10);
>> E10 = imerode(A, se);
>> subplot(221);imshow(A);axis off
>> subplot(222);imshow(E10);axis off
>> se = strel('disk', 5);
>> E5 = imerode(A, se);
>> subplot(223);imshow(E5);axis off
>> E20 = imerode(A, strel('disk', 20));
>> subplot(224);imshow(E20);axis off

3. Combining Dilation and Erosion
3.1 Opening and Closing
A被B形态学开操作表示为 A ∘ B A\circ B A∘B,定义为A被B腐蚀,然后再用B膨胀腐蚀结果: A ∘ B = ( A ⊖ B ) ⊕ B A\circ B=(A\ominus B)\oplus B A∘B=(A⊖B)⊕B。与开操作等价的数学表达式为: A ∘ B = ∪ { ( B ) z ∣ ( B ) z ⊆ A } A\circ B=\cup \{(B)_z|(B)_z\subseteq A\} A∘B=∪{(B)z∣(B)z⊆A}。(这里的 ∪ { ∙ } \cup\{ \bullet \} ∪{∙}表示所有集合的并集)
A被B形态学闭操作表示为 A ∙ B A\bullet B A∙B,指先膨胀再腐蚀: A ∙ B = ( A ⊕ B ) ⊖ B A\bull B=(A\oplus B)\ominus B A∙B=(A⊕B)⊖B。

对应的MATLAB中使用function imopen和function imclose实现:
C = imopen(A, B)
C = imclose(A, B)
Example Working with functions imopen and imclose:
>> f = imread('Fig0910(a).tif');
>> se = strel('square', 40);
>> fo = imopen(f, se);
>> imshow(fo)
>> fc = imclose(f, se);
>> figure, imshow(fc);
>> foc = imclose(fo, se);
>> figure, imshow(foc);

开操作和闭操作的顺序可用于噪声消除。
>> f = imread('Fig0911(a).tif');
>> se = strel('square', 3);
>> fo = imopen(f, se);
>> foc = imclose(fo, se);
>> subplot(131);imshow(f);title('Orginal image');axis off
>> subplot(132);imshow(fo);title('Opening');axis off
>> subplot(133);imshow(foc);title('Opening followed by closing');axis off
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-diWi3XSh-1647138958625)(/Users/wanghe/Library/Application Support/typora-user-images/image-20211213233746753.png)]
通过先开后闭的操作,可以消除一些噪音。但是我们发现会仍有代价,即部分指纹处引入的有一些缺口。
3.2 The Hit-or-Miss Transformation 击中或击不中变换
通常,能够匹配一幅图像中像素的特定结构是很有用的,比如孤立的前景像素或是线段的端点像素。击中或击不中变换对这类应用非常实用。
A被B击中或击不中变换表示为 A ⊗ B A\otimes B A⊗B。其中, B B B是结构元对,即 B = ( B 1 , B 2 ) B=(B_1,B_2) B=(B1,B2)。击中或击不中变换用两个结构元定义为: A ⊗ B = ( A ⊖ B 1 ) ∩ ( A C ⊖ B 2 ) A\otimes B=(A\ominus B_1)\cap(A^C\ominus B_2) A⊗B=(A⊖B1)∩(AC⊖B2)。
在MATLAB中使用function bwhitmiss实现:
C = bwhitmiss(A, B1, B2)
Example Using function bwhitmiss:
>> f = imread('Fig0913(a).tif');
>> B1 = strel([0 0 0; 0 1 1; 0 1 0]);
>> B2 = strel([1 1 1; 1 0 0; 1 0 0]);
>> g = bwhitmiss(f, B1, B2);
>> subplot(121);imshow(f);axis off;title('original image');
>> subplot(122);imshow(g);axis off;title('result of applying the hit-or-miss transformation');

bwhitmiss的替代语法可以把B1和B2组合成间隔矩阵。只要B1等于1或-1,B2等于1,间隔矩阵就等于1。对于不关心的像素,间隔矩阵等于0。
对于B1和B2的间隔矩阵为
>> interval = [-1 -1 -1; -1 1 1; -1 1 0]
interval =
-1 -1 -1
-1 1 1
-1 1 0
使用这个矩阵,function bwhitmiss中的B1和B2就可以替换为interval进行计算。
3.3 Using Lookup Tables 运用查询表
当击中或击不中结构元较小的时候,计算击中或击不中变换比较快速的方法是利用查找表(LUT)。
LUT是预先计算出每个可能邻域结构的输出像素,然后把这些值存入查找表以备后面使用。
例如:在二值图像中,有 2 9 = 512 2^9=512 29=512个不同的3x3像素值的结构。为了实际使用查找表,我们需要为每个可能的结构分配唯一的索引。在3x3情况下,可以创建如下的矩阵: [ 1 8 64 2 16 128 4 32 256 ] \begin{bmatrix}1&8&64\\2&16&128\\4&32&256\end{bmatrix} ⎣⎡1248163264128256⎦⎤,通过结构元素与该矩阵没个元素相乘累加得到结果。结果是在区域 [ 0 , 511 ] [0,511] [0,511]内对每个3x3邻域结构赋予唯一值。例如,邻域为 [ 1 1 0 1 0 1 1 0 1 ] \begin{bmatrix}1&1&0\\1&0&1\\1&0&1\end{bmatrix} ⎣⎡111100011⎦⎤时,对应的赋值应该为 1 × 1 + 2 × 1 + 4 × 1 + 8 × 1 + 16 × 0 + 32 × 0 + 64 × 0 + 128 × 1 + 256 × 1 = 399 1\times1+2\times1+4\times1+8\times1+16\times0+32\times0+64\times0+128\times1+256\times1=399 1×1+2×1+4×1+8×1+16×0+32×0+64×0+128×1+256×1=399。
在MATLAB中使用function makelut和function applylut。其中makelut来构造查询表,applylut用这个查询表处理二值图像。
为此,我们编写function endpoints,用makelut和applylut在二值图像中检测端点:
function g = endpoints(f)
%ENDPOINTS Computes end points of a binary image.
% G = ENDPOINTS (F) computes the end points of the binary image F
% and returns them in the binary image G.
persistent lut
if isempty(lut)
lut = makelut(@endpoint_fcn, 3);
end
g = applylut(f, lut);
% ----------------------------------------------------------------------
function is_end_point = endpoint_fcn(nhood)
%Determines if a pixel is an end point.
% IS_END_POINT = ENDPOINT_FCN(NHOOD) accepts a 3-by-3 binary
% neighborhood, NHOOD, and returns a 1 if the center element is an end
% point; otherwise it returns a 0.
interval1 = [0 1 0; -1 1 -1; -1 -1 -1];
interval2 = [1 -1 -1; -1 1 -1; -1 -1 -1];
% Use bwhitmiss to see if the input neighborhood matches either interval1
% or interval2, or any of their 90-degree rotations.
for k = 1:4
% rot90(A, K) rotates the matrix A by 90 degrees k times
C = bwhitmiss(nhood, rot90(interval1, k));
D = bwhitmiss(nhood, rot90(interval2, k));
if (C(2,2) == 1) || (D(2,2) == 1)
% Pixel neighborhood matches one of the end-point configurations,
% so return true
is_end_point = true;
return
end
end
% Pixel neighborhood did not match any of the end_point configurations, so
% return false.
is_end_point = false;
建立了名为lut的变量,并表明这个变量是永久的(persistent)。在函数间调用时,MATLAB会记住这个永久变量的值。在第一次调用该函数时,lut变量自动创建,之后便永久保存,不需要重新计算。
Example Playing Conway’s Game of Life using binary images and look-up-table-based computations:
Conway的遗传规则描述了如何从现在的一代计算下一代(或是下一个二值图像)。
- 每个拥有2个或3个相邻前景像素的前景像素可以产生下一代。
- 每个拥有0个、1个或至少4个相邻前景像素的前景像素会“死亡”(变成背景像素),因为已被“隔离”或是“数目过剩”。
- 每个背景像素与3个真正的前景邻域相连就可以“出生”像素,并且变成前景像素。
在计算描述下一代的下一个二值图像的过程中,所有出生和死亡同时发生。针对Conway游戏规则的方法如下:
function out = conwaylaws(nhood)
%CONWAYLAWS Applies Conway's genetic laws to a single pixel.
% OUT = CONWAYLAWS (NHOOD) applies Conway's genetic laws
% to a single pixel and its 3-by-3 neighborhood, NHOOD.
num_neighbors = sum(nhood(:)) - nhood(2,2);
if nhood(2,2) == 1
if num_neighbors <= 1
out = 0; % Pixel dies from isolation
elseif num_neighbors >= 1
out = 0; % Pixel dies from overpopulation
else
out = 1; % Pixel survives.
end
else
if num_neighbors == 3
out = 1; % Birth pixel.
else
out = 0; % Pixel remains empty.
end
end
使用具有函数句柄的makelut调用为conwayslaws建立查找表:
>> lut = makelut(@conwaylaws, 3);
各种初始图像已被设计用来演示在相继代中Conway规则的效果,例如,考虑“咧嘴笑猫”的初始图像:
>> bw1 = [0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 1 0 0 1 0 0 0
0 0 0 1 1 1 1 0 0 0
0 0 1 0 0 0 0 1 0 0
0 0 1 0 1 1 0 1 0 0
0 0 1 0 0 0 0 1 0 0
0 0 0 1 1 1 1 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0];
显示第三代结果:
>> subplot(131), imshow(bw1, 'InitialMagnification', 'fit'), title('Generation 1')
>> bw2 = applylut(bw1, lut);
>> subplot(132), imshow(bw2, 'InitialMagnification', 'fit'), title('Generation 2')
>> bw3 = applylut(bw2, lut);
>> subplot(133), imshow(bw3, 'InitialMagnification', 'fit'), title('Generation 3')

3.4 Function bwmorph
function bwmorph执行许多以膨胀、腐蚀和查找表运算相结合为基础的形态学操作。语法如下:
g = bwmorph(f, operation, n)
operation是指定所希望运算的字符串,n是制定重复次数。
Operation | Description |
---|---|
bothat | 用3x3结构元进行“底帽”操作;对其他结构元使用imbothat操作 |
bridge | 连接被单像素缝隙分割的像素 |
clean | 去掉孤立的前景像素 |
close | 用由1构成的3x3结构元进行闭操作;对其他结构元使用imclose操作 |
diag | 填充围绕对角线相连的前景像素 |
dilate | 用3x3结构元进行膨胀;对其他结构元使用imdilate操作 |
erode | 用3x3结构元进行腐蚀;对其他结构元使用imerode操作 |
fill | 填充单个像素的“洞”(被前景像素环绕的背景像索);使用imfl1操作来填充更大的洞 |
hbreak | 去掉H型相连的前景像素 |
majority | 如果 N 8 ( p ) N_8(p) N8(p)中至少有5个像索为前景像素,产生前景像素P;否则产生背景像素P |
open | 用由1组成的3x3结构元进行开操作;对其他结构使用imopen操作. |
remove | 去掉“内部”像素(没有背景邻域的前景像素) |
shrink | 将物体收缩成没有洞的点;将物体收缩成带洞的环形 |
skel | 骨骼化图像 |
spur | 去掉“毛刺”像素 |
thicken | 粗化物体,并且不加入不连贯的1 |
thin | 细化没有洞的物体到最低限度相连的笔画;将物体细化成带洞的环形 |
tophat | 用由1 组成的3x3结构元进行“顶帽”操作;对其他结构元使用imtophat操作 |
4. Labeling Connected Components 标记连通分量
连通分量:即可以根据某种特定的规则或连接方式,实现某一点到某一点的通路,这两个点就构成该规则下的连通分量。
在MATLAB中使用function bwlabel实现:
[L, num] = bwlabel(f, conn)
其中conn是希望的连接方式,即我们所说的规则。当conn省略时,默认为8连接,即中心p可以向四周八个方向移动。输出L叫做标记矩阵,num则给出找到的连通分量的总数。
Example Computing and displaying the center of mass of connected components:
我们计算并显示该图片中每个字母和图像的连通分量的质心,使用的是8连通分量:
>> f = imread('Fig0917(a).tif');
>> [L, n] = bwlabel(f);
>> [r, c] = find(L == 3);
>> rbar = mean(r);
>> cbar = mean(c);
>> imshow(f)
>> hold on
>> for k = 1:n
[r, c] = find(L == k);
rbar = mean(r);
cbar = mean(c);
plot(cbar, rbar, 'Marker', 'o', 'MarkerEdgeColor', 'k', 'MarkerFaceColor', 'k','MarkerSize', 15);
plot(cbar, rbar, 'Marker', '*', 'MarkerEdgeColor', 'w');
end

5. Morphological Reconstruction 形态学重建
重建是一种形态学变换,包含了两幅图像和一个结构元。一幅图像是标记,是变换的开始点;另一幅图像是模版,用于约束变换过程;结构元来定义连通性。
如果G为模版,F为标记,F必须是G的子集。从F重建G记做 R G ( F ) R_G(F) RG(F),迭代过程如下:
- 将标记图像F初始化 h 1 h_1 h1
- 建立结构元:B=ones(3)
- 重复 h k + 1 = ( h k ⊕ B ) ∩ G h_{k+1}=(h_k\oplus B)\cap G hk+1=(hk⊕B)∩G,直到 h k + 1 = h k h_{k+1}=h_k hk+1=hk
在MATLAB中可以使用function imreconstruct实现:
out = imreconstruct(marker, mask)
5.1 Opening by Reconstruction 通过重建进行开操作
在形态学开操作中,腐蚀典型地去除小的物体,且随后的膨胀趋向于恢复保留的物体形状。但是这种恢复的精确度取决于形状和结构元之间的相似性。
Example Opening by reconstruction:
通过进行开闭操作,我们将具有长垂直线的字母保留
>> f = imread('Fig0922(a).tif');
>> fe = imerode(f, ones(51, 1));
>> fo = imopen(f, ones(51, 1));
>> fobr = imreconstruct(fe, f);
>> subplot(221);imshow(f);axis off;title('Original Image')
>> subplot(222);imshow(fe);axis off;title('Image eroded with a vertical line')
>> subplot(223);imshow(fo);axis off;title('Image opened with a vertical line')
>> subplot(224);imshow(fobr);axis off;title('Image opened by reconstruction with a vertical line')
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-YOoexdcW-1647138958626)(/Users/wanghe/Library/Application Support/typora-user-images/image-20211215134542632.png)]
5.2 Filling Holes 填充孔洞
我们令
I
I
I表示二值图像,
F
F
F表示标记图像,除了图像边缘外,其余位置都为0,边缘部分值为
1
−
I
1-I
1−I:
F
(
x
,
y
)
=
{
1
−
I
(
x
,
y
)
if
(
x
,
y
)
is on the border of
I
0
otherwise
F(x,y)=\begin{cases}1-I(x,y)&\text{if}\space(x,y)\space\text{is on the border of} \space I\\0 &\text{otherwise}\end{cases}
F(x,y)={1−I(x,y)0if (x,y) is on the border of Iotherwise
H
=
[
R
I
C
(
F
)
]
C
H=[R_{I^C}(F)]^C
H=[RIC(F)]C是一幅相当于
I
I
I的填充了所有孔洞的二值图像。在MATLAB中使用function imfill:
g = imfill(f, 'holes')
5.3 Clearing Border Objects 清除外界物体
重建的另一-种应用是清除图像中与边缘相接触的物体。同样,关键任务仍然是选择合适的标记来达到希望的效果。假定定义标记图像F为:
F
(
x
,
y
)
=
{
I
(
x
,
y
)
if (x,y) is on the border of
I
0
otherwise
F(x,y)=\begin{cases}I(x,y) & \text{if (x,y) is on the border of }I\\0&\text{otherwise}\end{cases}
F(x,y)={I(x,y)0if (x,y) is on the border of Iotherwise
其中
I
I
I是原始图像。用
I
I
I作为模版图像,重建
H
=
R
I
(
F
)
H=R_I(F)
H=RI(F)得到一幅图像
H
H
H,即原图中与外界接触的物体。那么
I
−
H
I-H
I−H即原图中不与边缘接触的物体。
在MATLAB中,我们用function imclearborder产生 I − H I-H I−H,这样可以去除鼻周围物体更亮且与图像边界相连接的结构:
g = imclearborder(f, conn)
6. Gray-Scale Morphology 灰度级形态学
6.1 Dilation and Erosion
我们用结构元b对灰度图像f进行dilation operation,表示为
f
⊕
b
f\oplus b
f⊕b,定义如下:
(
f
⊕
b
)
(
x
,
y
)
=
m
a
x
{
f
(
x
−
x
′
,
y
−
y
′
)
+
b
(
x
′
,
y
′
)
∣
(
x
′
,
y
′
)
∈
D
b
}
(f\oplus b)(x,y) = max\{f(x-x',y-y')+b(x',y')|(x',y')\in D_b\}
(f⊕b)(x,y)=max{f(x−x′,y−y′)+b(x′,y′)∣(x′,y′)∈Db}
其中,
D
b
D_b
Db为结构元
b
b
b的域,
f
(
x
,
y
)
f(x,y)
f(x,y)假设在
f
f
f域之外为
−
∞
-\infin
−∞。执行过程类似于卷积,只不过这里是使用的mask是二值的,且需要计算最大值。
一般情况下,灰度膨胀通常用flat structuring element平坦结构元完成,其中 b b b的值在 D b D_b Db域中所有坐标均为0,即 b ( x ′ , y ′ ) = 0 b(x',y')=0 b(x′,y′)=0。此时,灰度膨胀运算公式就可以简化为: ( f ⊕ b ) ( x , y ) = m a x { f ( x − x ′ , y − y ′ ) ∣ ( x ′ , y ′ ) ∈ D b } (f\oplus b)(x,y) = max\{f(x-x',y-y')|(x',y')\in D_b\} (f⊕b)(x,y)=max{f(x−x′,y−y′)∣(x′,y′)∈Db}。
因此,平的灰度膨胀时局部最大值算子,其中,这个最大值是由 D b D_b Db中1值元素的空间形状决定的一系列邻域值得到。
对于non-flat structuring element非平坦结构元,可以通过两个矩阵使用function strel构造。两个矩阵分别为:
- 由结构元的域指定的0/1矩阵
- 由高度值制定的矩阵
同理,用结构元b对灰度图像f进行erosion operation,表示为
f
⊖
b
f\ominus b
f⊖b,定义如下:
(
f
⊖
b
)
(
x
,
y
)
=
m
i
n
{
f
(
x
−
x
′
,
y
−
y
′
)
−
b
(
x
′
,
y
′
)
∣
(
x
′
,
y
′
)
∈
D
b
}
(f\ominus b)(x,y) = min\{f(x-x',y-y')-b(x',y')|(x',y')\in D_b\}
(f⊖b)(x,y)=min{f(x−x′,y−y′)−b(x′,y′)∣(x′,y′)∈Db}
其中,
D
b
D_b
Db为结构元
b
b
b的域,
f
(
x
,
y
)
f(x,y)
f(x,y)假设在
f
f
f域之外为
−
∞
-\infin
−∞。其余信息与上述的gray-scale dilation operation一致。同理,灰度腐蚀也常用平坦结构元完成,其中
b
b
b的值在
D
b
D_b
Db域中所有坐标均为0,即
b
(
x
′
,
y
′
)
=
0
b(x',y')=0
b(x′,y′)=0。此时,灰度膨胀运算公式就可以简化为:
(
f
⊖
b
)
(
x
,
y
)
=
m
i
n
{
f
(
x
−
x
′
,
y
−
y
′
)
∣
(
x
′
,
y
′
)
∈
D
b
}
(f\ominus b)(x,y) = min\{f(x-x',y-y')|(x',y')\in D_b\}
(f⊖b)(x,y)=min{f(x−x′,y−y′)∣(x′,y′)∈Db}。
膨胀和腐蚀可以相互结合以获得各种效果。例如,从图像的膨胀结果中减去腐蚀过的图像可产生“Morphological gradient形态学梯度”,这是图像中局部灰度变化的一种度量。
>> f = imread('Fig0923(a).tif');
>> subplot(221);imshow(f);title('Original image');axis off
>> se = strel('square', 3);
>> gd = imdilate(f, se);
>> subplot(222);imshow(gd);title('Dilated image');axis off
>> ge = imerode(f, se);
>> subplot(223);imshow(ge);title('Eroded image');axis off
>> morph_grad = gd - ge;
>> subplot(224);imshow(morph_grad);title('Morphological gradient');axis off

6.2 Opening and Closing 开操作和闭操作
在灰度图像中,开操作和闭操作的表达式与二值图像拥有相同的形式。用结构元b对灰度图f进行开操作时,表达式为 f ∘ b = ( f ⊖ b ) ⊕ b f\circ b=(f\ominus b)\oplus b f∘b=(f⊖b)⊕b;进行闭操作的表达式为 f ∙ b = ( f ⊕ b ) ⊖ b f\bullet b = (f\oplus b)\ominus b f∙b=(f⊕b)⊖b。
灰度开/闭操作可以用几何解释。假设图像 f ( x , y ) f(x,y) f(x,y)用三维表面表示,即图像的Intensity为xy平面上的高度值。 b b b对 f f f的开操作可以在几何上解释为推动结构元 b b b使之沿表面 f f f的下沿平移,并移过整个 f f f的域。开操作的结果是寻找结构元滑过的 f f f下沿上所能达到的最高点,即波谷。
同理,闭操作在几何上解释为推动结构元 b b b使之沿表面 f f f的上沿平移,并移过整个 f f f的域。闭操作的结果是寻找结构元滑过的 f f f上沿上所能达到的最低点,即波峰。
也就是说,开操作是消除了小的亮点细节,闭操作是消除了小的黑点细节。
Example Morphological smoothing using opening and closing:
我们利用开操作和闭操作对图像中的亮暗区域进行平滑处理:
>> f = imread('Fig0925(a).tif');
>> se = strel('disk', 5);
>> fo = imopen(f, se);
>> foc = imclose(fo, se);
>> fc = imclose(f, se);
>> fco = imopen(fc, se);
>> for k = 2:5
se = strel('disk', k);
fasf = imclose(imopen(fasf, se), se);
imshow(fasf);
end

6.2.1 Tophat and Bottomhat Transformation
顶帽变换和底帽变换:
tophat transformation顶帽变换:灰度图像 f f f的顶帽变换定义为 f f f减去其开运算: T h a t ( f ) = f − ( f ∘ b ) T_{hat}(f)=f-(f\circ b) That(f)=f−(f∘b)。
bottomhat transformation底帽变换:灰度图像 f f f的顶帽变换定义为 f f f的闭运算减去 f f f: T h a t ( f ) = ( f ∙ b ) − f T_{hat}(f)=(f\bullet b)-f That(f)=(f∙b)−f。
对应的MATLAB中,使用function imtophat和function imbothat进行计算。
这些变换的主要应用之一是,在开运算或闭运算中用一个结构元从图像中删除目标,而不是拟合将被删除的目标。然后,差运算得到一幅仅保留已删除分量的图像。顶帽变换用于暗背景上的亮目标,而底帽变换则用于亮背景上的暗目标。因此,我们通常将这两个变换称为白顶帽变换和黑底帽变换。
顶帽变换的一个重要用途是校正不均匀光照的影响。
6.2.2 Granulometry 粒度测定
粒度测定是指确定图像中颗粒的大小分布。由于颗粒几乎无法整齐地分开,因此采用逐个识别颗粒的方法来计算颗粒的数量非常困难。形态学可间接估计颗粒的大小分布,而不需要识别和测量各个颗粒。
这种方法很简单。对于比背景亮且形状规则的颗粒,这种方法是用逐渐增大的结构元对图像进行开运算。
基本思想是,某个特殊大小的开运算应会对包含类似大小颗粒的输入图像的那些区域产生最大影响。对于开运算得到的每幅图像,我们计算像素值之和。这个和值称为表面区域,它随结构元大小的增大而减小,因为开运算会减小图像中的亮特征。这一过程得到的是一个一维阵列,阵列中的每个元素都是矩阵中对应该位置的结构元大小的开运算的像素之和。
6.3 Reconstruction
形态学中灰度图像的重建仍基于如下迭代过程:
如果G为模版,F为标记,F必须是G的子集。从F重建G记做 R G ( F ) R_G(F) RG(F),迭代过程如下:
- 将标记图像F初始化 h 1 h_1 h1
- 建立结构元:B=ones(3)
- 重复 h k + 1 = ( h k ⊕ B ) ∩ G h_{k+1}=(h_k\oplus B)\cap G hk+1=(hk⊕B)∩G,直到 h k + 1 = h k h_{k+1}=h_k hk+1=hk
完