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

已于 2023-05-18 19:24:19 修改
阅读量1.3k 收藏 11
点赞数
分类专栏： # 数字图像处理 OpenCV 文章标签： opencv python 计算机视觉 图像处理 算法
于 2022-02-19 07:30:00 首次发布
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/youcans/article/details/122991620
版权

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


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

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


文章目录

3. 频率域低通滤波器

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

3.1 频率域图像滤波基础

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

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

空间域滤波器和频率域滤波器也是相互对应的，二者形成一个傅里叶变换对：
f ( x , y ) ⊗ h ( x , y ) ⇔ F ( u , v ) H ( u , v ) f ( x , y ) h ( x , y ) ⇔ F ( u , v ) ⊗ H ( u , v ) f(x,y) \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) f(x,y)⊗h(x,y)⇔F(u,v)H(u,v)f(x,y)h(x,y)⇔F(u,v)⊗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(x,y) f(x,y) 经傅里叶变换为 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) 的频谱成分进行修改，最后通过傅里叶逆变换（IDFT）返回空间域，得到增强的图像 g ( x , y ) g(x,y) g(x,y)。

f ( x , y ) → D F T F ( u , v ) → 滤波 H ( u , v ) G ( u , v ) → I D F T g ( x , y ) f(x,y) \xrightarrow {DFT} F(u,v) \xrightarrow [滤波] {H(u,v)} G(u,v) \xrightarrow {IDFT} g(x,y) f(x,y)DFT ​F(u,v)H(u,v) 滤波​G(u,v)IDFT ​g(x,y)

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

（1）对原始图像 f ( x , y ) f(x,y) f(x,y) 进行傅里叶变换，得到 F ( u , v ) F(u,v) F(u,v)；
（2）将图像的傅里叶变换 F ( u , v ) F(u,v) F(u,v) 与传递函数 H ( u , v ) H(u,v) H(u,v) 进行卷积运算，得到滤波后的频谱 G ( u , v ) G(u,v) G(u,v)；
（3）对 G ( u , v ) G(u,v) G(u,v) 进行傅里叶逆变换，得到增强图像 g ( x , y ) g(x,y) 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 以理想低通滤波器为例，给出了频率域图像滤波的一般步骤。

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

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

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

利用二维高斯蒙板，可以滤去高频、保留低频。频域高斯函数的方差 σ 2 \sigma^2 σ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 D_0 D0​ 处的 n 阶巴特沃斯（Butterworth）低通滤波器的传递函数为：
H ( u , v ) = 1 1 + [ D ( u , v ) / D 0 ] 2 n H(u,v) = \frac {1}{1+[D(u,v) / D_0]^{2n}} H(u,v)=1+[D(u,v)/D0​]2n1​

当 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 D_0 D0​ 的高斯低通滤波器简单处理后的修复效果。当然这只是修复的一个步骤，接下来还要进行其它处理如阈值处理、细化处理，这些方法将在后续内容讨论。

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

    # 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 H P ( u , v ) = 1 − H L P ( u , v ) H_{HP}(u,v) = 1- H_{LP}(u,v) HHP​(u,v)=1−HLP​(u,v)
式中， H H P ( u , v ) H_{HP}(u,v) HHP​(u,v)、 H L P ( u , v ) H_{LP}(u,v) HLP​(u,v) 分别表示高通滤波器、低通滤波器的传递函数。

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


例程 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 L P ( x , y ) f_{LP}(x,y) fLP​(x,y) 表示低通滤波的平滑图像，则：
g m a s k ( x , y ) = f ( x , y ) − f L P ( x , y ) g ( x , y ) = f ( x , y ) + k ∗ g m a s k ( x , y ) ， k > 0 \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} gmask​(x,y)g(x,y)​=f(x,y)−fLP​(x,y)=f(x,y)+k∗gmask​(x,y)，k>0​​
当 k>1 时，实现高提升滤波；当 k=1 时，实现钝化掩蔽；k<1时，可以减弱钝化掩蔽的强度。

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

在频率域实现钝化掩蔽，高频强调滤波器传递函数为：
g ( x , y ) = J − 1 { [ 1 + k H H P ( u , v ) ] F ( u , v ) } g(x,y)= J^{-1} \{ [1+k H_{HP}(u,v)]F(u,v) \} g(x,y)=J−1{[1+kHHP​(u,v)]F(u,v)}
高频强调滤波的通用公式是：
g ( x , y ) = J − 1 { [ k 1 + k 2 H H P ( u , v ) ] F ( u , v ) } g(x,y)= J^{-1} \{ [k_1 + k_2 H_{HP}(u,v)]F(u,v) \} g(x,y)=J−1{[k1​+k2​HHP​(u,v)]F(u,v)}
式中， k 1 ≥ 0 k_1 \geq 0 k1​≥0 偏移传递函数的值， k 2 ≥ 0 k_2 \geq 0 k2​≥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）是最简单的各向同性导数算子（卷积核）：
∇ 2 f = ∂ 2 f ∂ x 2 + ∂ 2 f ∂ y 2 ∂ 2 f ∂ x 2 = f ( x + 1 , y ) − 2 f ( x , y ) + f ( x − 1 , y ) ∂ 2 f ∂ y 2 = f ( x , y + 1 ) − 2 f ( x , y ) + f ( x , y − 1 ) ∇ 2 f ( x , y ) = f ( x + 1 , y ) + f ( x − 1 , y ) + f ( x , y + 1 ) + f ( x , y − 1 ) − 4 f ( x , y ) \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} ∇2f∂x2∂2f​∂y2∂2f​∇2f(x,y)​=∂x2∂2f​+∂y2∂2f​=f(x+1,y)−2f(x,y)+f(x−1,y)=f(x,y+1)−2f(x,y)+f(x,y−1)=f(x+1,y)+f(x−1,y)+f(x,y+1)+f(x,y−1)−4f(x,y)​

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

K 1 = [ 0 1 0 1 − 4 1 0 1 0 ] ,   K 2 = [ 1 1 1 1 − 8 1 1 1 1 ] ,   K 3 = [ 0 − 1 0 − 1 4 − 1 0 − 1 0 ] ,   K 4 = [ − 1 − 1 − 1 − 1 8 − 1 − 1 − 1 − 1 ] K1= \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} K1= ​010​1−41​010​ ​, K2= ​111​1−81​111​ ​, K3= ​0−10​−14−1​0−10​ ​, K4= ​−1−1−1​−18−1​−1−1−1​ ​

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

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

于是，

∇ 2 f ( x , y ) = J − 1 { H ( u , v ) F ( u , v ) } g ( x , y ) = f ( x , y ) + c ∇ 2 f ( x , y ) = J − 1 { [ 1 + 4 π 2 D 2 ( u , v ) ] F ( 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} ∇2f(x,y)g(x,y)​=J−1{H(u,v)F(u,v)}=f(x,y)+c∇2f(x,y)=J−1{[1+4π2D2(u,v)]F(u,v)}​​

例程 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 ) = { 0 ,   ( C 0 − W / 2 ) ≤ D ( u , v ) ≤ ( C 0 + W / 2 ) 1 ,   e l s e H(u,v)=\begin{cases} 0,\ (C_0-W/2) \leq D(u,v) \leq (C_0+W/2)\\ 1,\ else \end{cases} H(u,v)={0, (C0​−W/2)≤D(u,v)≤(C0​+W/2)1, else​
高斯带阻滤波器（GBRF） 的传递函数为：
H ( u , v ) = 1 − e − [ D 2 ( u , v ) − C 0 2 D ( u , v ) W ] 2 H(u,v)=1-e^{-[ \frac {D^2(u,v) - C_0^2} {D(u,v)W}]^2} H(u,v)=1−e−[D(u,v)WD2(u,v)−C02​​]2
**巴特沃斯带阻滤波器（BBRF） 的传递函数为： **
H ( u , v ) = 1 1 + [ D ( u , v ) W D 2 ( u , v ) − C 0 2 ] 2 n H(u,v)= \frac {1} {1 +[ \frac {D(u,v)W} {D^2(u,v) - C_0^2}]^{2n}} H(u,v)=1+[D2(u,v)−C02​D(u,v)W​]2n1​


例程 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) (u0​,v0​) 的陷波滤波器传递函数在 ( − u 0 , − v 0 ) (-u_0, -v_0) (−u0​,−v0​) 位置必须有一个对应的陷波。

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

其中，滤波器的距离计算公式为：
D k ( u , v ) = ( u − M / 2 − u k ) 2 + ( v − N / 2 − v k ) 2 D − k ( u , v ) = ( 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} \\ D_{-k}(u,v) = \sqrt{(u-M/2+u_k)^2 + (v-N/2+v_k)^2} Dk​(u,v)=(u−M/2−uk​)2+(v−N/2−vk​)2 ​D−k​(u,v)=(u−M/2+uk​)2+(v−N/2+vk​)2 ​
例如，具有 3个陷波对的 n 阶巴特沃斯陷波带阻滤波器为：
H N R ( u , v ) = ∏ k = 1 3 [ 1 1 + [ D 0 k / D k ( u , v ) ] n ] [ 1 1 + [ D − k / D k ( u , v ) ] 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}] HNR​(u,v)=k=1∏3​[1+[D0k​/Dk​(u,v)]n1​][1+[D−k​/Dk​(u,v)]n1​]
陷波带通滤波器的传递函数可用陷波带阻滤波器得到：
H N P ( u , v ) = 1 − H N R ( u , v ) H_{NP}(u,v) = 1 - H_{NR}(u,v) HNP​(u,v)=1−HNR​(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. 频率域图像滤波（下）

opencv 写的频率域滤波（高斯，理想，巴特沃斯）
07-28
opencv 写的频率域滤波包含高斯，理想，巴特沃斯低通滤波，
数字图像处理 A_图像频率域滤波.pdf
01-29
3. 掌握MATLAB中与图像频率域处理相关的函数，例如fft2、ifft2、fftshift、ifftshift等函数的使用方法和它们在处理过程中的作用。这包括函数的输入输出数据类型和结构的理解，以及函数内部执行的数学运算。 4. 通过...
OpenCV —— 频率域滤波（傅里叶变换，低通和高通滤波，带通和带阻滤波，同态滤波）
热门推荐
学无止境
03-01 1万+
频率域滤波基本概念傅里叶变换二维离散的傅里叶变换快速傅里叶变换傅里叶幅度谱与相位谱谱残差显著性检测卷积与傅里叶变换的频率域滤波低通滤波和高通滤波带通和带阻滤波同态滤波 基本概念 频率域滤波 —— 百度百科 频率域滤波是对图像进行傅里叶变换，将图像由图像空间转换到频域空间，然后在频率域中对图像的频谱作分析处理，以改变图像的频率特征。 滤波: 狭义地说，滤波是指改变信号中各个频率分量的相对大小、或者分离出来加以抑制、甚至全部滤除某些频率分量的过程。广义地说，滤波是把某种信号处理成为另一种信号的过程。空间域滤波
【OpenCV】OpenCV基础教程（10）—— 频率域滤波
Quentin的博客
11-10 3456
10. 频率域滤波 10.1 概述及原理详解 频率域滤波的步骤如下： 10.2 低通滤波 低通滤波器越靠近中心点位置的值越接近于 1，越远离中心位置的值就越小于 1，与傅里叶变换相乘后，相当于保留了低频信息，消弱或者移除了高频信息。 10.2.1 理想低通滤波器 假设 𝐻、𝑊 分别代表图像快速傅里叶变换的高、宽。傅里叶谱的最大值（中心点）的位置在 (maxR, maxC)。radius 代表截断频率。𝐷(𝑟, 𝑐) 代表 (𝑟, 𝑐) 到中心位置的距离。 ilpFilter(r,c)={1D(r,c)≤r
opencv图像分析与处理（11）- 频率域滤波消除周期噪声
kevin_zhao_zl的博客
06-20 5945
本节为opencv数字图像处理（10）：图像复原与重建的第二小节，频率域滤波消除周期噪声，主要包括：带阻滤波器、带通滤波器、陷波滤波器以及最佳陷波滤波器的数学推导。
【OpenCV】Chapter6.频率域图像滤波
兴趣使然的创作者
08-29 1993
最近想对OpenCV进行系统学习，看到网上这份教程写得不错，于是跟着来学习实践一下。
paddedsize.rar_cardx9z_greatestbj9_paddedsize_paddedsize.m_频率域滤波
07-15
在IT领域，尤其是在信号处理和图像处理中，频率域滤波是一种常见的技术，它通过将信号或图像转换到频域来分析和修改其特性。在这个上下文中，"paddedsize.rar_cardx9z_greatestbj9_paddedsize_paddedsize.m_频率域...
image-process.rar_频率域滤波
09-19
在图像处理领域，频率域滤波是一种常见的图像增强技术，尤其在处理噪声和改善图像质量方面有着重要作用。本文将深入探讨这一主题，并结合提供的“image-process.rar”压缩包中的资源，如“图像增强滤波的PPT及代码....
数字图像处理实验报告-数字图像空间与频率滤波.docx
01-15
在数字图像处理领域，滤波是一种常见的操作，用于去除噪声、增强图像特征或者改变图像的局部特性。本实验报告主要探讨了数字图像的空间滤波和频率滤波，使用MATLAB软件进行实现。实验目的是让学生熟悉MATLAB环境，...
matlab数字图像处理课程作业-频率域滤波matlab源码.zip
12-28
matlab数字图像处理课程作业_频率域滤波matlab源码.zipmatlab数字图像处理课程作业_频率域滤波matlab源码.zipmatlab数字图像处理课程作业_频率域滤波matlab源码.zipmatlab数字图像处理课程作业_频率域滤波matlab源码...
OpenCV写的频率域低通滤波器
11-20
OpenCV写的频率率低通滤波器。
OpenCV笔记21 频率域滤波
weixin_39354845的博客
02-07 2169
一、概述 图像的傅里叶变换及其两个重要的度量：幅度谱和相位谱。了解两个重要的概念：低频和高频。低频指的是图 的傅里叶变换“中心位置”附近的区域。注意，如无特殊说明，后面所提到的图像的傅里叶变换都是中心化后的。高频随着到“中心位置”距离的增加而增加，即傅里叶变换中心位置的外围区域，这里的“中心位置”指的是傅里叶变换所对应的幅度谱最大值的位置。 频率域滤波器在程序或者数学运算中的呈现可以理解为一个矩阵，该矩阵的宽、高和图像的傅里叶变换的宽、高是相同的，下面所涉及的常用的低通、高通、带通、带阻 等滤...
opencv图像分析与处理（4）- 频率域滤波的基础概念
kevin_zhao_zl的博客
06-12 554
本节为opencv数字图像处理（4）：频率域滤波的第一小节，频率域滤波的基础概念，主要包括：复数定义的回顾、傅里叶级数的形式、冲激及其取样特性以及连续函数的傅里叶变换。
OpenCV实现图像基础频率域滤波
小学僧的博客
03-02 3703
频率域滤波
【OpenCV 例程200篇】78. 频率域图像滤波基础
youcans的博客
01-24 2749
空间域滤波器和频率域滤波器实际上是相互对应的，有些空间域滤波器在频率域通过傅里叶变换实现会更方便、更快速。 对信号或图像进行傅里叶变换后，可以得到信号或图像的低频信息和高频信息。低频信息对应图像中缓慢变化的灰度分量，高频信息则对应尖锐变化的灰度分量。 低频滤波器本质上就是构造一个矩阵，越靠近中心的位置越接近于 1，而远离中心位置的值则接近于 0。简单地，生成一个矩形窗口遮罩，在黑色（置 0）遮罩图像的中心开有白色（置 1）窗口，就得到一个低通滤波器。
OpenCV 笔记(34)：频域低通滤波——理想低通滤波器
Java与Android技术栈
06-25 1732
1. 频域低通滤波1.1 频域滤波频域滤波器是一种通过对图像进行傅里叶变换，然后在频域上选择性地滤除或增强某些频率成分来实现图像滤波的方法。频域滤波和空域滤波都有各自的优缺点和使用场景：特性频域滤波空域滤波原理基于频谱进行操作基于像素进行操作实现计算量较大计算量较小滤波效果灵活可控较简单应用去周期性噪声、增强细节、图像配准去随机噪声、平滑图像、锐化图像频域滤波器的基本步骤如下：对图像进行傅里叶变...
OpenCV 低通滤波
weixin_30467087的博客
08-22 400
#include<opencv2/core/core.hpp> #include<opencv2/highgui/highgui.hpp> #include<opencv2/imgproc/imgproc.hpp> using namespace std; using namespace cv; Mat I;//输入的图像矩阵 M...
opencv3.0 图像滤波方式代码实现
一只在拱菜的猪
03-13 1424
#include #include #include using namespace std; int main() { string imName = "lena.jpg"; cv::Mat im = cv::imread(imName); cv::imshow("原图",im); cv::Mat imGs = im.clone(); int kSize = 21; cv::Ma
c++ opencv数字图像处理：频率域滤波--低通滤波--巴特沃斯低通滤波
最新发布
生活需要深度
08-12 308
文章目录 前言 一、巴特沃斯低通滤波器（BLPF） 二、代码 三、说明 前言 数字图像处理c++ opencv（VS2019 opencv4.53）持续更新一、巴特沃斯低通滤波器（BLPF）D 2 ( u , v ) 为 距 离 中 心 点 的 距 离 平 方 ， D 0 为 设 置 半 径 ， 通 过 n 的 值 来 调 整 滤 波 器 D^2(u,v)为距离中心点的距离平方，D_0为设置半径，通过n的值来调整滤波器D 2 (u,v)为距离中心点的距离平方，D 0 ​ 为设置半径，通过n的值来调整
北京大学数字图像处理课件详解：频率域增强与滤波
本课件是关于数字图像处理的详细介绍，由彭宇新教授提供，适用于研究生级别的学习。内容涵盖了800多页的深度讲解，主要聚焦于频率域图像增强技术。频率域在图像处理中的重要性在于： 1. 频率域滤波的优势：在频率域...
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

youcans_

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付：¥1
获取中
扫码支付

您的余额不足，请更换扫码支付或充值

打赏作者

实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值