youcans 的 OpenCV 学习课—9.频率域图像滤波（下）

Posted 2022-02-24 小白YouCans

youcans 的 OpenCV 学习课—8.频率域图像滤波（下）

本系列面向 Python 小白，从零开始实战解说 OpenCV 项目实战。
图像滤波是在尽可能保留图像细节特征的条件下对目标图像的噪声进行抑制，是常用的图像预处理操作。
频率域图像处理先将图像进行傅里叶变换，然后在变换域进行处理，最后进行傅里叶反变换转换回空间域。
频率域滤波是滤波器传递函数与输入图像的傅里叶变换的对应像素相乘。频率域中的滤波概念更加直观，滤波器设计也更容易。使用快速傅里叶变换算法，比空间卷积运算速度快很多。
本文提供上述各种算法的完整例程和运行结果，并通过一个应用案例示范综合使用多种图像增强方法。

欢迎关注 『youcans 的 OpenCV 学习课』 系列，持续更新

文章目录

• youcans 的 OpenCV 学习课—8.频率域图像滤波（下）
• • 3. 频率域低通滤波器
  • • 3.1 频率域图像滤波基础
    • 例程 8.13：简单的频率域图像滤波
    • 3.2 频率域图像滤波的步骤
    • 例程 8.15：频率域图像滤波的一般步骤
    • 例程 8.16：频率域图像滤波的一般步骤（理想低通滤波器）
    • 3.3 频率域高斯低通滤波器（GLPF）
    • 例程 8.18：频率域高斯低通滤波器 (GLPF)
    • 3.4 频率域巴特沃斯低通滤波器（BLPF）
    • 例程 8.19：频率域巴特沃斯低通滤波器 (BLPF)
    • 3.5 频率域低通滤波案例：印刷文本字符修复
  • 4. 频率域高通滤波器
  • • 4.1 由低通滤波器得到高通滤波器
    • 例程 8.23 由低通滤波器得到高通滤波器
    • 例程 8.24 频率域高通滤波器的应用
    • 例程 8.25：指纹图像处理（高通滤波+阈值处理）
    • 4.2 频率域钝化掩蔽
    • 例程 8.26：频率域钝化掩蔽
    • 4.3 频率域拉普拉斯高通滤波
    • 例程 8.27：频率域图像锐化（拉普拉斯算子）
  • 5. 选择性滤波
  • • 5.1 带阻与带通
    • 例程 8.28 带阻滤波器的传递函数
    • 5.2 陷波滤波器 （Notch Filter）
    • 例程 8.29 使用陷波滤波删除数字化印刷图像中的莫尔模式

---

3. 频率域低通滤波器

图像变换是对图像信息进行变换，使能量保持但重新分配，以便于滤除噪声、加强感兴趣的部分或特征。

3.1 频率域图像滤波基础

傅里叶变换的目的是将图像从空间域转换到频率域，在频率域内实现对图像中特定信息的处理，在图像分析、图像增强、图像去噪、边缘检测、特征提取、图像压缩和加密中都有重要的应用。

空间取样和频率间隔是相互对应的，频域中样本之间的间隔，与空间样本之间的间隔及样本数量的乘积成反比。

空间域滤波器和频率域滤波器也是相互对应的，二者形成一个傅里叶变换对：
$$\begin{gathered} f(x,y)\otimes h(x,y)\iff F(u,v)H(u,v) \\ f(x,y)h(x,y)\iff F(u,v)\otimes H(u,v) \end{gathered}$$
也就是说，空间域滤波器和频率域滤波器实际上是相互对应的，有些空间域滤波器在频率域通过傅里叶变换实现会更方便、更快速。

对信号或图像进行傅里叶变换后，可以得到信号或图像的低频信息和高频信息。低频信息对应图像中缓慢变化的灰度分量，高频信息则对应尖锐变化的灰度分量。

低通滤波就是保留傅里叶变换的低频信息、削弱高频信息，而高通滤波则是保留傅里叶变换的高频信息、削弱低频信息。

低频滤波器本质上就是构造一个矩阵，越靠近中心的位置越接近于 1，而远离中心位置的值则接近于 0。简单地，生成一个矩形窗口遮罩，在黑色（置 0）遮罩图像的中心开有白色（置 1）窗口，就得到一个低通滤波器。

例程 8.13：简单的频率域图像滤波

    # 8.13：简单的频率域图像滤波(窗口遮罩低通滤波器)
    imgGray = cv2.imread("../images/imgLena.tif", flags=0)  # flags=0 读取为灰度图像
    height, width = imgGray.shape[:2]  # 图片的高度和宽度
    centerY, centerX = int(height/2), int(width/2)  # 图片中心

    # （1）首先对图像进行傅里叶变换
    imgFloat32 = np.float32(imgGray)  # 将图像转换成 float32
    dft = cv2.dft(imgFloat32, flags=cv2.DFT_COMPLEX_OUTPUT)  # 傅里叶变换
    dftShift = np.fft.fftshift(dft)  # 将低频分量移动到频域图像的中心

    d = [20, 40, 80]
    plt.figure(figsize=(9, 6))
    for i in range(3):
        # 构造低通滤波器矩形窗口遮罩 mask
        mask = np.zeros((height, width, 2), np.uint8)
        mask[centerY-d[i]:centerY+d[i], centerX-d[i]:centerX+d[i]] = 1  # 设置低通滤波矩形窗口遮罩，过滤高频
        maskAmp = np.uint8(np.sqrt(np.power(mask[:,:,0], 2) + np.power(mask[:,:,1], 2)))
        print("d=, maskAmp: max=, min=".format(d[i],maskAmp.max(), maskAmp.min()))

        # （2）然后在频率域修改傅里叶变换
        dftMask = dftShift * mask  # 修改傅里叶变换实现滤波

        # （3）最后通过傅里叶逆变换返回空间域
        iShift = np.fft.ifftshift(dftMask)  # 将低频逆转换回图像四角
        iDft = cv2.idft(iShift)  # 逆傅里叶变换
        imgRebuild = cv2.magnitude(iDft[:,:,0], iDft[:,:,1])  # 重建图像

        plt.subplot(2,3,i+1), plt.title("Mask (d=)".format(d[i])), plt.axis('off')
        plt.imshow(maskAmp, cmap='gray')
        plt.subplot(2,3,i+4), plt.title("LowPass (d=)".format(d[i])), plt.axis('off')
        plt.imshow(imgRebuild, cmap='gray')

    plt.tight_layout()
    plt.show()

程序说明：

本例程构造了不同尺寸的矩形窗口遮罩，中心低频置 1（白色）四周高频置 0（黑色），是一种低通滤波器。

低通滤波遮罩 mask 与图像傅里叶变换 dftShift 相乘，就使傅里叶变换的高频部分为 0，从而屏蔽原始图像中高频信号，实现了低通滤波。

类似地，将本例程中的低通滤波矩形窗口遮罩反向，改为中心高频置 0（黑色）四周低频置 1（白色），就是一种高通滤波器，可以实现图像锐化和边缘提取。

3.2 频率域图像滤波的步骤

上节例程中通过一个简单的低通滤波遮罩 mask 与图像傅里叶变换 dftShift 相乘，使傅里叶变换的高频部分为 0，从而屏蔽原始图像中高频信号，实现了低通滤波。

频率域图像滤波，首先对原图像 $f(x,y)$ 经傅里叶变换为 $F(u,v)$，然后用适当的滤波器函数 $H(u,v)$ 对傅里叶变换 $F(u,v)$ 的频谱成分进行修改，最后通过傅里叶逆变换（IDFT）返回空间域，得到增强的图像 $g(x,y)$。

$$f(x,y)\stackrel{DFT}{\to }F(u,v)\underset{滤波}{\overset{H(u,v)}{\to }}G(u,v)\stackrel{IDFT}{\to }g(x,y)$$

频率域图像滤波的基本步骤如下：

（1）对原始图像 $f(x,y)$ 进行傅里叶变换，得到 $F(u,v)$；
（2）将图像的傅里叶变换 $F(u,v)$ 与传递函数 $H(u,v)$ 进行卷积运算，得到滤波后的频谱 $G(u,v)$；
（3）对 $G(u,v)$ 进行傅里叶逆变换，得到增强图像 $g(x,y)$。

中心化是将频谱移频到原点，中心化的图像频谱是以原点为圆心，对称分布的。将频谱移频到圆心除了可以清晰地看出图像频率分布以外，还可以分离出有周期性规律的干扰信号。

例程 8.15：频率域图像滤波的一般步骤

    # # 8.15：频率域图像滤波的一般步骤 (本例对应于 冈萨雷斯 数字图像处理第四版 P183 图4.35)
    imgGray = cv2.imread("../images/Fig0431.tif", flags=0)  # flags=0 读取为灰度图像
    H, W = imgGray.shape[:2]  # 图片的高度(H)和宽度(W)
    P, Q = 2*H, 2*W
    # centerY, centerX = int(H/2), int(W/2)  # 图片中心
    fig = plt.figure(figsize=(10, 6))

    # 这里对原图像进行pad，以得到更好的显示
    imgDouble = np.ones([2*H, 2*W]) * 255
    imgDouble[:H, W:] = imgGray # 原始图像位于右上侧
    plt.subplot(241), plt.axis('off'), plt.title("Origin")
    # plt.imshow(imgDouble, cmap='gray')
    plt.imshow(imgGray, cmap='gray')

    fp = np.pad(imgGray, ((0,H), (0,W)), mode='constant')  # 零填充
    plt.subplot(242), plt.axis('off'), plt.title("ZeroFilled")
    plt.imshow(fp, cmap='gray')

    # 中心化
    mask = np.ones(fp.shape)
    mask[1::2, ::2] = -1
    mask[::2, 1::2] = -1
    fpCen = fp * mask
    fpCenShow = np.clip(fpCen, 0, 255)  # 截断函数，将数值限制在 [0,255]
    print("fpCenShow.shape", fpCenShow.shape, fpCenShow.min(), fpCenShow.max())
    plt.subplot(243), plt.axis('off'), plt.title("Centralization")
    plt.imshow(fpCenShow, cmap='gray')

    # cv2.dft 实现图像的傅里叶变换
    imgFloat32 = np.float32(fp)  # 将图像转换成 float32
    dft = cv2.dft(imgFloat32, flags=cv2.DFT_COMPLEX_OUTPUT)  # 傅里叶变换
    dftShift = np.fft.fftshift(dft)  # 将低频分量移动到频域图像的中心
    dftShiftAmp = cv2.magnitude(dftShift[:,:,0], dftShift[:,:,1])  # 幅度谱，中心化
    dftAmpLog = np.log(1 + dftShiftAmp)  # 幅度谱对数变换，以便于显示
    plt.subplot(244), plt.axis('off'), plt.title("DFTshift")
    plt.imshow(dftAmpLog, cmap='gray')

    # Gauss low_pass filter: 1/(2*pi*s2) * exp(-(x**2+y**2)/(2*s2))
    sigma2 = 0.005  # square of sigma
    x, y = np.mgrid[-1:1:1.0/H, -1:1:1.0/W]
    z = 1 / (2 * np.pi * sigma2) * np.exp(-(x**2 + y**2) / (2*sigma2))
    # maskGauss = np.uint8(255 * (z - z.min()) / (z.max() - z.min()))
    maskGauss = np.uint8(cv2.normalize(z, None, 0, 255, cv2.NORM_MINMAX))  # 归一化为 [0,255]
    print("maskGauss.shape", maskGauss.shape, maskGauss.min(), maskGauss.max())
    plt.subplot(245), plt.axis('off'), plt.title("Gauss low-pass")
    plt.imshow(maskGauss, cmap='gray')

    maskGauss2 = np.zeros((maskGauss.shape[0], maskGauss.shape[1], 2), np.uint8)
    maskGauss2[:, :, 0] = maskGauss
    maskGauss2[:, :, 1] = maskGauss
    dftTrans = dftShift * maskGauss2  # 修改傅里叶变换实现滤波
    dftTransAmp = cv2.magnitude(dftTrans[:,:,0], dftTrans[:,:,1])  # 幅度谱，中心化
    dftTransLog = np.log(1 + dftTransAmp)  # 幅度谱对数变换，以便于显示
    plt.subplot(246), plt.axis('off'), plt.title("DFT Trans")
    plt.imshow(dftTransLog, cmap='gray')

    # 通过傅里叶逆变换返回空间域
    ishift = np.fft.ifftshift(dftTrans)  # 将低频逆转换回图像四角
    idft = cv2.idft(ishift)  # 逆傅里叶变换
    idftMag = cv2.magnitude(idft[:,:,0], idft[:,:,1])  # 重建图像
    plt.subplot(247), pltyoucans 的 OpenCV 学习课12. 彩色图像的处理

 youcans 的 OpenCV 学习课1.2 编译生成带有 OpenCV_contrib 的 OpenCV 库
 youcans 的 OpenCV 学习课1.2 编译生成带有 OpenCV_contrib 的 OpenCV 库
 youcans 的 OpenCV 学习课—10. 图像复原与重建
 youcans 的 OpenCV 例程 200 篇103. 陷波带阻滤波器消除周期噪声干扰
 youcans 的图像处理学习课11. 形态学图像处理（中）