第4章 Python 数字图像处理(DIP) - 频率域滤波12 - 选择性滤波

选择性滤波

处理特定的频带的滤波器称为频带滤波器

带阻滤波器：
- 若某个频带中的频率被滤除
带通滤波器：
- 若某个频带中的频率被通过

处理小频率矩形区域的滤波器称为陷波滤波器

陷波带阻滤波器：
- 若某个频带中的频率被拒绝
陷波带通滤波器：
- 若某个频带中的频率被通过

带阻滤波器和带通滤波器

频率域中的带通和带阻滤波器传递函数，可通过组合低通和高通滤波器传递函数来构建。然后高通滤波器也是由低通滤波器推导而来，所以说低通滤波器传递函数是形成高通、带阻、带通滤波器传递函数的基础。

可以由带阻滤波器传递函数获得带通滤波器传递函数

HBP(u,v)=1−HBR(u,v)(4.148)H_{BP}(u, v) = 1 - H_{BR}(u, v) \tag{4.148}HBP(u,v)=1−HBR(u,v)(4.148)

高斯带阻滤波器传递函数

H(u,v)=1−e−[(D(u,v)−C0)2W2](4.149)H(u, v) = 1 - e^{-\big[\frac{(D(u, v) - C_0)^2}{W^2} \ \ \big]} \tag{4.149}H(u,v)=1−e−[W2(D(u,v)−C0)2 ](4.149)

低于C0C_0C0时，该函数表现为一个低通高斯函数；等于C0C_0C0时，始终为0；高于C0C_0C0时，表现为一个高通高斯函数。但该函数在原点关不总是1。可以修改为：

H(u,v)=1−e−[D2(u,v)−C02D(u,v)W]2(4.150)H(u, v) = 1 - e^{-\big[\frac{D^2(u, v) - C_0^2}{D(u, v) W} \ \ \big]^2} \tag{4.150}H(u,v)=1−e−[D(u,v)WD2(u,v)−C02 ]2(4.150)

巴特沃斯带阻滤波器传递函数

H(u,v)=11+[D(u,v)WD2(u,v)−C02]2nH(u, v) = \frac{1}{1 + \bigg[\frac{D(u, v) W}{D^2(u, v) - C_0^2} \bigg]^{2n}}H(u,v)=1+[D2(u,v)−C02D(u,v)W]2n1

def idea_band_resistant_filter(source, center, radius=10, w=5):"""create idea band resistant 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, 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 5return a [0, 1] value band resistant 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 = radiushalf_w = w / 2kernel_1 = D.copy()assert radius > half_w, "radius must greater than W/2"#==================piecewise================kernel = np.piecewise(kernel_1, [kernel_1 <= D0 + half_w, kernel_1 <= D0 - half_w], [1, 0])kernel = 1 - kernel#==================where==================
#     kernel = np.where(kernel_1 > D0 + half_w, 1, kernel_1)
#     kernel = np.where(kernel <= D0 - half_w, 1, kernel)
#     kernel = np.where(kernel != 1, 0, kernel)# =================公式法================
#     kernel_2 = D.copy()
#     kernel_1[D >  D0 + half_w] = 1
#     kernel_1[D <= D0 + half_w] = 0
#     kernel_2[D > D0 - half_w] = 1
#     kernel_2[D <= D0 - half_w] = 0
#     kernel = kernel_1 - kernel_2 return kernel

def gauss_band_resistant_149(source, center, radius=10, w=5):"""create gaussian band resistant filter, equation 4.149param: 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 5return a [0, 1] value gaussian band resistant filter"""    N, 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 = radiuskernel = 1 - np.exp(-(D - C0)**2 / (w**2))return kernel

def gauss_band_resistant_filter(source, center, radius=10, w=5):"""create gaussian 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 5return a [0, 1] value gaussian band resistant filter"""    N, 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 = radiuskernel = 1 - np.exp(-((D**2 - C0**2) / (D * w))**2)return kernel

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"""    N, 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)kernel = 1 / (1 + temp ** (2*n)) return kernel

# 带阻滤波器传递函数
img_temp = np.zeros([1000, 1000])
C0 = 100# 1理想带阻滤波器
IBRF = idea_band_resistant_filter(img_temp, img_temp.shape, radius=C0, w=100)
hx_i = IBRF[500:, 500].flatten()# 2由高斯低通和高斯高通滤波器函数相加形成的带阻传递函数，最小值不是0，并且与C0不重合
GHPF = gauss_high_pass_filter(img_temp, img_temp.shape, radius=C0)
GLPF = gauss_low_pass_filter(img_temp, img_temp.shape, radius=C0/2)
GLPF = GHPF + GLPF
hx_g = GLPF[500:, 500].flatten()# 3由式4.149得到的，原点处的值不是1
GBRF_149 = gauss_band_resistant_149(img_temp, img_temp.shape, radius=C0, w=100)
hx_g149 = GBRF_149[500:, 500].flatten()# 4由式4.150得到的
GBRF = gauss_band_resistant_filter(img_temp, img_temp.shape, radius=C0, w=100)
hx_gbrf = GBRF[500:, 500].flatten()fig = plt.figure(figsize=(16, 3))
ax_1 = fig.add_subplot(1, 4, 1)
ax_1.plot(hx_i), ax_1.set_yticks([0, 1.0]), ax_1.set_xticks([100, 500]), ax_1.set_ylim(0, 1.1), ax_1.set_xlim(0, 500)ax_2 = fig.add_subplot(1, 4, 2)
ax_2.plot(hx_g), ax_2.set_yticks([0.75, 1.0]), ax_2.set_xticks([100, 500]), ax_2.set_ylim(0.75, 1.1), ax_2.set_xlim(0, 500)ax_3 = fig.add_subplot(1, 4, 3)
ax_3.plot(hx_g149), ax_3.set_yticks([0, 1.0]), ax_3.set_xticks([100, 500]), ax_3.set_ylim(0, 1.1), ax_3.set_xlim(0, 500)ax_3 = fig.add_subplot(1, 4, 4)
ax_3.plot(hx_gbrf), ax_3.set_yticks([0, 1.0]), ax_3.set_xticks([100, 500]), ax_3.set_ylim(0, 1.1), ax_3.set_xlim(0, 500)plt.tight_layout()
plt.show()

# 理想、高斯、巴特沃斯带阻传递函数
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
from matplotlib import pyplot as plt
from matplotlib import cmimg_temp = np.zeros([512, 512])
center = img_temp.shaperadius = 128
w = 60
IBRF = idea_band_resistant_filter(img_temp, img_temp.shape, radius=radius, w=w)
GBFR = gauss_band_resistant_filter(img_temp, img_temp.shape, radius=radius, w=w)
BBRF = butterworth_band_resistant_filter(img_temp, img_temp.shape, radius=radius, w=w, n=1)filters = ['IBRF', 'GBFR', 'BBRF']
# 用来绘制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=(15, 15))for i in range(len(filters)):ax_1 = fig.add_subplot(3, 3, i*3 + 1, projection='3d')plot_3d(ax_1, u, v, eval(filters[i]))ax_2 = fig.add_subplot(3, 3, i*3 + 2)ax_2.imshow(eval(filters[i]),'gray'), ax_2.set_title(filters[i]), ax_2.set_xticks([]), ax_2.set_yticks([])h_1 = eval(filters[i])[img_temp.shape[0]//2:, img_temp.shape[1]//2]ax_3 = fig.add_subplot(3, 3, i*3 + 3)ax_3.plot(h_1), ax_3.set_xticks([0, radius//2]), ax_3.set_yticks([0, 1]), ax_3.set_xlim([0, 320]), ax_3.set_ylim([0, 1.2])
plt.tight_layout()
plt.show()

陷波滤波器

陷波滤波器是最有用的选择性滤波

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

每个滤波器的距离计算公式为
Dk(u,v)=[(u−M/2−uk)2+(v−N/2−vk)2]1/2(4.152)D_{k}(u, v) = \big[(u - M / 2 - u_{k})^2 + (v - N / 2 - v_{k})^2 \big]^{1/2} \tag{4.152}Dk(u,v)=[(u−M/2−uk)2+(v−N/2−vk)2]1/2(4.152)
D−k(u,v)=[(u−M/2+uk)2+(v−N/2+vk)2]1/2(4.153)D_{-k}(u, v) = \big[(u - M / 2 + u_{k})^2 + (v - N / 2 + v_{k})^2 \big]^{1/2} \tag{4.153}D−k(u,v)=[(u−M/2+uk)2+(v−N/2+vk)2]1/2(4.153)

nnn阶巴特沃斯带阴滤波器
HNR(u,v)=∏k=13[11+[D0k/Dk(u,v)]n][11+[D0k/D−k(u,v)]n](4.154)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{4.154}HNR(u,v)=k=1∏3[1+[D0k/Dk(u,v)]n1][1+[D0k/D−k(u,v)]n1](4.154)

常数D0kD_{0k}D0k对每对陷波是相同的，但对不同的陷波对，它可以不同。

陷波带通滤波器传递函数可用陷波带阻滤波器得到
HNP(u,v)=1−HNR(u,v)(4.155)H_{NP}(u, v) = 1 - H_{NR}(u, v) \tag{4.155}HNP(u,v)=1−HNR(u,v)(4.155)

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

# 巴特沃斯带阻陷波滤波器 BNRF
img_temp = np.zeros([512, 512])
BNF_1 = butterworth_notch_resistant_filter(img_temp, radius=10, uk=30, vk=40, n=3)
BNF_2 = butterworth_notch_resistant_filter(img_temp, radius=10, uk=30, vk=80, n=3)
BNF_3 = butterworth_notch_resistant_filter(img_temp, radius=10, uk=-30, vk=80, n=3)plt.figure(figsize=(16, 12))
plt.subplot(141), plt.imshow(BNF_1, 'gray'), plt.title('BNF_1')
plt.subplot(142), plt.imshow(BNF_2, 'gray'), plt.title('BNF_2')
plt.subplot(143), plt.imshow(BNF_3, 'gray'), plt.title('BNF_3')BNF_dst = BNF_1 * BNF_2 * BNF_3plt.subplot(144), plt.imshow(BNF_dst, 'gray'), plt.title('BNF_dst')plt.tight_layout()
plt.show()

# 使用陷波滤波删除数字化印刷图像中的莫尔模式
img_ori = cv2.imread("DIP_Figures/DIP3E_Original_Images_CH04/Fig0464(a)(car_75DPI_Moire).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 = spectrum_fft(fft)
# 对频谱做对数变换
spectrum_log = np.log(1 + spectrum)# 巴特沃斯陷波带阻滤波器
BNRF_1 = butterworth_notch_resistant_filter(fp, radius=9, uk=60, vk=80, n=4)
BNRF_2 = butterworth_notch_resistant_filter(fp, radius=9, uk=-60, vk=80, n=4)
BNRF_3 = butterworth_notch_resistant_filter(fp, radius=9, uk=60, vk=160, n=4)
BNRF_4 = butterworth_notch_resistant_filter(fp, radius=9, uk=-60, vk=160, n=4)BNRF = BNRF_1 * BNRF_2 * BNRF_3 * BNRF_4 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)fig = plt.figure(figsize=(10, 14))ax_1 = fig.add_subplot(2, 2, 1)
ax_1.imshow(img_ori, 'gray'), ax_1.set_title('Original'), ax_1.set_xticks([]), ax_1.set_yticks([])ax_2 = fig.add_subplot(2, 2, 2)
ax_2.imshow(spectrum_log, 'gray'), ax_2.set_title('Spectrum Before Filter'), ax_2.set_xticks([]), ax_2.set_yticks([])ax_3 = fig.add_subplot(2, 2, 3)
ax_3.imshow(spectrum_filter_log, 'gray'), ax_3.set_title('Spectrum After Filter'), ax_3.set_xticks([]), ax_3.set_yticks([])ax_4 = fig.add_subplot(2, 2, 4)
ax_4.imshow(img_new, 'gray'), ax_4.set_title('Denoising'), ax_4.set_xticks([]), ax_4.set_yticks([])plt.tight_layout()
plt.show()

使用陷波滤波去除周期干扰

def narrow_notch_filter(img, w=5, opening=10, vertical=True, horizontal=False):"""create narrow notch resistant filter, using opencvparam: 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"""      dst = np.ones(img.shape, dtype=np.uint8) * 255c_height, c_width = img.shape[0] // 2, img.shape[1] // 2if vertical:cv2.rectangle(dst, ((img.shape[1] - w)//2, 0), (c_width + w//2, img.shape[0]), (0), -1)cv2.rectangle(dst, (0, (img.shape[0] - opening)//2), (img.shape[1], c_height + opening//2), (255), -1)horizontal_ = np.ones(img.shape, dtype=np.uint8) * 255if horizontal:        cv2.rectangle(horizontal_, (0, (img.shape[0] - w)//2), (img.shape[1], c_height + w//2), (0), -1)cv2.rectangle(horizontal_, ((img.shape[1] - opening)//2, 0), (c_width + opening//2, img.shape[0]), (255), -1)dst = dst * horizontal_dst = dst / dst.max()return dst

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])N, M = 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

# NNF narrow_notch_filter
img_temp = np.zeros([512, 512])
NNF = narrow_notch_filter(img_temp, 5, 20, vertical=True, horizontal=False)
plt.figure(figsize=(10, 8))
plt.imshow(NNF,'gray'),plt.title('NNF')
plt.show()

# 使用陷波滤波去除周期干扰
img_ori = cv2.imread("DIP_Figures/DIP3E_Original_Images_CH04/Fig0465(a)(cassini).tif", 0)M, N = img_ori.shape[:2]# 填充为'constant'可得到跟书上一样的频谱，这里使用'reflect'
fp = pad_image(img_ori, mode='constant')
# 中心化
fp_cen = centralized_2d(fp)
# 正变换
fft = np.fft.fft2(fp_cen)# 频谱
spectrum = spectrum_fft(fft)
# 对频谱做对数变换
spectrum_log = np.log(1 + spectrum)# 巴特沃斯陷波带阻滤波器
NRF = narrow_notch_filter(fp, w=8, opening=20, vertical=True, horizontal=False)fft_filter = fft * NRF# 滤波后的频谱
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.uint8(normalize(img_new) * 255)fig = plt.figure(figsize=(10, 10))ax_1 = fig.add_subplot(2, 2, 1)
ax_1.imshow(img_ori, 'gray'), ax_1.set_title('Original'), ax_1.set_xticks([]), ax_1.set_yticks([])ax_2 = fig.add_subplot(2, 2, 2)
ax_2.imshow(spectrum_log, 'gray'), ax_2.set_title('Spectrum Before Filter'), ax_2.set_xticks([]), ax_2.set_yticks([])ax_3 = fig.add_subplot(2, 2, 3)
ax_3.imshow(spectrum_filter_log, 'gray'), ax_3.set_title('Spectrum After Filter'), ax_3.set_xticks([]), ax_3.set_yticks([])ax_4 = fig.add_subplot(2, 2, 4)
ax_4.imshow(img_new, 'gray'), ax_4.set_title('Denoising'), ax_4.set_xticks([]), ax_4.set_yticks([])plt.tight_layout()
plt.show()

# 周期干扰的空间模式
img_ori = cv2.imread("DIP_Figures/DIP3E_Original_Images_CH04/Fig0465(a)(cassini).tif", 0)
M, N = img_ori.shape[:2]# 填充为'constant'可得到跟书上一样的频谱，这里使用'reflect'
fp = pad_image(img_ori, mode='constant')
# 中心化
fp_cen = centralized_2d(fp)
# 正变换
fft = np.fft.fft2(fp_cen)# # 频谱
# spectrum = spectrum_fft(fft)
# # 对频谱做对数变换
# spectrum_log = np.log(1 + spectrum)# 巴特沃斯陷波带阻滤波器
NRF = narrow_notch_filter(fp, w=8, opening=20, vertical=True, horizontal=False)
NRF = 1 - NRF
fft_filter = fft * NRF# 滤波后的频谱
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.uint8(normalize(img_new) * 255)fig = plt.figure(figsize=(10, 10))# ax_1 = fig.add_subplot(2, 2, 1)
# ax_1.imshow(img_ori, 'gray'), ax_1.set_title('Original'), ax_1.set_xticks([]), ax_1.set_yticks([])# ax_2 = fig.add_subplot(2, 2, 2)
# ax_2.imshow(spectrum_log, 'gray'), ax_2.set_title('Spectrum Before Filter'), ax_2.set_xticks([]), ax_2.set_yticks([])ax_3 = fig.add_subplot(2, 2, 3)
ax_3.imshow(NRF, 'gray'), ax_3.set_title('Spectrum After Filter'), ax_3.set_xticks([]), ax_3.set_yticks([])ax_4 = fig.add_subplot(2, 2, 4)
ax_4.imshow(img_new, 'gray'), ax_4.set_title('Noise'), ax_4.set_xticks([]), ax_4.set_yticks([])plt.tight_layout()
plt.show()