图像频域增强：傅里叶变换

图像的傅里叶变换

图像的频域增强其实就是：先对图像进行变换，将图像转换到变换域，然后在变换域进行操作以实现图像增强。常用的变换域就是频域（频率域，傅里叶变换的结果）。图像的频域增强有直观的物理意义。例如，图像模糊是图像中高频分量不足的结果，在频域里增加高频分量或减少低频分量就能消除一些模糊。又如，图像有时会受到重复出现的有规律周期噪声的影响，周期噪声具有特定的频率，所以可以采取频域滤波的方法滤除相应噪声频率，从而消除周期噪声。

傅里叶变换(2D)及其性质

f(x,y)→F(u,v)f(x, y) \to F(u, v)f(x,y)→F(u,v)

F(u,v)=∑x=0M−1∑y=0N−1f(x,y)e−j2π(uxM+vyN)F(u, v) = \sum_{x=0}^{M-1} \sum_{y=0}^{N-1} f(x,y)e^{-j2\pi(\frac{ux}{M} + \frac{vy}{N})}F(u,v)=x=0∑M−1y=0∑N−1f(x,y)e−j2π(Mux+Nvy)

f(x,y)f(x, y)f(x,y) : 像素矩阵中x,yx, yx,y处的像素值
M,NM, NM,N : 像素矩阵的大小，行列
F(x,y)F(x, y)F(x,y) : 变换后u,vu, vu,v处的频域值，为复数。

F(u,v)→f(x,y)F(u, v) \to f(x, y)F(u,v)→f(x,y)

f(x,y)=1MN∑u=0M−1∑v=0N−1F(u,v)ej2π(uxM+vyN)f(x, y) = \frac{1}{MN} \sum_{u=0}^{M-1} \sum_{v=0}^{N-1} F(u, v)e^{j2\pi(\frac{ux}{M} + \frac{vy}{N})}f(x,y)=MN1u=0∑M−1v=0∑N−1F(u,v)ej2π(Mux+Nvy)

F(u,v)=R(u,v)+I(u,v)F(u, v) = R(u, v) + I(u, v)F(u,v)=R(u,v)+I(u,v)
频谱：
∣F(u,v)∣=(R(u,v)2+I(u,v)2)|F(u, v)| = \sqrt{({R(u, v)}^2 + {I(u, v)}^2)}∣F(u,v)∣=(R(u,v)2+I(u,v)2)
相位角：
ϕ=arctanI(u,v)R(u,v)\phi = arctan \frac{I(u, v)}{R(u, v)}ϕ=arctanR(u,v)I(u,v)
功率谱：
P(u,v)=∣F(u,v)∣2=R(u,v)2+I(u,v)2P(u, v) = |F(u, v)|^2 = {R(u, v)}^2 + {I(u, v)}^2P(u,v)=∣F(u,v)∣2=R(u,v)2+I(u,v)2

傅里叶变换的性质：可分离性和对称性。在进行2D傅里叶变换时，可以利用其性质来简化计算，一个2D傅里叶变换核可以分解为两个1D傅里叶变换。一个N×NN\times NN×N的2D傅里叶变换需要 N4N^4N4 次复数乘法运算和 N2(N2−1)N^2(N^2-1)N2(N2−1) 次复数加法运算，而进行一个长度为N的1D傅里叶变换只需进行N2N^2N2次复数乘法运算和N(N−1)N(N-1)N(N−1)次复数加法运算。

指数运算部分：
e−j2π(uxM+vyN)=e−j2πuxM+e−j2πvyNe^{-j2\pi(\frac{ux}{M} + \frac{vy}{N})} = e^{-j2\pi\frac{ux}{M}} + e^{-j2\pi\frac{vy}{N}}e−j2π(Mux+Nvy)=e−j2πMux+e−j2πNvy

ej2π(uxM+vyN)=ej2πuxM+ej2πvyNe^{j2\pi(\frac{ux}{M} + \frac{vy}{N})} = e^{j2\pi\frac{ux}{M}} + e^{j2\pi\frac{vy}{N}}ej2π(Mux+Nvy)=ej2πMux+ej2πNvy

列变换：
G(x,y)=∑y=0N−1f(x,y)e−j2πvyNG(x, y) = \sum_{y = 0}^{N-1} f(x, y) e^{\frac{-j2\pi vy}{N}}G(x,y)=y=0∑N−1f(x,y)eN−j2πvy
行变换：
F(u,v)=∑x=0M−1G(x,y)e−j2πuxMF(u, v) = \sum_{x = 0}^{M-1} G(x, y) e^{\frac{-j2\pi ux}{M}}F(u,v)=x=0∑M−1G(x,y)eM−j2πux

实现傅里叶变换

简单实现慢速傅里叶变换
加速版本
numpy.fft.fft2()

import cv2 as cv
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt

img = cv.imread('./images/guigu1_mini.png')
img = cv.cvtColor(img, cv.COLOR_BGR2RGB)
img_gray = cv.cvtColor(img, cv.COLOR_RGB2GRAY)

测试用例:

傅里叶变换与逆变换

def slowly_dft(img_gray):ret = np.zeros_like(img_gray).astype('complex')# ret[0]应包含零频率项，# ret[1:n//2]应该包含正频率项，# ret[n//2 + 1:]应该包含负频率项，从最负频率开始按升序排列。m, n = img_gray.shapefor u in range(m):for v in range(n):for x in range(m):for y in range(n):ret[u, v] += img_gray[x, y]*np.exp(-2j*np.pi*(u*x/m + v*y/n))return retdef slowly_idft(img_freq):img_gray = np.zeros_like(img_freq)m, n = img_freq.shapefor x in range(m):for y in range(n):for u in range(m):for v in range(n):img_gray[x, y] += img_freq[u, v]*np.exp(2j*np.pi*(u*x/m + v*y/n))/m*nimg_gray = np.clip(img_gray, 0, 255).astype(int)return img_graydef restore_img(idft_array):real_part = idft_array.realimg_gray = np.clip(real_part, 0, 255).astype("uint8")return img_graydef accelerate_dft(img_gray):ret = np.zeros_like(img_gray).astype('complex')m, n = img_gray.shapex_arr = np.concatenate([np.arange(m).reshape(m, 1)], axis=1)y_arr = np.concatenate([np.arange(n).reshape(1, n)], axis=0)for u in range(m):for v in range(n):ret[u, v] = np.sum(np.multiply(img_gray,np.exp(-2j*np.pi*(u*x_arr/m + v*y_arr/n))))return ret

%%time
img_freq_slow = slowly_dft(img_gray)

CPU times: total: 1min 3s
Wall time: 1min 3s

%%time
img_freq_acce = accelerate_dft(img_gray)

CPU times: total: 1.08 s
Wall time: 1.08 s

%%time
img_freq_fast = np.fft.fft2(img_gray)

CPU times: total: 0 ns
Wall time: 0 ns

快速傅里叶逆变换

img_1 = np.fft.ifft2(img_freq_slow)
img_2 = np.fft.ifft2(img_freq_fast)

还原和频谱图（log变换,0频率移动）

图像频域增强(OpenCV)

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

F(u,v)F(u, v)F(u,v) : 傅里叶变换后的图像频域矩阵
H(u,v)H(u, v)H(u,v) : 转移函数，作用于f(u,v)f(u, v)f(u,v)
G(u,v)G(u, v)G(u,v) : 增强后的频域矩阵

图像频域增强步骤：

傅里叶变换得到频域矩阵。
将其与一个转移函数（根据需要设计）相乘。
将结果进行傅里叶反变换以得到增强后的图像。

下通过一个简单转移操作，实现高通滤波，观察原始图像的变化：

img = cv.imread('./images/hj.png')
img = cv.cvtColor(img, cv.COLOR_BGR2RGB)
img_gray = cv.cvtColor(img, cv.COLOR_RGB2GRAY)

dft = cv.dft(img_gray.astype('float32'),flags = cv.DFT_COMPLEX_OUTPUT)           # 傅里叶变换(Opencv是用深度为2数组表示复数)
dft_shift = np.fft.fftshift(dft)                                                 # 移动零频分量
magnitude_spectrum = 20*np.log(cv.magnitude(dft_shift[:,:,0],dft_shift[:,:,1]))  # 幅值对数变换

测试用例（逃学威龙1）

从图像的幅值可以看出，中间最亮也代表着中间频率最低(幅值越大频率越小)。

1. 屏蔽中间10×1010 \times 1010×10的低频

crow = dft_shift.shape[0]//2                         # 中心点坐标
ccol = dft_shift.shape[1]//2
dft_shift[crow-5:crow+5, ccol-5:ccol+5] = 0          # 去除中心低频
dft_ishift = np.fft.ifftshift(dft_shift)             # 还原频谱图
img_ = cv.idft(dft_ishift)                           # 逆傅里叶变换
img_back1 = cv.magnitude (img_[:,:,0],img_[:,:,1])   # 还原图像

屏蔽部分低频后(高通滤波)，图像得到了类似边缘检测的结果。