第5章 Python 数字图像处理(DIP) - 图像复原与重建13 - 空间滤波 - 线性位置不变退化 - 退化函数估计、运动模糊函数

在这一节中,得到的结果,有些不是很好,我需要再努力多找资料,重新完成学习,如果大佬有相关资料推荐,不胜感激。

线性位置不变退化

# 巴特沃斯带阻陷波滤波器 BNRF
img_temp = np.zeros([512, 512])
BNF_1 = butterworth_notch_resistant_filter(img_temp, radius=20, uk=-80, vk=60)
BNF_2 = butterworth_notch_resistant_filter(img_temp, radius=10, uk=30, vk=80)
BNF_3 = butterworth_notch_resistant_filter(img_temp, radius=10, uk=-30, vk=80)

plt.figure(figsize=(16, 16))
plt.subplot(221), plt.imshow(BNF_1, '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')

BNF_dst = BNF_1 * BNF_2 * BNF_3

plt.subplot(224), plt.imshow(BNF_dst, 'gray'), plt.title('BNF_dst')

plt.tight_layout()
plt.show()

估计退化函数

In this section, I think I still got some problem have to sort out, when I have some more time or some more reading.

采用观察法估计退化函数

选择一个信号内容很强的区域(如一个高对比度区域)表示为 g ( x , y ) g(x, y) g(x,y),令 f ^ ( x , y ) \hat{f}(x, y) f^(x,y)表示为处理后的子图像,则有:
H s ( u , v ) = G s ( u , v ) F ^ s ( u , v ) (5.66) H_{s}(u, v) = \frac{G_{s}(u, v)}{\hat{F}_{s}(u, v)} \tag{5.66} Hs(u,v)=F^s(u,v)Gs(u,v)(5.66)

根据位置不变的假设来推断完整的退化函数 H ( u , v ) H(u, v) H(u,v)

采用试验法估计退化函数

一个冲激由一个亮点来模拟,这个点应亮到能降低噪声对可忽略值的影响。一个冲激的傅里叶变换是一个常量:
H ( u , v ) = G ( u , v ) A (5.67) H(u, v) = \frac{G(u, v)}{A} \tag{5.67} H(u,v)=AG(u,v)(5.67)

# 试验法估计退化函数
img_impulse = cv2.imread("DIP_Figures/DIP3E_Original_Images_CH05/Fig0524(a)(impulse).tif", 0)
img_blurred = cv2.imread("DIP_Figures/DIP3E_Original_Images_CH05/Fig0524(b)(blurred-impulse).tif", 0)

fig = plt.figure(figsize=(10, 5))
plt.subplot(1, 2, 1), plt.imshow(img_impulse, cmap='gray'), plt.xticks([]), plt.yticks([])
plt.subplot(1, 2, 2), plt.imshow(img_blurred, cmap='gray'), plt.xticks([]), plt.yticks([])
plt.tight_layout()
plt.show()

在这里插入图片描述
下面两个例子,你会看到傅立叶变换后,频谱图像的美。把频谱图像上了颜色后,更是美极啦

# 傅里叶变换
fp_impulse = pad_image(img_impulse)
impluse_cen = centralized_2d(fp_impulse)
fft_impulse = np.fft.fft2(impluse_cen)
impulse_spectrume = np.log(1 + spectrum_fft(fft_impulse))

fp_blurred = pad_image(img_blurred)
blurred_cen = centralized_2d(fp_blurred)
fft_blurred = np.fft.fft2(blurred_cen)
blurred_spectrum = np.log(1 + spectrum_fft(fft_blurred))

H = fft_blurred / fft_impulse

h_spectrum = np.log(1 + spectrum_fft(H))
h_spectrum = h_spectrum / h_spectrum.max()

fig = plt.figure(figsize=(15, 5))
plt.subplot(1, 3, 1), plt.imshow(impulse_spectrume, cmap='gray'), plt.xticks([]), plt.yticks([])
plt.subplot(1, 3, 2), plt.imshow(blurred_spectrum, cmap='gray'), plt.xticks([]), plt.yticks([])
plt.subplot(1, 3, 3), plt.imshow(h_spectrum, cmap='gray'), plt.xticks([]), plt.yticks([])
plt.tight_layout()
plt.show()

在这里插入图片描述

# 一些傅里叶变换
img_temp = np.zeros([256, 256])
# H = butterworth_low_pass_filter(img_temp, 10, 500)
H = 1 - butterworth_band_resistant_filter(img_temp, img_temp.shape, radius=50, w=5, n=5)
fp_blurred = pad_image(H)
blurred_cen = centralized_2d(fp_blurred)
fft_blurred = np.fft.fft2(blurred_cen)
blurred_spectrum = np.log(1 + spectrum_fft(fft_blurred))

fig = plt.figure(figsize=(15, 15))
plt.imshow(blurred_spectrum, cmap='PiYG'), plt.xticks([]), plt.yticks([])
# plt.savefig("bbrf_4.png", dpi=300, quality=100)
# plt.subplot(1, 3, 1), plt.imshow(impulse_spectrume, cmap='gray'), plt.xticks([]), plt.yticks([])
# plt.subplot(1, 3, 2), plt.imshow(blurred_spectrum, cmap='gray'), plt.xticks([]), plt.yticks([])
# # plt.subplot(1, 3, 3), plt.imshow(h_spectrum, cmap='gray'), plt.xticks([]), plt.yticks([])
plt.tight_layout()
plt.show()

在这里插入图片描述

采用建模法估计退化函数

H ( u , v ) = e − k ( u 2 + v 2 ) 5 6 (5.68) H(u,v) = e^{-k(u^2 + v^2)^{\frac{5}{6}}} \tag{5.68} H(u,v)=ek(u2+v2)65(5.68)

关于频率矩形的中心,可用如下函数
H ( u , v ) = e − k ( ( u − P / 2 ) 2 + ( v − Q / 2 ) 2    ) 5 6 H(u, v) = e^{-k((u - P/2)^2 + (v - Q/2)^2 \ \ )^{\frac{5}{6}}} H(u,v)=ek((uP/2)2+(vQ/2)2  )65

参加书上P247页,运动导的图像模糊的退化过程,是否用错?
这个问题已经得到解决啦,解决方案如下。

def modeling_degrade(img, k=1):
    """
    modeling degradation fuction, math: $$H(u,v) = e^{-k(u^2 +v^2)^{\frac{5}{6}}}$$
    param: img: input img
    param: k: 
    """
    N, M = img.shape[:2]
    
    u = np.arange(M)
    v = np.arange(N)
    
    u, v = np.meshgrid(u, v)
    
    temp = (u - M//2)**2 + (v - N//2)**2
    
    kernel = np.exp(-k * np.power(temp, 5/6))
    
    return kernel
# 不填充,结果与书上一致啦
def get_degenerate_image(img, img_deg):
    """
    不填充图像做傅里叶变换后与退化函数做乘积,再反傅里叶变换
    """

    # FFT--------------------------------------------
    fft = np.fft.fft2(img)

    # FFT * H(u, v)----------------------------------
    fft_huv = fft * img_deg

    # IFFT-------------------------------------------
    ifft = np.fft.ifft2(fft_huv)
    
    return ifft

img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH05/Fig0525(a)(aerial_view_no_turb).tif', 0)

# k = [1, 0.1, 0.01, 0.001, 0.0025, 0.00025]
k = [0.0025, 0.001, 0.00025]

fp_cen = centralized_2d(img_ori)

fig = plt.figure(figsize=(12, 12))
for i in range(len(k) + 1):
    ax = fig.add_subplot(2, 2, i+1, xticks=[], yticks=[])
    if i == 0:
        ax.imshow(img_ori, 'gray'), ax.set_title(f"Original")
    else:
        img_deg = modeling_degrade(fp_cen, k=k[i-1])
        ifft = get_degenerate_image(fp_cen, img_deg)
        img_new = centralized_2d(ifft.real)
        img_new = np.clip(img_new, 0, img_new.max())
        img_new = np.uint8(normalize(img_new) * 255)
        ax.imshow(img_new, 'gray')
        ax.set_title(f"k = {k[i-1]}")
plt.tight_layout()
plt.show()    

在这里插入图片描述

运动模糊函数

H ( u , v ) = T π ( u a + v b ) s i n [ π ( u a + v b ) ] e − j π ( u a + v b ) H(u,v) =\frac{T}{\pi(ua + vb)}sin[\pi(ua+vb)]e^{-j\pi(ua+vb)} H(u,v)=π(ua+vb)Tsin[π(ua+vb)]ejπ(ua+vb)

下面的代码可能比较混乱,因为实验过程,而得出的结果不太好,还没有整理。需要继续学习后,再完成整理。

img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH05/Fig0526(a)(original_DIP).tif', 0)
def motion_huv(img, a, b, T):
    
    eps = 1e-8
    M, N = img.shape[1], img.shape[0]
    
    u = np.arange(1, M+1)
    v = np.arange(1, N+1)
    
    u, v = np.meshgrid(u, v)
    
    temp = np.pi * (u * a + v * b)
    
    kernel = (T * np.sin(temp) * np.exp(-temp*1j) /(temp + eps))
    
    return kernel
# 对图片进行运动模糊
def make_blurred(img, PSF, eps):

    #=====================
#     fft = np.fft.fft2(img)
# #     fft_shift = np.fft.fftshift(fft)
    
#     fft_psf = fft * PSF
    
#     ifft = np.fft.ifft2(fft_psf)
# #     ifft_shift = np.fft.ifftshift(ifft)
#     blurred = abs(ifft.real)
    
    #=========================
    M, N = img.shape[:2]
    fp = pad_image(img, mode='constant')
    fp_cen = centralized_2d(fp)
    
    img_fft = np.fft.fft2(fp_cen)
    img_fft_psf = img_fft * PSF

    ifft = np.fft.ifft2(img_fft_psf)
    blurred = centralized_2d(ifft.real)[:N, :M]
# #     blurred = ifft.real[:N, :M]
    
    return blurred
def get_motion_dsf(image_size, motion_angle, motion_dis):
    PSF = np.zeros(image_size)  # 点扩散函数
    x_center = (image_size[0] - 1) / 2
    y_center = (image_size[1] - 1) / 2
 
    sin_val = np.sin(motion_angle * np.pi / 180)
    cos_val = np.cos(motion_angle * np.pi / 180)
 
    # 将对应角度上motion_dis个点置成1
    for i in range(motion_dis):
        x_offset = round(sin_val * i)
        y_offset = round(cos_val * i)
        PSF[int(x_center - x_offset), int(y_center + y_offset)] = 1
 
    return PSF / PSF.sum()    # 归一化
img_motion = get_motion_dsf((480, 480), 70, 200)

plt.figure(figsize=(10, 8))
plt.subplot(121), plt.imshow(img_motion,'gray'), plt.title('img_motion')
plt.show()

在这里插入图片描述

OpenCV Motion Blur
def motion_blur(image, degree=12, angle=45):
    """
    create motion blur using opencv
    param: image: input image
    param: degree: the size of the blurry
    param: angle: blur angle
    return uint8 image
    """
    image = np.array(image)

    # 这里生成任意角度的运动模糊kernel的矩阵, degree越大,模糊程度越高
    M = cv2.getRotationMatrix2D((degree / 2, degree / 2), angle, 1)
    motion_blur_kernel = np.diag(np.ones(degree))
    motion_blur_kernel = cv2.warpAffine(motion_blur_kernel, M, (degree, degree))

    motion_blur_kernel = motion_blur_kernel / degree
    blurred = cv2.filter2D(image, -1, motion_blur_kernel)

    # convert to uint8
    cv2.normalize(blurred, blurred, 0, 255, cv2.NORM_MINMAX)
    blurred = np.array(blurred, dtype=np.uint8)
    return blurred
# 运动模糊图像
img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH05/Fig0526(a)(original_DIP).tif', 0)

img_blur = motion_blur(img_ori, degree=75, angle=15)

plt.figure(figsize=(12, 8))
plt.subplot(121), plt.imshow(img_ori,'gray'), plt.title('img_deg')
plt.subplot(122), plt.imshow(img_blur,'gray'), plt.title('high_pass')
plt.show()

在这里插入图片描述

  • 0
    点赞
  • 20
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

jasneik

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值