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

最新推荐文章于 2024-05-10 10:23:26 发布

jasneik

最新推荐文章于 2024-05-10 10:23:26 发布

阅读量1.4k

点赞数 4

分类专栏：数字图像处理(DIP)--Python # 第4章 - 频率域滤波文章标签：图像处理 python numpy opencv 算法

本文链接：https://blog.csdn.net/jasneik/article/details/114045829

版权

数字图像处理(DIP)--Python 同时被 2 个专栏收录

57 篇文章 127 订阅

订阅专栏

第4章 - 频率域滤波

12 篇文章 6 订阅

订阅专栏

选择性滤波

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

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

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

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

带阻滤波器和带通滤波器

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

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

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

高斯带阻滤波器传递函数

$e^{-\big[\frac{(D(u, v) - C_0)^2}{W^2} \ \ \big]} \tag{4.149}$

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

$e^{-\big[\frac{D^2(u, v) - C_0^2}{D(u, v) W} \ \ \big]^2} \tag{4.150}$

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

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

def idea_band_resistant_filter(source, center, radius=10, w=5):
    """
    create idea band resistant 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, int, the radius of circle of the band pass filter, default is 10
    param: w:      input, int, the width of the band of the filter, default is 5
    return 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 = radius
    half_w = w / 2
    kernel_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.149
    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, int, the radius of circle of the band pass filter, default is 10
    param: w:      input, int, the width of the band of the filter, default is 5
    return 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 = radius
    
    kernel = 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.150
    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, int, the radius of circle of the band pass filter, default is 10
    param: w:      input, int, the width of the band of the filter, default is 5
    return 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 = radius
    
    kernel = 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.150
    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, int, the radius of circle of the band pass filter, default is 10
    param: w:      input, int, the width of the band of the filter, default is 5
    param: 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 = radius
    
    temp = (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 cm

img_temp = np.zeros([512, 512])
center = img_temp.shape

radius = 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) \tag{4.151}$

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

$n$ 阶巴特沃斯带阴滤波器
$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}$

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

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

def butterworth_notch_resistant_filter(img, uk, vk, radius=10, n=1):
    """
    create butterworth notch resistant filter, equation 4.155
    param: img:    input, source image
    param: uk:     input, int, center of the height
    param: vk:     input, int, center of the width
    param: radius: input, int, the radius of circle of the band pass filter, default is 10
    param: w:      input, int, the width of the band of the filter, default is 5
    param: 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 = radius
    kernel = (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_3

plt.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 opencv
    param: img:        input, source image
    param: w:          input, int, width of the resistant, value is 0, default is 5
    param: opening:    input, int, opening of the resistant, value is 1, default is 10
    param: 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) * 255
    c_height, c_width = img.shape[0] // 2, img.shape[1] // 2
    if 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) * 255
    if 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 filter
    param: img:        input, source image
    param: w:          input, int, width of the resistant, value is 0, default is 5
    param: opening:    input, int, opening of the resistant, value is 1, default is 10
    param: 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//2
    opening_half = opening//2
    img_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, :] = 0
        img_horizontal[:, N//2 - opening_half:N//2 + opening - opening_half] = 1
    if vertical:
        img_vertical[:, N//2 - w_half:N//2 + w - w_half] = 0
        img_vertical[M//2 - opening_half:M//2 + opening - opening_half, :] = 1
    img_dst = img_horizontal * img_vertical
    
    return 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()

在这里插入图片描述

jasneik

关注

4
点赞
踩
20

收藏

觉得还不错? 一键收藏
打赏
2
评论
第4章 Python 数字图像处理(DIP) - 频率域滤波12 - 选择性滤波 - 带阻

目录频率域滤波基础频率域的其他特性频率域滤波基础知识频率域滤波步骤小结空间域和频率域滤波之间的对应关系使用低通频率域滤波器平滑图像理想低通滤波器(ILPF)高斯低通滤波器(GLPF)巴特沃斯低通滤波器低通滤波的例子使用高通滤波器锐化图像由低通滤波器得到理想、高斯和巴特沃斯高通滤波器指纹增强频域中的拉普拉斯钝化掩蔽、高提升滤波和高频强调滤波同态滤波选择性滤波带阻滤波器和带通滤波器陷波滤波器频率域滤波基础频率域的其他特性频率域中的滤波过程如下：首先修改傅里叶变换以在到特定目的然后计算IDFT，返回
复制链接

扫一扫