基于卡方无偏估计的邻域收缩MRI去噪方法与流程

本发明属于图像处理技术领域。具体来说，涉及一种以去除MRI莱斯噪声，提高信噪比为目的的基于卡方无偏估计的邻域收缩MRI去噪方法。

背景技术：

现代医学中，磁共振成像(MRI)作为一项基本的医学成像技术，已经成为了医生诊断和治疗的重要辅助手段，磁共振图像在获取过程中由于硬件电路和人体的因素会引入噪声，噪声降低了MRI的质量，使得一些组织的边界变得模糊，细微结构难以辨别，增加了对图像细节的识别分析和处理的难度，影响医学诊断。因此，针对噪声MRI用软件去噪是一件实用方便且性价比高的事情。

常用的MRI获取方法是对频率域(k域)内感兴趣的目标进行采样，然后对频率域数据进行离散傅里叶反变换，得到包含实部和虚部的复数图像数据，再转换为幅值数据，在实部和虚部引入的加性高斯白噪声，转化为幅值图像后，呈现莱斯(Rician)分布。当信噪比较高时，莱斯分布倾向于高斯分布，当信噪比低时，倾向于瑞利分布，因此针对低信噪比的MRI，用传统的去除高斯噪声的方法进行噪声估计和去噪会产生较大的偏差。

目前，MRI去噪中一个重要问题是保留有用的边缘、细节信息，因为医学图像事关病患安全，任何微小的失误都是不允许的。MRI去噪通常有三种方式，一是在幅值域针对莱斯分布的模型去噪，二是在幅值的平方域去噪，也有对三维的MRI进行去噪的。1999年Nowak提出，平方使得噪声图像产生了一个与信号无关的、加性的偏置，经小波变换后该偏置在尺度系数中保留，直接减去偏置即可提升去噪效果，而且平方幅值图像服从两个自由度的非中心卡方分布。Nowak将平方图像的尺度系数去除偏置，对高频小波系数用卡方分布的性质计算信号方差，再应用类似维纳滤波的逐点比例收缩方法，得到优于传统小波收缩的去噪效果。Anand等人(2010年)在Nowak的基础上，将平方幅值图像经平稳小波变换后，高低频分别采用NeighShrinkSURE和双边滤波器进行去噪，对高频系数的处理不再是逐点处理，而是考虑了像素的邻域。在低频系数中，含有大部分的图像能量，信噪比高，而双边滤波器良好的保持边缘的能力，可以去模糊，该方法较好的滤除了噪声。Luisier等人(2012年)在论文“CURE for noisy Magnetic Resonance Imgaes：Chi-square unbiased risk estimation”针对平方MRI的非中心卡方分布模型，推导出针对卡方分布的均方误差的无偏估计表达式(CURE)，应用到两种线性阈值扩展去噪方法(LET)中，其一为图像域的CURE估计确定小波域UWT变换下的LET参数，其二为用小波域的CURE估计确定小波域的LET的参数，由于非中心卡方分布的特殊性质，小波域的CURE估计只能应用在未归一化的哈尔(Haar)小波变换。

技术实现要素：

为了去除噪声，提高MRI的信噪比，本发明推导了在小波域仅取决于小波系数的阈值收缩方法的CURE估计表达式，使其可以适用于邻域收缩去噪方法上，结合低频部分的双边滤波器处理和快速的循环移位技术，取得了较好的MRI去噪效果。

基于卡方无偏估计的邻域收缩MRI去噪方法，包括如下步骤：

(1)在MRI的背景区域估计噪声标准差σ；

(2)对噪声图像m进行平方，然后除以噪声标准差σ的平方，得图像y

(3)对y进行未归一化的平稳哈尔小波变换(SWT)，分解L层，得高频系数和低频系数，高频系数有L*3个子带；

(4)对低频系数用双边滤波器去模糊，对高频系数用基于卡方无偏估计的邻域收缩法(NeighShrinkCURE)去噪；

(5)对去噪后的小波系数连同低频系数进行循环移位、小波反变换(ISWT)，得到若干循环移位的去噪结果；

(6)将上述若干循环移位的去噪结果通过非线性映射、求平均得到最终的去噪图像。

所述的NeighShrinkCURE包括两个部分，其一是小波域的取决于小波系数的阈值收缩方法的CURE估计表达式的推导，其二是邻域收缩方法的CURE估计值的计算。

1、小波域关于小波系数的阈值收缩方法的CURE估计表达式的推导过程如下：

如果用表示一个N_x×N_y大小的图像的索引集合，用μ＝[μ_n]_n∈Ω∈C^N表示理想(无噪)的大小为N＝N_x×N_y的MR复数图像，m＝[m_n]_n∈Ω∈C^N表示含噪声的MR图像，去噪问题转换为从m中估计μ的问题。在MR数据获取过程中，实部和虚部会分别引入标准差为σ的高斯噪声，其统计特性可表述为：

表示取复数的实部，表示取复数的虚部。我们的目标是估计μ_n的幅值，定义x和y如下：

则y服从非中心的卡方分布，非中心参数为x，自由度K＝2，去噪问题又可以转换为从非中心卡方分布的数据y中估计非中心参数x的问题。在设计去噪函数f(y)来估计x时，一个自然的评价标准是均方误差(MSE)：

从式中可以看出，计算MSE需要无噪信号x，实际情况中这是无法预知的，当即y服从非中心参数为x，自由度为K时，Luisier等人根据非中心卡方分布的统计特性，推导出了MSE的无偏估计一卡方无偏估计(CURE)，只要f(y)满足连续或是平滑的限制，就可以由f(y)计算MSE的无偏估计CURE，并将其推广到小波域，提出未归一化的哈尔小波变换下的CURE估计。本方法只考虑未归一化哈尔小波变换下的CURE估计，未归一化的哈尔小波的意思是基本哈尔小波变换的尺度函数小波函数而未归一化的哈尔小波尺度函数[1，1]，小波函数[1，-1]，这是为了使小波系数仍然满足非中心卡方分布的性质，在小波反变换时，重建滤波器也要做相应的改变。在Luisier论文中的相关定理是这样表述的，对图像y做未归一化哈尔小波变换，得到第j层的尺度系数s^j和小波系数w^j，用θ(w，s)＝θ^j(w^j，s^j)表示对无噪信号x的小波变换系数ω＝ω^j的估计，则第j层小波系数的MSE可以由下式估计：

式中K_j＝2^jK，K是非中心卡方分布的自由度，在磁共振仪器用单线圈采集成像时，K＝2。去噪函数取决于噪声小波系数w、尺度系数s，而很多常用的去噪算法诸如软阈值收缩、二变量收缩都是不考虑尺度系数的，在不考虑尺度系数的情况下，本发明推导出仅考虑小波系数w的阈值收缩方法的CURE估计，即用θ(w)代替θ(w，s)估计无噪信号的小波系数。

为简化问题，我们把图像展成一维的信号，噪声信号y、理想信号x通过第j层小波分解所得尺度系数分别为s^j、高频系数分别为则根据未归一化哈尔小波变换的性质可得：

式中K_j＝2^jK＝2^j+1。

第j层小波系数的MSE可表示为：

考虑第一层分解j＝1，噪声信号y＝s^j-1，理想信号K＝K_j/2，则term2可表述为：

E{ω_nθ_n(w)}＝E{x_2nθ_n(w)}-E{x_2n-1θ_n(w)} (7)

Luisier已经证明，当f(y)连续可导时：

套用公式(6)求(5)：

代入term2：

套用公式(10)求term3：

最后将term1、term2、term3代入(8)可得仅取决于小波系数w的阈值收缩方法的小波域CURE估计：

其余层的分解亦可以得到类似的证明，而二维图像的小波变换是对水平和垂直方向分别做小波变换，上述公式亦可直接推广到二维图像，所做的改动只是K_j＝2^2·jK。

2、如上述在推导出仅考虑小波系数w的阈值方法的小波域CURE估计之后，应用到邻域收缩去噪方法中，主要问题是阈值函数的导数计算以及阈值函数的平滑。

邻域收缩法考虑当前待处理的小波系数w_ij为中心的平方邻域B_ij，给定平方邻域系数和设λ为当前小波子带的全局阈值，为了表述方便，用w_n表示w_ij，S_n表示S_ij，则阈值收缩公式可表述为：

式中，θ(w_n)为去噪后的小波系数，用CURE估计来确定最佳的阈值λ，即λ_opt＝argminCURE(λ)，窗口大小取3*3，CURE估计需要计算θ(w_n)的一阶和二阶微分，令g＝max(u，0)，通过分层求导步骤如下：

CURE估计对去噪函数有平滑的要求，g＝max(u，0)须做平滑，利用不完全贝塔函数做平滑，设平滑的去噪函数为g＝f(u)＝u.*I(u)，其中

式中，a、b用于控制平滑的程度，为了保证平滑的曲线中心对称，a、b须取值相同，a＝b＝80。

3、所述的循环移位(Cycle Spinning，CS)是为了消除由于平稳小波变换缺乏平移不变性而产生的伪吉布斯现象。传统的循环平移(Cycle Spinning)，即对含噪信号进行“循环平移-阈值去噪-逆向循环平移”。由于对每次平移后的信号进行阈值去噪会使伪吉布斯现象出现在不同的地方，因此针对行和列方向上的每组平移量都会得到一个不同的去噪结果μ_i，j，i、j分别为行和列方向上的平移量，对所有去噪结果进行线性平均将得到抑制伪吉布斯现象的去噪结果。然而，传统平移方法每次平移都要调用去噪函数，耗费时间长，为了加快运行速度，本发明考虑了只针对去噪后的小波系数的循环平移，这样可以大大节省时间，因为平稳小波变换每个子带的尺寸和原图相同，对图像平移与对去噪系数的平移影响不大，尤其在图像尺寸较大时，仅在边缘处会有所不同，事实证明在平稳小波变换下，这种快速循环平移方法和基本方法结果是一样的，运行时间却大大减少。

附图说明

下面结合附图对本发明做进一步的说明：

图1基于卡方无偏估计的邻域收缩MRI去噪方法流程图；

图2邻域收缩法的MSE和CURE估计对比图；

图3传统循环平移和快速循环平移示意图；

图4 T2w，7％噪声第40切片噪声、无噪和无噪图像局部放大图

图5 T2w，7％噪声第40切片去噪效果对比图，包括去噪图像，残差图像(噪声图减去去噪图)和局部放大图像。

图6 T2w，7％噪声第40到第100切片的去噪指标对比

图7 T2w第40切片模拟加不同强度Rician噪声去噪指标对比

具体实施方式

本发明提供了一种基于卡方无偏估计的邻域收缩MRI去噪方法，当需要对噪声MRI进行去噪处理时，如图1所示依次执行如下步骤：

步骤1 在MRI的背景区域估计噪声标准差σ：

μ是选中的背景区域像素值的均值。

步骤2 对噪声图像m进行平方，然后除以噪声标准差σ的平方，得图像y；

y＝m²/σ² (2)

步骤3 对y进行未归一化的平稳哈尔小波变换，分解L层，得高频系数和低频系数，高频系数有L*3个子带；

步骤4 对低频系数用双边滤波器去模糊；对每个子带在预定的阈值搜索范围内用NeighShrinkCURE方法去噪，在阈值范围内用CURE估计来确定最佳的阈值λ，即λ_opt＝arg min CURE(λ)，再用进行去噪；

步骤5 对去噪后的小波系数连同低频系数进行循环移位，CS＝9，即分别向右移i位，向下移j位，i∈[0，2]，j∈[0，2]，进行小波反变换，再反平移相应的位数得

步骤6 由得

步骤7 对得到的求平均的最终的去噪图像

下面通过三组实验来详细分析该方法的性能：

从MRI数据库(http：//brainweb.bic.mni.mcgill.ca/brainweb/)下载T1加权(T1w)、T2加权(T2w)，和质子密度加权(PD)的不同噪声等级的8比特MRI，原始数据为三维立体的人脑模拟MRI扫描数据，分辨率为181*217*181，将本发明的方法(NeighShrinkCURE+BF/CS＝9)和Anand等人的双边滤波器结合基于Stein无偏估计的邻域收缩法(NeighShrinkSURE+BF)、Luisier等人的CURE-LET haar/CS＝16方法做对比。

实验1：取T1w，5％噪声、PD，7％噪声、T2w，7％噪声的MRI的第40个切片进行去噪，该切片含有较完整的人脑边缘和细节信息，用上述三种算法进行去噪，得去噪指标峰值信噪比(PSNR)和结构相似度(SSIM)对比如表1所示。其中T2w，7％的噪声图像如图3所示，去噪效果对比图像如图4所示。

表1不同去噪方法指标对比

实验2：对T2w，7％噪声的MRI的第40到第100个切片用上述三种算法分别进行去噪，去噪指标PSNR和SSIM对比如图5所示。

实验3：对人为加噪声的T2w，第40切片进行试验，将噪声等级从3％逐渐增大到15％，间隔2％，得到不同去噪方法的去噪指标对比。如图6所示

通过上述三组实验说明本发明算法能够很好的实现图像去噪，通过对比，证明本发明的算法基本能得到最高的PSNR和SSIM，视觉效果好，且鲁棒性强，适用不同噪声级的MRI去噪，能清晰地保留MRI的边缘和细节信息。