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

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

陷波濾波深入介紹

零相移濾波器必須關于原點(頻率矩形中心)對稱，中以為$(u_0,v_0)$的陷波濾波器傳遞函數在$(?u_0,?v_0)$位置必須有一個對應的陷波。陷波帶阻濾波器傳遞函數可用中心被平移到陷波濾波中心的高通濾波器函數的乘積來產生

$$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)}$$

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

陷波帶通濾波器傳遞函數可用陷波帶阻濾波器得到
$$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)# 對頻譜做對數變換
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()

最優陷波濾波

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

首先提取干模式的主頻率分量，提取方法是在每個尖峰位置放一個陷波帶通濾波器傳遞函數$H_{NP}(u,v)$，則干擾噪聲模式的傅里葉變換為：
$$N(u,v)=H_{NP}(u,v)G(u,v)\text{\hspace{1em}(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)$減去噪聲模式$\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)$是一個加權函數或調制函數，這個方法的目的就是選取$w(x,y)$，以便以某種意義的方式來優化結果。一種方法是選擇$w(x,y)$，使$\hat{f}(x,y)$在每點$(x,y)$的規定鄰域上的方差最小。

$m\times n$(奇數)的鄰域$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}$的平均值

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

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

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

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

$\overline{\eta }$是鄰域$S_{xy}$中的平均值，所以式(5.43)變為：

$${\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)求關于$w(x,y)$的偏導數，并令為偏導數為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$中的每個點執行這一過程，可得到完全復原的圖像。

# 這里還沒有實現，遲點再弄吧
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()