图像处理基本概念

25spring BME3304 生物医学图像处理（I）与BME4301 计算机辅助手术与治疗技术的课程笔记

图像分割

图像分割是介于图像预处理（复原，增强）与更为复杂的图像理解（如识别分类理解等）之间的重要操作。图像分割基于的逻辑很简单：区域间具有不连续性（边缘检测），以及区域内具有相似性（区域生长）。

基于区域间不连续性的分割策略

区域间不连续性体现在有边缘，而有边缘进一步可以体现为存在灰度的突变，例如点，线，边沿。

点/线检测

比较常见的检测方式例如空间拉普拉斯算子。除此之外，在某些情况下我们也可能有对线的方向有一定要求的场景，这时我们可以采用一些其他的各向不同的算子，例如 $\begin{bmatrix} -1 & 2 & -1 \\ -1 & 2 & -1 \\ -1 & 2 & -1 \end{bmatrix}, \begin{bmatrix} 2 & -1 & -1 \\ -1 & 2 & -1 \\ -1 & -1 & 2 \end{bmatrix}$ 等。

边缘检测

边缘存在的形式可能有黑->白，也可能是黑->白->黑，不过边缘检测的原理都是导数会在边缘处取得最大值。之前已经介绍过的办法包括各类梯度算子（基础梯度算子，Roberts，Prewitt，Sobel 算子等，见前文）， 或者拉普拉斯算子（对噪声非常敏感不建议单独使用）。

LoG 算子

LoG（Laplacian of Gaussian）算子其可以有效地抑制噪声对边缘检测的影响。LoG 算子的模板可以通过将高斯模板与拉普拉斯模板进行卷积得到。顾名思义，LoG 算子就是先对其应用高斯滤波器，然后再计算拉普拉斯算子。其理论的算子公式：

（看着好复杂），实际使用的算子（$\sigma = 5$）如：

一般来说，在应用 LoG 算子之后通常还要通过检测零交叉来定位边缘，即看像素的左右或上下两对像素有没有出现零交叉（当然为了防止噪声影响也可以加一个阈值判定）

LoG 算子的高斯滤波阶段一定程度上克服了噪声的影响，继而增强了拉普拉斯算子的效果，但其仍然可能会产生假边缘，以及对一些曲线边缘的检测不够理想。

Canny 边缘检测

Canny 边缘检测是一种多阶段的边缘检测算法, 其目标是检测出图像中的真实边缘, 同时尽量减少噪声和假边缘的影响。其不是一个通用固定的算子。Canny 边缘检测的逻辑是提出了三个可以量化的边缘检测的损失函数, 边缘检测的目的就被转化为了求这三个函数的极值:

1. 信噪比越大，提取的边缘质量越高。信噪比 $\mathrm{SNR}$ 定义为：

   $$\mathrm{SNR} = \frac{\left| \int_{-W}^{+W} G(-x)h(x)dx \right|}{\sigma \sqrt{\int_{-W}^{+W} h^2(x)dx}}$$

   $G(x)$ 为表示边缘的函数，$h(x)$ 为滤波器的脉冲响应。假设边缘的中心在 $x=0$ 外，滤波器只在 $[-W, W]$ 范围内响应。

2. 标记为边缘的点与真实的边缘中心距离最小: 边缘定位精度定义为 $L$：

   $$L = \frac{\left| \int_{-W}^{+W} G'(-x)h'(x)dx \right|}{\sigma \sqrt{\int_{-W}^{+W} h'^2(x)dx}}$$

   L 越大表明分类精度越高

3. 单个边缘点响应: 真实边缘周围的局部最大数应该最小。 要保证对单边缘只有一个响应，边缘检测算子的脉冲响应的导数的零交叉点平均距离应满足：

   $$D_{\text{ZCA}}(f') = \pi \left( \frac{\int_{-\infty}^{+\infty} h'^2(x)dx}{\int_{-W}^{+W} h''(x)dx} \right)^{1/2}$$

实际的 Canny 算子实现步骤如下:

1. 使用高斯滤波器平滑图像以减少噪声。
2. 计算图像的梯度强度和方向, 通常使用 Sobel 算子。得到每个像素位置的梯度幅值 $M(x, y)$ 和方向 $\theta(x, y)$。
3. 非极大值抑制: 对梯度幅值图像进行非极大值抑制, 以细化边缘。对于每个像素, 检查其梯度方向上的邻域像素, 如果当前像素的梯度幅值不是局部最大值, 则将其设为 0。这个步骤可以使边缘变得更细更准确。
4. 双阈值处理: 应用两个阈值（高阈值和低阈值）将像素分为三类：强边缘像素（大于高阈值）、弱边缘像素（介于低阈值和高阈值之间）和非边缘像素（小于低阈值）。
5. 抑制孤立弱边缘: 从强边缘部分开始，沿着强边缘末端的梯度方向检查，若存在弱边缘像素，则将其也标记为边缘，直到边缘无法继续延伸，以连接断裂的边缘。

边缘连接

前述的边缘连接算法虽然部分会有连接边缘的功能, 但有时仍然会出现边缘断裂的情况. 这时我们可以通过一些后处理算法来连接这些断裂的边缘.

局部连接处理

首先，我们需要确定哪些是“真边缘”，这种方法认为，边缘中相邻的边缘点的梯度大小应该大于某阈值，且梯度方向相似。这样我们可以得到大量断续的边缘片段，然后通过连接这些边缘片段来形成完整的边缘。具体步骤如下：

随后逐行扫描图片，若每行中存在 <L 的非边缘片段则填充为边缘。

注意这里我们实际实现时，我们看的不是像素之间相对的梯度方向差异，而是一次看所有像素是否在一个方向范围内，因此，这样遍历一次只会识别冰并填充单个方向上的边缘，一般来说需要对图片进行一定角度的旋转反复重复这一行为。

全局连接处理（霍夫连接）

我们把之前找到的断断续续的边界点视为一个点集，那我们可以比较自然地想到，找到真实的边界（假设为直线）就是找到能经过这个点集里最多点的那一条。不过，直接在笛卡尔坐标系中寻找直线计算量相当大(o(n^2)), 因此我们可以将问题转换到极坐标系/参数空间中进行处理。

将问题转换到其他坐标系的核心逻辑是实现点/线的转换，根据点找过点线计算量大，但找线的交点是一件非常轻松的事情，以极坐标系为例，过点的直线簇 $\rho = x \cos \theta + y \sin \theta$ 在 $\theta$ 为 x 轴， $\rho$ 为 y 轴的坐标系中表现为一条正弦曲线，其上每个点都是一个过点的直线，因此我们只需要画出点集中每个点的直线簇所对应的正弦曲线然后找交点即可。

虽然上述描述只考虑了直线的情况，但实际上我们可以将绝大部分形状都转换为参数空间中的曲线簇，然后通过找交点的方式来确定真实边界的位置。其抗噪声能力显著优于局部连接处理等手段。不过，边界越复杂，其计算量会极大幅度上升。

基于区域内相似性的分割策略

阈值处理

边缘信息改进全局阈值处理

对于有些直接采取阈值处理效果不佳的情况，例如目标与背景面积相差较大等，我们知道，其实我们最关心的部分，就是边缘位置上的分割准确性，如果我们选择的这个阈值在边缘上能良好表现，那么其在背景与目标内部自然也能有不错的表现，所以我们可以利用拉普拉斯或者其他梯度算子来先提取边界，然后再只用边界来计算阈值，再将这个阈值应用到原图片上。

局部图像性质可变阈值处理

也就是不同分区手动设定不同的阈值/阈值计算方法。

基于移动平均的可变阈值处理

也就是滑动窗口计算阈值。

基于区域的分割

连通域标记法 [only in BME4301]

算法的核心是找到图像中所有满足指定条件（IC）的像素，并且将其分为若干独立连通区域。即所有选中的像素满足的是同一个条件，但其分属于不同的区域会导致其连通域标签不同。目标图像通常为已二值化，且定义了邻域连接方式（如 4/8 邻域）。：

1. 从上到下从左到右遍历，对每一个满足 IC 的像素 $P$，检查其已扫描的邻域像素（通常为上方和左侧）：

   • 新连通域生成： 若 $P$ 的所有邻域像素均不满足 IC（即均为背景），则判定 $P$ 为新连通区域的起始点，赋予其一个新的临时标签（Label）。
   • 标签继承： 若邻域中存在满足 IC 的像素，且仅有一种标签，则 $P$ 继承该标签，归入该连通域。
   • 冲突记录： 若邻域中存在多个满足 IC 的像素，且它们持有不同的标签（例如 $L_a$ 与 $L_b$），说明像素 $P$ 桥接了两个先前被认为独立的区域。此时，将 $P$ 标记为其中之一（通常取最小值），并同时在等价表中记录这两个标签的等价关系（$L_a \equiv L_b$）。

2. 等价关系解析（Resolve Equivalence）

   • 对等价标签对进行合并。
   • 将所有相互等价的临时标签归并为同一个等价类，并为每个等价类分配一个唯一的最终规范标签。

3. 再次遍历图像：

   • 查找每个像素当前的临时标签。
   • 根据步骤 2 建立的映射关系，将所有临时标签替换为其所属等价类的最终规范标签。

区域增长法

自下而上的区域分割方法，需要选定一个起始种子，设定生长规则和终止规则，然后不断地将符合规则的邻域像素加入区域，直到没有新的像素可以加入为止。

区域分离与聚合

更多是一种自上而下的分割策略，其可以解决需要手动指定种子的问题。其分为两种操作

• 基本逻辑
  • 分离（Split）： 如果一个区域内部的性质不均匀（比如方差太大，或者包含多种灰度），就把它切成4 个更小的区域（通常基于四叉树结构）。
  • 聚合（Merge）： 切分后，如果相邻的两个小区域性质非常相似，就把它们合并回去。

在具体实现层面，在初始状态时，我们可以把整幅图像看作一个大区域，而我们可以检查这个区域的性质，如果不均匀就分裂成 4 个象限，然后再对每个象限重复这个过程，直到所有区域都均匀为止。随后，我们再检查相邻区域，如果它们非常相似就合并它们。这个过程不断重复，直到无法再分裂也无法再合并为止。这个标准通常可以用简单的区域内方差来衡量，也可以引入一些更复杂的标准如阈值/目标背景等。

聚类分割

聚类的应用当然非常广泛，也包括图像分割方面。

K-means 聚类原理：随机选取 K 个初始聚类中心，然后计算并将每个样本归到离其最近的聚类中心，然后对每个簇把所有簇内样本的坐标均值视为新的中心，反复直到收敛/达到最大迭代次数。

基于频率域的图像增强

二维的频域与傅里叶变换

在一维连续信号中，傅里叶变换和逆傅里叶变换的定义如下：

那么很显然有一维的傅里叶变换对那也有二维的，具体公式如下：

值得一提的是，按照标准的傅里叶变换公式，二维离散傅里叶得到的结果如果画成变换幅度谱，其低频成分在四角，而高频成分在中间，这与我们通常的直觉有点不太一样，因此我们通常会对变换结果进行中心化处理，即将频谱的四个象限进行交换，使得低频成分在中间，高频成分在四角。具体实现时，我们可以通过对图像的每个像素乘以 $(-1)^{x+y}$ 来实现这一点。或者在变换后执行 fftshift 操作。

对二维傅里叶变换，我们关注了三种图谱，也就是其幅度谱，相位谱和功率谱，功率谱定义即为 $P(u, v) = |F(u, v)|^2 = R^2(u, v) + I^2(u, v)$。

与一维傅里叶变换类似，二维傅里叶变换也满足诸多性质：

• 线性性质

• 平移性质

• 旋转性质

• 可分性

• 卷积定理：$f(x, y) * h(x, y) \leftrightarrow F(u, v) \cdot H(u, v)$

频率域图像增强

频率域图像增强的基本思想是通过设计不同的频率域滤波器来增强图像的某些特征。对一些在时域上可能不是很好描述但在频域上比较明显的特征进行增强/减弱。

低通滤波器

理想低通滤波器

画成幅度谱就是一个圆形，圆形内保留，圆形外一律去除。这种滤波器的优点是计算简单，效果明显，缺点是会产生振铃效应（Gibbs 现象）。而且其无法在物理器件上直接实现。

巴特沃斯低通滤波器

其相对降噪效果较温和，其可能仍然会掺杂一些高频噪声，但有效地解决的振铃效应的问题。

高斯低通滤波器

相对于前两者，高斯低通滤波器的频率响应更加平滑，进一步抑制了振铃效应，同时其降噪效果也较为温和，不过相对而言其对低频信号也会有一定的抑制。常用于修复破碎文字，减少皱纹，去除扫描线等工作，直接用于降噪的情况较少。

高通滤波器

• 理想高通滤波器

• 巴特沃斯高通滤波器

• 高斯高通滤波器

• 频域拉普拉斯

同态滤波

同态滤波核心是想解决例如光照等问题。其基本思路是将图像的乘法模型（含光照影响的图=原图 × 光照）转换为加法模型，然后在频域上进行滤波处理，最后再转换回去。具体步骤如下：

1. 对数变换： 对输入图像 $f(x, y)$ 进行对数变换，得到 $g(x, y) = \ln(f(x, y))$。这样，乘法模型 $f(x, y) = i(x, y) \cdot r(x, y)$ 转换为加法模型 $g(x, y) = \ln(i(x, y)) + \ln(r(x, y))$。
2. 傅里叶变换： 对 $g(x, y)$ 进行二维傅里叶变换，得到频域表示 $G(u, v)$。
3. 设计滤波器： 设计一个同态滤波器 $H(u, v)$，其通常是一个高通滤波器，用于增强图像的高频成分（细节）并抑制低频成分（光照变化）。
4. 频域滤波： 将滤波器应用于频域表示，得到滤波后的频域表示 $S(u, v) = H(u, v) \cdot G(u, v)$。
5. 逆傅里叶变换： 对 $S(u, v)$ 进行逆傅里叶变换，得到空间域表示 $s(x, y)$。
6. 指数变换： 对 $s(x, y)$ 进行指数变换，得到最终增强后的图像 $f_{enhanced}(x, y) = e^{s(x, y)}$。

图像复原

图像增强是一个比较宽泛的概念，针对不同的场合和不同的图像，可能有不同的增强方法。因此，其更多也没有一个统一的评判标准，是一个比较主观的操作。而图像复原是一个相对比较客观的概念，且其目的也很直接，就是尽可能提高图像得到保真度。

我们可以将图像的退化/复原视为一个简单的系统，假设原图像为 $f(x, y)$，经过退化函数 $h(x, y)$ 卷积，再加上噪声 $\eta(x, y)$，得到观测图像 $g(x, y) = f(x, y) \ast h(x, y) + \eta(x, y)$，即

而复原的目标就是通过 $g(x, y)$ 来尽可能准确地估计 $f(x, y)$。其关键就在于获得退化函数 $h(x, y)$ 和噪声 $\eta(x, y)$ 。

去噪

常见噪声

常见的噪声包括高斯分布噪声，瑞利分布噪声，伽马分布噪声，均匀分布·噪声，（脉冲）椒盐噪声等。

• 高斯分布噪声：相当常见的一种噪声，来源于电子器件的热噪声等，电路内部相互影响，拍摄现场环境等 $p(z) = \frac{1}{\sigma \sqrt{2\pi}} e^{-\frac{(z - \mu)^2}{2\sigma^2}}$ 其遵守大部分（70%）噪声都落在 $\mu \pm \sigma$ 范围内，绝大部分（95%）噪声落在 $\mu \pm 2\sigma$ 范围内。
• 瑞利分布噪声：通常出现在雷达图像中 $p(z) = \frac{z - a}{b} e^{-\frac{(z - a)^2}{2b}}, z \geq a， 0 \text{ for others}$ $\mu = a + \sqrt{\pi b / 4}, \sigma^2 = b(4 - \pi)/4$
• 伽马分布噪声：通常出现在激光成像的图像中 $p(z) = \frac{a^bz^{b-1}}{(b-1)!}e^{-az}, z \geq a, 0 \text{ for others}$ $\mu = b/a, \sigma^2 = b/a^2$ 当 b=1 时，也称指数分布噪声。
• 均匀分布噪声：所有灰度值出现的概率相同 $p(z) = \frac{1}{b - a}, a \leq z \leq b, 0 \text{ for others}$ $\mu = (a + b)/2, \sigma^2 = (b - a)^2/12$
• 椒盐噪声：图像中随机出现的黑（椒）白（盐）点，表现在图像的短暂停留中，比如错误的开关操作。 $$p(z) = \begin{cases} P_a, & z = a \\ P_b, & z = b \\ 0, & \text{others} \end{cases}$$
• 周期噪声：周期性出现的噪声，通常来源于电磁干扰等，其在频域中可以比较清晰地观察到

只有加性噪声的图像复原比较简单，且往往实际上就是增强。复原的第一步是估计噪声种类，我们可以通过直观经验观察，先验知识获取，，理想图像获取，或者尝试通过观察退化图像获取。 从退化图像获取就是找到图像本身中比较平坦的区域，然后计算这些区域的均值和方差/灰度直方图，从而估计噪声的分布类型和参数。

常见滤波器

均值滤波器

• 算术均值滤波器：最简单的一种，之前我们已经见过了。

• 几何均值滤波器：对乘法噪声效果较好

  $$\hat{f}(x, y) = [\prod_{(s, t) \in S_{xy}} g(s, t)]^{1/|S_{xy}|}$$

  在大部分情况下，其平滑度与算术均值类似，但细节损失更少，但有个致命缺陷是如果有一个像素值为 0，整个窗口输出都是 0。

• 谐波均值滤波器：对盐噪声和高斯噪声效果，但不适用于椒噪声

  $$\hat{f}(x, y) = \frac{S_{xy}}{\sum_{(s, t) \in S_{xy}} \frac{1}{g(s, t)}}$$

• 逆谐波均值滤波器：当 Q 为正值时看一看消除胡椒噪声，当 Q 为负值时可以消除盐噪声

  $$\hat{f}(x, y) = \frac{\sum_{(s, t) \in S_{xy}} g(s, t)^{Q+1}}{\sum_{(s, t) \in S_{xy}} g(s, t)^Q}$$

  当 Q=0 时就是算术均值滤波器，当 Q=-1 时就是谐波均值滤波器。

统计排序滤波器

• 中值滤波器：对椒盐噪声效果非常好

  $$\hat{f}(x, y) = \text{median} \{g(s, t) | (s, t) \in S_{xy}\}$$

• 最大/小值滤波器：

  $$\hat{f}(x, y) = \max \{g(s, t) | (s, t) \in S_{xy}\}, \text{max filter} \\ \hat{f}(x, y) = \min \{g(s, t) | (s, t) \in S_{xy}\}, \text{min filter}$$

• 中点滤波器：对均匀/高斯噪声效果较好：

  $$\hat{f}(x, y) = \frac{1}{2} [\max \{g(s, t) | (s, t) \in S_{xy}\} + \min \{g(s, t) | (s, t) \in S_{xy}\}]$$

• 修正 α 均值滤波器：均值和中置混合滤波器，对椒盐噪声和高斯噪声均有较好效果

  $$\hat{f}(x, y) = \frac{1}{|S_{xy}| - 2\alpha} \sum_{(s, t) \in S_{xy}} g_r(s, t)$$

  其中 g_r(s, t)是去掉了窗口中最大和最小的 α 个像素值后的剩余像素值。

• 自适应中值滤波器：其通过两步来自适应调整窗口中心的像素

  • 它会先在一个小窗口（比如 3×3）里计算最大值、最小值和中值。
  • 判断中值是否可靠：它首先检查中值本身是不是噪声（即中值是否在最大值和最小值之间）。如果中值不是噪声，那么就继续下一步；否则，增大窗口尺寸，重新计算，直到找到一个可靠的中值或者达到设定的最大窗口限制。
  • 然后检查窗口中心的像素值是否是噪声：如果窗口中心的像素值在最大值和最小值之间，那么它就是一个正常的像素，保持不变；否则，将其替换为之前计算得到的可靠中值。

除了常规的时域滤波以外，正如上面所言，部分周期性噪声也可以用频率滤波消除。

退化

回忆图像退化的表达式 $g(x, y) = f(x, y) \ast h(x, y) + \eta(x, y)$

$G(u, v) = F(u, v) \cdot H(u, v) + N(u, v)$

上述我们只考虑了加性噪声也即是 $\eta(x, y)$，默认没有退化函数即为 1，但实际上图像在采集过程中往往会受到各种退化因素的影响，例如运动模糊，散焦模糊等，这些都可以通过退化函数 $h(x, y)$ 来描述。

估算退化函数有几种常见方法：

• 图像观察估计法：如果能手动还原一部分图片，我们可以直接同通过原子图像$G_s(u, v)$和重建后的子图像$F_s(u, v)$得到$H(u, v) = \frac{G_s(u, v)}{F_s(u, v)}$

• 实验估计法：通过实验手段获取退化函数，例如拍摄一个已知的图像（最简单的例如模拟冲激，其会在频域上显示为常数），然后通过分析拍摄结果来估计退化函数。

• 模型估计法：部分退化函数是由明确的物理过程产生的，例如运动模糊，我们可以通过分析这些物理过程来建立数学模型，从而得到退化函数。

常见退化模型

• 运动模糊模型：假设图像在曝光期间以恒定速度沿直线运动，则实际图像 f(x, y)与观测图像 g(x, y)之间的关系可以表示为：

  $g(x, y) = \int_{0}^{T} f(x - at, y - bt) dt + \eta(x, y)$

  根据傅里叶变换的相关性质，

  $$G(u, v) = F(u, v) \int_0^T e^{-j2\pi(ua + vb)t} dt$$

  退化函数可以表示为：

  $$H(u, v) = \int_0^T e^{-j2\pi(ua + vb)t} dt$$

  继而化简为

  $$H(u, v) = \frac{T}{\pi(ua + vb)} \sin[\pi(ua + vb)] e^{-j\pi(ua + vb)}$$

  其中 $a$ 和 $b$ 是运动的速度分量，$T$ 是曝光时间。

• 大气湍流模型：大气湍流会引起图像的模糊和失真，其退化函数可以表示为：

  $$H(u, v) = e^{-k(u^2 + v^2)^{5/6}}$$

  其中 $k$ 是与大气条件相关的常数，k 值越高，湍流越剧烈，模糊也就越严重。

复原

在求得退化函数和噪声模型后，我们就可以进行图像复原了。常见的复原方法包括逆滤波方法和维纳滤波方法。

逆滤波

即 $\hat{F} = \frac{G(u, v)}{H(u, v)} = F(u, v) + \frac{N(u, v)}{H(u, v)}$ 也就是说这里噪声频谱也会被 $H(u, v)$ 放大，因此当 $H(u, v)$ 很小时，噪声会被放大，从而导致复原图像质量降低，如果为 0 的时候更会直接导致运算错误，所以一般我们可以在分母中加入一个小的常数 k，即原式可被修改为。

维纳滤波

标准维纳滤波的频域表达式如下：

其中 $\frac{S_\eta(u, v)}{S_f(u, v)}$ 表示噪声与原图的功率谱密度比值。不过很显然大部分实际情况下我们的信噪比并不是已知的，因此我们可以将其视为一个常数 $K$，从而得到简化的维纳滤波器：

$\hat{F}(u, v) = [\frac{H^*(u, v)}{|H(u, v)|^2 + K}] G(u, v)$

其中 K 是一个可调参数

复原的质量评价

有很多评价方法，如平均绝对误差（MAE），均方误差（MSE），归一化均方误差（NMSE），峰值信噪比（PSNR）等，信噪比（SNR）等。