数字图像处理：图像复原应用

1.1 什么是图像复原？

图像复原是图像处理中的一个重要领域，旨在从退化（例如噪声、模糊等）图像中恢复出尽可能接近原始图像的结果。图像复原与图像增强不同，复原更多地依赖于图像退化的数学模型，通过这些模型恢复出原始图像的细节。
退化模型

1.2 点扩散函数（PSF）

点扩散函数（Point Spread Function, PSF） 是描述图像退化过程中的一个关键概念，它反映了成像系统对一个点光源的响应。简单来说，点扩散函数定义了一个单一的点（如一个像素）在经过成像系统后，如何在图像平面上扩散和模糊。

1.2 逆滤波

逆滤波是图像复原中的一种经典方法，旨在通过已知的退化模型（即退化函数）反转退化过程，以恢复图像的原始内容。逆滤波假设退化过程是线性、时不变的，并且噪声影响可以忽略。因此，逆滤波主要依赖于已知的退化系统（通常通过点扩散函数（PSF）来描述），来恢复原始图像。

1.3 维纳滤波

维纳滤波是一种基于最小均方误差（MSE, Mean Square Error）准则的图像复原方法。它的目标是通过最小化复原图像与真实图像之间的均方误差来提高图像质量。与逆滤波不同，维纳滤波不仅考虑了退化过程（即模糊），还将噪声的影响纳入考虑，因此它是一种更为有效的图像复原方法。
在图像复原中，维纳滤波是一种自适应滤波器，它能够根据退化图像的统计特性进行调整，从而在去除噪声的同时尽可能保留图像细节。

1.4 代码实现

【实验内容】
（1）选取实验用的图像，完成图像的读取和显示；
（2）对图像进行退化处理，添加不同类型的噪声（高斯噪声、椒盐噪声等）；
（3）使用逆滤波、维纳滤波和其他图像复原算法进行处理，比较复原效果；
（4）评估复原后图像的质量，使用峰值信噪比（PSNR）、结构相似性（SSIM）等指标进行分析。

# （1）选取实验用的图像，完成图像的读取和显示；
# （2）对图像进行退化处理，添加不同类型的噪声（高斯噪声、椒盐噪声等）；
# （3）使用逆滤波、维纳滤波和其他图像复原算法进行处理，比较复原效果；
# （4）评估复原后图像的质量，使用峰值信噪比（PSNR）、结构相似性（SSIM）等指标进行分析。

import numpy as np
import cv2 as cv
import matplotlib.pyplot as plt
import numpy.fft as fft
import random

# 读取和显示图像
img_dir = r'D:\Document\Experiment\data\image1.jpg'
gray = cv.imread(img_dir, 0)  # 以灰度模式读取图像
image_array = np.array(gray)  # 将图像转换为 NumPy 数组

# 添加高斯噪声
def add_gaussian_noise(image, mean=0, sigma=25):
    """向图像添加高斯噪声"""
    gaussian_noise = np.random.normal(mean, sigma, image.shape)  # 生成高斯噪声
    noisy_image = image + gaussian_noise  # 图像加上噪声
    return np.clip(noisy_image, 0, 255).astype(np.uint8)  # 将像素值限制在0-255之间

# 添加椒盐噪声
def add_salt_pepper_noise(img, percentage=0.05):
    """添加椒盐噪声，percentage表示噪声比例"""
    num = int(percentage * img.size)  # 噪声点的数量
    img2 = img.copy()
    
    for i in range(num):
        x = random.randint(0, img2.shape[0] - 1)
        y = random.randint(0, img2.shape[1] - 1)
        if random.randint(0, 1) == 0:  # 随机选择是盐噪声还是椒噪声
            img2[x, y] = 255  # 盐噪声，白色
        else:
            img2[x, y] = 0  # 椒噪声，黑色
    return img2

# 逆滤波
def inverse_filter(noisy_image):
    """应用逆滤波"""
    H = np.ones(noisy_image.shape)  # 假设退化函数 H 为 1（理想情况）
    f_transform = fft.fft2(noisy_image)
    H_inv = np.where(H == 0, 0, 1 / H)  # 构造 H 的逆滤波器
    restored_fft = f_transform * H_inv  # 频域乘法
    restored_image = np.abs(fft.ifft2(restored_fft))  # 逆傅里叶变换
    return np.clip(restored_image, 0, 255).astype(np.uint8)

# 维纳滤波
def wiener_filter(noisy_image, K=0.01):
    """应用维纳滤波"""
    noisy_fft = fft.fft2(noisy_image)
    H = np.ones(noisy_image.shape)  # 假设退化函数 H(u, v) 为 1（理想情况）
    S_nn = K  # 噪声功率
    S_xx = np.var(noisy_image)  # 原图像功率（使用噪声图像近似）
    
    wiener_filter_fft = (np.conj(H) / (np.abs(H) ** 2 + S_nn / S_xx)) * noisy_fft
    restored_image = np.abs(fft.ifft2(wiener_filter_fft))
    return np.clip(restored_image, 0, 255).astype(np.uint8)

# 高斯滤波
def gaussian_smoothing(noisy_image, kernel_size=5, sigma=1):
    """高斯滤波复原图像"""
    return cv.GaussianBlur(noisy_image, (kernel_size, kernel_size), sigma)

# 计算峰值信噪比 (PSNR)
def calculate_psnr(original_image, restored_image):
    mse = np.mean((original_image - restored_image) ** 2)
    if mse == 0:
        return float('inf')
    max_pixel = 255.0
    psnr = 20 * np.log10(max_pixel / np.sqrt(mse))
    return psnr

# 计算结构相似性 (SSIM)
def calculate_ssim(original_image, restored_image):
    C1 = (0.01 * 255) ** 2
    C2 = (0.03 * 255) ** 2

    original_image = original_image.astype(np.float64)
    restored_image = restored_image.astype(np.float64)

    mu_x = np.mean(original_image)
    mu_y = np.mean(restored_image)
    sigma_x = np.var(original_image)
    sigma_y = np.var(restored_image)
    sigma_xy = np.mean((original_image - mu_x) * (restored_image - mu_y))

    ssim = ((2 * mu_x * mu_y + C1) * (2 * sigma_xy + C2)) / ((mu_x ** 2 + mu_y ** 2 + C1) * (sigma_x + sigma_y + C2))
    return ssim

# 添加高斯噪声和椒盐噪声
gaussian_noisy_image = add_gaussian_noise(image_array, sigma=45)
salt_pepper_noisy_image = add_salt_pepper_noise(image_array, percentage=0.5)

# 高斯噪声应用逆滤波、维纳滤波、高斯平滑
restored_gaussian_inverse = inverse_filter(gaussian_noisy_image)
restored_gaussian_wiener = wiener_filter(gaussian_noisy_image)
restored_gaussian_smooth = gaussian_smoothing(gaussian_noisy_image)

# 椒盐噪声应用逆滤波、维纳滤波、高斯平滑
restored_salt_inverse = inverse_filter(salt_pepper_noisy_image)
restored_salt_wiener = wiener_filter(salt_pepper_noisy_image)
restored_salt_smooth = gaussian_smoothing(salt_pepper_noisy_image)

# 计算高斯噪声复原图像质量
psnr_gaussian_inverse = calculate_psnr(image_array, restored_gaussian_inverse)
psnr_gaussian_wiener = calculate_psnr(image_array, restored_gaussian_wiener)
psnr_gaussian_smooth = calculate_psnr(image_array, restored_gaussian_smooth)

ssim_gaussian_inverse = calculate_ssim(image_array, restored_gaussian_inverse)
ssim_gaussian_wiener = calculate_ssim(image_array, restored_gaussian_wiener)
ssim_gaussian_smooth = calculate_ssim(image_array, restored_gaussian_smooth)

# 计算椒盐噪声复原图像质量
psnr_salt_inverse = calculate_psnr(image_array, restored_salt_inverse)
psnr_salt_wiener = calculate_psnr(image_array, restored_salt_wiener)
psnr_salt_smooth = calculate_psnr(image_array, restored_salt_smooth)

ssim_salt_inverse = calculate_ssim(image_array, restored_salt_inverse)
ssim_salt_wiener = calculate_ssim(image_array, restored_salt_wiener)
ssim_salt_smooth = calculate_ssim(image_array, restored_salt_smooth)

# 展示图像和计算结果
plt.figure(figsize=(12, 8))

plt.subplot(3, 3, 1)
plt.imshow(image_array, cmap='gray')
plt.title('Original Image')

plt.subplot(3, 3, 2)
plt.imshow(gaussian_noisy_image, cmap='gray')
plt.title('Gaussian Noisy Image')

plt.subplot(3, 3, 3)
plt.imshow(salt_pepper_noisy_image, cmap='gray')
plt.title('Salt Noisy Image')

plt.subplot(3, 3, 4)
plt.imshow(restored_gaussian_inverse, cmap='gray')
plt.title(f'Gaussian_Noise_Inverse_Filter\nPSNR: {psnr_gaussian_inverse:.2f}, SSIM: {ssim_gaussian_inverse:.4f}')

plt.subplot(3, 3, 5)
plt.imshow(restored_gaussian_wiener, cmap='gray')
plt.title(f'Gaussian_Noise_Wiener_Filter\nPSNR: {psnr_gaussian_wiener:.2f}, SSIM: {ssim_gaussian_wiener:.4f}')

plt.subplot(3, 3, 6)
plt.imshow(restored_gaussian_smooth, cmap='gray')
plt.title(f'Gaussian_Noise_Gaussian_Smoothing\nPSNR: {psnr_gaussian_smooth:.2f}, SSIM: {ssim_gaussian_smooth:.4f}')

plt.subplot(3, 3, 7)
plt.imshow(restored_salt_inverse, cmap='gray')
plt.title(f'Salt_Noise_Inverse_Filter\nPSNR: {psnr_salt_inverse:.2f}, SSIM: {ssim_salt_inverse:.4f}')

plt.subplot(3, 3, 8)
plt.imshow(restored_salt_wiener, cmap='gray')
plt.title(f'Salt_Noise_Wiene_Filter\nPSNR: {psnr_salt_wiener:.2f}, SSIM: {ssim_salt_wiener:.4f}')

plt.subplot(3, 3, 9)
plt.imshow(restored_salt_smooth, cmap='gray')
plt.title(f'Salt_Noise_Smoothing\nPSNR: {psnr_salt_smooth:.2f}, SSIM: {ssim_salt_smooth:.4f}')

plt.tight_layout()
plt.show()