Richardson-Lucy (RL) 反卷积算法 —— 通过不断迭代更新图像估计值

一、RL反卷积算法

Richardson-Lucy（RL）反卷积算法：是一种基于最大似然估计（Maximum Likelihood Estimation, MLE）的图像去卷积方法，常用于去模糊和图像增强。通过不断地更新图像估计值，来最小化模糊图像和恢复图像之间的误差，逐步消除模糊，最终接近真实的原始图像，每次迭代都会根据当前估计的图像来修正下次迭代的结果。该算法最初由 Richardson 和 Lucy 在 1972 年提出，应用于天文学中的 X 射线成像，但随后被广泛应用于各类图像处理领域，尤其是在医学成像、显微镜图像分析、遥感图像恢复等领域取得了显著的应用成果。

（1）主要特点

• RL 算法在处理具有已知点扩散函数（Point Spread Function，PSF）的图像时表现较好，适用于泊松噪声或高斯噪声的图像。对于强噪声或未知 PSF 的图像，算法的效果可能受到影响。
  • （1）不仅能恢复模糊图像，还能处理一定类型的噪声 —— 最适合处理泊松噪声，对高斯噪声也有一定的适应性。
  • （2）点扩散函数（PSF）是描述成像系统模糊特性的函数，它决定了图像的模糊程度。RL 反卷积的性能与 PSF 的估计精度密切相关。RL 反卷积能够很好地恢复使用已知 PSF 的图像，并且可以处理某些复杂的模糊（如运动模糊和焦距模糊），但对未知或不准确估计的 PSF 较为敏感，恢复效果较差。
    • 已知的PSF（效果最优）
      • 适用情况：如果模糊过程中的 PSF 已知，RL 反卷积能够利用该信息进行精准的图像恢复。这是 RL 反卷积最常见的应用场景。
      • 效果：在这种情况下，RL 可以准确地去除图像中的模糊，恢复出接近原始图像的细节。
    • 未知的PSF（效果打折）
      • 适用情况：在实际应用中，PSF 可能并不完全已知，这时需要通过其他方法进行估计。常见的 PSF 估计方法包括自适应估计、图像梯度法等。
      • 效果：如果 PSF 的估计不准确，RL 反卷积的恢复效果会显著降低，因为不准确的 PSF 可能导致图像恢复过程中出现严重的伪影或不真实的细节。在这种情况下，RL 的效果会大打折扣，恢复出的图像质量可能不如预期。
    • 扩展的PSF（效果复杂）
      • 适用情况：在一些情况下，模糊过程可能是由多个不同的模糊源引起的。例如，图像可能同时经历了不同类型的运动模糊和焦距模糊，这些模糊可能由不同的 PSF 共同作用。
      • 效果：RL 算法能够处理一定程度的复合模糊（即多个 PSF 的组合），但需要足够的迭代次数和合适的初始估计，效果较为复杂且不一定完全精确。

（2）基本原理

• 基本思想：利用图像的模糊模型，通过迭代过程逐步逼近原始图像。
• 模糊过程：通常可以通过卷积（或乘积）来表示，即假设模糊图像 𝑔 是原始图像 𝑓 与点扩散函数（PSF）ℎ 的卷积操作∗：𝑔 = 𝑓 * ℎ + noise。
• 迭代更新：通过对 𝑔 进行迭代更新，逐步恢复 𝑓，即寻找一个最佳的 𝑓，使得通过卷积得到的结果最接近模糊图像 𝑔。

（3）关键步骤

• 初始化：首先需要选择一个初始图像，通常可以选择为全 1 图像或其他合适的初值。
• 迭代更新：根据最大似然估计，通过不断迭代更新图像估计值。在每次迭代中，根据当前图像估计与模糊图像的比值进行更新，直到达到停止条件（如达到最大迭代次数或收敛）。
• 停止条件：通常基于误差的收敛或达到预设的最大迭代次数来停止。

（4）优化算法

Richardson-Lucy 反卷积算法在图像去模糊和恢复过程中，可能会放大噪声，导致伪影和失真。为了提高恢复质量，可以采用以下优化策略：

• 噪声抑制：在 RL 迭代过程中，噪声可能会被放大，因此可以引入以下方法进行抑制：
  • Wiener 滤波：在 RL 迭代前或后使用 Wiener 滤波，以抑制高频噪声并平滑图像。
  • 约束最小二乘 (Constrained Least Squares, CLS)：在 RL 过程中加入约束条件，减少噪声放大，增强稳定性。
• 迭代终止准则：RL 反卷积的迭代次数过多可能导致过拟合噪声，可以采取以下策略动态终止：
  • 计算相邻两次迭代结果的差异，若变化幅度低于阈值，则停止迭代。
  • 根据噪声水平设定最大迭代次数，防止过度增强噪声。
• 正则化技术：在 RL 反卷积过程中，可加入正则化约束，以避免恢复过程中过度增强噪声或出现伪影：
  • Total Variation (TV) 正则化：通过抑制图像的总变差，减少高频噪声，同时保留边缘信息。TV 正则化可与 RL 结合，在每次迭代后进行 TV 去噪。
  • L2 正则化：通过约束恢复图像的幅值变化，减少震荡效应，增强鲁棒性。
• 多分辨率策略：对于复杂模糊情况，可以使用多分辨率（图像金字塔） 进行优化：
  • 低分辨率恢复 → 逐步细化至高分辨率：先在低分辨率图像上进行 RL 反卷积，然后逐步增加分辨率，细化细节。
  • 拉普拉斯金字塔分解：在不同尺度恢复低频信息，并逐层添加高频细节，避免直接高分辨率恢复带来的噪声放大问题。
• 结合深度学习
  • 使用深度神经网络（如 U-Net、GAN）辅助 RL 反卷积，提高去模糊效果。
  • 通过深度学习预测更合理的 PSF（点扩散函数），提升恢复质量。

二、项目实战

（1）RL 反卷积

import numpy as np
import cv2
import matplotlib.pyplot as plt
from scipy.ndimage import gaussian_filterdef generate_psf(size=5, sigma=3):"""更灵活的生成 PSF"""psf = np.zeros((size, size))psf[size//2, size//2] = 1  # 设置中心点psf = gaussian_filter(psf, sigma=sigma)  # 高斯模糊psf /= psf.sum()  # 归一化return psfdef richardson_lucy_deconvolution(image, psf, iterations=10):"""使用 Richardson-Lucy 算法进行图像去卷积（去模糊）。参数：:param image: 输入的模糊图像。:param psf: 点扩散函数（PSF）。:param iterations: 迭代次数。:return: 去模糊后的图像。"""# （1）初始化估计图像# estimate = np.copy(image).astype(np.float32)  # 使用输入图像，作为初始估计图像estimate = np.ones_like(image, dtype=np.float32)   # 使用全 1 的图像，作为初始估计图像# （2）获取 PSF 的共轭矩阵psf_conj = np.conj(psf)  # 共轭矩阵 = 矩阵的转置 + 虚数单位矩阵# （3）迭代过程for _ in range(iterations):convolved = cv2.filter2D(estimate, -1, psf)  # 计算卷积（估计图像与 PSF）relative_blur = image / (convolved + 1e-8)  # 计算比值（模糊图像与卷积结果） ———— 1e-8防止除零错误correction = cv2.filter2D(relative_blur, -1, psf_conj)  # 计算卷积（比值图像与 PSF 的共轭矩阵）estimate *= correction  # 更新估计图像return estimateif __name__ == '__main__':image = cv2.imread('YH220_235_gaussianBlur.png', cv2.IMREAD_GRAYSCALE)###################################################################################"""定义高斯PSF"""psf = cv2.getGaussianKernel(5, 3)  # 5x5 高斯核psf = psf @ psf.T  # 生成 2D 高斯核print("psf =", psf)"""##################################################################################函数说明：cv2.getGaussianKernel(ksize, sigma)参数说明：ksize（核大小）：生成的1D高斯核。如 5×1，通过 @ 矩阵乘法扩展成 5×5。核的尺寸一般为奇数（如 3、5、7...），以保证卷积中心对称。较小高斯核：仅影响邻近像素，去模糊能力较弱。较大高斯核：影响范围更大，适合更大尺度的模糊或去噪。sigma（标准差）：控制高斯分布的 扩散程度，值越大，高斯分布越平滑（即模糊范围更大）。较小sigma：高斯核集中，影响范围小（锐化效果更强）。较大sigma：高斯核扩展，影响范围大（模糊效果更强）。###################################################################################"""value = 1if value == 0:"""单次测试"""restored_image = richardson_lucy_deconvolution(image, psf, iterations=10)plt.figure(figsize=(12, 6))plt.subplot(1, 2, 1), plt.title("Blurred Image"), plt.imshow(image, cmap='gray'), plt.axis('off')plt.subplot(1, 2, 2), plt.title("Restored Image"), plt.imshow(restored_image, cmap='gray'), plt.axis('off')plt.show()elif value != 0:"""多次测试"""results = []for iters in range(12):results.append(richardson_lucy_deconvolution(image, psf, iterations=iters))# 可视化plt.figure(figsize=(12, 8))for i, (img, iters) in enumerate(zip(results, range(12))):plt.subplot(3, 4, i+1)plt.imshow(img, cmap='gray')plt.title(f'iterations: {iters}')plt.axis('off')plt.tight_layout()plt.show()

（2）优化：RL 反卷积 + Total Variation (TV) 正则化 —— 相比于RL反卷积，效果各有千秋

import numpy as np
import cv2
import matplotlib.pyplot as plt
from skimage.restoration import denoise_tv_chambolledef richardson_lucy_deconvolution(image, psf, iterations=10, tv_denoise=True):estimate = np.ones_like(image, dtype=np.float32)  # 初始化估计图像psf_conj = np.flip(psf)  # 取 PSF 的翻转（近似共轭）for i in range(iterations):convolved = cv2.filter2D(estimate, -1, psf)relative_blur = image / (convolved + 1e-8)correction = cv2.filter2D(relative_blur, -1, psf_conj)estimate *= correctionif tv_denoise:print(f"iterations: {i}, TV denoise: {tv_denoise}")estimate = denoise_tv_chambolle(estimate, weight=0.2)  # TV 去噪return estimateif __name__ == '__main__':# 读取输入图像（灰度图）image = cv2.imread('YH220_235_gaussianBlur.png', cv2.IMREAD_GRAYSCALE)if image is None:raise FileNotFoundError("图像文件未找到，请检查路径！")# 定义高斯 PSF（点扩散函数）psf = cv2.getGaussianKernel(5, 3)  # 5x1 高斯核psf = psf @ psf.T  # 生成 5x5 2D 高斯核# 进行 Richardson-Lucy 反卷积restored_image = richardson_lucy_deconvolution(image, psf, iterations=10)# 显示结果plt.figure(figsize=(12, 6))plt.subplot(1, 2, 1), plt.title("Blurred Image"), plt.imshow(image, cmap='gray'), plt.axis('off')plt.subplot(1, 2, 2), plt.title("Restored Image"), plt.imshow(restored_image, cmap='gray'), plt.axis('off')plt.show()

（3）优化：RL 反卷积 + 动态停止迭代 —— 用于设置较大迭代次数（如：100）

采用 PSNR、MSE 或 SSIM 作为衡量标准，动态决定何时停止迭代。
目前你的代码采用的是固定 10 次迭代，但在实际应用中，RL 迭代次数的选择很重要：

• 过少迭代：图像仍然模糊，恢复效果不理想。
• 过多迭代：可能放大噪声，导致伪影。

import numpy as np
import cv2
import matplotlib.pyplot as plt
from skimage.metrics import peak_signal_noise_ratio as psnrdef richardson_lucy_deconvolution(image, psf, iterations=50, tol=1e-3):estimate = np.ones_like(image, dtype=np.float32)psf_conj = np.flip(psf)prev_psnr = 0for i in range(iterations):convolved = cv2.filter2D(estimate, -1, psf)relative_blur = image / (convolved + 1e-8)correction = cv2.filter2D(relative_blur, -1, psf_conj)estimate *= correction# 计算 PSNR，若提升不明显则终止迭代current_psnr = psnr(image, estimate)print(f"Iteration {i + 1},"f"prev_psnr = {prev_psnr:.2f}, "f"current_psnr: {current_psnr:.2f}, "f"PSNR DIFF: {abs(current_psnr - prev_psnr):.2f}")if abs(current_psnr - prev_psnr) < tol:print(f"Early stopping at iteration {i + 1}")breakprev_psnr = current_psnrreturn estimateif __name__ == '__main__':# 读取输入图像（灰度图）image = cv2.imread('YH220_235_gaussianBlur.png', cv2.IMREAD_GRAYSCALE)if image is None:raise FileNotFoundError("图像文件未找到，请检查路径！")# 定义高斯 PSF（点扩散函数）psf = cv2.getGaussianKernel(5, 3)  # 5x1 高斯核psf = psf @ psf.T  # 生成 5x5 2D 高斯核# 进行 Richardson-Lucy 反卷积restored_image = richardson_lucy_deconvolution(image, psf, iterations=100, tol=0.5)# 显示结果plt.figure(figsize=(12, 6))plt.subplot(1, 2, 1), plt.title("Blurred Image"), plt.imshow(image, cmap='gray'), plt.axis('off')plt.subplot(1, 2, 2), plt.title("Restored Image"), plt.imshow(restored_image, cmap='gray'), plt.axis('off')plt.show()
"""
Iteration 1,prev_psnr = 0.00, current_psnr: 43.01, PSNR DIFF: 43.01
Iteration 2,prev_psnr = 43.01, current_psnr: 48.22, PSNR DIFF: 5.20
Iteration 3,prev_psnr = 48.22, current_psnr: 43.99, PSNR DIFF: 4.22
Iteration 4,prev_psnr = 43.99, current_psnr: 42.24, PSNR DIFF: 1.75
Iteration 5,prev_psnr = 42.24, current_psnr: 41.34, PSNR DIFF: 0.91
Iteration 6,prev_psnr = 41.34, current_psnr: 40.78, PSNR DIFF: 0.55
Iteration 7,prev_psnr = 40.78, current_psnr: 40.41, PSNR DIFF: 0.37
Early stopping at iteration 7
"""

（4）优化：RL 反卷积 + 多分辨率（分辨率层数越多，效果越好） —— 相比于输入图像略有提升，但不如RL反卷积

import numpy as np
import cv2
import matplotlib.pyplot as pltdef richardson_lucy_deconvolution(image, psf, iterations=10):"""使用 Richardson-Lucy 反卷积进行图像去模糊"""estimate = np.ones_like(image, dtype=np.float32)  # 初始估计图像psf_conj = np.flip(psf)  # PSF 共轭矩阵（对 PSF 进行翻转）for _ in range(iterations):convolved = cv2.filter2D(estimate, -1, psf)  # 计算卷积relative_blur = image / (convolved + 1e-8)  # 计算比值correction = cv2.filter2D(relative_blur, -1, psf_conj)  # 计算修正项estimate *= correction  # 更新估计图像return estimatedef build_gaussian_pyramid(image, levels):"""构建高斯金字塔（用于多尺度处理）"""pyramid = [image]for _ in range(levels - 1):image = cv2.pyrDown(image)  # 逐层降采样pyramid.append(image)return pyramiddef build_laplacian_pyramid(gaussian_pyramid):"""构建拉普拉斯金字塔（保存高频细节）"""laplacian_pyramid = []for i in range(len(gaussian_pyramid) - 1):upsampled = cv2.pyrUp(gaussian_pyramid[i + 1],dstsize=(gaussian_pyramid[i].shape[1], gaussian_pyramid[i].shape[0]))  # 调整大小high_freq = cv2.subtract(gaussian_pyramid[i], upsampled)  # 计算高频laplacian_pyramid.append(high_freq)laplacian_pyramid.append(gaussian_pyramid[-1])  # 顶层保持不变return laplacian_pyramiddef multi_scale_richardson_lucy(image, psf, iterations=10, levels=3):"""多尺度 Richardson-Lucy 反卷积"""# 1️⃣ 构建高斯金字塔gaussian_pyramid = [image.astype(np.float32)]for _ in range(levels - 1):image = cv2.pyrDown(image)gaussian_pyramid.append(image.astype(np.float32))# 2️⃣ 构建拉普拉斯金字塔（保存高频信息）laplacian_pyramid = build_laplacian_pyramid(gaussian_pyramid)# 3️⃣ 低分辨率层开始反卷积estimate = richardson_lucy_deconvolution(gaussian_pyramid[-1], psf, iterations)# 4️⃣ 从低分辨率逐层恢复for i in range(levels - 2, -1, -1):estimate = cv2.pyrUp(estimate, dstsize=(gaussian_pyramid[i].shape[1], gaussian_pyramid[i].shape[0]))  # 放大estimate = cv2.resize(estimate, (laplacian_pyramid[i].shape[1], laplacian_pyramid[i].shape[0]))  # 确保尺寸匹配# 🚀 **确保数据类型一致**estimate = estimate.astype(np.float32)laplacian_pyramid[i] = laplacian_pyramid[i].astype(np.float32)# 叠加高频信息estimate = cv2.add(estimate, laplacian_pyramid[i])return estimateif __name__ == '__main__':"""加载图像"""image = cv2.imread('YH220_235_gaussianBlur.png', cv2.IMREAD_GRAYSCALE)"""定义高斯PSF"""psf = cv2.getGaussianKernel(5, 3)  # 5×5 高斯核psf = psf @ psf.T  # 生成 2D 高斯核"""使用多尺度 RL 反卷积进行去模糊"""restored_image = multi_scale_richardson_lucy(image, psf, iterations=10, levels=50)"""可视化结果"""plt.figure(figsize=(12, 6))plt.subplot(1, 2, 1), plt.title("Blurred Image"), plt.imshow(image, cmap='gray'), plt.axis('off')plt.subplot(1, 2, 2), plt.title("Restored Image"), plt.imshow(restored_image, cmap='gray'), plt.axis('off')plt.show()