defidea_low_pass_filter(source, center, radius=5):"""create idea 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 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)D = np.sqrt((u - center[1]//2)**2+(v - center[0]//2)**2)D0 = radiuskernel = D.copy()kernel[D > D0]=0kernel[D <= D0]=1return kernel
defplot_3d(ax, x, y, z):ax.plot_surface(x, y, z, antialiased=True, shade=True)ax.view_init(20,60), ax.grid(b=False), ax.set_xticks([]), ax.set_yticks([]), ax.set_zticks([])# 理想低通濾波器 ILPFfrom mpl_toolkits.mplot3d import Axes3D
import numpy as np
from matplotlib import pyplot as plt
from matplotlib import cmcenter = img_ori.shape
ILPF = idea_low_pass_filter(img_ori, center, radius=50)# 用來繪制3D圖
M, N = img_ori.shape[1], img_ori.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, ILPF)ax_2 = fig.add_subplot(1,3,2)
ax_2.imshow(ILPF,'gray'), ax_2.set_xticks([]), ax_2.set_yticks([])h = ILPF[img_ori.shape[0]//2:, img_ori.shape[1]//2]
ax_3 = fig.add_subplot(1,3,3)
ax_3.plot(h), ax_3.set_xticks([0,50]), ax_3.set_yticks([0,1]), ax_3.set_xlim([0,320]), ax_3.set_ylim([0,1.2])plt.tight_layout()
plt.show()
defbutterworth_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
defplot_3d(ax, x, y, z):ax.plot_surface(x, y, z, antialiased=True, shade=True)ax.view_init(20,60), ax.grid(b=False), ax.set_xticks([]), ax.set_yticks([]), ax.set_zticks([])# 巴特沃斯低通濾波器 BLPFfrom mpl_toolkits.mplot3d import Axes3D
import numpy as np
from matplotlib import pyplot as plt
from matplotlib import cmcenter = img_ori.shapeBLPF_60_1 = butterworth_low_pass_filter(img_ori, center, radius=60, n=1)
h_1 = BLPF_60_1[img_ori.shape[0]//2:, img_ori.shape[1]//2]
BLPF_60_2 = butterworth_low_pass_filter(img_ori, center, radius=60, n=2)
h_2 = BLPF_60_2[img_ori.shape[0]//2:, img_ori.shape[1]//2]
BLPF_60_3 = butterworth_low_pass_filter(img_ori, center, radius=60, n=3)
h_3 = BLPF_60_3[img_ori.shape[0]//2:, img_ori.shape[1]//2]
BLPF_60_4 = butterworth_low_pass_filter(img_ori, center, radius=60, n=4)
h_4 = BLPF_60_4[img_ori.shape[0]//2:, img_ori.shape[1]//2]# 用來繪制3D圖
M, N = img_ori.shape[1], img_ori.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, BLPF_60_1)ax_2 = fig.add_subplot(1,3,2)
ax_2.imshow(BLPF_60_1,'gray'), ax_2.set_xticks([]), ax_2.set_yticks([])ax_3 = fig.add_subplot(1,3,3)
ax_3.plot(h_1, label='$n=1$'), ax_3.set_xticks([0,50]), ax_3.set_yticks([0,1]), ax_3.set_xlim([0,320]), ax_3.set_ylim([0,1.2])
ax_3.plot(h_2, label='$n=2$')
ax_3.plot(h_3, label='$n=3$')
ax_3.plot(h_4, label='$n=4$')
plt.legend(loc='best')
plt.tight_layout()
plt.show()defblpf_test(img_ori, mode='constant', radius=10, n=1):M, N = img_ori.shape[:2]# 填充fp = pad_image(img_ori, mode=mode)# 中心化fp_cen = centralized_2d(fp)# 正變換fft = np.fft.fft2(fp_cen)# 濾波器H = butterworth_low_pass_filter(fp, center=fp.shape, radius=radius, n=n)# 濾波HF = fft * H# 反變換ifft = np.fft.ifft2(HF)# 去中心化gp = centralized_2d(ifft.real)# 還回來與原圖像的大小g = gp[:M,:N]dst = np.uint8(normalize(g)*255)return dst
# 巴特沃斯低通濾波器在頻率域濾波的使用
img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH04/Fig0441(a)(characters_test_pattern).tif',-1)radius =[10,30,60,160,460]fig = plt.figure(figsize=(15,10))for i inrange(len(radius)+1):ax = fig.add_subplot(2,3, i+1)if i ==0:ax.imshow(img_ori,'gray'), ax.set_title('Original'), ax.set_xticks([]), ax.set_yticks([])else:img = blpf_test(img_ori, mode='reflect', radius=radius[i-1], n=2.25)ax.imshow(img,'gray'), ax.set_title("radius = "+str(radius[i-1])), ax.set_xticks([]), ax.set_yticks([])
plt.tight_layout()
plt.show()