標(biāo)題

• • 使用頻率域濾波降低周期噪聲
  • • 陷波濾波深入介紹
    • 最優(yōu)陷波濾波

本章陷波濾波器有部分得出的結(jié)果不佳，如果有更好的解決方案，請賜教，不勝感激。

使用頻率域濾波降低周期噪聲

陷波濾波深入介紹

零相移濾波器必須關(guān)于原點(頻率矩形中心)對稱，中以為$(u_0,v_0)$的陷波濾波器傳遞函數(shù)在$(?u_0,?v_0)$位置必須有一個對應(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}(5.33)}$$

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

3個$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}(5.36)}$$

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

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

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 def idea_notch_resistant_filter(source, uk, vk, radius=5):"""create idea notch resistant filter param: source: input, source imageparam: uk: input, int, center of the heightparam: vk: input, int, center of the widthparam: radius: input, the radius of the lowest value, greater value, bigger blocker out range, if the radius is 0, then allvalue is 0return a [0, 1] value filter"""M, N = source.shape[1], source.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 = radiusk_1 = DK.copy()k_2 = D_K.copy()k_1[DK > D0] = 1k_1[DK <= D0] = 0k_2[D_K > D0] = 1k_2[D_K <= D0] = 0kernel = k_1 * k_2return kernel def gauss_notch_resistant_filter(source, uk, vk, radius=5):"""create gauss low pass filter param: source: input, source imageparam: uk: input, int, center of the heightparam: vk: input, int, center of the widthparam: radius: input, the radius of the lowest value, greater value, bigger blocker out range, if the radius is 0, then allvalue is 0return a [0, 1] value filter""" M, N = source.shape[1], source.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 = radiusk_1 = 1 - np.exp(- (DK**2)/(D0**2)) k_2 = 1 - np.exp(- (D_K**2)/(D0**2)) kernel = k_1 * k_2return kernel def plot_3d(ax, x, y, z, cmap):ax.plot_surface(x, y, z, antialiased=True, shade=True, cmap=cmap)ax.view_init(20, -20), ax.grid(b=False), ax.set_xticks([]), ax.set_yticks([]), ax.set_zticks([]) # 理想、高斯、巴特沃斯陷波濾波器 from mpl_toolkits.mplot3d import Axes3D import numpy as np from matplotlib import pyplot as plt from matplotlib import cmimg_temp = np.zeros([256, 256])INRF = idea_notch_resistant_filter(img_temp, radius=20, uk=30, vk=80) GNRF = gauss_notch_resistant_filter(img_temp, radius=20, uk=30, vk=80) BNRF = butterworth_notch_resistant_filter(img_temp, radius=20, uk=30, vk=80, n=5)# 用來繪制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=(21, 7)) ax_1 = fig.add_subplot(1, 3, 1, projection='3d') plot_3d(ax_1, u, v, INRF, cmap=cm.plasma)ax_1 = fig.add_subplot(1, 3, 2, projection='3d') plot_3d(ax_1, u, v, GNRF, cmap=cm.PiYG)ax_1 = fig.add_subplot(1, 3, 3, projection='3d') plot_3d(ax_1, u, v, BNRF, cmap=cm.PiYG) plt.tight_layout() plt.show()

def centralized_2d(arr):"""centralized 2d array $f(x, y) (-1)^{x + y}$, about 4.5 times faster than index, 160 times faster than loop,"""def get_1_minus1(img):"""get 1 when image index is even, -1 while index is odd, same shape as input image, need this array to multiply with input imageto get centralized image for DFTParameter: img: input, here we only need img shape to create the arrayreturn such a [[1, -1, 1], [-1, 1, -1]] array, example is 3x3"""dst = np.ones(img.shape)dst[1::2, ::2] = -1dst[::2, 1::2] = -1return dst#==================中心化=============================mask = get_1_minus1(arr)dst = arr * maskreturn dst def pad_image(img, mode='constant'):"""pad image into PxQ shape, orginal is in the top left corner.param: img: input imageparam: mode: input, str, is numpy pad parameter, default is 'constant'. for more detail please refere to Numpy padreturn PxQ shape image padded with zeros or other values"""dst = np.pad(img, ((0, img.shape[0]), (0, img.shape[1])), mode=mode)return dst def add_sin_noise(img, scale=1, angle=0):"""add sin noise for imageparam: img: input image, 1 channel, dtype=uint8param: scale: sin scaler, smaller than 1, will enlarge, bigger than 1 will shrinkparam: angle: angle of the rotationreturn: output_img: output image is [0, 1] image which you could use as mask or any you want to"""height, width = img.shape[:2] # original image shape# convert all the angleif int(angle / 90) % 2 == 0:rotate_angle = angle % 90else:rotate_angle = 90 - (angle % 90)rotate_radian = np.radians(rotate_angle) # convert angle to radian# get new image height and widthnew_height = int(np.ceil(height * np.cos(rotate_radian) + width * np.sin(rotate_radian)))new_width = int(np.ceil(width * np.cos(rotate_radian) + height * np.sin(rotate_radian))) # if new height or new width less than orginal height or width, the output image will be not the same shape as input, here set it rightif new_height < height:new_height = heightif new_width < width:new_width = width# meshgridu = np.arange(new_width)v = np.arange(new_height)u, v = np.meshgrid(u, v)# get sin noise image, you could use scale to make some difference, better you could add some shift # noise = abs(np.sin(u * scale))noise = 1 - np.sin(u * scale)# here use opencv to get rotation, better write yourself rotation functionC1 = cv2.getRotationMatrix2D((new_width/2.0, new_height/2.0), angle, 1)new_img = cv2.warpAffine(noise, C1, (int(new_width), int(new_height)), borderValue=0)# ouput image should be the same shape as input, so caculate the offset the output image and the new image# I make new image bigger so that it will cover all output imageoffset_height = abs(new_height - height) // 2offset_width = abs(new_width - width) // 2img_dst = new_img[offset_height:offset_height + height, offset_width:offset_width+width]output_img = normalize(img_dst)return output_img def spectrum_fft(fft):"""return FFT spectrum"""return np.sqrt(np.power(fft.real, 2) + np.power(fft.imag, 2)) # 陷波濾波器處理周期噪聲 img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH05/Fig0507(a)(ckt-board-orig).tif', 0) #直接讀為灰度圖像# 正弦噪聲 noise = add_sin_noise(img_ori, scale=0.35, angle=-20) img = np.array(img_ori / 255, np.float32) img_noise = img + noise img_noise = np.uint8(normalize(img_noise)*255)# 頻率域中的其他特性 # FFT img_fft = np.fft.fft2(img_noise.astype(np.float32)) # 中心化 fshift = np.fft.fftshift(img_fft) # 將變換的頻率圖像四角移動到中心 # 中心化后的頻譜 spectrum_fshift = spectrum_fft(fshift) spectrum_fshift_n = np.uint8(normalize(spectrum_fshift) * 255)# 對頻譜做對數(shù)變換 spectrum_log = np.log(1 + spectrum_fshift)BNRF = butterworth_notch_resistant_filter(img_ori, radius=5, uk=25, vk=10, n=4)f1shift = fshift * (BNRF) f2shift = np.fft.ifftshift(f1shift) #對新的進(jìn)行逆變換 img_new = np.fft.ifft2(f2shift) img_new = np.abs(img_new)plt.figure(figsize=(15, 15)) plt.subplot(221), plt.imshow(img_noise, 'gray'), plt.title('With Sine noise'), plt.xticks([]),plt.yticks([]) plt.subplot(222), plt.imshow(spectrum_log, 'gray'), plt.title('Spectrum'), plt.xticks([]),plt.yticks([]) # 在圖像上加上箭頭 plt.arrow(180, 180, 25, 30, width=5,length_includes_head=True, shape='full') plt.arrow(285, 265, -25, -30, width=5,length_includes_head=True, shape='full')plt.subplot(223), plt.imshow(BNRF, 'gray'), plt.title('Spectrum'), plt.xticks([]),plt.yticks([]) # 在圖像上加上箭頭 plt.arrow(180, 180, 25, 30, width=5,length_includes_head=True, shape='full') plt.arrow(285, 265, -25, -30, width=5,length_includes_head=True, shape='full')plt.subplot(224), plt.imshow(img_new, 'gray'), plt.title('Spectrum'), plt.xticks([]),plt.yticks([])plt.tight_layout() plt.show()

# 陷波濾波器提取周期噪聲 img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH05/Fig0507(a)(ckt-board-orig).tif', 0) #直接讀為灰度圖像# 正弦噪聲 noise = add_sin_noise(img_ori, scale=0.35, angle=-20) img = np.array(img_ori / 255, np.float32) img_noise = img + noise img_noise = np.uint8(normalize(img_noise)*255)# 頻率域中的其他特性 # FFT img_fft = np.fft.fft2(img_noise.astype(np.float32)) # 中心化 fshift = np.fft.fftshift(img_fft) # 將變換的頻率圖像四角移動到中心 # 中心化后的頻譜 spectrum_fshift = spectrum_fft(fshift) spectrum_fshift_n = np.uint8(normalize(spectrum_fshift) * 255)# 對頻譜做對數(shù)變換 spectrum_log = np.log(1 + spectrum_fshift)BNRF = 1 - butterworth_notch_resistant_filter(img_ori, radius=5, uk=25, vk=10, n=4)f1shift = fshift * (BNRF) f2shift = np.fft.ifftshift(f1shift) #對新的進(jìn)行逆變換 img_new = np.fft.ifft2(f2shift) img_new = np.abs(img_new)plt.figure(figsize=(15, 15)) # plt.subplot(221), plt.imshow(img_noise, 'gray'), plt.title('With Sine noise'), plt.xticks([]),plt.yticks([]) # plt.subplot(222), plt.imshow(spectrum_log, 'gray'), plt.title('Spectrum'), plt.xticks([]),plt.yticks([]) # # 在圖像上加上箭頭 # plt.arrow(180, 180, 25, 30, width=5,length_includes_head=True, shape='full') # plt.arrow(285, 265, -25, -30, width=5,length_includes_head=True, shape='full')# plt.subplot(223), plt.imshow(BNRF, 'gray'), plt.title('Spectrum'), plt.xticks([]),plt.yticks([]) # # 在圖像上加上箭頭 # plt.arrow(180, 180, 25, 30, width=5,length_includes_head=True, shape='full') # plt.arrow(285, 265, -25, -30, width=5,length_includes_head=True, shape='full')plt.subplot(224), plt.imshow(img_new, 'gray'), plt.title('Sine pattern'), plt.xticks([]),plt.yticks([])plt.tight_layout() plt.show()

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""" epsilon = 1e-8N, 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) + epsilon)kernel = 1 / (1 + temp ** (2*n)) return kerneldef butterworth_low_pass_filter(img, center, radius=5, n=1):"""create butterworth low pass 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, the radius of the lowest value, greater value, bigger blocker out range, if the radius is 0, then allvalue is 0param: n: input, float, the order of the filter, if n is small, then the BLPF will be close to GLPF, and more smooth from lowfrequency to high freqency.if n is large, will close to ILPFreturn a [0, 1] value filter""" epsilon = 1e-8M, N = img.shape[1], img.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 = radiuskernel = (1 / (1 + (D / (D0 + epsilon))**(2*n)))return kernel # 陷波濾波器處理周期噪聲，用巴特沃斯低通濾波器得到的效果比目前的陷波濾波器效果還要好 img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH05/Fig0516(a)(applo17_boulder_noisy).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_fshift = spectrum_fft(fft) spectrum_log = np.log(1 + spectrum_fshift)# 濾波器 n = 15 r = 20 H = butterworth_low_pass_filter(fp, fp.shape, radius=100, n=4) # BNRF_1 = butterworth_notch_resistant_filter(fp, radius=r, uk=355, vk=0, n=n) # BNRF_2 = butterworth_notch_resistant_filter(fp, radius=r, uk=0, vk=355, n=n) # BNRF_3 = butterworth_notch_resistant_filter(fp, radius=r, uk=250, vk=250, n=n) # BNRF_4 = butterworth_notch_resistant_filter(fp, radius=r, uk=-250, vk=250, n=n) # BNRF = BNRF_1 * BNRF_2 * BNRF_3 * BNRF_4 * H BNRF = Hfft_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)plt.figure(figsize=(15, 12)) plt.subplot(221), plt.imshow(img_ori, 'gray'), plt.title('With noise'), plt.xticks([]),plt.yticks([]) plt.subplot(222), plt.imshow(spectrum_log, 'gray'), plt.title('Spectrum'), plt.xticks([]),plt.yticks([]) plt.subplot(223), plt.imshow(spectrum_filter_log, 'gray'), plt.title('Spectrum With Filter'), plt.xticks([]),plt.yticks([]) plt.subplot(224), plt.imshow(img_new, 'gray'), plt.title('IDFT'), plt.xticks([]),plt.yticks([]) plt.tight_layout() plt.show()

# 陷波濾波器提取周期噪聲 img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH05/Fig0516(a)(applo17_boulder_noisy).tif', 0) M, N = img_ori.shape[:2]fp = pad_image(img_ori, mode='constant') fp_cen = centralized_2d(fp) fft = np.fft.fft2(fp_cen)# 中心化后的頻譜 spectrum_fshift = spectrum_fft(fft) spectrum_log = np.log(1 + spectrum_fshift)# 濾波器 n = 15 r = 20 H = butterworth_low_pass_filter(fp, fp.shape, radius=100, n=3) # BNRF_1 = butterworth_notch_resistant_filter(fp, radius=r, uk=355, vk=0, n=n) # BNRF_2 = butterworth_notch_resistant_filter(fp, radius=r, uk=0, vk=355, n=n) # BNRF_3 = butterworth_notch_resistant_filter(fp, radius=r, uk=250, vk=250, n=n) # BNRF_4 = butterworth_notch_resistant_filter(fp, radius=r, uk=-250, vk=250, n=n) # BNRF = BNRF_1 * BNRF_2 * BNRF_3 * BNRF_4 * H BNRF = H fft_filter = fft * (1 - 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)plt.figure(figsize=(15, 12)) # plt.subplot(221), plt.imshow(img_ori, 'gray'), plt.title('With Sine noise'), plt.xticks([]),plt.yticks([]) # plt.subplot(222), plt.imshow(spectrum_log, 'gray'), plt.title('Spectrum'), plt.xticks([]),plt.yticks([]) # plt.subplot(223), plt.imshow(spectrum_filter_log, 'gray'), plt.title('Spectrum'), plt.xticks([]),plt.yticks([]) plt.subplot(224), plt.imshow(img_new, 'gray'), plt.title('Spectrum'), plt.xticks([]),plt.yticks([]) plt.tight_layout() plt.show()

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])M, N = 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 # 陷波濾波器處理周期噪聲，用巴特沃斯低通濾波器得到的效果比目前的陷波濾波器效果還要好 img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH05/Fig0519(a)(florida_satellite_original).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_fshift = spectrum_fft(fft) spectrum_log = np.log(1 + spectrum_fshift)# 濾波器 n = 15 r = 20 H = narrow_notch_filter(fp, w=10, opening=30, vertical=True, horizontal=False) # BNRF_1 = butterworth_notch_resistant_filter(fp, radius=r, uk=355, vk=0, n=n) # BNRF_2 = butterworth_notch_resistant_filter(fp, radius=r, uk=0, vk=355, n=n) # BNRF_3 = butterworth_notch_resistant_filter(fp, radius=r, uk=250, vk=250, n=n) # BNRF_4 = butterworth_notch_resistant_filter(fp, radius=r, uk=-250, vk=250, n=n) # BNRF = BNRF_1 * BNRF_2 * BNRF_3 * BNRF_4 * H BNRF = Hfft_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)plt.figure(figsize=(15, 16)) plt.subplot(221), plt.imshow(img_ori, 'gray'), plt.title('With noise'), plt.xticks([]),plt.yticks([]) plt.subplot(222), plt.imshow(spectrum_log, 'gray'), plt.title('Spectrum'), plt.xticks([]),plt.yticks([]) plt.subplot(223), plt.imshow(spectrum_filter_log, 'gray'), plt.title('Spectrum With Filter'), plt.xticks([]),plt.yticks([]) plt.subplot(224), plt.imshow(img_new, 'gray'), plt.title('IDFT'), plt.xticks([]),plt.yticks([]) plt.tight_layout() plt.show()

# 陷波濾波器提取周期噪聲，用巴特沃斯低通濾波器得到的效果比目前的陷波濾波器效果還要好 img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH05/Fig0519(a)(florida_satellite_original).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_fshift = spectrum_fft(fft) spectrum_log = np.log(1 + spectrum_fshift)# 濾波器 n = 15 r = 20 H = narrow_notch_filter(fp, w=10, opening=30, vertical=True, horizontal=False) # BNRF_1 = butterworth_notch_resistant_filter(fp, radius=r, uk=355, vk=0, n=n) # BNRF_2 = butterworth_notch_resistant_filter(fp, radius=r, uk=0, vk=355, n=n) # BNRF_3 = butterworth_notch_resistant_filter(fp, radius=r, uk=250, vk=250, n=n) # BNRF_4 = butterworth_notch_resistant_filter(fp, radius=r, uk=-250, vk=250, n=n) # BNRF = BNRF_1 * BNRF_2 * BNRF_3 * BNRF_4 * H BNRF = Hfft_filter = fft * (1 - 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)plt.figure(figsize=(15, 16)) # plt.subplot(221), plt.imshow(img_ori, 'gray'), plt.title('With noise'), plt.xticks([]),plt.yticks([]) # plt.subplot(222), plt.imshow(spectrum_log, 'gray'), plt.title('Spectrum'), plt.xticks([]),plt.yticks([]) # plt.subplot(223), plt.imshow(spectrum_filter_log, 'gray'), plt.title('Spectrum With Filter'), plt.xticks([]),plt.yticks([]) plt.subplot(224), plt.imshow(img_new, 'gray'), plt.title('IDFT'), plt.xticks([]),plt.yticks([]) plt.tight_layout() plt.show()

# 使用陷波帶阻濾波器濾波 img_florida = cv2.imread("DIP_Figures/DIP3E_Original_Images_CH05/Fig0519(a)(florida_satellite_original).tif", -1)#-------------------------------- fft = np.fft.fft2(img_florida) fft_shift = np.fft.fftshift(fft) amp_img = np.abs(np.log(1 + np.abs(fft_shift)))#-------------------------------- BNF = narrow_notch_filter(img_florida, w=5, opening=20, vertical=True, horizontal=False)fft_NNF = np.fft.fft2(BNF*255) fft_shift_NNF = np.fft.fftshift(fft_NNF) amp_img_NNF = np.abs(np.log(1 + np.abs(fft_shift_NNF)))#-------------------------------- f1shift = fft_shift * (BNF) f2shift = np.fft.ifftshift(f1shift) #對新的進(jìn)行逆變換 img_new = np.fft.ifft2(f2shift)#出來的是復(fù)數(shù)，無法顯示 img_new = np.abs(img_new)#調(diào)整大小范圍便于顯示 img_new = (img_new-np.amin(img_new))/(np.amax(img_new)-np.amin(img_new))fft_mask = amp_img * BNFplt.figure(figsize=(15, 16)) plt.subplot(221),plt.imshow(img_florida,'gray'),plt.title('Image with noise') plt.subplot(222),plt.imshow(amp_img,'gray'),plt.title('FFT') plt.subplot(223),plt.imshow(fft_mask,'gray'),plt.title('FFT with mask') plt.subplot(224),plt.imshow(img_new,'gray'),plt.title('Denoising') plt.tight_layout() plt.show()

最優(yōu)陷波濾波

這種濾波方法的過程如下：
首先分離干擾模式的各個主要貢獻(xiàn)，然后從被污染圖像中減去該模式的一個可變加權(quán)部分。

首先提取干模式的主頻率分量，提取方法是在每個尖峰位置放一個陷波帶通濾波器傳遞函數(shù)$H_{NP}(u,v)$，則干擾噪聲模式的傅里葉變換為：
$$N(u,v)=H_{NP}(u,v)G(u,v)\text{\hspace{1em}(5.38)}$$

則有噪聲模式：
$$\eta (x,y)={\mathfrak{J}}^{?}1\{H_{NP}(u,v)G(u,v)\}\text{\hspace{1em}(5.39)}$$

如果我們知道了噪聲模式，我們假設(shè)噪聲是加性噪聲，只可以用污染的噪聲$g(x,y)$減去噪聲模式$\eta (x,y)$可得到$\hat{f}(x,y)$，但通常這只是一個近似值。
$$\hat{f}(x,y)=g(x,y)?w(x,y)\eta (x,y)\text{\hspace{1em}(5.40)}$$
$w(x,y)$是一個加權(quán)函數(shù)或調(diào)制函數(shù)，這個方法的目的就是選取$w(x,y)$，以便以某種意義的方式來優(yōu)化結(jié)果。一種方法是選擇$w(x,y)$，使$\hat{f}(x,y)$在每點$(x,y)$的規(guī)定鄰域上的方差最小。

$m\times n$(奇數(shù))的鄰域$S_{xy}$。$\hat{f}(x,y)$的“局部”方差估計如下：
$${\sigma }^2(x,y)=\frac{1}{mn}\sum _{(r,c)\in S_{xy}}[\hat{f}(r,c)?\bar{\hat{f}}{]}^2\text{\hspace{1em}(5.41)}$$

$\bar{\hat{f}}$是鄰域$\hat{f}$的平均值，
$$\bar{\hat{f}}=\frac{1}{mn}\sum _{(r,c)\in S_{xy}}\hat{f}(r,c)\text{\hspace{1em}(5.42)}$$

將式(5.40)代入(5.41)，得
$${\sigma }^2(x,y)=\frac{1}{mn}\sum _{(r,c)\in S_{xy}}\left\{\right[g(r,c)?w(r,c)\eta (r,c)]?[\overline{g}?\overline{w\eta }]{\}}^2\text{\hspace{1em}(5.43)}$$

$\overline{g}$和$\overline{w\eta }$分別是$g$和$w\eta$在鄰域$S_{xy}$的平均值

若假設(shè)$w$在$S_{xy}$內(nèi)近似為常數(shù)，則可用該鄰域中心的$w$值來代替$w(r,c)$:

$$w(r,c)=w(x,y)\text{\hspace{1em}(5.44)}$$

因為$w(x,y)$在$S_{xy}$中被假設(shè)為常數(shù)，因此在$S_{xy}$中根據(jù)$\overline{w}=w(x,y)$有

$$\overline{w\eta }=w(x,y)\overline{\eta }\text{\hspace{1em}(5.45)}$$

$\overline{\eta }$是鄰域$S_{xy}$中的平均值，所以式(5.43)變?yōu)?#xff1a;

$${\sigma }^2(x,y)=\frac{1}{mn}\sum _{(r,c)\in S_{xy}}\left\{\right[g(r,c)?w(x,y)\eta (r,c)]?[\overline{g}?w(x,y)\overline{\eta }]{\}}^2\text{\hspace{1em}(5.44)}$$

要使得${\sigma }^2(x,y)$相對$w(x,y)$最小，我們可以對式(5.44)求關(guān)于$w(x,y)$的偏導(dǎo)數(shù)，并令為偏導(dǎo)數(shù)為0；

$$\frac{?{\sigma }^2(x,y)}{?w(x,y)}=0\text{\hspace{1em}(5.47)}$$

求得$w(x,y)$:
$$w(x,y)=\frac{\overline{g\eta }?\bar{g}\bar{\eta }}{\overline{{\eta }^2}?{\bar{\eta }}^2}\text{\hspace{1em}(5.48)}$$

把式(5.48)代入式(5.40)并在噪聲圖像$g$中的每個點執(zhí)行這一過程，可得到完全復(fù)原的圖像。

# 這里還沒有實現(xiàn)，遲點再弄吧 img_mariner = cv2.imread("DIP_Figures/DIP3E_Original_Images_CH05/Fig0520(a)(NASA_Mariner6_Mars).tif", 0) M, N = img_mariner.shape[:2]fp = pad_image(img_mariner, mode='reflect') fp_cen = centralized_2d(fp) fft = np.fft.fft2(fp_cen)# 中心化后的頻譜 spectrum_fshift = spectrum_fft(fft) spectrum_log = np.log(1 + spectrum_fshift)# 未中心化的頻譜 fft_fp = np.fft.fft2(fp) spectrum_fp = spectrum_fft(fft_fp) spectrum_fp_log = np.log(1 + spectrum_fp)# 濾波器 n = 15 r = 20 H = butterworth_band_resistant_filter(fp, fp.shape, radius=40, w=5, n=5) # BNRF_1 = butterworth_notch_resistant_filter(fp, radius=r, uk=355, vk=0, n=n) # BNRF_2 = butterworth_notch_resistant_filter(fp, radius=r, uk=0, vk=355, n=n) # BNRF_3 = butterworth_notch_resistant_filter(fp, radius=r, uk=250, vk=250, n=n) # BNRF_4 = butterworth_notch_resistant_filter(fp, radius=r, uk=-250, vk=250, n=n) # BNRF = BNRF_1 * BNRF_2 * BNRF_3 * BNRF_4 * Hfft_filter = fft_fp * (1 - H) ifft = np.fft.ifft2(fft_filter) img_new = ifft.real[:M, :N]# # show = spectrum_fp_log * H # 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)plt.figure(figsize=(15, 15)) plt.subplot(221), plt.imshow(img_mariner, 'gray'), plt.title('With noise'), plt.xticks([]),plt.yticks([]) plt.subplot(222), plt.imshow(spectrum_log, 'gray'), plt.title('Spectrum Centralied'), plt.xticks([]),plt.yticks([]) plt.subplot(223), plt.imshow(spectrum_fp_log, 'gray'), plt.title('Spectrum Not Centralized'), plt.xticks([]),plt.yticks([]) plt.subplot(224), plt.imshow(img_new, 'gray'), plt.title('IDFT'), plt.xticks([]),plt.yticks([]) plt.tight_layout() plt.show()

# 巴特沃斯帶阻陷波濾波器 BNRF img_dst = img_mariner - img_new plt.figure(figsize=(16, 16)) plt.subplot(221), plt.imshow(img_dst, 'gray'), plt.title('BNF_1') # plt.subplot(222), plt.imshow(BNF_2, 'gray'), plt.title('BNF_2') # plt.subplot(223), plt.imshow(BNF_3, 'gray'), plt.title('BNF_3') # plt.subplot(224), plt.imshow(BNF_dst, 'gray'), plt.title('BNF_dst') plt.tight_layout() plt.show()

總結(jié)

以上是生活随笔為你收集整理的第5章 Python 数字图像处理(DIP) - 图像复原与重建12 - 空间滤波 - 使用频率域滤波降低周期噪声 - 陷波滤波、最优陷波滤波的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

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