90. 频率域陷波滤波器

6. 频率域图像滤波索引

一、陷波滤波器 (Notch Filter)

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

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

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

  • 陷波带阻滤波器的传递函数可以用中心已被平移到陷波滤波器中心的高通滤波器的乘积来构造 :



    其中,滤波器的距离计算公式为:



    于是,3 阶巴特沃斯陷波带阻滤波器为:
  • 陷波带通滤波器的传递函数可用陷波带阻滤波器得到:


二、例程

  • 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()

三、资料

youcans_的博客:
https://blog.csdn.net/youcans/article/details/122781106
©著作权归作者所有,转载或内容合作请联系作者
平台声明:文章内容(如有图片或视频亦包括在内)由作者上传并发布,文章内容仅代表作者本人观点,简书系信息发布平台,仅提供信息存储服务。

推荐阅读更多精彩内容