失眠网 > 【OpenCV 例程 300 篇】103. 陷波带阻滤波器消除周期噪声干扰

【OpenCV 例程 300 篇】103. 陷波带阻滤波器消除周期噪声干扰

时间：2022-02-10 11:27:40

专栏地址：『youcans 的 OpenCV 例程 300篇 - 总目录』

【第 7 章：图像复原与重建】
102. 陷波带阻滤波器的传递函数
103. 陷波带阻滤波器消除周期噪声干扰

【youcans 的 OpenCV 例程 200 篇】103. 陷波带阻滤波器消除周期噪声干扰

4.1 陷波滤波器（Notch Filter）

陷波滤波器阻止或通过预定的频率矩形邻域中的频率，可以很好地复原被周期性噪声干扰的图像。

陷波滤波器可以在某一个频率点迅速衰减输入信号，以达到阻碍此频率信号通过的滤波效果的滤波器。

陷波带阻滤波器的传递函数是中心平移到陷波中心的各个高通滤波器的乘积：

HNR(u,v)=∏k=1QHk(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∏QHk(u,v)H−k(u,v)

其中，滤波器的距离计算公式为：

Dk(u,v)=(u−M/2−uk)2+(v−N/2−vk)2D−k(u,v)=(u−M/2+uk)2+(v−N/2+vk)2D_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)2D−k(u,v)=(u−M/2+uk)2+(v−N/2+vk)2

例如，具有 3个陷波对的 n 阶巴特沃斯陷波带阻滤波器为：

HNR(u,v)=∏k=13[11+[D0k/Dk(u,v)]n][11+[D0k/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_{0k}/D_{-k}(u,v)]^n}] HNR(u,v)=k=1∏3[1+[D0k/Dk(u,v)]n1][1+[D0k/D−k(u,v)]n1]

例程 9.17：陷波带阻滤波器消除周期噪声干扰

# 9.17: 陷波带阻滤波器消除周期噪声干扰def butterworthNRFilter(img, radius=10, uk=10, vk=10, n=2): # 巴特沃斯陷波带阻滤波器M, N = img.shape[1], img.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 = radiusn2 = 2 * nkernel = (1 / (1 + (D0 / (Dm + 1e-6))**n2)) * (1 / (1 + (D0 / (Dp + 1e-6))**n2))return kernel# (1) 读取原始图像img = cv2.imread("../images/Fig0505a.tif", flags=0) # flags=0 读取为灰度图像imgFloat32 = np.float32(img) # 将图像转换成 float32rows, cols = img.shape[:2] # 图片的高度和宽度fig = plt.figure(figsize=(9, 6))plt.subplot(231), plt.title("Original image"), plt.axis('off'), plt.imshow(img, cmap='gray')# (2) 中心化, centralized 2d array f(x,y) * (-1)^(x+y)mask = np.ones(img.shape)mask[1::2, ::2] = -1mask[::2, 1::2] = -1fImage = imgFloat32 * mask # f(x,y) * (-1)^(x+y)# (3) 快速傅里叶变换rPadded = cv2.getOptimalDFTSize(rows) # 最优 DFT 扩充尺寸cPadded = cv2.getOptimalDFTSize(cols) # 用于快速傅里叶变换dftImage = np.zeros((rPadded, cPadded, 2), np.float32) # 对原始图像进行边缘扩充dftImage[:rows, :cols, 0] = fImage # 边缘扩充，下侧和右侧补0cv2.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(232), plt.axis('off'), plt.title("DFT spectrum")plt.imshow(dftAmpNorm, cmap='gray')plt.arrow(445, 370, 25, 30, width=5, length_includes_head=True, shape='full') # 在图像上加上箭头plt.arrow(550, 490, -25, -30, width=5, length_includes_head=True, shape='full') # 在图像上加上箭头# (4) 构建陷波带阻滤波器传递函数BRFilter = butterworthNRFilter(dftImage, radius=15, uk=25, vk=16, n=3) # 巴特沃斯陷波带阻滤波器, 处理周期噪声plt.subplot(233), plt.axis('off'), plt.title("Butterworth notch resist filter")plt.imshow(BRFilter, cmap='gray')# (5) 在频率域修改傅里叶变换: 傅里叶变换点乘陷波带阻滤波器dftFilter = np.zeros(dftImage.shape, dftImage.dtype) # 快速傅里叶变换的尺寸(优化尺寸)for i in range(2):dftFilter[:rPadded, :cPadded, i] = dftImage[:rPadded, :cPadded, i] * BRFilter# 频域滤波傅里叶变换的傅里叶谱nrfDftAmp = cv2.magnitude(dftFilter[:, :, 0], dftFilter[:, :, 1]) # 傅里叶变换的幅度谱nrfDftAmpLog = np.log(1.0 + nrfDftAmp) # 幅度谱对数变换，以便于显示nrfDftAmpNorm = np.uint8(cv2.normalize(nrfDftAmpLog, None, 0, 255, cv2.NORM_MINMAX)) # 归一化为 [0,255]plt.subplot(234), plt.axis('off'), plt.title("BNRF DFT Spectrum")plt.imshow(nrfDftAmpNorm, cmap='gray')# (6) 对频域滤波傅里叶变换执行傅里叶逆变换，并只取实部idft = np.zeros(dftAmp.shape, 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(dftAmp.shape)mask2[1::2, ::2] = -1mask2[::2, 1::2] = -1idftCen = idft * mask2 # g(x,y) * (-1)^(x+y)plt.subplot(235), plt.axis('off'), plt.title("g(x,y)*(-1)^(x+y)")plt.imshow(idftCen, cmap='gray')# (8) 截取左上角，大小和输入图像相等idftCenClip = np.clip(idftCen, 0, 255) # 截断函数，将数值限制在 [0,255]imgFiltered = idftCenClip.astype(np.uint8)imgFiltered = imgFiltered[:rows, :cols]plt.subplot(236), plt.axis('off'), plt.title("BNRF filtered image")plt.imshow(imgFiltered, cmap='gray')plt.tight_layout()plt.show()print("image.shape:{}".format(img.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("dftFilter.shape:{}".format(dftFilter.shape))print("imgFiltered.shape:{}".format(imgFiltered.shape))

程序说明：

傅立叶变换的频谱反映能量分布。通过中心化移频到原点以后，傅里叶变换的频谱图是以原点为中心对称分布的。中心化不仅可以清晰地看出图像频率分布，还可以分离出有周期性规律的干扰信号。

本例程中的图像受到正弦噪声的干扰，从中心化的频谱图可以看出，除原点以外还存在一对对称分布的亮点（箭头指示处），这就干扰噪声产生的。因此，在该亮点位置设计陷波带阻滤波器，可以消除干扰正弦噪声。

（本节完）