文章目录
在各种数字图像处理的资料中,频繁出现如下的几个概念:
- 空间域
- 频域
- 时域
- 傅里叶变换 & 傅里叶级数
- DCT
- DFT
特别是几个概念搅和在一起,分分钟劝退。琢磨了半天,在在这个系列中插播一点理论方面的理解。还是那句老话,我是用比较通俗的话说出我自己对这些概念的理解,不一定很严谨,但是基本正确,能更好的理解。毕竟我们是程序员,不是专门做数学研究的。
什么是域
域是一个数学上的概念。还记得数学中的定义域和值域么?就是表示变量的数字在哪个范围里面。在数学中,一般表示某种映射关系的话,都可以使用X-Y坐标系来表示,X轴(横轴)表示自变量,Y轴(数轴)表示因变量。
表示某种数字信号或者函数与其因变量之间的映射关系。
空间域
空间域就是说,自变量是图像中的像素点位置,而值域就是针对这个位置所标记的像素的值进行一些操作和变换。这个就是空间域,这些操作也就叫做空间域操作。
y
=
f
(
x
,
y
)
y = f(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)
输出的结果是这样的一个数组:
[ 93 88 85 85 85 84 88 92 92 92 92 92 92 92 92 92 92 92 …]
用plt画出来的曲线是这样的:
图二
这个图像实际上可以描述这一行图像像素数据的灰度变换程度,某种程度上就可以被成为图像的频率,针对图像这样频率的处理的方式就可以称为频率域。
空间域到频域的转换
DCT(离散余弦转换)
一般的空间域到频域的转换都是从DFT,也就是离散傅里叶变换开始说起,但是我觉得DCT更好理解一点,所以先说说DCT。
DCT的基本逻辑:
-
和频率域的基本逻辑一样,把图像中的某一行拿出来,其中的每一个像素点可以被描述成(x, y) -> 像素值的映射关系,其中x还是固定的。如果用像素的y做横坐标,像素值做纵坐标,可以形成一个二维坐标系。这些像素数据或者形成的二维曲线可以被成为原始信号。
-
根据傅里叶变换的概念,任何函数曲线都可以用正弦或者余弦去模拟,那么上面这个图二肯定也是可以的。
-
傅里叶变化是可以用很多个cos或者sin函数去逼近的,那么这里用多少个呢。因为上面的曲线是由若干个点组成的,可以理解是对这个未知曲线做的一个采样。
-
而dct中的逻辑是使用采样点个cos余弦函数去逼近,也就是这个函数被认为是
f ( x ) = X 0 c o s ( 0 ) + X 1 c o s ( x ) + X 2 c o s ( 2 x ) + . . . + X N c o s ( N x ) f(x) = X_0cos(0)+X_1cos(x)+X_2cos(2x)+...+X_Ncos(Nx) f(x)=X0cos(0)+X1cos(x)+X2cos(2x)+...+XNcos(Nx)
要想确定这个函数,就是要找出这N个权重系数X,用N个余弦波是为了后面矩阵计算的方便与可行。
- 关于上面这一条,这里我看了一个视频,我对这个地方的理解是,第k个样本点可以替代频率为kx的余弦函数在这个函数中占的比例,总共N个样本点,所以是N个余弦函数。
放上视频的链接:
- 使用下面的dct公式,来通过样本点计算这些权重系数X。
X k = ∑ n = 0 N − 1 x n c o s [ ( 2 n + 1 ) π k 2 N ] X_k = \sum_{n=0}^{N-1}x_ncos[\frac{(2n+1)\pi k}{2N}] Xk=n=0∑N−1xncos[2N(2n+1)πk]
公式一
- x n x_n xn可以称作原始信号采样点,而cos可以被看做是这个样本点代表余弦频率的余弦波的采样点。如果把上面的公式一写成向量形式:
X = C x X=Cx 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_0=1 + 1 + 1 + 1 + 1 = 5
X0=1+1+1+1+1=5
因为所有的
x
n
x_n
xn都是1,所以公式就只剩下余弦部分了。
X
1
=
c
o
s
[
(
2
∗
0
+
1
)
∗
1
π
10
]
+
c
o
s
[
(
2
∗
1
+
1
)
∗
1
π
10
]
+
c
o
s
[
(
2
∗
2
+
1
)
∗
1
π
10
]
+
c
o
s
[
(
2
∗
3
+
1
)
∗
1
π
10
]
+
c
o
s
[
(
2
∗
4
+
1
)
∗
1
π
10
]
=
c
o
s
(
π
10
)
+
c
o
s
(
3
π
10
)
+
c
o
s
(
π
2
)
+
c
o
s
(
7
π
10
)
+
c
o
s
(
9
π
10
)
=
0
X_1 = cos[\frac{(2*0+1)*1\pi}{10}] + cos[\frac{(2*1+1)*1\pi}{10}] \\ + cos[\frac{(2*2+1)*1\pi}{10}] + cos[\frac{(2*3+1)*1\pi}{10}] + cos[\frac{(2*4+1)*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
X1=cos[10(2∗0+1)∗1π]+cos[10(2∗1+1)∗1π]+cos[10(2∗2+1)∗1π]+cos[10(2∗3+1)∗1π]+cos[10(2∗4+1)∗1π]=cos(10π)+cos(103π)+cos(2π)+cos(107π)+cos(109π)=0
其他三个X大家有兴趣也可以算一算,都是等于0。所以输出就是(5,0,0,0,0),当然opencv输出的shape是一个(5, 0)的数组,而且opencv在公式一的基础上还增加了一点东西,矩阵
C
C
C加了一个前缀:
C
i
j
(
N
)
=
α
j
/
N
c
o
s
(
π
(
2
k
+
1
)
j
2
N
)
α
0
=
1
,
α
j
=
2
,
j
>
0
C_{ij}^{(N)}=\sqrt{\alpha_j/N}cos(\frac{\pi (2k+1)j}{2N}) \\ \alpha_0=1, \alpha_j=2, j>0
Cij(N)=αj/Ncos(2Nπ(2k+1)j)α0=1,αj=2,j>0
所以输出值不是5,没有去细究了,大体的逻辑肯定是一样的。
line = (1, 1, 1, 1, 1)
dft = cv2.dft(np.float32(line))
print(dft.shape)
print(dft)
输出的dft如下:
- 物理含义,假设(1, 1, 1, 1, 1)是像素值,那么这张图片的像素之间是没有任何变化的,也可以理解没有频率,和直流电一样,没有变化,最终变换到频率域之后,就只有常量 X 0 X_0 X0上有作用,很好的表示了图像在频率上的描述。
- 而且常量值可以很好的表达像素值,如果像素值比较大,这个常量值 X 0 X_0 X0也会比较大。
- 在上面的向量公式中,矩阵C中的不同行向量之间是互相正交的,也就是点积为0。在这个小例子中,原始信号就是直流信号,只有 X 0 X_0 X0会对其影响,所以其他的行向量进行点乘后,结果都为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
C的逆矩阵
C
−
1
C^{-1}
C−1,再用这个矩阵去乘以权重向量X就可以得到原来的数组了。
x
=
C
−
1
X
x = C^{-1}X
x=C−1X
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)
就可以得到原来的数据了。
二维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)
打开图片一看,是这么个东西。
我的理解如下:回忆一下二维DCT的过程,在X方向上,从左往右的数据分别是
c
o
s
(
0
x
)
,
c
o
s
(
x
)
,
c
o
s
(
2
x
)
,
c
o
s
(
3
x
)
.
.
.
.
.
.
c
o
s
(
N
x
)
cos(0x), cos(x), cos(2x), cos(3x) ...... cos(Nx)
cos(0x),cos(x),cos(2x),cos(3x)......cos(Nx)
的权重系数
X
0
,
X
1
,
X
2
.
.
.
.
.
.
X
N
X_0,X_1,X_2......X_N
X0,X1,X2......XN
而在图像中,只有左上角是比较黑的,所以可以说这幅图像在X方向上的频率(也就是像素值的变化率)集中在低频区。换言之,在Y轴方向上的频率也是集中在低频区。这样就可以采用一些低通,高通等不同的滤波算法进行区分和操作了。
当然,我个人觉得这里的频率变换可能也是相对的,我感觉最右边的可能也到了 c o s ( 100 x ) cos(100x) cos(100x)的高频了吧。
DFT,离散傅里叶变换
在图像处理中,用的更多的是这个DFT,也就是离散傅里叶变换。应该说DCT是DFT的一种特例,DCT说完了,就再说说DFT。
从DFT的公式说起:
X
k
=
∑
n
=
0
N
−
1
x
n
e
−
2
π
j
k
n
N
X_k = \sum_{n=0}^{N-1}x_ne^{-2\pi j\frac{kn}{N}}
Xk=n=0∑N−1xne−2πjNkn
然后根据欧拉公式:
e
j
w
=
c
o
s
(
w
)
−
j
s
i
n
(
w
)
e^{jw}=cos(w)-jsin(w)
ejw=cos(w)−jsin(w)
所以DFT公式可以展开为:
X
k
=
∑
n
=
0
N
−
1
x
n
[
c
o
s
(
2
π
k
n
/
N
)
−
j
s
i
n
(
2
π
k
n
/
N
)
]
X_k = \sum_{n=0}^{N-1}x_n[cos(2\pi kn/N) - jsin(2\pi kn/N)]
Xk=n=0∑N−1xn[cos(2πkn/N)−jsin(2πkn/N)]
相对于DCT,DFT有两个部分,实数部分与虚数部分,其中实数部分是和DCT一样的。那么虚数部分代表啥呢,图像中的像素是一个一个真实的像素点,而且是0-255之间的一个整数。
其实我的理解是在opencv中,傅里叶变换就是把这两个部分做了图像两个维度的变换,一个叫频率谱(实数部分),一个叫相位谱(虚数部分)。不要把它和数学中的虚数一一对应,就是两个维度的数字就好。
其中频率谱部分在DCT的部分已经说过了,就是代表了图像两个方向的像素变换情况。
那么相位谱又代表了什么呢?我在网上找了蛮久的资料,没有找到准确合适的描述,下面这段话是我自己的理解,去理解图像的相位谱到底代表啥意思。
- 首先,图像的数据是离散的,在信号处理领域,如果信号上的点是离散的,相当于对信号进行采样:
图中有8个采样点,每个点的相位就是他在X轴上的位置,比如第一个点的相位就是 π / 8 \pi /8 π/8。
- 那么回到图像中,一行图像的数据可以被认为是某种信号,这个上面说过了。
那么这个信号可以拆成若干个频率不同的余弦波。那么我的理解是
傅里叶的虚数部分,也就是相位谱,就是这个采样点在频率为x或者2x的不同频率中的波信号中所对应的相位。
比如第一个点,可能在 cos(x)这个信号中的相位是 π / 9 \pi /9 π/9,cos(2x)这个频率的信号中的相位是 π / 4 \pi /4 π/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()
输出图像为: