如何使用Python实现各向同性线性扩散平滑图像去噪

你提到的常数扩散系数去噪滤波器本质是求解二维热传导方程，形式为：$\frac{\partial u}{\partial t} = d \cdot \Delta u$，初始条件为 $u(x,y,0) = I(x,y)$，其中$d=1$为扩散系数，$I(x,y)$为输入含噪图像，$u(x,y,t)$为扩散时长$t$对应的输出图像。

依赖库 #

• numpy：负责数值计算与有限差分迭代
• opencv-python：负责图像读写与预处理
• matplotlib：负责结果可视化

实现逻辑 #

我们采用有限差分法求解方程：

1. 图像预处理：将输入图像转为灰度图，归一化到[0,1]区间，避免迭代过程中数值溢出
2. 迭代计算：使用5点拉普拉斯差分模板计算每个像素的二阶微分，按dt时间步长更新像素值（dt取0.25满足CFL稳定条件，避免数值发散）
3. 结果后处理：将迭代完成的数值映射回[0,255]的8位图像格式

完整可运行代码 #

import cv2
import numpy as np
import matplotlib.pyplot as plt

def linear_diffusion_denoise(noisy_img, d=1, t=5, dt=0.25):
    # 输入归一化到0-1区间
    u = noisy_img.astype(np.float32) / 255.0
    # 计算总迭代次数
    iter_num = int(t / dt)
    # 5点拉普拉斯差分卷积核
    laplacian_kernel = np.array([[0, 1, 0],
                                 [1, -4, 1],
                                 [0, 1, 0]], dtype=np.float32)
    for _ in range(iter_num):
        # 计算当前步拉普拉斯值
        laplacian = cv2.filter2D(u, -1, laplacian_kernel)
        # 更新扩散结果
        u += dt * d * laplacian
        # 截断数值避免越界
        u = np.clip(u, 0, 1)
    # 映射回0-255的8位图像
    denoised_img = (u * 255).astype(np.uint8)
    return denoised_img

# 调用示例
if __name__ == "__main__":
    # 读入含噪灰度图，替换为你的本地图像路径即可
    noisy_img = cv2.imread("test_noisy_img.png", 0)
    # 测试指定的三个扩散时长
    t_list = [5, 10, 30]
    denoise_results = []
    for t in t_list:
        res = linear_diffusion_denoise(noisy_img, d=1, t=t)
        denoise_results.append(res)
    
    # 结果可视化
    plt.figure(figsize=(16, 4))
    plt.subplot(1, 4, 1)
    plt.imshow(noisy_img, cmap="gray")
    plt.title("含噪原图")
    plt.axis("off")
    for idx, t in enumerate(t_list):
        plt.subplot(1, 4, idx+2)
        plt.imshow(denoise_results[idx], cmap="gray")
        plt.title(f"扩散时长t={t}")
        plt.axis("off")
    plt.show()

注意事项 #

• 该方法为各向同性扩散，扩散时长越大去噪效果越强，但图像边缘也会越模糊，若需要保留边缘可改用非线性PM扩散模型
• 如果输入是彩色图像，拆分R/G/B三个通道分别调用函数处理后再合并即可
• 时间步长dt不要超过0.25，否则会出现数值不稳定，输出异常噪点

内容的提问来源于stack exchange，提问作者pythn langu