【youcans 的图像处理学习课】9. 频率域图像滤波（下）

youcans_

已于 2023-05-18 19:24:19 修改

阅读量1.1k

点赞数

分类专栏： # 图像处理学习课文章标签： opencv python 计算机视觉图像处理算法

于 2022-02-19 07:30:00 首次发布

本文链接：https://blog.csdn.net/youcans/article/details/122991620

版权

图像处理学习课专栏收录该内容

23 篇文章 302 订阅

订阅专栏

专栏地址：『youcans 的图像处理学习课』
文章目录：『youcans 的图像处理学习课 - 总目录』

【youcans 的图像处理学习课】9. 频率域图像滤波（下）

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

文章目录

【youcans 的图像处理学习课】9. 频率域图像滤波（下）

3. 频率域低通滤波器

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

3.1 频率域图像滤波基础

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

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

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

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

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

低频滤波器本质上就是构造一个矩阵，越靠近中心的位置越接近于 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)$ 。

$\xrightarrow {DFT} F(u,v) \xrightarrow [滤波] {H(u,v)} G(u,v) \xrightarrow {IDFT} 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), plt.axis('off'), plt.title("IDFT rebuild")
    plt.imshow(idftMag, cmap='gray')

    imgRebuild = idftMag[:H, :W]
    plt.subplot(248), plt.axis('off'), plt.title("DFT filtered")
    plt.imshow(imgRebuild, cmap='gray')

    plt.tight_layout()
    plt.show()

在这里插入图片描述

例程 8.16：频率域图像滤波的一般步骤（理想低通滤波器）

    # 8.16：频率域图像滤波的一般步骤 (本例对应于 张平 OpenCV算法精解 P383 低通滤波的 Python 实现)
    # (1) 读取原始图像
    # imgGray = cv2.imread("../images/imgLena.tif", flags=0)  # flags=0 读取为灰度图像
    imgGray = cv2.imread("../images/Fig0431.tif", flags=0)  # flags=0 读取为灰度图像
    imgFloat32 = np.float32(imgGray)  # 将图像转换成 float32
    rows, cols = imgGray.shape[:2]  # 图片的高度和宽度

    fig = plt.figure(figsize=(10, 6))
    plt.subplot(241), plt.axis('off'), plt.title("Origin image")
    plt.imshow(imgGray, cmap='gray')

    # (2) 中心化, centralized 2d array f(x,y) * (-1)^(x+y)
    mask = np.ones(imgGray.shape)
    mask[1::2, ::2] = -1
    mask[::2, 1::2] = -1
    fImage = imgFloat32 * mask  # f(x,y) * (-1)^(x+y)
    plt.subplot(242), plt.axis('off'), plt.title("f(x,y)*(-1)^(x+y)")
    plt.imshow(fImage, cmap='gray')

    # (3) 快速傅里叶变换
    rPadded = cv2.getOptimalDFTSize(rows)  # 最优 DFT 扩充尺寸
    cPadded = cv2.getOptimalDFTSize(cols)  # 用于快速傅里叶变换
    dftImage = np.zeros((rPadded, cPadded, 2), np.float32)  # 对原始图像进行边缘扩充
    dftImage[:rows, :cols, 0] = fImage  # 边缘扩充，下侧和右侧补0
    cv2.dft(dftImage, dftImage, cv2.DFT_COMPLEX_OUTPUT)  # 快速傅里叶变换

    dftAmp = cv2.magnitude(dftImage[:,:,0], dftImage[:,:,1])  # 傅里叶变换的幅度谱 (rPad, cPad)
    dftAmpLog = np.log(1.0 + dftAmp)  # 幅度谱对数变换，以便于显示
    dftAmpNorm = np.uint8(cv2.normalize(dftAmpLog, None, 0, 255, cv2.NORM_MINMAX))  # 归一化为 [0,255]
    plt.subplot(243), plt.axis('off'), plt.title("DFT spectrum")
    plt.imshow(dftAmpNorm, cmap='gray')

    # (4) 构建低通滤波器 传递函数
    # 找到傅里叶谱最大值的位置
    minValue, maxValue, minLoc, maxLoc = cv2.minMaxLoc(dftAmp)
    rPadded, cPadded = dftImage.shape[:2]  # 快速傅里叶变换的尺寸, 原始图像尺寸优化
    u, v = np.mgrid[0:rPadded:1, 0:cPadded:1]
    Dsquare = np.power((u-maxLoc[1]), 2.0) + np.power((v-maxLoc[0]), 2.0)
    D0 = 50  # radius
    D = np.sqrt(Dsquare)
    lpFilter = np.zeros(dftImage.shape[:2], np.float32)
    lpFilter[D <= D0] = 1  # 理想低通滤波 (Idea low pass filter)
    plt.subplot(244), plt.axis('off'), plt.title("LP Filter")
    plt.imshow(lpFilter, cmap='gray')

    # (5) 在频率域修改傅里叶变换: 傅里叶变换 点乘 低通滤波器
    # rPadded, cPadded = dftImage.shape[:2]  # 快速傅里叶变换的尺寸, 原始图像尺寸优化
    dftLPfilter = np.zeros(dftImage.shape, dftImage.dtype)  # 快速傅里叶变换的尺寸(优化尺寸)
    for i in range(2):
        dftLPfilter[:rPadded, :cPadded, i] = dftImage[:rPadded, :cPadded, i] * lpFilter

    # 低通傅里叶变换的傅里叶谱
    lpDftAmp = cv2.magnitude(dftLPfilter[:, :, 0], dftLPfilter[:, :, 1])  # 傅里叶变换的幅度谱
    lpDftAmpLog = np.log(1.0 + lpDftAmp)  # 幅度谱对数变换，以便于显示
    lpDftAmpNorm = np.uint8(cv2.normalize(lpDftAmpLog, None, 0, 255, cv2.NORM_MINMAX))  # 归一化为 [0,255]
    plt.subplot(245), plt.axis('off'), plt.title("lpfDft Spectrum")
    plt.imshow(lpDftAmpNorm, cmap='gray')

    # (6) 对低通傅里叶变换 执行傅里叶逆变换，并只取实部
    idft = np.zeros(dftAmp.shape, np.float32)  # 快速傅里叶变换的尺寸(优化尺寸)
    cv2.dft(dftLPfilter, idft, cv2.DFT_REAL_OUTPUT + cv2.DFT_INVERSE + cv2.DFT_SCALE)
    plt.subplot(246), plt.axis('off'), plt.title("IDFT Spectrum")
    plt.imshow(idft, cmap='gray')

    # (7) 中心化, centralized 2d array g(x,y) * (-1)^(x+y)
    mask2 = np.ones(dftAmp.shape)
    mask2[1::2, ::2] = -1
    mask2[::2, 1::2] = -1
    idftCen = idft * mask2  # g(x,y) * (-1)^(x+y)
    plt.subplot(247), plt.axis('off'), plt.title("g(x,y)*(-1)^(x+y)")
    plt.imshow(idftCen, cmap='gray')

    # (8) 截取左上角，大小和输入图像相等
    result = np.clip(idftCen, 0, 255)  # 截断函数，将数值限制在 [0,255]
    lpResult = result.astype(np.uint8)
    lpResult = lpResult[:rows, :cols]

    plt.subplot(248), plt.axis('off'), plt.title("DFT filtered")
    plt.imshow(lpResult, cmap='gray')
    plt.tight_layout()
    plt.show()

    print("image.shape:{}".format(imgGray.shape))
    print("imgFloat32.shape:{}".format(imgFloat32.shape))
    print("dftImage.shape:{}".format(dftImage.shape))
    print("dftAmp.shape:{}".format(dftAmp.shape))
    print("idft.shape:{}".format(idft.shape))
    print("lpFilter.shape:{}".format(lpFilter.shape))
    print("result.shape:{}".format(result.shape))

在这里插入图片描述
程序说明：

（1）使用 cv2.getOptimalDFTSize() 获得快速傅里叶变换的优化尺寸，对原始图像进行了边缘扩充和补 0 填充，因此与原始图像尺寸不一定相同。频域滤波后截取左上角，得到与原始图像大小一致的滤波图像。
（2）注意滤波器的尺寸要与（快速）傅里叶变换的尺寸相同，与原始图像的大小不一定相同。

3.3 频率域高斯低通滤波器（GLPF）

例程 8.16 以理想低通滤波器为例，给出了频率域图像滤波的一般步骤。

高斯滤波器的滤波效果随着样值点到傅里叶变换中心的距离的变化而变化，当距离增大时滤波效果变好。

频率域高斯低通滤波器也是一个掩模蒙板：

$e^{-(u^2+v^2)/2 \sigma^2}\\ D(u,v) = \sqrt {(u-P/2)^2+(v-Q/2)^2}$
式中 $D_0=\sigma$ 是截止频率，当 $D(u,v) = D_0$ 时 GLPF 下降到最大值的 0.607处。

利用二维高斯蒙板，可以滤去高频、保留低频。频域高斯函数的方差 $\sigma^2$ 越小，高斯函数越狭窄，滤除的高频成分（图像细节）越多，图像越模糊。

空间域高斯滤波器的计算量随着模板的增大而增大，而频率域高斯滤波的计算量独立于滤波函数。

例程 8.18：频率域高斯低通滤波器 (GLPF)

    # 8.18：频率域高斯低通滤波器 (GLPF)
    imgGray = cv2.imread("../images/imgLena.tif", flags=0)  # flags=0 读取为灰度图像

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

    # # 幅度谱
    # imgDFT = 10 * np.log(cv2.magnitude(dftShift[:,:,0], dftShift[:,:,1]))
    # plt.figure(figsize=(9, 6))
    # plt.subplot(121), plt.axis('off'), plt.title("Original"), plt.imshow(imgGray, 'gray')
    # plt.subplot(122), plt.axis('off'), plt.title("FFt"), plt.imshow(imgDFT, cmap='gray')
    # plt.show()

    plt.figure(figsize=(9, 6))
    rows, cols = imgGray.shape[:2]  # 图片的高度和宽度
    sigma2 = [0.5, 0.09, 0.01]  # square of sigma
    for i in range(3):
        # 构造高斯滤波器遮罩：# Gauss = 1/(2*pi*s2) * exp(-(x**2+y**2)/(2*s2))
        x, y = np.mgrid[-1:1:2.0 / rows, -1:1:2.0 / cols]
        z = 1 / (2 * np.pi * sigma2[i]) * np.exp(-(x ** 2 + y ** 2) / (2 * sigma2[i]))
        # maskGauss = np.uint8(cv2.normalize(z, None, 0, 255, cv2.NORM_MINMAX))  # 归一化为 [0,255]
        zNorm = np.uint8(cv2.normalize(z, None, 0, 255, cv2.NORM_MINMAX))  # 归一化为 [0,255]
        maskGauss = np.zeros((rows, cols, 2), np.uint8)
        maskGauss[:,:,0] = zNorm
        maskGauss[:,:,1] = zNorm
        print(maskGauss.shape, maskGauss.max(), maskGauss.min())

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

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

        plt.subplot(2,3,i+1), plt.title("Mask (s^2={})".format(sigma2[i])), plt.axis('off')
        plt.imshow(maskGauss[:,:,0], cmap='gray')
        plt.subplot(2,3,i+4), plt.title("DFT GLPF (s^2={})".format(sigma2[i])), plt.axis('off')
        plt.imshow(imgRebuild, cmap='gray')

    plt.tight_layout()
    plt.show()

在这里插入图片描述

3.4 频率域巴特沃斯低通滤波器（BLPF）

截止频率位于距频率中心 $D_0$ 处的 n 阶巴特沃斯（Butterworth）低通滤波器的传递函数为：
$\frac {1}{1+[D(u,v) / D_0]^{2n}}$

当 n 较大时，巴特沃斯低通滤波器 BLPF 可以逼近理想低通滤波器 ILPF 的特性；而当 n 较小时，巴特沃斯低通滤波器 BLPF 可以逼近高斯低通滤波器 GLPF 的特性，同时提供从低频到高频的平滑过渡。

利用二维高斯蒙板，可以滤去高频、保留低频。频域高斯函数的越小，高斯函数越狭窄，滤除的高频成分（图像细节）越多，图像越模糊。

巴特沃斯滤波器的特点是通频带内的频率响应曲线最大限度平坦，没有起伏，而在阻频带则逐渐下降为零。
在振幅的对数对角频率的波特图上，从某一边界角频率开始，振幅随着角频率的增加而逐步减小，趋向负无穷大。巴特沃斯滤波器的频率特性曲线，无论在通带内还是阻带内都是频率的单调函数。因此，当通带的边界处满足指标要求时，通带内肯定会有裕量。所以，更有效的设计方法应该是将精确度均匀的分布在整个通带或阻带内，或者同时分布在两者之内。

例程 8.19：频率域巴特沃斯低通滤波器 (BLPF)

    # 8.19：频率域巴特沃斯低通滤波器 (BLPF)
    # (1) 读取原始图像
    # imgGray = cv2.imread("../images/imgLena.tif", flags=0)  # flags=0 读取为灰度图像
    imgGray = cv2.imread("../images/Fig0431.tif", flags=0)  # flags=0 读取为灰度图像
    imgFloat32 = np.float32(imgGray)  # 将图像转换成 float32
    rows, cols = imgGray.shape[:2]  # 图片的高度和宽度

    # (2) 中心化, centralized 2d array f(x,y) * (-1)^(x+y)
    mask = np.ones(imgGray.shape)
    mask[1::2, ::2] = -1
    mask[::2, 1::2] = -1
    fImage = imgFloat32 * mask  # f(x,y) * (-1)^(x+y)

    # (3) 快速傅里叶变换
    # dftImage = fft2Image(fImage)  # 快速傅里叶变换 (rPad, cPad, 2)
    rPadded = cv2.getOptimalDFTSize(rows)  # 最优 DFT 扩充尺寸
    cPadded = cv2.getOptimalDFTSize(cols)  # 用于快速傅里叶变换
    dftImage = np.zeros((rPadded, cPadded, 2), np.float32)  # 对原始图像进行边缘扩充
    dftImage[:rows, :cols, 0] = fImage  # 边缘扩充，下侧和右侧补0
    cv2.dft(dftImage, dftImage, cv2.DFT_COMPLEX_OUTPUT)  # 快速傅里叶变换

    dftAmp = cv2.magnitude(dftImage[:,:,0], dftImage[:,:,1])  # 傅里叶变换的幅度谱 (rPad, cPad)
    dftAmpLog = np.log(1.0 + dftAmp)  # 幅度谱对数变换，以便于显示
    dftAmpNorm = np.uint8(cv2.normalize(dftAmpLog, None, 0, 255, cv2.NORM_MINMAX))  # 归一化为 [0,255]
    minValue, maxValue, minLoc, maxLoc = cv2.minMaxLoc(dftAmp)  # 找到傅里叶谱最大值的位置

    plt.figure(figsize=(9, 6))
    # rows, cols = imgGray.shape[:2]  # 图片的高度和宽度
    u, v = np.mgrid[0:rPadded:1, 0:cPadded:1]
    D = np.sqrt(np.power((u-maxLoc[1]), 2) + np.power((v-maxLoc[0]), 2))
    D0 = [20, 40, 80]  # cut-off frequency
    n = 2
    for k in range(3):
        # (4) 构建低通滤波器 传递函数
        # 巴特沃斯低通滤波 (Butterworth low pass filter)
        epsilon = 1e-8  # 防止被 0 除
        lpFilter = 1.0 / (1.0 + np.power(D / (D0[k] + epsilon), 2*n))

        # (5) 在频率域修改傅里叶变换: 傅里叶变换 点乘 低通滤波器
        dftLPfilter = np.zeros(dftImage.shape, dftImage.dtype)  # 快速傅里叶变换的尺寸(优化尺寸)
        for j in range(2):
            dftLPfilter[:rPadded, :cPadded, j] = dftImage[:rPadded, :cPadded, j] * lpFilter

        # (6) 对低通傅里叶变换 执行傅里叶逆变换，并只取实部
        idft = np.zeros(dftAmp.shape, np.float32)  # 快速傅里叶变换的尺寸(优化尺寸)
        cv2.dft(dftLPfilter, idft, cv2.DFT_REAL_OUTPUT + cv2.DFT_INVERSE + cv2.DFT_SCALE)

        # (7) 中心化, centralized 2d array g(x,y) * (-1)^(x+y)
        mask2 = np.ones(dftAmp.shape)
        mask2[1::2, ::2] = -1
        mask2[::2, 1::2] = -1
        idftCen = idft * mask2  # g(x,y) * (-1)^(x+y)

        # (8) 截取左上角，大小和输入图像相等
        result = np.clip(idftCen, 0, 255)  # 截断函数，将数值限制在 [0,255]
        imgBLPF = result.astype(np.uint8)
        imgBLPF = imgBLPF[:rows, :cols]

        plt.subplot(2,3,k+1), plt.title("BLPF mask(D0={})".format(D0[k])), plt.axis('off')
        plt.imshow(lpFilter[:,:], cmap='gray')
        plt.subplot(2,3,k+4), plt.title("BLPF rebuild(D0={})".format(D0[k])), plt.axis('off')
        plt.imshow(imgBLPF, cmap='gray')

    plt.tight_layout()
    plt.show()

在这里插入图片描述

3.5 频率域低通滤波案例：印刷文本字符修复

低通滤波技术主要应用于印刷和出版业。本节给出频率域低通滤波的实际应用：印刷文本的字符修复，案例来自冈萨雷斯《数字图像处理（第四版）》P194 图 4.48。

原始图像是一幅低分辨率文本图像的样本，我们会在传真、复印或历史文档中遇到这样的文本。图片中的文本不存在污迹、折痕和撕裂，但文本中的字符因分辨率不足而产生的失真，其中许多字符是断开的。解决这个问题的一种方法是，通过模糊输入图像来桥接这些小缝隙。

例程显示了使用不同 $D_0$ 的高斯低通滤波器简单处理后的修复效果。当然这只是修复的一个步骤，接下来还要进行其它处理如阈值处理、细化处理，这些方法将在后续内容讨论。

为了方便今后的使用，例程将最优扩充的快速傅立叶变换、高斯低通滤波器封装为函数。

    # 8.22：GLPF 低通滤波案例：印刷文本字符修复 (Text character recognition)
    def gaussLowPassFilter(shape, radius=10):  # 高斯低通滤波器
        # 高斯滤波器：# Gauss = 1/(2*pi*s2) * exp(-(x**2+y**2)/(2*s2))
        u, v = np.mgrid[-1:1:2.0/shape[0], -1:1:2.0/shape[1]]
        D = np.sqrt(u**2 + v**2)
        D0 = radius / shape[0]
        kernel = np.exp(- (D ** 2) / (2 *D0**2))
        return kernel

    def dft2Image(image):  #　最优扩充的快速傅立叶变换
        # 中心化, centralized 2d array f(x,y) * (-1)^(x+y)
        mask = np.ones(image.shape)
        mask[1::2, ::2] = -1
        mask[::2, 1::2] = -1
        fImage = image * mask  # f(x,y) * (-1)^(x+y)

        # 最优 DFT 扩充尺寸
        rows, cols = image.shape[:2]  # 原始图片的高度和宽度
        rPadded = cv2.getOptimalDFTSize(rows)  # 最优 DFT 扩充尺寸
        cPadded = cv2.getOptimalDFTSize(cols)  # 用于快速傅里叶变换

        # 边缘扩充(补0), 快速傅里叶变换
        dftImage = np.zeros((rPadded, cPadded, 2), np.float32)  # 对原始图像进行边缘扩充
        dftImage[:rows, :cols, 0] = fImage  # 边缘扩充，下侧和右侧补0
        cv2.dft(dftImage, dftImage, cv2.DFT_COMPLEX_OUTPUT)  # 快速傅里叶变换
        return dftImage


    # (1) 读取原始图像
    imgGray = cv2.imread("../images/Fig0419a.tif", flags=0)  # flags=0 读取为灰度图像
    rows, cols = imgGray.shape[:2]  # 图片的高度和宽度

    plt.figure(figsize=(10, 6))
    plt.subplot(2, 3, 1), plt.title("Original"), plt.axis('off'), plt.imshow(imgGray, cmap='gray')

    # (2) 快速傅里叶变换
    dftImage = dft2Image(imgGray)  # 快速傅里叶变换 (rPad, cPad, 2)
    rPadded, cPadded = dftImage.shape[:2]  # 快速傅里叶变换的尺寸, 原始图像尺寸优化
    print("dftImage.shape:{}".format(dftImage.shape))


    D0 = [10, 30, 60, 90, 120]  # radius
    for k in range(5):
        # (3) 构建 高斯低通滤波 (Gauss low pass filter)
        lpFilter = gaussLowPassFilter((rPadded, cPadded), radius=D0[k])

        # (5) 在频率域修改傅里叶变换: 傅里叶变换 点乘 低通滤波器
        dftLPfilter = np.zeros(dftImage.shape, dftImage.dtype)  # 快速傅里叶变换的尺寸(优化尺寸)
        for j in range(2):
            dftLPfilter[:rPadded, :cPadded, j] = dftImage[:rPadded, :cPadded, j] * lpFilter

        # (6) 对低通傅里叶变换 执行傅里叶逆变换，并只取实部
        idft = np.zeros(dftImage.shape[:2], np.float32)  # 快速傅里叶变换的尺寸(优化尺寸)
        cv2.dft(dftLPfilter, idft, cv2.DFT_REAL_OUTPUT + cv2.DFT_INVERSE + cv2.DFT_SCALE)

        # (7) 中心化, centralized 2d array g(x,y) * (-1)^(x+y)
        mask2 = np.ones(dftImage.shape[:2])
        mask2[1::2, ::2] = -1
        mask2[::2, 1::2] = -1
        idftCen = idft * mask2  # g(x,y) * (-1)^(x+y)

        # (8) 截取左上角，大小和输入图像相等
        result = np.clip(idftCen, 0, 255)  # 截断函数，将数值限制在 [0,255]
        imgLPF = result.astype(np.uint8)
        imgLPF = imgLPF[:rows, :cols]

        plt.subplot(2,3,k+2), plt.title("GLPF rebuild(n={})".format(D0[k])), plt.axis('off')
        plt.imshow(imgLPF, cmap='gray')

    print("image.shape:{}".format(imgGray.shape))
    print("lpFilter.shape:{}".format(lpFilter.shape))
    print("dftImage.shape:{}".format(dftImage.shape))

    plt.tight_layout()
    plt.show()

在这里插入图片描述

4. 频率域高通滤波器

图像边缘化其它灰度的急剧变化与高频分量有关，因此可以在频率域通过高通滤波实现图像锐化。高通滤波衰减傅里叶变换中的低频分量而不干扰高频信息。

4.1 由低通滤波器得到高通滤波器

简单地，在频率域中用 1 减去低通滤波器的传递函数，就可以得到相应的高通滤波器传递函数：
$H_{HP}(u,v) = 1- H_{LP}(u,v)$
式中， $H_{HP}(u,v)$ 、 $H_{LP}(u,v)$ 分别表示高通滤波器、低通滤波器的传递函数。

理想高通滤波器（IHPF）的传递函数为：
$H(u,v)=\begin{cases} 0,\ D(u,v) \leq D_0\\ 1,\ D(u,v)>D_0 \end{cases}$
高斯高通滤波器（GHPF）的传递函数为：
$H(u,v)=1-e^{-D^2 (u,v)/2D_0^2}$
巴特沃斯高通滤波器（BHPF）的传递函数为：
$\frac {1} {1 + [D_0/D(u,v)]^{2n}}$

例程 8.23 由低通滤波器得到高通滤波器

    # 8.23：频率域高通滤波器
    def ideaHighPassFilter(shape, radius=10):  # 理想高通滤波器
        u, v = np.mgrid[-1:1:2.0/shape[0], -1:1:2.0/shape[1]]
        D = np.sqrt(u**2 + v**2)
        D0 = radius / shape[0]
        kernel = np.ones(shape)
        kernel[D <= D0] = 0  # 理想低通滤波 (Idea low pass filter)
        return kernel

    def gaussHighPassFilter(shape, radius=10):  # 高斯高通滤波器
        # 高斯滤波器：# Gauss = 1/(2*pi*s2) * exp(-(x**2+y**2)/(2*s2))
        u, v = np.mgrid[-1:1:2.0/shape[0], -1:1:2.0/shape[1]]
        D = np.sqrt(u**2 + v**2)
        D0 = radius / shape[0]
        kernel = 1 - np.exp(- (D ** 2) / (2 *D0**2))
        return kernel

    def butterworthHighPassFilter(shape, radius=10, n=2):  # 巴特沃斯高通滤波
        u, v = np.mgrid[-1:1:2.0/shape[0], -1:1:2.0/shape[1]]
        epsilon = 1e-8
        D = np.sqrt(u**2 + v**2)
        D0 = radius / shape[0]
        kernel = 1.0 / (1.0 + np.power(D0/(D + epsilon), 2*n))

        return kernel

    # 理想、高斯、巴特沃斯高通传递函数
    shape = [128, 128]
    radius = 32
    IHPF = ideaHighPassFilter(shape, radius=radius)
    GHPF = gaussHighPassFilter(shape, radius=radius)
    BHPF = butterworthHighPassFilter(shape, radius=radius)

    filters = ['IHPF', 'GHPF', 'BHPF']
    u, v = np.mgrid[-1:1:2.0/shape[0], -1:1:2.0/shape[1]]
    fig = plt.figure(figsize=(10, 8))
    for i in range(3):
        hpFilter = eval(filters[i]).copy()

        ax1 = fig.add_subplot(3, 3, 3*i+1)
        ax1.imshow(hpFilter, 'gray')
        ax1.set_title(filters[i]), ax1.set_xticks([]), ax1.set_yticks([])

        ax2 = plt.subplot(3,3,3*i+2, projection='3d')
        ax2.set_title("transfer function")
        ax2.plot_wireframe(u, v, hpFilter , rstride=2, linewidth=0.5, color='c')
        ax2.set_xticks([]), ax2.set_yticks([]), ax2.set_zticks([])

        ax3 = plt.subplot(3,3,3*i+3)
        profile = hpFilter[shape[0]//2:, shape[1]//2]
        ax3.plot(profile), ax3.set_title("profile"), ax3.set_xticks([]), ax3.set_yticks([])

    plt.show()

在这里插入图片描述

例程 8.24 频率域高通滤波器的应用

    # 8.24：频率域高通滤波器的应用
    def ideaHighPassFilter(shape, radius=10):  # 理想高通滤波器
        u, v = np.mgrid[-1:1:2.0/shape[0], -1:1:2.0/shape[1]]
        D = np.sqrt(u**2 + v**2)
        D0 = radius / shape[0]
        kernel = np.ones(shape)
        kernel[D <= D0] = 0  # 理想低通滤波 (Idea low pass filter)
        return kernel

    def gaussHighPassFilter(shape, radius=10):  # 高斯高通滤波器
        # 高斯滤波器：# Gauss = 1/(2*pi*s2) * exp(-(x**2+y**2)/(2*s2))
        u, v = np.mgrid[-1:1:2.0/shape[0], -1:1:2.0/shape[1]]
        D = np.sqrt(u**2 + v**2)
        D0 = radius / shape[0]
        kernel = 1 - np.exp(- (D ** 2) / (2 *D0**2))
        return kernel

    def butterworthHighPassFilter(shape, radius=10, n=2):  # 巴特沃斯高通滤波
        u, v = np.mgrid[-1:1:2.0/shape[0], -1:1:2.0/shape[1]]
        epsilon = 1e-8
        D = np.sqrt(u**2 + v**2)
        D0 = radius / shape[0]
        kernel = 1.0 / (1.0 + np.power(D0/(D + epsilon), 2*n))

        return kernel

    def dft2Image(image):  #　最优扩充的快速傅立叶变换
        # 中心化, centralized 2d array f(x,y) * (-1)^(x+y)
        mask = np.ones(image.shape)
        mask[1::2, ::2] = -1
        mask[::2, 1::2] = -1
        fImage = image * mask  # f(x,y) * (-1)^(x+y)

        # 最优 DFT 扩充尺寸
        rows, cols = image.shape[:2]  # 原始图片的高度和宽度
        rPadded = cv2.getOptimalDFTSize(rows)  # 最优 DFT 扩充尺寸
        cPadded = cv2.getOptimalDFTSize(cols)  # 用于快速傅里叶变换

        # 边缘扩充(补0), 快速傅里叶变换
        dftImage = np.zeros((rPadded, cPadded, 2), np.float32)  # 对原始图像进行边缘扩充
        dftImage[:rows, :cols, 0] = fImage  # 边缘扩充，下侧和右侧补0
        cv2.dft(dftImage, dftImage, cv2.DFT_COMPLEX_OUTPUT)  # 快速傅里叶变换
        return dftImage


    # (1) 读取原始图像
    imgGray = cv2.imread("../images/Fig0515a.tif", flags=0)  # flags=0 读取为灰度图像
    # imgGray = cv2.imread("../images/imgLena.tif", flags=0)  # flags=0 读取为灰度图像
    rows, cols = imgGray.shape[:2]  # 图片的高度和宽度

    plt.figure(figsize=(10, 6))
    plt.subplot(2, 3, 1), plt.title("Original"), plt.axis('off'), plt.imshow(imgGray, cmap='gray')

    # (2) 快速傅里叶变换
    dftImage = dft2Image(imgGray)  # 快速傅里叶变换 (rPad, cPad, 2)
    rPadded, cPadded = dftImage.shape[:2]  # 快速傅里叶变换的尺寸, 原始图像尺寸优化
    print("dftImage.shape:{}".format(dftImage.shape))


    D0 = [20, 40, 80, 120, 160]  # radius
    for k in range(5):
        # (3) 构建高通滤波器
        # hpFilter = ideaHighPassFilter((rPadded, cPadded), radius=D0[k])  # 理想高通滤波器
        # hpFilter = gaussHighPassFilter((rPadded, cPadded), radius=D0[k])  # 高斯高通滤波器
        hpFilter = butterworthHighPassFilter((rPadded, cPadded), radius=D0[k])  # 巴特沃斯高通滤波器

        # (5) 在频率域修改傅里叶变换: 傅里叶变换 点乘 低通滤波器
        dftHPfilter = np.zeros(dftImage.shape, dftImage.dtype)  # 快速傅里叶变换的尺寸(优化尺寸)
        for j in range(2):
            dftHPfilter[:rPadded, :cPadded, j] = dftImage[:rPadded, :cPadded, j] * hpFilter

        # (6) 对高通傅里叶变换 执行傅里叶逆变换，并只取实部
        idft = np.zeros(dftImage.shape[:2], np.float32)  # 快速傅里叶变换的尺寸(优化尺寸)
        cv2.dft(dftHPfilter, idft, cv2.DFT_REAL_OUTPUT + cv2.DFT_INVERSE + cv2.DFT_SCALE)

        # (7) 中心化, centralized 2d array g(x,y) * (-1)^(x+y)
        mask2 = np.ones(dftImage.shape[:2])
        mask2[1::2, ::2] = -1
        mask2[::2, 1::2] = -1
        idftCen = idft * mask2  # g(x,y) * (-1)^(x+y)

        # (8) 截取左上角，大小和输入图像相等
        result = np.clip(idftCen, 0, 255)  # 截断函数，将数值限制在 [0,255]
        imgHPF = result.astype(np.uint8)
        imgHPF = imgLPF[:rows, :cols]

        plt.subplot(2,3,k+2), plt.title("HPFilter rebuild(n={})".format(D0[k])), plt.axis('off')
        plt.imshow(imgHPF, cmap='gray')

    print("image.shape:{}".format(imgGray.shape))
    print("hpFilter.shape:{}".format(hpFilter.shape))
    print("dftImage.shape:{}".format(dftImage.shape))

    plt.tight_layout()
    plt.show()

在这里插入图片描述

例程 8.25：指纹图像处理（高通滤波+阈值处理）

    # 8.25：指纹图像处理（高通滤波+阈值处理）
    def gaussHighPassFilter(shape, radius=10):  # 高斯高通滤波器
        # 高斯滤波器：# Gauss = 1/(2*pi*s2) * exp(-(x**2+y**2)/(2*s2))
        u, v = np.mgrid[-1:1:2.0/shape[0], -1:1:2.0/shape[1]]
        D = np.sqrt(u**2 + v**2)
        D0 = radius / shape[0]
        kernel = 1 - np.exp(- (D ** 2) / (2 *D0**2))
        return kernel

    def dft2Image(image):  #　最优扩充的快速傅立叶变换
        # 中心化, centralized 2d array f(x,y) * (-1)^(x+y)
        mask = np.ones(image.shape)
        mask[1::2, ::2] = -1
        mask[::2, 1::2] = -1
        fImage = image * mask  # f(x,y) * (-1)^(x+y)

        # 最优 DFT 扩充尺寸
        rows, cols = image.shape[:2]  # 原始图片的高度和宽度
        rPadded = cv2.getOptimalDFTSize(rows)  # 最优 DFT 扩充尺寸
        cPadded = cv2.getOptimalDFTSize(cols)  # 用于快速傅里叶变换

        # 边缘扩充(补0), 快速傅里叶变换
        dftImage = np.zeros((rPadded, cPadded, 2), np.float32)  # 对原始图像进行边缘扩充
        dftImage[:rows, :cols, 0] = fImage  # 边缘扩充，下侧和右侧补0
        cv2.dft(dftImage, dftImage, cv2.DFT_COMPLEX_OUTPUT)  # 快速傅里叶变换
        return dftImage


    def imgHPFilter(image, D0=50):  #　图像高通滤波
        rows, cols = image.shape[:2]  # 图片的高度和宽度
        # 快速傅里叶变换
        dftImage = dft2Image(image)  # 快速傅里叶变换 (rPad, cPad, 2)
        rPadded, cPadded = dftImage.shape[:2]  # 快速傅里叶变换的尺寸, 原始图像尺寸优化

        # 构建 高斯高通滤波器 (Gauss low pass filter)
        hpFilter = gaussHighPassFilter((rPadded, cPadded), radius=D0)  # 高斯高通滤波器

        # 在频率域修改傅里叶变换: 傅里叶变换 点乘 高通滤波器
        dftHPfilter = np.zeros(dftImage.shape, dftImage.dtype)  # 快速傅里叶变换的尺寸(优化尺寸)
        for j in range(2):
            dftHPfilter[:rPadded, :cPadded, j] = dftImage[:rPadded, :cPadded, j] * hpFilter

        # 对高通傅里叶变换 执行傅里叶逆变换，并只取实部
        idft = np.zeros(dftImage.shape[:2], np.float32)  # 快速傅里叶变换的尺寸(优化尺寸)
        cv2.dft(dftHPfilter, idft, cv2.DFT_REAL_OUTPUT + cv2.DFT_INVERSE + cv2.DFT_SCALE)

        # 中心化, centralized 2d array g(x,y) * (-1)^(x+y)
        mask2 = np.ones(dftImage.shape[:2])
        mask2[1::2, ::2] = -1
        mask2[::2, 1::2] = -1
        idftCen = idft * mask2  # g(x,y) * (-1)^(x+y)

        # 截取左上角，大小和输入图像相等
        result = np.clip(idftCen, 0, 255)  # 截断函数，将数值限制在 [0,255]
        imgHPF = result.astype(np.uint8)
        imgHPF = imgHPF[:rows, :cols]
        return imgHPF


    imgGray = cv2.imread("../images/Fig0457a.tif", flags=0)  # flags=0 读取为灰度图像
    rows, cols = imgGray.shape[:2]  # 图片的高度和宽度
    imgHPF = imgHPFilter(imgGray, D0=50)
    imgThres = np.clip(imgHPF, 0, 1)

    plt.figure(figsize=(10, 5))
    plt.subplot(131), plt.imshow(imgGray, 'gray'), plt.title('origial'), plt.xticks([]), plt.yticks([])
    plt.subplot(132), plt.imshow(imgHPF, 'gray'), plt.title('GaussHPF'), plt.xticks([]), plt.yticks([])
    plt.subplot(133), plt.imshow(imgThres, 'gray'), plt.title('Threshold'), plt.xticks([]), plt.yticks([])
    plt.tight_layout()
    plt.show()

在这里插入图片描述

4.2 频率域钝化掩蔽

简单地，从原始图像中减去一幅平滑处理的钝化图像，也可以实现图像锐化效果，称为钝化掩蔽。

令 $f_{LP}(x,y)$ 表示低通滤波的平滑图像，则：
$\begin{align} g_{mask} (x,y) &= f(x,y) - f_{LP}(x,y) \\ g(x,y) &= f(x,y) + k * g_{mask}(x,y)， k>0 \end{align}$
当 k>1 时，实现高提升滤波；当 k=1 时，实现钝化掩蔽；k<1时，可以减弱钝化掩蔽的强度。

原图减去模糊图的结果为模板，输出图像等于原图加上加权后的模板，当权重为1得到非锐化掩蔽，当权重大于1成为高提升滤波。

在频率域实现钝化掩蔽，高频强调滤波器传递函数为：
$g(x,y)= J^{-1} \{ [1+k H_{HP}(u,v)]F(u,v) \}$
高频强调滤波的通用公式是：
$g(x,y)= J^{-1} \{ [k_1 + k_2 H_{HP}(u,v)]F(u,v) \}$
式中， $k_1 \geq 0$ 偏移传递函数的值， $k_2 \geq 0$ 控制高频的贡献。

例程 8.26：频率域钝化掩蔽

    # 8.26：频率域钝化掩蔽
    def gaussHighPassFilter(shape, radius=10):  # 高斯高通滤波器
        # 高斯滤波器：# Gauss = 1/(2*pi*s2) * exp(-(x**2+y**2)/(2*s2))
        u, v = np.mgrid[-1:1:2.0/shape[0], -1:1:2.0/shape[1]]
        D = np.sqrt(u**2 + v**2)
        D0 = radius / shape[0]
        kernel = 1 - np.exp(- (D ** 2) / (2 *D0**2))
        return kernel

    def dft2Image(image):  #　最优扩充的快速傅立叶变换
        # 中心化, centralized 2d array f(x,y) * (-1)^(x+y)
        mask = np.ones(image.shape)
        mask[1::2, ::2] = -1
        mask[::2, 1::2] = -1
        fImage = image * mask  # f(x,y) * (-1)^(x+y)

        # 最优 DFT 扩充尺寸
        rows, cols = image.shape[:2]  # 原始图片的高度和宽度
        rPadded = cv2.getOptimalDFTSize(rows)  # 最优 DFT 扩充尺寸
        cPadded = cv2.getOptimalDFTSize(cols)  # 用于快速傅里叶变换

        # 边缘扩充(补0), 快速傅里叶变换
        dftImage = np.zeros((rPadded, cPadded, 2), np.float32)  # 对原始图像进行边缘扩充
        dftImage[:rows, :cols, 0] = fImage  # 边缘扩充，下侧和右侧补0
        cv2.dft(dftImage, dftImage, cv2.DFT_COMPLEX_OUTPUT)  # 快速傅里叶变换
        return dftImage


    # 高频强调滤波 + 直方图均衡化
    image = cv2.imread("../images/Fig0459a.tif", flags=0)  # flags=0 读取为灰度图像
    rows, cols = image.shape[:2]  # 图片的高度和宽度
    print(rows, cols)

    # 快速傅里叶变换
    dftImage = dft2Image(image)  # 快速傅里叶变换 (rPad, cPad, 2)
    rPadded, cPadded = dftImage.shape[:2]  # 快速傅里叶变换的尺寸, 原始图像尺寸优化
    # 构建 高斯高通滤波器 (Gauss low pass filter)
    hpFilter = gaussHighPassFilter((rPadded, cPadded), radius=10)  # 高斯高通滤波器
    # 在频率域修改傅里叶变换: 傅里叶变换 点乘 低通滤波器
    dftHPfilter = np.zeros(dftImage.shape, dftImage.dtype)  # 快速傅里叶变换的尺寸(优化尺寸)
    for j in range(2):
        dftHPfilter[:rPadded, :cPadded, j] = dftImage[:rPadded, :cPadded, j] * hpFilter
    # 对高通傅里叶变换 执行傅里叶逆变换，并只取实部
    idft = np.zeros(dftImage.shape[:2], np.float32)  # 快速傅里叶变换的尺寸(优化尺寸)
    cv2.dft(dftHPfilter, idft, cv2.DFT_REAL_OUTPUT + cv2.DFT_INVERSE + cv2.DFT_SCALE)
    # 中心化, centralized 2d array g(x,y) * (-1)^(x+y)
    mask2 = np.ones(dftImage.shape[:2])
    mask2[1::2, ::2] = -1
    mask2[::2, 1::2] = -1
    idftCen = idft * mask2  # g(x,y) * (-1)^(x+y)
    # 截取左上角，大小和输入图像相等
    result = np.clip(idftCen, 0, 255)  # 截断函数，将数值限制在 [0,255]
    imgHPF = result.astype(np.uint8)
    imgHPF = imgHPF[:rows, :cols]

    # # =======高频增强滤波===================
    k1 = 0.5
    k2 = 0.75
    # 在频率域修改傅里叶变换: 傅里叶变换 点乘 低通滤波器
    hpEnhance = np.zeros(dftImage.shape, dftImage.dtype)  # 快速傅里叶变换的尺寸(优化尺寸)
    for j in range(2):
        hpEnhance[:rPadded, :cPadded, j] = dftImage[:rPadded, :cPadded, j] * (k1 + k2*hpFilter)
    # 对高通傅里叶变换 执行傅里叶逆变换，并只取实部
    idft = np.zeros(dftImage.shape[:2], np.float32)  # 快速傅里叶变换的尺寸(优化尺寸)
    cv2.dft(hpEnhance, idft, cv2.DFT_REAL_OUTPUT + cv2.DFT_INVERSE + cv2.DFT_SCALE)
    # 中心化, centralized 2d array g(x,y) * (-1)^(x+y)
    mask2 = np.ones(dftImage.shape[:2])
    mask2[1::2, ::2] = -1
    mask2[::2, 1::2] = -1
    idftCen = idft * mask2  # g(x,y) * (-1)^(x+y)
    # 截取左上角，大小和输入图像相等
    result = np.clip(idftCen, 0, 255)  # 截断函数，将数值限制在 [0,255]
    imgHPE= result.astype(np.uint8)
    imgHPE = imgHPE[:rows, :cols]

    # =======直方图均衡===================
    imgEqu = cv2.equalizeHist(imgHPE)  # 使用 cv2.qualizeHist 完成直方图均衡化变换

    plt.figure(figsize=(9, 6))
    plt.subplot(221), plt.imshow(image, 'gray'), plt.title("Origin"), plt.xticks([]), plt.yticks([])
    plt.subplot(222), plt.imshow(imgHPF, 'gray'), plt.title("Gauss high-pass filter"), plt.xticks([]), plt.yticks([])
    plt.subplot(223), plt.imshow(imgHPE, 'gray'), plt.title("High frequency emphasis"), plt.xticks([]), plt.yticks([])
    plt.subplot(224), plt.imshow(imgEqu, 'gray'), plt.title("Histogram of equalized"), plt.xticks([]), plt.yticks([])
    plt.tight_layout()
    plt.show()

在这里插入图片描述

4.3 频率域拉普拉斯高通滤波

拉普拉斯算子（Laplace）是导数算子，会突出图像中的急剧灰度变化，抑制灰度缓慢变化区域，往往会产生暗色背景下的灰色边缘和不连续图像。将拉普拉斯图像与原图叠加，可以得到保留锐化效果的图像。

拉普拉斯算子（Laplace）是最简单的各向同性导数算子（卷积核）：
$\begin{aligned} \nabla ^2 f &= \dfrac{\partial ^2 f}{\partial x ^2} + \dfrac{\partial ^2 f}{\partial y ^2} \\ \dfrac{\partial ^2 f}{\partial x ^2} &= f(x+1,y) - 2f(x,y) + f(x-1,y) \\ \dfrac{\partial ^2 f}{\partial y ^2} &= f(x,y+1) - 2f(x,y) + f(x,y-1) \\ \nabla ^2 f(x,y) &= f(x+1,y) + f(x-1,y) + f(x,y+1) + f(x,y-1) - 4f(x,y) \end{aligned}$

于是可以得到空间域拉普拉斯核 K1，考虑对角项后可以得到拉普拉斯核 K2。

$\begin{bmatrix} 0 & 1 &0\\ 1 & -4 &1\\ 0 & 1 &0\\ \end{bmatrix}, \ K2= \begin{bmatrix} 1 & 1 &1\\ 1 & -8 &1\\ 1 & 1 &1\\ \end{bmatrix}, \ K3= \begin{bmatrix} 0 & -1 &0\\ -1 & 4 &-1\\ 0 & -1 &0\\ \end{bmatrix}, \ K4= \begin{bmatrix} -1 & -1 &-1\\ -1 & 8 &-1\\ -1 & -1 &-1\\ \end{bmatrix}$

在频率域中拉普拉斯算子可以用传递函数描述：

$\pi ^2 (u^2 + v^2)$
采用 D(u,v) 表示到频率中心的距离函数，可以进一步表示为：
$\pi ^2 [(u-P/2)^2 + (v-Q/2)^2] = - 4 \pi ^2 D^2(u,v)$

于是，

$\begin{align} \nabla ^2 f(x,y) &= J^{-1} \{H(u,v) F(u,v)\}\\ g(x,y) &= f(x,y) + c \nabla ^2 f(x,y) \\ &= J^{-1} \{[1+4 \pi ^2 D^2(u,v)] F(u,v)\} \end{align}$

例程 8.27：频率域图像锐化（拉普拉斯算子）

    # 8.27：频率域图像锐化（拉普拉斯算子）
    def LaplacianFilter(shape):  # 频域 Laplacian 滤波器
        u, v = np.mgrid[-1:1:2.0/shape[0], -1:1:2.0/shape[1]]
        D = np.sqrt(u**2 + v**2)
        kernel = -4 * np.pi**2 * D**2
        return kernel

    def imgHPfilter(image, lpTyper="Laplacian"):  #　频域高通滤波
        # (1) 中心化, centralized 2d array f(x,y) * (-1)^(x+y)
        mask = np.ones(image.shape)
        mask[1::2, ::2] = -1
        mask[::2, 1::2] = -1
        fImage = image * mask  # f(x,y) * (-1)^(x+y)

        # (2) 最优 DFT 扩充尺寸, 快速傅里叶变换的尺寸扩充
        rows, cols = image.shape[:2]  # 原始图片的高度和宽度
        rPadded = cv2.getOptimalDFTSize(rows)  # 最优 DFT 扩充尺寸
        cPadded = cv2.getOptimalDFTSize(cols)  # 用于快速傅里叶变换

        # (3) 边缘扩充(补0), 快速傅里叶变换
        dftImage = np.zeros((rPadded, cPadded, 2), np.float32)  # 对原始图像进行边缘扩充
        dftImage[:rows, :cols, 0] = fImage  # 边缘扩充，下侧和右侧补0
        cv2.dft(dftImage, dftImage, cv2.DFT_COMPLEX_OUTPUT)  # 快速傅里叶变换 (rPad, cPad, 2)

        # (4) 构建 频域滤波器传递函数: 以 Laplacian 为例
        LapFilter = LaplacianFilter((rPadded, cPadded))  # 拉普拉斯滤波器

        # (5) 在频率域修改傅里叶变换: 傅里叶变换 点乘 滤波器传递函数
        dftFilter = np.zeros(dftImage.shape, dftImage.dtype)  # 快速傅里叶变换的尺寸(优化尺寸)
        for j in range(2):
            dftFilter[:rPadded, :cPadded, j] = dftImage[:rPadded, :cPadded, j] * LapFilter

        # (6) 对高通傅里叶变换 执行傅里叶逆变换，并只取实部
        idft = np.zeros(dftImage.shape[:2], np.float32)  # 快速傅里叶变换的尺寸(优化尺寸)
        cv2.dft(dftFilter, idft, cv2.DFT_REAL_OUTPUT + cv2.DFT_INVERSE + cv2.DFT_SCALE)

        # (7) 中心化, centralized 2d array g(x,y) * (-1)^(x+y)
        mask2 = np.ones(dftImage.shape[:2])
        mask2[1::2, ::2] = -1
        mask2[::2, 1::2] = -1
        idftCen = idft * mask2  # g(x,y) * (-1)^(x+y)

        # (8) 截取左上角，大小和输入图像相等
        result = np.clip(idftCen, 0, 255)  # 截断函数，将数值限制在 [0,255]
        imgFilter = result.astype(np.uint8)
        imgFilter = imgFilter[:rows, :cols]
        return imgFilter


    # (1) 读取原始图像
    img = cv2.imread("../images/Fig0338a.tif", flags=0)  # NASA 月球影像图
    rows, cols = img.shape[:2]  # 图片的高度和宽度
    print(rows, cols)

    # (2) 空间域 拉普拉斯算子 (Laplacian)
    # 使用函数 filter2D 实现 Laplace 卷积算子
    kernLaplace = np.array([[0, 1, 0], [1, -4, 1], [0, 1, 0]])  # Laplacian kernel
    imgLaplace1 = cv2.filter2D(img, -1, kernLaplace, borderType=cv2.BORDER_REFLECT)
    # 使用 cv2.Laplacian 实现 Laplace 卷积算子
    imgLaplace2 = cv2.Laplacian(img, -1, ksize=3)
    imgLapReSpace = cv2.add(img, imgLaplace2)  # 恢复原图像

    # (3) 频率域 拉普拉斯算子 (Laplacian)
    imgLaplace = imgHPfilter(img, "Laplacian")  # 调用自定义函数 imgHPfilter()
    imgLapRe = cv2.add(img, imgLaplace)  # 恢复原图像

    plt.figure(figsize=(10, 6))
    plt.subplot(131), plt.imshow(img, 'gray'), plt.title("Origin from NASA"), plt.xticks([]), plt.yticks([])
    plt.subplot(132), plt.imshow(imgLapReSpace, 'gray'), plt.title("Spatial Lapalacian"), plt.xticks([]), plt.yticks([])
    plt.subplot(133), plt.imshow(imgLapRe, 'gray'), plt.title("Freauency Lapalacian"), plt.xticks([]), plt.yticks([])
    plt.tight_layout()
    plt.show()

在这里插入图片描述

5. 选择性滤波

5.1 带阻与带通

空间域和频率域线性滤波器可以分为四类：低通滤波器、高通滤波器、带通滤波器和带阻滤波器。高通滤波和低通滤波都是在整个频率矩形上操作，带通滤波和带阻滤波则是对特定频带处理，属于选择性滤波。

频率域的高通滤波器可以由低通滤波器推导而来。类似地，频率域中的带通和带阻滤波器的传递函数，可以通过低通滤波器和高通滤波器的组合来构建。

理想带阻滤波器（IBRF）的传递函数为：
$H(u,v)=\begin{cases} 0,\ (C_0-W/2) \leq D(u,v) \leq (C_0+W/2)\\ 1,\ else \end{cases}$
高斯带阻滤波器（GBRF）的传递函数为：
$H(u,v)=1-e^{-[ \frac {D^2(u,v) - C_0^2} {D(u,v)W}]^2}$
**巴特沃斯带阻滤波器（BBRF）的传递函数为： **
$\frac {1} {1 +[ \frac {D(u,v)W} {D^2(u,v) - C_0^2}]^{2n}}$

例程 8.28 带阻滤波器的传递函数

    # 例程 8.28 带阻滤波器的传递函数
    def ideaBondResistFilter(shape, radius=10, w=5):  # 理想带阻滤波器
        u, v = np.meshgrid(np.arange(shape[1]), np.arange(shape[0]))
        D = np.sqrt((u - shape[1]//2)**2 + (v - shape[0]//2)**2)
        D0 = radius
        halfW = w/2
        kernel = np.piecewise(D, [D<=D0+halfW, D<=D0-halfW], [1, 0])
        kernel = 1 - kernel  # 带阻
        return kernel

    def gaussBondResistFilter(shape, radius=10, w=5):  # 高斯带阻滤波器
        # 高斯滤波器：# Gauss = 1/(2*pi*s2) * exp(-(x**2+y**2)/(2*s2))
        u, v = np.meshgrid(np.arange(shape[1]), np.arange(shape[0]))
        D = np.sqrt((u - shape[1]//2)**2 + (v - shape[0]//2)**2)
        C0 = radius
        kernel = 1 - np.exp(-(D-C0)**2 / (w**2))
        return kernel

    def butterworthBondResistFilter(shape, radius=10, w=5, n=1):  # 巴特沃斯带阻滤波
        u, v = np.meshgrid(np.arange(shape[1]), np.arange(shape[0]))
        D = np.sqrt((u - shape[1]//2)**2 + (v - shape[0]//2)**2)
        C0 = radius
        epsilon = 1e-8  # 防止被 0 除
        kernel = 1.0 / (1.0 + np.power(D*w/(D**2-C0**2+epsilon), 2*n))
        return kernel

    # 理想、高斯、巴特沃斯带阻滤波器传递函数
    shape = [128, 128]
    radius = 32
    IBRF = ideaBondResistFilter(shape, radius=radius)
    GBRF = gaussBondResistFilter(shape, radius=radius)
    BBRF = butterworthBondResistFilter(shape, radius=radius)

    filters = ["IBRF", "GBRF", "BBRF"]
    u, v = np.mgrid[-1:1:2.0/shape[0], -1:1:2.0/shape[1]]
    fig = plt.figure(figsize=(10, 8))
    for i in range(3):
        hpFilter = eval(filters[i]).copy()

        ax1 = fig.add_subplot(3, 3, 3*i+1)
        ax1.imshow(hpFilter, 'gray')
        ax1.set_title(filters[i]), ax1.set_xticks([]), ax1.set_yticks([])

        ax2 = plt.subplot(3,3,3*i+2, projection='3d')
        ax2.set_title("transfer function")
        ax2.plot_wireframe(u, v, hpFilter, rstride=2, linewidth=0.5, color='c')
        ax2.set_xticks([]), ax2.set_yticks([]), ax2.set_zticks([])

        ax3 = plt.subplot(3,3,3*i+3)
        profile = hpFilter[shape[0]//2:, shape[1]//2]
        ax3.plot(profile), ax3.set_title("profile"), ax3.set_xticks([]), ax3.set_yticks([])

    plt.show()

在这里插入图片描述

5.2 陷波滤波器（Notch Filter）

陷波滤波器阻止或通过预定的频率矩形邻域中的频率，是重要的选择性滤波器。

陷波滤波器可以在某一个频率点迅速衰减输入信号，以达到阻碍此频率信号通过的滤波效果的滤波器。陷波滤波器属于带阻滤波器的一种，它的阻带非常狭窄，起阶数必须是二阶（含二阶）以上。

零相移滤波器必须给予原点（频率矩形中心）对称，因此中心为 $u_0,v_0)$ 的陷波滤波器传递函数在 $u_0, -v_0)$ 位置必须有一个对应的陷波。

陷波带阻滤波器的传递函数可以用中心已被平移到陷波滤波器中心的高通滤波器的乘积来构造：
$H_{NR}(u,v) = \prod_{k=1}^Q H_k(u,v) H_{-k}(u,v)$

其中，滤波器的距离计算公式为：
$D_k(u,v) = \sqrt{(u-M/2-u_k)^2 + (v-N/2-v_k)^2} \\ D_{-k}(u,v) = \sqrt{(u-M/2+u_k)^2 + (v-N/2+v_k)^2}$
例如，具有 3个陷波对的 n 阶巴特沃斯陷波带阻滤波器为：
$H_{NR}(u,v) = \prod_{k=1}^3 [\frac{1}{1+[D_{0k}/D_k(u,v)]^n}] [\frac{1}{1+[D_{-k}/D_k(u,v)]^n}]$
陷波带通滤波器的传递函数可用陷波带阻滤波器得到：
$H_{NP}(u,v) = 1 - H_{NR}(u,v)$

例程 8.29 使用陷波滤波删除数字化印刷图像中的莫尔模式

本例使用陷波滤波降低数字化印刷图像中的莫尔波纹。

摩尔纹是差拍原理的表现。两个频率接近的等幅正弦波叠加后，信号幅度将按照两个频率之差变化。如果感光元件像素的空间频率与影像中条纹的空间频率接近，就会产生摩尔纹。

陷波滤波是选择性的修改DFT的局部区域。典型的处理方法是交互式操作，直接用鼠标在傅里叶频谱中选择矩形区域，找出最大值点作为 (uk,vk)。为了简化程序，本例程删除了鼠标交互部分，只保留了给出 (uk,vk) 后的滤波处理过程。

    # 例程 8.29 使用陷波滤波删除数字化印刷图像中的莫尔模式

    def gaussLowPassFilter(shape, radius=10):  # 高斯低通滤波器
        u, v = np.mgrid[-1:1:2.0 / shape[0], -1:1:2.0 / shape[1]]
        D = np.sqrt(u ** 2 + v ** 2)
        D0 = radius / shape[0]
        kernel = np.exp(- (D ** 2) / (2 * D0 ** 2))
        return kernel

    def butterworthNRFilter(shape, radius=9, uk=60, vk=80, n=2):  # 巴特沃斯陷波带阻滤波器
        M, N = shape[1], shape[0]
        u, v = np.meshgrid(np.arange(M), np.arange(N))
        Dm = np.sqrt((u - M // 2 - uk) ** 2 + (v - N // 2 - vk) ** 2)
        Dp = np.sqrt((u - M // 2 + uk) ** 2 + (v - N // 2 + vk) ** 2)
        D0 = radius
        n2 = n * 2
        kernel = (1 / (1 + (D0 / (Dm + 1e-6)) ** n2)) * (1 / (1 + (D0 / (Dp + 1e-6)) ** n2))
        return kernel

    def imgFrequencyFilter(img, lpTyper="GaussLP", radius=10):
        normalize = lambda x: (x - x.min()) / (x.max() - x.min() + 1e-8)

        # (1) 边缘填充
        imgPad = np.pad(img, ((0, img.shape[0]), (0, img.shape[1])), mode="reflect")
        # (2) 中心化: f(x,y) * (-1)^(x+y)
        mask = np.ones(imgPad.shape)
        mask[1::2, ::2] = -1
        mask[::2, 1::2] = -1
        imgPadCen = imgPad * mask
        # (3) 傅里叶变换
        fft = np.fft.fft2(imgPadCen)

        # (4) 构建 频域滤波器传递函数:
        if lpTyper == "GaussLP":
            print(lpTyper)
            freFilter = gaussLowPassFilter(imgPad.shape, radius=60)
        elif lpTyper == "GaussHP":
            freFilter = gaussLowPassFilter(imgPad.shape, radius=60)
        elif lpTyper == "ButterworthNR":
            print(lpTyper)
            freFilter = butterworthNRFilter(imgPad.shape, radius=9, uk=60, vk=80, n=2)  # 巴特沃斯陷波带阻滤波器
        elif lpTyper == "MButterworthNR":
            print(lpTyper)
            BNRF1 = butterworthNRFilter(imgPad.shape, radius=9, uk=60, vk=80, n=2)  # 巴特沃斯陷波带阻滤波器
            BNRF2 = butterworthNRFilter(imgPad.shape, radius=9, uk=-60, vk=80, n=2)
            BNRF3 = butterworthNRFilter(imgPad.shape, radius=9, uk=60, vk=160, n=2)
            BNRF4 = butterworthNRFilter(imgPad.shape, radius=9, uk=-60, vk=160, n=2)
            freFilter = BNRF1 * BNRF2 * BNRF3 * BNRF4
        else:
            print("Error of unknown filter")
            return -1

        # (5) 在频率域修改傅里叶变换: 傅里叶变换 点乘 滤波器传递函数
        freTrans = fft * freFilter
        # (6) 傅里叶反变换
        ifft = np.fft.ifft2(freTrans)
        # (7) 去中心化反变换的图像
        M, N = img.shape[:2]
        mask2 = np.ones(imgPad.shape)
        mask2[1::2, ::2] = -1
        mask2[::2, 1::2] = -1
        ifftCenPad = ifft.real * mask2
        # (8) 截取左上角，大小和输入图像相等
        imgFilter = ifftCenPad[:M, :N]
        imgFilter = np.clip(imgFilter, 0, imgFilter.max())
        imgFilter = np.uint8(normalize(imgFilter) * 255)
        return imgFilter


    # 使用陷波滤波删除数字化印刷图像中的莫尔模式
    # (1) 读取原始图像
    img = cv2.imread("../images/Fig0464a.tif", flags=0)  # flags=0 读取为灰度图像
    fig = plt.figure(figsize=(10, 5))
    plt.subplot(141), plt.title("Original"), plt.axis('off'), plt.imshow(img, cmap='gray')

    # (2) 图像高斯低通滤波
    imgGLPF = imgFrequencyFilter(img, lpTyper="GaussLP", radius=30)  # 图像高斯低通滤波
    plt.subplot(142), plt.title("GaussLP filter"), plt.axis('off'), plt.imshow(imgGLPF, cmap='gray')

    # (3) 图像巴特沃斯陷波带阻滤波
    imgBNRF = imgFrequencyFilter(img, lpTyper="ButterworthNR", radius=9)
    plt.subplot(143), plt.title("ButterworthNR filter"), plt.axis('off'), plt.imshow(imgBNRF, cmap='gray')

    # (4) 叠加巴特沃斯陷波带阻滤波
    imgSBNRF = imgFrequencyFilter(img, lpTyper="MButterworthNR", radius=9)
    plt.subplot(144), plt.title("Superimposed BNRF"), plt.axis('off'), plt.imshow(imgSBNRF, cmap='gray')

    plt.tight_layout()
    plt.show()

在这里插入图片描述

注：本文中部分原始图片来自 Rafael C. Gonzalez “Digital Image Processing, 4th.Ed.”，特此致谢。

youcans@xupt 原创作品，转载必须标注原文链接：(https://blog.csdn.net/youcans/article/details/122991620)
Copyright 2022 youcans, XUPT
欢迎关注『youcans 的 OpenCV 学习课』系列，持续更新
【youcans 的 OpenCV 学习课】1. 安装与环境配置
 【youcans 的 OpenCV 学习课】2. 图像读取与显示
 【youcans 的 OpenCV 学习课】3. 图像的创建与修改
 【youcans 的 OpenCV 学习课】4. 图像的叠加与混合
 【youcans 的 OpenCV 学习课】5. 图像的几何变换
 【youcans 的 OpenCV 学习课】6. 灰度变换与直方图处理
 【youcans 的 OpenCV 学习课】7. 空间域图像滤波
 【youcans 的 OpenCV 学习课】8. 频率域图像滤波（上）
【youcans 的 OpenCV 学习课】9. 频率域图像滤波（下）

youcans_

关注

0
点赞
踩
9

收藏

觉得还不错? 一键收藏
打赏
0
评论
【youcans 的图像处理学习课】9. 频率域图像滤波（下）

频率域滤波是滤波器传递函数与输入图像的傅里叶变换的对应像素相乘。频率域中的滤波概念更加直观，滤波器设计也更容易。本文提供上述各种算法的完整例程和运行结果，并通过一个应用案例示范综合使用多种图像增强方法。1. 安装与环境配置2. 图像读取与显示3. 图像的创建与修改4. 图像的叠加与混合5. 图像的几何变换6. 灰度变换与直方图处理7. 空间域图像滤波8. 频率域图像滤波...
复制链接

扫一扫