图像频域增强：陷波滤波器

理想陷波滤波器

陷波滤波器可以阻止或允许以某个频率为中心的邻域里的频率通过，所以其在本质上仍然是带阻滤波器或带通滤波器，并且可分别称为
带阻陷波滤波器和带通陷波滤波器。

理想带阻陷波滤波器转移函数为:

H(u,v)={0,D(u,v)<=D01,D(u,v)>D0H(u, v) = \left\{\begin{aligned}&0, D(u, v) <= D_0 \\ &1, D(u, v) > D_0 \end{aligned}\right.H(u,v)={0,D(u,v)<=D01,D(u,v)>D0

D(u,v)=[(u−u0)2+(v−v0)2]12D(u, v) = [(u - u_0)^2 + (v - v_0)^2]^{\frac{1}{2}}D(u,v)=[(u−u0)2+(v−v0)2]21

带阻陷波滤波器一般是两两工作的，有两个频率中心：

H(u,v)={0,D1(u,v)<=D0orD2(u,v)<=D01,OtherH(u, v) = \left\{\begin{aligned}&0, D_1(u, v) <= D_0 or D_2(u, v) <= D_0 \\ &1, Other \end{aligned}\right.H(u,v)={0,D1(u,v)<=D0orD2(u,v)<=D01,Other

D1(u,v)=[(u−u0)2+(v−v0)2]12D_1(u, v) = [(u - u_0)^2 + (v - v_0)^2]^{\frac{1}{2}}D1(u,v)=[(u−u0)2+(v−v0)2]21
D2(u,v)=[(u−u1)2+(v+v1)2]12D_2(u, v) = [(u - u_1)^2 + (v + v_1)^2]^{\frac{1}{2}}D2(u,v)=[(u−u1)2+(v+v1)2]21

取反后得到陷波带通滤波器。

Numpy陷波滤波：

import cv2 as cv
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt

img = cv.imread('./images/R-C.png')
img = cv.cvtColor(img, cv.COLOR_BGR2RGB)
img_gray = cv.cvtColor(img, cv.COLOR_RGB2GRAY)

测试用例:

理想陷波带通、带阻滤波：

def notch_filter(img_gray, u0=0, v0=0, d0=50, ftype='pass'):# 以频谱左上角为坐标原点dft = cv.dft(img_gray.astype('float32'),flags = cv.DFT_COMPLEX_OUTPUT) dft_shift = np.fft.fftshift(dft)                                    m, n, _ = dft_shift.shapemask = np.zeros_like(dft_shift)x_arr = np.concatenate([np.arange(m).reshape(m, 1)], axis=1)y_arr = np.concatenate([np.arange(n).reshape(1, n)], axis=0)dist1 = np.sqrt((x_arr - u0)**2 + (y_arr - v0)**2)u1 = m - u0v1 = n - v0dist2 = np.sqrt((x_arr - u1)**2 + (y_arr - v1)**2)mask[dist1 <= d0] = 1mask[dist2 <= d0] = 1if ftype != 'pass':mask = 1 - maskbpf_dft_shift = dft_shift * mask         magnitude_spectrum = cv.magnitude(bpf_dft_shift[:,:,0], bpf_dft_shift[:,:,1])log_magnitude_spectrum = 20*np.log(magnitude_spectrum+1)bpf_dft = np.fft.ifftshift(bpf_dft_shift)         img_ = cv.idft(bpf_dft)                           img_bpf = cv.magnitude (img_[:,:,0],img_[:,:,1])return img_bpf, log_magnitude_spectrum

img_1, log_1 = notch_filter(img_gray, u0=0, v0=0, d0=150, ftype='pass')
img_2, log_2 = notch_filter(img_gray, u0=0, v0=0, d0=150, ftype='stop')

img_3, log_3 = notch_filter(img_gray, u0=150, v0=250, d0=150, ftype='pass')
img_4, log_4 = notch_filter(img_gray, u0=150, v0=250, d0=150, ftype='stop')

巴特沃斯陷波带阻滤波器

H(u,v)=11+[D02D1(u,v)D2(u,v)]NH(u, v) = \frac{1}{1 + \left[\frac{{D_0}^2}{D_1(u, v)D_2(u, v)}\right]^N}H(u,v)=1+[D1(u,v)D2(u,v)D02]N1

D1(u,v)=[(u−u0)2+(v−v0)2]12D_1(u, v) = [(u - u_0)^2 + (v - v_0)^2]^{\frac{1}{2}}D1(u,v)=[(u−u0)2+(v−v0)2]21
D2(u,v)=[(u+u0)2+(v+v0)2]12D_2(u, v) = [(u + u_0)^2 + (v + v_0)^2]^{\frac{1}{2}}D2(u,v)=[(u+u0)2+(v+v0)2]21

def bw_notch_filter(img_gray, u0=0, v0=0, d0=50, N=1, ftype='pass'):# 以频谱左上角为坐标原点dft = cv.dft(img_gray.astype('float32'),flags = cv.DFT_COMPLEX_OUTPUT) dft_shift = np.fft.fftshift(dft)                                    m, n, _ = dft_shift.shapex_arr = np.concatenate([np.arange(m).reshape(m, 1)], axis=1)y_arr = np.concatenate([np.arange(n).reshape(1, n)], axis=0)dist1 = np.sqrt((x_arr - u0)**2 + (y_arr - v0)**2)dist2 = np.sqrt((x_arr + u0)**2 + (y_arr + v0)**2)mask = 1 / (1. + ((d0**2)/(np.multiply(dist1, dist2)+0.00001))**N)if ftype == 'pass':mask = 1 - maskbpf_dft_shift = dft_shift * mask.reshape(m, n, 1)         magnitude_spectrum = cv.magnitude(bpf_dft_shift[:,:,0], bpf_dft_shift[:,:,1])log_magnitude_spectrum = 20*np.log(magnitude_spectrum+1)bpf_dft = np.fft.ifftshift(bpf_dft_shift)         img_ = cv.idft(bpf_dft)                           img_bpf = cv.magnitude (img_[:,:,0],img_[:,:,1])return img_bpf, log_magnitude_spectrum

测试效果：

高斯带阻陷波滤波器

H(u,v)=1−e[−12×D1(u,v)D2(u,v)D02]H(u, v) = 1 - e^{\left[\frac{-1}{2}\times \frac{D_1(u,v)D_2(u,v)}{D_0^2}\right]}H(u,v)=1−e[2−1×D02D1(u,v)D2(u,v)]

def gaussian_notch_filter(img_gray, u0=0, v0=0, d0=50, ftype='pass'):# 以频谱左上角为坐标原点dft = cv.dft(img_gray.astype('float32'),flags = cv.DFT_COMPLEX_OUTPUT) dft_shift = np.fft.fftshift(dft)                                    m, n, _ = dft_shift.shapex_arr = np.concatenate([np.arange(m).reshape(m, 1)], axis=1)y_arr = np.concatenate([np.arange(n).reshape(1, n)], axis=0)dist1 = np.sqrt((x_arr - u0)**2 + (y_arr - v0)**2)dist2 = np.sqrt((x_arr + u0)**2 + (y_arr + v0)**2)mask = 1 - np.exp(np.multiply(dist1, dist2)/(d0**2)*(-0.5))if ftype == 'pass':mask = 1 - maskbpf_dft_shift = dft_shift * mask.reshape(m, n, 1)         magnitude_spectrum = cv.magnitude(bpf_dft_shift[:,:,0], bpf_dft_shift[:,:,1])log_magnitude_spectrum = 20*np.log(magnitude_spectrum+1)bpf_dft = np.fft.ifftshift(bpf_dft_shift)         img_ = cv.idft(bpf_dft)                           img_bpf = cv.magnitude (img_[:,:,0],img_[:,:,1])return img_bpf, log_magnitude_spectrum

img_1, log_1 = gaussian_notch_filter(img_gray, u0=400, v0=400, d0=100, ftype='pass')
img_2, log_2 = gaussian_notch_filter(img_gray, u0=400, v0=400, d0=100, ftype='stop')