图像去噪

0. 项目使用

本项目实现了TV-L1/2 ADMM去噪算法和PnP-ADMM去噪算法。

部分运行效果：

使用方法如下：

git clone https://github.com/BugNotFoundX/image_denoise.git
cd image_denoise
conda create -n denoise python=3.12
conda activate denoise
pip3 install torch torchvision torchaudio --index-url https://download.pytorch.org/whl/cu126
pip install -r requirements.txt
python main.py

文件目录

.
├── main.py  // 主程序入口（可选不同方法，详见代码）
├── PnP.py  // PnP-ADMM算法实现（BM3D和DnCNN两种去噪器）
├── TotalVariation.py  // TV-L1/2 ADMM去噪算法实现
├── utils.py  // 工具函数和
├── Method.py // Method基类
├── model.pth // DnCNN模型
├── model.py  // DnCNN模型类
├── README.md  // 项目说明文档
├── image.png  // 测试图片
├── .gitignore  // Git忽略文件
├── denoised_comparison.png // 去噪结果对比图（运行main.py后更新）
└── requirements.txt  // python依赖

Tip：

• 代码中的命名一般使用公式中的符号，便于理解。
• github的公式渲染可能不支持所有公式，建议使用本地Markdown编辑器查看，本项目中使用VSCode编写和查看（typora也可正常查看）。

下文介绍了本项目实现的算法及其推导过程，符号含义如下表：

符号	含义
$u$	去噪后的图像
$f$	输入图像
$\lambda$	正则化参数
$\nabla u$	图像的梯度
$\nabla^T$	图像的散度
$F$	傅里叶变换
$F^{-1}$	傅里叶逆变换
$L$	拉格朗日函数
$F(\nabla^2)$	拉普拉斯算子的傅里叶变换
$|x|_1$	L1范数（曼哈顿范数）
$|x|_2$	L2范数（欧几里得范数）

1. TV-L1去噪算法

TV-L1去噪算法是一种基于总变差（Total Variation, TV）和L1范数的图像去噪方法。其主要思想是通过最小化图像的总变差来保留图像的边缘信息，同时使用L1范数来减少噪声。

该算法的目标是通过交替方向乘子法（ADMM）来求解以下优化问题[1]：

$$ \min_{u} |u - f|_1 + \lambda |\nabla u|_1 $$

其中，$f$ 是输入图像，$u$ 是去噪后的图像，$\lambda$ 是正则化参数，$\nabla u$ 表示图像的梯度。

ADMM 形式：

$$ \min_{u,z,y} |z|_1 + \lambda|y|_1 \\ \text{s.t. } z = u - f, y = \nabla u $$

对应的拉格朗日增广函数为（通过配方吸收一次项）：

$$ L_{\rho}(u,z,y,a,b)=|z|{1}+\lambda|y|{1}+\frac{\beta}{2}|z-u+f+\frac{a}{\beta}|^{2}+\frac{\alpha}{2}|y-\nabla u+\frac{b}{\alpha}|^{2} $$

参数更新：

1、u 更新：

$$ \frac{\partial L_{\rho}(u^{k},z^{k},y^{k},a^{k},b^{k})}{\partial u}=-\beta(z^{k}-u^{k}+f+\frac{a^{k}}{\beta})-\alpha\nabla^{T}(y^{k}-\nabla u^{k}+\frac{b^{k}}{\alpha})=0 $$

$$ =(z^{k}-u^{k}+f+\frac{a^{k}}{\beta})+\frac{\alpha}{\beta}\nabla^{T}(y^{k}-\nabla u^{k}+\frac{b^{k}}{\alpha})=0 $$

$$ =(z^{k}+f+\frac{a^{k}}{\beta})+\frac{\alpha}{\beta}\nabla^{T}(y^{k}+\frac{b^{k}}{\alpha})=u^{k}+\frac{\alpha}{\beta}\nabla^{2}u^{k} $$

$\rightarrow$ Fourier Transform :

$$ F(z^{k}+f+\frac{a^{k}}{\beta})+\frac{\alpha}{\beta}F(\nabla^{T}(y^{k}+\frac{b^{k}}{\alpha}))=F(u^{k}+\frac{\alpha}{\beta}\nabla^{2}u^{k}) $$

$$ \therefore u^{k+1}=F^{-1}(\frac{F(z^{k}+f+\frac{a^{k}}{\beta})+\frac{\alpha}{\beta}F(\nabla^{T}(y^{k}+\frac{b^{k}}{\alpha}))}{1+F(\frac{\alpha}{\beta}\nabla^{2})}) $$

2、z 更新：

(1) if $z>0$

$$ \frac{\partial L_{\rho}(u^{k+1},z^{k},y^{k},a^{k},b^{k})}{\partial z}=1+\beta(z^{k}-u^{k+1}+f+\frac{a^{k}}{\beta})=0 $$

$$ \rightarrow z^{k+1}=-\frac{1}{\beta}+u^{k+1}-f-\frac{a^{k}}{\beta}=v_{1}-\frac{1}{\beta},v_{1}=u^{k+1}-f-\frac{a^{k}}{\beta} $$

(2) if $z<0$

$$ \frac{\partial L_{\rho}(u^{k+1},z^{k},y^{k},a^{k},b^{k})}{\partial z}=-1+\beta(z^{k}-u^{k+1}+f+\frac{a^{k}}{\beta})=0 $$

$$ \rightarrow z^{k+1}=\frac{1}{\beta}+u^{k+1}-f-\frac{a^{k}}{\beta}=v_{1}+\frac{1}{\beta},v_{1}=u^{k+1}-f-\frac{a^{k}}{\beta} $$

(3) if $z=0$

(1)+(2)+(3): 结果为软阈值函数

$$ \therefore z^{k+1}=max(|v_{1}|-\frac{1}{\beta},0)sign(v_{1}),v_{1}=u^{k+1}-f-\frac{a}{\beta} $$

3、y 更新：

(1) if $y>0:$

$$ \frac{\partial L_{\rho}(u^{k+1},z^{k+1},y^{k},a^{k},b^{k})}{\partial y}=\lambda+\alpha(y^{k}-\nabla u^{k+1}+\frac{b^{k}}{\alpha})=0 $$

$$ \rightarrow y^{k+1}=-\frac{\lambda}{\alpha}+\nabla u^{k+1}-\frac{b^{k}}{\alpha}=v_{2}-\frac{\lambda}{\alpha},v_{2}=\nabla u^{k+1}-\frac{b^{k}}{\alpha} $$

(2) if $y<0:$

$$ \frac{\partial L_{\rho}(u^{k+1},z^{k+1},y^{k},a^{k},b^{k})}{\partial y}=-\lambda+\alpha(y^{k}-\nabla u^{k+1}+\frac{b^{k}}{\alpha})=0 $$

$$ \rightarrow y^{k+1}=\frac{\lambda}{\alpha}+\nabla u^{k+1}-\frac{b^{k}}{\alpha}=v_{2}+\frac{\lambda}{\alpha}, v_{2}=\nabla u^{k+1}-\frac{b^{k}}{\alpha} $$

(3) if $y=0$

(1)+(2)+(3)：结果为软阈值函数

$$ \therefore y^{k+1}=max(|v_{2}|-\frac{\lambda}{\alpha},0)sign(v_{2}), v_{2}=\nabla u^{k+1}-\frac{b^{k}}{\alpha} $$

4、a 更新：

$$ a^{k+1}=a^{k}+\beta(z^{k+1}-u^{k+1}+f) $$

5、b 更新：

$$ b^{k+1}=b^{k}+\alpha(y^{k+1}-\nabla u^{k+1}) $$

6、 $\alpha$ 更新：

$$ \alpha^{k+1}=\alpha^{k}+\rho $$

7、 $\beta$ 更新：

$$ \beta^{k+1}=\beta^{k}+\rho $$

一般的 $\rho > 0$ ，本项目取 $\rho=0.07$ 。

2. TV-L2去噪算法

TV-L2去噪算法是一种基于总变差（Total Variation, TV）和L2范数的图像去噪方法。其主要思想是通过最小化图像的总变差来保留图像的边缘信息，同时使用L2范数来减少噪声。

该算法的目标是通过交替方向乘子法（ADMM）来求解以下优化问题：

$$ \min_{u} \frac{1}{2}|u - f|_2^2 + \lambda |\nabla u|_1 $$

ADMM 形式：

$$ \min_{u,z,y} \frac{1}{2}|z|_2^2 + \lambda|y|_1 \\ \text{s.t. } z = u - f, y = \nabla u $$

对应的拉格朗日增广函数为（通过配方吸收一次项）：

$$ L_{\rho}(u,z,y,a,b)=\frac{1}{2}|z|_2^2 + \lambda|y|_1 + \frac{\beta}{2}|z-u+f+\frac{a}{\beta}|^2 + \frac{\alpha}{2}|y-\nabla u+\frac{b}{\alpha}|^2 $$

参数更新： 1、u 更新：

$$ \frac{\partial L_{\rho}(u^{k},z^{k},y^{k},a^{k},b^{k})}{\partial u}=-\beta(z^{k}-u^{k}+f+\frac{a^{k}}{\beta})-\alpha\nabla^{T}(y^{k}-\nabla u^{k}+\frac{b^{k}}{\alpha})=0 $$

$$ =(z^{k}-u^{k}+f+\frac{a^{k}}{\beta})+\frac{\alpha}{\beta}\nabla^{T}(y^{k}-\nabla u^{k}+\frac{b^{k}}{\alpha})=0 $$

$$ =(z^{k}+f+\frac{a^{k}}{\beta})+\frac{\alpha}{\beta}\nabla^{T}(y^{k}+\frac{b^{k}}{\alpha})=u^{k}+\frac{\alpha}{\beta}\nabla^{2}u^{k} $$

$\rightarrow$ Fourier Transform :

$$ F(z^{k}+f+\frac{a^{k}}{\beta})+\frac{\alpha}{\beta}F(\nabla^{T}(y^{k}+\frac{b^{k}}{\alpha}))=F(u^{k}+\frac{\alpha}{\beta}\nabla^{2}u^{k}) $$

$$ \therefore u^{k+1}=F^{-1}(\frac{F(z^{k}+f+\frac{a^{k}}{\beta})+\frac{\alpha}{\beta}F(\nabla^{T}(y^{k}+\frac{b^{k}}{\alpha}))}{1+F(\frac{\alpha}{\beta}\nabla^{2})}) $$

2、z 更新：

$$ \frac{\partial L_{\rho}(u^{k+1},z^{k},y^{k},a^{k},b^{k})}{\partial z}=z^{k}+\beta(z^{k}-u^{k+1}+f+\frac{a^{k}}{\beta})=0 $$

$$ \rightarrow (1+\beta)z^{k} = \beta(u^{k+1}-f-\frac{a^{k}}{\beta}) $$

$$ \rightarrow z^{k+1} = \frac{\beta}{1+\beta}(u^{k+1}-f-\frac{a^{k}}{\beta}) $$

令 $v_1 = u^{k+1}-f-\frac{a^{k}}{\beta}$, $z^{k+1} = \frac{\beta}{1+\beta}v_1$.

3、y 更新：

(1) if $y>0:$

$$ \frac{\partial L_{\rho}(u^{k+1},z^{k+1},y^{k},a^{k},b^{k})}{\partial y}=\lambda+\alpha(y^{k}-\nabla u^{k+1}+\frac{b^{k}}{\alpha})=0 $$

$$ \rightarrow y^{k+1}=-\frac{\lambda}{\alpha}+\nabla u^{k+1}-\frac{b^{k}}{\alpha}=v_{2}-\frac{\lambda}{\alpha},v_{2}=\nabla u^{k+1}-\frac{b^{k}}{\alpha} $$

(2) if $y<0:$

$$ \frac{\partial L_{\rho}(u^{k+1},z^{k+1},y^{k},a^{k},b^{k})}{\partial y}=-\lambda+\alpha(y^{k}-\nabla u^{k+1}+\frac{b^{k}}{\alpha})=0 $$

$$ \rightarrow y^{k+1}=\frac{\lambda}{\alpha}+\nabla u^{k+1}-\frac{b^{k}}{\alpha}=v_{2}+\frac{\lambda}{\alpha}, v_{2}=\nabla u^{k+1}-\frac{b^{k}}{\alpha} $$

(3) if $y=0$

(1)+(2)+(3)：结果为软阈值函数

$$ \therefore y^{k+1}=max(|v_{2}|-\frac{\lambda}{\alpha},0)sign(v_{2}), v_{2}=\nabla u^{k+1}-\frac{b^{k}}{\alpha} $$

4、a 更新：

$$ a^{k+1}=a^{k}+\beta(z^{k+1}-u^{k+1}+f) $$

5、b 更新：

$$ b^{k+1}=b^{k}+\alpha(y^{k+1}-\nabla u^{k+1}) $$

6、 $\alpha$ 更新：

$$ \alpha^{k+1}=\alpha^{k}+\rho $$

7、 $\beta$ 更新：

$$ \beta^{k+1}=\beta^{k}+\rho $$

一般的 $\rho > 0$，本项目取 $\rho=0.07$ 。

3. PnP-ADMM算法

PnP-ADMM（Plug-and-Play Alternating Direction Method of Multipliers）算法是一种结合了ADMM和深度学习的图像去噪方法。该算法通过引入深度学习模型作为先验知识，来提高去噪效果。

在TV-L2的基础上，PnP-ADMM算法将深度学习模型作为一个黑盒函数$g$，用于更新去噪图像。其优化问题可以表示为[2]：

$$ \min_{u, v} \frac{1}{2}|u - f|_2^2 + \lambda g(v) \\ \text{s.t. } u - v = 0 $$

该式的Lagrangian函数为：

$$ L_{\rho}(u,v,y)=\frac{1}{2}|u - f|_2^2 + \lambda g(v) + y^T(u - v) + \frac{\rho}{2}|u - v|_2^2 $$

根据ADMM更新参数：

$$ u^{k+1} = \arg\min_u \frac{1}{2}|u - f|_2^2 + \frac{\rho}{2}|u - v^k + \frac{y^k}{\rho}|_2^2 $$

$$ v^{k+1} = \arg\min_v \lambda g(v) + \frac{\rho}{2}|v - u^{k+1} - \frac{y^{k}}{\rho}|_2^2 $$

$$ y^{k+1} = y^k + \rho(u^{k+1} - v^{k+1}) $$

$u$ 更新含有正向模型 $\frac{1}{2}|u - f|_2^2$ ，可看作反演步骤。由于 $v$ 更新涉及到先验知识 $g$，所以可以看作是去噪步骤。$v$的更新可以进一步写为：

$$ v^{k+1} = \arg\min_v g(v) + \frac{1}{2 \sigma^2}|v - \bar{v}^{k}|_2^2 $$

其中 $\bar{v}^{k} = u^{k+1} + \frac{y^{k}}{\rho}$ ，$\sigma^2 = \lambda / \rho$ 是噪声的方差。如果将 $\bar{v}^{k}$ 视为一个噪声图像，则后项就是噪声图像与原始图像的残差，再结合先验信息 $g(v)$ ，这其实就是一个去噪模型。上式可直接写为：

$$ v^{k+1} = D_\sigma(\bar{v}^{k}) $$

其中 $D_\sigma$ 是一个去噪模型。

此外，可以将 $\rho$ 采用变量形式加速收敛：

$$ \rho^{k+1} = \gamma^k \rho^k, \gamma^k > 1 $$

$u$ 更新策略：

$$ u^{k+1} = \arg\min_u \frac{1}{2}|u - f|_2^2 + \frac{\rho^k}{2}|u - v^k + \frac{y^k}{\rho^k}|_2^2 $$

$$ \frac{\partial L_{\rho}(u,v^{k},y^{k})}{\partial u} = u - f + \rho^k\left(u - v^k + \frac{y^k}{\rho^k}\right) = 0 $$

$$ \Rightarrow u^{k+1} = \frac{f + v^k\rho^k - y^k}{1 + \rho^k} $$

关于惩罚参数 $\rho_k$ 的更新策略中，$\gamma_k$ 的选择也有两种方案：

1. 对于所有的迭代 $k$，$\gamma_k$ 都取同一个常量，即 $\gamma_k = \gamma$，其中 $\gamma > 1$。
2. 自适应地选取 $\gamma_k$，选取规则可以是：

$$ \gamma_k = \begin{cases} \gamma, & \Delta_{k+1} \ge \eta \Delta_k \ 1, & \Delta_{k+1} < \eta \Delta_k \end{cases} $$

其中，$\eta \in (0,1)$，$\gamma > 1$，

$$ \Delta_{k+1} = (1/\sqrt{n}) (|u^{k+1} - u^{k}|_2 + |v^{k+1} - v^{k}|_2 + |y^{k+1} - y^{k}|_2) $$

$n$ 是当前的循环次数

综上所述，该 PnP-ADMM 的迭代模型（使用 scaled dual variable $u$ 的形式，对应图片中的 $u^{k}$ 为对偶变量，图片中的 $x^{k}$ 为主变量 $u$，图片中的 $v^{k}$ 为辅助变量 $v$可以总结为（工程上可将将变量 $\frac{y^{k}}{\rho_k}$ 视作整体，及令 $y^k=\frac{y^{k}}{\rho_k}$）：

$$ \begin{cases} u^{k+1} = (f + \rho^k(v^k - y^k))/(1 + \rho^k) \\ v^{k+1} = D_{\sigma_k}(u^{k+1} + y^{k}) \\ y^{k+1} = y^{k} + u^{k+1} - v^{k+1} \\ \rho_{k+1} = \gamma_k \rho_k \end{cases} $$

代码中PnP_BMD3D和PnP_CNN的去噪模型分别为BM3D[3]和DnCNN[4]

Tips：

• 关于 $\gamma$ 的选取：该参数影响CNN去噪强度，图片的初始噪声越大，初始值应该越大，本项目中当噪声的 $\sigma=0.07$ ，（ $\rho=0.07,\lambda=1$ ）， $\gamma_0$ 较好的取值大约为 $10$ 左右。
• 因为不明原因使用BM3D该方法的表现并不好，论文[5]中似乎效果也不理想
• 增大 $\gamma$ 的初始值也可以增大去噪强度（图片看起来也会越油腻，对TV方法也生效）

代码参考

• DnCNN：https://github.com/cszn/DnCNN
• TV_L1_Denoising_via_ADMM：https://github.com/Git-H2K/TV_L1_Denoising_via_ADMM
• pypher：https://github.com/aboucaud/pypher

参考文献

[1] QIN Z T, GOLDFARB D, MA S. An alternating direction method for total variation denoising[J/OL]. Optimization Methods Software, 2015, 30(3): 594-615. DOI:10.1080/10556788.2014.955100.

[2] CHAN S H, WANG X, ELGENDY O A. Plug-and-Play ADMM for Image Restoration: Fixed Point Convergence and Applications[A/OL]. arXiv, 2016[2025-06-13]. http://arxiv.org/abs/1605.01710. DOI:10.48550/arXiv.1605.01710.

[3] DABOV K, FOI A, KATKOVNIK V, 等. Image Denoising by Sparse 3-D Transform-Domain Collaborative Filtering[J/OL]. IEEE Transactions on Image Processing, 2007, 16(8): 2080-2095. DOI:10.1109/TIP.2007.901238.

[4] ZHANG K, ZUO W, CHEN Y, 等. Beyond a Gaussian Denoiser: Residual Learning of Deep CNN for Image Denoising[J/OL]. IEEE Transactions on Image Processing, 2017, 26(7): 3142-3155. DOI:10.1109/TIP.2017.2662206.

[5] 魏俊豪. 基于ADMM的图像去噪算法研究[D]. 黑龙江:哈尔滨工程大学,2020.