数字图像处理-小波域维纳滤波

问题描述

编程作业:
以 lena 图像为例，编程实现小波域维纳滤波

编程实现:

step1: 导入图片

step2: 构造添加噪声的函数得到有噪声的图片

def get_noise(img, sigma):
    row, col = img.shape
    noise = np.random.randn(row,col)*sigma
    img_noise = img + noise
    img_noise = np.where(img_noise <= 0, 0, img_noise)
    img_noise = np.where(img_noise > 255, 255, img_noise)
    return np.uint8(img_noise)

噪声方差取值为 10，得到的噪声图片如下所示

step3: 构造维纳滤波的方差估计函数，其公式为:

def wiener_filter(img,HH):
    std_noise = np.median(np.abs(HH))/0.6745
    sigma_noise = std_noise**2
    sigma_img = (np.sum(img*img)-sigma_noise)/(img.shape[0]*img.shape[1])
    img_no_noise = sigma_img/(sigma_img+sigma_noise)*img
    return img_no_noise

step4: 由于小波域变换后的维纳滤波是一个递归过程，LL(低频分量)要不断的被分割为多尺度，又要被反变换回图片。故将小波域的变换和反变换过程写到一个函数里面，进行递归过程。采用的 bior4.4 小波基。

def mutial_dwt2(img, index):
    img = pywt.dwt2(img, 'bior4.4')
    LL, (LH, HL, HH) = img
    LH = wiener_filter(LH, HH)
    HL = wiener_filter(HL, HH)
    HH = wiener_filter(HH, HH)
    if index > 1:
        LL = mutial_dwt2(LL, index-1)
        row, col = LL.shape
        #小波重构可能会改变矩阵维数，进行统一
        d1 = row - HH.shape[0]
        d2 = col - HH.shape[1]
        if d1 > 0 or d2 > 0:
            d1 = row - np.arange(d1) - 1
            d2 = col - np.arange(d2) - 1
            LL = np.delete(LL, d1, axis=0)
            LL = np.delete(LL, d2, axis=1)
    img1 = pywt.idwt2((LL, (LH, HL, HH)), 'bior4.4')
    return img1

step5: 运行函数，其中分解和重构次数为 5，得到的去燥图片如下

step6: 图片对比，由于噪声值较小，所以进行了归一化显示

从结果可以看出，去出了不少的噪声，已经能看出原图像的大概轮廓。

运行说明:

环境
Numpy == 1.19.2
OpenCV-python == 4.4.0.44
PyWavelets == 1.1.1
全代码

import pywt
import cv2
import numpy as np
def normalization(img):
    img = (img - np.min(img))/(np.max(img)-np.min(img))
    return img

def get_noise(img, sigma):
    row, col = img.shape
    noise = np.random.randn(row,col)*sigma
    img_noise = img + noise
    img_noise = np.where(img_noise <= 0, 0, img_noise)
    img_noise = np.where(img_noise > 255, 255, img_noise)
    return np.uint8(img_noise)


def wiener_filter(img,HH):
    std_noise = np.median(np.abs(HH))/0.6745
    sigma_noise = std_noise**2
    sigma_img = (np.sum(img*img)-sigma_noise)/(img.shape[0]*img.shape[1])
    img_no_noise = sigma_img/(sigma_img+sigma_noise)*img
    return img_no_noise


def mutial_dwt2(img, index):
    img = pywt.dwt2(img, 'bior4.4')
    LL, (LH, HL, HH) = img
    LH = wiener_filter(LH, HH)
    HL = wiener_filter(HL, HH)
    HH = wiener_filter(HH, HH)
    if index > 1:
        LL = mutial_dwt2(LL, index-1)
        row, col = LL.shape
        #小波重构可能会改变矩阵维数，进行统一
        d1 = row - HH.shape[0]
        d2 = col - HH.shape[1]
        if d1 > 0 or d2 > 0:
            d1 = row - np.arange(d1) - 1
            d2 = col - np.arange(d2) - 1
            LL = np.delete(LL, d1, axis=0)
            LL = np.delete(LL, d2, axis=1)
    img1 = pywt.idwt2((LL, (LH, HL, HH)), 'bior4.4')
    return img1



if __name__ == '__main__':
    img = cv2.imread('lena512color.tiff',cv2.IMREAD_GRAYSCALE)
    cv2.imshow('original', img)
    img_ = get_noise(img, 10)
    cv2.imshow('noise', img_)
    img1 = mutial_dwt2(img_, 5)
    img1 = np.where(img1 < 0, 0, img1)
    img1 = np.where(img1 > 255, 255, img1)
    img1 = np.uint8(img1)
    cv2.imshow('qz', img1)
    cv2.imshow('o_n', normalization(np.uint8(np.abs(img_.astype('int16')-img.astype('int16')))))
    cv2.imshow('o_q', normalization(np.uint8(np.abs(img1.astype('int16') - img.astype('int16')))))
    cv2.waitKey(0)