数字图像的频率域滤波

引言

“任何”周期信号都可以用一系列成谐波关系的正弦曲线来表示。
—— 让·巴普蒂斯·约瑟夫·傅里叶

由于频率和与时间变化率直接相关，因此容易想到 傅里叶变换中的频率与图像中的灰度值变化模式 也可以关联起来。

空间域　英文： spatial domain。释义：又称图像空间(image space)。由图像像元组成的空间。在图像空间中以长度(距离)为自变量直接对像元值进行处理称为空间域处理。

频率域 英文： spatial frequency domain。释义：以频率(即波数)为自变量描述图像的特征,可以将一幅图像像元值在空间上的变化分解为具有不同振幅、空间频率和相位的简振函数的线性叠加,图像中各种频率成分的组成和分布称为空间频谱。这种对图像的频率特征进行分解、处理和分析称为频率域处理或波数域处理。

二者关系：

空间域与频率域可互相转换。在频率域中可以引用已经很成熟的频率域技术,处理的一般步骤为：
①对图像施行二维离散傅立叶变换或小波变换,将图像由图像空间转换到频域空间。
②在频率域中对图像的频谱作分析处理,以改变图像的频率特征。即设计不同的数字滤波器,对图像的频谱进行滤波。
③通过反变换将图像从频率域变换到空间域。

频率域处理主要用于与图像空间频率有关的处理中。如图像恢复、图像重建、辐射变换、边缘增强、图像锐化、图像平滑、噪声压制、频谱分析、纹理分析等处理和分析中。

1.图像傅里叶变换的物理意义

频谱图(点击查看代码)

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

# plt 显示中文
import matplotlib as mpl
mpl.rcParams['font.family'] = 'SimHei'
plt.rcParams['axes.unicode_minus'] = False


img = cv2.imread("640_400.jpg", 0)

# 傅里叶变换
f = np.fft.fft2(img)
fshift = np.fft.fftshift(f)
f_magnitude_spectrum = np.log(np.abs(f))
fshift_magnitude_spectrum = np.log(np.abs(fshift))

# opencv 显示
# magnitude_spectrum = magnitude_spectrum.astype("int8")
# cv2.imshow("test", magnitude_spectrum)
# cv2.waitKey(0)

# plt 显示
plt.subplot(131)
plt.imshow(img, cmap='gray')
plt.title("原图"), plt.xticks([]), plt.yticks([])

plt.subplot(132)
plt.imshow(f_magnitude_spectrum, cmap='gray')
plt.title("二维FFT"), plt.xticks([]), plt.yticks([])

plt.subplot(133)
plt.imshow(fshift_magnitude_spectrum, cmap='gray')
plt.title("零频率分量移至频谱中心"), plt.xticks([]), plt.yticks([])

plt.show()

# 1. np.fft.fft2() # Numpy 提供的实现傅里叶变换的函数
#  - 这里需要注意的是，参数“原始图像”的类型是灰度图像，函数的返回值是一个复数数组（complex ndarray）。
#  - 经过该函数的处理，就能得到图像的频谱信息，此时，图像频谱中的零频率分量位于频谱图像（频域图像）的左上角
# 2. np.fft.fftshift() # 将零频率成分移动到频域图像的中心位置
# 3. np.log(np.abs(频谱值)) # 对图像进行傅里叶变换后，得到的是一个复数数组。为了显示为图像，需要将它们的值调整到 [0,255] 的灰度空间内

说明：

图像的频率是表征图像中灰度变化剧烈程度的指标，是灰度在平面空间上的梯度。
如：大面积的沙漠在图像中是一片灰度变化缓慢的区域，对应的频率值很低；而对于地表属性变换剧烈的边缘区域在图像中是一片灰度变化剧烈的区域，对应的频率值较高。
傅里叶变换在实际中有非常明显的物理意义，从物理效果看，傅里叶变换是将图像从空间域转换到频率域，其逆变换是将图像从频率域转换到空间域。换句话说，傅里叶变换的物理意义是将图像的灰度分布函数变换为图像的频率分布函数。
为了便于频域和频谱分析，常在变换后进行频谱中心化，即对调频谱的四个象限。
频谱中心化后，中间最亮的点是低频率，属于直流分量，越往外频率越高。

低频对应图像中变化缓慢的灰度分量
高频对应图像中灰度急剧变化的物体边缘或其他分量

低通滤波器(点击查看代码)

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

# plt 显示中文
import matplotlib as mpl
mpl.rcParams['font.family'] = 'SimHei'
plt.rcParams['axes.unicode_minus'] = False


img = cv2.imread("640_400.jpg", 0)

# 傅里叶变换
f = np.fft.fft2(img)
fshift = np.fft.fftshift(f)

# 获取原图长宽和中心点
(h, w) = img.shape
ch, cw = h//2, w//2

# 低通滤波提取高频分量
mask = np.zeros((h, w), np.uint8)
mask[ch-30:ch+30, cw-30:cw+30] = 1 # 低通滤波器

fshift = fshift * mask
f_ishift = np.fft.ifftshift(fshift)
img_back = np.fft.ifft2(f_ishift)
img_back = np.abs(img_back)

plt.subplot(131)
plt.imshow(img, cmap= "gray")
plt.title("原图"), plt.xticks([]), plt.yticks([])

plt.subplot(132)
plt.imshow(mask * 255, cmap= "gray")
plt.title("低通滤波器"), plt.xticks([]), plt.yticks([])

plt.subplot(133)
plt.imshow(img_back, cmap="gray")
plt.title("低通滤波后（低频分量）"), plt.xticks([]), plt.yticks([])

plt.show()

高通滤波器（点击查看代码）

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

# plt 显示中文
import matplotlib as mpl
mpl.rcParams['font.family'] = 'SimHei'
plt.rcParams['axes.unicode_minus'] = False


img = cv2.imread("640_400.jpg", 0)

# 傅里叶变换
f = np.fft.fft2(img)
fshift = np.fft.fftshift(f)

# 获取原图长宽和中心点
(h, w) = img.shape
ch, cw = h//2, w//2

# 低通滤波提取高频分量
mask = np.zeros((h, w), np.uint8)
mask[ch-30:ch+30, cw-30:cw+30] = 1 # 低通滤波器

fshift = fshift * mask
f_ishift = np.fft.ifftshift(fshift)
img_back = np.fft.ifft2(f_ishift)
img_back = np.abs(img_back)

plt.subplot(131)
plt.imshow(img, cmap= "gray")
plt.title("原图"), plt.xticks([]), plt.yticks([])

plt.subplot(132)
plt.imshow(mask * 255, cmap= "gray")
plt.title("低通滤波器"), plt.xticks([]), plt.yticks([])

plt.subplot(133)
plt.imshow(img_back, cmap="gray")
plt.title("低通滤波后（低频分量）"), plt.xticks([]), plt.yticks([])

plt.show()

2. 频率滤波步骤

频率滤波步骤

输入图像f(x,y)
做二维FFT得到 F(u,v)
频谱中心化
构建一个实对称滤波器传递函数 H(u,v)
采用对应像素相乘得到G(u,v)=H(u,v)F(u,v)
对 G(u,v) 进行 IDFT 得到滤波后的图像

3. 常见滤波器

3.1 低通滤波器

常见低通滤波器（点击查看代码）

import matplotlib.pyplot as plt
import numpy as np


def cal_distance(pa, pb):
    from math import sqrt
    dis = sqrt((pa[0] - pb[0]) ** 2 + (pa[1] - pb[1]) ** 2)
    return dis


# 滤波器模板
class Filter:

    @classmethod
    def generate_filter(cls, d, shape, *args, **kwargs):
        # 这里要反过来 shape 的两个维度
        transfer_matrix = np.zeros((shape[0], shape[1]))
        center_point = tuple(map(lambda x: (x - 1) // 2, shape))
        for i in range(transfer_matrix.shape[0]):
            for j in range(transfer_matrix.shape[1]):
                dist = cal_distance(center_point, (i, j))
                transfer_matrix[i, j] = cls.get_one(d, dist, *args, **kwargs)
        return transfer_matrix

    @classmethod
    def get_one(cls, d, dist, *args, **kwargs) -> float:
        return 1

# 理想低通滤波器
class ILPFFilter(Filter):

    @classmethod
    def get_one(cls, d, dist, *args, **kwargs) -> float:
        if dist <= d:
            return 1
        else:
            return 0

# 高斯低通滤波器
class GLPFFilter(Filter):

    @classmethod
    def get_one(cls, d, dist, *args, **kwargs) -> float:
        return np.exp(-(dist ** 2) / (2 * (d ** 2)))


# 巴特沃斯低通滤波器
class BLPFFilter(Filter):

    @classmethod
    def get_one(cls, d, dist, *args, **kwargs) -> float:
        n = kwargs["n"]
        return 1 / ((1 + dist / d) ** (2 * n))

if __name__ == "__main__":
    # plt 显示中文
    import matplotlib as mpl
    mpl.rcParams['font.family'] = 'SimHei'
    plt.rcParams['axes.unicode_minus'] = False

    ILPFFilter = ILPFFilter.generate_filter(10, (100, 100), w=5)
    GLPFFilter = GLPFFilter.generate_filter(10, (100, 100), w=5)
    BLPFFilter = BLPFFilter.generate_filter(10, (100, 100), w=5, n=1)

    plt.subplot(131)
    plt.imshow(ILPFFilter, cmap="gray"), plt.xticks([]), plt.yticks([])
    plt.title("理想低通滤波器")
    plt.subplot(132)
    plt.imshow(GLPFFilter, cmap="gray"), plt.xticks([]), plt.yticks([])
    plt.title("高斯低通滤波器")
    plt.subplot(133)
    plt.imshow(BLPFFilter, cmap="gray"), plt.xticks([]), plt.yticks([])
    plt.title("巴特沃斯低通滤波器 n=1")

    plt.show()

理想低通滤波器

传递函数为:

\[H(u,v)= \left\{\begin{matrix} 1, D(u,v)\leqslant D_0 \\ 0, D(u,v)> D_0 \end{matrix}\right. \]

式中，\(D_0\)是一个正常数，\(D(u,v)\)是频率域中点\((u,v)\)到\(P×Q\)频率矩形中心的距离，即

\[D(u,v) = [(u-\frac{P}{2})^{2}+(v-\frac{Q}{2})^{2}]^{\frac{1}{2}} \]

高斯低通滤波器

传递函数为:

\[H(u,v)=e^{-D^{2}(u,v)/2D_{0}^{2}} \]

巴特沃斯低通滤波器

传递函数为:

\[H(u,v)=\frac{1}{1+[D(u,v)/D_0)]^{2n}]} \]

巴特沃斯滤波器的形状由滤波器阶数 \(n\)控制，当阶数取值较大时，巴特沃斯滤波器接近于理想滤波器，取值较小时，接近于高斯滤波器。

3.2 高通滤波器

常见高通滤波器（点击查看代码）

# 理想高通滤波器
class IHPFFilter(Filter):

    @classmethod
    def get_one(cls, d, dist, *args, **kwargs) -> float:
        if dist <= d:
            return 0
        else:
            return 1

# 高斯高通滤波器
class GHPFFilter(Filter):

    @classmethod
    def get_one(cls, d, dist, *args, **kwargs) -> float:
        return 1 - np.exp(-(dist ** 2) / (2 * (d ** 2)))

# 巴特沃斯高通滤波器
class BHPFFilter(Filter):

    @classmethod
    def get_one(cls, d, dist, *args, **kwargs) -> float:
        n = kwargs["n"]
        return 1 - 1 / ((1 + dist / d) ** (2 * n))

# 拉普拉斯
class LPFilter(Filter):

    @classmethod
    def get_one(cls, d, dist, *args, **kwargs) -> float:
        return 4 * math.pi**2 * dist**2

在频率域中用1减去低通滤波器传递函数，会得到对应的高通滤波器传递函数。

3.3 选择性滤波

处理特定频带的滤波器称为频带滤波器。

若某个频带中的频率被滤除，则称该频带滤波器为带阻滤波器；
若某个频带中的频率被通过，则称该频带滤波器为带通滤波器；
可处理小频率区域的滤波器称为陷波滤波器，根据陷波区域中的频率被滤除还是被通过，可将陷波滤波器进一步分为带阻陷波滤波器和带通陷波滤波器。

常见选择性滤波（点击查看代码）

# 巴特沃斯带通滤波
class BBPFFilter(Filter):

    @classmethod
    def get_one(cls, d, dist, *args, **kwargs) -> float:
        w = kwargs["w"]
        n = kwargs["n"]
        if dist == d:
            return 1
        else:
            return 1- 1 / (1 + ((dist * w) / (dist ** 2 - d ** 2)) ** (2 * n))

# 巴特沃斯带阻滤波
class BBRFFilter(Filter):

    @classmethod
    def get_one(cls, d, dist, *args, **kwargs) -> float:
        w = kwargs["w"]
        n = kwargs["n"]
        if dist == d:
            return 0
        else:
            return 1 / (1 + ((dist * w) / (dist ** 2 - d ** 2)) ** (2 * n))

# 高斯带通滤波
class GBPFFilter(Filter):

    @classmethod
    def get_one(cls, d, dist, *args, **kwargs) -> float:
        w = kwargs["w"]
        return np.exp(-((dist ** 2 - d ** 2) / (d * w)) ** 2)

# 高斯带阻滤波
class GBRFFilter(Filter):

    @classmethod
    def get_one(cls, d, dist, *args, **kwargs) -> float:
        w = kwargs["w"]
        return 1-np.exp(-((dist ** 2 - d ** 2) / (d * w)) ** 2)


# 理想带通滤波器
class IBPFFilter(Filter):

    @classmethod
    def get_one(cls, d, dist, *args, **kwargs) -> float:
        w = kwargs["w"]
        if d - w / 2 <= dist <= d + w / 2:
            return 1
        else:
            return 0

# 理想带阻滤波器
class IBRFFilter(Filter):

    @classmethod
    def get_one(cls, d, dist, *args, **kwargs) -> float:
        w = kwargs["w"]
        if d - w / 2 <= dist <= d + w / 2:
            return 0
        else:
            return 1

空间域和频率域滤波之间的对应关系

空间卷积运算非常适合采用硬件和/或固件的小型核。但在使用普通计算机时，使用快速傅里叶变换（FFT）算法计算DFT的频率域方法，要比空间卷积运算的速度快几百倍，具体取决于所用核的大小。

频率域中的滤波概念更加直观，且频率域中的滤波器设计也更容易。充分利用频率域和空间域性质的一种方法是，在频率域中规定一个滤波器，计算其IDFT，然后利用生成的全尺寸空间核的性质，指导构建较小的核。

\[h(x,y) = F ^{-1}[H(u,v)] \]

验证1：使用低通频率域滤波器平滑图像

理想低通滤波器
高斯低通滤波器
巴特沃斯低通滤波器

验证2：使用高通滤波器锐化图像

理想高通滤波器
高斯高通滤波器
巴特沃斯高通滤波器
拉普拉斯
同态滤波

总结

使用低通频率域滤波器平滑图像，常见的低通滤波器有理想低通滤波器、高斯低通滤波器、巴特沃斯低通滤波器。
使用高通滤波器锐化图像，常见的高通滤波器有理想高通滤波器、高斯高通滤波器、巴特沃斯高通滤波。
理想滤波器非常尖锐；高斯滤波器非常平滑；巴特沃斯滤波器的形状由滤波器阶数控制，当阶数取值较大时，巴特沃斯滤波器接近于理想滤波器，取值较小时，接近于高斯滤波器。
理想滤波器因为具有振铃效应，并不实用；频率域高斯函数的傅里叶反变化也是高斯函数，说明通过反傅里叶变换后得到的空间高斯滤波器核没有振铃效应；针对巴特沃斯滤波器，滤波器阶数为1时，没有振铃效应，阶数越大，振铃效应越明显。
在频率域中用1减去低通滤波器传递函数，会得到对应的高通滤波器传递函数。
陷波滤波器是最有用的选择性滤波器，可去除周期干扰、印刷物图像中的莫尔模式。
空间域滤波核频率域滤波之间的纽带是卷积定理。
频率域卷积比空间域卷积快几百倍（具体取决于所用核的大小）。

参考：

《数字图像处理》第四版
数字图像处理与Python实现笔记之频域滤波
File:Fourier transform time and frequency domains.gif
numpy.fft.fft2
numpy.fft.fftshift
《数字图像处理(第四版)》阅读随笔 ch4 (频率域滤波)
https://github.com/YuiCtwo/Learn-DIP

数字图像处理

数字图像的频率域滤波

引言

1.图像傅里叶变换的物理意义

2. 频率滤波步骤

3. 常见滤波器

3.1 低通滤波器

理想低通滤波器

高斯低通滤波器

巴特沃斯低通滤波器

3.2 高通滤波器

3.3 选择性滤波

空间域和频率域滤波之间的对应关系

验证1：使用低通频率域滤波器平滑图像

验证2：使用高通滤波器锐化图像

总结

参考：

相关

matlab数字图像处理(3)——数字形态学处理

数字图像处理整理（2022）

标签