第5章 Python 数字图像处理(DIP) - 图像复原与重建12 - 空间滤波 - 使用频率域滤波降低周期噪声

标题

使用频率域滤波降低周期噪声
- 陷波滤波深入介绍
- 最优陷波滤波

本章陷波滤波器有部分得出的结果不佳，如果有更好的解决方案，请赐教，不胜感激。

使用频率域滤波降低周期噪声

陷波滤波深入介绍

零相移滤波器必须关于原点(频率矩形中心)对称，中以为(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)(5.33)H_{NR}(u, v) = \prod_{k=1}^Q H_k(u, v) H_{-k}(u, v) \tag{5.33}HNR(u,v)=k=1∏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/2+uk)2+(v−N/2+vk)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/2+uk)2+(v−N/2+vk)2]1/2(5.35)

3个nnn阶巴特沃斯带阻滤波器
HNR(u,v)=∏k=13[11+[D0k/Dk(u,v)]n][11+[D0k/D−k(u,v)]n](5.36)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{5.36}HNR(u,v)=k=1∏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, 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)# 对频谱做对数变换
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) #对新的进行逆变换
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)# 对频谱做对数变换
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) #对新的进行逆变换
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) #对新的进行逆变换
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)w=w(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, 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()