CN104200441B - 基于高阶奇异值分解的磁共振图像去噪方法 - Google Patents

Abstract

一种基于高阶奇异值分解的磁共振图像去噪方法，包括：（1）对原始图像进行方差稳定变换；（2）对变换后的图像进行第一次高阶奇异值分解去噪，具体是：（a）对每一个标的像素对应的参考块，逐个寻找该参考块的相似块组成高维数组；（b）对高维数组进行高阶奇异值分解变换得到系数和自适应基；（c）通过第一阈值操作；（d）进行高阶奇异值分解逆变换；（3）加权平均进行像素合并；（4）加权平均得到第一次加权平均后的图像；（5）对第一次加权平均后的图像进行第二次高阶奇异值分解去噪获得第二次去噪后的图像；（6）对第二次去噪后的图像进行方差稳定逆变换得到滤波图像。本发明的方法能够有效地抑制磁共振图像噪声，提高图像质量。

Description

基于高阶奇异值分解的磁共振图像去噪方法

技术领域

本发明属于医学图像处理的技术领域，具体地说属于一种基于高阶奇异值分解的磁共振图像去噪方法。

背景技术

磁共振成像是当前临床医学影像学的重要检查手段之一。然而，由于成像机制的影响，成像过程中不可避免的会引入噪声。图像中的噪声会大大降低图像的质量，使得图像边缘变得模糊，细微结构难以辨认，从而影响临床诊断，降低分析任务的可信度，例如图像配准，图像分割和一些相关组织的参数测量(如灌注图像和弛豫时间的相关参数)。因此，对于临床诊断和图像分析来说，降低噪声是十分必要的。

降低噪声可以通过两种方法来实现，一种是增加采集次数，平均多次采集的数据来提高图像的信噪比。这种方法是以增加采集时间来提高信噪比的。另一种方法是通过图像后处理滤波技术来提高信噪比。图像滤波技术不需要增加采集时间，因而在临床研究中得到广泛应用。

近期，有研究者采用高阶奇异值分解方法(Rajwade，A.，Rangarajan，A.，Banerjee，A..Image denoising using the higher order singular valuedecomposition.IEEE Trans Pattern Anal Mach 2013，Intell35，849-862)对加有高斯噪声的自然图像进行去噪。该方法是一种简单的，基于块匹配，机器学习的稀疏去噪方法。该方法是将相似块组成高维数组，利用相似块之间的稀疏性进行去噪。该方法针对不同的高维数组，通过训练学习得到不同的自适应基，能够更好地表达图像内容。不同于奇异值分解，高阶奇异值分解不需要将高维数组展开成矩阵来进行分解，不会破坏块内部的相关信息。然而，这种通过图像本身学习得到的基，易受图像中噪声的影响尤其是当噪声严重的时候。

因此，针对现有技术不足，提供一种对加有莱斯噪声的磁共振图像进行去噪处理的方法以克服现有技术不足甚为必要。

发明内容

本发明的目的在于针对现有技术不足，提供一种基于高阶奇异值分解的磁共振图像去噪方法，具有去噪性好，且能够提高去噪处理后的图像质量。

本发明的上述目的通过如下技术手段实现。

一种基于高阶奇异值分解的磁共振图像去噪方法，依次包括如下步骤：

(1)对加有莱斯噪声的原始磁共振图像进行方差稳定变换，将依赖信号分布的莱斯噪声变成独立于信号分布的噪声，得到变换后的噪声图像；

(2)将变换后的噪声图像中的每一个像素作为标的像素或者按照一定距离间隔取像素作为标的像素，取标的像素及其周围的像素作为参考块，执行以下操作进行第一次高阶奇异值分解去噪，具体是：

(a)通过k最近邻方法依次逐个寻找参考块的相似块，组成高维数组；

(b)对高维数组进行高阶奇异值分解变换，得到对应的系数和自适应基；

(c)通过第一阈值操作，将绝对值小于第一阈值的系数置零得到修正后的系数；

(d)根据自适应基与修正后的系数，进行高阶奇异值分解逆变换，将逆变换后的结果作为该高维数组中所有图像块去除噪声后的估计值；

(3)对步骤(2)得到的估计值加权平均进行像素合并，得到像素合并后的估计值，对每一个像素进行像素合并构成第一次去噪后的图像；

(4)对方差稳定变换后的噪声图像和第一次去噪后的图像进行加权平均，得到第一次加权平均后的图像；

(5)对第一次加权平均后的图像进行第二次高阶奇异值分解去噪，具体是：

(5.1)将第一次加权平均后的图像中的每一个像素作为二次标的像素或者按照一定距离取像素作为二次标的像素，取二次标的像素及其周围的像素作为二次参考块，执行以下操作：

(a)以第一次去噪后的图像进行相似块判断，通过k最近邻方法依次逐个寻找二次参考块的二次相似块，组成二次高维数组；

(b)对二次高维数组进行高阶奇异值分解变换，得到对应的二次系数和二次自适应基；

(c)通过二次阈值操作，将绝对值小于二次阈值的二次系数置零得到修正后的二次系数；

(d)根据二次自适应基与修正后的二次系数，进行高阶奇异值分解逆变换，将逆变换后的结果作为该高维数组中所有图像块去除噪声后的二次估计值；

(5.2)对步骤(5.1)得到的二次估计值加权平均进行像素合并，得到像素合并后的二次估计值，对每一个像素进行像素合并构成第二次去噪后的图像；

(6)对第二次去噪后的图像进行方差稳定逆变换，得到最终的滤波图像。

优选的，上述步骤(2)中(c)分步骤的第一阈值大小其中p为块大小，d为图像维度，K为相似块的个数，σ为方差稳定变换后图像的标准差。

优选的，上述步骤(4)具体是采用迭代正则化方法将方差稳定变换后的噪声图像和第一次去噪后的图像进行加权平均。

优选的，上述步骤(4)对方差稳定变换后的噪声图像和第一次去噪后的图像进行加权平均具体为：

其中表示第一次去噪后的图像，Y表示方差稳定变化后的噪声图像，ξ∈[0，1]是松弛参数，表示噪声图像Y对图像Y_w的贡献。

优选的，上述步骤(5)中(c)分步骤的二次阈值大小其中γ∈[0，1]为尺度参数，控制图像Y_w的标准差σ_w大小。

优选的，上述步骤(2)的(a)中，通过k最近邻方法依次逐个寻找参考块的相似块，相似块的块大小和相似块个数的取值根据噪声水平而定；

当噪声水平不大于2％时，相似块的块大小为3，相似块个数为30-70；

当噪声水平大于2％而小于等于7％时，相似块的块大小为4，相似块个数为70-80；

当噪声水平大于7％而小于等于15％时，相似块的块大小为5，相似块个数为100-130。

优选的，上述步骤(5)的(a)中，通过k最近邻方法依次逐个寻找二次参考块的二次相似块，二次相似块的块大小和二次相似块个数的取值根据噪声水平而定；

当噪声水平不大于2％时，二次相似块的块大小为3，二次相似块个数为30-70；

当噪声水平大于2％而小于等于7％时，二次相似块的块大小为4，二次相似块个数为70-80；

当噪声水平大于7％而小于等于15％时，二次相似块的块大小为5，二次相似块个数为100-130。

本发明基于高阶奇异值分解的磁共振图像去噪方法，该方法通过两步高阶奇异值分解去噪。第一步去噪后，对方差稳定变换后的噪声图像和第一次去噪后的图像进行加权平均，对加权平均后的图像进行第二步高阶奇异值分解去噪。与原始噪声图像相比，加权平均后的图像含有的噪声变小，第二步去噪可以减少噪声对高阶奇异值分解所得基的影响，使得高阶奇异值分解后的基能够更好地表达图像内容。故本发明方法在去噪的同时可以很好的保留图像细节，提高图像质量。

附图说明

利用附图对本发明作进一步的说明，但附图中的内容不构成对本发明的任何限制。

图1为本发明基于高阶奇异值分解的磁共振图像去噪方法的流程示意图；

图2为加入噪声水平为15％莱斯噪声的T1加权、T2加权和质子密度(PD)加权图像使用本发明方法去噪后的结果示意图；

图3为对噪声水平为3％的真实数据，采用本发明方法去噪后的结果示意图；

图4是对应图3的局部放大图；

图5为对噪声水平为4.5％的真实数据，采用本发明方法去噪后的结果示意图；

图6是对应图5的局部放大图。

具体实施方式

结合以下实施例对本发明作进一步描述。

实施例1。

如图1所示，一种基于高阶奇异值分解的磁共振图像去噪方法，依次包括如下步骤。

(1)对加有莱斯噪声的原始磁共振图像进行方差稳定变换，将依赖信号分布的莱斯噪声变成独立于信号分布的噪声，得到变换后的噪声图像。

(2)将变换后的噪声图像中的每一个像素作为标的像素或者按照一定距离间隔取像素作为标的像素，取标的像素及其周围的像素作为参考块，执行以下操作进行第一次高阶奇异值分解去噪。

具体是：

(a)通过k最近邻方法依次逐个寻找参考块的相似块，组成高维数组。块大小p和相似块个数K的取值根据噪声水平不同而不同，两个参数在不同噪声水平下的取值通常遵循如下规律：

当噪声水平不大于2％时，相似块的块大小为3，相似块个数为30-70；当噪声水平大于2％而小于等于7％时，相似块的块大小为4，相似块个数为70-80；当噪声水平大于7％而小于等于15％时，相似块的块大小为5，相似块个数为100-130。

(b)对高维数组进行高阶奇异值分解变换，得到对应的系数和自适应基。

(c)通过第一阈值操作，将绝对值小于第一阈值的系数置零得到修正后的系数。其中，第一阈值大小其中p为块大小，d为图像维度，K为相似块的个数，σ为方差稳定变换后图像的标准差。

(d)根据自适应基与修正后的系数，进行高阶奇异值分解逆变换，将逆变换后的结果作为该高维数组中所有图像块去除噪声后的估计值。

(3)像素合并，因为可能存在某些像素出现在不同图像块中而具有多个估计值，通过对这些估计值加权进行像素合并，具体是对步骤(2)得到的估计值加权平均进行像素合并，得到像素合并后的估计值，对每一个像素进行像素合并构成第一次去噪后的图像。

(4)对方差稳定变换后的噪声图像和第一次去噪后的图像进行加权平均，得到第一次加权平均后的图像。

步骤(4)是采用迭代正则化方法将方差稳定变换后的噪声图像和第一次去噪后的图像进行加权平均，具体为：

其中表示第一次去噪后的图像，Y表示方差稳定变化后的噪声图像，ξ∈[0，1]是松弛参数，表示噪声图像Y对图像Y_w的贡献。

(5)对第一次加权平均后的图像进行第二次高阶奇异值分解去噪，具体是：

(5.1)将第一次加权平均后的图像中的每一个像素作为二次标的像素或者按照一定距离取像素作为二次标的像素，取二次标的像素及其周围的像素作为二次参考块，执行以下操作：

(a)通过k最近邻方法依次逐个寻找二次参考块的二次相似块，组成二次高维数组。为了提高寻找相似块的准确度，这里是用第一次去噪后的图像进行相似块判断。二次相似块的块大小和二次相似块个数的取值根据噪声水平而定，两个参数在不同噪声水平下的取值通常遵循如下规律：

当噪声水平不大于2％时，二次相似块的块大小为3，二次相似块个数为30-70；

当噪声水平大于2％而小于等于7％时，二次相似块的块大小为4，二次相似块个数为70-80；

当噪声水平大于7％而小于等于15％时，二次相似块的块大小为5，二次相似块个数为100-130。

(b)对二次高维数组进行高阶奇异值分解变换，得到对应的二次系数和二次自适应基。

(c)通过二次阈值操作，将绝对值小于二次阈值的二次系数置零得到修正后的二次系数。二次阈值大小其中γ∈[0，1]为尺度参数，控制图像Y_w的标准差σ_w大小。

(d)根据二次自适应基与修正后的二次系数，进行高阶奇异值分解逆变换，将逆变换后的结果作为该高维数组中所有图像块去除噪声后的二次估计值。

(5.2)对步骤(5.1)得到的二次估计值加权平均进行像素合并，得到像素合并后的二次估计值，对每一个像素进行像素合并构成第二次去噪后的图像。

(6)对第二次去噪后的图像进行方差稳定逆变换，得到最终的滤波图像。

本发明的方法通过两步高阶奇异值分解去噪。第一步去噪后，对方差稳定变换后的噪声图像和第一次去噪后的图像进行加权平均，对加权平均后的图像进行第二步高阶奇异值分解去噪。与原始噪声图像相比，加权平均后的图像含有的噪声变小，第二步去噪可以减少噪声对高阶奇异值分解所得基的影响，使得高阶奇异值分解后的基能够更好地表达图像内容。故本发明方法在去噪的同时可以很好的保留图像细节，提高图像质量。

实施例2。

为了验证本发明方法的效果，以模拟数据为去噪对象，进行本发明基于高阶奇异值分解的磁共振去噪方法处理过程。

本发明基于高阶奇异值分解的磁共振去噪方法，具体实施包括如下步骤：

(1)输入加有莱斯噪声的磁共振图像，通过方差稳定变换，将依赖信号分布的莱斯噪声变成独立于信号分布的噪声。

(2)对变换后的图像中的每一个像素(为提高计算效率，本实施例中隔一定距离取一个像素作为目标像素，每个目标像素的距离为p-1，p为块大小)，取其周围的像素作为参考块，执行以下操作：

(a)通过k最近邻方法找到参考块的相似块，然后把它们聚集到一个高维数组，块大小p和相似块个数K的取值根据噪声水平不同而不同，两个参数在不同噪声水平下的取值按照如下规律选取：

当噪声水平为1％时，相似块的块大小为3，相似块个数为35；当噪声水平为3％、5％和7％时，相似块的块大小为4，相似块个数为75；当噪声水平为9％、11％、13％时15％时，相似块的块大小为5，相似块个数为120。

(b)对该高维数组进行高阶奇异值分解，得到对应的系数和基。

(c)通过第一阈值操作，将绝对值小于第一阈值的那些系数置零来减弱噪声，第一阈值大小其中p为块大小，d为图像维度，K为相似块个数，σ为方差稳定变换后图像的标准差，一般稳定到1。

(d)根据自适应基与修正后的系数，进行高阶奇异值分解逆变换，将逆变换后的结果作为该高维数组中所有图像块去除噪声后的估计值。

(3)像素合并，因为可能存在某些像素出现在不同图像块中而具有多个估计值，通过对这些估计值加权进行像素合并，得到像素合并后的估计值，对每一个像素进行像素合并构成第一次去噪后的图像。

(4)对方差稳定变换后的噪声图像和第一次去噪后的图像进行加权平均：

其中表示第一次去噪后的图像，Y表示方差稳定变化后的噪声图像。ξ∈[0，1]是一个松弛参数，表示噪声图像Y对图像Y_w的贡献。本实验中参数ξ取值为0.38。

(5)对加权平均后的图像Y_w进行第二次高阶奇异值分解去噪。具体是：

(5.1)将加权平均后的图像按照一定距离取像素作为二次标的像素，取二次标的像素及其周围的像素作为二次参考块，执行以下操作：

(a)通过k最近邻方法依次逐个寻找二次标的像素的二次相似块，组成二次高维数组。为了提高寻找相似块的准确度，这里是用第一次去噪后的图像进行相似块判断。二次相似块的块大小和二次相似块个数的取值根据噪声水平而定，两个参数在不同噪声水平下的取值通常遵循如下规律：当噪声水平为1％时，相似块的块大小为3，相似块个数为35；当噪声水平为3％、5％和7％时，相似块的块大小为4，相似块个数为75；当噪声水平为9％、11％、13％时15％时，相似块的块大小为5，相似块个数为120。

(b)对二次高维数组进行高阶奇异值分解变换，得到对应的二次系数和二次自适应基。

(c)通过二次阈值操作，将绝对值小于二次阈值的二次系数置零得到修正后的二次系数。二次阈值大小其中γ∈[0，1]为尺度参数，控制图像Y_w的标准差σ_w大小。本实施例中参数γ取值0.65。

(d)根据二次自适应基与修正后的二次系数，进行高阶奇异值分解逆变换，将逆变换后的结果作为该高维数组中所有图像块去除噪声后的二次估计值。

(5.2)对步骤(5.1)得到的二次估计值加权平均进行像素合并，得到像素合并后的二次估计值，对每一个像素进行像素合并构成第二次去噪后的图像。

(6)为了得到图像的无偏估计，对第二次去噪后的图像进行方差稳定逆变换，得到最终的滤波图像。

将无噪声模拟数据T1加入15％莱斯噪声，采用本发明的方法对该噪声水平为15％的噪声图像进行去噪处理。模拟数据T1的图像、与T1对应的噪声水平为15％的噪声图像及采用本发明的方法去噪后的图像分别按照从左到右侧的顺序示意于图2(a)中。从图2(a)可以看出，通过本发明的方法去噪后的图像大大降低图像噪声的同时，保留了图像细节，大大提高了图像质量。

将无噪声模拟数据T2加入15％莱斯噪声，采用本发明的方法对该噪声水平为15％的噪声图像进行去噪处理。模拟数据T2的图像、与T2对应噪声水平为15％的噪声图像及采用本发明的方法去噪后的图像分别按照从左到右侧的顺序示意于图2(b)中。从图2(b)可以看出，通过本发明的方法去噪后的图像大大降低图像噪声的同时，保留了图像细节，大大提高了图像质量。

将PD加权的无噪声模拟数据(简称模拟数据PD)加入15％莱斯噪声，采用本发明的方法对该噪声水平为15％的噪声图像进行去噪处理。模拟数据PD的图像、与PD对应噪声水平为15％的噪声图像及采用本发明的方法去噪后的图像分别按照从左到右侧的顺序示意于图2(c)中。从图2(c)可以看出，通过本发明的方法去噪后的图像大大降低图像噪声的同时，保留了图像细节，大大提高了图像质量。

采用本发明的方法对如图3第一行所示的噪声水平为3％的真实数据进行去噪处理，去噪后的对应的图像如图3第二行所示，对应的残差图像如图3第三行所示。残差图像是原始噪声图像减去去噪后的图像。图4是与图3对应的局部放大图，第一行表示原始噪声图像的局部放大图，第二行表示使用本发明方法去噪后的图像的局部放大图。从图3中可以看出，通过本发明的方法所去除的噪声当中几乎不含结构信息，说明去噪的同时没有将结构信息滤除。从图4的局部放大图可以看出，通过本发明方法进行去噪，原本模糊不清的结构信息变得清晰可见。可见，通过本发明的方法去噪后的图像大大降低图像噪声的同时，保留了图像细节，大大提高了图像质量。

采用本发明的方法对如图5第一行所示的噪声水平为4.5％的真实数据进行去噪处理，去噪后的对应的图像如图5第二行所示，对应的残差图像如图5第三行所示。残差图像是原始噪声图像减去去噪后的图像。图6是与图5对应的局部放大图，第一行表示原始噪声图像的局部放大图，第二行表示使用本发明方法去噪后的图像的局部放大图。从图5中可以看出，通过本发明的方法所去除的噪声当中几乎不含结构信息，说明去噪的同时没有将结构信息滤除。从图6的局部放大图可以看出，通过本发明方法进行去噪，原本模糊不清的结构信息变得清晰可见。可见，通过本发明的方法去噪后的图像大大降低图像噪声的同时，保留了图像细节，大大提高了图像质量。

以上实施例仅用以说明本发明的技术方案而非对本发明保护范围的限制，尽管参照较佳实施例对本发明作了详细说明，本领域的普通技术人员应当理解，可以对本发明的技术方案进行修改或者等同替换，而不脱离本发明技术方案的实质和范围。

Claims (7)

1.一种基于高阶奇异值分解的磁共振图像去噪方法，依次包括如下步骤：

(1)对加有莱斯噪声的原始磁共振图像进行方差稳定变换，将依赖信号分布的莱斯噪声变成独立于信号分布的噪声，得到变换后的噪声图像；

(2)将变换后的噪声图像中的每一个像素作为标的像素或者按照一定距离间隔取像素作为标的像素，取标的像素及其周围的像素作为参考块，执行以下操作进行第一次高阶奇异值分解去噪，具体是：

(a)通过k最近邻方法依次逐个寻找参考块的相似块，组成高维数组；

(b)对高维数组进行高阶奇异值分解变换，得到对应的系数和自适应基；

(c)通过第一阈值操作，将绝对值小于第一阈值的系数置零得到修正后的系数；

(d)根据自适应基与修正后的系数，进行高阶奇异值分解逆变换，将逆变换后的结果作为该高维数组中所有图像块去除噪声后的估计值；

(3)对步骤(2)得到的估计值加权平均进行像素合并，得到像素合并后的估计值，对每一个像素进行像素合并得到第一次去噪后的图像；

(4)对方差稳定变换后的噪声图像和第一次去噪后的图像进行加权平均，得到第一次加权平均后的图像；

(5)对第一次加权平均后的图像进行第二次高阶奇异值分解去噪，具体是：

(5.1)将第一次加权平均后的图像中的每一个像素作为二次标的像素或者按照一定距离取像素作为二次标的像素，取二次标的像素及其周围的像素作为二次参考块，执行以下操作：

(a)以第一次去噪后的图像进行相似块判断，通过k最近邻方法依次逐个寻找二次参考块的二次相似块，组成二次高维数组；

(b)对二次高维数组进行高阶奇异值分解变换，得到对应的二次系数和二次自适应基；

(c)通过二次阈值操作，将绝对值小于二次阈值的二次系数置零得到修正后的二次系数；

(d)根据二次自适应基与修正后的二次系数，进行高阶奇异值分解逆变换，将逆变换后的结果作为该高维数组中所有图像块去除噪声后的二次估计值；

(5.2)对步骤(5.1)得到的二次估计值加权平均进行像素合并，得到像素合并后的二次估计值，对每一个像素进行像素合并得到第二次去噪后的图像；

(6)对第二次去噪后的图像进行方差稳定逆变换，得到最终的滤波图像。

2.根据权利要求1所述的基于高阶奇异值分解的磁共振图像去噪方法，其特征在于：所述步骤(2)中(c)分步骤的第一阈值大小其中p为块大小，d为图像维度，K为相似块的个数，σ为方差稳定变换后图像的标准差。

3.根据权利要求2所述的基于高阶奇异值分解的磁共振图像去噪方法，其特征在于：所述步骤(4)具体是采用迭代正则化方法将方差稳定变换后的噪声图像和第一次去噪后的图像进行加权平均。

4.根据权利要求3所述的基于高阶奇异值分解的磁共振图像去噪方法，其特征在于：所述步骤(4)对方差稳定变换后的噪声图像和第一次去噪后的图像进行加权平均具体为：

$Y_w={\hat{X}}_{HOSVD1}+\mathrm{\ξ}(Y-{\hat{X}}_{HOSVD1})=(1-\mathrm{\ξ}){\hat{X}}_{HOSVD1}+\mathrm{\ξ}Y\mathrm{......}\left(I\right);$

其中表示第一次去噪后的图像，Y表示方差稳定变化后的噪声图像，ξ∈[0，1]是松驰参数，表示噪声图像Y对图像Y_w的贡献。

5.根据权利要求4所述的基于高阶奇异值分解的磁共振图像去噪方法，其特征在于：

所述步骤(5)中(c)分步骤的二次阈值大小

其中γ∈[0，1]为尺度参数，控制图像Y_w的标准差σ_w大小。

6.根据权利要求1至5任意一项所述的基于高阶奇异值分解的磁共振图像去噪方法，其特征在于：

所述步骤(2)的(a)中，通过k最近邻方法依次逐个寻找参考块的相似块，相似块的块大小和相似块个数的取值根据噪声水平而定；

当噪声水平不大于2％时，相似块的块大小为3，相似块个数为30—70；

当噪声水平大于2％而小于等于7％时，相似块的块大小为4，相似块个数为70—80；

当噪声水平大于7％而小于等于15％时，相似块的块大小为5，相似块个数为100—130。

7.根据权利要求1至5任意一项所述的基于高阶奇异值分解的磁共振图像去噪方法，其特征在于：

所述步骤(5)的(a)中，通过k最近邻方法依次逐个寻找二次标的像素参考块的二次相似块，二次相似块的块大小和二次相似块个数的取值根据噪声水平而定；

当噪声水平不大于2％时，二次相似块的块大小为3，二次相似块个数为30—70；

当噪声水平大于2％而小于等于7％时，二次相似块的块大小为4，二次相似块个数为70—80；

当噪声水平大于7％而小于等于15％时，二次相似块的块大小为5，二次相似块个数为100—130。