目录
选择性滤波带阻滤波器和带通滤波器陷波滤波器选择性滤波
处理特定的频带的滤波器称为频带滤波器
带阻滤波器：
若某个频带中的频率被滤除
带通滤波器：
若某个频带中的频率被通过
处理小频率矩形区域的滤波器称为陷波滤波器
陷波带阻滤波器：
若某个频带中的频率被拒绝
陷波带通滤波器：
若某个频带中的频率被通过
带阻滤波器和带通滤波器
频率域中的带通和带阻滤波器传递函数，可通过组合低通和高通滤波器传递函数来构建。然后高通滤波器也是由低通滤波器推导而来，所以说低通滤波器传递函数是形成高通、带阻、带通滤波器传递函数的基础。
可以由带阻滤波器传递函数获得带通滤波器传递函数
HBP(u,v)=1−HBR(u,v)(4.148)H_{BP}(u, v) = 1 - H_{BR}(u, v) \tag{4.148}HBP(u,v)=1−HBR(u,v)(4.148)
高斯带阻滤波器传递函数
H(u,v)=1−e−[(D(u,v)−C0)2W2](4.149)H(u, v) = 1 - e^{-\big[\frac{(D(u, v) - C_0)^2}{W^2} \ \ \big]} \tag{4.149}H(u,v)=1−e−[W2(D(u,v)−C0)2](4.149)
低于C0C_0C0时，该函数表现为一个低通高斯函数；等于C0C_0C0时，始终为0；高于C0C_0C0时，表现为一个高通高斯函数。但该函数在原点关不总是1。可以修改为：
H(u,v)=1−e−[D2(u,v)−C02D(u,v)W]2(4.150)H(u, v) = 1 - e^{-\big[\frac{D^2(u, v) - C_0^2}{D(u, v) W} \ \ \big]^2} \tag{4.150}H(u,v)=1−e−[D(u,v)WD2(u,v)−C02]2(4.150)
巴特沃斯带阻滤波器传递函数
H(u,v)=11+[D(u,v)WD2(u,v)−C02]2nH(u, v) = \frac{1}{1 + \bigg[\frac{D(u, v) W}{D^2(u, v) - C_0^2} \bigg]^{2n}}H(u,v)=1+[D2(u,v)−C02D(u,v)W]2n1
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
# 带阻滤波器传递函数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由高斯低通和高斯高通滤波器函数相加形成的带阻传递函数，最小值不是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 + GLPFhx_g = GLPF[500:, 500].flatten()# 3由式4.149得到的，原点处的值不是1GBRF_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()
# 理想、高斯、巴特沃斯带阻传递函数from mpl_toolkits.mplot3d import Axes3Dimport numpy as npfrom matplotlib import pyplot as pltfrom matplotlib import cmimg_temp = np.zeros([512, 512])center = img_temp.shaperadius = 128w = 60IBRF = 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()
陷波滤波器
陷波滤波器是最有用的选择性滤波
零相移滤波器必须关于原点(频率矩形中心)对称，中以为(u0,v0)(u_0, v_0)(u0,v0)的陷波滤波器传递函数在(−u0,−v0)(-u_0, -v_0)(−u0,−v0)位置必须有一个对应的陷波。陷波带阻滤波器传递函数可用中心被平移到陷波滤波中心的高通滤波器函数的乘积来产生
HNR(u,v)=∏k=1QHk(u,v)H−k(u,v)(4.151)H_{NR}(u, v) = \prod_{k=1}^Q H_k(u, v) H_{-k}(u, v) \tag{4.151}HNR(u,v)=k=1∏QHk(u,v)H−k(u,v)(4.151)
每个滤波器的距离计算公式为
Dk(u,v)=[(u−M/2−uk)2+(v−N/2−vk)2]1/2(4.152)D_{k}(u, v) = \big[(u - M / 2 - u_{k})^2 + (v - N / 2 - v_{k})^2 \big]^{1/2} \tag{4.152}Dk(u,v)=[(u−M/2−uk)2+(v−N/2−vk)2]1/2(4.152)
D−k(u,v)=[(u−M/2+uk)2+(v−N/2+vk)2]1/2(4.153)D_{-k}(u, v) = \big[(u - M / 2 + u_{k})^2 + (v - N / 2 + v_{k})^2 \big]^{1/2} \tag{4.153}D−k(u,v)=[(u−M/2+uk)2+(v−N/2+vk)2]1/2(4.153)
nnn阶巴特沃斯带阴滤波器
HNR(u,v)=∏k=13[11+[D0k/Dk(u,v)]n][11+[D0k/D−k(u,v)]n](4.154)H_{NR}(u, v) = \prod_{k=1}^3\bigg[ \frac{1}{1 + [D_{0k}/D_{k}(u,v)]^n} \bigg] \bigg[ \frac{1}{1 + [D_{0k}/D_{-k}(u,v)]^n} \bigg] \tag{4.154}HNR(u,v)=k=1∏3[1+[D0k/Dk(u,v)]n1][1+[D0k/D−k(u,v)]n1](4.154)
常数D0kD_{0k}D0k对每对陷波是相同的，但对不同的陷波对，它可以不同。
陷波带通滤波器传递函数可用陷波带阻滤波器得到
HNP(u,v)=1−HNR(u,v)(4.155)H_{NP}(u, v) = 1 - H_{NR}(u, v) \tag{4.155}HNP(u,v)=1−HNR(u,v)(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
# 巴特沃斯带阻陷波滤波器 BNRFimg_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()
# 使用陷波滤波删除数字化印刷图像中的莫尔模式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)# 对频谱做对数变换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_filterimg_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)# 对频谱做对数变换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)# # 对频谱做对数变换# spectrum_log = np.log(1 + spectrum)# 巴特沃斯陷波带阻滤波器NRF = narrow_notch_filter(fp, w=8, opening=20, vertical=True, horizontal=False)NRF = 1 - NRFfft_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()
