形态学图像处理
说明：“形态学”通常指的是生物学的某个分支，常用来处理动物和植物的形状和结构。现在，我们也将这个词用于图像处理中，作为提取图像分量的一种工具，这些分量在表示和描述区域形状、形态学滤波、细化、粗化等是非常有用的。
（一）基础知识
（1）集合论中的基本概念
如果（x,y）是来自Z2的整数，f是分配给每个不同坐标对（x,y）的亮度值的映射，那么函数f(x,y)被称为数字图像。令A是Z2中的集合，A中的元素是坐标（x,y）处的像素。
若w=(x,y)属于 A中的元素，则表示为：
若w 不属于 A中的元素，表示为：
若满足一定条件的像素坐标集合B，表示为：
若像素坐标的集合均不属于集合A，则表示为Ac：
若集合C同属于集合A和集合B，交集表示为：
集合A和集合B的并集C，表示为：
集合A和集合B的差集（属于A但不属于B的所有像素坐标），表示为：
集合B的反射表示为：
点z=(z1,z2)对集合A的平移表示为：
（2）二值图像、集合及逻辑算子
二值图像是指在图像中，灰度等级只有两种，通俗来说，就是图像中的任何像素不是0就是1，再无其他过渡的灰度值。
编写实验代码：
f = imread('C:\Users\Public\Pictures\Sample Pictures\Fig0627(b).tif');
g = imread('C:\Users\Public\Pictures\Sample Pictures\Fig0627(c).tif');
subplot(231),imshow(f),title('f');
subplot(232),imshow(g),title('g');
a = f&g;
subplot(233),imshow(a,[]),title('f和g的交集');
b = f|g;
subplot(234),imshow(b),title('f和g的并集');
c = f&~g;
subplot(235),imshow(c),title('f和g的差集');
d = ~f;
subplot(236),imshow(d),title('f的补集');
（二）膨胀与腐蚀
（1）膨胀
顾名思义，膨胀就是使图像中的目标“变粗”的操作。这种特殊的方式和变粗程度由一种被称为结构元的形状来控制的。A被B膨胀，表示为A⊕B，集合操作表示为：
其中，A是图像，B是结构元，结构元通常比图像小得多。膨胀满足结合律和交换律。
在matlab中使用函数reflect来计算结构元的反射。
使用函数imdilate进行膨胀，语法为：
D = imdilate(A,B) %A是图像，B是结构元
编写实验代码：
A = imread('C:\Users\Public\Pictures\Sample Pictures\dingdangmao.jpg');
B = [0 1 0;1 1 1;0 1 0];
subplot(121),imshow(A),title('原图像')
D = imdilate(A,B);
subplot(122),imshow(D),title('执行膨胀后')
B = [1 1 1;0 0 0;1 1 1];更换结构元后的图像
实验分析： 当设置结构元B后，对A图像进行膨胀，变换后的图像是将内部进行膨胀，因而看到叮当猫的黑色轮廓变细了、眼睛变小了、胡须也变细了。是由于实验所设计的结构元B是内部为1，而边缘部分为0所造成的。当设置不同的结构元，操作膨胀的结果是不同的，可以根据需要调整结构元。
（2）结构元的分解
结构元B可以描述为结构元B1和B2的膨胀：B=B1⊕B2。由于膨胀满足结合律，所以A⊕B=A⊕(B1⊕B2)=(A⊕B1)⊕B2，这样B就可以分解为B1和B2两个结构元。
进行结构元的分解是使膨胀的速度增加，减少其他不必要的开销。 比如说：
这样的矩阵操作膨胀需要对25个元素进行膨胀，若将其分解成如下形式，
这样的操作速度就会从膨胀25个元素变成操作10个元素即可。
（3）strel函数
strel函数是用于构造各种形状和大小的结构元，语法为：
se = strel(shape,parameters) %shape是用于指定希望形状的字符串，parameters是形状信息的参数
用函数strel编写实验代码：
se = strel('diamond',5) %构造菱形结构元，水平和垂直均扩展5个像素
decomp = getsequence(se); %提取并检查分解中单独的结构元
whos
into decomp:
decomp(1).Neighborhood
decomp(2).Neighborhood
decomp(3).Neighborhood
decomp(4).Neighborhood
（4）腐蚀
腐蚀是膨胀的反操作，是将图像中的目标进行“细化”的操作，收缩的方法和程度也由结构元控制。A被B腐蚀表示为AΘB，集合操作表示为：
在matlab工具中，运用函数imerode进行腐蚀操作，编写代码如下：
A = imread('C:\Users\Public\Pictures\Sample Pictures\dingdangmao.jpg');
subplot(221),imshow(A),title('原始图像')
se = strel('disk',10);
E10 = imerode(A,se);
subplot(222),imshow(E10),title('半径为10进行腐蚀')
se = strel('disk',5);
E5 = imerode(A,se);
subplot(223),imshow(E5),title('半径为5进行腐蚀')
se = strel('disk',20);
E20 = imerode(A,se);
subplot(224),imshow(E20),title('半径为20进行腐蚀')
实验分析： 腐蚀相当于膨胀的反操作，将内部图像进行腐蚀，黑色边缘就变得越来越粗，当我们使用半径不同的结构元时，可以得到不一样的腐蚀结果。腐蚀可以用于弥补小裂缝的。假设把叮当猫的眼睛当成孔洞，用半径为20的圆形结构元进行腐蚀，结果显示可以将腐蚀操作可以眼睛填充，作用弥合裂缝孔洞等。
（三）膨胀与腐蚀的结合
用膨胀与腐蚀有三种组合：开操作、闭操作、击中或击不中变换。
（1）开操作和闭操作
开操作是指A被B腐蚀，然后再用B膨胀的结果，开操作表示为：
开操作主要用于平滑目标的轮廓，断开了细的连接部分，去掉了细的突出。
闭操作指的是A被B膨胀，再进行腐蚀的操作，表示为：
闭操作主要用于连接窄的断裂并填满细长的“沟壑”,填满比结构元小的孔洞。
在matlab中用函数imopen表示开操作，函数imclose表示闭操作：
C = imopen(A,B) %开操作，B是0和1的矩阵
C = imclose(A,B) %闭操作，B是0和1的矩阵
编写实验代码：
f = imread('C:\Users\Public\Pictures\Sample Pictures\dingdangmao.jpg');
subplot(221),imshow(f),title('原始图像')
se = strel('square',10);
fo = imopen(f,se);
subplot(222),imshow(fo),title('开操作');
fc = imclose(f,se);
subplot(223),imshow(fc),title('闭操作');
foc = imclose(fo,se);
subplot(224),imshow(foc),title('先开操作，后闭操作');
se = strel('square',10); 修改方形结构元为10后得到的结果
se = strel('square',3); 修改方形结构元为3后得到的结果
实验分析： 开操作可以用在去除噪声的实验中，这并没有改变图像的位置与基本轮廓，当实验只关心物体的位置和个数的时候，腐蚀的计算速度比单独过滤噪声更具有优势。另外说明，开闭运算是对偶的，对于某图像多次应用开运算或者闭运算，和只进行一次运算的效果相同。正如上面说明的那样，开操作和闭操作中的腐蚀和膨胀与结构元的状态紧密相关，当选取的结构元半径比较小的时候，基本看不出变换，反之选取的半径过大，会将原始图像变成另一种模样，因此选择膨胀和腐蚀的结构元时要选择合适的结构元的结构状态。
（2）击中和或击不中变换
击中击不中变换是形态学形状检测的基本工具，是用来查找像素局部模式的形态学运算符。‘局部’是指‘结构元素’的大小。
A被B击中或击不中变换表示为：
其中，B是结构元对，表示如下：
击中或击不中变换的输出图像由所有在B1中匹配的像素（击中）和未在B2中匹配的像素（击不中）组成。
使用函数bwhitmiss实现，语法为：
C = bwhitmiss(A,B1,B2)
编写代码如下：
f = imread('C:\Users\Public\Pictures\Sample Pictures\Fig0913(a).tif');
subplot(121),imshow(f),title('原始图像')
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(122),imshow(g,[]),title('击中或击不中变换')
interval = [-1 -1 -1;-1 1 1;-1 1 0] %间隔矩阵
执行击中或击不中变换的图像，仔细一看，是将原图像中符合匹配的点提取出来，显示在上面，由于图像比较暗，所以不是很明显。击中或击不中变换就是当符合的点提取，用来查找局部内容的形态学操作。
（3）运用查询表
当击中或击不中结构元较小的时候，计算击中或击不中变换比较快速的方法是利用查询表（LUT）。这种方法预先计算出每个可能邻域结构的输出像素值，然后把这些值存入查询表中以备后面使用。
比如说，有3×3矩阵系数和像素值结构：
（矩阵系数）
（像素值结构）
然后将对应相乘并相加起来，得到最后的结果值：1（1）+2（1）+4（1）+8（1）+16（0）+32（0）+64（0）+128（1）+256（1）=399
工具箱中提供了两个函数makelut和函数applylut用于查询表，函数makelut用于构造查询表，函数applylut用这个查询表处理二值图像。我们编写函数endpoints.m用于makelut和applylut在二值图像中检测端点，使用如下：persistent lut，建立名为lut的变量，并标明这个变量是永久的。MATLAB在函数间调用时会记住这个变量的值。函数endpoints第一次被调用时，变量lut自动初始化为空矩阵([ ])。当lut为空时，调用函数makelut用它发送句柄给函数endpoint_fc。之后，用函数applylut查询表寻找端点时，查询表被存入永久变量lut中，以后调用就不用重新计算了。
函数endpoints的M-文件：
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.
% Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
% Digital Image Processing Using MATLAB, Prentice-Hall, 2004
% $Revision: 1.3 $ $Date: 2003/04/22 01:18:18 $
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.
is_end_point = nhood(2,2) & (sum(nhood(:)) == 2);
编写实验代码：
persistent lut
lut = makelut(@conwaylaws,3)
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];
imshow(bw1,'InitialMagnification','fit'),title('Generation 1')
bw2 = applylut(bw1,lut);
figure,imshow(bw2,'InitialMagnification','fit'),title('Generation 2')
bw3 = applylut(bw2,lut);
figure,imshow(bw3,'InitialMagnification','fit'),title('Generation 3')
（4）bwmorph函数
函数bwmorph h用于执行许多以膨胀、腐蚀和查找表运算相结合为基础的形态学操作，语法是：
g = bwmorph(f,operation,n) %operation是指定所希望运算的字符串，n是指定重复次数的正整数
以上操作，主要关注的是细化和骨骼化。细化就是指将图像中的目标形状减少为单个像素宽度的“笔画”。
f = imread('C:\Users\Public\Pictures\Sample Pictures\Fig0911(a).tif');
g1 = bwmorph(f,'thin',1);
g2 = bwmorph(f,'thin',2);
subplot(231),imshow(f),title('原始图像');
subplot(232),imshow(g1),title('一次细化');
subplot(233),imshow(g2),title('二次细化');
ginf = bwmorph(f,'thin',Inf);
subplot(234),imshow(ginf),title('细化到稳定状态')
fs = bwmorph(f,'skel',Inf);
subplot(235),imshow(fs),title('骨骼化图像')
fs = bwmorph(fs,'spur',5); %与骨骼化效果大概一致
subplot(236),imshow(fs),title('除去“毛刺”图像')
实验分析： 进行去除“毛刺”操作，得到的图像与上面骨骼化后的图像基本一致。虽然所运用的算法操作不同，但是结果会可能基本一致。这个去除“毛刺”的操作被称为裁剪，可以使用endpoints来获得，其原因是进行反复确认并去除端点。通过5次去除端点的迭代，变成了骨骼化的结果。
实验分析： 用指纹图像再执行一次算法，说明算法的稳定性。将第五幅图和第六幅图进行对比，可以明显的看出，去除“毛刺”操作与骨骼化是大概一致，但附近的噪声点发现在第六幅图中得到了滤除。一次细化后比原始图像的噪声点少一些，而且指纹也变细了，细化到稳定状态时，指纹基本上只剩下一个像素的轮廓，说明已经达到了最细的状态。
（四）标记连通分量
连通分量概念：
无向图的极大连通子图称为的连通分量( Connected Component)。任何连通图的连通分量只有一个，即是其自身，非连通的无向图有多个连通分量（来自百度百科）。
邻域定义：
坐标为(x,y)的像素p有两个水平和两个垂直相邻像素，它们的坐标分别为(x+1,y)、(x-1,y)、(x,y+1)、(x,y-1)，这4个相邻像素的集合表示为
。p的4个对角线相邻像素的坐标分别是(x-1,y-1)、(x+1,y-1)、(x-1,y++1)、(x+1,y+1)，它们称为
。上面两种邻域的并集是p的8相邻像素，表示为
。
邻接定义：
如果
，则q和p就被称为4邻接的；同样，q若属于8邻域、D邻域内，则q和p就分别是8邻接和D邻接。
如果在像素p和q之间存在一条4连接通路，就称这两个前景像素为4连接的；同理，如果在像素p和q之间存在一条8连接通路，就称这两个前景像素为8连接的。连通分量用通路来定义，而通路的定义按照所用邻接的类型而定。通俗来说，连通分量的性质决定我们所选择的邻接方式，常见的为4邻接和8邻接，同时邻接方式还决定图像中连通分量的数目。
在matlab中用函数bwlabel计算二值图像中所有的连通分量，基本语法为：
[L,num] = belabel(f,conn) %L是标记矩阵，num是找出的连通分量的总数，f是二值图像，conn是连接方式（默认为8）
编写实验代码：
f = imread('C:\Users\Public\Pictures\Sample Pictures\Fig0219(a).tif');
[L,n] = bwlabel(f);
imshow(f),title('原始图像');
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',10);
plot(cbar,rbar,'Marker','*','MarkerEdgeColor','w'),title('相应连通分量上的质心显示为*');
end
此操作为在连通分量的质心上标记*，以此来说明连通分量的位置及个数。
（五）形态学重建
重建是一种形态学变换，包括两幅图像和一个结构元。一幅图像是标记，是变换的开始点；另一幅图是模板，用于约束变换过程。结构元时用来定义连通性，默认是8连接。
基本迭代过程是（模板是G，标记是F）：
- 将标记图像F初始化为h1；
- 建立结构元：B = onese(3);
- 重复hk+1=(hk⊕B)∩G，直到hk+1=hk。
以上是理论的方法，但实际会使用更快的计算方法，在matlab中使用函数imreconstruct用于快速混合重建算法，编写为：
out = imreconstruct(marker,mask)
（1）通过重建进行开操作
在形态学开操作中，腐蚀典型地去除小的物体，且随后的膨胀趋向于恢复保留的物体形状。但这个恢复的精确度取决于形状和结构元之间的相似性，说明并不是所有情况都是适用的。现在，用一种比较好的方式来解决这个问题，那就是通过重建进行开操作能准确地恢复腐蚀之后的物体形状。
其集合表示为：
编写实验代码：
f = imread('C:\Users\Public\Pictures\Sample Pictures\Fig0917(a).tif');
fe = imerode(f,ones(51,1));
fo = imopen(f,ones(51,1));
fobr = imreconstruct(fe,f);
subplot(221),imshow(f),title('原始图像');
subplot(222),imshow(fe),title('用竖线腐蚀后图像');
subplot(223),imshow(fo),title('用竖线进行开操作图像');
subplot(224),imshow(fobr),title('用竖线重建后图像');
（2）填充孔洞
孔洞填充： 先选取十字形结构作为结构元；接着生成一个由0组成的阵列X，大小与二值图像A相同，随机选取A中一个灰度值为1的点，将X中相应位置的灰度值也设置为1；之后根据填充的公式多次迭代，直到A不再变化，就停止迭代；最后对A和X求并集则得到了填充的图像。
数学表达式表示为：假定F是标记图像，边缘部分为1-I(x,y)，
然后进行
H是一幅相当于I的填充了所有孔洞的二值图像。
在MATLAB工具箱中可以使用函数imfill将所有的孔洞进行填充：
g = imfill(f,'holes')
观察图像发现，这个操作将原来是封闭部分的图像填充起来，完成孔洞填充的应用效果。但实验结果表明，要是孔洞不是完全封闭的，则该操作并不是将其填充起来。孔洞填充的前提是，图像是完全封闭的，才可以进行该操作。
（3）清除边界物体
重建的又一应用是清除图像中与边缘相接触的物体。
假定定义标记图像F为：
其中I是原始图像，然后以I为模板图像，重建
得到一幅图像H，清除边缘物体的图像。
在MATLAB工具箱中可以使用函数imclearborder将图像边缘接触的物体清除：
g = imclearborder(f,4) %conn的值默认是8，也可以选择4
（六）灰度级形态学
（1）膨胀和腐蚀
用结构元b对灰度图像f的灰度进行膨胀，表示为：
可以将灰度膨胀简化为如下形式，这是使用平的结构元进行计算的：
在实验中，非平的结构元通过两个矩阵使用strel函数进行构造。这两个矩阵包括由结构元的域指定的0和1的矩阵和由高度值指定的矩阵，编写代码：
b = strel([1 1 1],[1 2 1])
然而平坦的结构元也是通过函数strel构造的，编写实验代码：
f = imread('C:\Users\Public\Pictures\Sample Pictures\Fig0906(a).tif');
subplot(121),imshow(f),title('原始图像');
se = strel('square',3);
gd = imdilate(f,se);
subplot(122),imshow(gd),title('膨胀后图像');
用结构元b对灰度图像f的灰度进行腐蚀，表示为：
同理，通常灰度腐蚀使用的结构元时平的结构元，简化为：
编写实验代码：
f = imread('C:\Users\Public\Pictures\Sample Pictures\Fig0917(a).tif');
subplot(121),imshow(f),title('原始图像');
ge = imerode(f,se);
subplot(122),imshow(ge),title('腐蚀后图像');
morph_grad = gd-ge; %用于计算形态学的梯度
imshow(morph_grad),title('形态学梯度')
（2）开操作和闭操作
在灰度图像中，开操作和闭操作的表达式与二值图像拥有一样的形式，表示为：
（开操作）
（闭操作）
几何解释：假设图像函数f(x,y)用三维表面表示，然后b对f的开操作可以在几何上解释为推动结构元b，使之沿表面f的下沿平移，移过整个f的域。可以用下面的图像进一步理解：
开操作和闭操作应用于平滑图像、非均匀背景的补偿、粒度的测定：
编写实验一代码：平滑图像
f = imread('C:\Users\Public\Pictures\Sample Pictures\Fig0925(a).tif');
subplot(131),imshow(f),title('原始图像');
se = strel('disk',5);
fo = imopen(f,se);
foc = imclose(fo,se);
subplot(132),imshow(fo),title('用半径5的圆盘进行开操作图像');
subplot(133),imshow(foc),title('再进行闭操作的图像');
f = imread('C:\Users\Public\Pictures\Sample Pictures\Fig0925(a).tif');
subplot(131),imshow(f),title('原始图像');
se = strel('disk',5);
fc = imclose(f,se);
fco = imopen(fc,se);
subplot(132),imshow(fc),title('用半径5的圆盘进行开操作图像');
subplot(133),imshow(fco),title('再进行闭操作的图像');
fasf = f; %交替顺序滤波
for k = 2:5
se = strel('disk',k);
fasf = imclose(imopen(fasf,se),se);
end
进行交替顺序滤波后得到的结果：对比三幅图，看出交替顺序滤波后的图像更加平滑，图像也变得明亮一些。
编写实验二代码：非均匀背景补偿
f = imread('C:\Users\Public\Pictures\Sample Pictures\Fig0926(a).tif');
subplot(131),imshow(f),title('原始图像');
se = strel('disk',10);
fo = imopen(f,se);
subplot(132),imshow(fo),title('开操作后图像');
f2 = f - fo;
f2 = imtophat(f,se);
subplot(133),imshow(f2),title('阈值处理后图像');
se = strel('disk',3);
g = f + imtophat(f,se)-imbothat(f,se); %顶帽变换减去底帽变化
imshow(g)
编写实验三代码：粒度测定
f = imread('C:\Users\Public\Pictures\Sample Pictures\Fig0926(a).tif');
imshow(f),title('原始图像')
sumpixels = zeros(1,36);
for k = 0:35
se = strel('disk',k)
fo = imopen(f,se);
sumpixels(k+1) = sum(fo(:));
end
plot(0:35,sumpixels),xlabel('k'),ylabel('Surface area')
plot(-diff(sumpixels))
xlabel('k')
ylabel('Surface area reduction')
（相对于结构元半径的表面区域）
（相对于结构元半径的表面区域的减少）
（3）重建
灰度重建和二值图像的重建是理论基本一样的，不同在于灰度图像中通常使用的结构是平坦结构元。形态学重建可以看作是对一个图像进行连续的膨胀操作，这个图像称之为标识图像marker。当连续的膨胀操作使得标识图像的轮廓匹配了另一个图像，重建过程结束，另一个图像便是掩模图像mark。
编写实验一代码：使用形态学开操作重建图像
f = imread('C:\Users\Public\Pictures\Sample Pictures\Fig0309(a).tif');
subplot(131),imshow(f),title('原始图像');
se =strel('disk',5);
fe = imerode(f,se);
subplot(132),imshow(fe),title('腐蚀过后图像');
fobr = imreconstruct(fe,f);
subplot(133),imshow(fobr),title('腐蚀后再重建的图像');
fobrc = imcomplement(fobr);
subplot(131),imshow(fobrc),title('补集的图像');
fobrce = imerode(fobrc,se);
subplot(132),imshow(fobrce),title('补集后腐蚀的图像');
fobrcbr = imcomplement(imreconstruct(fobrce,fobrc));
subplot(133),imshow(fobrcbr),title('补集后腐蚀再补集的图像');
编写实验二代码：使用重建移除复制的背景
f = imread('C:\Users\Public\Pictures\Sample Pictures\Fig0917(a).tif');
subplot(131),imshow(f),title('原始图像');
f_obr = imreconstruct(imerode(f,ones(1,71)),f); %重建的开操作
subplot(132),imshow(f_obr),title('重建后图像');
f_o = imopen(f,ones(1,71)); %标准的开操作
subplot(133),imshow(f_o),title('开操作图像');
f_thr = f - f_obr; %从原始图像中减去重建开操作
f_th = f -f_o; %从原始图像中减去标准开操作
g_obr = imreconstruct(imerode(f_thr,ones(1,11)),f_thr); %用水平线对f_thr的重建开操作
g_obrd = imdilate(g_obr,ones(1,21)); %用水平线对g_obr进行膨胀处理
f2 = imreconstruct(min(g_obrd,f_thr),f_thr); %最后重建步骤
subplot(221),imshow(f_th),title('变换后图像');
subplot(222),imshow(g_obr),title('用水平线对f_thr的重建开操作后图像');
subplot(223),imshow(g_obrd),title('用水平线对g_obr进行膨胀处理后图像');
subplot(224),imshow(f2),title('再重建后图像');
3695
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
