第4章 Python 数字图像处理(DIP) - 频率域滤波10 - 使用低通频率域滤波器平滑图像 - 理想、高斯、巴特沃斯低通滤波器

使用低通频率域滤波器平滑图像

理想低通滤波器(ILPF)

在以原点为中心的一个圆内无衰减地通过所有频率，而在这个圆外“截止”所有的频率的二维低通滤波器。

$$H(u,v)=\{\begin{array}{ll}1,& D(u,v)\le D_0\\ 0,& D(u,v)>D_0\end{array}\text{\hspace{1em}(4.111)}$$

$$D(u,v)=[(u-P/2{)}^2+(v-Q/2{)}^2{]}^{1/2}\text{\hspace{1em}(4.12)}$$

$D_0$是一个正常数，控制圆的大小

def idea_low_pass_filter(source, center, radius=5):
    """
    create idea low pass filter 
    param: source: input, source image
    param: 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 image
    param: radius: input, the radius of the lowest value, greater value, bigger blocker out range, if the radius is 0, then all
                   value is 0
    return 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 = radius
    kernel = D.copy()
    
    kernel[D > D0] = 0
    kernel[D <= D0] = 1
    
    return kernel

def plot_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([])

# 理想低通滤波器 ILPF
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
from matplotlib import pyplot as plt
from matplotlib import cm

center = 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()

总图像能量$P_T$是对填充零后图像的功率谱的各个分量在点$(u,v)$处求和得到。
$$P_T=\sum _{u=0}^{P-1}\sum _{v=0}^{Q-1}P(u,v)\text{\hspace{1em}(4.113)}$$

半径为$D_0$的圆将包含$a$%的功率
$$a=100[\sum _{u=0}\sum _{v=0}P(u,v)/P_T]\text{\hspace{1em}(4.114)}$$

def mask_ring(img_ori, d=10):
    ILPF_1 = idea_low_pass_filter(img_ori, img_ori.shape, radius=d)
    ILPF_2 = idea_low_pass_filter(img_ori, img_ori.shape, radius=d-1)
    ILPF = ILPF_1 - ILPF_2
    return ILPF

# 测试模型
img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH04/Fig0441(a)(characters_test_pattern).tif', -1)
M, N = img_ori.shape[:2]

# 填充
fp = pad_image(img_ori, mode='reflect')

radius = [10, 30, 60, 160, 460]
mask = np.zeros(fp.shape)
for i in range(len(radius)):
    mask += mask_ring(fp, d=radius[i])

plt.figure(figsize=(16, 8))
plt.subplot(1, 2, 1), plt.imshow(img_ori, cmap='gray'), plt.xticks([]), plt.yticks([])
plt.subplot(1, 2, 2), plt.imshow(mask, cmap='gray'), plt.xticks([]), plt.yticks([])
plt.tight_layout()
plt.show()

注：功率计算不太正确

#功率计算不太正确
# 填充
fp = pad_image(img_ori, mode='constant')
# 中心化
fp_cen = centralized_2d(fp)
# 正变换
fft = np.fft.fft2(fp_cen)
spectrum = spectrum_fft(fft)
# spectrum = np.log(1 + spectrum)
PT = spectrum.sum()

ILPF = idea_low_pass_filter(fp, fp.shape, D0=10)

fh = fft * ILPF
spectrum = spectrum_fft(fh)
# spectrum = np.log(1 + spectrum)
P = spectrum.sum()
print(f"Power is -> {100*(P / PT)}")

Power is -> 3.269204251436661

def ilpf_test(img_ori, mode='constant', radius=10):
    M, N = img_ori.shape[:2]
    # 填充
    fp = pad_image(img_ori, mode='reflect')
    # 中心化
    fp_cen = centralized_2d(fp)
    # 正变换
    fft = np.fft.fft2(fp_cen)
    # 滤波器
    H = idea_low_pass_filter(fp, center=fp.shape, radius=radius)
    # 滤波
    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 in range(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 = ilpf_test(img_ori, radius=radius[i-1])
        ax.imshow(img, 'gray'), ax.set_title("radius = " + str(radius[i-1])), ax.set_xticks([]), ax.set_yticks([])
plt.tight_layout()
plt.show()

# 频率域ILPF传递函数对应的空间核函数
img_temp = np.zeros([1000, 1000])
ILPF = idea_low_pass_filter(img_temp, img_temp.shape, radius=15)
ifft = np.fft.ifft2(ILPF)
ifft = np.fft.ifftshift(ifft)
space = ifft.real * 1200
space_s = abs(space)
# space_s = np.clip(space, 0, space.max())
space_s = normalize(space_s)

hx = space[:, 500]
hx = centralized_2d(hx.reshape(1, -1)).flatten()

fig = plt.figure(figsize=(15, 5))
ax_1 = fig.add_subplot(1, 3, 1)
ax_1.imshow(ILPF, 'gray'), ax_1.set_xticks([]), ax_1.set_yticks([])

ax_2 = fig.add_subplot(1, 3, 2)
ax_2.imshow(space_s, 'gray'), ax_2.set_xticks([]), ax_2.set_yticks([])

ax_3 = fig.add_subplot(1, 3, 3)
ax_3.plot(hx), ax_3.set_xticks([]), ax_3.set_yticks([])

plt.tight_layout()
plt.show()

高斯低通滤波器(GLPF)

$$H(u,v)=e^{-D^2(u,v)/2D_0^2}\text{\hspace{1em}(4.116)}$$
$D_0$是截止频率。当$D(u,v)=D_0$时，GLPF传递函数下降到其最大值1.0的0.607。

def gauss_low_pass_filter(source, center, radius=5):
    """
    create gauss low pass filter 
    param: source: input, source image
    param: 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 image
    param: radius: input, the radius of the lowest value, greater value, bigger blocker out range, if the radius is 0, then all
                   value is 0
    return 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 = radius
    kernel = np.exp(- (D**2)/(2*D0**2))
    
    return kernel

def plot_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([])

# 高斯低通滤波器 GLPF
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
from matplotlib import pyplot as plt
from matplotlib import cm

center = img_ori.shape

GLPF_10 = gauss_low_pass_filter(img_ori, center, radius=10)
h_10 = GLPF_10[img_ori.shape[0]//2:, img_ori.shape[1]//2]
GLPF_20 = gauss_low_pass_filter(img_ori, center, radius=20)
h_20 = GLPF_20[img_ori.shape[0]//2:, img_ori.shape[1]//2]
GLPF_40 = gauss_low_pass_filter(img_ori, center, radius=40)
h_40 = GLPF_40[img_ori.shape[0]//2:, img_ori.shape[1]//2]
GLPF_60 = gauss_low_pass_filter(img_ori, center, radius=60)
h_60 = GLPF_60[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, GLPF_60)

ax_2 = fig.add_subplot(1, 3, 2)
ax_2.imshow(GLPF_60,'gray'), ax_2.set_xticks([]), ax_2.set_yticks([])

ax_3 = fig.add_subplot(1, 3, 3)
ax_3.plot(h_10, label='$D_0=10$'), 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_20, label='$D_0=20$')
ax_3.plot(h_40, label='$D_0=40$')
ax_3.plot(h_60, label='$D_0=60$')
plt.legend(loc='best')
plt.tight_layout()
plt.show()

def glpf_test(img_ori, mode='constant', radius=10):
    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 = gauss_low_pass_filter(fp, center=fp.shape, radius=radius)
    # 滤波
    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 in range(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 = glpf_test(img_ori, mode='reflect', radius=radius[i-1])
        ax.imshow(img, 'gray'), ax.set_title("radius = " + str(radius[i-1])), ax.set_xticks([]), ax.set_yticks([])
plt.tight_layout()
plt.show()

# 频率域GLPF传递函数对应的空间核函数
img_temp = np.zeros([1000, 1000])
GLPF = gauss_low_pass_filter(img_temp, img_temp.shape, radius=15)
ifft = np.fft.ifft2(GLPF)
ifft = np.fft.ifftshift(ifft)
space = ifft.real * 1200
space_s = abs(space)
space_s = normalize(space_s)

hx = space[:, 500]
hx = centralized_2d(hx.reshape(1, -1)).flatten()

fig = plt.figure(figsize=(15, 5))
ax_1 = fig.add_subplot(1, 3, 1)
ax_1.imshow(GLPF, 'gray'), ax_1.set_xticks([]), ax_1.set_yticks([])

ax_2 = fig.add_subplot(1, 3, 2)
ax_2.imshow(space_s, 'gray'), ax_2.set_xticks([]), ax_2.set_yticks([])

ax_3 = fig.add_subplot(1, 3, 3)
ax_3.plot(hx), ax_3.set_xticks([]), ax_3.set_yticks([])

plt.tight_layout()
plt.show()

# 不使用传统方法
import cv2
import numpy as np
import matplotlib.pyplot as plt
 
img_ic = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH04/Fig0429(a)(blown_ic).tif', 0) #直接读为灰度图像

plt.figure(figsize=(15, 12))
plt.subplot(221),plt.imshow(img_ic,'gray'),plt.title('origial'), plt.xticks([]), plt.yticks([])

#--------------------------------
fft = np.fft.fft2(img_ic)
fft_shift = np.fft.fftshift(fft)
amp_img = np.abs(np.log(1 + np.abs(fft_shift)))
plt.subplot(222),plt.imshow(amp_img,'gray'),plt.title('IC FFT'), plt.xticks([]), plt.yticks([])

#--------------------------------
glpf = gauss_low_pass_filter(img_ic, img_ic.shape, radius=20)
plt.subplot(223),plt.imshow(glpf,'gray'),plt.title('mask'), plt.xticks([]), plt.yticks([])

#--------------------------------
f1shift = fft_shift * glpf
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))
plt.subplot(224),plt.imshow(img_new,'gray'),plt.title('GLPF'), plt.xticks([]), plt.yticks([])

plt.tight_layout()
plt.show()

巴特沃斯低通滤波器

$$H(u,v)=\frac{1}{1+[D(u,v)/D_0{]}^{2n}}\text{\hspace{1em}(4.117)}$$
$$D(u,v)=[(u-M/2{)}^2+(v-N/2{)}^2{]}^{1/2}$$

• 特点
  • 较高的$n$值来控制这个BLPF函数可逼近ILPF的特性
  • 较低的$n$值来控制这个BLPF函数可逼近GLPF的特性，同时提供从低频到高频的平滑过渡。
  • 可用BLPF以小得多的振铃效应来逼近ILPF函数的清晰度

def butterworth_low_pass_filter(img, center, radius=5, n=1):
    """
    create butterworth low pass filter 
    param: source: input, source image
    param: 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 image
    param: radius: input, the radius of the lowest value, greater value, bigger blocker out range, if the radius is 0, then all
                   value is 0
    param: n: input, float, the order of the filter, if n is small, then the BLPF will be close to GLPF, and more smooth from low
              frequency to high freqency.if n is large, will close to ILPF
    return a [0, 1] value filter
    """  
    epsilon = 1e-8
    M, 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 = radius
    kernel = (1 / (1 + (D / (D0 + epsilon))**(2*n)))
    
    return kernel

def plot_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([])

# 巴特沃斯低通滤波器 BLPF
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
from matplotlib import pyplot as plt
from matplotlib import cm

center = img_ori.shape

BLPF_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()

def blpf_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 in range(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()

空间域的一阶巴特沃斯没有振铃效应。在2阶和3阶滤波器中，振铃效应通常难以察觉，但更高阶滤波器中的振铃效应很明显。

# 频率域GLPF传递函数对应的空间核函数
img_temp = np.zeros([1000, 1000])
BLPF = butterworth_low_pass_filter(img_temp, img_temp.shape, radius=15, n=25)
ifft = np.fft.ifft2(BLPF)
ifft = np.fft.ifftshift(ifft)
space = ifft.real * 1200
space_s = abs(space)
space_s = normalize(space_s)

hx = space[:, 500]
hx = centralized_2d(hx.reshape(1, -1)).flatten()

fig = plt.figure(figsize=(15, 5))
ax_1 = fig.add_subplot(1, 3, 1)
ax_1.imshow(GLPF, 'gray'), ax_1.set_xticks([]), ax_1.set_yticks([])

ax_2 = fig.add_subplot(1, 3, 2)
ax_2.imshow(space_s, 'gray'), ax_2.set_xticks([]), ax_2.set_yticks([])

ax_3 = fig.add_subplot(1, 3, 3)
ax_3.plot(hx), ax_3.set_xticks([]), ax_3.set_yticks([])

plt.tight_layout()
plt.show()

def frequen2spatial(filter):
    ifft = np.fft.ifft2(filter)
    ifft = np.fft.ifftshift(ifft)
    spatial = ifft.real * 1200
    spatial_s = abs(spatial)
    spatial_s = normalize(spatial_s)
    return spatial, spatial_s

# 频率域BLPF传递函数对应的空间核函数
img_temp = np.zeros([1000, 1000])
n = [1, 2, 5, 20]

fig = plt.figure(figsize=(20, 10))

for i in range(len(n)):
    # 这是显示空间域的核
    ax = fig.add_subplot(2, 4, i+1)
    BLPF = butterworth_low_pass_filter(img_temp, img_temp.shape, radius=15, n=n[i])
    spatial, spatial_s = frequen2spatial(BLPF)
    ax.imshow(spatial_s, 'gray'), ax.set_xticks([]), ax.set_yticks([])
    
    # 这里显示是对应的空间域核水平扫描线的灰度分布
    ax = fig.add_subplot(2, 4, i+5)
    hx = spatial[:, 500]
    hx = centralized_2d(hx.reshape(1, -1)).flatten()
    ax.plot(hx, label=f'n = {n[i]}'), ax.set_xticks([]), ax.set_yticks([])
    ax.legend(loc='best', fontsize=14)
plt.tight_layout()
plt.show()

低通滤波的例子

出下图，我们可以清晰看到不同的截止频率的核对图像的平滑效果。我们杺选择合适的核，以平滑图像，再加上其它图像处理技术，以达到想要的效果。
如下文字处理的例子，我们可以利用$D_0=20$，来平滑图像，再经过阈值处理，可以得到文字的蒙板。

# 高斯低通滤波器在印刷和出版业的应用
img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH04/Fig0419(a)(text_gaps_of_1_and_2_pixels).tif', -1)

radius = [10, 30, 60, 90, 120]

fig = plt.figure(figsize=(17, 10))
for i in range(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 = glpf_test(img_ori, mode='reflect', radius=radius[i-1])
        ax.imshow(img, 'gray'), ax.set_title("radius = " + str(radius[i-1])), ax.set_xticks([]), ax.set_yticks([])
plt.tight_layout()
plt.show()

# 高斯低通滤波器在印刷和出版业的应用
img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH04/Fig0419(a)(text_gaps_of_1_and_2_pixels).tif', -1)

radius = [10, 30, 60, 90, 120]

fig = plt.figure(figsize=(17, 10))
ax = fig.add_subplot(1, 3, 1)
ax.imshow(img_ori, 'gray'), ax.set_title('Original'), ax.set_xticks([]), ax.set_yticks([])
ax = fig.add_subplot(1, 3, 2)
img = glpf_test(img_ori, mode='reflect', radius=20)
ax.imshow(img, 'gray'), ax.set_title("radius = " + str(20)), ax.set_xticks([]), ax.set_yticks([])
ax = fig.add_subplot(1, 3, 3)
ret, img_thred = cv2.threshold(img, 0, 255, cv2.THRESH_OTSU + cv2.THRESH_BINARY_INV)
ax.imshow(img_thred, 'gray'), ax.set_title("Thred"), ax.set_xticks([]), ax.set_yticks([])
plt.tight_layout()
plt.show()

# 高斯低通滤波器在印刷和出版业的应用，平滑后的图像看上去更柔和、更美观
img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH04/Fig0427(a)(woman).tif', -1)

radius = [10, 50, 80, 130, 150]

fig = plt.figure(figsize=(17, 10))
for i in range(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 = glpf_test(img_ori, mode='reflect', radius=radius[i-1])
        ax.imshow(img, 'gray'), ax.set_title("radius = " + str(radius[i-1])), ax.set_xticks([]), ax.set_yticks([])
plt.tight_layout()
plt.show()

# 高斯低通滤波器在卫星图像的应用，这里对图像的滤波的目的是尺可能模糊更多的细节，而保留可识别的大特征。
img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH04/Fig0451(a)(satellite_original).tif', -1)

radius = [10, 20, 50, 80, 100]

fig = plt.figure(figsize=(17, 10))
for i in range(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 = glpf_test(img_ori, mode='reflect', radius=radius[i-1])
        ax.imshow(img, 'gray'), ax.set_title("radius = " + str(radius[i-1])), ax.set_xticks([]), ax.set_yticks([])
plt.tight_layout()
plt.show()