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

图像的傅里叶变换

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

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

f

(

x

,

y

)

→

F

(

u

,

v

)

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

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

F

(

u

,

v

)

=

∑

x

=

0

M

−

1

∑

y

=

0

N

−

1

f

(

x

,

y

)

e

−

j

2

π

(

u

x

M

+

v

y

N

)

F(u, v) = sum_{x=0}^{M-1} sum_{y=0}^{N-1} f(x,y)e^{-j2pi(frac{ux}{M} + frac{vy}{N})}

F(u,v)=x=0∑M−1​y=0∑N−1​f(x,y)e−j2π(Mux​+Nvy​)

• $f(x,y)$
• $M,N$
• $F(x,y)$

F

(

u

,

v

)

→

f

(

x

,

y

)

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

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

f

(

x

,

y

)

=

1

M

N

∑

u

=

0

M

−

1

∑

v

=

0

N

−

1

F

(

u

,

v

)

e

j

2

π

(

u

x

M

+

v

y

N

)

f(x, y) = frac{1}{MN} sum_{u=0}^{M-1} sum_{v=0}^{N-1} F(u, v)e^{j2pi(frac{ux}{M} + frac{vy}{N})}

f(x,y)=MN1​u=0∑M−1​v=0∑N−1​F(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)

​
相位角：

ϕ

=

a

r

c

t

a

n

I

(

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

)

2

P(u, v) = |F(u, v)|^2 = {R(u, v)}^2 + {I(u, v)}^2

P(u,v)=∣F(u,v)∣2=R(u,v)2+I(u,v)2

    傅里叶变换的性质：可分离性和对称性。在进行2D傅里叶变换时，可以利用其性质来简化计算，一个2D傅里叶变换核可以分解为两个1D傅里叶变换。一个

N

×

N

Ntimes N

N×N的2D傅里叶变换需要

N

4

N^4

N4 次复数乘法运算和

N

2

(

N

2

−

1

)

N^2(N^2-1)

N2(N2−1) 次复数加法运算，而进行一个长度为N的1D傅里叶变换只需进行

N

2

N^2

N2次复数乘法运算和

N

(

N

−

1

)

N(N-1)

N(N−1)次复数加法运算。

指数运算部分：

e

−

j

2

π

(

u

x

M

+

v

y

N

)

=

e

−

j

2

π

u

x

M

+

e

−

j

2

π

v

y

N

e^{-j2pi(frac{ux}{M} + frac{vy}{N})} = e^{-j2pifrac{ux}{M}} + e^{-j2pifrac{vy}{N}}

e−j2π(Mux​+Nvy​)=e−j2πMux​+e−j2πNvy​

e

j

2

π

(

u

x

M

+

v

y

N

)

=

e

j

2

π

u

x

M

+

e

j

2

π

v

y

N

e^{j2pi(frac{ux}{M} + frac{vy}{N})} = e^{j2pifrac{ux}{M}} + e^{j2pifrac{vy}{N}}

ej2π(Mux​+Nvy​)=ej2πMux​+ej2πNvy​

列变换：

G

(

x

,

y

)

=

∑

y

=

0

N

−

1

f

(

x

,

y

)

e

−

j

2

π

v

y

N

G(x, y) = sum_{y = 0}^{N-1} f(x, y) e^{frac{-j2pi vy}{N}}

G(x,y)=y=0∑N−1​f(x,y)eN−j2πvy​
行变换：

F

(

u

,

v

)

=

∑

x

=

0

M

−

1

G

(

x

,

y

)

e

−

j

2

π

u

x

M

F(u, v) = sum_{x = 0}^{M-1} G(x, y) e^{frac{-j2pi ux}{M}}

F(u,v)=x=0∑M−1​G(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.shape
    for 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 ret
    
def slowly_idft(img_freq):
    img_gray = np.zeros_like(img_freq)
    m, n = img_freq.shape
    for 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*n
    img_gray = np.clip(img_gray, 0, 255).astype(int)
    return img_gray

def restore_img(idft_array):
    real_part = idft_array.real
    img_gray = np.clip(real_part, 0, 255).astype("uint8")
    return img_gray

def accelerate_dft(img_gray):
    ret = np.zeros_like(img_gray).astype('complex')
    m, n = img_gray.shape
    x_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)$
• $H(u,v)$
• $G(u,v)$

图像频域增强步骤：

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

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

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\times 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])   # 还原图像

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