2024数字图像处理

😀 上大学宛如坐牢,而大三则是长达一年的禁闭

Chapter 1:Intro

图像表示

就是像素,在计算机程序中是用二维数组表示存储的:

$$f(x,y)=\left[\begin{array}{cccc}f(0,0)& f(0,1)& \cdots & f(0,N-1)\\ f(1,0)& f(1,1)& \cdots & f(1,N-1)\\ ⋮& ⋮& \ddots & ⋮\\ f(M-1,0)& f(M-1,1)& \cdots & f(M-1,N-1)\end{array}\right]$$

其中,$M$是图像的宽度，$N$是图像的高度。 图像的灰度级的计算方式是$L=2^K$,其中$K$是图像的位深度。

位深度其实就是图像的颜色深度,是一个像素颜色的位数,通常是8位或者16位.

因此存储一幅图像所需要的位数$b=M\ast N\ast K$

图像存储

交叉存储(Interlaced)

$$\left[\begin{array}{cccc}[B,G,R]& [B,G,R]& \cdots & [B,G,R]\\ [B,G,R]& [B,G,R]& \cdots & [B,G,R]\\ [B,G,R]& [B,G,R]& \cdots & [B,G,R]\\ ⋮& ⋮& \ddots & ⋮\\ [B,G,R]& [B,G,R]& \cdots & [B,G,R]\\ [B,G,R]& [B,G,R]& \cdots & [B,G,R]\\ [B,G,R]& [B,G,R]& \cdots & [B,G,R]\end{array}\right]$$

如上,也就是说是按照R→G→B→R→G→B的顺序进行存储的

连续存储(Sequential)

$$\left[\begin{array}{cccccc}[B]& [B]& [B]& [B]& \cdots & [B]\\ [B]& [B]& [B]& [B]& \cdots & [B]\\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ [G]& [G]& [G]& [G]& \cdots & [G]\\ [G]& [G]& [G]& [G]& \cdots & [G]\\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ [R]& [R]& [R]& [R]& \cdots & [R]\\ [R]& [R]& [R]& [R]& \cdots & [R]\\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \end{array}\right]$$

如上，也就是说是按照B→…→B→G→…→G→R→…→R的顺序进行存储的

例题1

一幅8位单通道图像中，某像素的内存地址为p，图像的步长为step（每行所占的字节数为step），其八邻域像素的内存地址是?

例题2

一副32位单通道图像，每一行的字节数为step，共有weight行，high列。（0,0）的存储位置是p。 1.求(u,v)的存储位置: $location=p+4\ast (v\ast step+u)$(乘4是因为32位四字节) 2.求(u,v)的八邻域存储位置:见上 3.(x,y)为浮点数坐标，求其双线性拟合值，用p表示: u = floor(x) v = floor(y),dx = x - u dy = y - v

$$\begin{gathered} p=(1-dx)\ast (1-dy)\ast img[v\ast step+u]+dx\ast (1-dy)\ast img[v\ast step+u+1]+ \\ (1-dx)\ast dy\ast img[(v+1)\ast step+u]+dx\ast dy\ast img[(v+1)\ast step+u+1] \end{gathered}$$

图像数字化

采样 图像在空间上的离散化称为采样，也就是用空间上部分点的灰度值代表图像，这些点称为采样点. 离散化就涉及到一个脉冲函数与原函数求交,最终离散结果的还原程度很大程度上取决于采样的脉冲函数频率(在图像上反映出来是密度)有多高. 课堂PPT的原句为:“Spatial sampling determines the resolution of an image(what level of detail can be seen)“采样越细即采样点越多,像素越小,越能精细地表现图像

量化 把采样后所得的各像素的灰度变换成离散的整数值称为图像灰度的量化. 同采样一样,也是一个离散化的过程,因此对于精细度也有类似的影响. 量化越细致即设置灰度值范围越大，图像灰度表现越丰富

Chapter 2:Fundamentals

亮度相对性

第一个出现的重要知识点就挺抽象的

背景亮度为I，前景和背景亮度差距的绝对值为ΔI，所以能用 $\Delta I/I$ 来表示前景和背景的亮度差距。对于这个例子来说，因为B的前景相对亮背景而言差距更小，因此B看起来更亮。

be like:

灰度插值

也就是对灰度值图像的像素点进行差值,一般的使用双线性插值就足够了

步骤：

1. 对上端的两个顶点进行线性插值，可得：$f(x,0)=f(0,0)+x[f(1,0)-f(0,0)]$
2. 对下端的两个顶点进行线性插值，可得：$f(x,1)=f(0,1)+x[f(1,1)-f(0,1)]$
3. 对上述两个结果进行线性插值，可得：$f(x,y)=f(x,0)+y[f(x,1)-f(x,0)]$

简单来讲就是横竖各插值一次,插值结果再插值一次

Resampling

对于放大缩小的图像(尤其是放大后的)产生的新像素的值需要由原始图像获取,但是如果新的像素反向映射到原图像中不是整数位置则需要差值.主要有三种方式:

1. 最近邻插值:其实没插,找距离最近的像素直接使用它的值
2. 双线性插值:见上文
3. 双三次插值:比双线性插值多一层,考虑周围16个像素的值,四个一组双线性插值一次,插值出的四个值作为四个”像素值”再双线性差值一次

由此可以看出计算量是几何倍上升的,质量也会越来越好

区域属性

区域属性是区域被标记(分割)后进行区域分析(或测量)的重要属性.

例如周长和面积.

给定一个简单多边形，其边界上的点坐标都是整数坐标（即网格点），用 Pick’s Formula 可以计算该多边形的面积$A(p)$:

$$A(p)=n_I+\frac{n_B}{2}+1$$

其中：

• $A(p)$ 是多边形的面积；
• $n_I$ 是多边形内部的网格点数（即多边形内部的整数点数）；
• $n_B$ 是多边形边界上的网格点数（即多边形边界上的整数点数）。

Chapter 3:Enhancement

翻译过来是空域上的图像增强

灰度变换

其中，对数变换的公式为：

$$s=clog(1+r)$$

幂律变换（Power-Law Transformations）的公式为：

$$s=c{r}^{\gamma }$$

其中,c和γ是正常数(正的常数),r是原始图像的灰度值,s是变换后的灰度值

明显对数变换和幂律变换的操作是相反的,但是最终的效果都在某种程度上做到了增大图像对比度

• 对数变换明显的,拓展了低的灰度值范围,压缩高灰度值;因此可以增强暗部细节,同时可以压缩图像的动态范围,尤其是高灰度值范围(在傅里叶变换展示频谱图时很有用)

幂律变换又称伽马校正，作用是避免灰度失真，校正显示屏的非线性变换。显然，随着γ的增大，曲线向下凹，图像的对比度增大.

• 幂律变换与对数相反,拓展高灰度值,压缩低灰度值;这种变换首先可以进行伽马校正(用于补偿显示设备的非线性特性，确保图像在不同设备上显示的一致性)相对于对数,可以增强亮部细节

灰度变换的输入输出当然可以是任意函数，因此可以随意调整，比如高亮某个灰度范围等等.

位平面切片

位平面切片（Bit-Plane Slicing）是一种灰度变换，它将图像的每个像素的二进制表示中的每一位作为一个灰度级，从而将图像分解为8个位平面($255=2^8-1$).

下图中展示了将lena图分解后，各个平面图的样子，其中（a）为lena的原图，从(b)到（i)展示了从最低位平面到最高位平面，其中我们可以很明显发现第八张子图(i)最接近lena的原图.

能看出:越是高位，其位平面越是低频，越保留图像的总体样子；越是低位，其位平面越是高频，越保留图像的细节.

这样看来,某种程度上,位平面其实也能算是一种基函数分割,某种傅里叶变换((

Q：如何从一幅图像中抽出某一张位平面切片?

A：拿二进制串（如00001000）去和灰度值二进制按位求与（AND）.

二值化

1.Isodata算法

Isodata算法是一种自动确定阈值的方法，它是一种迭代算法，其基本思想是：首先假设阈值为$T$，然后将图像中的所有像素分为两类：灰度值大于$T$的像素为一类，灰度值小于$T$的像素为另一类。然后计算这两类像素的平均灰度值$m_1$和$m_2$并计算新的阈值$T=(m_1+m_2)/2$。重复上述过程，直到$T$不再变化为止。

2.Otsu算法

Otsu算法基于聚类的思想，将图像按照灰度值分为两类，使得类内方差最小，类间方差最大。通过计算方差来寻找合适的阈值。

直方图

真好啊真好,数媒就是个小概统

直方图（Histogram）是一种统计报告，它对图像中每个灰度级的像素数进行计数.

1.直方图均衡化

直方图均衡化（Histogram Equalization）是一种用于增强图像对比度的图像处理技术.它通过重新分布图像的像素值，使得图像的直方图变得更加均匀，从而增强图像的视觉效果.

直方图均衡化的基本思想是将图像的像素值映射到一个新的值，使得映射后的直方图具有均匀分布.这样做的目的是使得图像中的像素值范围更广，从而增加图像的对比度

下面是直方图均衡化的一般步骤：

1. 计算原始图像的灰度直方图：将图像转换为灰度图像，并统计每个灰度级别的像素数目，得到原始图像的直方图。
2. 计算累积分布函数（CDF）：对原始图像的直方图进行归一化，得到每个灰度级别的累积概率。CDF表示了每个灰度级别在图像中的累积分布情况。
3. 映射像素值：根据CDF，计算每个灰度级别的新的映射值。这可以通过将CDF值乘以最大灰度级别来实现。
4. 应用映射：将原始图像中的每个像素值替换为其对应的映射值。这样就完成了直方图均衡化的过程。

直方图均衡化可以使得图像的像素值范围更广，从而增加图像的对比度.它在许多图像处理应用中都有广泛的应用，例如图像增强、图像匹配和图像分割等。然而，直方图均衡化也可能导致图像的噪声增加，因为它会放大图像中的细微差异.

2.直方图规定化

直方图均衡不允许交互式图像增强，并且仅生成一个结果：均匀直方图的近似值。但有时，我们需要能够指定能够突出显示某些灰度范围的特定直方图形状。为此，我们需要直方图规定化（Histogram Matching）。

步骤:

1. 均衡化原图:$s=T(r_k)={\sum }_{j=0}^k\frac{n_j}{n}$其中,n是像素总数,$n_j$是灰度级别j的像素数(其实就是计算$n_k$的CDF的过程)
2. 指定所需的密度函数并获得变换函数$G(z)$.$v=G(z)={\sum }_0^k{p}_z(w)\approx {\sum }_{j=0}^z\frac{n_j}{n}$其中,$p_z(w)$是所需的密度函数,$n_j$同上
3. 对步骤一中的结果进行反变换:$z=G^{-1}(s)=G^{-1}[T(r_k)]$,其中,$G^{-1}$是$G$的逆变换

反变换函数主要有两种：

1. SML（Single Mapping Law）
   1. 有偏
   2. 自原始累积直方图到匹配累积直方图寻找对应关系
   3. 简单，但有时有较大的取整误差
2. GML（Group Mapping Law）
   1. 无偏
   2. 自匹配累积直方图到原始累积直方图寻找对应关系
   3. GML产生的误差远小于SML

局部增强

对局部增强，我们需要根据每个像素邻域的灰度分布设计变换函数，基于像素邻域对灰度的贡献设计变换函数。

说人话就是前面的方法不要全局应用了，划个小区域做。定义一个正方形（或矩形）邻域，并将该区域的中心从一个像素移动到另一个像素。

两种方式：Pixel-to-pixel的（类似于滤波核）和 Non-overlapping region的（就是说划分成多个不重叠的矩形区域）

空域滤波

直接上才艺!!

均值滤波:

$$g(x,y)=\frac{1}{9}\sum _{s=-1}^1\sum _{t=-1}^{1}f(x+s,y+t)$$

最大值滤波:

$$g(x,y)=\underset{s,t}{\max }f(x+s,y+t)$$

最小值滤波:

$$g(x,y)=\underset{s,t}{\min }f(x+s,y+t)$$

线性滤波(也就是加权均值滤波):

$$g(x,y)=\sum _{s=-1}^1\sum _{t=-1}^{1}w(s,t)f(x+s,y+t)$$

其中, $w(s,t)$ 是权重函数, 满足 ${\sum }_{s=-1}^1{\sum }_{t=-1}^{1}w(s,t)=1$

有一些滤波器是行列可分离的，在二维分离意味着，卷积核可以分解成两个一维向量的乘积

$$\begin{array}{rl}g(x,y)& =\sum _{s=-a}^a\sum _{t=-b}^{b}w(s,t)f(x+s,y+t)\\ & =\sum _{s=-a}^a{w}_1(x)\sum _{t=-b}^b{w}_2(y)f(x+s,y+t)\end{array}$$

二维高斯函数这种可分离为两个一维向量(纵横)相乘的滤波器就可以

二维高斯函数卷积可分两步进行：首先将图像与一维高斯函数卷积，然后将卷积结果与方向垂直的相同一维高斯函数卷积。因此，二维高斯函数的计算量随模板宽度成线性增长而不是成平方增长，所以较大尺寸的高斯滤波器可以得以有效的实现

(假设图像大小为mn,滤波器大小为ab,直接二维计算量为mnab,而分离行列后则是(m+n)ab)

低通滤波

均值滤波就可以，不过不如Gaussian滤波好。

Gaussian滤波：

$$g(x,y)=\frac{1}{16}\left[\begin{array}{lll}1& 2& 1\\ 2& 4& 2\\ 1& 2& 1\end{array}\right]\ast f(x,y)$$

高通滤波

滤波核中心为正，边缘为负，这样就可以突出差异，可以说非常高频了。

一个简单的高通滤波器：

$$g(x,y)=\left[\begin{array}{ccc}-1& -1& -1\\ -1& 8& -1\\ -1& -1& -1\end{array}\right]\ast f(x,y)$$

来点 Laplacian 滤波器：

$$g(x,y)=\left[\begin{array}{ccc}0& -1& 0\\ -1& 4& -1\\ 0& -1& 0\end{array}\right]\ast f(x,y)$$

这是一个4-邻域的 Laplacian 滤波器，如果加上原图像，就可以锐化边缘，即：

$$\left[\begin{array}{ccc}0& -1& 0\\ -1& 4& -1\\ 0& -1& 0\end{array}\right]+\left[\begin{array}{lll}0& 0& 0\\ 0& 1& 0\\ 0& 0& 0\end{array}\right]=\left[\begin{array}{ccc}0& -1& 0\\ -1& 5& -1\\ 0& -1& 0\end{array}\right]$$

8-邻域的 Laplacian 滤波器：

$$g(x,y)=\left[\begin{array}{ccc}-1& -1& -1\\ -1& 8& -1\\ -1& -1& -1\end{array}\right]\ast f(x,y)$$

注意滤波器中的系数和为 0

微分滤波器

1. 一阶导数在图像中产生较厚的边缘：
   1. 一阶导数滤波器对图像中的边缘具有较强的响应。当应用一阶导数滤波器时，边缘区域的像素值变化较大，导致滤波器的响应较强，从而产生较厚的边缘。
2. 二阶导数对精细细节有更强的响应：
   1. 二阶导数滤波器对图像中的细节变化具有更敏感的响应。细节通常包含更高频率的变化，而二阶导数滤波器对高频信号有更强的响应，因此它能够捕捉到图像中的精细细节。
3. 一阶导数对灰度阶跃有更强的响应：
   1. 当图像中出现灰度阶跃（像素值突然变化）时，一阶导数滤波器能够产生较强的响应。这是因为阶跃边缘的像素值变化较大，导致一阶导数滤波器的输出较大。
4. 二阶导数在灰度阶跃变化时产生双重响应：
   1. 与一阶导数不同，二阶导数滤波器在灰度阶跃变化时会产生双重响应。这是因为二阶导数滤波器对图像中的边缘变化有两个极值点（正负），因此在阶跃边缘的两侧都会产生响应。
5. 二阶导数对线段的响应强于阶跃的响应，对点的响应强于线段的响应：
   1. 二阶导数滤波器对不同类型的特征有不同的响应强度。它对线段的响应较强，因为线段具有较大的二阶导数值。相比之下，阶跃边缘的二阶导数值较小。对于点（如孤立的像素），二阶导数滤波器的响应最强，因为点处的二阶导数为零，而在点周围的像素处有较大的二阶导数值。

我们约定一个3*3领域如下:

$$\left[\begin{array}{ccc}{z}_1& z_2& z_3\\ z_4& z_5& z_6\\ z_7& z_8& z_9\end{array}\right]$$

梯度算子:

$$\Delta f=\left[\begin{array}{c}G(x)\\ G(y)\end{array}\right]=\left[\begin{array}{c}\frac{\partial f}{\partial x}\\ \frac{\partial f}{\partial y}\end{array}\right]$$

平均类似于积分，会导致模糊，因此差分预计会产生相反的结果并使图像锐化

一元微分

$$\begin{gathered} f^{\prime }=f(x+1)-f(x) \\ {f}^{\prime \prime }=f(x+1)+f(x-1)-2f(x) \end{gathered}$$

Roberts算子(以上文3*3领域为例)

$$\begin{gathered} G_x=(z_9-z_5) \\ {G}_y=(z_8-z_6) \end{gathered}$$

也就对角线做差,右下减左上,左下减右上

近似（Roberts 交叉梯度算子）:

$$\Delta f\approx |z_9-z_5|+|z_8-z_6|$$

Roberts 算子模板:

• $G_x$:

  $$\left[\begin{array}{cc}-1& 0\\ 0& 1\end{array}\right]$$

• $G_y$:

  $$\left[\begin{array}{cc}0& -1\\ 1& 0\end{array}\right]$$

Prewitt算子

• 近似(Prewitt交叉梯度算子):

$$\nabla f=|(z_7+z_8+z_9)-(z_1+z_2+z_3)|+|(z_3+z_6+z_9)-(z_1+z_4+z_7)|$$

• $G_x$:

  $$\left[\begin{array}{ccc}-1& -1& -1\\ 0& 0& 0\\ 1& 1& 1\end{array}\right]$$

• $G_y$:

  $$\left[\begin{array}{ccc}-1& 0& 1\\ -1& 0& 1\\ -1& 0& 1\end{array}\right]$$

sobel算子

• 近似(sobel交叉梯度算子)

  $$\nabla f=|(z_7+2z_8+z_9)-(z_1+2z_2+z_3)|+|(z_3+2z_6+z_9)-(z_1+2z_4+z_7)|$$

• $G_x$:

  $$\left[\begin{array}{ccc}-1& -2& -1\\ 0& 0& 0\\ 1& 2& 1\end{array}\right]$$

• $G_y$:

  $$\left[\begin{array}{ccc}-1& 0& 1\\ -2& 0& 2\\ -1& 0& 1\end{array}\right]$$

二阶微分

$${\nabla }^{2}f=[f(x+1,y)+f(x-1,y)+f(x,y+1)+f(x,y-1)]-4f(x,y)$$

二阶微分只说 Laplacian:

$$g(x,y)=\{\begin{array}{ll}f(x,y)-{\nabla }^{2}f(x,y)& \text{如果滤波器中心是负数}\\ f(x,y)+{\nabla }^{2}f(x,y)& \text{如果滤波器中心是正数}\end{array}$$

例题:根据滤波器的可分离性，说明为什么prewitt、sobel算子抗噪性比简单梯度滤波器好

首先是二者的自身特性,不难看出两个都哟加权平均的成分,这一点与高斯滤波器很相似,有效抑制高频噪声 再者,根据可分离性,分离后的滤波器相当于窗口内所有像素都有参与,而并非只有梯度方向

总结

平滑

• 正值
• 和为 1
• 平滑量与模板尺寸成正比
• 过低频，滤高频

锐化

• 有正有负
• 和为0
• 锐化效果和滤波核中心权重有关
• 滤低频,过高频

微分

• 相反的符号（就是说有正有负）用于在高对比度区域获得高响应
• 和为 0，这是为了不响应平坦区域
• 高对比度点的绝对值高

Chapter 4:Frequency

终于来到这部分了.

空域表示：一幅图像是一个二维函数$f(x,y)$,其中xx和yy是空间坐标，幅值$f$在任何一点$(x,y)$处的值称为该点的强度或灰度.

类似于 Taylor 级数,任何连续周期信号都可以表示为一系列正弦曲线的和.频率从低到高排列,一些成分为直线是因为其振幅为零.也就是说,为了组成特殊的信号,有些正弦波成分是不需要的.

不同频率的正弦波被称为频率分量因为正弦波是一个圆周运动在一条直线的投影,所以频域的基本单元也可以理解成一个始终在旋转的圆.

在频域中，频率越大说明原始信号变化速度越快；频率越小说明原始信号越平缓。当频率为 0 时，表示直流信号，没有变化。因此，频率的大小反应了信号的变化快慢。且高频分量解释信号的突变部分，而低频分量决定信号的整体形象.

• 空域：图像的原始表示，即图像的像素矩阵.
• 频域：通过 Fourier 变换将图像从空域转换到频域。处理的是图像的频率分量，而不是图像的像素值.
• 时域：图像灰度的分布函数。所谓空域通过 Fourier 变换转频域，倒不如说是时域通过 Fourier 变换转频域.

在开始和奇妙又倒霉催的Fourier打交道之前,我们先约定符号,我们希望在整章中,公式中的符号总是一致的.

• $f(x,y)$：一般表示图像的空域表示，也就是图像的像素矩阵。其实，倒不如说是图像的时域表示。
• $F(u,v)$：图像的频域表示，也就是图像的频率分量。
• F：Fourier 变换。
• 和幻灯片中不一样，我们直接使用i表示虚数单位，说实话使用j也一样，没什么区别.

Fourier变换

$$F(u)={\int }_{-\infty }^{\infty }f(x)e^{-i2\pi ux}dx$$

其中,$F(u)$是频率为$u$的正弦波在$\ast f(x)\ast$中的振幅.

我们也很关心滤波函数$H(u)$，它处理 Fourier 变换后的图像（也就是频域），然后再进行反变换回到空域

一维连续傅里叶变换

傅里叶:

$$F(u)={\int }_{-\infty }^{\infty }f(x)e^{-i2\pi ux}dx$$

逆傅里叶变换:

$$f(x)={\int }_{-\infty }^{\infty }F(x)e^{i2\pi ux}du$$

可以这样理解：Fourier 逆变换中，$e$的指数是 Fourier 变换中的相反数，这样，$e$的指数相加就是 0，如果将一个 Fourier 变换和一个 Fourier 逆变换连续进行，那么积分中的$e$的指数相加就是 0，积分结果就是1，也就是说，连续进行 Fourier 变换和 Fourier 逆变换，最后结果不变。

Fourier 变换的物理意义：将图像的灰度分布函数变换为图像的频率分布函数。

Fourier 逆变换：将图像的频率分布函数变换为图像的灰度分布函数。

基函数

基函数：

$${\varphi }_u(x)=e^{-i2\pi ux}$$

欧拉公式:

$$e^{i\theta }=\cos \theta +i\sin \theta$$

因此傅里叶变换套欧拉公式为:

$$\begin{array}{rl}F(u)& ={\int }_{-\infty }^{\infty }f(x)e^{-i2\pi ux}\mathrm{d}x\\ & ={\int }_{-\infty }^{\infty }f(x)\left[\cos (2\pi ux)-i\sin (2\pi ux)\right]\mathrm{d}x\end{array}$$

其中，$u$体现了频率性质.

一维离散傅里叶

将离散 Fourier 变换应用采样函数，有DFT：

$$F(u)=\frac{1}{M}\sum _{x=0}^{M-1}f(x)e^{-i2\pi ux/M}$$

其中，MM 是采样点数，$u=0,1,2,\cdots ,M-1$。请注意，$u$的取值范围是$[0,M-1]$，每个都要取,下面引述 DFT 公式的时候不再赘述.

同理,再看一维离散Fourier 逆变换:

$$f(x)=\sum _{u=0}^{M-1}F(u)e^{i2\pi ux/M}$$

采样

采样增量在空域和频域的关系:

$$\begin{gathered} \Delta x=\frac{1}{M\Delta u} \\ \Delta y=\frac{1}{N\Delta v} \end{gathered}$$

其中，MM 和 NN 是采样点数，Δx和Δy是空域采样增量，Δu和Δv是频域采样增量。

这样，DFT 可以写成:

$$\begin{array}{rl}F(u)& =\frac{1}{M}\sum _{x=0}^{M-1}f(x)e^{-i2\pi ux/M}\\ & =\frac{1}{M}\sum _{x=0}^{K-1}f(x)A{e}^{-i2\pi ux/M}\\ & =\frac{A}{M}\sum _{x=0}^{K-1}f(x)e^{-i2\pi ux/M}\end{array}$$

其中：

• A：表示一个常数或者系数，用来调整 DFT 的幅度。在公式中，A是一个常数，可以用于缩放 DFT 的结果。通常情况下，A的取值为 1，表示不进行幅度调整。
• M：表示离散傅里叶变换（DFT）的长度或者采样点数。它决定了 DFT 的频率分辨率和计算的复杂度。在公式中，M代表了 DFT 的长度，也就是采样点的个数。M通常是一个正整数，且必须大于等于K。
• K：表示输入序列的长度或者有效数据点数。在公式中，K代表了输入序列的长度，也就是有效数据点的个数。K通常是一个正整数，且必须小于等于M。

需要注意的是，M和K的取值关系是K≤M，即有效数据点数K不能超过DFT的长度M。当K<M时，DFT公式中的求和项只对前K个数据点进行求和，忽略了后面的M−K个数据点。这样可以在计算 DFT 时只考虑有效数据，提高计算效率。

为什么要区分K和M呢？这是因为在 DFT 中，我们通常将输入序列看作是周期为M的周期信号，其中只有前K个数据点是有效的，后面的M−K个数据点被认为是零值或者被忽略。

而F(0)自然就有：

对于图像来说,F(0)不就是整个图像的平均亮度吗(有效点占比)

二维傅里叶

连续二维傅里叶:

$$F(u,v)={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }f(x,y)e^{-i2\pi (ux+vy)}\mathrm{d}x\mathrm{d}y$$

连续二维傅里叶逆:

$$f(x,y)={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }F(u,v)e^{i2\pi (ux+vy)}\mathrm{d}u\mathrm{d}v$$

其中，u和v是频率空间的变量.

离散二维傅里叶:

$$F(u,v)=\frac{1}{MN}\sum _{x=0}^{M-1}\sum _{y=0}^{N-1}f(x,y)e^{-i2\pi (ux/M+vy/N)}$$

其中,M和N是采样点数。正如一维离散Fourier变换一样，N其实就是第二个维度上的M。

二维离散 Fourier 逆变换为

$$f(x,y)=\sum _{u=0}^{M-1}\sum _{v=0}^{N-1}F(u,v)e^{i2\pi (ux/M+vy/N)}$$

Fourier频谱

$$|F(u,v)|=\sqrt{{\mathrm{Re}}^{2}F(u,v)+{\mathrm{Im}}^{2}F(u,v)}$$

其中,Re F(u,v)和Im F(u,v)分别是F(u,v)的实部与虚部

---

相位谱（Phase Spectrum）：

$$\varphi (u,v)=\arctan \left(\frac{\mathrm{Im}F(u,v)}{\mathrm{Re}F(u,v)}\right)$$

---

功率谱（Power Spectrum）：

$$\begin{array}{rl}P(u,v)& =|F(u,v){|}^2\\ & ={\mathrm{Re}}^{2}F(u,v)+{\mathrm{Im}}^{2}F(u,v)\end{array}$$

傅里叶频谱移中

在实际图像频谱分析中，由于低频分量区域较小，并且分散在四角，因此不利于对它进行分析。根据图像频谱的周期性和共轭对称性对频谱坐标移位，将所有的低频分量集中到频谱图中心，同时高频分量分散在四周。

所谓移中处理，就是将$F(0,0)$移动到$F(M/2,N/2)$处。

在证明这个事情之前，我们首先要证明 Fourier 变换的线性性质，即

$\mathfrak{F}[af(x,y)+bg(x,y)]=a\mathfrak{F}[f(x,y)]+b\mathfrak{F}[g(x,y)]$

证明如下:

$$\begin{array}{rl}\mathfrak{F}[af(x,y)+bg(x,y)]& ={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }[af(x,y)+bg(x,y)]e^{-i2\pi (ux+vy)}\mathrm{d}x\mathrm{d}y\\ & =a{\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }f(x,y)e^{-i2\pi (ux+vy)}\mathrm{d}x\mathrm{d}y+b{\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }g(x,y)e^{-i2\pi (ux+vy)}\mathrm{d}x\mathrm{d}y\\ & =a\mathfrak{F}[f(x,y)]+b\mathfrak{F}[g(x,y)]\end{array}$$

说白了就是积分的线性性质。积分是可加的.

卷积定理又说:在时域和频域之间，卷积操作和乘积操作是互相对偶的。也就是说，时域的卷积在频域中是乘积，时域的乘积在频域中是卷积.则有:

$$\begin{array}{rl}\mathfrak{F}\left[f(x,y)\cdot g(x,y)\right]& =\frac{1}{2\pi }\mathfrak{F}[f(x,y)]\ast \mathfrak{F}[g(x,y)]\\ & =\frac{1}{2\pi }F(u,v)\ast G(u,v)\end{array}$$

这个是频域卷积定理,等式左侧是数乘操作,右侧是卷积. (如果是时域卷积定理的话,等式左侧是卷积,右侧是数乘,且没有1/2π系数)

因此可得:

$$\begin{array}{rl}\mathfrak{F}\left[f(x,y)(-1{)}^{x+y}\right]& =\sum _{x=0}^{M-1}\sum _{y=0}^{N-1}f(x,y)(-1{)}^{x+y}{e}^{-i2\pi (ux/M+vy/N)}\\ & =\frac{1}{MN}\sum _{x=0}^{M-1}\sum _{y=0}^{N-1}f(x,y)e^{-i2\pi ((u-M/2)x/M+(v-N/2)y/N)}\\ & =F(u-M/2,v-N/2)\end{array}$$

$\mathfrak{F}$表示Fourier变换，$f(x,y)(-1{)}^{x+y}$表示将图像的原点移动到图像的中心。之所以要乘以$(-1{)}^{x+y}$，是因为这样可以使得图像的中心在原点，而不是在图像的四角。上面的式子证明了这个移中处理的正确性，它使得低频分量集中在频谱图的中心，而高频分量分散在四周.

解释

To be honest这个应该反向推理,使用了傅里叶变换的频移属性:

$$\begin{gathered} f(t)=\frac{1}{2\pi }{\int }_{-\infty }^{\infty }F(j(\omega -{\omega }_0))e^{j\omega t}d\omega \stackrel{\omega ={\omega }^{\prime }+{\omega }_0}{\to }\frac{1}{2\pi }{\int }_{-\infty }^{\infty }F(j{\omega }^{\prime })e^{j({\omega }^{\prime }+{\omega }_0)t}d{\omega }^{\prime } \\ =e^{j{\omega }_{0}t}\cdot \frac{1}{2\pi }{\int }_{-\infty }^{\infty }F(j\omega )e^{j\omega t}d\omega =e^{j{\omega }_{0}t}f(t) \end{gathered}$$

傅里叶频谱应用

Fourier 谱本身的动态范围很大，只有最亮的部分频谱可见，这会导致直接输出的 Fourier 谱图难以观察。因此，我们需要对 Fourier 谱图进行对数变换，以便于观察。

我们一般会走一个对数变换，将其变换到一个合适的范围内。

$D(u,v)=clog(1+\mid F(u,v)\mid )$

这个在之前的对数变换中有提到过.其中，c是一个常数，用来调整对数变换的幅度；之所以要加1，是为了避免对数变换中出现 0，导致无法计算.

傅里叶频谱图结论

回顾 Fourier 频谱图:

$$|F(u,v)|=\sqrt{{\mathrm{Re}}^{2}F(u,v)+{\mathrm{Im}}^{2}F(u,v)}$$

其中，Re⁡F(u,v)和 Im⁡F(u,v)分别 F(u,v)的实部和虚部。

移中后，图像中心表示低频分量，图像四角表示高频分量。所以通常来说，Fourier频谱图都是中心亮，四角暗的。

对于 Fourier 频谱图来说，显然，频谱图中心即为 F(0,0)F(0,0)，而图像四角则是 F(M/2,N/2)。F(0,0)表示图像的直流分量(也就是 Fourier 变换后的最低频率的那条直线)，也就是图像的平均灰度。

傅里叶描述器

这个东西有点抽象且不知道怎么用,大概可以理解为对复域的傅里叶变换与逆变换.

如果将一个坐标(x,y)看作是一个复数，然后对一系列点施加Fourier变换，那么就可以得到一个复数序列，这个序列就是Fourier描述器。

首先，我们将一个坐标(x,y)看作是一个复数：

$$z=x+iy$$

然后，对一系列点施加离散Fourier变换

$$F(u)=\frac{1}{K}\sum _{k=0}^{K-1}(x+iy)e^{-i2\pi uk/K}$$

其中K是点个数,重建原始序列:

$$\begin{array}{rl}& \hat{z}=\frac{1}{K}\sum _{u=0}^{K-1}F(u)e^{i2\pi uk/K}\\ & =\frac{1}{K}\sum _{u=0}^{K-1}\left[\frac{1}{K}\sum _{k=0}^{K-1}(x+iy)e^{-i2\pi uk/K}\right]e^{i2\pi uk/K}\\ & =\frac{1}{K}\sum _{u=0}^{K-1}\sum _{k=0}^{K-1}(x+iy)e^{-i2\pi uk/K}{e}^{i2\pi uk/K}\\ & =\frac{1}{K}\sum _{u=0}^{K-1}\sum _{k=0}^{K-1}(x+iy)\\ & =\frac{1}{K}K(x+iy)\\ & =z\end{array}$$

其实就是逆傅里叶变化吧,对一个复数域的点序列.

取不同的K，可以得到不同的 Fourier 描述器。K越大，描述器越精确，但相应地，对噪声也越敏感，计算量也越大；K越小，描述器越不精确，但相应地，鲁棒性更强，计算量也越小

二维傅里叶变换性质

其实很多傅里叶变换的性质也就是其根本:积分的性质

1.线性性质

线性性质在上面提到过,系数可以提到傅里叶变换外面

2.分配率

积分的分配率,令线性性质中两个系数都为1

$$\mathfrak{F}[f(x,y)+g(x,y)]=\mathfrak{F}[f(x,y)]+\mathfrak{F}[g(x,y)]$$

注意,只对加法成立,乘法就涉及到卷积定理的内容了,自然是不满足分配率的

3.平移

在上面的移中操作中我们提到了一些平移性质,那部分是频移,也就是频域函数平移

$$\begin{array}{rl}f(x,y)e^{i2\pi (u_{0}x/M+v_{0}y/N)}& ⟺F(u-u_0,v-v_0)\\ f(x-x_0,y-y_0)& ⟺F(u,v)e^{-i2\pi (u{x}_0/M+v{y}_0/N)}\end{array}$$

两个表达式都比较明显,重点记住平移出的系数.

注意，平移性质当然不会影响频率振幅

4.缩放

缩放性质描述如下：

$$\begin{gathered} af(x,y)⟺aF(u,v) \\ f(ax,by)⟺\frac{1}{|ab|}F\left(\frac{u}{a},\frac{v}{b}\right) \end{gathered}$$

第一个式子表明：对空域的幅度的缩放，将导致频域的幅度的缩放

第二个式子表明:对空域的空间尺度缩放(就是对xy缩放),会导致频域的空间尺度反方向缩放,还多了个系数表示频域空间幅度也会进行缩放.也就是说,如果$a<1,b<1$,那么频域的空间尺度和幅度都会增大

其实，还能注意到，空域中波形传播的方向和频域中亮点排列的方向是一样的。

5.旋转

在旋转中用极坐标来表示:

$$\begin{gathered} x=r\cos \theta \\ y=r\sin \theta \\ u=\omega \cos \varphi \\ v=\omega \sin \varphi \end{gathered}$$

这样就把直角坐标系转到极坐标系下,旋转操作就可以变为线性加和

带入傅里叶变换:

$$f(r,\theta +{\theta }_0)⟺F(\omega ,\varphi +{\theta }_0)$$

即，对$f(r,\theta )$旋转${\theta }_0$，将导致$F(\omega ,\varphi )$同样旋转${\theta }_0$。

与上文中”其实，还能注意到，空域中波形传播的方向和频域中亮点排列的方向是一样的”的方向性描述是一致的,或者说，正是因为方向总是一致，所以空域旋转多少角度，频域就要旋转多少角度.

6.周期性与共轭对称性

对实函数$f(x,y)$，其傅里叶变换结果$F(u,v)$具有共轭对称性

$$F(u,v)=F^{\ast }(-u,-v)$$

其中,F∗表示F的共轭复数（我们简称为复共轭）.当然了，写成这样也是可以的：

$$|F(u,v)|=|F(-u,-v)|$$

本质上，这个性质来自于Fourier变换的周期性.三角函数可不就有周期性吗？

然后，我们以一维Fourier变换为例，来说明这个性质。F(u)本身就有周期N，所以F(u)=F(u+N)，或者说 ∣F(u)∣=∣F(−u)∣.第二个式子成立是因为F(u)的对称中心是原点.

7.可分离性与卷积定理

DFT 可以表示为可分离的形式:

$$\begin{gathered} F(u,v)=\frac{1}{M}\sum _{x=0}^{M-1}F(x,v)e^{-i2\pi ux/M} \\ F(x,v)=\frac{1}{N}\sum _{y=0}^{N-1}f(x,y)e^{-i2\pi vy/N} \end{gathered}$$

换句话说，一个二维 DFT 可以分解为：先对每一行进行一维 DFT，然后对每一列进行一维 DFT。

另外，公式表明，对于任意值x，F(x,v)都是一维 DFT 的结果。因此，我们可以将二维 DFT 看作是一维 DFT 的一系列结果

其实和二维滤波核分解为两个一维的相当类似不是吗?

卷积定理就不用多说了,再上文多次提到.

FFT(快速傅里叶变换)

包快的

fft是一种快速计算DFT的散发，计算复杂度为$O（NlogN）$而朴素的DFT计算复杂度为$O（N^2）$。

朴素的DFT就不说了上文都有，细说在求和内的公式用欧拉公式展开后其实是f(x)和一个复数相乘。

那么FFT是怎么个F法呢，课上讲的FFT是最常见的“Cooley-Tukey“算法，但是FFT的实现方法不只有这一种

FFT算法是基于逐次加倍方法.通过推导将原始傅里叶变换转化成两个递推公式：

$$\begin{gathered} F(u)=\frac{1}{M}\sum _{x=0}^{M-1}f(x)e^{-i2\pi ux/M} \\ F(u)=\frac{1}{2}\left[F_{\mathrm{even}}(u)+F_{\mathrm{odd}}(u)W_{2k}^u\right] \\ F(u+k)=\frac{1}{2}\left[F_{\mathrm{even}}(u)-F_{\mathrm{odd}}(u)W_{2k}^u\right] \end{gathered}$$

看起来嘎嘎复杂,其实一点也不嘎嘎,但真的复杂.

我们先规定，$M=2K$。然后，$W_{2k}^u$是一个复数旋转因子，很学术是不是？其实它的定义是$e^{-i2\pi u/M}$，就是Fourier变换的基函数。

可以这样看：$\ast W_B^A=e^{-i2\pi A/B}\ast$，说A是旋转因子的指数，B是旋转因子的周期。我们称其为“B次单位根”。

在开始之前，我们应该介绍一下单位根的性质：

1. 周期性。那当然了，怎么想都是。$\ast W_N^{k+N}=W_N^k\ast$。
2. 对称性。$\ast W_N^{k+N/2}=-W_N^k\ast$。
3. 有 $\ast W_N^{mkn}=W_{N/m}^{kn}\ast$，如果m是N的因子的话。

现在，为了运用这个递推公式，我们首先要定义$F_{even}(u)$和$F_{odd}(u)$：

$$\begin{array}{rl}& F_{\mathrm{even}}(u)=\frac{1}{K}\sum _{x=0}^{K-1}f(2x)W_K^{ux}\\ & F_{\mathrm{odd}}(u)=\frac{1}{K}\sum _{x=0}^{K-1}f(2x+1)W_K^{ux}\end{array}$$

很显然这里的u是0,1,2,⋯,K−1。这里的Feven(u)和Fodd(u)分别是f(x)中偶数项和奇数项的 DFT.

这就是说,前K个F(u)可以通过奇部和偶部之和得到:

$$F(u)=\frac{1}{2}\left[F_{\text{even}}(u)+F_{\text{odd}}(u)W_{2k}^u\right]$$

这就是说，后k个F(u+k)可以通过奇部和偶部之差得到:

$$F(u+k)=\frac{1}{2}\left[F_{\text{even}}(u)-F_{\text{odd}}(u)W_{2k}^u\right]$$

Cooley-Tukey 算法的特性：

• 一个 M 点的 DFT 可以分解为两个 M/2 点的 DFT
• 通过计算两个 M/2 点的 DFT，得到Feven​(u) 和 Fodd​(u)
• 奇部和偶部之和得到前 k 个 F(u)
• 奇部和偶部之差得到后 k 个 F(u+k)

所以总的来说，其思想是：

• 计算两个单点的 DFT，然后合并成一个双点的 DFT
• 计算两个双点的 DFT，然后合并成一个四点的 DFT
• 依此类推，在 nlogn 内做完 FFT

频域滤波

兜兜转转在原理这么久,终于来到了最想看的东西了

之所以我们能在频域滤波，那多亏了卷积定理。

在空域（或者说时域）中，滤波器本质上是在做卷积。根据卷积定理，我们知道，在频域上做乘积等价于在空域上做卷积。所以，我们可以在频域上做乘积，然后再做逆变换，就得到了在空域上做卷积的结果。

滤波这件事，我们可以表述为：

$$g(x,y)=f(x,y)\ast h(x,y)$$

其中，f(x,y)是输入图像，h(x,y)是滤波器（也就是说，一个线性不变算子），g(x,y)是输出图像

根据卷积定理，我们有:

$$G(u,v)=F(u,v)H(u,v)$$

这几个都是通过 Fourier 变换到频域中的。

所以总的来说，做频域滤波的流程就是：

1. 转换到频域：对输入图像 f(x,y) 做 Fourier 变换，得到 F(u,v)
2. 转换到频域：对滤波器 h(x,y) 做 Fourier 变换，得到 H(u,v)
3. 在频域应用转移函数：对 F(u,v) 和 H(u,v) 做乘积，得到 G(u,v)
4. 转换回空域：对 G(u,v) 做 Fourier 逆变换，得到 g(x,y)

举个例子来说:

$$h(x)=\sqrt{2\pi }\sigma A{e}^{-2{\pi }^2{\sigma }^2{x}^2}$$

这是一维高斯滤波的时域表示,对其进行傅里叶变换:

$$H(u)=A{e}^{-u^2/2{\sigma }^2}$$

而高通滤波器的 Fourier 变换为:

$$\begin{array}{rl}h(x)& =\sqrt{2\pi }{\sigma }_{1}A{e}^{-2{\pi }^2{\sigma }_1^2{x}^2}-\sqrt{2\pi }{\sigma }_{2}B{e}^{-2{\pi }^2{\sigma }_2^2{x}^2}\\ & H(u)=A{e}^{-u^2/2{\sigma }_1^2}-B{e}^{-u^2/2{\sigma }_2^2}\end{array}$$

其中，σ1和σ2分别是两个 Gaussian 函数的标准差，并且σ1≥σ2,A≥B

低通滤波器

简而言之，留低频（低频通过，称为低通），滤高频。效果是使图像变得模糊（平滑）。

ILPF

理想低通滤波器（Ideal Low Pass Filter）的转移函数为:

$$H(u,v)=\{\begin{array}{llc}1& \mathrm{if}D(u,v)\le D_0& \\ 0& \mathrm{if}D(u,v)>D_0& \end{array}$$

这就是说，以截频D0为半径的圆内的所有频率都能无损通过，而以截频D0为半径的圆外的所有频率都被完全阻断

振铃效应

这就自然引出了振铃效应。振铃效应是指，当我们使用 ILPF 时，图像的边缘会出现一系列明暗相间的环带。

原因：

输入一个尖峰信号，相当于一个快速的输入变化，则输出会有明显的振荡。因为 Fourier 变换的频谱是周期性的，所以这种振荡会在频谱中反映出来。

环的半径反比于截止频率D0，也就是说：

• D0较小，h(x,y) 产生数量少但是较宽的环
• D0较大，h(x,y) 产生数量多但是较窄的环

Butterworth 低通滤波器

Butterworth 低通滤波器的转移函数为：

$$H(u,v)=\frac{1}{1+{\left[\frac{D(u,v)}{D_0}\right]}^{2n}}$$

其中，n是滤波器的阶数，D0是截止频率。

Butterworth 低通滤波器又称为最大平坦滤波器。和 ILPF 不同，它的通带和阻带之间没有明显的不连续性。也就是说，在通带和阻带之间有一个平滑的过渡带。

通常把 H(u,v)下降到某一值的那点定位为截止频率D0。譬如说罢，当D(u,v)=D0时，H(u,v)下降到1/2，那么我们就说D0是 H(u,v)的截止频率。

同 ILPF 的结果相比，Butterworth 低通滤波器的结果更加平滑（因此模糊程度减少），没有振铃效应；相应的，因为H(u,v)的下降速度更慢，它的尾部有不少高频成分.

Gaussian 低通滤波器

Gaussian 低通滤波器的转移函数为:

$$H(u,v)=e^{-D^2(u,v)/2D_0^2}$$

其中，D0是截止频率

高通滤波器

简而言之，留高频（高频通过，称为高通），滤低频。效果是使图像变得锐化.一般的,在频域中所有的频率减去低通的就是高通:

$$H_{\mathrm{HP}}(u,v)=1-H_{\mathrm{LP}}(u,v)$$

IHPF

$$H(u,v)=\{\begin{array}{llc}0& \mathrm{if}D(u,v)\le D_0& \\ 1& \mathrm{if}D(u,v)>D_0& \end{array}$$

即,仅保留远离中心的值.当然了,IHPF 因为生生截断，所以也会产生振铃效应.

Butterworth 高通滤波器

引入高频强调：在高通滤波器中添加一个常数以保留低频分量.

$$H(u,v)=\frac{1}{1+{\left[\frac{D_0}{D(u,v)}\right]}^{2n}}$$

同样的，$H_{BHPF}=1-H_{BLPF}$。

该滤波器没有明显的不连续性，在通过频率和滤波频率之间建立了清晰的截止点。这一特性和 BLPF 是一致的。

定义截止频率：在 H(u,v)H(u,v) 降至其最大值的特定数值处，约定为截止频率。

Gaussian 高通滤波器

$$H(u,v)=1-e^{-D^2(u,v)/2D_0^2}$$

同样的，$H_{GHPF}=1-H_{GLPF}$

同态滤波器

同态滤波是一种在频域中同时将图像亮度范围进行压缩和将图像对比度进行增强的方法，即把频率过滤和灰度变换结合起来。

滤波器性质证明

这个应该是必考的点之一了:如何利用傅里叶变换证明一个滤波器是高通还是低通?

很多空域滤波器确实一眼看不出来是高通滤波还是低通滤波(个人而言第一眼确实没看出来微分滤波器是高通还是低通)

空域下对高通/低通的描述就很泛泛:增强(锐化)边缘/模糊图像,但是在频域下就有很明确的数学特性:

• 低通滤波器:在低频区域(接近原点)具有较高的增益(函数值)
• 高通滤波器:在高频区域(远离原点)具有较高的增益(函数值)

这样就把一个模糊的问题具体化成了证明函数增减的问题了.

用我一开始没看明白的微分滤波器举个例子:

$$h(x,y)=\left[\begin{array}{ccc}-1& 0& 1\\ -2& 0& 2\\ -1& 0& 1\end{array}\right]$$

还记得这个吗.这个是sobel算子,对其应用离散傅里叶变换与欧拉公式后应该有个对应的频域函数$H(x,y)$,如何知道这个函数的增益呢?我们需要计算频域函数的频谱(幅度):

• 如果$|H(x,y)|$在高频区域比较大,在低频区域比较小则为高通
• 如果$|H(x,y)|$在低频区域比较大,在高频区域比较小则为低通

一般来讲,应用了欧拉公式的频域函数幅度都可以进行三角恒等变换为一个三角函数,这样甚至都不用求函数增减,可以直接根据图像获得结论

Chapter 5:Restoration

典型的图像恢复是根据图像退化的先验知识建立一个退化模型，以此模型为基础，采用各种逆退化处理方法进行恢复，得到质量改善的图像。

过程：找退化原因 —— 建立退化模型 —— 反向推演 —— 恢复图像

图像恢复和图像增强

相同之处：改进输入图像的视觉质量

不同之处：

• 图像增强借助人的视觉系统特性，以取得较好的视觉结果（不考虑退化原因）
• 图像恢复根据相应的退化模型和知识重建恢复原始的图像（考虑退化原因）

图像恢复的问题描述 (Problem Formulation)：

物体的理想图像应为$f(x,y)$，但由于成像系统、采集系统不理想或成像环境不理想等原因，实际得到的是退化（变质）图像$g(x,y)$。

图像恢复则是根据退化图像$g(x,y)$和退化模型$H(x,y)$，重建出理想图像$\hat{f}(x,y)$。

数学描述

空域退化模型:

$$g(x,y)=f(x,y)\ast h(x,y)+\eta (x,y)$$

频域退化模型:

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

只有噪声的空域恢复

均值滤波

算数均值（虽能降噪，但也模糊了图像）：

$$\hat{f}(x,y)=\frac{1}{mn}\sum _{s=-a}^a\sum _{t=-b}^{b}g(x+s,y+t)$$

几何均值（和算数均值差不多平滑，但丢失更少细节）：

$$\hat{f}(x,y)={\left[\prod _{s=-a}^a\prod _{t=-b}^{b}g(x+s,y+t)\right]}^{\frac{1}{mn}}$$

谐波均值（对盐噪声效果更好，但不适用于胡椒噪声，并且善于处理 Gaussian 噪声等其他噪声）：

$$\hat{f}(x,y)=\frac{mn}{{\sum }_{s=-a}^a{\sum }_{t=-b}^b\frac{1}{g(x+s,y+t)}}$$

逆谐波均值（适用于椒盐噪声，但不能同时消除）：

$$\hat{f}(x,y)=\frac{{\sum }_{s=-a}^a{\sum }_{t=-b}^{b}g(x+s,y+t{)}^{Q+1}}{{\sum }_{s=-a}^a{\sum }_{t=-b}^{b}g(x+s,y+t{)}^Q}$$

其中，Q能改变滤波器的性质。当Q>0时，滤波器对胡椒噪声更敏感；当Q<0时，滤波器对盐噪声更敏感。当Q=0时，退化为算数均值滤波器；当Q=−1时，退化为谐波均值滤波器。

统计排序（中值滤波）：

$$\hat{f}(x,y)=\text{median}\left\{g(x+s,y+t)\right\}$$

其中,$s=-a,\cdots ,a$,$t=-b,\cdots ,b$。中值滤波对多种随机噪声有良好的去噪能力，引起的模糊更少，尤其对脉冲噪声非常有效。

最大（小）值滤波：

$$\begin{gathered} \hat{f}(x,y)=\max \left\{g(x+s,y+t)\right\} \\ \hat{f}(x,y)=\min \left\{g(x+s,y+t)\right\} \end{gathered}$$

最大值滤波可发现图像中的亮点，消除胡椒噪声；最小值滤波可发现图像中的暗点，消除盐噪声。

中点滤波器：

$$\hat{f}(x,y)=\frac{1}{2}\left[\max \left\{g(x+s,y+t)\right\}+\min \left\{g(x+s,y+t)\right\}\right]$$

中点滤波器结合了统计排序和求平均操作，使得其对于高斯和均匀随机分布的噪声有较好的效果。

Alpha-Trimmed 均值滤波器：

$$\hat{f}(x,y)=\frac{1}{mn-d}\sum _{s=-a}^a\sum _{t=-b}^b{g}_r(x+s,y+t)$$

其中，d可取0,1,⋯,mn−1，r=(d−1)/2。当d=0时，退化为算数均值滤波器；当d=mn−1时，退化为中值滤波器。取其他值时，适合包含多种噪声的图像处理，如高斯和椒盐噪声的混合情况

自适应中值滤波

中值滤波器的尺寸如果太小，可以较好地保护图像的某些细节，但往往遗漏噪声；多次滤波造成特征丧失、图像模糊。反之，如果太大的话，可以加强噪声抑制能力，但图像会很模糊，有时候会滤去图像中的一些细线、尖锐边角等重要细节，从而破坏图像的几何结构。

因此，先采用小尺寸的中值滤波；当发现滤波器不够好时，扩大其尺寸，再次滤波；再不好；再扩大，直到达到最大尺寸。此谓之“自适应”。

AMF 的三个主要目的：

1. 去除椒盐(脉冲)噪声

2. 平滑其他非脉冲噪声

3. 减少诸如物体边界细化或粗化等失真

   

这个要看一下伪代码:

def adaptiveMedianFilter(img, int row, int col, int kernelSize, int maxSize)
{
		#对窗口图像排序
		window = getwindow(img,kernelSize)
		sort(window.begin(),window.end())
    auto min = pixels[0];
    auto max = pixels[kernelSize*kernelSize - 1];
    auto med = pixels[kernelSize*kernelSize / 2];
    auto zxy = img.at<uchar>(row, col);
    if (med > min && med < max)#看中值是不是噪声,是则出错,扩大窗口大小(level A)
    {
        #先看被滤波的像素是不是噪声,要不要滤掉
        if (zxy > min && zxy < max)#level B
            return zxy;
        else
            return med;
    }
    else
    {
        kernelSize += 2;
        if (kernelSize <= maxSize)
            return adaptiveMedianFilter(img, row, col, kernelSize, maxSize);#增大窗口尺寸，继续A过程。
        else
            return med;
    }
}

逆向滤波

如果噪声可忽略，并且退化函数已知或已经获得，那么我们可以直接将退化函数除掉：

$$\hat{F}(u,v)=\frac{G(u,v)}{H(u,v)}$$

我们把$\frac{1}{H(u,v)}$称为逆滤波器。

逆滤波复原过程：

1.对退化图像$g(x,y)$做Fourier变换，得到$G(u,v)$

2.计算系统点扩散函数$h(x,y)$的Fourier变换，得到$H(u,v)$

1. 计算逆滤波$F(u,v)=\frac{G(u,v)}{H(u,v)}$

4.对$F(u,v)$做Fourier逆变换，得到$\hat{f}(x,y)$

由于$G(u,v)=H(u,v)F(u,v)+N(u,v)$,所以：$\hat{F}(u,v)=F(u,v)+\frac{N(u,v)}{H(u,v)}$。因此，当$H(u,v)$的值很小的时候，噪声将占据主导地位。

倘若不存在噪声，那么逆滤波就是一个完美的恢复过程。但是，由于噪声的存在，逆滤波很可能会放大噪声，导致结果很差。除非我们知道噪声的分布情况(事实上，这也很难知道),逆滤波几乎不可用。

解决策略：忽略掉$H(u,v)$中的小值，或者用一个截止频率$D_0$来限制$H(u,v)$的值。

$$M(u,v)=\{\begin{array}{ll}1/H(u,v)& \text{if}{u}^2+v^2\le w^2\\ 1& \text{if}{u}^2+v^2>w^2\end{array}$$

Chapter 6：Morphological

形态学,某种程度上感觉有点像人的减肥(

基础运算

1.膨胀

膨胀（Dilation）是数学形态学中的一种基本操作，用于对图像或集合进行形态学变换。在二值图像中，膨胀可以扩展或增大图像中的对象，使其更加粗壮或更接近原始形状.

$$A\oplus B=\left\{z\mid (\hat{B}{)}_z\cap A\subseteq A\right\}$$

B 的反射进行平移与A的交集是 A 的子集。在给定的定义中，$(\hat{B}{)}_z$表示结构元素B的反射 (reflection)进行平移操作。这个操作将结构元素B沿着水平和垂直方向翻转，并将其在空间中进行平移，使得结构元素的原点(通常是结构元素的中心)与元素$z$对齐

平移操作是指将一个对象在平面上移动一定的距离，保持其形状和方向不变。在这种情况下，结构元素B的反射进行平移是指对结构元素B进行翻转并移动，使得结构元素的原点与元素$z$对齐.

然后，$(\hat{B}{)}_z\cap A$表示结构元素B的反射进行平移后与原始集合A的交集。这个交集表示了在结构元素B的形状和位置经过反射平移后，与原始集合A重叠的部分.

最后，定义中指出如果$(\hat{B}{)}_z\cap A$是$A的子集，即所有在结构元素B的反射平移后与A重叠的部分都包含在A中，那么元素$z$就被包含在膨胀操作的结果中

这个定义的意思是，膨胀操作将结构元素B的反射平移到原始集合A上，然后检查反射平移后的结构元素与A的交集。如果这个交集的所有元素都包含在A中，那么元素$z$就被包含在膨胀操作的结果中。换句话说，膨胀操作通过将结构元素B的形状扩展到原始集合A上来增大集合的大小。

对于这种全1结构元,膨胀结果可以简称为”胖一圈”,就像过完年的你

2.腐蚀

腐蚀(Erosion):$A\ominus B=\{z\mid (B{)}_z\subseteq A\}$

腐蚀 (Erosion) 是数学形态学中的另一种基本操作，用于对图像或集合进行形态学变换。腐蚀操作可以使图像中的对象变得更加细小或更接近原始形状。

在给定的定义中，A和B是两个集合，其中A是原始集合，B是结构元素(也称为核)。腐蚀操作$A\ominus B$的结果是一个新的集合，其中包含满足以下条件的所有元素$z$:

1. 将结构元素B的原点(通常是结构元素的中心)与元素$z$对齐。
2. 检查结构元素B的所有非零元素是否都包含在原始集合A中。
3. 如果结构元素B的所有非零元素都包含在原始集合A中，则元素z被包含在腐蚀操作的结果中。

换句话说，腐蚀操作将结构元素B在原始集合A上滑动，并检查每个位置上的包含关系。如果结构元素B的所有非零元素都被包含在原始集合A中，那么该位置就被包含在腐蚀操作的结果中。

腐蚀操作可以用于图像处理中的多种应用，例如去除图像中的小噪点、分离相连的对象或缩小对象的边界。它也可以与其他形态学操作 (如膨胀、开运算和闭运算)结合使用，以实现更复杂的图像处理任务。

3.开/闭运算

开运算 (Opening):$A\circ B=(A\ominus B)\oplus B$ 在不改变形状的前提下，使图像的轮廓变得光滑，能够断开狭窄的间断，消除细的突出物。 开操作的性质：

• 平移。$O(A+x,B)=O(A,B)+x$
• 反广延性(Antiextensivity)。$O(A,B)\subseteq A$
• 增加单调性 (Increasing Monotonicity).

$$A_1\subseteq A_2\Rightarrow (A_1\circ B)\subseteq (A_2\circ B)$$

• 幂等性 (ldempotence)。$A\circ B\circ B=A\circ B$

开运算从几何意义角度上看，可以说是结构元在目标图像边界内滑动时，所能到达的最远点

闭运算 (Closing): $A\bullet B=(A\oplus B)\ominus B$ 在不明显改变面积前提下，使图像的轮廓变得光滑。可以弥合狭窄的间断和细长的鸿沟，消除小的孔洞，填补轮廓线中的裂痕。闭运算的性质：

• 平移。$C(A+x,B)=C(A,B)+x$
• 广延性 (Extensivity)。$A\subseteq C(A,B)$
• 增加单调性 (Increasing Monotonicity)。$A_1\subseteq A_2\Rightarrow (A_1\bullet B)\subseteq (A_2\bullet B)$
• 幂等性 (ldempotence)。$A\bullet B\bullet B=A\bullet B$

同样的，闭运算的几何意义是结构元在目标图像边界外滑动，所能到达的最远位置

边界提取

边界提取 (Boundary Extraction):$\beta (A)=A-(A\ominus B)$

把腐蚀后的图像从原图像中减去，得到的就是边界叻

提取连通分量

令Y表示一个包含于集合A中的连通分量，并假设Y中的一个点p是已知的。用下列迭代式生成Y的所有元素：

$$Xk=Xk-1\oplus B\cap A$$

其中，$X_0=p，k=1,2,\cdots$。当$X_k=X_{k-1}$时，迭代停止，$X_k$就是Y的所有元素.

啥意思呢,其实用代码的形式看更加简单明了:

def morphological_region_growing(image, seed_point, kernel):
    rows, cols = image.shape
    mask = np.zeros_like(image)
    mask[seed_point] = 255
 
    prev_mask = np.zeros_like(mask)
    
    while not np.array_equal(mask, prev_mask):
        prev_mask = mask.copy()
        
        # 膨胀
        dilated = dilate(mask, kernel)
        
        # 与原图像相交
        mask = np.logical_and(dilated, image).astype(np.uint8) * 255
 
    return mask

• 在每一轮中prev_mask是上一轮留下的mask
• 对mask进行一次膨胀操作得到dilated，最终mask为膨胀和原图相交后的结果，这自然就有两种情况：
  1. mask比原图连通区域小，相交结果为mask
  2. mask比原图连通区域大，相交结果为原图连通空间
• while比较上一轮mask和上上轮mask（prev_mask）是否相同，很明显的，如果mask最后覆盖了连通区域，那根据第二步求交，mask就不会再变保持连通区域大小，这样mask和prev_mask就相同了