Digital Image Processing Basics

像素 \((x_1,y_1)\) 和 \((x_2,y_2)\) 的距离有不同的度量方法，如：
一，欧氏距离(Euclidean Distance)，可记为 \(D_E = \sqrt{(x_1 - x_2)^2 + (y_1 - y_2)^2}\)
二，城区距离(City-Block Distance)，两点间只允许横向和纵向的移动，可记为 \(D_4 = \vert x_1 - x_2 \vert + \vert y_1 - y_2 \vert\)
三，棋盘距离(Chessboard Distance)，两点间还允许对角线方向的移动（国际象棋里国王从一处移动到另一处的步数），可记为 \(D_8 = \max\{\vert x_1 - x_2 \vert, \vert y_1 - y_2 \vert\}\)

距离为 1 的两点是邻接的（采用城区距离时称为 4-邻接，采用棋盘距离时称为 8-邻接）。

从一个简单的例子开始讨论。假设一幅大小为 \(500 \times 500\) 像素的图像要放大 1.5 倍到 \(750 \times 750\) 像素。一种简单的放大方法是创建一个 \(750 \times 750\) 网格，它与原始图像有相同的间隔，然后将其收缩，使它准确地与原图像匹配。显然，收缩后的 \(750 \times 750\) 网格的像素间隔要小于原图像的像素间隔。为了对覆盖的每一个点赋以灰度值，我们在原图像中寻找最接近的像素，并把该像素的灰度赋给 \(750 \times 750\) 网格中的新像素。当完成对网格中覆盖的所有点的灰度赋值后，就把图像扩展到原来规定的大小，得到放大后的图像。

上面描述的方法称为最近邻插值(Nearest neighbor interpolation)，因为这种方法把原图像中最近邻的灰度赋给了每个新位置。

前面讨论的最近邻插值算法比较简单，但其效果不是很好，如可能导致直边缘严重失真。双线性插值是一种更实用的方法。

在介绍双线性插值前，先了解一下线性插值(Linear interpolation)。线性插值非常简单，使用连接两个已知量的直线来确定在这两个已知量之间的一个未知量的值的方法。线性插值的一个实例如图 3 所示。

双线性插值(Bilinear interpolation)是线性插值算法在二维空间的扩展。 双线性插值利用 4 个最近邻点的灰度值来估计点的灰度值，如图 4 中，已知 4 个最近邻点的灰度值分别为 91,210,162,95，需要估计点 \((20.2, 14.5)\) 处的灰度值。
双线性插值的具体步骤如下：
第 1 步，用线性插值法求得点 \((20, 14,5)\) 和点 \((21, 14.5)\) 的估计灰度值，即：
\[\begin{aligned} I_{(20, 14.5)} = \frac{15-14.5}{15-14} \cdot 91 + \frac{14.5 - 14}{15 -14} \cdot 210 = 150.5 \\ I_{(21, 14.5)} = \frac{15-14.5}{15-14} \cdot 162 + \frac{14.5-14}{15-14} \cdot 95 = 128.5 \end{aligned}\]
第 2 步，再次用线性插值法求得得 \((20.2, 14.5)\) 的估计灰度值，即：
\[I_{(20.2, 14.5)} = \frac{21-20.2}{21-20} \cdot 150.5 + \frac{20.2-20}{21-20} \cdot 128.5 = 146.1\]

双三次插值(Bicubic interpolation)通过 16 个最近邻点来估计点的灰度值，它比双线性插值(通过 4 个最近邻点来估计点的灰度值)得到的效果更好。双三次插值是商业图像编辑程序（比如 Adobe Photoshop 和 Corel Photopaint）的标准插值方法。

这里不详细介绍双三次插值的具体算法。

对于灰度级范围为 \([0, L-1]\) 的图像，其反转图像由下式给出：
\[s = L-1-r\]
其中， \(r\) 是输入灰度级， \(s\) 是输出灰度级。

这种类型的处理特别适用于增强嵌入在一幅图像的暗区域中的白色或灰色细节，特别是当黑色面积在尺寸上占主导地位时。 比如，图 6 左边图像是乳房 X 射线照片，显然分析乳房组织时使用反转图像会容易得多。

对数变换的通用形式为：
\[s = c \log(1+r)\]
其中， \(c \ge 0\) 是一个常数。r是输入灰度级，s是输出灰度级。

从图 5 中曲线的形状可知，对数变换会将输入中“范围较窄的低灰度值”映射为输出中“范围较宽的灰度值”，这样较暗的区域对比度将得到提升，从而增强图像较暗区域的细节。反对数变换的作用与此相反（它能增强图像较亮区域的细节）。不过，后文介绍的“幂律（伽马）变换”也能实现灰度级的扩展/压缩，而且“幂律（伽马）变换”更为通用。

对数变换的主要应用场景在于压缩像素值变化较大的图像的动态范围。 像素值有较大动态范围的一个典型应用是傅里叶频谱，频谱的动态范围可从 0 到 \(10^6\) ，甚至更高。图像的显示系统通常无法真实地再现如此大范围的灰度值。当作线性地缩放显示时，最亮的像素将支配该显示，频谱中的低值（恰恰是重要的）将损失掉。因而，最终结果是许多重要的灰度细节在傅里叶频谱的显示中丢失（无法察觉）了。 通过对数变换可以压缩像素值变化较大的图像的动态范围，通过这种非线性地变换，图像的一些细节将展示地更清楚。 图 7 中左边是傅里叶频谱，右边是应用对数变换(c=1)后的结果。

幂律（伽马）变换的基本形式为：
\[s = c r^{\gamma}\]
其中， \(c\) 和 \(\gamma\) 为正常数。r是输入灰度级，s是输出灰度级。有时考虑到偏移量（即输入为 0 时的一个可度量输出），上式也写为 \(s=c(r+\varepsilon)^{\gamma}\) ，不过，偏移量一般是显示标定(display calibration)问题。

不同伽马 \(\gamma\) 值（假设所有情形 c=1）的 \(s = c r^{\gamma}\) 曲线如图 8 所示。

从图 8 中分析易知， 当 \(\gamma\) 值小于 1 时，伽马变换把“范围较窄的低灰度值”转换为“范围较宽的灰度值”，这样可以增强图像较暗区域的对比度。当 \(\gamma\) 值大于 1 时，伽马变换把“范围较窄的高灰度值”转换为“范围较宽的灰度值”，这样可以增强图像较亮区域的对比度。

图像的直方图，表示了其灰度分布的特性。 直方图均衡(Histogram Equalization)的目的是使图像的灰度分布平均化（即各个灰度值的像素的个数差不多相等）。也就是说，直方图均衡就是把给定图像的直方图分布改变成“均匀”分布直方图分布。

直方图均衡对于处理背景和前景太亮或者太暗的图像非常有用，这种方法尤其是可以显示 X 射线图像中更好的骨骼结构显示以及曝光过度或者曝光不足照片中更好的细节。图 12 是直方图均衡的一个实例图。

图 13 是直方图均衡的另一个例子。

直方均衡可以得到一个灰度直方图分布平均的图像，直方图均衡的结果是唯一的。这是优点，也可以说是缺点。优点在于，无需多余的参数，就可以实现对比度的增强。缺点在于，没有参数，对于结果无法调整。

图 15 是一个“直方图均衡”处理后效果不是很好的例子。

在上面例子中，我们发现原图中的一个细节已经显示出来了，只是整个图像偏白而已，如果能够把处理后的图像的直方图向左移一些，则效果可能会好很多（不会这么偏白）。但“直方图均衡”无法调整，为了解决这个不足，就有了直方图匹配（规定化）。

产生处理后的图像具有特殊直方图的方法称为直方图匹配（或称为直方图规定化）。

这里不详细介绍“直方图匹配”过程。不过，很多时候“直方图匹配”都是试凑的过程。

令 \(r\) 表示在区间 \([0, L-1]\) 上代表灰度值的一个离散随机变量，并令 \(p(r_i)\) 表示对应于 \(r_i\) 值的归一化直方图分量。图像的灰度均值（即平均灰度）为：
\[m = \sum_{i=0}^{L-1} r_i p(r_i)\]
灰度方差为：
\[\sigma^2 = \sum_{i=0}^{L-1} (r_i - m)^2 p(r_i)\]
灰度方差是图像对比度的度量。

对于均值和方差，通常直接从取样值来估计它们，而不必计算直方图。这些估计值称为“取样均值”和“取样方差”。它们可以根据基本的统计学知识由下面的形式给出：
\[m = \frac{1}{MN} \sum_{x=0}^{M-1} \sum_{y=0}^{N-1} f(x,y)\]
和
\[\sigma^2 = \frac{1}{MN} \sum_{x=0}^{M-1} \sum_{y=0}^{N-1} [ f(x,y) - m ]^2\]

前面提到过，空间滤波器由下面两部分组成：
(1)一个邻域（典型地是一个较小的矩形）；
(2)对该邻域包围的图像像素执行的预定义操作。

如果在图像像素上执行的预定义操作是线性操作，则该滤波器称为“线性空间滤波器”。否则，滤波器就称为“非线性空间滤波器”。

在执行线性空间滤波时，必须清楚地理解两个相近的概念：一个是相关(Correlation)，另一个是卷积(Convolution)。如上面描述的那样， 相关是滤波器模板移过图像并计算每个位置乘积之和的处理（前文已经介绍过）。卷积的机理相似，但滤波器首先要旋转 \(180^{\circ}\) 。 下面将详细介绍“相关”和“卷积”处理过程。

一个大小为 \(m \times n\) 的滤波器 \(w(x,y)\) 与一幅图像 \(f(x,y)\) 做“相关”操作，可表示为 \(w(x,y) \mathbin{☆} f(x,y)\) ，在上一节中已经给出过其公式，为方便起见重写如下：

\begin{equation} w(x,y) \mathbin{☆} f(x,y) = \sum_{s=-a}^{a} \sum_{t=-b}^{b} w(s,t) f(x+s, y+t) \tag{Correlation Formula} \end{equation}

这一等式对所有位移变量 \(x\) 和 \(y\) 求值，以便 \(w\) 的所有元素访问 \(f\) 的每一个像素，其中假设 \(f\) 已被适当地填充（比如在边缘进行零填充，后文将介绍）。式中 \(a=(m-1)/2, b = (n-1)/2\) ，且假设 \(m,n\) 是奇整数。

类似地， \(w(x,y)\) 和 \(f(x,y)\) 的“卷积”表示为 \(w(x,y) \star f(x,y)\) ，它的定义由下式给出：

\begin{equation} w(x,y) \star f(x,y) = \sum_{s=-a}^{a} \sum_{t=-b}^{b} w(s,t) f(x-s, y-t) \tag{Convolution Formula} \end{equation}

其中，等式右侧的减号表示翻转 \(f\) （即旋转 \(180^{\circ}\) ）。为简化符号表示，我们遵循惯例，翻转和移位 \(w\) 而不是 \(f\) ，得到的结果是一样的（因为“卷积”满足交换律）。

显然， 如果滤波器矩阵是对称矩阵，则“卷积”操作和“相关”操作将得到相同的结果。

ImageMagick 和 GIMP 支持定制线性空间滤波的模板。

参考：
ImageMagick v6 Examples -- Convolution of Images
GIMP Convolution Matrix

平滑线性空间滤波器使用滤波模板确定的邻域内像素的平均灰度值代替图像中每个像素的值，这种处理的结果降低了图像灰度的“尖锐”变化。这些滤波器有时也称为均值滤波器(Averaging Filters)。显然，它们属于低通滤波器(Lowpass Filters)。

图 23 是两个 \(3 \times 3\) 平滑线性滤波器的示例。左边的平滑线性滤波器的所有系数都相等（都为 \(1/9\) ），有时也称为盒状滤波器(Box Filter)。右边的平滑线性滤波器实现是的“加权平均”，突出了处于模板中心位置的像素的重要性比其它像素的重要性要更大。不过，实践中，由于模板在一幅图像中任何一个位置所跨过的区域很小，通常很难看出使用图 23 中两个模板或类似方式进行平滑处理后的图像的区别。

图 24 是使用不同尺寸的模板对图像进行平滑处理后的效果。

高斯模糊(Gaussian blur)是图像处理常用滤波器，可以用它来 减少图像噪声，降低细节层次。

二维高斯函数为：
\[G(x,y) = \frac{1}{2 \pi \sigma^2} e^{- \frac{x^2 + y^2}{2 \sigma^2}}\]
它是各向同性的(isotropic)，函数图形如 25 所示。

图 26 是高斯滤波器（标准差 \(\sigma=1\) ，核大小为 \(5\times 5\) ）的一个离散近似。

下面是高斯滤波器（标准差 \(\sigma=0.84089642\) ，核大小为 \(7\times 7\) ）的另一个离散近似。

0.00000067	0.00002292	0.00019117	0.00038771	0.00019117	0.00002292	0.00000067
0.00002292	0.00078634	0.00655965	0.01330373	0.00655965	0.00078633	0.00002292
0.00019117	0.00655965	0.05472157	0.11098164	0.05472157	0.00655965	0.00019117
0.00038771	0.01330373	0.11098164	0.22508352	0.11098164	0.01330373	0.00038771
0.00019117	0.00655965	0.05472157	0.11098164	0.05472157	0.00655965	0.00019117
0.00002292	0.00078633	0.00655965	0.01330373	0.00655965	0.00078633	0.00002292
0.00000067	0.00002292	0.00019117	0.00038771	0.00019117	0.00002292	0.00000067

显然，高斯滤波器的输出是“邻域像素的加权平均”，同时离中心越近的像素权重越高。因此， 相对于均值滤波高斯滤波器的平滑效果更柔和，而且边缘保留的也更好。

计算高斯函数的离散近似时，在 \(3\sigma\) 距离之外的离散值可以认为不起作用，所以， 用户不指定核大小时，核大小可以由 \(\sigma\) 自动决定，为 \((6\sigma + 1) \times (6\sigma + 1)\) 。

参考：
http://homepages.inf.ed.ac.uk/rbf/HIPR2/gsmooth.htm

统计排序滤波器是一种非线性空间滤波器，这种滤波器的响应以滤波器包围的图像区域中所包含的像素的排序为基础，然后使用统计排序结果决定的值代替中心像素的值。这一类中最知名的滤波器是中值滤波器(Median Filter)，正如其名暗示的那样，它是将像素邻域内灰度的中值代替该像素的值。

数字函数的微分描述的是“差别”，有不同的方式来定义这些“差别”。 不过，对于一阶微分的任何定义都必须满足以下几点：(1) 在恒定灰度区域的微分值为零；(2) 在灰度台阶或斜坡的开始位置微分值非零；(3) 沿着斜坡的微分值为非零。类似地，任何二阶微分的定义必须满足以下几点：(1) 在恒定区域微分值为零；(2) 在灰度台阶或斜坡的开始和结束位置微分值非零；(3) 沿着斜坡的微分值为零。

对于一维函数 \(f(x)\) ，其一阶微分的基本定义是差值：
\[\frac{\partial f}{\partial x} = f(x+1) - f(x)\]
对上式关于 \(x\) 再次微分，得到二阶微分：
\[\begin{aligned} \frac{\partial^2 f}{\partial x^2} & = f'(x+1) - f'(x) \\ & = f(x+2) - f(x+1) - f(x+1) + f(x) \\ & = f(x+2) - 2f(x+1) + f(x) \\ \end{aligned}\]
但这一展开是关于点 \(x+1\) 的，我们感兴趣的是关于点 \(x\) 的二阶微分，故将上式中变量都减 1。即，我们将一维函数的二阶微分定义为如下差值：
\[\frac{\partial^2 f}{\partial x^2} = f(x+1) + f(x-1) - 2f(x)\]

注：这里使用偏导数符号仅仅是为了和二维函数在形式上一致。

容易验证这两个定义都满足前面所说的条件。

图 28 是一幅图像中一段水平灰度剖面对应一维数字函数的一阶微分和二阶微分。

在上面例子中，在灰度台阶的过渡中，二阶微分出现了“零交叉”，这对于边缘定位是非常有用的。

一个二维图像函数 \(f(x,y)\) 的拉普拉斯算子定义为：
\[\nabla^2 f = \frac{\partial^2 f}{\partial x^2} + \frac{\partial^2 f}{\partial y^2}\]
由于任意阶微分都是线性操作，所以拉普拉斯算子也是一个线性算子。为了以离散形式描述这一公式，和前面介绍的一维函数二阶微分的定义类似，对二维函数二阶微分可以这样定义，在 \(x\) 方向上：
\[\frac{\partial^2 f}{\partial x^2} = f(x+1,y) + f(x-1,y) - 2f(x,y)\]
在 \(y\) 方向上：
\[\frac{\partial^2 f}{\partial y^2} = f(x,y+1) + f(x,y-1) - 2f(x,y)\]
由上面三个公式，从而得 两个变量的离散拉普拉斯算子为：
\[\nabla^2 f(x,y) = f(x+1,y) + f(x-1,y) + f(x,y+1) + f(x,y-1) - 4f(x,y)\]
图 29 所示模板是上面离散拉普拉斯算子的滤波模板实现。

通过拉普拉斯算子得到的图像（称为拉普拉斯图像），体现了原图像灰度的突变。那些原图像中灰度恒定的区域在拉普拉斯图像中对应区域的灰度会为零（表现为黑色）。

图 30 是拉普拉斯图像的一个例子。

在实践中，除了图 29 所示的拉普拉斯模板外，还有其它三个拉普拉斯模板。图 31 是常用的拉普拉斯滤波模板实现的总结。

将拉普拉斯滤波后的图像和原图像合并（相加或相减）时，可以达到锐化原图像的效果。至于到底是相加还是相减，取决于拉普拉斯算子中心系数的符号（当中心系数的符号为负时，两图像相减；当中心系数的符号为正时，两图像相加）。即，用拉普拉斯滤波模板进行锐化处理的方法可用下式表示：
\[g(x,y) = \begin{cases} f(x,y) - \nabla^2 f(x,y) \qquad \text{如果拉普拉斯模板的中心系数为负} \\ f(x,y) + \nabla^2 f(x,y) \qquad \text{如果拉普拉斯模板的中心系数为正} \end{cases}\]

图 32 是用拉普拉斯滤波模板进行锐化处理的实例。从图中可知，右图锐化效果更加明显，因为它使用的是拉普拉斯“扩展模板”，在对角方向上产生了额外的锐化。

非锐化掩蔽(Unsharp Masking)是在印刷和出版界已用了多年的图像锐化处理技术。非锐化掩蔽的处理过程如下：
第 1 步：模糊原图像；
第 2 步：从原图像中减去模糊图像（产生的差值图像称为模板）；
第 3 步：将模板加到原图像上，将得到锐化图像。

非锐化掩蔽的原理如图 33 所示。

令 \(\overline{f}(x,y)\) 表示模糊后的图像，非锐化掩蔽用公式描述如下。首先，得到模板：
\[g_{\text{mask}}(x,y) = f(x,y) - \overline{f}(x,y)\]
然后，在原图像上加上该模板的一个权重部分：
\[g(x,y) = f(x,y) + k \times g_{\text{mask}}(x,y) \qquad (k \ge 0)\]
上式中，当 \(k=1\) 时，就是前面介绍的非锐化掩蔽。当 \(k < 1\) 时，该处理过程不强调钝化模板的贡献。 当 \(k > 1\) 时，该处理过程称为高提升滤波(Highboost Filtering)。

图像处理中的一阶微分是用梯度幅值来实现的。对于函数 \(f(x,y)\) ， \(f\) 在坐标 \((x,y)\) 处的梯度定义为二维列向量：
\[\nabla f \equiv \text{grad}(f) \equiv \begin{bmatrix} \frac{\partial f}{\partial x} \\ \frac{\partial f}{\partial y} \end{bmatrix}\]
该向量具有重要的几何特性，即它指出了在位置 \((x,y)\) 处 \(f\) 的最大变化率的方向（函数值增加最快的方向）。
向量 \(\nabla f\) 的幅度值(Magnitude)表示为 \(M(x,y)\) ，即：
\[M(x,y) = \sqrt{ \left(\frac{\partial f}{\partial x} \right)^2 + \left(\frac{\partial f}{\partial y} \right)^2}\]
它是梯度向量方向变化率在 \((x,y)\) 处的值。 当 \(x\) 和 \(y\) 在 \(f\) 中所有像素变化时计算每个 \(M(x,y)\) 得到与原图像大小相同的图像，该图像通常称为梯度图像（当含义很清楚时，可简称为梯度）。
在某些实现中，用绝对值来近似平方和平方根操作更适合计算，即：
\[M(x,y) \approx |\frac{\partial f}{\partial x}| + |\frac{\partial f}{\partial y}| \]

和拉普拉斯滤波的情况类似，我们需要对上面公式定义一个离散近似，并由此形成合适的滤波模板。
为了简化下面的讨论，我们将使用图 34 中的符号来表示一个 \(3 \times 3\) 区域内的图像点的灰度。使用这些符号，中心点 \(z_5\) 表示任意位置 \((x,y)\) 处的 \(f(x,y)\) ， \(z_1\) 表示为 \(f(x-1,y-1)\) ，等等。

由节 2.5.1 中介绍的知识可知，一阶微分的最简近似是 \(\frac{\partial f}{\partial x} = z_8 - z_5\) 和 \(\frac{\partial f}{\partial y} = z_6 - z_5\) 。在早期数字图像处理的研究中，由 Roberts 提出的其他两个定义使用了交叉差分： \(\frac{\partial f}{\partial x} = z_9 - z_5\) 和 \(\frac{\partial f}{\partial y} = z_8 - z_6\) ，用线性滤波器模板表示，如图 35 所示，它们被称为罗伯特交叉梯度算子(Roberts cross gradient operators)。

这时，对应的梯度图像为：
\[M(x,y) = \sqrt{(z_9 - z_5)^2 + (z_8 - z_6)^2}\]
或者近似为：
\[M(x,y) = |z_9 - z_5| + |z_8 - z_6|\]

前面介绍的微分算子（如拉普拉斯算子，Sobel 算子），其模板系数相加都为零。这体现了节 2.5.1 所描述的定义微分的条件（恒定灰度区域的微分值为零）。

基本集合操作包括：
(1) 子集。如果集合 A 中的每个元素又是另一个集合 B 中的一个元素，则称 A 为 B 的子集，表示为：
\[A \subseteq B\]
(2) 并集。如果某个集合 C 包含集合 A 和 B 中的所有元素，则称 C 为集合 A 和 B 的并集，表示为：
\[C = A \cup B\]
(3) 交集。如果某个集合 D 包含的元素同时属于集合 A 和 B，则称 D 为集合 A 和 B 的交集，表示为：
\[D = A \cap B\]
(4) 补集。集合 A 的补集是不包含于集合 A 的元素所组成的集合，表示为：
\[A^c = \{ w \mid w \notin A \}\]
(5) 差集。集合 A 和 B 的差集是属于集合 A，而不属于集合 B 的元素组成的集合，表示为：
\[A-B = \{ w \mid w \in A, w \notin B \} = A \cap B^c\]

除了上面这些基本的操作外，集合的“反射(Reflection)”和“平移(Translation)”的概念在形态学中应用非常广泛。
一个集合 B 的反射(Reflection)可记为 \(\hat{B}\) ，定义如下：
\[\hat{B} = \{ w \mid w = -b, b \in B \}\]
如果 B 是描述图像中物体的像素（二维点）的集合，则 \(\hat{B}\) 是 B 中 \((x,y)\) 坐标被 \((-x,-y)\) 替代的点的集合。

集合 B 按照点 \(z=(z_1, z_2)\) 的平移(Translation)可记为 \((B)_z\) ，定义如下：
\[(B)_z = \{ c \mid c = b + z, b \in B \}\]
如果 B 是描述图像中物体的像素（二维点）的集合，则 \((B)_z\) 是 B 中 \((x,y)\) 坐标被 \((x+z_1, y+z_2)\) 替代的点的集合。

集合反射和平移例子如图 37 所示。

在形态学中， 结构元(Structuring Element)是一个小集合（子图像），用它来研究图像中感兴趣的特征。

结构元都有一个“原点”，一般用实心黑点表示（如果结构元是对称的，且未明确画出原点时，常常假定结构元原点就是结构元的对称中心）。

当对图像操作时，为便于操作，一般要求结构元是矩形阵列，这可以通过“添加最小可能数量”的背景元素形成一个矩阵阵列来实现。

用结构元对图像进行处理时，图像边界处也往往需要添加适当的背景元素。

如何使用结构元对图像进行处理呢？我们通过一个例子来说明。
假定我们定义一个用结构元 B 在集合 A 上的操作如下：通过让 B 在 A 上运行，以便结构元 B 的原点访问 A 的每一个元素，来创建一个新集合。在 B 的每个原点位置，如果 B 完全被 A 包含，则将该位置标记为新集合的一个成员；否则，将该位置标记为非新集合的成员。 图 39 显示了这一操作结果。我们注意到，当 B 的原点位于 A 的边界元素上时，B的一部分将不再包含在 A 中，按照约定的规则，从而 B 处在中心位置的点没有成为新集合的成员。 最终结果是集合 A 的边界被腐蚀掉了。

关于“腐蚀”的数字公式表述在后文将介绍。

A 和 B 是 \(\mathbb{Z}^2\) 中的集合，集合 B 对集合 A 的“腐蚀”定义为：
\[A \ominus B=\{z \mid (B)_{z} \subseteq A\}\]
即 B 对 A 的腐蚀是一个用 z 平移的 B，那些满足“B 包含在 A 中”的所有的点 z 的集合。
腐蚀还有其他一些等价的定义形式，如：
\[A \ominus B=\{z \mid (B)_z \cap A^c = \varnothing \}\]

图 40 是用结构元 B 对图像 A 进行腐蚀的两个例子。

图 41 是用形态相同但大小不同的结构元对同一图像进行腐蚀的例子。

A 和 B 是 \(\mathbb{Z}^2\) 中的集合，集合 B 对集合 A 的“膨胀”定义为：
\[A \oplus B=\{z \mid (\hat{B})_{z} \cap A \neq \varnothing \}\]
集合 B 对集合 A 的“膨胀”就是用 z 平移的 B 关于原点的映像 \(\hat{B}\) ，那些满足“ \(\hat{B}\) 和 A 至少有一个元素是重叠的”的所有的点 z 的集合。

图 42 是用结构元 B（它是关于原点对称的，所以 \(\hat{B}=B\) ）对图像 A 进行膨胀的两个例子。

可以证明，腐蚀和膨胀彼此关于集合求补运算和反射运算是对偶的，即：
\[\begin{aligned} (A \ominus B)^c = A^c \oplus \hat{B} \\ (A \oplus B)^c = A^c \ominus \hat{B} \end{aligned}\]

上面两式表明，B对 A 的腐蚀是 \(\hat{B}\) 对 \(A^c\) 的膨胀的补；B对 A 的膨胀是 \(\hat{B}\) 对 \(A^c\) 的腐蚀的补。

当结构元 B 关于其原点对称时（通常如此），有 \(\hat{B}=B\) ，这时对偶性特别有用。这样，我们可以用相同的结构元简单地使用 B 膨胀图像 A 的背景（即膨胀 \(A^c\) ），并对该结果求补就可得到 B 对图像 A 的腐蚀。

结构元 B 对集合 A 的开操作，表示为 \(A\circ B\) ，其定义为：
\[A\circ B=(A\ominus B) \oplus B\]
即开操作是“先腐蚀后膨胀”。

开操作有一个简单的几何解释。假设我们把结构元 B 视为一个“转球”， \(A\circ B\) 的边界由 B 中的点建立： 当 B 在 A 的边界内侧滚动时，B所能达到的 A 的边界的最远点组成了 \(A\circ B\) 的边界。 图 43 是开操作的一个实例。

结构元 B 对集合 A 的开操作还可以表示为：
\[A\circ B= \bigcup \{ (B)_z \mid (B)_z \subseteq A \}\]
其中，记号 \(\bigcup \{ \cdot \}\) 表示大括号中所有集合的并集。

结构元 B 对集合 A 的闭操作，表示为 \(A\bullet B\) ，其定义为：
\[A\bullet B=(A\oplus B) \ominus B\]
即闭操作是“先膨胀后腐蚀”。

和开操作类似，闭操作也有一个简单的几何解释（不同在于开操作是沿“内侧边界”滚动，而闭操作是沿“外侧边界”滚动）。如图 44 所示。

和腐蚀和膨胀类似，开操作和闭操作彼此关于集合求补运算和反射运算是对偶的，即：
\[\begin{aligned} (A \circ B)^c = A^c \bullet \hat{B} \\ (A \bullet B)^c = A^c \circ \hat{B} \end{aligned}\]

图 45 是形态学基本操作（腐蚀、膨胀、开操作、闭操作）的实例。

当结构元 b 的原点位于 \((x,y)\) 处时，用一个平坦结构元 b 在 \((x,y)\) 处对图像 f 的“腐蚀”的定义为：图像 f 中与 b 重合区域的最小值。灰度图像的“腐蚀”用公式表示为：
\[[f \ominus b] (x,y) = \min_{(s,t) \in b} \{ f(x+s, y+t) \}\]

平坦结构元 b 在 \((x,y)\) 处对图像 f 的“膨胀”的定义为：图像 f 中与 \(\hat{b}\) （即 \(b(-x, -y)\) ）重合区域的最大值。灰度图像的“膨胀”用公式表示为：
\[[f \oplus b] (x,y) = \max_{(s,t) \in b} \{ f(x-s, y-t) \}\]

灰度图像的开操作和闭操作的表达式与二值图像的对应操作具有相同的形式。
结构元 b 对图像 f 的开操作为：
\[f \circ b = (f \ominus b) \oplus b\]
结构元 b 对图像 f 的闭操作为：
\[f \bullet B = (f \oplus b) \ominus b\]

灰度图像的开操作和闭操作同样具有简单的几何解释。假设我们将一个图像 \(f(x,y)\) 看成是一个三维表面，也就是说，它的灰度值可解释为 \(xy\) 平面上的高度值。
b 对 f 的开操作的几何解释为： 在 b 的原点位置，从 f 的下表面垂直向上推动结构元 b 时，b的任何部分所达到的最高值。当 b 的原点访问 f 的每一个坐标 \((x,y)\) 所得到的所有值组成的集合就是 b 对 f 的开操作的结果。 闭操作的几何解释和开操作类似，只是“向上推”改为“向下推”，能达到的“最高值”改为“最低值”。

图 46 以一维的形式说明了灰度图像的开操作和闭操作的几何解释。

显然， 开操作可去除“比结构元小”的明亮细节，而闭操作可去除“比结构元小”的暗部细节。

使用拉普拉斯算子（参考 2.5.2 ），如果在某点处该滤波器（模板）的响应的绝对值超过了一个指定的阈值，那么可以认为模板中心位置 \((x,y)\) 处的该点已经被检测到了。 在输出图像中，这样的点被标注为 1（表示找到孤立点），而所有其他点则被标注为 0，从而产生一幅二值图像。

一般来说，我们对特定方向的线感兴趣。图 51 所示的几个模板可以用来检测特定方向的线。

一阶导数和二阶导数都可以找到边缘，如图 52 所示。

在灰度斜坡开始处并沿着整个灰度斜坡，一阶导数都不为零，这会导致一阶导数产生比较“粗”的边缘。

对于图像中的每条边缘，二阶导数生成两个值（被称为双线效应，Double-line effect，这是一个副作用，需要处理）； 二阶导数的正负符号可以用于确定边缘的方向（从亮到暗，还是从暗到亮）；二阶导数的“零交叉点”可以用于定位边缘的中心。

边缘是灰度变化较大的位置。用一阶导数（即梯度）可以检测灰度变化。相关基本知识可参考 2.5.4 。

这里列举一些常见的边缘检测梯度算子，如图 53 所示。

节 2.5.2 中介绍过，使用二阶微分可以提取出图像中灰度变化较快的部分。但是，直接对图像使用拉普拉斯算子 \(\nabla^2 f = \frac{\partial^2 f}{\partial x^2} + \frac{\partial^2 f}{\partial y^2}\) ，得到的结果中图像的噪声很可能被放大。

一个更好的作法是先使用“高斯模糊滤波器（参考2.4.2）”抑制图像中的噪声，再使用二阶微分得到其拉普拉斯图像，最后寻找拉普拉斯图像中的零交叉点即可得到图像边缘。这就是 Marr-Hildreth 边缘检测算法的思想。

Marr-Hildreth 边缘检测算法 总结如下：
第一步：用“高斯模糊滤波器 \(G(x,y)\) ”对输入图像 \(f(x,y)\) 进行平滑处理。
第二步：计算去噪后图像（第一步的结果）的拉普拉斯图像：
\[g(x,y) = \nabla^2 [G(x,y) \star f(x,y)]\]
其中， \(\star\) 是卷积运算（由于高斯模糊滤波器是对称的，“相关”和“卷积”的结果是相同的，一般提到 Marr-Hildreth 算法时，习惯于使用卷积一词）。
第三步：计算 \(g(x,y)\) 的“零交叉点”，即为图像中物体边缘。

注：第一、二步往往合并为一步，后文会介绍。
第三步中，如何寻找 \(g(x,y)\) 中的零交叉点呢？ 由于噪声或计算不准确，试图通过查找满足 \(g(x,y)=0\) 的坐标都找零交叉点是不实际的。 寻找零交叉点的一种可行方法描述如下：
在滤波后的图像 \(g(x,y)\) 的任意像素 \(p\) 处，观察以 \(p\) 为中心的一个 \(3 \times 3\) 邻域，如果 \(p\) 为零交叉点，则意味着至少存在“相对的邻域像素的正负符号不同”，这有 4 种情况需要测试：左/右、上/下、和两个对角。存在一种情况“相对的邻域像素的正负符号不同”，并且它们的数值差的绝对值还超过某个阈值（这很重要），则我们称 \(p\) 为一个零交叉像素。

Canny edge detector是目前讨论过的边缘检测算法中最优秀的。

Canny 算法分为如下几个步骤：
(1) 用高斯滤波器平滑输入图像；
(2) 计算梯度幅值图像和梯度角度图像；
(3) 对梯度幅值图像应用“非最大抑制(Non-maximum suppression)”算法；
(4) 用“双阈值处理(Double threshold)”和来检测边缘；
(5) 用特征综合方法收集来自多尺度的最终的边缘信息（不过，通常的算法实现都省略了这一步）。

由第(2)步知，Canny 算法是利用一阶导数来检测边缘的，但由于一阶导数得到的边缘往往会较粗； 第(3)步使用非最大抑制(Non-maximum suppression)算法的目的是“细化上一步得到的边缘”。 后文将介绍非最大抑制(Non-maximum suppression)算法。

“双阈值处理”的目的是减少伪边缘点，其具体过程为：设置两个阈值，上限和下限。如果一个像素的梯度大于上限阈值，则认为是边缘像素；如果低于下限阈值，则不认为是边缘像素，被抛弃；如果介于二者之间，只有当其“与高于上限阈值的像素连接”时才被接受。

图 55 是用不同算法对同一图像进行边缘检测的例子。从这个例子可知 Canny 算法的效果最好。

霍夫变换(Hough Transform)可以用来寻找图像中的直线，扩展后可以用来寻找其他的形状（如圆形等等）。

下面用检测直线来说明霍夫变换的基本思想。假设直线由两点 \(A=(x_1,y_1)\) 和 \(B=(x_2,y_2)\) 定义（如图 57 (a)所示）。 通过点 \(A\) 的所有直线可用 \(y_1 = k x_1 + q\) 表示。 同一个方程可以解释为参数空间 \(k,q\) 的方程，因此通过点 \(A\) 的所有直线可以表示为方程 \(q = -x_1 k +y_1\) （如图 57 (b)所示），类似地， 通过点 \(B\) 的所有直线可以表示为方程 \(q = -x_2 k +y_2\) 。在参数空间 \(k,q\) 中，两条直线的唯一公共点是在原图像空间中表示连接点 \(A\) 和 \(B\) 的唯一存在的直线。这样，我们在参数空间中检测相交的点就相当于在原图像空间中检测直线，后文将介绍这个过程。

不过，当直线逼近垂直方向时，直线的斜率 \(k\) 会趋于无限大，比较难处理，这时可以把直线表达为下面形式（这时参数空间为 \(\rho, \theta\) ）： \(x \cos \theta + y \sin \theta = \rho\) ，经过原图像空间某点的所有直线，在参数空间 \(\rho, \theta\) 中是一条正弦曲线，如图 58 所示。

一般来说，由于噪声或计算误差，原图像空间中的直线不会精确地被变换到参数空间中的一个点，而是一个点群，点群中心可以作为原图像空间直线的参数。我们在参数空间寻找点群的方法是：把 \(\rho, \theta\) 空间划分为累加单元（如图 58 (c)所示），哪个累加单元的点非常多，就说明找到了原图像空间的直线。

详情可参考：《图像处理、分析与机器视觉（第 3 版，艾海舟等译）》，6.2.6 霍夫变换

假设某图像 \(f(x,y)\) 由暗色的背景和较亮的物体组成，其灰度直方图如图 60 (a)所示，从这个图像的背景中提取物体的一种明显方法是选择一个合适的阈值 \(T\) ，使得背景点都满足 \(f(x,y) \le T\) ，而物体像素点都满足 \(f(x,y) > T\) 。换句话说，分割后的图像 \(g(x,y)\) 由下式给出：
\[g(x,y) = \begin{cases} 1, & f(x,y) > T \\ 0, & f(x,y) \le T \end{cases}\]

如果 \(T\) 是一个适用于整个图像的常数，则上面公式给出的处理称为 “全局阈值处理” 。当 \(T\) 值在一幅图像上改变时，称为 “可变阈值处理” （有时也称为“局部阈值处理”或者“区域阈值处理”）。可以把阈值从单个扩展到多个（如图 60 (b)所示），这称为“多阈值处理”。

当物体和背景像素的灰度分布相差明显时，可以用全局阈值处理来分割背景和物体。但如何找到一个合适的阈值呢？
可以使用下面的迭代算法来寻找全局阈值 \(T\) ：
(1) 为全局阈值 \(T\) 选择一个初始估计值；
(2) 用阈值 \(T\) 分割该图像。这将产生两组像素： \(G_1\) 由灰度值大于 \(T\) 的所有像素组成， \(G_2\) 由灰度值小于等于 \(T\) 的所有像素组成。
(3) 对 \(G_1\) 和 \(G_2\) 的像素分别计算平均灰度值，记为 \(m_1\) 和 \(m_2\) 。
(4) 计算一个新的阈值：
\[T = \frac{1}{2} (m_1 + m_2)\]
(5) 重复步骤 2 到步骤 4，直到连续迭代中的 \(T\) 值间的差小于一个预定义的参数 \(\Delta T\) 为止。

图 63 是全局阈值处理的一个例子。

Otsu方法 是由日本学者 Nobuyuki Otsu 于 1979 年提出的一种对图像进行全局阈值处理的算法。

Otsu 方法的基本思想为：按图像的灰度特性，将图像分成背景和物体二个部分时，背景和目标之间的方差越大，说明图像的这二个部分的差别会越大，分类效果越好。 但怎么描述背景和目标之间的方差（称为类间方差）呢？请看下文说明。

介绍“类间方差”前，先介绍一些基本知识。
假设图像阈值化处理后为两类 \(C_1\) 和 \(C_2\) ，其中 \(C_1\) 由图像中灰度值在范围 \([0,k]\) 内的所有像素组成， \(C_2\) 由灰度值在范围 \([k+1, L-1]\) 内的所有像素组成。则，像素被分到类 \(C_1\) 中的概率 \(P_1(k)\) 为：
\[P_1(k) = \sum_{i=0}^{k} p_i\]
其中， \(p_i\) 是归一化直方图中灰度为 \(i\) 时的值。
像素被分到类 \(C_2\) 中的概率 \(P_2(k)\) 为：
\[P_2(k) = 1 - P_1(k)\]

记 \(m_1(k)\) 和 \(m_2(k)\) 分别为 \(C_1\) 和 \(C_2\) 的灰度均值， \(m_G\) 为整个图像的灰度均值（它和 \(m_1(k)\) 和 \(m_2(k)\) 之间存在关系： \(P_1(k) m_1(k) + P_2(k) m_2(k) = m_G\) ）。

由节 2.2.3 知识可知，整个图像的灰度方差为：
\[\sigma_G^2 = \sum_{i=0}^{L-1} (i - m_G)^2 p_i\]

为了使记号不过于累赘，下面的讨论我们暂时省略 \(k\) 。
记 \(\sigma_B^2\) 为类间方差(between-class variance)，其定义为：
\[{\color{red}{\sigma_B^2 = P_1 (m_1 - m_G)^2 + P_2 (m_2 - m_G)^2}}\]
利用关系式 \(P_1 m_1 + P_2 m_2 = m_G\) ，类间方差还可以化为下面形式：
\[\sigma_B^2 = P_1 P_2 (m_1 - m_2)^2\]
从上式可以看出， 两个类的灰度均值 \(m_1\) 和 \(m_2\) 彼此相隔越远时， \(\sigma_B^2\) 会越大，这表明类间方差 \(\sigma_B^2\) 是类之间的可分性度量。

\(C_1\) 和 \(C_2\) 类是以阈值 \(k\) 为分类依据的， \(k\) 取不同的值时，可得到不同的类间方差 \(\sigma_B^2(k)\) ， Otsu 方法的核心 就是求当 \(\sigma_B^2(k)\) 取得最大值时，对应的阈值 \(k^{*}\) 作为最终的阈值。用公式表达就是：
\[\sigma_B^2(k^*) = \max_{0 \le k \le L-1} \sigma_B^2(k)\]

图 65 (a)是一个被噪声污染得比较严重的图像，其直方图中已经找不到不同的波峰，用基本全局阈值得或 Otsu 方法都无法成功分割背景和物体。这时，可以尝试在全局阈值处理前，先对原图像进行平滑处理。如图 65 (c)是对(a)图使用 \(5\times5\) 均值模板进行平滑处理后的结果，这时直方图已经出现了不同的波峰，可以使用全局阈值处理来分割背景和物体了。

节 4.2.4 中介绍的利用平滑方法并不能改善所有的全局阈值处理，如图 66 所示，背景中的物体太小，经过平滑处理后，没有产生新的波峰，还是无法使用全局阈值处理。

这时怎么办呢？改进直方图形状的另一种思路是： 仅考虑那些位于或靠近“物体和背景的过渡位置（即边缘）”的像素，这些部分像素构成一个子图，用全局阈值处理（如 Ostu 方法）从这个子图中得到一个阈值，再利用这个阈值处理整个原始图像。

但，完整的边缘信息在用阈值处理进行分割操作期间是不可用的，因为寻找边缘正好是分割所要做的工作。利用梯度幅度图像或拉普拉斯绝对值图像可以得到边缘的基本信息。

下面通过一个例子来说明如何利用边缘来改善全局阈值处理。如图 67 所示。

可变阈值处理（又称为“局部阈值处理”或者“区域阈值处理”）是指阈值 \(T\) 在整个图像中是可变的。当阈值 \(T\) 取决于空间坐标 \((x,y)\) 时，可变阈值处理通常称为“动态阈值处理”，或者“自适应阈值处理”。

分水岭算法常常由于噪声和梯度的局部不规则性造成“过度分割”，这可能会使算法得到的结果变得毫无用处。如图 73 所示。

使用标记是控制过度分割的一种方法。一个标记是属于一幅图像的一个连通分量。把和感兴趣物体相联系的标记称为“内部标记”，与背景相关的标记称为“外部标记”。标记属于图像的先验知识（一般需要人为指定）。

利用标记改善图 73 的分水岭处理结果如图 74 所示。