目錄

• • 選擇性濾波
  • • 帶阻濾波器和帶通濾波器
    • 陷波濾波器

選擇性濾波

處理特定的頻帶的濾波器稱為頻帶濾波器

• 帶阻濾波器：

  • 若某個(gè)頻帶中的頻率被濾除

• 帶通濾波器：

  • 若某個(gè)頻帶中的頻率被通過

處理小頻率矩形區(qū)域的濾波器稱為陷波濾波器

• 陷波帶阻濾波器：

  • 若某個(gè)頻帶中的頻率被拒絕

• 陷波帶通濾波器：

  • 若某個(gè)頻帶中的頻率被通過

帶阻濾波器和帶通濾波器

頻率域中的帶通和帶阻濾波器傳遞函數(shù)，可通過組合低通和高通濾波器傳遞函數(shù)來構(gòu)建。然后高通濾波器也是由低通濾波器推導(dǎo)而來，所以說低通濾波器傳遞函數(shù)是形成高通、帶阻、帶通濾波器傳遞函數(shù)的基礎(chǔ)。

可以由帶阻濾波器傳遞函數(shù)獲得帶通濾波器傳遞函數(shù)

$$H_{BP}(u,v)=1?H_{BR}(u,v)\text{\hspace{1em}(4.148)}$$

高斯帶阻濾波器傳遞函數(shù)

$$H(u,v)=1?e^{?\left[\frac{(D(u,v)?C_0{)}^2}{W^2}\right]}\text{\hspace{1em}(4.149)}$$

低于$C_0$時(shí)，該函數(shù)表現(xiàn)為一個(gè)低通高斯函數(shù)；等于$C_0$時(shí)，始終為0；高于$C_0$時(shí)，表現(xiàn)為一個(gè)高通高斯函數(shù)。但該函數(shù)在原點(diǎn)關(guān)不總是1。可以修改為：

$$H(u,v)=1?e^{?[\frac{D^2(u,v)?C_0^2}{D(u,v)W}{]}^2}\text{\hspace{1em}(4.150)}$$

巴特沃斯帶阻濾波器傳遞函數(shù)

$$H(u,v)=\frac{1}{1+[\frac{D(u,v)W}{D^2(u,v)?C_0^2}{]}^{2n}}$$

def idea_band_resistant_filter(source, center, radius=10, w=5):"""create idea band resistant filter param: source: input, source imageparam: center: input, the center of the filter, where is the lowest value, (0, 0) is top left corner, source.shape[:2] is center of the source imageparam: radius: input, int, the radius of circle of the band pass filter, default is 10param: w: input, int, the width of the band of the filter, default is 5return a [0, 1] value band resistant filter""" M, N = source.shape[1], source.shape[0]u = np.arange(M)v = np.arange(N)u, v = np.meshgrid(u, v)D = np.sqrt((u - center[1]//2)**2 + (v - center[0]//2)**2)D0 = radiushalf_w = w / 2kernel_1 = D.copy()assert radius > half_w, "radius must greater than W/2"#==================piecewise================kernel = np.piecewise(kernel_1, [kernel_1 <= D0 + half_w, kernel_1 <= D0 - half_w], [1, 0])kernel = 1 - kernel#==================where================== # kernel = np.where(kernel_1 > D0 + half_w, 1, kernel_1) # kernel = np.where(kernel <= D0 - half_w, 1, kernel) # kernel = np.where(kernel != 1, 0, kernel)# =================公式法================ # kernel_2 = D.copy() # kernel_1[D > D0 + half_w] = 1 # kernel_1[D <= D0 + half_w] = 0 # kernel_2[D > D0 - half_w] = 1 # kernel_2[D <= D0 - half_w] = 0 # kernel = kernel_1 - kernel_2 return kernel def gauss_band_resistant_149(source, center, radius=10, w=5):"""create gaussian band resistant filter, equation 4.149param: source: input, source imageparam: center: input, the center of the filter, where is the lowest value, (0, 0) is top left corner, source.shape[:2] is center of the source imageparam: radius: input, int, the radius of circle of the band pass filter, default is 10param: w: input, int, the width of the band of the filter, default is 5return a [0, 1] value gaussian band resistant filter""" N, M = source.shape[:2]u = np.arange(M)v = np.arange(N)u, v = np.meshgrid(u, v)D = np.sqrt((u - center[1]//2)**2 + (v - center[0]//2)**2)C0 = radiuskernel = 1 - np.exp(-(D - C0)**2 / (w**2))return kernel def gauss_band_resistant_filter(source, center, radius=10, w=5):"""create gaussian band resistant filter, equation 4.150param: source: input, source imageparam: center: input, the center of the filter, where is the lowest value, (0, 0) is top left corner, source.shape[:2] is center of the source imageparam: radius: input, int, the radius of circle of the band pass filter, default is 10param: w: input, int, the width of the band of the filter, default is 5return a [0, 1] value gaussian band resistant filter""" N, M = source.shape[:2]u = np.arange(M)v = np.arange(N)u, v = np.meshgrid(u, v)D = np.sqrt((u - center[1]//2)**2 + (v - center[0]//2)**2)C0 = radiuskernel = 1 - np.exp(-((D**2 - C0**2) / (D * w))**2)return kernel def butterworth_band_resistant_filter(source, center, radius=10, w=5, n=1):"""create butterworth band resistant filter, equation 4.150param: source: input, source imageparam: center: input, the center of the filter, where is the lowest value, (0, 0) is top left corner, source.shape[:2] is center of the source imageparam: radius: input, int, the radius of circle of the band pass filter, default is 10param: w: input, int, the width of the band of the filter, default is 5param: n: input, int, order of the butter worth fuction, return a [0, 1] value butterworth band resistant filter""" N, M = source.shape[:2]u = np.arange(M)v = np.arange(N)u, v = np.meshgrid(u, v)D = np.sqrt((u - center[1]//2)**2 + (v - center[0]//2)**2)C0 = radiustemp = (D * w) / (D**2 - C0**2)kernel = 1 / (1 + temp ** (2*n)) return kernel # 帶阻濾波器傳遞函數(shù) img_temp = np.zeros([1000, 1000]) C0 = 100# 1理想帶阻濾波器 IBRF = idea_band_resistant_filter(img_temp, img_temp.shape, radius=C0, w=100) hx_i = IBRF[500:, 500].flatten()# 2由高斯低通和高斯高通濾波器函數(shù)相加形成的帶阻傳遞函數(shù)，最小值不是0，并且與C0不重合 GHPF = gauss_high_pass_filter(img_temp, img_temp.shape, radius=C0) GLPF = gauss_low_pass_filter(img_temp, img_temp.shape, radius=C0/2) GLPF = GHPF + GLPF hx_g = GLPF[500:, 500].flatten()# 3由式4.149得到的，原點(diǎn)處的值不是1 GBRF_149 = gauss_band_resistant_149(img_temp, img_temp.shape, radius=C0, w=100) hx_g149 = GBRF_149[500:, 500].flatten()# 4由式4.150得到的 GBRF = gauss_band_resistant_filter(img_temp, img_temp.shape, radius=C0, w=100) hx_gbrf = GBRF[500:, 500].flatten()fig = plt.figure(figsize=(16, 3)) ax_1 = fig.add_subplot(1, 4, 1) ax_1.plot(hx_i), ax_1.set_yticks([0, 1.0]), ax_1.set_xticks([100, 500]), ax_1.set_ylim(0, 1.1), ax_1.set_xlim(0, 500)ax_2 = fig.add_subplot(1, 4, 2) ax_2.plot(hx_g), ax_2.set_yticks([0.75, 1.0]), ax_2.set_xticks([100, 500]), ax_2.set_ylim(0.75, 1.1), ax_2.set_xlim(0, 500)ax_3 = fig.add_subplot(1, 4, 3) ax_3.plot(hx_g149), ax_3.set_yticks([0, 1.0]), ax_3.set_xticks([100, 500]), ax_3.set_ylim(0, 1.1), ax_3.set_xlim(0, 500)ax_3 = fig.add_subplot(1, 4, 4) ax_3.plot(hx_gbrf), ax_3.set_yticks([0, 1.0]), ax_3.set_xticks([100, 500]), ax_3.set_ylim(0, 1.1), ax_3.set_xlim(0, 500)plt.tight_layout() plt.show()

# 理想、高斯、巴特沃斯帶阻傳遞函數(shù) from mpl_toolkits.mplot3d import Axes3D import numpy as np from matplotlib import pyplot as plt from matplotlib import cmimg_temp = np.zeros([512, 512]) center = img_temp.shaperadius = 128 w = 60 IBRF = idea_band_resistant_filter(img_temp, img_temp.shape, radius=radius, w=w) GBFR = gauss_band_resistant_filter(img_temp, img_temp.shape, radius=radius, w=w) BBRF = butterworth_band_resistant_filter(img_temp, img_temp.shape, radius=radius, w=w, n=1)filters = ['IBRF', 'GBFR', 'BBRF'] # 用來繪制3D圖 M, N = img_temp.shape[1], img_temp.shape[0] u = np.arange(M) v = np.arange(N) u, v = np.meshgrid(u, v)fig = plt.figure(figsize=(15, 15))for i in range(len(filters)):ax_1 = fig.add_subplot(3, 3, i*3 + 1, projection='3d')plot_3d(ax_1, u, v, eval(filters[i]))ax_2 = fig.add_subplot(3, 3, i*3 + 2)ax_2.imshow(eval(filters[i]),'gray'), ax_2.set_title(filters[i]), ax_2.set_xticks([]), ax_2.set_yticks([])h_1 = eval(filters[i])[img_temp.shape[0]//2:, img_temp.shape[1]//2]ax_3 = fig.add_subplot(3, 3, i*3 + 3)ax_3.plot(h_1), ax_3.set_xticks([0, radius//2]), ax_3.set_yticks([0, 1]), ax_3.set_xlim([0, 320]), ax_3.set_ylim([0, 1.2]) plt.tight_layout() plt.show()

陷波濾波器

陷波濾波器是最有用的選擇性濾波

零相移濾波器必須關(guān)于原點(diǎn)(頻率矩形中心)對(duì)稱，中以為$(u_0,v_0)$的陷波濾波器傳遞函數(shù)在$(?u_0,?v_0)$位置必須有一個(gè)對(duì)應(yīng)的陷波。陷波帶阻濾波器傳遞函數(shù)可用中心被平移到陷波濾波中心的高通濾波器函數(shù)的乘積來產(chǎn)生

$$H_{NR}(u,v)=\prod _{k=1}^Q{H}_k(u,v)H_{?k}(u,v)\text{\hspace{1em}(4.151)}$$

每個(gè)濾波器的距離計(jì)算公式為
$$D_k(u,v)=[(u?M/2?u_k{)}^2+(v?N/2?v_k{)}^2{]}^{1/2}\text{\hspace{1em}(4.152)}$$
$$D_{?k}(u,v)=[(u?M/2+u_k{)}^2+(v?N/2+v_k{)}^2{]}^{1/2}\text{\hspace{1em}(4.153)}$$

$n$階巴特沃斯帶陰濾波器
$$H_{NR}(u,v)=\prod _{k=1}^3\left[\frac{1}{1+[D_{0k}/D_k(u,v){]}^n}\right]\left[\frac{1}{1+[D_{0k}/D_{?k}(u,v){]}^n}\right]\text{\hspace{1em}(4.154)}$$

常數(shù)$D_{0k}$對(duì)每對(duì)陷波是相同的，但對(duì)不同的陷波對(duì)，它可以不同。

陷波帶通濾波器傳遞函數(shù)可用陷波帶阻濾波器得到
$$H_{NP}(u,v)=1?H_{NR}(u,v)\text{\hspace{1em}(4.155)}$$

def butterworth_notch_resistant_filter(img, uk, vk, radius=10, n=1):"""create butterworth notch resistant filter, equation 4.155param: img: input, source imageparam: uk: input, int, center of the heightparam: vk: input, int, center of the widthparam: radius: input, int, the radius of circle of the band pass filter, default is 10param: w: input, int, the width of the band of the filter, default is 5param: n: input, int, order of the butter worth fuction, return a [0, 1] value butterworth band resistant filter""" M, N = img.shape[1], img.shape[0]u = np.arange(M)v = np.arange(N)u, v = np.meshgrid(u, v)DK = np.sqrt((u - M//2 - uk)**2 + (v - N//2 - vk)**2)D_K = np.sqrt((u - M//2 + uk)**2 + (v - N//2 + vk)**2)D0 = radiuskernel = (1 / (1 + (D0 / (DK+1e-5))**n)) * (1 / (1 + (D0 / (D_K+1e-5))**n))return kernel # 巴特沃斯帶阻陷波濾波器 BNRF img_temp = np.zeros([512, 512]) BNF_1 = butterworth_notch_resistant_filter(img_temp, radius=10, uk=30, vk=40, n=3) BNF_2 = butterworth_notch_resistant_filter(img_temp, radius=10, uk=30, vk=80, n=3) BNF_3 = butterworth_notch_resistant_filter(img_temp, radius=10, uk=-30, vk=80, n=3)plt.figure(figsize=(16, 12)) plt.subplot(141), plt.imshow(BNF_1, 'gray'), plt.title('BNF_1') plt.subplot(142), plt.imshow(BNF_2, 'gray'), plt.title('BNF_2') plt.subplot(143), plt.imshow(BNF_3, 'gray'), plt.title('BNF_3')BNF_dst = BNF_1 * BNF_2 * BNF_3plt.subplot(144), plt.imshow(BNF_dst, 'gray'), plt.title('BNF_dst')plt.tight_layout() plt.show()

# 使用陷波濾波刪除數(shù)字化印刷圖像中的莫爾模式 img_ori = cv2.imread("DIP_Figures/DIP3E_Original_Images_CH04/Fig0464(a)(car_75DPI_Moire).tif", 0)M, N = img_ori.shape[:2]# 填充 fp = pad_image(img_ori, mode='reflect') # 中心化 fp_cen = centralized_2d(fp) # 正變換 fft = np.fft.fft2(fp_cen)# 頻譜 spectrum = spectrum_fft(fft) # 對(duì)頻譜做對(duì)數(shù)變換 spectrum_log = np.log(1 + spectrum)# 巴特沃斯陷波帶阻濾波器 BNRF_1 = butterworth_notch_resistant_filter(fp, radius=9, uk=60, vk=80, n=4) BNRF_2 = butterworth_notch_resistant_filter(fp, radius=9, uk=-60, vk=80, n=4) BNRF_3 = butterworth_notch_resistant_filter(fp, radius=9, uk=60, vk=160, n=4) BNRF_4 = butterworth_notch_resistant_filter(fp, radius=9, uk=-60, vk=160, n=4)BNRF = BNRF_1 * BNRF_2 * BNRF_3 * BNRF_4 fft_filter = fft * BNRF# 濾波后的頻譜 spectrum_filter = spectrum_fft(fft_filter) spectrum_filter_log = np.log(1 + spectrum_filter)# 傅里葉反變換 ifft = np.fft.ifft2(fft_filter)# 去中心化反變換的圖像，并取左上角的圖像 img_new = centralized_2d(ifft.real)[:M, :N] img_new = np.clip(img_new, 0, img_new.max()) img_new = np.uint8(normalize(img_new) * 255)fig = plt.figure(figsize=(10, 14))ax_1 = fig.add_subplot(2, 2, 1) ax_1.imshow(img_ori, 'gray'), ax_1.set_title('Original'), ax_1.set_xticks([]), ax_1.set_yticks([])ax_2 = fig.add_subplot(2, 2, 2) ax_2.imshow(spectrum_log, 'gray'), ax_2.set_title('Spectrum Before Filter'), ax_2.set_xticks([]), ax_2.set_yticks([])ax_3 = fig.add_subplot(2, 2, 3) ax_3.imshow(spectrum_filter_log, 'gray'), ax_3.set_title('Spectrum After Filter'), ax_3.set_xticks([]), ax_3.set_yticks([])ax_4 = fig.add_subplot(2, 2, 4) ax_4.imshow(img_new, 'gray'), ax_4.set_title('Denoising'), ax_4.set_xticks([]), ax_4.set_yticks([])plt.tight_layout() plt.show()

使用陷波濾波去除周期干擾

def narrow_notch_filter(img, w=5, opening=10, vertical=True, horizontal=False):"""create narrow notch resistant filter, using opencvparam: img: input, source imageparam: w: input, int, width of the resistant, value is 0, default is 5param: opening: input, int, opening of the resistant, value is 1, default is 10param: vertical: input, boolean, whether vertical or not, default is "True"param: horizontal: input, boolean, whether horizontal or not, default is "False"return a [0, 1] value butterworth band resistant filter""" dst = np.ones(img.shape, dtype=np.uint8) * 255c_height, c_width = img.shape[0] // 2, img.shape[1] // 2if vertical:cv2.rectangle(dst, ((img.shape[1] - w)//2, 0), (c_width + w//2, img.shape[0]), (0), -1)cv2.rectangle(dst, (0, (img.shape[0] - opening)//2), (img.shape[1], c_height + opening//2), (255), -1)horizontal_ = np.ones(img.shape, dtype=np.uint8) * 255if horizontal: cv2.rectangle(horizontal_, (0, (img.shape[0] - w)//2), (img.shape[1], c_height + w//2), (0), -1)cv2.rectangle(horizontal_, ((img.shape[1] - opening)//2, 0), (c_width + opening//2, img.shape[0]), (255), -1)dst = dst * horizontal_dst = dst / dst.max()return dst def narrow_notch_filter(img, w=5, opening=10, vertical=True, horizontal=False):"""create narrow notch resistant filterparam: img: input, source imageparam: w: input, int, width of the resistant, value is 0, default is 5param: opening: input, int, opening of the resistant, value is 1, default is 10param: vertical: input, boolean, whether vertical or not, default is "True"param: horizontal: input, boolean, whether horizontal or not, default is "False"return a [0, 1] value butterworth band resistant filter""" assert w > 0, "W must greater than 0"w_half = w//2opening_half = opening//2img_temp = np.ones(img.shape[:2])N, M = img_temp.shape[:]img_vertical = img_temp.copy()img_horizontal = img_temp.copy()if horizontal:img_horizontal[M//2 - w_half:M//2 + w - w_half, :] = 0img_horizontal[:, N//2 - opening_half:N//2 + opening - opening_half] = 1if vertical:img_vertical[:, N//2 - w_half:N//2 + w - w_half] = 0img_vertical[M//2 - opening_half:M//2 + opening - opening_half, :] = 1img_dst = img_horizontal * img_verticalreturn img_dst # NNF narrow_notch_filter img_temp = np.zeros([512, 512]) NNF = narrow_notch_filter(img_temp, 5, 20, vertical=True, horizontal=False) plt.figure(figsize=(10, 8)) plt.imshow(NNF,'gray'),plt.title('NNF') plt.show()

# 使用陷波濾波去除周期干擾 img_ori = cv2.imread("DIP_Figures/DIP3E_Original_Images_CH04/Fig0465(a)(cassini).tif", 0)M, N = img_ori.shape[:2]# 填充為'constant'可得到跟書上一樣的頻譜，這里使用'reflect' fp = pad_image(img_ori, mode='constant') # 中心化 fp_cen = centralized_2d(fp) # 正變換 fft = np.fft.fft2(fp_cen)# 頻譜 spectrum = spectrum_fft(fft) # 對(duì)頻譜做對(duì)數(shù)變換 spectrum_log = np.log(1 + spectrum)# 巴特沃斯陷波帶阻濾波器 NRF = narrow_notch_filter(fp, w=8, opening=20, vertical=True, horizontal=False)fft_filter = fft * NRF# 濾波后的頻譜 spectrum_filter = spectrum_fft(fft_filter) spectrum_filter_log = np.log(1 + spectrum_filter)# 傅里葉反變換 ifft = np.fft.ifft2(fft_filter)# 去中心化反變換的圖像，并取左上角的圖像 img_new = centralized_2d(ifft.real)[:M, :N] img_new = np.uint8(normalize(img_new) * 255)fig = plt.figure(figsize=(10, 10))ax_1 = fig.add_subplot(2, 2, 1) ax_1.imshow(img_ori, 'gray'), ax_1.set_title('Original'), ax_1.set_xticks([]), ax_1.set_yticks([])ax_2 = fig.add_subplot(2, 2, 2) ax_2.imshow(spectrum_log, 'gray'), ax_2.set_title('Spectrum Before Filter'), ax_2.set_xticks([]), ax_2.set_yticks([])ax_3 = fig.add_subplot(2, 2, 3) ax_3.imshow(spectrum_filter_log, 'gray'), ax_3.set_title('Spectrum After Filter'), ax_3.set_xticks([]), ax_3.set_yticks([])ax_4 = fig.add_subplot(2, 2, 4) ax_4.imshow(img_new, 'gray'), ax_4.set_title('Denoising'), ax_4.set_xticks([]), ax_4.set_yticks([])plt.tight_layout() plt.show()

# 周期干擾的空間模式 img_ori = cv2.imread("DIP_Figures/DIP3E_Original_Images_CH04/Fig0465(a)(cassini).tif", 0) M, N = img_ori.shape[:2]# 填充為'constant'可得到跟書上一樣的頻譜，這里使用'reflect' fp = pad_image(img_ori, mode='constant') # 中心化 fp_cen = centralized_2d(fp) # 正變換 fft = np.fft.fft2(fp_cen)# # 頻譜 # spectrum = spectrum_fft(fft) # # 對(duì)頻譜做對(duì)數(shù)變換 # spectrum_log = np.log(1 + spectrum)# 巴特沃斯陷波帶阻濾波器 NRF = narrow_notch_filter(fp, w=8, opening=20, vertical=True, horizontal=False) NRF = 1 - NRF fft_filter = fft * NRF# 濾波后的頻譜 spectrum_filter = spectrum_fft(fft_filter) spectrum_filter_log = np.log(1 + spectrum_filter)# 傅里葉反變換 ifft = np.fft.ifft2(fft_filter)# 去中心化反變換的圖像，并取左上角的圖像 img_new = centralized_2d(ifft.real)[:M, :N] img_new = np.uint8(normalize(img_new) * 255)fig = plt.figure(figsize=(10, 10))# ax_1 = fig.add_subplot(2, 2, 1) # ax_1.imshow(img_ori, 'gray'), ax_1.set_title('Original'), ax_1.set_xticks([]), ax_1.set_yticks([])# ax_2 = fig.add_subplot(2, 2, 2) # ax_2.imshow(spectrum_log, 'gray'), ax_2.set_title('Spectrum Before Filter'), ax_2.set_xticks([]), ax_2.set_yticks([])ax_3 = fig.add_subplot(2, 2, 3) ax_3.imshow(NRF, 'gray'), ax_3.set_title('Spectrum After Filter'), ax_3.set_xticks([]), ax_3.set_yticks([])ax_4 = fig.add_subplot(2, 2, 4) ax_4.imshow(img_new, 'gray'), ax_4.set_title('Noise'), ax_4.set_xticks([]), ax_4.set_yticks([])plt.tight_layout() plt.show()

總結(jié)

以上是生活随笔為你收集整理的第4章 Python 数字图像处理(DIP) - 频率域滤波12 - 选择性滤波 - 带阻的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。