DIP Crash Course 数字图像处理速成笔记
将所有的开放问题的 solution 当作一个普通的解决方法来看待, 它们绝不是最优或者唯一的解决方案.
每个变换背后尽可能想象一个图像的画面变化.
1 Basic Conventions 基本约定
Image 默认为 sampling 和 quantization 之后的数字图像
- Spacial resolution 空间分辨率: 长宽方向的采样点数的乘积.
- 规定图片的左上角像素坐标为 \((0,0)\).
- Intensity resolution 强度分辨率: 量化的精度, 一般用每个像素的 bit 数表示. 默认为 \(8\) bit (\(\mathbb{Z}_{256}\), 或者不严格地, \([0, 255]\)).
- Spacial resolution 空间分辨率: 长宽方向的采样点数的乘积.
Color space 颜色空间: 单色可见光的波长在 \([380, 780] \text{ nm}\), 每一束光都可以有任意强度的单色光分量 (即「光谱」), 因此所有光构成一个不可数无穷维拓扑向量空间.
- 但是人眼有三种类型的视锥细胞 (L, M, S) 对不同的频率成分有 (不会变化的) 不同的敏感度, 它们自动将所有光线性映射到一个三维子空间 \(\Omega\) 中 (想想为什么是线性的, 因为是积分). 即使两个光的光谱不同, 但它们在这个三维子空间的投影相同, 人眼也无法区分它们.
- 不同的应用场景在 \(\Omega\) 中选择不同的基底:
- RGB (Red, Green, Blue): 显示器, 图像储存.
- CMY/CMYK (Cyan, Magenta, Yellow, K: Black ink): 打印机, K 成分是多余的, 就当作是为了节省墨水 lol.
- HSI/HSV (Hue, Saturation, Intensity/Value): 图像编辑软件.
默认情况下图片丢失拓扑结构
view()成一个 vector. \(r\) 代表处理前的图像 vector, \(s\) 代表处理后的图像 vector.RGB2Gray 转换方法: 可以任选一通道 / 取最大 / 平均值 / 加权平均.
Contrast 对比度: 灰度图像中像素值最大值与最小值之差 (去掉 outlier).
Binary Image 二值图像 \(f\): \(0\) 表示 background, \(1\) 表示 foreground.
1.1 Notations
- 图像空间 \(\text{Img}\): 所有不同大小的灰度图像构成的 Vector Space.
- 这是最大的图像空间, 只要是图像都在里面.
- 彩色图像空间记为 \(\text{\textcolor{red}{I}\textcolor{green}{m}\textcolor{blue}{g}}\).
2 Kernel-based processing 基于核的处理
本人没有按照 Spatial domain 和 Frequency domain 来划分, 因为他们本质上是等价的, 无法将二者严格区分.
Kernel-based method 就是将原图片与一个 kernel 做滑动内积 (纠结是内积还是卷积毫无意义). 重点在于 kernel 的设计, 我们通过 kernel 设计的思维路径将其划分为 Spatial filtering 和 Frequency filtering 两大类.
2.1 Spatial Filtering 空间滤波
Spatial filtering 指 kernel 的设计过程没有经过频域分析.
Padding 填充
NCC (Normalized Cross-Correlation): 图像 vector 之间的夹角余弦值.
Common filters 常用滤波器
- Smoothing filters 平滑滤波器 (低通滤波器)
- Order-statistic filters 次序统计滤波器: 将某框内的像素排序后取中位数, 一般可去除 salt-and-pepper (impulse) noise.
- Sharpening (derivative) filters 锐化滤波器 (高通滤波器): 将梯度、Laplacian 算子移植到图片上; 由于这些算子是线性的, 可等价于与一个 kernel 做卷积; 这些 kernel 也可以在形式上被强化来达到更强的锐化效果.
3 Non-kernel-based processing 非基于核的处理
Non-kernel-based method 指处理过程无法被等价为与一个确定的 kernel 做内积的过程.
3.1 Intensity Transformation
强度变换仅仅是单灰度像素的操作, 即用一个函数 \(s = T(r)\) 来独立地映射所有像素值, 一般情况下
- Monoid \(\mathcal{T}\): 一张图片 \(r\) 上的所有 Intensity Transformation 构成幺半群 \(\mathcal{T}\).
- Group \(\mathcal{T}_{\text{bij}}\): 从 \(\mathcal{T}\) 中取出所有 bijective 映射, 它构成一个群. 其中包括:
- Image Negatives 图像反相: \(s = L - 1 - r\).
- Group \(\mathcal{T}_\uparrow\): \(\mathbb{Z}_{256}\) 的保序同胚群 (所有保序映射 (严格单增) 的 \(T\)), 是 \(\mathcal{T}_{\text{bij}}\) 的子群. 引入此的意义是一般图像处理都只希望深颜色的像素在处理后还是比较深的颜色, 而不会变白了, 所以实际上几乎只会考虑 \(\mathcal{T}_\uparrow\) 中的变换. 常用的有:
- Contrast stretching 对比度拉伸
- Log transformation 对数变换: \(s = c \cdot \log(1 + r)\).
- Exponential transformation 指数变换
- Power-Law (Gamma) transformation 幂律变换: \(s = c \cdot r^{\gamma}\).
- Group \(\mathcal{T}_{\text{bij}}\): 从 \(\mathcal{T}\) 中取出所有 bijective 映射, 它构成一个群. 其中包括:
3.2 Histogram Processing 直方图处理
直接想出一个 Intensity Transformation 比较难, 我们引入分析灰度图像的另一种方法: 直方图 (Histogram). 注意这只是辅助我们设计 \(T\) 的中间步骤, 我们会看到用 Histogram 可以快速、有根据地帮我们得出一个 \(T\).
- Histogram 图像直方图: 每个图像以 Figure 2 的方式给 \(\mathbb{Z}_{256}\) 数轴引入了权重, 即 measure \(\#\). Histogram 定义为测度空间 \((\mathbb{Z}_{256}, \mathcal{P}(\mathbb{Z}_{256}), \#)\).
- 可以引入 \(\text{cdf}(\#)\).
为后文表达方便我们引入这些冗余概念:
- \(y\) 归一化直方图 \((\mathbb{Z}_{256}, \mathcal{P}(\mathbb{Z}_{256}), \mathbb{P})\): 纵坐标除以了像素总数 (Figure 3), 相当于概率.
- 可以引入 \(\text{cdf}(\mathbb{P})\).
- \(x\) 归一化直方图 \(([0,1], \mathcal{B}([0,1]), \#')\): 在归一化直方图的基础上横坐标也归一化到 \([0,1]\) 上 (见 Figure 2).
- 可以引入 \(\text{cdf}(\#')\).
- \(x,y\) 归一化直方图 \(([0,1], \mathcal{B}([0,1]), \mu)\): 横纵坐标都归一化到 \([0,1]\) 上 (见 Figure 3).
- 可以引入 \(\text{cdf}(\mu)\).
下文都默认在 \(x,y\) 归一化直方图的测度空间 \(([0,1], \mathcal{B}([0,1]), \mu)\) 上讨论! 以防到处出现恶心的 \(255\).1
1 PPT 上就没有用 \(x,y\) 归一化的方式表达, 非常 annoying!
3.2.1 Histogram Equalization 直方图均衡化
Motivation: 如果图像整体看起来比较一致, 对比度不高, 我们希望不同深浅的像素都要有一些 (尽量均匀一点), 这样图像看起来更清晰.
- 目标: 设计一个 \(T\) 使得输出图像的 histogram 均匀分布在 \([0,1]\) 上.
- 即: Find \(T: [0,1] \to [0,1]\) s.t., \(\mathbb{\mu}\) 的 pushforward 测度 \(T_* \mu\) 是 uniform 的 (Lebesgue 测度), i.e., \[T_* \mu = \lambda.\]
- Solution: Surprisingly!!! 这个 \(T\) 就是 \(\mu\) 的累积分布函数 (c.d.f.)! (见 Lemma 1) 仔细想想看, 这正是导数的定义!!
3.2.2 Histogram Matching 直方图匹配
Motivation: 我想让这张图的整体明暗感觉看起来像那张图, 请设计一个 \(T\) 来达到这个目的.
- 目标: 刚好是 Histogram Equalization 的 generalization, 即设计一个 \(T\) 使得输出图像的 histogram 是一个指定的分布 \(\nu\).
- 即: Find \(T: [0,1] \to [0,1]\) s.t., \(\mu\) 的 pushforward 测度 \(T_* \mu = \nu.\)
- Solution: 很显然, 我们可以得到一个均匀分布 \(\lambda\), 再用这个均匀分布 #lem-equalization 得到我们想要的分布 \(\nu\).
4 Image Restoration 图像复原
图像就是在已知图像是如何损坏 (Degradation model)的情况下, 修复图像的过程.
4.1 Degradation/Restoration Model 损坏/复原模型
显然有:
\[ \begin{aligned} r &= f * h + n, &\hat{f} = r * k \\ R &= F \cdot H + N, &\hat{F} = R \cdot K \end{aligned} \]
4.1.1 只有噪声 \(n(x,y)\)
- Noise Models 噪声模型
- Gaussian noise 高斯噪声
- Rayleigh noise 瑞利噪声: 一般出现在雷达图像中.
- Uniform noise 均匀噪声
- Salt-and-pepper (impulse) noise 椒盐噪声
- Periodic noise 周期噪声: 一般是收到了电磁干扰 (比如宇宙图像).
4.1.2 只有模糊 \(h(x,y)\)
或者噪声非常小可以忽略.
- Inverse filtering 逆滤波: 直接用 \[K=\frac{1}{H}\] 来复原图像.
- Limitation 局限性: 当 \(H\) 在某些频率上接近 \(0\) 时, 逆滤波会极大地放大噪声, 导致复原效果很差! 比如下面的 EXAMPLE:
4.1.3 同时有模糊 \(H\) 和噪声 \(n\)
5 Image Segmentation 图像分割
一个非常自然的需求是把一张图像分割成若干个有意义的区域. 由于「有意义」无法准确定义, 因此图像分割问题甚至没有标准答案.
- Neighbors of a pixel 二值图像像素的邻域: 有不同的定义! 4-neighbors (上下左右), 8-neighbors, diagonal-neighbors (对角线).
- 领域不一定是等价关系!
- Path 路径: 一个像素的有序序列, 后一个像素都是前一个像素的 neighbor.
- Connected set 连通集: 当邻域是等价关系 \(\sim\) 时, 二值图像 \(b/\sim\) 中的每个等价类就是一个连通集.
- Boundary 边界: 跟拓扑学里的定义一样, 边界点的邻域包含区域外的像素.
5.1 Edge Detection 边缘检测
5.1.1 空间一阶导
梯度大的地方很可能就是「边缘」, 因此很自然的想法是计算图片每个点处的梯度 (先各算 \(x,y\) 方向, 然后用平方加起来根号), 然后把梯度大于某个阈值的点标记为边缘. 而对于离散图像, 梯度有 Central difference (\(\delta(x,y)\)), Forward difference \(\Delta(x,y)\), Backward difference \(\nabla(x,y)\) [5]. 由于它们都是线性运算, 因此可以被等价为与一个 derivative kernel 做相关 (“[ ]” 表示相关出来的结果填在该位置):
Central Forward Backward \(x\) direction kernel \(\begin{pmatrix}-1 & [0] & 1\end{pmatrix}\) \(\begin{pmatrix}[-1] & 1\end{pmatrix}\) \(\begin{pmatrix}-1 & [1]\end{pmatrix}\) \(y\) direction kernel \(\begin{pmatrix}-1 \\ [0] \\ 1\end{pmatrix}\) \(\begin{pmatrix}[-1] \\ 1\end{pmatrix}\) \(\begin{pmatrix}-1 \\ [1]\end{pmatrix}\) 当然「用梯度来检测」是一个拍脑袋的想法, 在这个基础上我们可以将梯度算子一般化, 只要能检测出两侧的差异就行呗, 于是我们得到了 (当然完全不限于这些, 你可以自己设计):
Sobel Prewitt Roberts \(x\) direction kernel \(\begin{pmatrix}-1 & 0 & 1 \\ -2 & [0] & 2 \\ -1 & 0 & 1\end{pmatrix}\) \(\begin{pmatrix}-1 & 0 & 1 \\ -1 & [0] & 1 \\ -1 & 0 & 1\end{pmatrix}\) \(\begin{pmatrix}[0] & 1 \\ -1 & 0\end{pmatrix}\) (\(\nearrow\) direction) \(y\) direction kernel \(\begin{pmatrix}-1 & -2 & -1 \\ 0 & [0] & 0 \\ 1 & 2 & 1\end{pmatrix}\) \(\begin{pmatrix}-1 & -1 & -1 \\ 0 & [0] & 0 \\ 1 & 1 & 1\end{pmatrix}\) \(\begin{pmatrix}[1] & 0 \\ 0 & -1\end{pmatrix}\) (\(\nwarrow\) direction) 最后我们设一个阈值, 得到一个二值图像:
Figure 10: 以 Roberts operator 为例, 设置不同阈值后的效果 [6].
5.1.2 空间二阶导
上面的各种算子已经够 generalized 的了, 但是它们的本质还是 “first-order-like”, 也许空间的二阶导也会给我们一些关于边缘的信息? 我们将会设计一些 “second-order-like” 的 kernel.
空间二阶导就是 Laplacian operator \[\nabla^2 f = \frac{\partial^2 f}{\partial x^2} + \frac{\partial^2 f}{\partial y^2}.\] 这些偏导数也有 Central, Forward, Backward difference 的离散形式, 而且两次偏导数用的还可以不一样 (比如 Central + Forward), 但是下面只考虑两次用的都一样. 由于是直接将 \(x,y\) 方向上的数值相加, 只需要一个 kernel 就行了 (实际情况一般就用 Central, 因为另外两个边缘的位置会偏):
Central Forward Backward Laplacian kernel \(\begin{pmatrix}0 & -1 & 0 \\ -1 & [4] & -1 \\ 0 & -1 & 0\end{pmatrix}\) \(\begin{pmatrix}[2] & -2 & 1 \\ -2 & 0 & 0 \\ 1 & 0 & 0\end{pmatrix}\) \(\begin{pmatrix}0 & 0 & 1 \\ 0 & 0 & -2 \\ 1 & -2 & [2]\end{pmatrix}\)
Figure 12: 直接用 Laplacian \(=0\) 会发现全是 edge. 于是我们先用一个低通滤波器 (比如对原图进行高斯模糊), 再用 Laplacian (Laplacian of Gaussian, LoG) [6].
5.1.3 Canny Edge Detector
Canny edge detector 是一个在 CV 中广泛使用的方法, 他综合了前面两者的优点. 当然它也很拍脑袋, 你可以想象成是你一个很聪明的同学想出来的.
- Canny 的步骤很简单:
- Gaussian: 首先高斯滤波.
- Compute Gradient: 然后计算每个像素处梯度大小和其方向.
- Thresholding: 选出梯度大小在一个预设的范围内的像素点. 假设 Figure 13 中的蓝色像素被选出了.
- Directional second-derivative: 在所有被选出的像素上沿着梯度方向做一个一维的二阶导数, 找到其零点并在原图上标记出来. 比如 Figure 13 中蓝色像素的梯度方向是右上方, 橙色像素的二阶导值最接近 \(0\), 保留橙色像素, 丢弃蓝色像素 (想一想, 这样会让边缘更细、更精确). 对所有蓝色的像素都做这个操作, 最终得到 Figure 13 中的所有橙色像素, 这些就是 Canny 检测到的边缘.
Figure 13: Canny 步骤 [7]. - Edge linking 边缘连接: 边缘可能不连续, 可以利用梯度方向等信息把边缘连接起来.
5.2 Edge-based object detection 基于边缘的目标检测
5.2.1 Hough Transform
Hough Transform 目标不再仅仅是找到边缘, 而是基于 Edge Detector 的二值边缘图像标记出特定形状 (比如直线, 圆) 的位置以及它们的数量2.
2 我觉得 Hough Transform 的思想非常优美, 我在想如果直接用机器学习的方法训练一个网络来找直线位置以及数量, 是否在神经网络的某个角落里就隐藏着 Hough Transform 这个算法呢? 怎么找到呢? 或者更一般地, 如何提取神经网络学到的东西?
看一遍 This video 就懂了.
5.2.2 Distance Transform (DT) 距离变换
DT 的目标也是基于 Edge Detector 的二值边缘图像, 再给定一个模版的二值边缘图像, 目标是在原图中找到与模版边缘最相似的区域. 一个很自然的想法是直接计算二值模版图像在原二值图像上的归一化内积 (NCC), 然后找到最大的位置. 但是模版可能会有缩放, 导致匹配的结果很差. 因此我们引入 DT.
- 操作步骤:
5.3 Non-edge-based segmentation 非基于边缘的分割
这种分割图像的方法不从「变化」来先找边界, 而是从「相似性」出发直接找区域 (Region-based segmentation).
5.3.1 Thresholding
根据阈值是否随空间变化可分为 Global thresholding 和 Variable thresholding. 注意这两种方法都可以涉及阈值迭代的过程, 即有一种机制不断更新 threshold 直到收敛.
- Global thresholding:
- 对光照变化大的图像效果不好.
- Variable thresholding:
- 可以自动适应光照变化.
5.3.2 Region Growing
- 给定待分割的灰度图像 \(r\), 给定一个 seed binary mask \(m\), 计算 \(r\) 在 \(m\) 覆盖的区域的均值 \(\mu\), 然后设置一个 threshold \(T\), 把 \(r\) 中所有在 \([\mu - T, \mu + T]\) 范围内且与 \(m\) 相邻的像素点加入 \(m\), 重复这个过程即可.
- 当然算均值 + threshold 也可以换成别的相似性度量方法 (similarity criteria).
5.3.3 K-Means Clustering
- 假设我们想对一张图片进行基于颜色的 segmentation, 我们可以把每个像素的颜色 (比如 RGB) 看作一个三维空间中的点, 然后在 RGB 空间中进行 K-Means Clustering 来把这些点分成 \(K\) 类:
- 随机选择 \(K\) 个初始中心点 (centroids).
- 将每个点分配给距离它最近的中心点所属的类. (注意每个点都会被分配而不是留到以后再分).
- 重新计算每个类的中心点 (centroids).
- 重复上面两步直到收敛 (中心点不再变化, 或其变化小于某个阈值 \(\varepsilon\)).
- K-means 也可以用在别的 feature space 上! 比如它甚至可以用于将图像本身分类. 假设我们有很多图像, 每张图像在某个 feature space 中都是一个点. 我们可以在这个 feature space 上用 K-means 将这些图像分成 \(K\) 类.
6 Morphological 形态学图像处理
这里主要涉及 PS 软件中的各种基础操作的原理.
6.1 Binary erosion / dilation 二值图像腐蚀/膨胀
我们希望 foreground binary image 里面的小点消失 (Dilation)、快连上的线能连上 (Erosion).
- 想象一个 binary image \(f\) 的前景是凹下去的模具, 想象 SE (Structuring Element) 像俄罗斯方块一样掉落在 \(f\) 上, 记录下:
- Erosion: SE 能整体掉进去的地方, 这些就是腐蚀 \(f \ominus se\) 的结果. 腐蚀过的前景 (\(f \ominus se\)) 只会比 \(f\) 一样大或者更小.
- Dilation: 任何一个 SE 方块能掉进去的地方, 这些就是膨胀 \(f \oplus se\) 的结果. 膨胀过的前景 (\(f \oplus se\)) 只会比 \(f\) 一样大或者更大.
- Erosion 和 Dilation 互为 dual.
- 我们可以通过给原二值图像的 DT 图像加 threshold 来计算腐蚀和膨胀 (见 Figure 18).
6.2 Opening / Closing 开运算/闭运算
Erosion / Dilation 虽然比较有效, 但是有一个副作用就是前景的粗细都会顺带一起被改变: 小点消失的同时整体的线都变细了, 即将要连上的线连起来了但整体的线也变粗了. 为了解决这个问题, 我们引入 Opening / Closing.
- Opening 开运算: 先腐蚀后膨胀, \(b \circ se \equiv (b \ominus se) \oplus se\).
- 消失的小点不会因为膨胀而再出现.
- Closing 闭运算: 先膨胀后腐蚀, \(b \bullet se \equiv (b \oplus se) \ominus se\).
- 连上的线一般不会因为腐蚀而再断开.
7 PCA 主成分分析
8 Image Compression 图像压缩
信息永远不能被压缩, 只是被压缩的部分变成编解码规则了.
8.1 Basic Concepts 基本概念
- 若图片 \(f\) 压缩前可以被 \(100\) bit 表示, 压缩后可以被 \(75\) bit 表示, 则:
- Compression ratio 压缩比 \(C = 100/75\).
- Relative data redundancy 相对数据冗余 \(R = 1-0.75\).
8.2 Lossless Compression 无损压缩
- Huffman coding 霍夫曼编码
- Huffman Tree 的构建: 将所有符号的频次 (或概率) 从小到达排序, 选两个频次最小的创建他们的公共父节点, 频次为两者之和, 然后将这个新节点插回列表中重复步骤, 直到只剩一个节点为止.
- Arithmetic coding 算术编码
- LZW 编码
8.3 Lossy Compression 有损压缩
- JPEG:
- 颜色空间变换 (RGB \(\to\) YCbCr): 利用人眼对亮度更敏感的特点.
- 分块 (\(8\times 8\) blocks): 计算较快.
- DCT (离散余弦变换): 就是 DFT.
- Quantization (量化,唯一真正「有损」的一步): 在频域被选择性地 mask 掉一些系数.
- Zig-zag 扫描
- Entropy coding (Huffman / Arithmetic)
8.4 Watermarking 数字水印
Visible watermark 可见水印: 直接把水印叠加在图片上.
Invisible watermark 隐形水印
- Robustness 水印必须能抵抗以下操作: 压缩, 缩小/放大, 裁剪, 加噪声, 滤波, 镜像, etc.
- 简单来说, 可以将水印直接叠加在频域系数上, 再反变换得到水印图像.