python的opencv操作记录(五) - 空间域与频域转换

​ 在各种数字图像处理的资料中，频繁出现如下的几个概念：

• 空间域
• 频域
• 时域
• 傅里叶变换 & 傅里叶级数
• DCT
• DFT

​ 特别是几个概念搅和在一起，分分钟劝退。琢磨了半天，在在这个系列中插播一点理论方面的理解。还是那句老话，我是用比较通俗的话说出我自己对这些概念的理解，不一定很严谨，但是基本正确，能更好的理解。毕竟我们是程序员，不是专门做数学研究的。

什么是域

​ 域是一个数学上的概念。还记得数学中的定义域和值域么？就是表示变量的数字在哪个范围里面。在数学中，一般表示某种映射关系的话，都可以使用X-Y坐标系来表示，X轴(横轴)表示自变量，Y轴(数轴)表示因变量。

​ 表示某种数字信号或者函数与其因变量之间的映射关系。

空间域

​ 空间域就是说，自变量是图像中的像素点位置，而值域就是针对这个位置所标记的像素的值进行一些操作和变换。这个就是空间域，这些操作也就叫做空间域操作。
$$y=f(x,y)$$

​ 其中的x, y就表示图像中像素的x坐标和y坐标，y就表示这个像素的像素值，f表示这样一种映射关系。因为这个映射关系的自变量是x和y坐标，这是一个空间上的概念，所以这样一种映射关系就可以称作图像的空间域表达。

频域

​ 如果把图像中的一行抽取出来，比如下面这张图(512 * 512的尺寸)：

​ 图一

​ 我们提取出其中的第100行数据：

img = cv2.imread("e:\\lena.jpg", 0)
line = img[100]
print(line)
1
2
3

​ 输出的结果是这样的一个数组：

​ [ 93 88 85 85 85 84 88 92 92 92 92 92 92 92 92 92 92 92 …]

​ 用plt画出来的曲线是这样的：

​ 图二

​ 这个图像实际上可以描述这一行图像像素数据的灰度变换程度，某种程度上就可以被成为图像的频率，针对图像这样频率的处理的方式就可以称为频率域。

空间域到频域的转换

DCT(离散余弦转换)

​ 一般的空间域到频域的转换都是从DFT，也就是离散傅里叶变换开始说起，但是我觉得DCT更好理解一点，所以先说说DCT。

DCT的基本逻辑：

1. 和频率域的基本逻辑一样，把图像中的某一行拿出来，其中的每一个像素点可以被描述成(x, y) -> 像素值的映射关系，其中x还是固定的。如果用像素的y做横坐标，像素值做纵坐标，可以形成一个二维坐标系。这些像素数据或者形成的二维曲线可以被成为原始信号。

   2. 根据傅里叶变换的概念，任何函数曲线都可以用正弦或者余弦去模拟，那么上面这个图二肯定也是可以的。

   3. 傅里叶变化是可以用很多个cos或者sin函数去逼近的，那么这里用多少个呢。因为上面的曲线是由若干个点组成的，可以理解是对这个未知曲线做的一个采样。

   4. 而dct中的逻辑是使用采样点个cos余弦函数去逼近，也就是这个函数被认为是

   $$f(x)=X_{0}cos(0)+X_{1}cos(x)+X_{2}cos(2x)+...+X_{N}cos(Nx)$$

   要想确定这个函数，就是要找出这N个权重系数X，用N个余弦波是为了后面矩阵计算的方便与可行。

   5. 关于上面这一条，这里我看了一个视频，我对这个地方的理解是，第k个样本点可以替代频率为kx的余弦函数在这个函数中占的比例，总共N个样本点，所以是N个余弦函数。

   放上视频的链接：

   6. 使用下面的dct公式，来通过样本点计算这些权重系数X。

   $$X_k=\sum _{n=0}^{N-1}{x}_{n}cos[\frac{(2n+1)\pi k}{2N}]$$

   ​ 公式一

   7. $x_n$可以称作原始信号采样点，而cos可以被看做是这个样本点代表余弦频率的余弦波的采样点。如果把上面的公式一写成向量形式：

   $$X=Cx$$

   实际上就是表示第k个样本点所代表的kx的余弦波的N个采样点的值向量C，然后由N个这样的行向量组成了矩阵C，这个矩阵与原始信号向量x进行点乘，形成最后的权重向量X。

   可以参考下面两张图。

简单一维数组的DCT DEMO

假设有一个一维数组：(1, 1, 1, 1, 1)

那么根据dct的公式：

因为k=0，所以代入dct公式中，cos部分都等于1了，所以结果就是5。
$$X_0=1+1+1+1+1=5$$
​ 因为所有的$x_n$都是1，所以公式就只剩下余弦部分了。
$$\begin{gathered} X_1=cos[\frac{(2\ast 0+1)\ast 1\pi }{10}]+cos[\frac{(2\ast 1+1)\ast 1\pi }{10}] \\ +cos[\frac{(2\ast 2+1)\ast 1\pi }{10}]+cos[\frac{(2\ast 3+1)\ast 1\pi }{10}]+cos[\frac{(2\ast 4+1)\ast 1\pi }{10}] \\ =cos(\frac{\pi }{10})+cos(\frac{3\pi }{10})+cos(\frac{\pi }{2})+cos(\frac{7\pi }{10})+cos(\frac{9\pi }{10}) \\ =0 \end{gathered}$$
​ 其他三个X大家有兴趣也可以算一算，都是等于0。所以输出就是（5,0,0,0,0），当然opencv输出的shape是一个(5, 0)的数组，而且opencv在公式一的基础上还增加了一点东西，矩阵$C$加了一个前缀：
$$\begin{gathered} C_{ij}^{(N)}=\sqrt{{\alpha }_j/N}cos(\frac{\pi (2k+1)j}{2N}) \\ {\alpha }_0=1,{\alpha }_j=2,j>0 \end{gathered}$$
​ 所以输出值不是5，没有去细究了，大体的逻辑肯定是一样的。

line = (1, 1, 1, 1, 1)
dft = cv2.dft(np.float32(line))

print(dft.shape)
print(dft)
1
2
3
4
5

​ 输出的dft如下：

• 物理含义，假设(1, 1, 1, 1, 1)是像素值，那么这张图片的像素之间是没有任何变化的，也可以理解没有频率，和直流电一样，没有变化，最终变换到频率域之后，就只有常量$X_0$上有作用，很好的表示了图像在频率上的描述。
• 而且常量值可以很好的表达像素值，如果像素值比较大，这个常量值$X_0$也会比较大。
• 在上面的向量公式中，矩阵C中的不同行向量之间是互相正交的，也就是点积为0。在这个小例子中，原始信号就是直流信号，只有$X_0$会对其影响，所以其他的行向量进行点乘后，结果都为0。

换一组频率变化大和变化小的一维数组

​ 换一个一维数组：(1, 100, 200, 2, 4, 250, 23, 198, 2, 251)。

​ 这个数组的输出值为：

​ 再看看变化小的一个数组：(1, 2, 3, 4, 5, 6, 7, 8, 9, 10)

​ 从上面两个数组的dct输出可以看出，频率变化大和变化小的像素信号经过DCT转换后由很大的区别，那么就很好对这样的一些信号去进行区分和处理了。

DCT逆变换

​ 根据上面提到的DCT的向量表达信息，实际上只需要求出矩阵$C$的逆矩阵$C^{-1}$，再用这个矩阵去乘以权重向量X就可以得到原来的数组了。
$$x=C^{-1}X$$

demo_line=(1, 2, 3, 4, 5, 6, 7, 8, 9, 10)
dct_demo = cv2.dct(np.float32(demo_line))
iLine = cv2.idct(dct_demo)
print(iLine)
1
2
3
4

​ 就可以得到原来的数据了。

二维DCT

​ 从一维往二维拓展一下，比如一个8 * 8的矩阵，先对每一行做一次DCT操作，最终可以得到一个 8 * 8的DCT矩阵，然后再针对每一列做一次DCT，还是获得一个8 * 8的DCT矩阵。如下图所示。

真实图像

​ 我针对上面的图一用opencv的dct函数做了一次转换：

dct_img = cv2.imread("e:\\lena.jpg", 0)
dct = cv2.dct(np.float32(dct_img), flags=cv2.DFT_COMPLEX_OUTPUT)

cv2.imwrite("e:\\aaa.jpg", dct)
1
2
3
4

​ 打开图片一看，是这么个东西。

​ 我的理解如下：回忆一下二维DCT的过程，在X方向上，从左往右的数据分别是
$$cos(0x),cos(x),cos(2x),cos(3x)......cos(Nx)$$
​ 的权重系数
$$X_0,X_1,X_2......X_N$$
​ 而在图像中，只有左上角是比较黑的，所以可以说这幅图像在X方向上的频率(也就是像素值的变化率)集中在低频区。换言之，在Y轴方向上的频率也是集中在低频区。这样就可以采用一些低通，高通等不同的滤波算法进行区分和操作了。

​ 当然，我个人觉得这里的频率变换可能也是相对的，我感觉最右边的可能也到了$cos(100x)$的高频了吧。

DFT，离散傅里叶变换

​ 在图像处理中，用的更多的是这个DFT，也就是离散傅里叶变换。应该说DCT是DFT的一种特例，DCT说完了，就再说说DFT。

​ 从DFT的公式说起：
$$X_k=\sum _{n=0}^{N-1}{x}_n{e}^{-2\pi j\frac{kn}{N}}$$

​ 然后根据欧拉公式：
$$e^{jw}=cos(w)-jsin(w)$$
​ 所以DFT公式可以展开为：
$$X_k=\sum _{n=0}^{N-1}{x}_n[cos(2\pi kn/N)-jsin(2\pi kn/N)]$$
​ 相对于DCT，DFT有两个部分，实数部分与虚数部分，其中实数部分是和DCT一样的。那么虚数部分代表啥呢，图像中的像素是一个一个真实的像素点，而且是0-255之间的一个整数。

​ 其实我的理解是在opencv中，傅里叶变换就是把这两个部分做了图像两个维度的变换，一个叫频率谱(实数部分)，一个叫相位谱(虚数部分)。不要把它和数学中的虚数一一对应，就是两个维度的数字就好。

​ 其中频率谱部分在DCT的部分已经说过了，就是代表了图像两个方向的像素变换情况。

​ 那么相位谱又代表了什么呢？我在网上找了蛮久的资料，没有找到准确合适的描述，下面这段话是我自己的理解，去理解图像的相位谱到底代表啥意思。

• 首先，图像的数据是离散的，在信号处理领域，如果信号上的点是离散的，相当于对信号进行采样：

图中有8个采样点，每个点的相位就是他在X轴上的位置，比如第一个点的相位就是$\pi /8$。

• 那么回到图像中，一行图像的数据可以被认为是某种信号，这个上面说过了。

那么这个信号可以拆成若干个频率不同的余弦波。那么我的理解是

傅里叶的虚数部分，也就是相位谱，就是这个采样点在频率为x或者2x的不同频率中的波信号中所对应的相位。

比如第一个点，可能在 cos(x)这个信号中的相位是$\pi /9$，cos(2x)这个频率的信号中的相位是$\pi /4$，类似这种的。

• 那么在图像中的物理意义，我是这么理解，频率谱是表达了像素之间的变化，而相位包含了在某个频率信号中的位置，这个位置可以得到这个信号中具体的取值，可以对应到具体的信号幅度，某种意义上来说可以代表在图像中的像素值。

实例

DFT

python代码：

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

# 1 读取图像

img = cv.imread('e:\\lena.jpg',0)

# 2 傅里叶变换

# 2.1 正变换

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

# 2.2 频谱中心化

dft_shift = np.fft.fftshift(dft)

# 2.3 计算频谱和相位谱

mag, angle = cv.cartToPolar(dft_shift[:,:,0], dft_shift[:,:,1], angleInDegrees=True)
mag=20*np.log(mag)

# 3 傅里叶逆变换

# 3.1 反变换

img_back = cv.idft(dft)

# 3.2 计算灰度值

img_back = cv.magnitude(img_back[:,:,0],img_back[:,:,1])

# 4 图像显示

plt.figure(figsize=(10,8))
plt.subplot(221),plt.imshow(img, cmap = 'gray') 
plt.title('输入图像'), plt.xticks([]), plt.yticks([]) 
plt.subplot(222),plt.imshow(mag, cmap = 'gray')
plt.title('频谱'), plt.xticks([]), plt.yticks([]) 
plt.subplot(223),plt.imshow(angle, cmap = 'gray')
plt.title('相位谱'), plt.xticks([]), plt.yticks([])
plt.subplot(224),plt.imshow(img_back, cmap = 'gray')
plt.title('逆变换结果'), plt.xticks([]), plt.yticks([])
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45

输出图像为：