修复（去噪与反卷积） - 计算机视觉算法与实践-从入门到精通

本节简要回顾图像修复：利用噪声/模糊模型和正则化，从退化的观测中按原则恢复清晰图像。

概述¶

修复将图像恢复问题表述为一个反问题。给定观测模型（如噪声、点扩散函数）和先验/正则化项，我们估计能够最好地解释观测数据的图像。

去噪¶

去噪在于减少图像中的噪声。注意，通常无法完全消除噪声。本节我们首先列出最常见的噪声模型，然后介绍一些去噪方法。

# 配置中文字体以避免 Matplotlib 中文缺字警告
import matplotlib.pyplot as plt

# 按优先级尝试不同的中文字体（同 fourier.ipynb 方法）
plt.rcParams['font.sans-serif'] = [
    'Noto Sans CJK SC', 'Noto Sans CJK JP', 'SimHei',
    'Microsoft YaHei', 'WenQuanYi Micro Hei', 'DejaVu Sans'
]
plt.rcParams['axes.unicode_minus'] = False

# Prepare a base image and noisy variants so downstream cells can run standalone
import numpy as np
import matplotlib.pyplot as plt
from skimage import data
from skimage.util import random_noise

# Load a sample grayscale image
image = data.camera()

# Create noisy versions used in later examples
gaussian_noisy = random_noise(image, mode='gaussian', var=0.05)
poisson_noisy = random_noise(image, mode='poisson')
salt_pepper_noisy = random_noise(image, mode='s&p', amount=0.05)

噪声来源¶

数字图像中的主要噪声来源是在采集过程中（收集到的光子数量太少、传感器温度……）或在任何传输过程中（无线通信中的回声和大气失真）。在某些情况下，噪声也被用来模拟图像形成数学模型中的不准确性，因为后者与现实相比必然存在差异，就像任何物理模型一样！

噪声本质上是一种随机现象，通过表示其强度分布的概率密度进行建模。

在下文中，我们用 $y$ 表示观测到的（带噪声的）图像， $b$ 表示噪声， $x$ 表示无噪声的图像。

加性高斯白噪声¶

加性高斯白噪声（AWGN）将观测图像 $y$ 的每个像素 $(m,n)$ 建模为无噪声图像 $x$ 的像素 $(m,n)$ 与噪声 $b$ 的一个像素之和：

\forall\, m,n \quad y(m,n) = x(m,n) + b(m,n)

(1)

其中 $b(m,n) \sim \mathcal{N}(0,\sigma^2)$ 。

这个模型很简单，便于计算。它被用于包括摄影在内的大多数应用中。

泊松噪声¶

泊松噪声（也称散粒噪声）模拟光子在光敏元件上的采集过程。光子的数量是随机的，并取决于光照强度。相应的泊松过程的均值等于光照强度。观测图像 $y$ 的每个像素 $(m,n)$ 的强度为：

\forall\, m,n \qquad y(m,n) \sim \mathcal{P}\big(x(m,n)\big).

(2)

该模型用于光子数量较少的采集情况，例如在天文学中。

椒盐噪声¶

椒盐噪声，也不太诗意地称为脉冲噪声，模拟了饱和或失效的像素（由于光敏元件故障或饱和）。

\forall\, m,n \quad y(m,n) = \begin{cases} x_\mathrm{min} &\text{概率为}\,p_\mathrm{min}, \\ x_\mathrm{max} &\text{概率为}\,p_\mathrm{max}, \\ x(m,n) &\text{概率为}\,1\!-\!p_\mathrm{min}\!-\!p_\mathrm{max}. \end{cases}

(4)

其中 $x_\mathrm{min}$ 和 $x_\mathrm{max}$ 是强度的最小值和最大值。

# 用于设置依赖和显示噪声示例的代码块

# 如果尚未安装scikit-image，则安装
try:
    import skimage
except ImportError:
    # %pip install scikit-image
    print("scikit-image not installed. Please install it to run this cell.")

import numpy as np
import matplotlib.pyplot as plt
from skimage import data
from skimage.util import random_noise

# 加载示例图像
image = data.camera()
# --- 向图像添加噪声 ---
# 高斯噪声
gaussian_noisy = random_noise(image, mode='gaussian', var=0.05)
# 泊松噪声
poisson_noisy = random_noise(image, mode='poisson')
# 椒盐噪声
salt_pepper_noisy = random_noise(image, mode='s&p', amount=0.05)
# --- 显示结果 ---
fig, axes = plt.subplots(1, 4, figsize=(20, 5))
ax = axes.ravel()

ax[0].imshow(image, cmap='gray')
ax[0].set_title('原始图像')
ax[0].axis('off')

ax[1].imshow(gaussian_noisy, cmap='gray')
ax[1].set_title('高斯噪声')
ax[1].axis('off')

ax[2].imshow(poisson_noisy, cmap='gray')
ax[2].set_title('泊松噪声')
ax[2].axis('off')

ax[3].imshow(salt_pepper_noisy, cmap='gray')
ax[3].set_title('椒盐噪声')
ax[3].axis('off')

plt.tight_layout()
plt.show()

噪声功率¶

Figure 2 说明了先前噪声在同一图像上的效果。可以观察到：

对于高斯噪声，整个图像受到同样方式的噪声影响，
对于泊松噪声，较亮的部分比较暗的部分噪声更大，
对于脉冲噪声，只有少数像素被修改，并被替换为黑色或白色像素。

信噪比（SNR）是衡量噪声水平的指标。它定义为无噪声图像的功率与噪声功率之比，其中图像 $x$ 的功率定义为：

P_x = \frac{1}{M \times N} \sum_{m,n} x(m,n)^2

(5)

因为SNR通常以对数刻度（单位：分贝）表示，所以它也定义为：

\text{SNR} = 10 \log_{10} \left( \frac{\sum_{m,n} x(m,n)^2}{\sum_{m,n} b(m,n)^2} \right)

(6)

对于加性噪声，存在另一个度量：峰值信噪比（PSNR）是无噪声图像的平方动态范围（最大强度与最小强度之差）与噪声功率之比：

\text{PSNR} = 10 \log_{10} \left( \frac{\Delta x^{\;2}}{\frac{1}{M \times N} \sum_{m,n} b(m,n)^2} \right)

(7)

Figure 3 表示在不同SNR和PSNR下，被加性高斯白噪声污染的同一图像。如您所见，当SNR或PSNR增加时，噪声会减小！

from skimage.metrics import peak_signal_noise_ratio

def signal_to_noise_ratio(true_image, noisy_image):
    """计算信噪比。"""
    # 确保图像为浮点类型以便计算
    true_image = true_image.astype(np.float64)
    noisy_image = noisy_image.astype(np.float64)
    
    # 计算噪声
    noise = true_image - noisy_image
    
    # 计算信号和噪声的功率
    signal_power = np.sum(true_image**2)
    noise_power = np.sum(noise**2)
    
    # 计算信噪比
    if noise_power == 0:
        return float('inf')
    
    snr = 10 * np.log10(signal_power / noise_power)
    return snr

# --- 计算并打印带噪图像的SNR和PSNR ---
# 注意：scikit-image的random_noise会将图像缩放到[0, 1]或[-1, 1]
# 因此我们需要注意数据范围。我们将使用skimage的PSNR
# 来处理这个问题。对于SNR，我们的自定义函数是合适的。
# 原始图像是uint8，random_noise会将其转换为[0,1]范围内的float64。
image_float = image / 255.0

psnr_gaussian = peak_signal_noise_ratio(image_float, gaussian_noisy)
snr_gaussian = signal_to_noise_ratio(image_float, gaussian_noisy)

# 对于泊松噪声和椒盐噪声，噪声不完全是加性的，但我们仍然可以计算这些指标。
psnr_poisson = peak_signal_noise_ratio(image_float, poisson_noisy)
snr_poisson = signal_to_noise_ratio(image_float, poisson_noisy)

psnr_sp = peak_signal_noise_ratio(image_float, salt_pepper_noisy)
snr_sp = signal_to_noise_ratio(image_float, salt_pepper_noisy)

print(f"高斯噪声 -> SNR: {snr_gaussian:.2f} dB, PSNR: {psnr_gaussian:.2f} dB")
print(f"泊松噪声  -> SNR: {snr_poisson:.2f} dB, PSNR: {psnr_poisson:.2f} dB")
print(f"椒盐噪声      -> SNR: {snr_sp:.2f} dB, PSNR: {psnr_sp:.2f} dB")

高斯噪声 -> SNR: 9.52 dB, PSNR: 14.21 dB
泊松噪声  -> SNR: 22.37 dB, PSNR: 27.06 dB
椒盐噪声      -> SNR: 13.04 dB, PSNR: 17.73 dB

均值滤波器¶

均值滤波器是一种非常简单的去噪方法。去噪图像 $\widehat{x}$ 的每个像素 $(m,n)$ 是带噪图像 $y$ 中 $(m,n)$ 周围像素的平均值：

\forall m,\,n \quad \widehat{x}(m,n) = \frac{1}{|V_{m,n}|} \sum_{(u,v)\in V_{m,n}} y(u,v)

(8)

其中

$V_{m,n}$ 是邻域，即 $(m,n)$ 周围的像素集；
$|V_{m,n}|$ 是 $V_{m,n}$ 的基数，即邻域中的像素数。

Figure 4 说明了均值滤波器对于不同邻域大小的效果。如果邻域增大，则噪声减小，但同时图像变得更模糊。

均值滤波器可以用卷积积表示。确实，考虑邻域是大小为 $N\times N$ 像素的正方形的情况，那么均值滤波器的定义给出：

\widehat{x}(m,n) := \frac{1}{N^2} \sum_{(u,v)\in V_{m,n}} y(u,v) := \sum_{u,v} g(m-u,n-v) y(u,v)

(9)

其中

g(u,v) = \begin{cases} 1/N^2 &\text{if}\, u\in\left\{-\frac{N}{2},\dots,\frac{N}{2}\right\} \,\text{and}\, v\in\left\{-\frac{N}{2},\dots,\frac{N}{2}\right\} \\ 0 &\text{otherwise} \end{cases}

(10)

这个定义可以扩展到任何类型的核 $g$ ！例如，Figure 5 给出了两种不同核的结果。

总结：

均值滤波器计算邻域中像素的平均值，
均值滤波器可以写成卷积，
通过平均强度来减少噪声，但图像会变模糊。

from skimage.filters import rank
from skimage.morphology import disk
from skimage.util import img_as_ubyte

# rank 滤波器期望 8-bit 输入；将浮点噪声图转为 uint8
gaussian_noisy_u8 = img_as_ubyte(gaussian_noisy)

# 应用不同核大小的均值滤波器（uint8 输入/输出）
mean_denoised_3 = rank.mean(gaussian_noisy_u8, footprint=disk(3))
mean_denoised_7 = rank.mean(gaussian_noisy_u8, footprint=disk(7))

# --- 显示结果 ---
fig, axes = plt.subplots(1, 3, figsize=(15, 5))
ax = axes.ravel()

ax[0].imshow(gaussian_noisy, cmap='gray')
ax[0].set_title('带高斯噪声的图像')
ax[0].axis('off')

ax[1].imshow(mean_denoised_3, cmap='gray')
ax[1].set_title('均值滤波器 (3x3 核)')
ax[1].axis('off')

ax[2].imshow(mean_denoised_7, cmap='gray')
ax[2].set_title('均值滤波器 (7x7 核)')
ax[2].axis('off')

plt.tight_layout()
plt.show()

# --- 计算并打印去噪图像的PSNR ---
mean3_float = mean_denoised_3.astype(np.float64) / 255.0
mean7_float = mean_denoised_7.astype(np.float64) / 255.0
psnr_mean_3 = peak_signal_noise_ratio(image_float, mean3_float, data_range=1.0)
psnr_mean_7 = peak_signal_noise_ratio(image_float, mean7_float, data_range=1.0)
print(f"PSNR (均值 3x3): {psnr_mean_3:.2f} dB")
print(f"PSNR (均值 7x7): {psnr_mean_7:.2f} dB")

PSNR (均值 3x3): 23.27 dB
PSNR (均值 7x7): 21.83 dB

中值滤波器¶

一组数的中值是集合中的元素 $m$ ，使得小于 $m$ 的数和大于 $m$ 的数一样多。例如， $\{1,\,2,\,4,\,8,\,16\}$ 的中值是 4。

中值滤波器定义为：

\forall m,\,n \quad \widehat{x}(m,n) = \mathrm{median}\big(\{y(u,v) \mid (u,v)\in V_{m,n}\}\big)

(11)

中值滤波器非常适合用于去除椒盐噪声，因为它不会像均值滤波器那样使图像模糊。

尽管名为滤波器，但中值滤波器不是一个线性滤波器，因为它不满足线性性质。因此它不能写成卷积的形式。

# --- 比较均值和中值滤波器对椒盐噪声的效果 ---
from skimage.util import img_as_ubyte

# rank.mean 期望 uint8；先将输入转换为 uint8
sp_u8 = img_as_ubyte(salt_pepper_noisy)
mean_denoised_sp = rank.mean(sp_u8, footprint=disk(3))

# 中值滤波器也接受 uint8；保持不变
median_denoised_sp = rank.median(sp_u8, footprint=disk(3))

# --- 显示结果 ---
fig, axes = plt.subplots(1, 4, figsize=(20, 5))
ax = axes.ravel()

ax[0].imshow(image, cmap='gray')
ax[0].set_title('原始图像')
ax[0].axis('off')

ax[1].imshow(salt_pepper_noisy, cmap='gray')
ax[1].set_title('带椒盐噪声的图像')
ax[1].axis('off')

ax[2].imshow(mean_denoised_sp, cmap='gray')
ax[2].set_title('均值滤波器去噪')
ax[2].axis('off')

ax[3].imshow(median_denoised_sp, cmap='gray')
ax[3].set_title('中值滤波器去噪')
ax[3].axis('off')

plt.tight_layout()
plt.show()

# --- 计算并打印去噪图像的PSNR ---
mean_sp_float = mean_denoised_sp.astype(np.float64) / 255.0
median_sp_float = median_denoised_sp.astype(np.float64) / 255.0
psnr_mean_sp = peak_signal_noise_ratio(image_float, mean_sp_float, data_range=1.0)
psnr_median_sp = peak_signal_noise_ratio(image_float, median_sp_float, data_range=1.0)
print(f"PSNR (均值滤波器对椒盐噪声):   {psnr_mean_sp:.2f} dB")
print(f"PSNR (中值滤波器对椒盐噪声): {psnr_median_sp:.2f} dB")

PSNR (均值滤波器对椒盐噪声):   25.03 dB
PSNR (中值滤波器对椒盐噪声): 27.63 dB

周期性噪声滤波¶

周期性噪声的特点是在傅里叶变换中呈现为结构性图案。这些结构可以通过在傅里叶域中应用一个掩模来消除相应的系数，从而被移除。去噪后的图像随后通过逆傅里叶变换得到。

from scipy import fftpack

# --- 创建周期性噪声 ---
# 创建一个带周期性噪声的合成图像
rows, cols = image.shape
crow, ccol = rows // 2 , cols // 2
u, v = 30, 50 # 噪声的频率分量

# 创建正弦噪声
noise = np.zeros((rows, cols), dtype=np.float64)
for r in range(rows):
    for c in range(cols):
        noise[r, c] = 20 * (np.sin(2 * np.pi * u * r / rows) + np.sin(2 * np.pi * v * c / cols))

# 将噪声添加到图像中
periodic_noisy = np.clip(image_float + noise / 255.0, 0, 1)

# --- 傅里叶变换与滤波 ---
# 计算傅里叶变换
f_transform = fftpack.fft2(periodic_noisy)
f_shift = fftpack.fftshift(f_transform)
magnitude_spectrum = 20 * np.log(np.abs(f_shift))

# 创建一个掩模以移除周期性噪声
mask = np.ones_like(f_shift, dtype=bool)
# 噪声尖峰以水平和垂直轴为中心呈对称分布，因此我们将它们掩蔽掉
mask[crow - u - 2 : crow - u + 2, ccol - 2 : ccol + 2] = False
mask[crow + u - 2 : crow + u + 2, ccol - 2 : ccol + 2] = False
mask[crow - 2 : crow + 2, ccol - v - 2 : ccol - v + 2] = False
mask[crow - 2 : crow + 2, ccol + v - 2 : ccol + v + 2] = False

# 应用掩模
f_shift_filtered = f_shift * mask
magnitude_spectrum_filtered = 20 * np.log(np.abs(f_shift_filtered) + 1) # 加1以避免log(0)

# --- 逆傅里叶变换 ---
f_ishift = fftpack.ifftshift(f_shift_filtered)
image_filtered = fftpack.ifft2(f_ishift)
image_filtered = np.abs(image_filtered)


# --- 显示结果 ---
fig, axes = plt.subplots(2, 2, figsize=(10, 10))
ax = axes.ravel()

ax[0].imshow(periodic_noisy, cmap='gray')
ax[0].set_title('带周期性噪声的图像')
ax[0].axis('off')

ax[1].imshow(magnitude_spectrum, cmap='gray')
ax[1].set_title('幅度谱')
ax[1].axis('off')

ax[2].imshow(magnitude_spectrum_filtered, cmap='gray')
ax[2].set_title('滤波后的幅度谱')
ax[2].axis('off')

ax[3].imshow(image_filtered, cmap='gray')
ax[3].set_title('去噪后的图像')
ax[3].axis('off')

plt.tight_layout()
plt.show()

TV正则化¶

从某种角度看，去噪的目标是获得一个图像 $\widehat{x}$ ，它不仅像素间的强度变化小，而且接近观测值 $y$ 。 TV正则化（total variation） [Rudin et al. 1992, Chambolle 2004] 是一种去噪方法，它通过数学函数（即所谓的“准则”）来描述这两个目标。

希望去噪图像 $\widehat{x}$ 接近观测值 $y$ 的愿望，产生了一个“数据拟合”准则，该准则衡量 $x$ 和 $y$ 之间的差异。一个经典的选择是最小二乘准则：
$E(x,y) = \sum_{m,n} \left(y(m,n)-x(m,n)\right)^2$
(12)
确实，为了使 $E$ 减小，我们必须找到一个接近 $y$ 的图像 $x$ 。
希望图像强度变化小的愿望，产生了一个“正则化”准则，该准则衡量 $x$ 的相邻像素之间的差异。一个简单的选择是“总变分”：
$R(x) = \sum_{m,n} \left|x(m+1,n)-x(m,n)\right| + \sum_{m,n} \left|x(m,n+1)-x(m,n)\right|$
(13)
因此，为了使 $R$ 减小，无论是行还是列，两个连续像素之间的差异都必须很小。这意味着图像 $x$ 的强度几乎是恒定的。

目标是找到一个既能最小化数据拟合又能最小化正则化的图像 $x$ 。在数学上，我们寻找最小化 $E(x,y) + \lambda R(x)$ 的图像 $x$ ，其中 $\lambda$ 是“正则化参数”，用于调整两个准则之间的权衡。 $\lambda$ 的值由用户选择。在数学上，我们写成：

\widehat{x} = \arg\min_x E(x,y) + \lambda R(x)

(14)

这变成了一个优化问题，有大量的算法可以最小化 $E(x,y) + \lambda R(x)$ 。

使用TV正则化去噪，对应两个不同的 \lambda 值。 — Figure 8:使用TV正则化去噪，对应两个不同的 $\lambda$ 值。

from skimage.restoration import denoise_tv_chambolle

# --- 应用总变分 (TV) 去噪 ---
# 我们使用带高斯噪声的图像作为输入
tv_denoised_low = denoise_tv_chambolle(gaussian_noisy, weight=0.1)
tv_denoised_high = denoise_tv_chambolle(gaussian_noisy, weight=0.2)


# --- 显示结果 ---
fig, axes = plt.subplots(1, 4, figsize=(20, 5))
ax = axes.ravel()

ax[0].imshow(image, cmap='gray')
ax[0].set_title('原始图像')
ax[0].axis('off')

ax[1].imshow(gaussian_noisy, cmap='gray')
ax[1].set_title('带高斯噪声的图像')
ax[1].axis('off')

ax[2].imshow(tv_denoised_low, cmap='gray')
ax[2].set_title('TV 去噪 (weight=0.1)')
ax[2].axis('off')

ax[3].imshow(tv_denoised_high, cmap='gray')
ax[3].set_title('TV 去噪 (weight=0.2)')
ax[3].axis('off')

plt.tight_layout()
plt.show()

# --- 计算并打印去噪后图像的PSNR ---
psnr_tv_low = peak_signal_noise_ratio(image_float, tv_denoised_low)
psnr_tv_high = peak_signal_noise_ratio(image_float, tv_denoised_high)
print(f"PSNR (TV 去噪, weight=0.1): {psnr_tv_low:.2f} dB")
print(f"PSNR (TV 去噪, weight=0.2): {psnr_tv_high:.2f} dB")

PSNR (TV 去噪, weight=0.1): 22.21 dB
PSNR (TV 去噪, weight=0.2): 22.54 dB

反卷积¶

通常，由视觉系统采集的图像会受到退化，这种退化可以建模为卷积。例如，一些图像呈现出相机抖动效果（Figure 9）或由于对焦不良导致的模糊（Figure 10）。反卷积的目标是消除卷积的影响。

退化现象的模型如Figure 11所示：观测图像 $y$ 受到与PSF $h$ 的卷积以及可能的噪声 $b$ （假设为加性噪声）的退化。

y = h*x + b

(15)

从观测图像 $y$ 估计原始图像 $\widehat{x}$ 的过程称为反卷积。这是一个经典的 反问题。

挑战在于卷积和加性噪声都破坏了原始信号。恢复 $x$ 的难度在很大程度上取决于我们已知的信息：

非盲反卷积：我们假设PSF $h$ 是已知的。
- 如果 没有噪声 ( $b=0$ )，问题简化为 $y = h*x$ 。理论上，我们可以通过“反转”卷积来恢复 $x$ 。这就是 逆滤波器 背后的原理。
- 如果 存在噪声 ( $b \neq 0$ )，简单地反转卷积会急剧放大噪声，导致结果很差。这需要更复杂的方法，如 维纳滤波器，它在去模糊和噪声抑制之间取得平衡。
盲反卷积：这是最具挑战性的情况，原始图像 $x$ 和 PSF $h$ 都是未知的。这超出了我们当前的讨论范围，但它是一个活跃的研究领域。

在接下来的部分中，我们将专注于非盲反卷积，假设PSF $h$ 是已知的或可以被估计，并且我们有一个噪声模型。

逆滤波器¶

逆滤波器是最简单的反卷积方法。由于退化模型为 $y = h*x + b$ ，该方程在傅里叶域中变为：

Y = HX + B

(16)

因此我们可以写出：

X = \frac{Y-B}{H}.

(17)

我们通过计算前一个表达式的逆傅里叶变换来得到 $x$ ：

x = \mathcal{F}^{-1} \left[ \frac{Y-B}{H} \right].

(18)

由于噪声（及其频谱 $B$ ）是未知的，我们可以通过在前一个表达式中将 $B$ 置零来近似 $x$ ，从而得到反卷积图像：

\widehat{x} = \mathcal{F}^{-1} \left[ \frac{Y}{H} \right]

(19)

逆滤波器应用于图像的结果在 Figure 12 中给出。结果是不可用的，然而观察到的图像模糊很小，噪声也很少！

逆滤波器的灾难性结果是由于假设噪声为零。确实，根据 $\widehat{x}$ 的定义并考虑 $Y = HX + B$ ，则：

\widehat{x} := \mathcal{F}^{-1} \left[ \frac{Y}{H} \right] := \mathcal{F}^{-1} \left[ X + \frac{B}{H} \right] := x + \mathcal{F}^{-1} \left[ \frac{B}{H} \right]

(20)

因此，反卷积图像 $\widehat{x}$ 对应于 $x$ 加上一个附加项，即 $B/H$ 的逆傅里叶变换。 PSF $H$ 通常是一个低通滤波器，因此对于高频 $(m,n)$ ， $H(m,n)$ 的值趋于 0。因为 $H$ 在分母中，这会急剧放大噪声的高频，然后 $B/H$ 项很快就主导了 $X$ 。这解释了 Figure 12 的结果。

一个解决方案是只考虑 $Y/H$ 的低频部分。这相当于在计算逆傅里叶变换之前，通过将高频置零来截断逆滤波器的结果。反卷积的结果变得可以接受得多，如 Figure 13 所示，尽管结果仍然远非完美（物体周围存在许多强度变化，例如树干）...

维纳滤波器¶

维纳滤波器，记为 $g$ （其傅里叶变换为 $G$ ），应用于观测值 $y$ 使得：

\widehat{x} = g * y \qquad\Leftrightarrow\qquad \widehat{X} = GY.

(21)

该滤波器是在统计框架下建立的：图像 $x$ 和噪声 $b$ 被视为随机变量。它们也被假定为统计上独立的。因此，观测值 $y$ 和估计值 $\widehat{x}$ 也是随机变量。

为简单起见，计算在傅里叶域中进行（因为卷积变为乘法）。维纳滤波器的目标是找到最接近 $X = \mathcal{F}[x]$ 的图像 $\widehat{X} = \mathcal{F}[\widehat{x}]$ ，在均方误差 $\mathrm{MSE} = \mathbb{E}\left[(\widehat{X}-X)^2\right]$ 的意义上。因此：

\begin{align} \mathrm{MSE} &= \mathbb{E}\left[(\widehat{X}-X)^2\right] \\ &= \mathbb{E}\left[(GY-X)^2\right] \\ &= \mathbb{E}\left[\big(G(HX+B)-X\big)^2\right] \\ &= \mathbb{E}\left[\big((GH-I)X+GB\big)^2\right] \end{align}

(22)

其中 $I$ 是一个全为1的图像。所以：

\mathrm{MSE} = \mathbb{E}\Big[ (GH-1)^*(GH-1)X^*X + (GH-1)^*GX^*B + G^*(GH-1)B^*X + G^*GB^*B \Big]

(23)

其中 $\cdot^*$ 表示变量的共轭。由于期望 $\mathbb{E}$ 是线性的，并且只有 $X$ 和 $B$ 是随机变量，我们可以将前一个表达式分解为四个项：

\begin{align} \mathrm{MSE} &= (GH-I)^*(GH-I) \,\mathbb{E}\big[X^*X\big]\\ &\quad + (GH-I)^*G \,\mathbb{E}\big[X^*B\big]\\ &\quad + G^*(GH-I) \,\mathbb{E}\big[B^*X\big]\\ &\quad + G^*G \,\mathbb{E}\big[B^*B\big]. \end{align}

(24)

由于 $X$ 和 $B$ 是独立的，协方差 $\mathbb{E}\big[X^*B\big]$ 和 $\mathbb{E}\big[B^*X\big]$ 为零。此外， $\mathbb{E}\big[X^*X\big]$ 和 $\mathbb{E}\big[B^*B\big]$ 是功率谱密度，分别记为 $S_x$ 和 $S_b$ （功率谱密度是傅里叶变换模的平方的期望）。所以均方误差简化为：

\mathrm{MSE} = (GH-1)^*(GH-1) S_x + G^*G S_b

(25)

我们寻找最小化MSE的滤波器 $G$ ，或者等价地，使MSE的导数等于零：

\begin{align} \frac{\partial \mathrm{MSE}}{\partial G} &= (GH-1)^*H S_x + G^* S_b = 0 \\ \Leftrightarrow\quad G^*H^*H S_x - H S_x + G^* S_b &= 0 \\ \Leftrightarrow\quad G^* ( H^*H S_x + S_b) &= H S_x \\ \Leftrightarrow\quad G^* &= \frac{H S_x}{H^*H S_x + S_b} \\ \Leftrightarrow\quad G &= \frac{H^* S_x}{H^*H S_x + S_b} \\ \Leftrightarrow\quad G &= \frac{H^* S_x}{|H|^2 S_x + S_b} \end{align}

(26)

我们得到了维纳滤波器 $G$ 的表达式！最后，反卷积图像是 $GY$ 的逆傅里叶变换：

\widehat{x} = \mathcal{F}^{-1} \Bigg[\frac{H^* S_x}{|H|^2 S_x + S_b} Y \Bigg]

(27)

我们可以认为功率谱密度 $S_x$ 和 $S_b$ 是常数（对于 $S_b$ ，需要假设是白噪声）。因此，维纳滤波器的表达式可以写成

\widehat{x} = \mathcal{F}^{-1} \Bigg[ \frac{H^*}{|H|^2 + S_b/S_x} Y \Bigg]

(28)

而 $S_b/S_x$ 项被一个常数 $K$ 代替，它成为该方法的参数，由用户设置。

两点说明：

当 $H$ 趋于零时（通常在高频），不再出现像逆滤波器那样的噪声放大问题，因为滤波器趋于0，
此外，如果图像中没有噪声，则 $S_b = 0$ ，维纳滤波器变回逆滤波器：
$\widehat{x} = \mathcal{F}^{-1} \Bigg[ \frac{H^*}{|H|^2} Y \Bigg] = \mathcal{F}^{-1} \Bigg[ \frac{Y}{H} \Bigg]$
(29)

维纳滤波器的结果在 Figure 14 中展示：它明显优于逆滤波器，甚至是其截断版本！

维纳滤波器的结果
（\lambda 被设置为使得估计在MSE方面是最佳的）。 — Figure 14:维纳滤波器的结果（ $\lambda$ 被设置为使得估计在MSE方面是最佳的）。

import numpy as np
import matplotlib.pyplot as plt

# ====================== 主脚本 ======================

def main():
    """
    主函数，用于运行图像恢复流程。
    它会加载一张图像，应用不同程度的运动模糊和噪声，
    然后尝试使用逆滤波器和维纳滤波器恢复图像。
    最后，它会计算并显示性能指标（MSE, PSNR, SSIM），
    并展示处理后的图像。
    """
    # ---------- 1. 加载并准备图像 ----------
    try:
        # 使用在前一个notebook单元格中预加载的'image'变量。
        img = image
        img = np.array(img, dtype=np.float64)
    except NameError:
        print("未找到图像'image'！请确保它已在前一个单元格中加载。")
        return
    except Exception as e:
        print(f"加载图像时发生错误: {e}")
        return

    # ---------- 2. 转换为灰度图并归一化 ----------
    # 如果图像是彩色的（3个维度），则使用标准的亮度公式将其转换为灰度图。
    if img.ndim == 3:
        R, G, B = img[:, :, 0], img[:, :, 1], img[:, :, 2]
        I = 0.2989 * R + 0.5870 * G + 0.1140 * B
    else:
        I = img
    # 将图像强度归一化到[0, 1]范围，以便进行一致的处理。
    I = (I - I.min()) / (I.max() - I.min())
    M, N = I.shape

    # ---------- 3. 定义测试用例 ----------
    # 定义了四种场景，用于在不同的模糊（L: 运动长度）和噪声（noise_sigma）
    # 条件下测试恢复算法。
    cases = [
        {"name": "低噪声, 低模糊",  "L": 5,  "theta": 30, "noise_sigma": 0.003},
        {"name": "低噪声, 高模糊", "L": 25, "theta": 30, "noise_sigma": 0.003},
        {"name": "高噪声, 低模糊", "L": 5,  "theta": 30, "noise_sigma": 0.02},
        {"name": "高噪声, 高模糊","L": 25, "theta": 30, "noise_sigma": 0.02},
    ]

    # ---------- 4. 处理每一种情况 ----------
    results = []  # 存储用于最终报告的指标。
    plt.figure(figsize=(15, 12))
    plot_idx = 1

    for case in cases:
        # --- a. 为当前情况设置参数 ---
        L = case["L"]                # 运动模糊的长度
        theta = case["theta"]            # 运动模糊的角度
        noise_sigma = case["noise_sigma"]  # 高斯噪声的标准差
        inv_eps = 1e-3               # 用于正则化逆滤波器的Epsilon

        # 估算维纳滤波器的噪声信号功率比（K）。
        img_var = np.var(I)
        noise_var = noise_sigma ** 2
        K = noise_var / max(img_var, 1e-10) # 避免对平坦图像做除零操作

        # --- b. 创建运动模糊传递函数 H ---
        # H_c 是运动模糊的中心化频率响应。
        H_c = motion_H(M, N, L, theta)
        # H 经过移位以与numpy的FFT兼容，后者期望零频分量位于角落。
        H = np.fft.ifftshift(H_c)

        # --- c. 应用模糊并添加噪声 ---
        F = np.fft.fft2(I)      # 原始图像的FFT
        G_clean = H * F         # 在频域中应用模糊
        g_clean = np.real(np.fft.ifft2(G_clean)) # 空间域中的模糊图像
        # 向模糊图像添加合成的高斯噪声。
        g = g_clean + noise_sigma * np.random.randn(M, N)

        # --- d. 使用逆滤波器进行恢复 ---
        G = np.fft.fft2(g)      # 退化图像的FFT
        Fhat_inv = inverse_filter(G, H, inv_eps) # 应用逆滤波器
        I_inv = np.real(np.fft.ifft2(Fhat_inv))   # 转换回空间域
        I_inv = np.clip(I_inv, 0, 1)          # 将值裁剪到有效的[0, 1]范围

        # --- e. 使用维纳滤波器进行恢复 ---
        Fhat_wiener = wiener_filter(G, H, K)    # 应用维纳滤波器
        I_wiener = np.real(np.fft.ifft2(Fhat_wiener)) # 转换回空间域
        I_wiener = np.clip(I_wiener, 0, 1)       # 将值裁剪到有效的[0, 1]范围

        # --- f. 计算性能指标 ---
        # 将原始图像与退化和恢复后的版本进行比较。
        mse_blur = mse(I, g)
        mse_inv = mse(I, I_inv)
        mse_wiener = mse(I, I_wiener)

        psnr_blur = psnr(I, g)
        psnr_inv = psnr(I, I_inv)
        psnr_wiener = psnr(I, I_wiener)

        ssim_blur = ssim_manual(I, g)
        ssim_inv = ssim_manual(I, I_inv)
        ssim_wiener = ssim_manual(I, I_wiener)

        # 存储此案例的所有指标。
        results.append({
            "case": case["name"],
            "mse_blur": mse_blur, "psnr_blur": psnr_blur, "ssim_blur": ssim_blur,
            "mse_inv": mse_inv, "psnr_inv": psnr_inv, "ssim_inv": ssim_inv,
            "mse_wiener": mse_wiener, "psnr_wiener": psnr_wiener, "ssim_wiener": ssim_wiener
        })

        # --- g. 绘制结果以进行视觉比较 ---
        # 对于每种情况，显示：退化图像、逆滤波器结果、维纳滤波器结果。
        plt.subplot(4, 3, plot_idx); plt.imshow(g, cmap='gray', vmin=0, vmax=1)
        plt.title(f"{case['name']}\n模糊+噪声"); plt.axis('off')
        plt.subplot(4, 3, plot_idx+1); plt.imshow(I_inv, cmap='gray', vmin=0, vmax=1)
        plt.title("逆滤波器"); plt.axis('off')
        plt.subplot(4, 3, plot_idx+2); plt.imshow(I_wiener, cmap='gray', vmin=0, vmax=1)
        plt.title("维纳滤波器"); plt.axis('off')
        plot_idx += 3

    plt.tight_layout()
    plt.show()

    # ---------- 5. 以表格形式打印量化指标 ----------
    print("\n" + "=" * 80)
    print(" 所有四种情况的恢复性能指标 ")
    print("=" * 80)
    print(f"\n{'情况':<28}{'方法':<18}{'MSE':<12}{'PSNR(dB)':<12}{'SSIM':<10}")
    print("-" * 80)
    for r in results:
        print(f"{r['case']:<28}{'模糊+噪声':<18}{r['mse_blur']:<12.5f}{r['psnr_blur']:<12.3f}{r['ssim_blur']:<10.4f}")
        print(f"{'':<28}{'逆滤波器':<18}{r['mse_inv']:<12.5f}{r['psnr_inv']:<12.3f}{r['ssim_inv']:<10.4f}")
        print(f"{'':<28}{'维纳滤波器':<18}{r['mse_wiener']:<12.5f}{r['psnr_wiener']:<12.3f}{r['ssim_wiener']:<10.4f}")
        print("-" * 80)


# ====================== 辅助函数 ======================

def motion_H(M, N, L, theta):
    """
    计算线性运动模糊滤波器的频率响应。
    参数:
        M, N (int): 图像的尺寸。
        L (float): 运动模糊的长度（以像素为单位）。
        theta (float): 运动的角度（以度为单位）。
    返回:
        H_c (ndarray): 运动模糊的中心化二维频率响应。
    """
    a = L * np.cos(np.deg2rad(theta))
    b = L * np.sin(np.deg2rad(theta))
    # 创建一个中心化的频率网格。
    U, V = freqgrid_centered(M, N)
    P = a * U + b * V
    # 使用sinc函数作为运动模糊模型。
    S = np.ones_like(P)
    nz = P != 0 # 避免在原点处除以零。
    S[nz] = np.sin(np.pi * P[nz]) / (np.pi * P[nz])
    # 包括相移。
    H_c = S * np.exp(-1j * np.pi * P)
    return H_c

def freqgrid_centered(M, N):
    """
    生成一个二维中心化频率网格。
    网格范围从-0.5到0.5。
    参数:
        M, N (int): 网格的尺寸。
    返回:
        U, V (ndarray): 代表频率坐标的二维数组。
    """
    u = np.arange(-M//2, M//2) / M
    v = np.arange(-N//2, N//2) / N
    V, U = np.meshgrid(v, u)
    return U, V

def inverse_filter(G, H, eps):
    """
    应用正则化的逆滤波器。
    F_hat = G / H，但避免了H中的小值导致的除法问题。
    参数:
        G (ndarray): 退化图像的FFT。
        H (ndarray): 退化（模糊）的传递函数。
        eps (float): 防止除以零的阈值。
    返回:
        Fhat (ndarray): 估计的恢复图像的FFT。
    """
    Fhat = np.zeros_like(G)
    # 仅在H的幅度大于阈值的地方进行除法。
    mask = np.abs(H) > eps
    Fhat[mask] = G[mask] / H[mask]
    return Fhat

def wiener_filter(G, H, K):
    """
    应用简化的维纳滤波器。
    F_hat = (H* / (|H|^2 + K)) * G
    参数:
        G (ndarray): 退化图像的FFT。
        H (ndarray): 退化（模糊）的传递函数。
        K (float): 估计的噪声信号功率比。
    返回:
        Fhat (ndarray): 估计的恢复图像的FFT。
    """
    # 频域中维纳滤波器的公式。
    denom = np.abs(H)**2 + K
    W = np.conj(H) / denom
    Fhat = W * G
    return Fhat

# ---------- 指标 ----------
def mse(I1, I2):
    """计算两幅图像之间的均方误差。"""
    return np.mean((I1.astype(np.float64) - I2.astype(np.float64)) ** 2)

def psnr(I1, I2):
    """计算两幅图像之间的峰值信噪比。"""
    m = mse(I1, I2)
    if m == 0:
        return float('inf')
    # 假设像素值归一化到[0, 1]范围。
    return 10 * np.log10(1.0 / m)

def ssim_manual(X, Y):
    """
    计算两幅图像之间的结构相似性指数（SSIM）。
    这个简化版本使用全局统计数据。
    """
    # 用于稳定分母较弱的除法的常数。
    C1, C2 = 0.01**2, 0.03**2

    # 图像的均值、方差和协方差。
    mu_x = np.mean(X)
    mu_y = np.mean(Y)
    sigma_x = np.var(X)
    sigma_y = np.var(Y)
    sigma_xy = np.mean((X - mu_x) * (Y - mu_y))

    # SSIM 公式。
    numerator = (2 * mu_x * mu_y + C1) * (2 * sigma_xy + C2)
    denominator = (mu_x**2 + mu_y**2 + C1) * (sigma_x + sigma_y + C2)
    return numerator / denominator

# ====================== 运行 ======================
# 该脚本设计为在Jupyter环境中运行，其中'main()'被显式调用。
if __name__ == "__main__":
    main()


================================================================================
 所有四种情况的恢复性能指标 
================================================================================

情况                          方法                MSE         PSNR(dB)    SSIM      
--------------------------------------------------------------------------------
低噪声, 低模糊                    模糊+噪声             0.00529     22.767      0.9680    
                            逆滤波器              0.01882     17.255      0.8945    
                            维纳滤波器             0.00159     27.999      0.9906    
--------------------------------------------------------------------------------
低噪声, 高模糊                    模糊+噪声             0.01702     17.691      0.8922    
                            逆滤波器              0.07051     11.518      0.6543    
                            维纳滤波器             0.00691     21.606      0.9597    
--------------------------------------------------------------------------------
高噪声, 低模糊                    模糊+噪声             0.00568     22.456      0.9657    
                            逆滤波器              0.20589     6.864       0.2393    
                            维纳滤波器             0.00797     20.988      0.9535    
--------------------------------------------------------------------------------
高噪声, 高模糊                    模糊+噪声             0.01739     17.596      0.8901    
                            逆滤波器              0.27327     5.634       0.1019    
                            维纳滤波器             0.01092     19.616      0.9356    
--------------------------------------------------------------------------------

结论¶

噪声模型催生了简单的滤波器（均值/中值）和变分去噪（TV）。对于模糊，反卷积需要一个PSF；维纳滤波是超越朴素反演的稳健基线。

参考文献¶

A. Buades, B. Coll, J.-M. Morel, “A non-local algorithm for image denoising”, CVPR, 2005.
A. Chambolle, “An Algorithm for Total Variation Minimization and Applications”, Journal of Mathematical Imaging and Vision, vol. 20, p. 89--97, 2004.
L.B. Lucy, “An iterative technique for the rectification of observed distributions”, Astronomical Journal, vol. 79, no 6, p. 745--754, 1974.
W.H. Richardson, “Bayesian-based iterative method of image restoration”, Journal of the Optical Society of America, vol. 62, no 1, p. 55--59, 1972.
L.I. Rudin, S. Osher, E. Fatemi, “Nonlinear total variation based noise removal algorithms”, Physica D, vol. 60, no. 1--4, p. 259--268, 1992.