数字图像处理

目标：对图像进行相对底层的基本操作

Introduction

图像采样与量化：

• 采样：图像空间连续坐标的离散化，决定图像的空间分辨率。
• 量化：图像函数值（幅值）的数字化，决定图像的幅度（灰度级）分辨率。

图像表示

图像文件格式：

• Vector images(.ai,.eps,.ps,…)：矢量图，放缩时不会走样和模糊，但难以获得。
• Bitmap(.bmp,.jpg,.png,.gif,…)：位图，易于获得，但放缩时会走样和模糊。
  • GIF(Graphics Interchange Format)
    • 8位索引，最多可保存256种颜色
    • 具有"抖动(dither)"选项：可以混合两种不同可用颜色的像素来创建另一种颜色的像素
    • 可以动画；透明
    • 小内存、高品质
  • JPEG(Joint Photographic Experts Group)
    • 16位，能够同时显示数百万种颜色而不会抖动
    • 约60%的压缩设置实现质量和文件大小的最佳平衡
  • PNG(Portable Network Graphics)
    • 基于ZIP的无损压缩
    • 可以透明（4通道图像）
  • BMP(Windows Bitmap)
    • 简单未压缩
    • 可以索引或不索引
    • DIB(Device Independent Bitmap)/DDB(Device Dependent Bitmap)
• 矢量图通过光栅化（rasterize）得到位图

图像（位图）可以表示成在二维域上定义的函数 $f(x, y)$，每一个元素称为像素（pixel）

图像在内存中的表示：

本部分的其余内容可在“计算机视觉”的课程复习博客中的Introduction部分查看

Fundamental

本部分的内容可在“计算机视觉”的课程复习博客中的Image Processing部分查看

Spatial Filtering

本部分的内容可在“计算机视觉”的课程复习博客中的Spatial Filtering部分查看

Statistics

图像直方图是表示图像中灰度分布的直方图，描述了各个灰度上像素的数量关系。

以灰度级为x轴，像素数量为y轴，较暗的图像像素会集中在直方图的左侧；较亮的图像像素会集中在直方图的右侧；对比度低的图像像素会挤在直方图的中间；对比度高的图像像素会分散在直方图的各个灰度级上。

直方图均衡化
将灰度级 $r$ 归一化到 $[0, 1]$
灰度变换函数：$s_k=T(r_k)=\sum_{j=0}^k\frac{n_j}{n}=\sum_{j=0}^kp_r(r_j)$，$n$为总像素数，$n_j$为灰度级为$j$的像素数，$p(r)$为概率密度函数。该公式表示将灰度为$r_k$的像素变换到灰度$s_k$上

全黑(0)图像直方图均衡化后为全白(255)图像；全白(255)图像直方图均衡化后还为全白(255)图像；一半黑一半白的图像直方图均衡化后，黑变灰(127.5)，白不变

直方图应用

• 图像检索：使用每幅图像的直方图作为索引
• 目标跟踪：搜索图像，找到与前一帧目标区域的直方图最匹配的区域
• 图像分割：利用直方图模拟图像的前景和背景的颜色分布，分割出前景和背景

高斯混合和直方图：

• 都可以用于像素亮度/颜色分布的描述
• 优缺点（以交互式图像分割为例）
  • 直方图对交互分割
    • 优点：计算速度快，易于实现，能够捕捉到数据的整体分布。
    • 缺点：无法捕捉到数据的局部特征，对异常值敏感，可能会出现数据过拟合现象。
  • 高斯混合模型对交互分割
    • 优点：可以捕捉到数据的局部特征，对异常值不敏感，具有较高的准确率。
    • 缺点：计算量较大，需要手动确定高斯分布的个数，对初始化参数敏感，可能会出现局部最优解。

Structure

基于图的图像处理方法

图像可以被看作是一张图，每个像素为一个结点，像素之间根据邻接关系和颜色值等以边相连

• 连通域：同一连通域的任意两结点都能以路径相连
  • 区域增长/种子填充：
  • 快速连通域算法：一次扫描+合并等价类
• 最短路：
  • 单源最短路：所有结点到某个结点的最短路
    • Dijkstra
  • 多源最短路：所有结点到某几个结点的最短路
    • 再设置一个节点，连接该节点和几个目标结点，边权为0，求该节点的单源最短路
  • 定义相邻像素之间的距离为其颜色差的函数，颜色差越大，距离越小：$d_{ij}\propto e^{-\beta\parallel I_i-I_j\parallel^2}$
• 图像距离场：所有像素到种子像素的最短距离
  • 距离变换：
    • 
    • 
• 图割和最小割
  • 边的容量：每一条边允许通过的最大流量
  • 图的切割：能将源和目标之间所有路径切断的图的剖分
  • 最小割：所有图的切割中所切断的边的容量之和最小的一个，对应于网络流的瓶劲，即最大流
  • 应用：图像分割，根据像素颜色建图，将源点和目标点分别设置在前景和背景上，通过求解最小割来得到最佳分割

形态学与二值图结构

形态学图像处理

• 膨胀：使图像扩大
  • $A$ 和 $B$ 是 $Z^2$ 中的集合，$A$ 被 $B$ 膨胀定义为：$A\oplus B=\{z|(\hat B)_z\cap A\neq \varnothing\}$ 上式表示 $B$ 的反射进行平移与 $A$ 的交集不为空。$B$ 的反射：相对于自身原点的映象，$B$ 的平移：对 $B$ 的反射进行位移。也可以写成$A\oplus B=\{z|(\hat B)_z\cap A\subseteq A\}$
• 腐蚀：使图像缩小
  • $A$ 和 $B$ 是 $Z^2$ 中的集合，$A$ 被 $B$ 腐蚀定义为：$A\ominus B=\{z|(B)_z\subseteq A\}$ 上式表示 $B$ 进行平移后包含于 $A$
• 开操作：在不改变形状的前提下，使图像的轮廓变得光滑；断开狭窄的间断；消除细的突出物
  • 使用结构元素 $B$ 对集合 $A$ 进行开操作，定义为$A\circ B=(A\ominus B)\oplus B$ 先用 $B$ 对 $A$ 腐蚀，然后用 $B$ 对结果膨胀，另一个定义形：$A\circ B=\cup\{(B)_z|(B)_z\subseteq A\}$
  • 几何解释：$B$ 在 $A$ 的边界内转动时，$B$ 中的点所能达到的 $A$ 的边界的最远点
  • 性质
    • $A\circ B\subseteq A$
    • 若 $C\subseteq D$，则 $C\circ B\subseteq D\circ B$
    • $(A\circ B)\circ B=A\circ B$
• 闭操作：在不明显改变面积前提下，使图像的轮廓变得光滑；消除小的孔洞；消除狭窄的间断；细长的鸿沟；填补轮廓线中的裂痕
  • 使用结构元素 $B$ 对集合 $A$ 进行闭操作，定义为：$A\bullet B=(A\oplus B)\ominus B$ 先用 $B$ 对 $A$ 膨胀，然后用 $B$ 对结果腐蚀
  • 几何解释：$B$ 在 $A$ 的边界外侧转动时，$B$ 中的点所能达到的离 $A$ 的边界的最近点
  • 性质
    • $A\subseteq A\bullet B$
    • 若 $C\subseteq D$，则 $C\bullet B\subseteq D\bullet B$
    • $(A\bullet B)\bullet B=A\bullet B$

Frequency

频率(Frequency)即信号进行周期性变化的速率
图像的频率指图像亮度/颜色在水平/垂直方向上周期性变化的速率

傅里叶变换

图像从空间域到频率域的转换：

• 确定“某种频率”：选择一组具有特定频率的信号（信号的基），且通过这组信号的组合可以表示其它任何的信号
• 计算“频率的成分”：对任意的输入信号，以及某给定的基信号（频率已知），计算该基信号所占的成分（系数）

一维连续傅里叶变换：$F(u)=\int_{x=-\infty}^{\infty}f(x)e^{-j2\pi ux}\,{\rm d}x$
一维离散傅里叶变换：$F(u)=\frac{1}{M}\sum_{x=0}^{M-1}f(x)e^{-j2\pi ux/M},u=0,1,2,\dots,M-1$
以上公式中，$F(u)$表示输入信号占基信号的成分，$f(x)$为输入信号，$e^{-j2\pi ux}$为基信号。

$F(u)$为什么表示输入信号占基信号的成分？
由欧拉公式 $e^{j\theta}=\cos\theta +j\sin\theta$，可得 $F(u)=\frac{1}{M}\sum_{x=0}^{M-1}f(x)(\cos2\pi ux/M -j\sin2\pi ux/M)$
先暂时忽略虚数部分，得：$F^{real}(u)=\frac{1}{M}\sum_{x=0}^{M-1}f(x)\cos2\pi ux/M$

令$\vec{f(x)}=[f(0),f(1),f(2),\dots ,f(M-1)]$，$\vec{\cos u\delta x}=[1,\cos\delta x,\cos 2\delta x,\dots,\cos(M-1)\delta x],\delta=\frac{2\pi}{M}$
则$F^{real}(u)=\frac{1}{M}\vec{f(x)}\cdot\vec{\cos u\delta x}$
$F(u)$是原信号与基信号的内积，$F(u)$与原信号在基信号上投影的长度成正比，是原信号与基信号相似性的描述。傅里叶变换将原信号分解为基信号的组合，$F(u)$为组合的系数。

用这组基信号能否表示所有的其它信号？
傅里叶变换是矩阵（线性）变换
$F^{real}(u)=\frac{1}{M}\sum_{x=0}^{M-1}f(x)\cos2\pi ux/M,u=0,1,2,\dots,M-1$

$\frac{1}{M} \begin{bmatrix} 1 & 1 & 1 & \cdots & 1 \\ 1 & \cos\delta x & \cos 2\delta x & \cdots & \cos(M-1)\delta x \\ 1 & \cos 2\delta x & \cos 4\delta x & \cdots & \cos 2(M-1)\delta x \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \cos(M-1)\delta x & \cos 2(M-1)\delta x & \cdots & \cos (M-1)^2\delta x \\ \end{bmatrix} \begin{bmatrix} f(0) \\ f(1) \\ f(2) \\ \vdots \\ f(M-1) \\ \end{bmatrix} = \begin{bmatrix} F(0) \\ F(1) \\ F(2) \\ \vdots \\ F(M-1) \\ \end{bmatrix}$

$\because 1,\cos x,\sin x,\cos 2x,\sin 2x,\dots,\cos nx,\sin nx,\dots$ 任意两个不同函数之积在$[0,2\pi]$内积分为 $0$
$\therefore$ 它们在$[0,2\pi]$内相互正交
$\therefore$ 它们是$[0,2\pi]$区间所有函数的一组基，其他函数都可以表示为这组函数的线性组合
$\therefore$ 用这一组基信号可以表示所有其他信号

为什么是复函数？

$\begin{aligned} F(u) &= \frac{1}{M}\sum_{x=0}^{M-1}f(x)e^{-j2\pi ux/M} \\ &= \frac{1}{M}\sum_{x=0}^{M-1}f(x)(\cos2\pi ux/M -j\sin2\pi ux/M) \\ &= \frac{1}{M}\vec{f(x)}\cdot(\vec{\cos 2\pi ux/M}-j\vec{\sin 2\pi ux/M}) \end{aligned}$

$\cos 2\pi ux/M$ 和 $\sin 2\pi ux/M$ 的频率相同，但相位不同。引入复数变换的目的是为了刻画信号的相位。

$R(u)=\frac{1}{M}\vec{f(x)}\cdot \vec{\cos 2\pi ux/M},I(u)=-\frac{1}{M}\vec{f(x)}\cdot \vec{\sin 2\pi ux/M}$
$F(u)=R(u)+jI(u)$
幅度(magnitude)：$|F(u)|=[R^2(u)+I^2(u)]^{1/2}$
相角(phase angle)：$\phi(u)=\arctan[\frac{I(u)}{R(u)}]$

傅里叶逆变换

一维连续傅里叶变换

• 正向：$F(u)=\int_{x=-\infty}^{\infty}f(x)e^{-j2\pi ux}\,{\rm d}x$
• 逆向：$f(x)=\int_{x=-\infty}^{\infty}F(u)e^{j2\pi ux}\,{\rm d}u$

一维离散傅里叶变换

• 正向：$F(u)=\frac{1}{M}\sum_{x=0}^{M-1}f(x)e^{-j2\pi ux/M},u=0,1,2,\dots,M-1$
• 逆向：$f(x)=\sum_{x=0}^{M-1}F(u)e^{j2\pi ux/M},u=0,1,2,\dots,M-1$

为什么逆变换是这种形式？
$\because$ 傅里叶变换是矩阵（线性）变换，傅里叶基是正交基
$\therefore$ 逆变换即求矩阵的逆
$F(u)=\frac{1}{M}\sum_{x=0}^{M-1}f(x)(\cos2\pi ux/M -j\sin2\pi ux/M)$

$F= \begin{bmatrix} 1 & 1 & 1 & \cdots & 1 \\ 1 & \cos\delta x-j\sin \delta x & \cos 2\delta x-j\sin 2\delta x & \cdots & \cos(M-1)\delta x-j\sin(M-1)\delta x \\ 1 & \cos 2\delta x-j\sin 2\delta x & \cos 4\delta x-j\sin 4\delta x & \cdots & \cos 2(M-1)\delta x-j\sin 2(M-1)\delta x \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \cos(M-1)\delta x-j\sin(M-1)\delta x & \cos 2(M-1)\delta x-j\sin 2(M-1)\delta x & \cdots & \cos (M-1)^2\delta x-j\sin (M-1)^2\delta x \\ \end{bmatrix}$

$F^{-1}=\frac{1}{M}F^*$，$F^*$ 是 $F$ 的共轭转置

$F^{-1}=\frac{1}{M} \begin{bmatrix} 1 & 1 & 1 & \cdots & 1 \\ 1 & \cos\delta x+j\sin \delta x & \cos 2\delta x+j\sin 2\delta x & \cdots & \cos(M-1)\delta x+j\sin(M-1)\delta x \\ 1 & \cos 2\delta x+j\sin 2\delta x & \cos 4\delta x+j\sin 4\delta x & \cdots & \cos 2(M-1)\delta x+j\sin 2(M-1)\delta x \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \cos(M-1)\delta x+j\sin(M-1)\delta x & \cos 2(M-1)\delta x+j\sin 2(M-1)\delta x & \cdots & \cos (M-1)^2\delta x+j\sin (M-1)^2\delta x \\ \end{bmatrix}$

$\therefore$ 一维离散傅里叶变换可以写成

• 正向：$F(u)=\alpha\sum_{x=0}^{M-1}f(x)e^{-j2\pi ux/M}$
• 逆向：$f(x)=\beta\sum_{x=0}^{M-1}F(u)e^{j2\pi ux/M}$
• $\alpha\beta=\frac{1}{M}$

二维离散傅里叶变换

• 正向：$F(u,v)=\frac{1}{MN}\sum_{x=0}^{M-1}\sum_{y=0}^{N-1}f(x,y)e^{-j2\pi(ux/M+vy/N)}$
• 逆向：$f(x,y)=\sum_{u=0}^{M-1}\sum_{v=0}^{N-1}F(u,v)e^{j2\pi(ux/M+vy/N)}$
• 行列可分离

傅里叶频率图像的移中
给输入信号 $f(x,y)$ 乘上 $(-1)^{(x+y)}$，再进行傅里叶变换，可以将$F(0,0)$移到$F(M/2,N/2)$

$\begin{aligned} & \frac{1}{MN}\sum_{x=0}^{M-1}\sum_{y=0}^{N-1}f(x,y)(-1)^{x+y}e^{-j2\pi(ux/M+vy/N)} \\ = & \frac{1}{MN}\sum_{x=0}^{M-1}\sum_{y=0}^{N-1}f(x,y)e^{j\pi(x+y)}e^{-j2\pi(ux/M+vy/N)} \\ = & \frac{1}{MN}\sum_{x=0}^{M-1}\sum_{y=0}^{N-1}f(x,y)e^{j\pi(x+y)-j2\pi(ux/M+vy/N)} \\ = & \frac{1}{MN}\sum_{x=0}^{M-1}\sum_{y=0}^{N-1}f(x,y)e^{-j2\pi((u-M/2)x/M+(v-N/2)y/N)} \\ = & F(u-M/2,v-N/2) \end{aligned}$

实际上就是将图像分成四个象限，然后将对角位置的区域互换

傅里叶谱(The Fourier Spectrum)：$|F(u,v)|$

• $F(0,0)$ 是整个图像的平均灰度，也即最低频率
• $F(u,v)$ 在靠近$(u,v)=(0,0)$附近包含了低频信息：代表的是灰度变化缓慢的光滑区域
• 距离$(u,v)=(0,0)$ 越远，其频率越大
• $F(u,v)$在远离$(u,v)=(0,0)$区域包含了高频信息：代表的是灰度变化剧烈的区域，如边缘

通过移中操作将频率域的低频信息移到图像中央

快速傅里叶变换

（离散）傅里叶变换的复杂度，对长度为 $M$ 的信号，总共需要 $M^2$ 次求和运算和乘法运算。
快速傅里叶变换(FFT)只需要 $M\log_2 M$ 次运算。

FFT算法基本思想：基于逐次加倍(successive doubling)的方法，通过推导将原始傅里叶转换成两个递推公式。
$F(u)=\frac{1}{M}\sum_{x=0}^{M-1}f(x)e^{-j2\pi ux/M},u=0,1,2,\dots,M-1$

$\begin{cases} F(u)=\frac{1}{2}[F_{even}(u)+F_{odd}(u)W_{2K}^u] \\ F(u+K)=\frac{1}{2}[F_{even}(u)-F_{odd}(u)W_{2K}^u] \end{cases} u=0,1,2,\dots,M-1$

其中，$M=2K,W_M=e^{-j2\pi /M}$，$F_{even}(u),F_{odd}(u)$ 分别是偶数和奇数位置上 $K$ 个点的傅里叶值

频率域的图像滤波

卷积理论
傅里叶变换对的卷积理论：
空域上的卷积等价于频域上的乘积
频域上的卷积等价于空域上的乘积

空域和频域的滤波函数也构成傅里叶变换对，利用该性质可以在频域上设计滤波器，再转换到空域进行实现（通常效率更高）

振铃现象：对一幅图像进行滤波处理，若选用的频域滤波器具有陡峭的变化，则会使滤波图像产生“振铃”，所谓“振铃”，就是指输出图像的灰度剧烈变化处产生的震荡，就好像钟被敲击后产生的空气震荡。
出现的原因：接近窗函数的滤波器，逆傅里叶变换之后，其空域函数两边的余波将对图像产生振铃现象。