网站公司销售,电子购物网站,哪个浏览器可以做网站,地方网站发展标题使用频率域滤波降低周期噪声陷波滤波深入介绍最优陷波滤波本章陷波滤波器有部分得出的结果不佳#xff0c;如果有更好的解决方案#xff0c;请赐教#xff0c;不胜感激。
使用频率域滤波降低周期噪声
陷波滤波深入介绍
零相移滤波器必须关于原点(频率矩形中心)对称如果有更好的解决方案请赐教不胜感激。
使用频率域滤波降低周期噪声
陷波滤波深入介绍
零相移滤波器必须关于原点(频率矩形中心)对称中以为(u0,v0)(u_0, v_0)(u0,v0)的陷波滤波器传递函数在(−u0,−v0)(-u_0, -v_0)(−u0,−v0)位置必须有一个对应的陷波。陷波带阻滤波器传递函数可用中心被平移到陷波滤波中心的高通滤波器函数的乘积来产生
HNR(u,v)∏k1QHk(u,v)H−k(u,v)(5.33)H_{NR}(u, v) \prod_{k1}^Q H_k(u, v) H_{-k}(u, v) \tag{5.33}HNR(u,v)k1∏QHk(u,v)H−k(u,v)(5.33)
每个滤波器的距离计算公式为 Dk(u,v)[(u−M/2−uk)2(v−N/2−vk)2]1/2(5.34)D_{k}(u, v) \big[(u - M / 2 - u_{k})^2 (v - N / 2 - v_{k})^2 \big]^{1/2} \tag{5.34}Dk(u,v)[(u−M/2−uk)2(v−N/2−vk)2]1/2(5.34) D−k(u,v)[(u−M/2uk)2(v−N/2vk)2]1/2(5.35)D_{-k}(u, v) \big[(u - M / 2 u_{k})^2 (v - N / 2 v_{k})^2 \big]^{1/2} \tag{5.35}D−k(u,v)[(u−M/2uk)2(v−N/2vk)2]1/2(5.35)
3个nnn阶巴特沃斯带阻滤波器 HNR(u,v)∏k13[11[D0k/Dk(u,v)]n][11[D0k/D−k(u,v)]n](5.36)H_{NR}(u, v) \prod_{k1}^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{5.36}HNR(u,v)k1∏3[1[D0k/Dk(u,v)]n1][1[D0k/D−k(u,v)]n1](5.36)
常数D0kD_{0k}D0k对每对陷波是相同的但对不同的陷波对它可以不同。
陷波带通滤波器传递函数可用陷波带阻滤波器得到 HNP(u,v)1−HNR(u,v)(5.37)H_{NP}(u, v) 1 - H_{NR}(u, v) \tag{5.37}HNP(u,v)1−HNR(u,v)(5.37)
def butterworth_notch_resistant_filter(img, uk, vk, radius10, n1):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 / (DK1e-5))**n)) * (1 / (1 (D0 / (D_K1e-5))**n))return kerneldef idea_notch_resistant_filter(source, uk, vk, radius5):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 filterM, 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 kerneldef gauss_notch_resistant_filter(source, uk, vk, radius5):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 kerneldef plot_3d(ax, x, y, z, cmap):ax.plot_surface(x, y, z, antialiasedTrue, shadeTrue, cmapcmap)ax.view_init(20, -20), ax.grid(bFalse), 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, radius20, uk30, vk80)
GNRF gauss_notch_resistant_filter(img_temp, radius20, uk30, vk80)
BNRF butterworth_notch_resistant_filter(img_temp, radius20, uk30, vk80, n5)# 用来绘制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, projection3d)
plot_3d(ax_1, u, v, INRF, cmapcm.plasma)ax_1 fig.add_subplot(1, 3, 2, projection3d)
plot_3d(ax_1, u, v, GNRF, cmapcm.PiYG)ax_1 fig.add_subplot(1, 3, 3, projection3d)
plot_3d(ax_1, u, v, BNRF, cmapcm.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 3x3dst np.ones(img.shape)dst[1::2, ::2] -1dst[::2, 1::2] -1return dst#中心化mask get_1_minus1(arr)dst arr * maskreturn dstdef pad_image(img, modeconstant):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 valuesdst np.pad(img, ((0, img.shape[0]), (0, img.shape[1])), modemode)return dst def add_sin_noise(img, scale1, angle0):add sin noise for imageparam: img: input image, 1 channel, dtypeuint8param: 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 toheight, 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)), borderValue0)# 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_widthwidth]output_img normalize(img_dst)return output_img def spectrum_fft(fft):return FFT spectrumreturn 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, scale0.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)# 对频谱做对数变换
spectrum_log np.log(1 spectrum_fshift)BNRF butterworth_notch_resistant_filter(img_ori, radius5, uk25, vk10, n4)f1shift fshift * (BNRF)
f2shift np.fft.ifftshift(f1shift) #对新的进行逆变换
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, width5,length_includes_headTrue, shapefull)
plt.arrow(285, 265, -25, -30, width5,length_includes_headTrue, shapefull)plt.subplot(223), plt.imshow(BNRF, gray), plt.title(Spectrum), plt.xticks([]),plt.yticks([])
# 在图像上加上箭头
plt.arrow(180, 180, 25, 30, width5,length_includes_headTrue, shapefull)
plt.arrow(285, 265, -25, -30, width5,length_includes_headTrue, shapefull)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, scale0.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)# 对频谱做对数变换
spectrum_log np.log(1 spectrum_fshift)BNRF 1 - butterworth_notch_resistant_filter(img_ori, radius5, uk25, vk10, n4)f1shift fshift * (BNRF)
f2shift np.fft.ifftshift(f1shift) #对新的进行逆变换
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, width5,length_includes_headTrue, shapefull)
# plt.arrow(285, 265, -25, -30, width5,length_includes_headTrue, shapefull)# plt.subplot(223), plt.imshow(BNRF, gray), plt.title(Spectrum), plt.xticks([]),plt.yticks([])
# # 在图像上加上箭头
# plt.arrow(180, 180, 25, 30, width5,length_includes_headTrue, shapefull)
# plt.arrow(285, 265, -25, -30, width5,length_includes_headTrue, shapefull)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, radius10, w5, n1):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, radius5, n1):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, modereflect)
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, radius100, n4)
# BNRF_1 butterworth_notch_resistant_filter(fp, radiusr, uk355, vk0, nn)
# BNRF_2 butterworth_notch_resistant_filter(fp, radiusr, uk0, vk355, nn)
# BNRF_3 butterworth_notch_resistant_filter(fp, radiusr, uk250, vk250, nn)
# BNRF_4 butterworth_notch_resistant_filter(fp, radiusr, uk-250, vk250, nn)
# 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, modeconstant)
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, radius100, n3)
# BNRF_1 butterworth_notch_resistant_filter(fp, radiusr, uk355, vk0, nn)
# BNRF_2 butterworth_notch_resistant_filter(fp, radiusr, uk0, vk355, nn)
# BNRF_3 butterworth_notch_resistant_filter(fp, radiusr, uk250, vk250, nn)
# BNRF_4 butterworth_notch_resistant_filter(fp, radiusr, uk-250, vk250, nn)
# 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, w5, opening10, verticalTrue, horizontalFalse):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 Trueparam: horizontal: input, boolean, whether horizontal or not, default is Falsereturn a [0, 1] value butterworth band resistant filter assert w 0, W must greater than 0w_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, modereflect)
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, w10, opening30, verticalTrue, horizontalFalse)
# BNRF_1 butterworth_notch_resistant_filter(fp, radiusr, uk355, vk0, nn)
# BNRF_2 butterworth_notch_resistant_filter(fp, radiusr, uk0, vk355, nn)
# BNRF_3 butterworth_notch_resistant_filter(fp, radiusr, uk250, vk250, nn)
# BNRF_4 butterworth_notch_resistant_filter(fp, radiusr, uk-250, vk250, nn)
# 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, modereflect)
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, w10, opening30, verticalTrue, horizontalFalse)
# BNRF_1 butterworth_notch_resistant_filter(fp, radiusr, uk355, vk0, nn)
# BNRF_2 butterworth_notch_resistant_filter(fp, radiusr, uk0, vk355, nn)
# BNRF_3 butterworth_notch_resistant_filter(fp, radiusr, uk250, vk250, nn)
# BNRF_4 butterworth_notch_resistant_filter(fp, radiusr, uk-250, vk250, nn)
# 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, w5, opening20, verticalTrue, horizontalFalse)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) #对新的进行逆变换
img_new np.fft.ifft2(f2shift)#出来的是复数无法显示
img_new np.abs(img_new)#调整大小范围便于显示
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()最优陷波滤波
这种滤波方法的过程如下 首先分离干扰模式的各个主要贡献然后从被污染图像中减去该模式的一个可变加权部分。
首先提取干模式的主频率分量提取方法是在每个尖峰位置放一个陷波带通滤波器传递函数HNP(u,v)H_{NP}(u, v)HNP(u,v)则干扰噪声模式的傅里叶变换为 N(u,v)HNP(u,v)G(u,v)(5.38)N(u, v) H_{NP}(u, v)G(u, v) \tag{5.38}N(u,v)HNP(u,v)G(u,v)(5.38)
则有噪声模式 η(x,y)J−1{HNP(u,v)G(u,v)}(5.39)\eta(x, y) \mathfrak{J}^-1 \{ H_{NP}(u, v)G(u, v) \} \tag{5.39}η(x,y)J−1{HNP(u,v)G(u,v)}(5.39)
如果我们知道了噪声模式我们假设噪声是加性噪声只可以用污染的噪声g(x,y)g(x, y)g(x,y)减去噪声模式η(x,y)\eta(x, y)η(x,y)可得到f^(x,y)\hat{f}(x, y)f^(x,y)但通常这只是一个近似值。 f^(x,y)g(x,y)−w(x,y)η(x,y)(5.40)\hat{f}(x, y) g(x, y) - w(x, y)\eta(x, y) \tag{5.40}f^(x,y)g(x,y)−w(x,y)η(x,y)(5.40) w(x,y)w(x, y)w(x,y)是一个加权函数或调制函数这个方法的目的就是选取w(x,y)w(x, y)w(x,y)以便以某种意义的方式来优化结果。一种方法是选择w(x,y)w(x, y)w(x,y)使f^(x,y)\hat{f}(x, y)f^(x,y)在每点(x,y)(x, y)(x,y)的规定邻域上的方差最小。
m×nm\times{n}m×n(奇数)的邻域SxyS_{xy}Sxy。f^(x,y)\hat{f}(x, y)f^(x,y)的“局部”方差估计如下 σ2(x,y)1mn∑(r,c)∈Sxy[f^(r,c)−f^ˉ]2(5.41)\sigma^2(x, y) \frac{1}{mn} \sum_{(r,c)\in S_{xy}} \Big[ \hat{f}(r, c) - \bar{\hat{f}} \Big]^2 \tag{5.41}σ2(x,y)mn1(r,c)∈Sxy∑[f^(r,c)−f^ˉ]2(5.41)
f^ˉ\bar{\hat{f}}f^ˉ是邻域f^\hat{f}f^的平均值 f^ˉ1mn∑(r,c)∈Sxyf^(r,c)(5.42)\bar{\hat{f}} \frac{1}{mn} \sum_{(r,c)\in S_{xy}} \hat{f}(r, c) \tag{5.42}f^ˉmn1(r,c)∈Sxy∑f^(r,c)(5.42)
将式(5.40)代入(5.41)得 σ2(x,y)1mn∑(r,c)∈Sxy{[g(r,c)−w(r,c)η(r,c)]−[g‾−wη‾]}2(5.43)\sigma^2(x, y) \frac{1}{mn} \sum_{(r,c)\in S_{xy}} \Big\{\big[g(r, c) - w(r, c) \eta(r, c)\big] - \big[\overline{g} - \overline{w\eta} \big ] \Big\}^2\tag{5.43}σ2(x,y)mn1(r,c)∈Sxy∑{[g(r,c)−w(r,c)η(r,c)]−[g−wη]}2(5.43)
g‾\overline{g}g和wη‾\overline{w\eta}wη分别是ggg和wηw\etawη在邻域SxyS_{xy}Sxy的平均值
若假设www在SxyS_{xy}Sxy内近似为常数则可用该邻域中心的www值来代替w(r,c)w(r, c)w(r,c):
w(r,c)w(x,y)(5.44)w(r, c) w(x, y) \tag{5.44}w(r,c)w(x,y)(5.44)
因为w(x,y)w(x, y)w(x,y)在SxyS_{xy}Sxy中被假设为常数因此在SxyS_{xy}Sxy中根据w‾w(x,y)\overline{w} w(x, y)ww(x,y)有
wη‾w(x,y)η‾(5.45)\overline{w\eta} w(x, y) \overline{\eta} \tag{5.45}wηw(x,y)η(5.45)
η‾\overline{\eta}η是邻域SxyS_{xy}Sxy中的平均值所以式(5.43)变为
σ2(x,y)1mn∑(r,c)∈Sxy{[g(r,c)−w(x,y)η(r,c)]−[g‾−w(x,y)η‾]}2(5.44)\sigma^2(x, y) \frac{1}{mn} \sum_{(r,c)\in S_{xy}} \Big\{\big[g(r, c) - w(x, y) \eta(r, c)\big] - \big[\overline{g} - {w(x, y)}\overline{\eta} \big ] \Big\}^2\tag{5.44}σ2(x,y)mn1(r,c)∈Sxy∑{[g(r,c)−w(x,y)η(r,c)]−[g−w(x,y)η]}2(5.44)
要使得σ2(x,y)\sigma^2(x, y)σ2(x,y)相对w(x,y)w(x, y)w(x,y)最小我们可以对式(5.44)求关于w(x,y)w(x, y)w(x,y)的偏导数并令为偏导数为0
∂σ2(x,y)∂w(x,y)0(5.47)\frac{\partial{\sigma^2(x, y)}}{\partial{w(x, y)}} 0 \tag{5.47}∂w(x,y)∂σ2(x,y)0(5.47)
求得w(x,y)w(x, y)w(x,y): w(x,y)gη‾−gˉηˉη2‾−ηˉ2(5.48)w(x, y) \frac{\overline{g\eta} - \bar{g}\bar{\eta}}{\overline{\eta^2} - \bar{\eta}^2}\tag{5.48}w(x,y)η2−ηˉ2gη−gˉηˉ(5.48)
把式(5.48)代入式(5.40)并在噪声图像ggg中的每个点执行这一过程可得到完全复原的图像。
# 这里还没有实现迟点再弄吧
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, modereflect)
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, radius40, w5, n5)
# BNRF_1 butterworth_notch_resistant_filter(fp, radiusr, uk355, vk0, nn)
# BNRF_2 butterworth_notch_resistant_filter(fp, radiusr, uk0, vk355, nn)
# BNRF_3 butterworth_notch_resistant_filter(fp, radiusr, uk250, vk250, nn)
# BNRF_4 butterworth_notch_resistant_filter(fp, radiusr, uk-250, vk250, nn)
# 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()
