2021年10月21日10:24:05
占坑，这周五图像处理课讲完PPT后填。
2021年10月22日21:07:55 -> 讲课完毕

数字图像处理、图像恢复，噪声去除PPT

图像处理课课件，包含图像恢复、噪声去除部分内容及相关注解
链接: https://pan.baidu.com/s/1wBfXM0ZrSw3O221GUShVtw  密码: 424w
--来自百度网盘超级会员V5的分享

分割线

说明：在准备课件上台讲解的时候，为去除周期性噪声，需在频域上进行操作；所以在网上查看了许多视频博客知乎等，有了初步的理解，现在将收集的资料整理出来，并发布到博客，在网上参考的地方，我会一一给出参考，方便其它同学快速理解将图像从空域到频域转换的这个过程，也提供一些比较好的博文及视频，希望能对你有参考作用！

以下内容来自BILIBILI的一个非常棒的视频：《P04 图像的快速傅里叶变换及其应用，作者：flybirdAI》

说明：图像和公式来自此视频中的截图，此视频讲的非常的好，清晰明了，对于空域到频域的转换，不涉及代码的话看前28分钟就足够了！

频域变换一般过程：

原图->DFT/FFT(正变换)->中心化->频域显示(处理)->去中心化->IDFT/IFFT(反变换)->原图

图像高频：亮度或灰度变化激烈的部分，如边缘和轮廓
图像低频：变化不大或变化较为平缓的地方，如大块的结构

DFT（离散傅里叶变换）：Discrete Fourier transform
FFT（快速傅里叶变换）：Fast Fourier transform

二维图像（M行N列）
正变换：
$$F(u,v)=\sum _{x=0}^{M-1}\sum _{y=0}^{N-1}f(x,y)e^{-j2\pi (ux/M+vy/N)}$$

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

其中下式为欧拉公式：
$$e^{j\theta }=\cos (\theta )+j\sin (\theta )$$

频域分析就是对F(u，v)进行操作：

由于共轭对称性：可对FFT变换的结果再进行FFT变换得到IFFT（I指的是反变换，IFFT反快速傅里叶变换，将图像从频域转换到空域）

注：傅里叶谱通过对log(1+F|(u,v)|)表示，另外自行频率分析时，也往往将频率窗口移动窗口中心处。

FFT还可以做卷积：

1. 将图像和卷积模板都进行正变换到频率域
2. 将两者在再频域直接进行相乘运算
3. 再反变换到空间域作为输出

频域变换过程时刻牢记下图：

如果深入理解后，可以观看3Blue1Brown的视频：【官方双语】形象展示傅里叶变换

代码：【原文博客：python的numpy库和cv2库实现图像傅里叶变换】

# writer:wojianxinygcl@163.com

# date  : 2020.3.30 

# 原文地址 : https://www.cnblogs.com/wojianxin/p/12530172.html


import cv2
import numpy as np
from matplotlib import pyplot as plt


def fft_cv(img):
    """
    用cv进行傅里叶变换
    :param img:
    :return:
    """

    # 傅里叶变换

    dft = cv2.dft(np.float32(img), flags=cv2.DFT_COMPLEX_OUTPUT)

    dftshift = np.fft.fftshift(dft)

    res1 = 20 * np.log(cv2.magnitude(dftshift[:, :, 0], dftshift[:, :, 1]))

    # 傅里叶逆变换

    ishift = np.fft.ifftshift(dftshift)

    iimg = cv2.idft(ishift)

    res2 = cv2.magnitude(iimg[:, :, 0], iimg[:, :, 1])

    # 显示图像

    plt.subplot(131), plt.imshow(img, 'gray'), plt.title('Original Image')

    plt.axis('off')

    plt.subplot(132), plt.imshow(res1, 'gray'), plt.title('Fourier Image')

    plt.axis('off')

    plt.subplot(133), plt.imshow(res2, 'gray'), plt.title('Inverse Fourier Image')

    plt.axis('off')

    plt.show()


def fft_numpy(img):
    """
    用numpy对图像进行傅里叶变换
    """
    # 傅里叶变换

    f = np.fft.fft2(img)

    fshift = np.fft.fftshift(f)

    res = np.log(np.abs(fshift))

    # 傅里叶逆变换

    ishift = np.fft.ifftshift(fshift)

    iimg = np.fft.ifft2(ishift)

    iimg = np.abs(iimg)

    # 展示结果

    plt.subplot(131), plt.imshow(img, 'gray'), plt.title('Original Image')

    plt.axis('off')

    plt.subplot(132), plt.imshow(res, 'gray'), plt.title('Fourier Image')

    plt.axis('off')

    plt.subplot(133), plt.imshow(iimg, 'gray'), plt.title('Inverse Fourier Image')

    plt.axis('off')

    plt.show()


if __name__ == '__main__':
    # 读取图像
    img = cv2.imread('LENA256.BMP', 0)  # 读取灰度图
    fft_numpy(img)  # 使用numpy进行变换
    # fft_cv(img)   # 使用cv进行变化

Lena标准测试图: https://pan.baidu.com/s/1FGAYFsNaPR7XuKbL6aTnhQ  密码: m19a
--来自百度网盘超级会员V5的分享

结果：

参考以下博客：

1. 【傅里叶分析之掐死教程（完整版）更新于2014.06.06】
2. 【通俗讲解：图像傅里叶变换】
3. 【傅立叶变换频谱图怎么看？cccw的问题回答】
   

我的总结：

1.图2是转换到频域二维可视化的结果，称作频谱图

2.由于频谱动态范围很大，所以用log变换降低频谱的动态范围

3.图1为空域：f(x, y) , 图2为频域：H(u, v)，中心化之后的结果，频谱图的正中心点为H（0，0），图像大小与原图相同

4.图1和图2相同位置的点没啥直接的对应关系，右边的每一个点可以理解为所描述的是一种平面波，它是构成左边图无数个波面中的一个，它与左边的每一个点都相关

一个曲面可以有很多个正弦平面波组成，二维正弦平面波，可以用以下4个参数来描述它：频率、幅度、相位、方向

所有波都可以用多个正弦波叠加；波又可以通过频率、幅值和方向来表示（频谱图中相位信息被舍弃）

5.图2中的每一个点：
它到中点的距离描述的是频率
中点到它的方向是构成平面波的方向
那一点的灰度值描述的是它的幅值

6.频谱图里面描述的信息有：频率f，振幅A，方向n

7.每一个频谱图中的点它不是实数，是一个虚数，a+bi，为二维可视化，计算的是它的绝对值

8.点点越亮代表这个数值越大，点点越暗代表这个数值越小（幅值的大小）

带通/带阻滤波器

如下图所示，左图是被周期性正弦噪声所污染的图像，右图是去噪后的结果，想要在空域去除这样的噪声是很困难的，但是在频域上确能很轻松的进行操作。

带阻滤波器可以移除图像中特定的频率范围，从而很轻松的去除周期性噪声。

简单解释带阻滤波器以及带通滤波器：

• 带阻滤波器：带阻，阻就是阻止，阻止特定的频率范围通过
• 带通滤波器：带通，通就是通过，只允许特定的频率范围通过

由上简单的定义可知，带阻滤波器和带通滤波器可以相互转换，是一种互补的关系，带阻滤波器取反就是带通滤波器，带通滤波器取反就是带阻滤波器。

理想的带阻滤波器定义如下：

$$H(u,v)=\{\begin{array}{ll}1& \text{ if }D(u,v)<D_0-\frac{W}{2}\\ 0& \text{ if }{D}_0-\frac{W}{2}\le D(u,v)\le D_0+\frac{W}{2}\\ 1& \text{ if }D(u,v)>D_0+\frac{W}{2}\end{array}$$
【参考博客】：数字图像处理与python实现-带通滤波器
带阻滤波器，如下图3（第二排第一个），理想的带阻滤波器的定义如上，D(u,v)为频域中任意一点到中心点的距离，D0代表带阻滤波器的半径（从中心到带阻滤波器第一个内圈和第一个外圈之间的中间位置），W代表带阻滤波器的带宽（内圈到外圈的距离），上述公式代表的就是带阻滤波器，含义就是带阻滤波器带宽范围内的部分：不允许它通过，带阻滤波器带宽范围外的地方保持不变，可视化的结果如下图3。

图像从上到下，从左到右分别表示图1：周期性正弦噪声污染的图像；图2：图1转换的频域的结果；图3：理想带阻滤波器；图4：在频域去除周期性噪声后转换到空域的结果

如何使用带阻滤波器

1. 将图像1转换到频域，如图2，图1周期性正弦噪声在频域可视化的结果，表现在图2上就是围绕中心一圈分布均匀亮亮的8个小点点
2. 去掉周期性噪声，就是去掉这些小点点
3. 设计合适的带阻滤波器，如图3，黑色部分设为0，白色部分设为1，与图2的大小相同，两幅图像进行哈达玛积（对应点相乘）->得到结果（带阻滤波器带宽上的点都被置为0了，周期性噪声点在带宽内也被置为0，消除掉了）
4. 然后将第三步去噪后的结果，通过IFFT反快速傅里叶变换还原到空域，得到最终结果图4

1. 傅里叶变换理论层面理解

数学意义:傅里叶变换将一个任意的周期函数分解成为无穷个正弦函数的和的形式
物理效果:傅里叶变换实现了将信号从空间域到频率域的转换
信号分析: 一维傅里叶变换（将杂乱的信号由时域转化到频域中）一维傅里叶变化是将信号分解为正弦波的和的形式

时域
横轴是时间,纵轴是振幅

频域
横轴是频率,纵轴是振幅

对一个信号做时域到频域的变换，能够清除看到该信号主要由两个正弦波以及一些噪声混合，如下图所示：

将有效的信号频率提取并分离拟合出信号信息，然后将噪声过滤掉，得到滤波结果

频谱图: 二维傅里叶变换 (原图中的像素值是x,y坐标轴下的（即空间域），而傅立叶变换后的像素值是u,v坐标轴下的（即频域)),二维傅里叶变换可以将一个二维信号（图像）分解为三角平面波之和的形式

二维信号的离散傅里叶变换所得到的结果的频率成分的分布如下所示：
在经过频谱居中后的频谱中，中间最亮的点是最低频率，属于直流分量（DC分量）。越往边外走，频率越高。所以，频谱图中的四个角和X,Y轴的尽头都是高频(未经过频谱居中的相反）

二维离散傅里叶变换的作用：可以将空间域（二维灰度数表）的图像转换到频域（频率数表），使得更直观地观察和处理图像，也更有利于进行频域滤波等操作。

一副图像（不论是灰度的图像还是彩色图像）所提供的信息都是显而易见的（排除一些本身就很抽象的图片，不搞艺术哈哈哈）。然而，一副图像的傅里叶频谱图，却常常让人难以理解，满脸问号，一脸蒙B。

2. 常见频谱图

单色图
只有一种颜色

正弦波

分析：频率越低对称的频点越靠近频谱中心，随着频率增加远离中心（因为中心是低频，两端为高频）

灰度变化

分析：图像灰度变化方向体现在频谱图的方向中

根据上面两个信息可以得到下面结果：

常见几何形状

分析：还是和上面讲解的一样，频谱图的绘制规律是朝着梯度变化的方向改变的

这边可以这样理解傅里叶变换：
第一种理解: 二维图像进行傅里叶变化得到的频谱图(图像梯度的分布图)，当然频谱图上的各点与原图像各点并不存在一一对应关系，即使在不移频的情况下也不对应。傅里叶频谱图上看到的明亮不一的亮点，实际上图像上某一点的像素灰度值与它的邻域点差异的强弱，即梯度的大小，也即该点的频率的大小（也可以这么理解，低频部分指低梯度的点，高频部分指高梯度的点）。
另一种理解：图像二维频谱图通过对输入图像进行水平和竖直两个方向的所有扫描线的一维傅立叶变换进行叠加得到,用来表示输入图像的频率分布。
频谱图以图像的中心为圆心,圆的相位对应原图中频率分量的相位,半径对应频率高低。低频部分半径小,高频对应的半径较大,中心为直流分量（图像灰度的平均值）,某点的灰度值对应该频率的能量高低。

3. 傅里叶变换在图像中的应用

Opencv中的傅里叶变换实现
实现步骤：

• getOptimalDFTSize()函数得到DFT变换后结果的最优尺寸大小
• 根据得到的尺寸大小，使用copyMakeBorder()函数填充图像，得到填充后的Mat
• 根据新生成的Mat，使用merge()函数得到一个双通道的Mat，命名为planes
• 使用dft()函数进行傅里叶变换，得到通道1为实部，通道2为虚部

实现代码如下所示：

image.convertTo(image, CV_32F);
	vector<Mat> channels;
	split(image, channels);  //分离RGB通道
	Mat image_B = channels[0];
	//选取最适合做fft的宽和高
	int m1 = getOptimalDFTSize(image_B.rows);
	int n1 = getOptimalDFTSize(image_B.cols);
	Mat padded;
	//填充
	copyMakeBorder(image_B, padded, 0, m1 - image_B.rows, 0, n1 - image_B.cols, BORDER_CONSTANT, Scalar::all(0));
	Mat planes[] = { Mat_<float>(padded), Mat::zeros(padded.size(), CV_32F) };
	Mat complexI;
	merge(planes, 2, complexI);  //planes[0], planes[1]是实部和虚部

	dft(complexI, complexI, DFT_SCALE | DFT_COMPLEX_OUTPUT);
	split(complexI, planes);

	//定义幅度谱和相位谱
	Mat ph, mag, idft;
	phase(planes[0], planes[1], ph);
	magnitude(planes[0], planes[1], mag);  //由实部planes[0]和虚部planes[1]得到幅度谱mag和相位谱ph

	//重新排列傅里叶图像中的象限，使得原点位于图像中心
	int cx = mag.cols / 2;
	int cy = mag.rows / 2;
	Mat q0(mag, Rect(0, 0, cx, cy));       //左上角图像划定ROI区域
	Mat q1(mag, Rect(cx, 0, cx, cy));      //右上角图像
	Mat q2(mag, Rect(0, cy, cx, cy));      //左下角图像
	Mat q3(mag, Rect(cx, cy, cx, cy));     //右下角图像

	//变换左上角和右下角象限
	Mat tmp;
	q0.copyTo(tmp);
	q3.copyTo(q0);
	tmp.copyTo(q3);

	//变换右上角和左下角象限
	q1.copyTo(tmp);
	q2.copyTo(q1);
	tmp.copyTo(q2);

	imshow("mag", mag);

实现结果如下所示：

傅里叶逆变换：

//傅里叶逆变换
	polarToCart(mag, ph, planes[0], planes[1]);  //由幅度谱mag和相位谱ph恢复实部planes[0]和虚部planes[1]
	merge(planes, 2, idft);
	dft(idft, idft, DFT_INVERSE | DFT_REAL_OUTPUT);
	image_B = idft(Rect(0, 0, image.cols & -2, image.rows & -2));
	image_B.copyTo(channels[0]);

	merge(channels, image);
	image.convertTo(image, CV_8U);

傅里叶逆变换实现结果：

低通滤波
实现机制：将高频的信号舍去；(通低频，阻高频）

	//低通滤波
 	for (int i = 0; i < mag.cols; i++){
 		for (int j = 0; j < mag.rows; j++){
 			if (abs(i - mag.cols / 2) > mag.cols / 10 || abs(j - mag.rows / 2) > mag.rows / 10)
 				mag.at<float>(j, i) = 0;
 		}
 	}

低通滤波实现结果：

高通滤波
实现机制：将低频的信号舍去；(通高频，阻低频）

	//高通滤波
 	for (int i = 0; i < mag.cols; i++){
 		for (int j = 0; j < mag.rows; j++){
 			if (abs(i - mag.cols / 2) < mag.cols / 30 && abs(j - mag.rows / 2) < mag.rows / 30)
 				mag.at<float>(j, i) = 0;
 		}
 	}

高通滤波实现结果