選擇性濾波

處理特定的頻帶的濾波器稱為頻帶濾波器

• 帶阻濾波器：

  • 若某個頻帶中的頻率被濾除

• 帶通濾波器：

  • 若某個頻帶中的頻率被通過

處理小頻率矩形區域的濾波器稱為陷波濾波器

• 陷波帶阻濾波器：

  • 若某個頻帶中的頻率被拒絕

• 陷波帶通濾波器：

  • 若某個頻帶中的頻率被通過

帶阻濾波器和帶通濾波器

頻率域中的帶通和帶阻濾波器傳遞函數，可通過組合低通和高通濾波器傳遞函數來構建。然后高通濾波器也是由低通濾波器推導而來，所以說低通濾波器傳遞函數是形成高通、帶阻、帶通濾波器傳遞函數的基礎。

可以由帶阻濾波器傳遞函數獲得帶通濾波器傳遞函數

$$H_{BP}(u,v)=1?H_{BR}(u,v)\text{\hspace{1em}(4.148)}$$

高斯帶阻濾波器傳遞函數

$$H(u,v)=1?e^{?\left[\frac{(D(u,v)?C_0{)}^2}{W^2}\right]}\text{\hspace{1em}(4.149)}$$

低于$C_0$時，該函數表現為一個低通高斯函數；等于$C_0$時，始終為0；高于$C_0$時，表現為一個高通高斯函數。但該函數在原點關不總是1。可以修改為：

$$H(u,v)=1?e^{?[\frac{D^2(u,v)?C_0^2}{D(u,v)W}{]}^2}\text{\hspace{1em}(4.150)}$$

巴特沃斯帶阻濾波器傳遞函數

$$H(u,v)=\frac{1}{1+[\frac{D(u,v)W}{D^2(u,v)?C_0^2}{]}^{2n}}$$

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

陷波濾波器

陷波濾波器是最有用的選擇性濾波

零相移濾波器必須關于原點(頻率矩形中心)對稱，中以為$(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}(4.151)}$$

每個濾波器的距離計算公式為
$$D_k(u,v)=[(u?M/2?u_k{)}^2+(v?N/2?v_k{)}^2{]}^{1/2}\text{\hspace{1em}(4.152)}$$
$$D_{?k}(u,v)=[(u?M/2+u_k{)}^2+(v?N/2+v_k{)}^2{]}^{1/2}\text{\hspace{1em}(4.153)}$$

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

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

陷波帶通濾波器傳遞函數可用陷波帶阻濾波器得到
$$H_{NP}(u,v)=1?H_{NR}(u,v)\text{\hspace{1em}(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()