OpenCV 完整例程90. 频率域陷波滤波器

Posted 2022-02-17 小白YouCans

【OpenCV 完整例程】90. 频率域陷波滤波器

欢迎关注 『OpenCV 完整例程 100 篇』 系列，持续更新中
欢迎关注 『Python小白的OpenCV学习课』 系列，持续更新中

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)$$

其中，滤波器的距离计算公式为：
$$\begin{gathered} 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} \end{gathered}$$
于是，3 阶巴特沃斯陷波带阻滤波器为：
$$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) 后的滤波处理过程。

# OpenCVdemo08.py
# Demo08 of OpenCV
# 8. 图像的频率域滤波
# Copyright 2021 Youcans, XUPT
# Crated：2021-12-30  
   
    # 例程 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")以上是关于OpenCV 完整例程90. 频率域陷波滤波器的主要内容，如果未能解决你的问题，请参考以下文章
 
 youcans 的 OpenCV 例程 200 篇103. 陷波带阻滤波器消除周期噪声干扰
 OpenCV 完整例程79. 频率域图像滤波的基本步骤
 OpenCV 完整例程81. 频率域高斯低通滤波器
 OpenCV 完整例程85. 频率域高通滤波器的应用
 OpenCV 完整例程82. 频率域巴特沃斯低通滤波器
 OpenCV 完整例程80. 频率域图像滤波详细步骤