数字图像处理形态学的空域与频域实现

王玉文 胡顺波

摘要：数字图像处理课程的形态学主要是从几何学的角度分析图像，不同于课程前面章节从空域和频域角度分析图像，因而难以被学生理解和接受。文章探讨形态学、空间域和频率域运算三个主要课程内容之间的联系，阐述形态学基本操作的空间域卷积和频率域乘积的实现步骤。仿真实验结果表明，形态学操作可以通过空间域卷积和频率域乘积实现，它们是密切联系在一起的。

关键词：数字图像处理;形态学;腐蚀;膨胀;卷积;频率

中图分类号：TP393        文献标识码：A

文章编号：1009-3044（2022）18-0074-03

开放科学（资源服务）标识码（OSID）：

1 引言

新工科自2016年提出以来，形成了三个纲领性文件：“复旦共识”“天大行动”和“北京指南”，揭示了新工科的特质：引领性、通宽性、前瞻性、交叉性、开放性和实践性[1-2]。新工科是主动应对新一轮科技革命与产业变革的战略发展新需求，对应的是新兴产业，如人工智能、智能制造、机器人、云计算等，当然也包括传统工科专业的升级改造[3]。数字图像处理作为计算机、通信、信息、人工智能等专业方向的重要课程，在新工科的智能处理方向建设中起重要作用。融会贯通数字图像处理的理论框架，串联形态学、空间域处理和频率域处理的关系[4-5]，有助于数字图像处理课程的教学工作，促进学生对这门课程的学习和掌握。

形态学处理在数字图像处理课程中，常单独作为一个章节。主要从几何学的角度，定义数学形态学的一些算法，执行处理任务[6]。从理论上讲，形态学好像和其他章节关系不大，相对难以被学生理解和接受。但这部分内容可以和整个数字图像处理的理论框架（空间域处理和频率域处理）密切相关。空（间）域处理的核心是卷积运算，频（率）域处理的核心是乘积操作，卷积定理是他们之间的桥梁和纽带。本文将先说明形态学的基本概念：膨胀和腐蚀，然后再探讨其“空域卷积，频域乘积”实现，最后给出具体实现步骤和实验结果。

2 形态学中腐蚀和膨胀的基本概念

假设[A]表示二值图像，[B]表示一个二值模板，[C]是图像[A]中数值为1（目标）的像素集合，结构元[D]是模板[B]中数值为1的像素集合，则[B]对[A]的腐蚀操作定义如下：

[A?B=z]                              （1）

z表示像素。[Dz]表示当模板[B]的中心平移到z位置处，[A]在[B]为1位置处（结构元[D]）的像素集合。在像素z处，如果满足条件[Dz?C]，即[A]在结构元[D]对应的所有[D]个位置处的像素值都为1，则腐蚀图像在像素z位置处取值1;否则取值0。[D]表示集合[D]中的元素总数目。

同理，非反转结构元版本[7]的膨胀操作的数学定义如下：

[A⊕B=z]                                       （2）

当模板[B]的中心平移到z位置处，[A]在结构元对应位置处的所有[D]个像素中，只要有一个像素值为1，则膨胀图像在像素z位置处取值为1;否则为0。

3 腐蚀和膨胀的卷积实现

设二值图像[A]的大小为[M×N]，模板[B]的大小为[m×n]，[m=2a+1]，[n=2b+1]。则他们之间（线性）卷积的定义如下：

[A?Bx，y=s=-aat=-bbBs，tA（x-s，y-t）]       （3）

上述卷积操作可以用八个字来概括：反转，平移，相（点）乘，累加。首先，绕模板[B]的中心原点反转;然后，把反转后的模板平移到[x，y]处;对应[m×n]个位置处的像素值相乘;最后，把[m×n]个乘积项累加到一起，得到卷积在[x，y]处的输出值。本文中，“线性卷积”简称“卷积”。

从卷积、腐蚀和膨胀的定义可见，共同点是都需要平移操作。其不同点体现在如下两个方面：1）卷积运算对模板[B]先进行反转，再进行平移，而腐蚀和膨胀运算对模板[B]直接进行平移操作，不需要反转;2）腐蚀和膨胀使用了集合关系进行定义，而卷积使用了代数运算进行定义。所以，针对这两个不同点，对卷积运算和运算后的输出图像做如下设计：1）[B=imrotateB，180?，center]，即模板[B]和[B]互為反转关系，这将保证卷积运算与腐蚀、膨胀运算在进行平移操作时候，使用相同的模板;2）对卷积的结果[A?B]进行二值化处理，得到二值图像[E]：

[Ex，y=1A?Bx，y≥T0A?Bx，y

如果阈值[T]取数值1，则二值图像[E]等于膨胀图像，[E=A⊕B]。如果阈值[T]取数值[D]（相对于图像[A]，无卷积出界点的区域），或在图像[A]四个边界附近取[D∩A]（相对于图像[A]，存在卷积出界点的区域），则二值图像[E]等于腐蚀图像，[E=A?B]。这样，通过卷积运算和后续处理，实现了形态学的腐蚀和膨胀。

综上，通过空间域卷积运算，实现腐蚀和膨胀的主要步骤如下：

1）输入二值图像[A]和模板[B];

2）对模板[B]进行反转，得到反转后的模板[B];

3）对[A]和[B]进行卷积操作，得到卷积结果[Gp]=[conv2A，B];

4）对[Gp]进行裁剪，保留大小和[A]相同的中间部分，记为[G];

5）类似公式（4），对[G]进行二值化处理，得到二值图像[E]。当阈值分别取1和[D]（当模板出界，[D∩A]）时，对应[E]分别是膨胀图像和腐蚀图像。

裁剪操作目的是得到和输入二值图像[A]同样大小的图像。如果省略该步骤，则得到全尺寸为[M+m-1×N+n-1]的输出图像。步骤4）和5）交换顺序，可得到相同结果。

4 腐蚀和膨胀的频率域实现

二维离散卷积定理说明周期卷积可以通过频率域乘积来实现。而在线性移不变系统中，系统的真实输出，是输入图像与单位冲击响应的线性卷積。所以如何使图像和模板的周期卷积等于线性卷积，是频率域乘积求线性卷积的关键。但是，由于离散傅里叶变换具有周期性，导致在使用卷积定理时，二维数字图像在空间域不再认定为单张图片，而是认定为周期性出现的。进而，导致了使用‘频域乘积，再进行离散傅里叶反变换计算线性卷积的时候，会发生空间域的周期缠绕错误。为了避免这种空间域的周期缠绕错误，需要对输入图像和模板，通过0填充进行周期延拓。当填充足够多的0时，空间域图像就无法缠绕在一起。此时，通过频域乘积计算的周期卷积，其有效部分就等于线性卷积。所以，我们可以通过频率域乘积，来实现空间域的线性卷积，进而实现腐蚀和膨胀。

图1给出线性卷积、周期卷积和补0后周期卷积的计算对比示意图。[4×4]黑色数字部分表示原始图像，[3×3]红色数字部分表示对称模板。图中，模板的中心都对应原始图像的右下角值为7的位置，此时其三种卷积的输出值分别如下：1）线性卷积，1?6+1?4+1?5+1?7+（1?0）?5=22; 2）周期卷积，1?6+1?4+1?5 +1?7 +1?8 +1?2 + 1?3+ 1?1=36;3）补0后周期卷积，1?6+1?4+1?5+1?7+（1?0）?5=22。由图1可见，由于离散傅里叶变换引入了周期性，导致了周期卷积的缠绕错误，在图像边界处的周期卷积和线性卷积结果不一致。通过在原图像的行和列进行补0（紫色）操作，消除了周期卷积的缠绕错误，从而使补0后的周期卷积和线性卷积结果相同。

综上所述，通过频率域乘积，实现腐蚀和膨胀运算的主要步骤如下：

1）输入二值图像[A]和模板[B];

2）对模板[B]进行反转，得到反转后的模板[B];

3）对[A]和[B]进行0填充，得到[Ap]和[Bp]：

[Apx，y=Ax，y0≤x≤M-1 和 0≤y≤N-10M≤x≤P 或 N≤y≤Q ]，

[Bpx，y=Bx，y0≤x≤m-1 和 0≤y≤n-10m≤x≤P 或 n≤y≤Q]。

其中，[P≥M+m-1，  Q≥N+n-1];

4）分别对[Ap]和[Bp]进行二维离散傅里叶变换，得到[FAp=fft2Ap]，[FBp=fft2Bp];

5）频域乘积，对[FAp]和[FBp]进行对应元素相乘，即哈达玛乘积，得到频域滤波结果[FG p]=[FAp?FBp];

6）对[FGp]进行二维离散傅里叶反变换，得到[Gp]=[ ifft2FGp]。

7）对[Gp]进行裁剪，保留大小和输入图像[A]相同的中间部分，记为[G]。

8）同上节步骤5）。

步骤7）和8）可以交换顺序，就是先进行二值化，再进行裁剪，也能得到同样的结果。在实际的频率域实现过程中，考虑离散傅里叶正变换和反变换过程中的数值离散化误差和计算精度误差，阈值可以取稍微小一点的数值，[T=T-ε]，[ε]是很小的正数。

5 实验

使用Matlab语言编程，分别对具有简单边缘的仿真图像和具有复杂边缘的人物图像进行实验，实验结果分别见图2和图3。实验代码在Github上公开[8]。

在图2中，（a）是原始彩色图像，大小为[385×512]，含有六种简单图形目标;（b）是对称模板[B]，大小为[25×25]，每个元素值都为1，[D=625];（c）是由（a）生成的二值图像[A]，即先通过rgb2gray函数把彩色图像（a）转换为灰度图像，然后通过二值化处理获得;（d）是腐蚀图像;（e）是膨胀图像。

因为对于定义法、空域卷积法和频域乘积法，都得到了完全一样的实验结果[8]，所以腐蚀图像（d）和膨胀图像（e）仅画出一次，代表三类实现方案的结果。当使用定义法的时候，腐蚀运算、膨胀运算分别用imerode、imdilate函数;当使用空域卷积法的时候，主要用带有‘same参数的conv2函数来实现;对于频域乘积法，主要使用了fft2、矩阵点乘运算、ifft2函数来编程实现。

对实际人物的二值图像，我们也做了类似实验，实验结果如图3所示。在图3中，（a）是原始二值图像[A]，大小为[269×221]，具有相对复杂的轮廓;（b）是对称模板[B]，大小为[11×11]，为1的像素占据中间菱形区域，剩余四角部分的像素值为0，[D=61]。因为[A]下边界附近存在白色目标区域，所以，当使用空域卷积法和频域乘积法的进行腐蚀操作时候，存在卷积出界点的腐蚀区域，阈值采用前文中的出界点策略计算方法，实现代码见链接 [8]。（c）是腐蚀图像，（d）是膨胀图像。同样地，三类方法得到相同的腐蚀图像和膨胀图像[8]。

以上两个实验表明，无论对腐蚀运算，还是膨胀运算，三类实现方案都可以得到同样的结果，所以这三类实现方案可以相互转换。形态学、空间域卷积和频率域乘积三者之间是密切联系在一起的。这样，把数字图像处理课程中凌乱的三块知识点统一起来。此外，卷积是现在流行的深度学习中卷积神经网络的基本模块，把这三个概念进行外延，就和深度学习和人工智能相联系。

6 结论

本文对形态学腐蚀和膨胀运算，给出了三类实现方案：定义法、空域卷积法、频域乘积法。这三类实现方案可以从理论上和算法上相互转换，并通过实验，验证了三类方案能得到同样结果。因而，从概念和算法中理顺数字图像处理中三大板块（空间域处理、频率域处理、形态学处理）的区别和联系，洞悉腐蚀运算和膨胀运算的空域卷积和频域乘积的实现细节。通过腐蚀运算或膨胀运算，可以来定义其他更复杂的形态学运算，例如开、闭运算、击中-击不中变换等。这些复杂形态学运算也可以通过空域卷积实现，或频域乘积来实现。

参考文献：

[1] 周端明，沈燕培.习近平科技创新重要论述指引新工科建设的方向[J].高等工程教育研究，2021（4）：1-7.

[2] 赵继，谢寅波.新工科建设与工程教育创新[J].高等工程教育研究，2017（5）：13-17，41.

[3] 徐艳，朱孔伟.“数字信号处理”课程思政教学的融入点探索[J].教育教学论坛，2021（7）：117-120.

[4] 杨森泉，丁凡，陈景华，刘建军，胡松喜.新工科背景下数字图像处理课程实训改革与实践[J].物联网技术，2021，11（7）：126-127+130.

[5] 唐晓芬，郭辉.以创新能力培养为目标的数字图像处理课程教学改革研究[J].计算机教育，2021（8）：112-115，120.

[6] Rafael C. Gonzalez， Richard E. Woods. 数字图像处理（第四版）[M]. 阮秋琦，阮宇智，译. 北京：电子工业出版社， 2020：454-467.

[7] Dilate image. Matlab imdilate[EB/OL]. https：//www.mathworks.com/help/images/ref/ imdilate.html.

[8] Hu Shunbo. dilation_erosion_conv_freq[EB/OL]. https：//github.com/hushunbo/dilation_ erosion_conv_freq.

【通联编辑：王力】