基于重加权相关全变分模型的医学图像去噪方法

Medical Image Denoising Method Based on Reweighted Correlation Total Variation Model

摘要: 图像去噪一直是一个受到大量研究者关注的问题，并成功应用到医学等领域。典型的图像去噪方法是利用图像中存在的先验信息，例如低秩先验，局部光滑先验等。但是，目前的图像去噪方法并没有充分利用到图像的这些先验信息，导致去噪效果不是很理想。针对上述问题，本文提出了基于重加权相关全变分正则项的图像去噪模型。该模型利用重加权核范数的方式对相关全变分正则项进行约束，来保证更加充分地利用图像中的低秩先验和局部光滑先验，以此来提升图像恢复效果。我们应用该方法到医学图像中，并和常见的几种图像去噪方法进行比较，实验结果显示，该方法所得到的图像质量得到了明显的提升。

Abstract: Image denoising has always been a problem that attracts the attention of a large number of researchers, and has been successfully applied to medical fields. Typical image denoising methods use the prior information existing in the image, such as low-rank prior, local smooth prior, etc. However, the current image denoising methods do not make full use of the prior information of the image, resulting that the denoising effect is not very ideal. To solve the above problems, this paper proposed an image denoising model based on reweighted correlation total variation regularization term. In this model, the reweighted nuclear norm method is used to constrain the correlated total variation regularization term to ensure that the low-rank prior and local smooth prior in the image are fully utilized, so as to improve the image restoration effect. We applied this method to medical images and compared it with several common image denoising methods. The experimental results show that the image quality obtained by this method has been significantly improved.

1. 引言

图像作为基本的视觉信息，占据着人类获得信息的主要部分。然而，在采集、传输和存储图像的过程中，图像不可避免地会受到噪声污染。噪声的存在不仅降低了图像的质量，而且还干扰了信息传输和信息获取的准确性，这影响了后续图像处理目标检测 [1] 、数据分类 [2] 和图像分割 [3] 等下游任务。尤其在分析和阅读医学图像时，图像对比度、边缘特征和信噪比对医学诊断的正确性至关重要。因此，医学图像去噪是图像去噪的重要任务之一。

一般来讲，基于模型的图像去噪方法通常是在数学模型中设计不同的先验项来反映图像的先验信息，以此来恢复图像。所以，设计准确的先验项可以极大地利用数据丰富的先验信息，从退化的图像中恢复包含更多特征信息的高质量图像。在现实生活中，大多数图像都具有低秩先验，基于此先验信息，一个代表性的图像去噪方法 [4] 是鲁棒主成分分析(RPCA)。它的优化模型如下：

$\begin{array}{l}\underset{X,S}{\min }{\Vert X\Vert }_{\ast }+\lambda {\Vert S\Vert }_{{}_1}\\ \text{s}\text{.t}.M=X+S,\end{array}$ (1)

其中 $M\in R^{n_1\times n_2}$ 是观测矩阵，X和S分别对应低秩和稀疏矩阵； ${\Vert ·\Vert }_{*}={\displaystyle {\sum }_{i=1}^{\min \left(n_1,n_2\right)}{\sigma }_i}$ 表示矩阵的核范数， ${\sigma }_i$ 代表矩阵的第i个最大的奇异值； ${\Vert \Vert }_1$ 表示矩阵的 $l_1$ 范数。但除了低秩先验之外，局部光滑先验也是图像的一个代表性先验。近年来，许多有关考虑低秩先验和局部光滑先验的图像去噪方法被提出。这类去噪方法通常通过线性组合对全变分(TV)正则项和低秩正则项进行建模。例如，Wang [5] 等人将TV正则项加入低秩矩阵恢复模型，可以有效地将干净的低秩矩阵与高密度稀疏噪声分开，同时保持图像结构的平滑性。然而，两个正则项会给模型增加一个平衡参数，导致比较复杂的算法，并且可能导致在图像去噪过程中一些信息的损失。为了改善这一方式，最近Peng [6] 等人提出了联合考虑图像低秩先验和局部光滑先验的相关全变分正则项(CTV) ${\Vert X\Vert }_{\text{CTV}}={\displaystyle {\sum }_{k=1}^2{\Vert {\nabla }_{k}X\Vert }_{*}}$ ，其中 ${\nabla }_k$ 表示对矩阵的第k维度的微分。该正则项将低秩先验和局部光滑先验合并为一个正则项，能够较明显提升图像的恢复效果，但是核范数最小化(NNM)有一个明显的缺点，即同时最小化所有的奇异值，这可能没有充分利用图像的先验信息，因此，图像的恢复效果还有待提升。

基于上述分析，本文我们提出一个基于重加权相关全变分模型的去噪模型(WCTV)，尽可能地充分利用图像的先验信息。具体来说，我们对相关全变分正则项进行重加权操作，这样既可以在联合考虑图像低秩先验和局部光滑先验的同时，又能尽可能地减少图像信息的流失。我们采用著名的交替方向乘子法(ADMM)来解决这个图像去噪模型。然后，我们将提出的去噪模型应用于不同噪声强度条件下的医学图像，并与不同类型的去噪方法进行比较实验。实验结果表明，我们的方法在几乎所有的噪声情况下都能得到比其他比较方法更好的图像恢复效果。

2. WCTV模型的建立

2.1. 相关全变分图像去噪模型

Peng [6] 等人提出的相关全变分图像去噪模型，将低秩先验和局部光滑先验联合考虑在一个正则项中。该模型与其他图像去噪模型相比，在既利用更多图像先验信息的同时，又避免了算法的复杂性，很好地提升了目标图像的恢复效果。该模型的表达式如下所示：

$\begin{array}{l}\underset{X,S}{\min }{\displaystyle \underset{k=1}{\overset{2}{\sum }}{\Vert {\nabla }_{k}X\Vert }_{*}}+2\lambda {\Vert S\Vert }_1\\ \text{s}\text{.t}.M=X+S\end{array}$ (2)

我们假设 $G_k={\nabla }_k\left(X\right)$ ，则该模型的等价模型为：

$\begin{array}{l}\underset{X,S}{\min }{\displaystyle \underset{k=1}{\overset{2}{\sum }}{\Vert G_k\Vert }_{\ast }+2\lambda {\Vert S\Vert }_1}\\ \text{s}\text{.t}.M=X+S\text{,}{G}_k={\nabla }_k\left(X\right),k=1,2\end{array}$ (3)

2.2. 加权核范数方式的选择

为了在图像去噪工作中提高NNM的灵活性，Gu [7] [8] 等人提出了加权核范数的去噪模型(WNNM)，其中核范数的加权策略为： ${\Vert X\Vert }_{\omega ,\ast }={\displaystyle {\sum }_{i=1}^{\min \left(n_1,n_2\right)}{\omega }_i{\sigma }_i}\left(X\right)$ ，其中 ${\omega }_i=C/\left({\sigma }_i\left(X\right)+\epsilon \right)$ ，C是一个折衷常数且大于0， $\epsilon >0$ 是一个小常数，目的是为防止分母为0。与WNNM相比，Xie [9] 等人提出了一个更加灵活的核范数加权方式，即加权Schatten-p范数(WSNM)， ${\Vert X\Vert }_{\omega ,S_P}={\left({\displaystyle {\sum }_{i=1}^{\min \left(m,n\right)}{\omega }_i{\sigma }_i^p\left(X\right)}\right)}^{1/p}$ ，其中 $0<p\le 1$ 。WSNM不仅对原始矩阵的秩给出了更好的近似，而且还考虑不同奇异值的重要性。然而该方式的求解问题是非凸的，很难直接进行求解。

最近，Huang [10] 等人提出了一个广义的迭代重加权核范数的加权方式，即：

${\omega }_i=C/{\left({\sigma }_i\left(X\right)+\epsilon \right)}^p$ (4)

与WSNM相比，这种加权方式在求解时是一个凸问题，它保留了真实自然图像中常见的关于低秩矩阵奇异值重尾分布的更好先验信息，这使得其在图像去噪中对噪声和异常值具有很好的鲁棒性 [11] 。因此，对比以上几种加权方式，本文选择将重加权核范数的加权方式应用到相关全变分去噪模型当中，以此更加充分地利用图像的低秩先验和局部光滑先验。

2.3. WCTV图像去噪模型

经过上述两节的分析，我们提出了重加权相关全变分图像去噪模型：

$\begin{array}{l}\underset{X,S}{\min }{\displaystyle \underset{k=1}{\overset{2}{\sum }}{\Vert G_k\Vert }_{\omega ,\ast }+2\lambda {\Vert S\Vert }_1}\\ \text{s}\text{.t}.M=X+S,G_k={\nabla }_k\left(X\right),k=1,2\end{array}$ (5)

其中的核范数加权策略我们采用(4)中的重加权方式。

3. WCTV模型的求解

本章我们使用著名的交替方向乘子法 [12] [13] 对模型(5)进行求解。首先，我们写出它的增广拉格朗日函数如下所示：

$\begin{array}{l}L\left(X,S,{\left\{{\Gamma }_k\right\}}_{k=1}^3,{\left\{G_k\right\}}_{k=1}^2\right)\\ ={\displaystyle \underset{k=1}{\overset{2}{\sum }}\left({\Vert G_k\Vert }_{\omega ,\ast }+\frac{\mu }{2}{\Vert {\nabla }_k\left(X\right)-G_k+\frac{{\Gamma }_k}{\mu }\Vert }_F^2\right)}+2\lambda {\Vert S\Vert }_{{}_1}+\frac{\mu }{2}{\Vert M-X-S+\frac{{\Gamma }_3}{\mu }\Vert }_F^2\end{array}$ (6)

其中 $\mu$ 是惩罚参数； ${\Gamma }_k,k=1,2,3$ 是拉格朗日乘子。在ADMM中，需要优化一个变量的同时固定其他变量。接下来，我们将讨论如何解决每个相关变量的子问题。

求解 $G_k$ 的子问题：固定其他变量，式(6)可以变为：

$\underset{G_k}{\min }{\Vert G_k\Vert }_{\omega ,\ast }+\frac{\mu }{2}{\Vert G_k-\left({\nabla }_k\left(X^t\right)+\frac{{\Gamma }_k^t}{\mu }\right)\Vert }_F^2$ (7)

我们定义 $Y_L={\nabla }_k\left(X^t\right)+\frac{{\Gamma }_k^t}{\mu }$ ，那么 $Y_L$ 的奇异值分解为 $Y_L=U\sum V^{\text{T}}$ ，其中 $\sum =diag\left\{{\sigma }_i\left(Y_L\right)\right\}$ 。根据文献 [10] 中的引理1，我们可以得到该子问题的解：

$G_k^{t+1}=U{\mathcal{S}}_{{\omega }_i/\mu }\left(\Sigma \right)V^{\text{T}}$ (8)

求解X的子问题：固定其他变量，式(6)可以变为：

$\underset{X}{\min }{\displaystyle \underset{k=1}{\overset{2}{\sum }}\frac{\mu }{2}{\Vert {\nabla }_k\left(X\right)-G_k^{t+1}+\frac{{\Gamma }_k^t}{\mu }\Vert }_F^2}+\frac{\mu }{2}{\Vert M-X-S^t+\frac{{\Gamma }_3^t}{\mu }\Vert }_F^2$ (9)

优化上述问题可视为求解以下线性问题：

$\left(\mu I+\mu {\displaystyle \underset{k=1}{\overset{2}{\sum }}{\nabla }_k^{\text{T}}{\nabla }_k}\right)\left(X\right)=\mu \left(M-S^t\right)+{\Gamma }_3^t+\mu {\displaystyle \underset{k=1}{\overset{2}{\sum }}{\nabla }_k^{\text{T}}}\left(G_k^{t+1}\right)-{\nabla }_k^T\left({\Gamma }_k^t\right)$ (10)

其中 ${\nabla }_k^{\text{T}}(·)$ 表示 ${\nabla }_k(·)$ 的转置算子； ${\nabla }_k^{\text{T}}{\nabla }_k$ 对应矩阵的块循坏结构，它可以被快速傅里叶变换(FFT)矩阵对角化。事实上，与文献 [14] 类似， $X^{t+1}$ 的闭式解可以通过对式(10)的两边进行傅里叶变换来快速确定：

$\{\begin{array}{l}H={\displaystyle \underset{k=1}{\overset{2}{\sum }}\mathcal{F}{\left(D_k\right)}^{*}\odot \mathcal{F}\left(\mu G_k^{t+1}-{\Gamma }_k^t\right),}\\ T_x={\left|\mathcal{F}\left(D_1\right)\right|}^2+{\left|\mathcal{F}\left(D_2\right)\right|}^2,\\ X^{t+1}={\mathcal{F}}^{-1}\left(\frac{\mathcal{F}\left(\mu M-\mu S^t+{\Gamma }_3^t\right)+H}{\mu 1+\mu T_x}\right)\end{array}$ (11)

其中1表示所有元素为1的张量， $\odot$ 表示对应元素的乘积， $\mathcal{F}(·)$ 表示傅里叶变换， ${|·|}^2$ 表示元素的平方运算。

求解S的子问题：固定其他变量，式(6)可以变为：

$\underset{S}{\min }2\lambda {\Vert S\Vert }_{{}_1}+\frac{\mu }{2}{\Vert S-\left(M-X^{t+1}+\frac{{\Gamma }_3^t}{\mu }\right)\Vert }_F^2$ (12)

该问题可以通过文献 [15] 中的软阈值操作来解决：

$S^{t+1}={\mathcal{S}}_{2\lambda /\mu }\left(M-X^{t+1}+\frac{{\Gamma }_3^t}{\mu }\right)$ (13)

其中 $\mathcal{S}$ 表示软阈值算子。

求解 ${\Gamma }_k$ 的子问题：根据ADMM的原理，乘子的更新如下公式：

$\{\begin{array}{l}{\Gamma }_k^{t+1}={\Gamma }_k^t+\mu \left({\nabla }_k\left(X^{t+1}\right)-G_k^{t+1}\right),k=1,2,\\ {\Gamma }_3^{t+1}={\Gamma }_3^t+\mu \left(M-X^{t+1}-S^{t+1}\right),\\ \mu =\mu \rho \end{array}$ (14)

其中 $\rho$ 是一个大于1的常数。我们将上述迭代过程总结为如下算法1。

4. 实验结果及分析

在本章中，我们应用我们所提出的方法到医学图像中，并于其他五种出色的图像去噪方法进行比较，以评估我们方法的恢复性能。这五种对比方法分别为：RPCA [4] ，WNNM [7] ，TNNR [16] ，WSTV [17] 和CTV [6] 。这些对比方法中的所有相关参数都通过给定的默认设置或遵循相应论文中的规则进行微调，以确保它们的最佳性能。实验数据我们选取了4个大小为256 × 256的医学图像，原始图像如图1所示。其中三个显示人的大脑，另外一个显示人的脚部。对于实验评估指标，我们采用图像处理领域常用的两个指标信噪比(PSNR)和结构相似性(SSIM)指标 [18] 来评估被各模型恢复图像的质量。

Figure 1. The 4 test images, all are numbered in order from 1 to 4, from left to right

图1. 4张测试图，图片从左到右依次编号为1至4

4.1. 实验设置

在WCTV模型中，我们需要提前设定的参数是平衡参数 $\lambda$ 和加权指数p。对于 $\lambda$ ，我们仍然延续CTV模型中的原始设定，即 $\lambda =2/\sqrt{n_1}$ 。至于加权指数p的设定，我们使用网格搜索的方法来确定其的最佳值。在图2中我们以数据图像1分别在稀疏噪声强度为0.1，0.15和0.2时作为测试实验，列出了不同p值对PSNR的影响。从图中可以看出，在不同噪声强度下，p值在0.1~0.2之间时，恢复性能较为稳定。因此我们决定在下面的实验中，统一将p取值为0.2。

Figure 2. Effect of different values of p on the PSNR of the test image

图2. 不同p的取值对测试图像PSNR的影响

4.2. 医学图像去噪

医学图像在现代疾病诊断中发挥着相当重要的作用，但难免在传输和存储等过程中受到噪声的影响，从而影响医师对诊断结果的判断 [19] 。本节为了验证我们重加权相关全变分模型的有效性与正确性，我们将我们的算法应用到四张测试医学图像数据中，并于五种对比算法进行比较。其中噪声强度我们施加了三种不同强度的稀疏噪声，分别为0.1，0.15和0.2。表1和表2分别显示了本文提出的去噪方法与另外5种方法在4张测试图像上实验的PSNR和SSIM指标的实验结果，为避免数据偶然性，表中的数值均为10次实验后的均值。具体来说，我们的方法无论在不同噪声强度和不同数据图像下均可以达到所有对比算法中最好的数值结果。与CTV方法进行比较，我们方法的PSNR值在每个数据图像中几乎都可以高出CTV方法的1 dB以上，同时SSIM值也有0.01以上的提高。因此，从评估数据指标来说，我们的方法具有提升恢复图像质量的性能。

Table 1. Comparison of PSNR values of each algorithm after noise removal under different noise intensities

表1. 不同噪声强度下各个算法去除噪声后的PSNR值比较

Table 2. Comparison of SSIM values of each algorithm after noise removal under different noise intensities

表2. 不同噪声强度下各个算法去除噪声后的SSIM值比较

为了将实验效果更加直观化，我们选取数据3在稀疏噪声强度为0.2和数据4在稀疏噪声强度为0.1时的实验效果图进行展示，如图3和图4。从图示的恢复效果来看，WCTV和CTV的去噪效果明显优于其他对比算法。而我们的方法与CTV相比，WCTV在恢复医学图像的局部纹理和细节方面更有优势，这种细节和纹理的变化对医生诊断病人的病情有很大的帮助。综上所述，WCTV方法比其他方法具有更好的性能。

(a) Truth (b) Noisy (c) RPCA (d) WNNM(e) TNNR (f) WSTV (g) CTV (h) WCTV

Figure 3. Results of various algorithmic noise restoration for medical image data 3 with sparse noise intensity of 0.2. (a) Truth; (b) Noisy: PSNR = 10.80 dB, SSIM = 0.10; (c) RPCA:PSNR = 23.01 dB, SSIM = 0.83; (d) WNNM:PSNR = 24.01 dB, SSIM = 0.77; (e) TNNR:PSNR = 25.04 dB, SSIM = 0.82; (f) WSTV:PSNR = 27.61 dB, SSIM = 0.87; (g) CTV:PSNR = 29.25 dB, SSIM = 0.93; (h) WCTV:PSNR = 30.50 dB, SSIM = 0.95

图3. 医学图像数据1在稀疏噪声强度为0.2的条件下各种算法噪声恢复结果。(a) 原图；(b) 噪声图：PSNR = 10.80 dB，SSIM = 0.10；(c) RPCA：PSNR = 23.01 dB，SSIM = 0.83；(d) WNNM：PSNR = 24.01 dB，SSIM = 0.77；(e) TNNR：PSNR = 25.04 dB，SSIM = 0.82；(f) WSTV：PSNR = 27.61 dB，SSIM = 0.87；(g) CTV：PSNR = 29.25 dB，SSIM = 0.93；(h) WCTV：PSNR = 30.50 dB，SSIM = 0.95

(a) Truth (b) Noisy (c) RPCA (d) WNNM(e) TNNR (f) WSTV (g) CTV (h) WCTV

Figure 4. Results of various algorithmic noise restoration for medical image data 4 with sparse noise intensity of 0.1. (a) Truth; (b) Noisy: PSNR = 13.64 dB, SSIM = 0.24; (c) RPCA:PSNR = 20.31 dB, SSIM = 0.84; (d) WNNM:PSNR = 23.92 dB, SSIM = 0.85; (e) TNNR:PSNR = 22.59 dB, SSIM = 0.85; (f) WSTV:PSNR = 24.72 dB, SSIM = 0.87; (g) CTV:PSNR = 25.13 dB, SSIM = 0.91; (h) WCTV:PSNR = 26.63 dB, SSIM = 0.92

图4. 医学图像数据4在稀疏噪声强度为0.1的条件下各种算法噪声恢复结果。(a) Truth；(b) Noisy：PSNR = 13.64 dB，SSIM = 0.24；(c) RPCA：PSNR = 20.31 dB，SSIM = 0.84；(d) WNNM：PSNR = 23.92 dB，SSIM = 0.85；(e) TNNR：PSNR = 22.59 dB，SSIM = 0.85；(f) WSTV：PSNR = 24.72 dB，SSIM = 0.87；(g) CTV：PSNR = 25.13 dB，SSIM = 0.91；(h) WCTV：PSNR = 26.63 dB，SSIM = 0.92

5. 总结与展望

在本文中，我们提出了基于重加权相关全变分的医学图像去噪模型。该去噪方法通过重加权的方式可以更好地利用图像的低秩先验和局部光滑先验，以此完成更高的图像恢复效果。通过医学图像去噪实验证明，我们的方法无论在数值结果还是视觉效果上都要优于其他对比去噪方法。在未来的工作中，我们将考虑如何自适应参数的问题以及在不同类噪声条件下如何保持更佳的图像恢复效果。