第5章 Python 数字图像处理(DIP) - 图像复原与重建12 - 空间滤波 - 使用频率域滤波降低周期噪声 - 陷波滤波、最优陷波滤波

最新推荐文章于 2023-07-20 13:23:37 发布

jasneik

最新推荐文章于 2023-07-20 13:23:37 发布

阅读量2.3k

点赞数 3

分类专栏：数字图像处理(DIP)--Python # 第5章图像复原与重建文章标签：图像处理 python numpy opencv 图像识别

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

版权

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

57 篇文章 127 订阅

订阅专栏

第5章图像复原与重建

17 篇文章 22 订阅

订阅专栏

标题

- 使用频率域滤波降低周期噪声
- - 陷波滤波深入介绍
  - 最优陷波滤波

本章陷波滤波器有部分得出的结果不佳，如果有更好的解决方案，请赐教，不胜感激。

使用频率域滤波降低周期噪声

陷波滤波深入介绍

零相移滤波器必须关于原点(频率矩形中心)对称，中以为 $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{5.33}$

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

3个 $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{5.36}$

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

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

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

def idea_notch_resistant_filter(source, uk, vk, radius=5):
    """
    create idea notch resistant filter 
    param: source: input, source image
    param: uk:     input, int, center of the height
    param: vk:     input, int, center of the width
    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)
    
    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
    k_1 = DK.copy()
    k_2 = D_K.copy()
    
    k_1[DK > D0] = 1
    k_1[DK <= D0] = 0
    
    k_2[D_K > D0] = 1
    k_2[D_K <= D0] = 0
    
    kernel = k_1 * k_2
    
    return kernel

def gauss_notch_resistant_filter(source, uk, vk, radius=5):
    """
    create gauss low pass filter 
    param: source: input, source image
    param: uk:     input, int, center of the height
    param: vk:     input, int, center of the width
    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)
    
    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
    
    k_1 = 1 - np.exp(- (DK**2)/(D0**2))   
    k_2 = 1 - np.exp(- (D_K**2)/(D0**2))   
    
    kernel = k_1 * k_2
    return kernel

def plot_3d(ax, x, y, z, cmap):
    ax.plot_surface(x, y, z, antialiased=True, shade=True, cmap=cmap)
    ax.view_init(20, -20), ax.grid(b=False), ax.set_xticks([]), ax.set_yticks([]), ax.set_zticks([])

# 理想、高斯、巴特沃斯陷波滤波器
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
from matplotlib import pyplot as plt
from matplotlib import cm

img_temp = np.zeros([256, 256])

INRF = idea_notch_resistant_filter(img_temp, radius=20, uk=30, vk=80)
GNRF = gauss_notch_resistant_filter(img_temp, radius=20, uk=30, vk=80)
BNRF = butterworth_notch_resistant_filter(img_temp, radius=20, uk=30, vk=80, n=5)

# 用来绘制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=(21, 7))
ax_1 = fig.add_subplot(1, 3, 1, projection='3d')
plot_3d(ax_1, u, v, INRF, cmap=cm.plasma)

ax_1 = fig.add_subplot(1, 3, 2, projection='3d')
plot_3d(ax_1, u, v, GNRF, cmap=cm.PiYG)

ax_1 = fig.add_subplot(1, 3, 3, projection='3d')
plot_3d(ax_1, u, v, BNRF, cmap=cm.PiYG)
plt.tight_layout()
plt.show()

在这里插入图片描述

def centralized_2d(arr):
    """
    centralized 2d array $f(x, y) (-1)^{x + y}$, about 4.5 times faster than index, 160 times faster than loop,
    """
    
    def get_1_minus1(img):
        """
        get 1 when image index is even, -1 while index is odd, same shape as input image, need this array to multiply with input image
        to get centralized image for DFT
        Parameter: img: input, here we only need img shape to create the array
        return such a [[1, -1, 1], [-1, 1, -1]] array, example is 3x3
        """
        dst = np.ones(img.shape)
        dst[1::2, ::2] = -1
        dst[::2, 1::2] = -1
        return dst

    #==================中心化=============================
    mask = get_1_minus1(arr)
    dst = arr * mask

    return dst

def pad_image(img, mode='constant'):
    """
    pad image into PxQ shape, orginal is in the top left corner.
    param: img: input image
    param: mode: input, str, is numpy pad parameter, default is 'constant'. for more detail please refere to Numpy pad
    return PxQ shape image padded with zeros or other values
    """
    dst = np.pad(img, ((0, img.shape[0]), (0, img.shape[1])), mode=mode)
    return dst

def add_sin_noise(img, scale=1, angle=0):
    """
    add sin noise for image
    param: img: input image, 1 channel, dtype=uint8
    param: scale: sin scaler, smaller than 1, will enlarge, bigger than 1 will shrink
    param: angle: angle of the rotation
    return: output_img: output image is [0, 1] image which you could use as mask or any you want to
    """
    
    height, width = img.shape[:2]  # original image shape
    
    # convert all the angle
    if int(angle / 90) % 2 == 0:
        rotate_angle = angle % 90
    else:
        rotate_angle = 90 - (angle % 90)
    
    rotate_radian = np.radians(rotate_angle)    # convert angle to radian
    
    # get new image height and width
    new_height = int(np.ceil(height * np.cos(rotate_radian) + width * np.sin(rotate_radian)))
    new_width = int(np.ceil(width * np.cos(rotate_radian) + height * np.sin(rotate_radian))) 
    
    # if new height or new width less than orginal height or width, the output image will be not the same shape as input, here set it right
    if new_height < height:
        new_height = height
    if new_width < width:
        new_width = width
    
    # meshgrid
    u = np.arange(new_width)
    v = np.arange(new_height)
    u, v = np.meshgrid(u, v)
    
    # get sin noise image, you could use scale to make some difference, better you could add some shift
#     noise = abs(np.sin(u * scale))
    noise = 1 - np.sin(u * scale)
    
    # here use opencv to get rotation, better write yourself rotation function
    C1 = cv2.getRotationMatrix2D((new_width/2.0, new_height/2.0), angle, 1)
    new_img = cv2.warpAffine(noise, C1, (int(new_width), int(new_height)), borderValue=0)
      
    # ouput image should be the same shape as input, so caculate the offset the output image and the new image
    # I make new image bigger so that it will cover all output image
    
    offset_height = abs(new_height - height) // 2
    offset_width = abs(new_width - width) // 2
    img_dst = new_img[offset_height:offset_height + height, offset_width:offset_width+width]
    output_img = normalize(img_dst)
    
    return output_img

def spectrum_fft(fft):
    """
    return FFT spectrum
    """
    return np.sqrt(np.power(fft.real, 2) + np.power(fft.imag, 2))

# 陷波滤波器处理周期噪声
img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH05/Fig0507(a)(ckt-board-orig).tif', 0) #直接读为灰度图像

# 正弦噪声
noise = add_sin_noise(img_ori, scale=0.35, angle=-20)
img = np.array(img_ori / 255, np.float32)
img_noise = img + noise
img_noise = np.uint8(normalize(img_noise)*255)

# 频率域中的其他特性
# FFT
img_fft = np.fft.fft2(img_noise.astype(np.float32))
# 中心化
fshift = np.fft.fftshift(img_fft)            # 将变换的频率图像四角移动到中心
# 中心化后的频谱
spectrum_fshift = spectrum_fft(fshift)
spectrum_fshift_n = np.uint8(normalize(spectrum_fshift) * 255)

# 对频谱做对数变换
spectrum_log = np.log(1 + spectrum_fshift)

BNRF = butterworth_notch_resistant_filter(img_ori, radius=5, uk=25, vk=10, n=4)

f1shift = fshift * (BNRF)
f2shift = np.fft.ifftshift(f1shift) #对新的进行逆变换
img_new = np.fft.ifft2(f2shift)
img_new = np.abs(img_new)

plt.figure(figsize=(15, 15))
plt.subplot(221), plt.imshow(img_noise, 'gray'), plt.title('With Sine noise'), plt.xticks([]),plt.yticks([])
plt.subplot(222), plt.imshow(spectrum_log, 'gray'), plt.title('Spectrum'), plt.xticks([]),plt.yticks([])
# 在图像上加上箭头
plt.arrow(180, 180, 25, 30, width=5,length_includes_head=True, shape='full')
plt.arrow(285, 265, -25, -30, width=5,length_includes_head=True, shape='full')

plt.subplot(223), plt.imshow(BNRF, 'gray'), plt.title('Spectrum'), plt.xticks([]),plt.yticks([])
# 在图像上加上箭头
plt.arrow(180, 180, 25, 30, width=5,length_includes_head=True, shape='full')
plt.arrow(285, 265, -25, -30, width=5,length_includes_head=True, shape='full')

plt.subplot(224), plt.imshow(img_new, 'gray'), plt.title('Spectrum'), plt.xticks([]),plt.yticks([])

plt.tight_layout()
plt.show()

在这里插入图片描述

# 陷波滤波器提取周期噪声
img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH05/Fig0507(a)(ckt-board-orig).tif', 0) #直接读为灰度图像

# 正弦噪声
noise = add_sin_noise(img_ori, scale=0.35, angle=-20)
img = np.array(img_ori / 255, np.float32)
img_noise = img + noise
img_noise = np.uint8(normalize(img_noise)*255)

# 频率域中的其他特性
# FFT
img_fft = np.fft.fft2(img_noise.astype(np.float32))
# 中心化
fshift = np.fft.fftshift(img_fft)            # 将变换的频率图像四角移动到中心
# 中心化后的频谱
spectrum_fshift = spectrum_fft(fshift)
spectrum_fshift_n = np.uint8(normalize(spectrum_fshift) * 255)

# 对频谱做对数变换
spectrum_log = np.log(1 + spectrum_fshift)

BNRF = 1 - butterworth_notch_resistant_filter(img_ori, radius=5, uk=25, vk=10, n=4)

f1shift = fshift * (BNRF)
f2shift = np.fft.ifftshift(f1shift) #对新的进行逆变换
img_new = np.fft.ifft2(f2shift)
img_new = np.abs(img_new)

plt.figure(figsize=(15, 15))
# plt.subplot(221), plt.imshow(img_noise, 'gray'), plt.title('With Sine noise'), plt.xticks([]),plt.yticks([])
# plt.subplot(222), plt.imshow(spectrum_log, 'gray'), plt.title('Spectrum'), plt.xticks([]),plt.yticks([])
# # 在图像上加上箭头
# plt.arrow(180, 180, 25, 30, width=5,length_includes_head=True, shape='full')
# plt.arrow(285, 265, -25, -30, width=5,length_includes_head=True, shape='full')

# plt.subplot(223), plt.imshow(BNRF, 'gray'), plt.title('Spectrum'), plt.xticks([]),plt.yticks([])
# # 在图像上加上箭头
# plt.arrow(180, 180, 25, 30, width=5,length_includes_head=True, shape='full')
# plt.arrow(285, 265, -25, -30, width=5,length_includes_head=True, shape='full')

plt.subplot(224), plt.imshow(img_new, 'gray'), plt.title('Sine pattern'), plt.xticks([]),plt.yticks([])

plt.tight_layout()
plt.show()

在这里插入图片描述

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

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

# 陷波滤波器处理周期噪声，用巴特沃斯低通滤波器得到的效果比目前的陷波滤波器效果还要好
img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH05/Fig0516(a)(applo17_boulder_noisy).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_fshift = spectrum_fft(fft)
spectrum_log = np.log(1 + spectrum_fshift)

# 滤波器
n = 15
r = 20
H = butterworth_low_pass_filter(fp, fp.shape, radius=100, n=4)
# BNRF_1 = butterworth_notch_resistant_filter(fp, radius=r, uk=355, vk=0, n=n)
# BNRF_2 = butterworth_notch_resistant_filter(fp, radius=r, uk=0, vk=355, n=n)
# BNRF_3 = butterworth_notch_resistant_filter(fp, radius=r, uk=250, vk=250, n=n)
# BNRF_4 = butterworth_notch_resistant_filter(fp, radius=r, uk=-250, vk=250, n=n)
# BNRF = BNRF_1 * BNRF_2 * BNRF_3 * BNRF_4  * H
BNRF = H

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)

plt.figure(figsize=(15, 12))
plt.subplot(221), plt.imshow(img_ori, 'gray'), plt.title('With noise'), plt.xticks([]),plt.yticks([])
plt.subplot(222), plt.imshow(spectrum_log, 'gray'), plt.title('Spectrum'), plt.xticks([]),plt.yticks([])
plt.subplot(223), plt.imshow(spectrum_filter_log, 'gray'), plt.title('Spectrum With Filter'), plt.xticks([]),plt.yticks([])
plt.subplot(224), plt.imshow(img_new, 'gray'), plt.title('IDFT'), plt.xticks([]),plt.yticks([])
plt.tight_layout()
plt.show()

在这里插入图片描述

# 陷波滤波器提取周期噪声
img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH05/Fig0516(a)(applo17_boulder_noisy).tif', 0)
M, N = img_ori.shape[:2]

fp = pad_image(img_ori, mode='constant')
fp_cen = centralized_2d(fp)
fft = np.fft.fft2(fp_cen)

# 中心化后的频谱
spectrum_fshift = spectrum_fft(fft)
spectrum_log = np.log(1 + spectrum_fshift)

# 滤波器
n = 15
r = 20
H = butterworth_low_pass_filter(fp, fp.shape, radius=100, n=3)
# BNRF_1 = butterworth_notch_resistant_filter(fp, radius=r, uk=355, vk=0, n=n)
# BNRF_2 = butterworth_notch_resistant_filter(fp, radius=r, uk=0, vk=355, n=n)
# BNRF_3 = butterworth_notch_resistant_filter(fp, radius=r, uk=250, vk=250, n=n)
# BNRF_4 = butterworth_notch_resistant_filter(fp, radius=r, uk=-250, vk=250, n=n)
# BNRF = BNRF_1 * BNRF_2 * BNRF_3 * BNRF_4 * H
BNRF = H
fft_filter = fft * (1 - 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)

plt.figure(figsize=(15, 12))
# plt.subplot(221), plt.imshow(img_ori, 'gray'), plt.title('With Sine noise'), plt.xticks([]),plt.yticks([])
# plt.subplot(222), plt.imshow(spectrum_log, 'gray'), plt.title('Spectrum'), plt.xticks([]),plt.yticks([])
# plt.subplot(223), plt.imshow(spectrum_filter_log, 'gray'), plt.title('Spectrum'), plt.xticks([]),plt.yticks([])
plt.subplot(224), plt.imshow(img_new, 'gray'), plt.title('Spectrum'), plt.xticks([]),plt.yticks([])
plt.tight_layout()
plt.show()

在这里插入图片描述

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

# 陷波滤波器处理周期噪声，用巴特沃斯低通滤波器得到的效果比目前的陷波滤波器效果还要好
img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH05/Fig0519(a)(florida_satellite_original).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_fshift = spectrum_fft(fft)
spectrum_log = np.log(1 + spectrum_fshift)

# 滤波器
n = 15
r = 20
H = narrow_notch_filter(fp, w=10, opening=30, vertical=True, horizontal=False)
# BNRF_1 = butterworth_notch_resistant_filter(fp, radius=r, uk=355, vk=0, n=n)
# BNRF_2 = butterworth_notch_resistant_filter(fp, radius=r, uk=0, vk=355, n=n)
# BNRF_3 = butterworth_notch_resistant_filter(fp, radius=r, uk=250, vk=250, n=n)
# BNRF_4 = butterworth_notch_resistant_filter(fp, radius=r, uk=-250, vk=250, n=n)
# BNRF = BNRF_1 * BNRF_2 * BNRF_3 * BNRF_4  * H
BNRF = H

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)

plt.figure(figsize=(15, 16))
plt.subplot(221), plt.imshow(img_ori, 'gray'), plt.title('With noise'), plt.xticks([]),plt.yticks([])
plt.subplot(222), plt.imshow(spectrum_log, 'gray'), plt.title('Spectrum'), plt.xticks([]),plt.yticks([])
plt.subplot(223), plt.imshow(spectrum_filter_log, 'gray'), plt.title('Spectrum With Filter'), plt.xticks([]),plt.yticks([])
plt.subplot(224), plt.imshow(img_new, 'gray'), plt.title('IDFT'), plt.xticks([]),plt.yticks([])
plt.tight_layout()
plt.show()

在这里插入图片描述

# 陷波滤波器提取周期噪声，用巴特沃斯低通滤波器得到的效果比目前的陷波滤波器效果还要好
img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH05/Fig0519(a)(florida_satellite_original).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_fshift = spectrum_fft(fft)
spectrum_log = np.log(1 + spectrum_fshift)

# 滤波器
n = 15
r = 20
H = narrow_notch_filter(fp, w=10, opening=30, vertical=True, horizontal=False)
# BNRF_1 = butterworth_notch_resistant_filter(fp, radius=r, uk=355, vk=0, n=n)
# BNRF_2 = butterworth_notch_resistant_filter(fp, radius=r, uk=0, vk=355, n=n)
# BNRF_3 = butterworth_notch_resistant_filter(fp, radius=r, uk=250, vk=250, n=n)
# BNRF_4 = butterworth_notch_resistant_filter(fp, radius=r, uk=-250, vk=250, n=n)
# BNRF = BNRF_1 * BNRF_2 * BNRF_3 * BNRF_4  * H
BNRF = H

fft_filter = fft * (1 - 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)

plt.figure(figsize=(15, 16))
# plt.subplot(221), plt.imshow(img_ori, 'gray'), plt.title('With noise'), plt.xticks([]),plt.yticks([])
# plt.subplot(222), plt.imshow(spectrum_log, 'gray'), plt.title('Spectrum'), plt.xticks([]),plt.yticks([])
# plt.subplot(223), plt.imshow(spectrum_filter_log, 'gray'), plt.title('Spectrum With Filter'), plt.xticks([]),plt.yticks([])
plt.subplot(224), plt.imshow(img_new, 'gray'), plt.title('IDFT'), plt.xticks([]),plt.yticks([])
plt.tight_layout()
plt.show()

在这里插入图片描述

# 使用陷波带阻滤波器滤波
img_florida = cv2.imread("DIP_Figures/DIP3E_Original_Images_CH05/Fig0519(a)(florida_satellite_original).tif", -1)

#--------------------------------
fft = np.fft.fft2(img_florida)
fft_shift = np.fft.fftshift(fft)
amp_img = np.abs(np.log(1 + np.abs(fft_shift)))

#--------------------------------
BNF = narrow_notch_filter(img_florida, w=5, opening=20, vertical=True, horizontal=False)

fft_NNF = np.fft.fft2(BNF*255)
fft_shift_NNF = np.fft.fftshift(fft_NNF)
amp_img_NNF = np.abs(np.log(1 + np.abs(fft_shift_NNF)))


#--------------------------------
f1shift = fft_shift * (BNF)
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))

fft_mask = amp_img * BNF

plt.figure(figsize=(15, 16))
plt.subplot(221),plt.imshow(img_florida,'gray'),plt.title('Image with noise')
plt.subplot(222),plt.imshow(amp_img,'gray'),plt.title('FFT')
plt.subplot(223),plt.imshow(fft_mask,'gray'),plt.title('FFT with mask')
plt.subplot(224),plt.imshow(img_new,'gray'),plt.title('Denoising')
plt.tight_layout()
plt.show()

在这里插入图片描述

最优陷波滤波

这种滤波方法的过程如下：
首先分离干扰模式的各个主要贡献，然后从被污染图像中减去该模式的一个可变加权部分。

首先提取干模式的主频率分量，提取方法是在每个尖峰位置放一个陷波带通滤波器传递函数 $H_{NP}(u, v)$ ，则干扰噪声模式的傅里叶变换为：
$H_{NP}(u, v)G(u, v) \tag{5.38}$

则有噪声模式：
$\eta(x, y) = \mathfrak{J}^-1 \{ H_{NP}(u, v)G(u, v) \} \tag{5.39}$

如果我们知道了噪声模式，我们假设噪声是加性噪声，只可以用污染的噪声 $g (x, y)$ 减去噪声模式 $\eta(x, y)$ 可得到 $\hat{f}(x, y)$ ，但通常这只是一个近似值。
$\hat{f}(x, y) = g(x, y) - w(x, y)\eta(x, y) \tag{5.40}$
$w (x, y)$ 是一个加权函数或调制函数，这个方法的目的就是选取 $w (x, y)$ ，以便以某种意义的方式来优化结果。一种方法是选择 $w (x, y)$ ，使 $\hat{f}(x, y)$ 在每点 $(x, y)$ 的规定邻域上的方差最小。

$m\times{n}$ (奇数)的邻域 $S_{xy}$ 。 $\hat{f}(x, y)$ 的“局部”方差估计如下：
$\sigma^2(x, y) = \frac{1}{mn} \sum_{(r,c)\in S_{xy}} \Big[ \hat{f}(r, c) - \bar{\hat{f}} \Big]^2 \tag{5.41}$

$\bar{\hat{f}}$ 是邻域 $\hat{f}$ 的平均值，
$\bar{\hat{f}} = \frac{1}{mn} \sum_{(r,c)\in S_{xy}} \hat{f}(r, c) \tag{5.42}$

将式(5.40)代入(5.41)，得
$\sigma^2(x, y) = \frac{1}{mn} \sum_{(r,c)\in S_{xy}} \Big\{\big[g(r, c) - w(r, c) \eta(r, c)\big] - \big[\overline{g} - \overline{w\eta} \big ] \Big\}^2\tag{5.43}$

$\overline{g}$ 和 $\overline{w\eta}$ 分别是 $g$ 和 $w\eta$ 在邻域 $S_{xy}$ 的平均值

若假设 $w$ 在 $S_{xy}$ 内近似为常数，则可用该邻域中心的 $w$ 值来代替 $w (r, c)$ :

$\tag{5.44}$

因为 $w (x, y)$ 在 $S_{xy}$ 中被假设为常数，因此在 $S_{xy}$ 中根据 $\overline{w} = w(x, y)$ 有

$\overline{w\eta} = w(x, y) \overline{\eta} \tag{5.45}$

$\overline{\eta}$ 是邻域 $S_{xy}$ 中的平均值，所以式(5.43)变为：

$\sigma^2(x, y) = \frac{1}{mn} \sum_{(r,c)\in S_{xy}} \Big\{\big[g(r, c) - w(x, y) \eta(r, c)\big] - \big[\overline{g} - {w(x, y)}\overline{\eta} \big ] \Big\}^2\tag{5.44}$

要使得 $\sigma^2(x, y)$ 相对 $w (x, y)$ 最小，我们可以对式(5.44)求关于 $w (x, y)$ 的偏导数，并令为偏导数为0；

$\frac{\partial{\sigma^2(x, y)}}{\partial{w(x, y)}} = 0 \tag{5.47}$

求得 $w (x, y)$ :
$\frac{\overline{g\eta} - \bar{g}\bar{\eta}}{\overline{\eta^2} - \bar{\eta}^2}\tag{5.48}$

把式(5.48)代入式(5.40)并在噪声图像 $g$ 中的每个点执行这一过程，可得到完全复原的图像。

# 这里还没有实现，迟点再弄吧
img_mariner = cv2.imread("DIP_Figures/DIP3E_Original_Images_CH05/Fig0520(a)(NASA_Mariner6_Mars).tif", 0)
M, N = img_mariner.shape[:2]

fp = pad_image(img_mariner, mode='reflect')
fp_cen = centralized_2d(fp)
fft = np.fft.fft2(fp_cen)

# 中心化后的频谱
spectrum_fshift = spectrum_fft(fft)
spectrum_log = np.log(1 + spectrum_fshift)

# 未中心化的频谱
fft_fp = np.fft.fft2(fp)
spectrum_fp = spectrum_fft(fft_fp)
spectrum_fp_log = np.log(1 + spectrum_fp)

# 滤波器
n = 15
r = 20
H = butterworth_band_resistant_filter(fp, fp.shape, radius=40, w=5, n=5)
# BNRF_1 = butterworth_notch_resistant_filter(fp, radius=r, uk=355, vk=0, n=n)
# BNRF_2 = butterworth_notch_resistant_filter(fp, radius=r, uk=0, vk=355, n=n)
# BNRF_3 = butterworth_notch_resistant_filter(fp, radius=r, uk=250, vk=250, n=n)
# BNRF_4 = butterworth_notch_resistant_filter(fp, radius=r, uk=-250, vk=250, n=n)
# BNRF = BNRF_1 * BNRF_2 * BNRF_3 * BNRF_4  * H

fft_filter = fft_fp * (1 - H)
ifft = np.fft.ifft2(fft_filter)
img_new = ifft.real[:M, :N]

# # show = spectrum_fp_log * H
# 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)

plt.figure(figsize=(15, 15))
plt.subplot(221), plt.imshow(img_mariner, 'gray'), plt.title('With noise'), plt.xticks([]),plt.yticks([])
plt.subplot(222), plt.imshow(spectrum_log, 'gray'), plt.title('Spectrum Centralied'), plt.xticks([]),plt.yticks([])
plt.subplot(223), plt.imshow(spectrum_fp_log, 'gray'), plt.title('Spectrum Not Centralized'), plt.xticks([]),plt.yticks([])
plt.subplot(224), plt.imshow(img_new, 'gray'), plt.title('IDFT'), plt.xticks([]),plt.yticks([])
plt.tight_layout()
plt.show()

在这里插入图片描述

# 巴特沃斯带阻陷波滤波器 BNRF
img_dst = img_mariner - img_new
plt.figure(figsize=(16, 16))
plt.subplot(221), plt.imshow(img_dst, 'gray'), plt.title('BNF_1')
# plt.subplot(222), plt.imshow(BNF_2, 'gray'), plt.title('BNF_2')
# plt.subplot(223), plt.imshow(BNF_3, 'gray'), plt.title('BNF_3')
# plt.subplot(224), plt.imshow(BNF_dst, 'gray'), plt.title('BNF_dst')
plt.tight_layout()
plt.show()

在这里插入图片描述

jasneik

关注

3
点赞
踩
36

收藏

觉得还不错? 一键收藏
打赏
2
评论
第5章 Python 数字图像处理(DIP) - 图像复原与重建12 - 空间滤波 - 使用频率域滤波降低周期噪声 - 陷波滤波、最优陷波滤波

标题使用频率域滤波降低周期噪声陷波滤波深入介绍最优陷波滤波本章陷波滤波器有部分得出的结果不佳，如果有更好的解决方案，请赐教，不胜感激。使用频率域滤波降低周期噪声陷波滤波深入介绍零相移滤波器必须关于原点(频率矩形中心)对称，中以为(u0,v0)(u_0, v_0)(u0,v0)的陷波滤波器传递函数在(−u0,−v0)(-u_0, -v_0)(−u0,−v0)位置必须有一个对应的陷波。陷波带阻滤波器传递函数可用中心被平移到陷波滤波中心的高通滤波器函数的乘积来产生HNR(u,v)=∏k=1QHk
复制链接

扫一扫