基于改进小波分析的图像降噪方法（Python）

哥廷根数学学派

于 2024-07-26 07:49:49 发布

阅读量40

点赞数

分类专栏：降噪图像处理小波分析文章标签： python 开发语言人工智能深度学习信号处理

本文链接：https://blog.csdn.net/weixin_39402231/article/details/140704840

版权

小波分析同时被 3 个专栏收录

154 篇文章 34 订阅

订阅专栏

图像处理

62 篇文章 1 订阅

订阅专栏

降噪

45 篇文章 0 订阅

订阅专栏

import matplotlib.pyplot as plt


from skimage.restoration import (denoise_wavelet, estimate_sigma)
from skimage import data, img_as_float
from skimage.util import random_noise
from skimage.metrics import peak_signal_noise_ratio




original = img_as_float(data.chelsea()[100:250, 50:300])


sigma = 0.12
noisy = random_noise(original, var=sigma**2)


fig, ax = plt.subplots(nrows=2, ncols=3, figsize=(8, 5),
                       sharex=True, sharey=True)


plt.gray()


# Estimate the average noise standard deviation across color channels.
sigma_est = estimate_sigma(noisy, channel_axis=-1, average_sigmas=True)
# Due to clipping in random_noise, the estimate will be a bit smaller than the
# specified sigma.
print(f'Estimated Gaussian noise standard deviation = {sigma_est}')


im_bayes = denoise_wavelet(noisy, channel_axis=-1, convert2ycbcr=True,
                           method='BayesShrink', mode='soft',
                           rescale_sigma=True)
im_visushrink = denoise_wavelet(noisy, channel_axis=-1, convert2ycbcr=True,
                                method='VisuShrink', mode='soft',
                                sigma=sigma_est, rescale_sigma=True)


# VisuShrink is designed to eliminate noise with high probability, but this
# results in a visually over-smooth appearance.  Repeat, specifying a reduction
# in the threshold by factors of 2 and 4.
im_visushrink2 = denoise_wavelet(noisy, channel_axis=-1, convert2ycbcr=True,
                                 method='VisuShrink', mode='soft',
                                 sigma=sigma_est/2, rescale_sigma=True)
im_visushrink4 = denoise_wavelet(noisy, channel_axis=-1, convert2ycbcr=True,
                                 method='VisuShrink', mode='soft',
                                 sigma=sigma_est/4, rescale_sigma=True)


# Compute PSNR as an indication of image quality
psnr_noisy = peak_signal_noise_ratio(original, noisy)
psnr_bayes = peak_signal_noise_ratio(original, im_bayes)
psnr_visushrink = peak_signal_noise_ratio(original, im_visushrink)
psnr_visushrink2 = peak_signal_noise_ratio(original, im_visushrink2)
psnr_visushrink4 = peak_signal_noise_ratio(original, im_visushrink4)


ax[0, 0].imshow(noisy)
ax[0, 0].axis('off')
ax[0, 0].set_title(f'Noisy\nPSNR={psnr_noisy:0.4g}')
ax[0, 1].imshow(im_bayes)
ax[0, 1].axis('off')
ax[0, 1].set_title(
    f'Wavelet denoising\n(BayesShrink)\nPSNR={psnr_bayes:0.4g}')
ax[0, 2].imshow(im_visushrink)
ax[0, 2].axis('off')
ax[0, 2].set_title(
    'Wavelet denoising\n(VisuShrink, $\\sigma=\\sigma_{est}$)\n'
     'PSNR=%0.4g' % psnr_visushrink)
ax[1, 0].imshow(original)
ax[1, 0].axis('off')
ax[1, 0].set_title('Original')
ax[1, 1].imshow(im_visushrink2)
ax[1, 1].axis('off')
ax[1, 1].set_title(
    'Wavelet denoising\n(VisuShrink, $\\sigma=\\sigma_{est}/2$)\n'
     'PSNR=%0.4g' % psnr_visushrink2)
ax[1, 2].imshow(im_visushrink4)
ax[1, 2].axis('off')
ax[1, 2].set_title(
    'Wavelet denoising\n(VisuShrink, $\\sigma=\\sigma_{est}/4$)\n'
     'PSNR=%0.4g' % psnr_visushrink4)
fig.tight_layout()


plt.show()

# load package
import matplotlib.pyplot as plt
import skimage
import matplotlib.image as mpig
from skimage.restoration import (denoise_wavelet, estimate_sigma)
from skimage import data, img_as_float
from skimage.util import random_noise
from skimage.metrics import peak_signal_noise_ratio,structural_similarity
from skimage import io
from numpy import *
from numpy import linalg as la
from math import log10, sqrt
import numpy as np
import matplotlib.pyplot as plt
import scipy.io as sci
from sklearn.decomposition import PCA
import numpy as np
def sep_RGB(noisy_or):
    # R
    red=noisy_or[:,:,0]
    # G
    green=noisy_or[:,:,1]
    # B
    blue=noisy_or[:,:,2]


    fig,axes = plt.subplots(1,3,figsize=(12,7))
    axes[0].imshow(noisy_or[:,:,0])
    axes[0].set_title('R')
    axes[1].imshow(noisy_or[:,:,1])
    axes[1].set_title('G')
    axes[2].imshow(noisy_or[:,:,2])
    axes[2].set_title('B')


    plt.show()


    return red,green,blue
def pic(original,x,method):
    '''
    showing the pic, PSNR and SSIM
    ---
    input: pic matrix, denoising method
    output: plot
    '''
    psnr_noisy = peak_signal_noise_ratio(original, x)
    ssm_noisy=structural_similarity(original, x,channel_axis=-1)
    plt.imshow(x)
    plt.title(f'{method}\nPSNR={psnr_noisy:0.4g}'
         '\nSSIM=%0.4g' % ssm_noisy)
    plt.savefig(f'{method}.png')
def PCA_denoising(red,green,blue,n):
    '''
    Doing pca on RGB paths repectively
    ---
    input: noisy pic
    output: denoise pic
    '''
    import matplotlib.pyplot as plt
    import scipy.io as sci
    from sklearn.decomposition import PCA
    import numpy as np
    #initialize PCA with first  principal components
    pca = PCA(n,random_state=1)
    #Applying to red channel and then applying inverse transform to transformed array.
    red_transformed = pca.fit_transform(red)
    red_inverted = pca.inverse_transform(red_transformed)


    #Applying to Green channel and then applying inverse transform to transformed array.
    green_transformed = pca.fit_transform(green)
    green_inverted = pca.inverse_transform(green_transformed)


    #Applying to Blue channel and then applying inverse transform to transformed array.
    blue_transformed = pca.fit_transform(blue)
    blue_inverted = pca.inverse_transform(blue_transformed)


    x= np.dstack((red_inverted, green_inverted, blue_inverted))


    return x
def svd(img):
    '''
    Doing SVD on each 2D data
    ---
    input: one of RGB path
    output: one of RGB SVD path
    '''
    u, sigma, vt = la.svd(img)
    h, w = img.shape[:2]
    h1 = int(h * 0.1) #The first 10% singular value is used to reconstruct the image
    sigma1 = diag(sigma[:h1],0) #Generate a diagonal matrix with singular values
    u1 = zeros((h,h1), float)
    u1[:,:] = u[:,:h1]
    vt1 = zeros((h1,w), float)
    vt1[:,:] = vt[:h1,:]
    return u1.dot(sigma1).dot(vt1)


def SVD_denoising(red,green,blue):
    '''
    Doing SVD on pic
    ---
    input:RGB matrix
    output:pic after svd
    '''
    r=svd(red)
    g=svd(green)
    b=svd(blue)
    noisy_svd= np.dstack((r, g, b))
    return noisy_svd
def sigmaest(noisy):
    '''
    Estimate the average noise standard deviation across color channels.
    Due to clipping in random_noise, the estimate will be a bit smaller than the specified sigma.
    ---
    input: noisy pic
    output: estimated sigma
    '''
    sigma_est = estimate_sigma(noisy, channel_axis=-1, average_sigmas=True)
    print(f'Estimated Gaussian noise standard deviation = {sigma_est}')
    return sigma_est
def wavelet_denoise(noisy,sigma_est,method,original,picname):
    '''
    Doing Wavelet denoise
    ---
    input: noisy pic,estimated sigma of original noisy pic, dimension reduction method,pic name
    output: plot result, bayes denoise matrix, visu denoise matrix
    '''
    fig, ax = plt.subplots(nrows=2, ncols=3, figsize=(12, 8),sharex=True, sharey=True)
    plt.gray()
    #convert2ycbcr=True(for colorful pic)
    im_bayes = denoise_wavelet(noisy,channel_axis=-1,method='BayesShrink', mode='soft', wavelet_levels=3, wavelet='sym8',rescale_sigma=True)
    im_visushrink = denoise_wavelet(noisy, channel_axis=-1,method='VisuShrink', mode='soft',sigma=sigma_est, rescale_sigma=True)


    # VisuShrink is designed to eliminate noise with high probability, but this
    # results in a visually over-smooth appearance.  Repeat, specifying a reduction
    # in the threshold by factors of 2 and 4.
    im_visushrink2 = denoise_wavelet(noisy,channel_axis=-1, method='VisuShrink', mode='soft',sigma=sigma_est/2, rescale_sigma=True)
    im_visushrink4 = denoise_wavelet(noisy, channel_axis=-1,method='VisuShrink', mode='soft',sigma=sigma_est/4, rescale_sigma=True)


    # Compute PSNR and SSIM as indications of image quality
    psnr_noisy = peak_signal_noise_ratio(original, noisy, data_range=1)
    psnr_bayes = peak_signal_noise_ratio(original, im_bayes, data_range=1)
    psnr_visushrink = peak_signal_noise_ratio(original, im_visushrink,data_range=1)
    psnr_visushrink2 = peak_signal_noise_ratio(original, im_visushrink2, data_range=1)
    psnr_visushrink4 = peak_signal_noise_ratio(original, im_visushrink4, data_range=1)
    ssm_noisy=structural_similarity(original, noisy,channel_axis=-1, data_range=1)
    ssm_bayes = structural_similarity(original, im_bayes,channel_axis=-1, data_range=1)
    ssm_visushrink = structural_similarity(original, im_visushrink,channel_axis=-1, data_range=1)
    ssm_visushrink2 = structural_similarity(original, im_visushrink2,channel_axis=-1, data_range=1)
    ssm_visushrink4 = structural_similarity(original, im_visushrink4,channel_axis=-1, data_range=1)


    ax[0, 1].imshow(noisy)
    ax[0, 1].axis('off')
    ax[0, 1].set_title(f'{method} Denoisy\nPSNR={psnr_noisy:0.4g}\nSSIM={ssm_noisy:0.4g}')
    ax[0, 2].imshow(im_bayes)
    ax[0, 2].axis('off')
    ax[0, 2].set_title(
        f'Wavelet denoising\n(BayesShrink)\nPSNR={psnr_bayes:0.4g}'
        '\nSSIM=%0.4g' % ssm_bayes)
    ax[1, 0].imshow(im_visushrink)
    ax[1, 0].axis('off')
    ax[1, 0].set_title(
        'Wavelet denoising\n(VisuShrink, $\\sigma=\\sigma_{est}$)\n'
         f'PSNR={psnr_visushrink:0.4g}\nSSIM={ssm_visushrink:0.4g}')
    ax[0, 0].imshow(original)
    ax[0, 0].axis('off')
    ax[0, 0].set_title('Original')
    ax[1, 1].imshow(im_visushrink2)
    ax[1, 1].axis('off')
    ax[1, 1].set_title(
        'Wavelet denoising\n(VisuShrink, $\\sigma=\\sigma_{est}/2$)\n'
         f'PSNR={psnr_visushrink2:0.4g}\nSSIM={ssm_visushrink2:0.4g}')
    ax[1, 2].imshow(im_visushrink4)
    ax[1, 2].axis('off')
    ax[1, 2].set_title(
        'Wavelet denoising\n(VisuShrink, $\\sigma=\\sigma_{est}/4$)\n'
         f'PSNR={psnr_visushrink4:0.4g}\nSSIM={ssm_visushrink4:0.4g}')
    fig.tight_layout()
    plt.savefig(f'{method} Denoising_{picname}.png')
    plt.show()


    return im_bayes,im_visushrink4
def proc3D(I,n,picname):
    '''
    Doing whole process on 3D pic
    '''
    plt.imshow(I)
    original=I
    sigma = 0.15
    noisy_or = random_noise(original, var=sigma**2)
    red,green,blue =sep_RGB(noisy_or)
    noisy_pca=PCA_denoising(red,green,blue,n)
    noisy_svd=SVD_denoising(red,green,blue)
    sigma_est=sigmaest(noisy_or)
    no_b,no_visu=wavelet_denoise(noisy_or,sigma_est,'No',original,picname)
    pca_b,pca_visu=wavelet_denoise(noisy_pca,sigma_est,'PCA',original,picname)
    svd_b,svd_visu=wavelet_denoise(noisy_svd,sigma_est,'SVD',original,picname)
    return no_b,no_visu,pca_b,pca_visu,svd_b,svd_visu
psnr_noisy=[]
ssm_noisy=[]
i=0
x=[]
original=I
sigma = 0.15
best_p=0
best_s=0
n_best=0
n_bests=0
noisy_or = random_noise(original, var=sigma**2)
red,green,blue =sep_RGB(noisy_or)
for n in range(20,201,10):
    x.append(n)
    pca_no=PCA_denoising(red,green,blue,n)
    p=peak_signal_noise_ratio(original, pca_no,data_range=1)
    s=structural_similarity(original, pca_no,channel_axis=-1,data_range=1)
    psnr_noisy.append(p)
    ssm_noisy.append(s)
    i+=1
    if best_p<p:
        best_p=p
        n_best=n
    if best_s<s:
        best_s=s
        n_bests=n
plt.plot(x,psnr_noisy,'s-',color = 'r',label="PSNR")#s-:方形
plt.xlabel("PCA ncomponent")
plt.ylabel("value")
plt.legend(loc = "best")
print(f'best ncomponents={n_best}')

best ncomponents=20

plt.plot(x,ssm_noisy,'d-',color = 'b',label="SSIM")#s-:方形
plt.xlabel("PCA ncomponent")
plt.ylabel("value")
plt.legend(loc = "best")
print(f'best ncomponents={n_best}')
best ncomponents=20

I=mpig.imread('wave line.jpeg')/255
no_b,no_visu,pca_b,pca_visu,svd_b,svd_visu=proc3D(I,50)

担任《Mechanical System and Signal Processing》审稿专家，担任《中国电机工程学报》，《控制与决策》等EI期刊审稿专家，擅长领域：现代信号处理，机器学习，深度学习，数字孪生，时间序列分析，设备缺陷检测、设备异常检测、设备智能故障诊断与健康管理PHM等。
知乎学术咨询：https://www.zhihu.com/consult/people/792359672131756032?isMe=1