CN113724158B - 一种动态对比增强磁共振成像的降噪方法 - Google Patents

Abstract

本发明公开了一种动态对比增强磁共振成像的降噪方法：(1)获得添加噪声后的动态对比增强磁共振图像数据I_noise；(2)KIBF利用I_noise中待滤波像素点与其周围相邻像素点的信号强度差距得到信号强度加权矩阵、空间距离大小得到空间距离加权矩阵，两者相乘再归一化得整体加权矩阵，再与I_noise中周围相邻像素点信号对应相乘加和后获得降噪后的I_noise；(3)根据降噪前后的均方误差确定KIBF的L、σ_s和σ_z，得到优化后的KIBF；(4)采用动态增强磁共振成像实验数据和优化后的KIBF分别替换(2)中的I_noise和KIBF，再进行(2)，获得降噪后的动态增强磁共振成像实验数据。该降噪方法不仅可以去除动态对比增强磁共振成像中的噪声，同时保留图像的药代动力学信息和精细的结构细节。

Description

一种动态对比增强磁共振成像的降噪方法

技术领域

本发明涉及磁共振成像的技术领域，特别涉及一种动态对比增强磁共振成像的降噪方法。

背景技术

动态对比度增强磁共振成像可以用于测量动物及人体组织各成分占比、微血管渗透性、血流动力学等灌注参数。该技术在注射造影剂之前，期间和之后进行连续成像，进而利用药代动力学模型对数据进行分析，从而得到定量生理参数，包括血管密度、通透性等。它通常用于增强不同组织，特别是正常和病灶之间的对比度。动态对比度增强磁共振成像的临床应用包括肿瘤的检测和表征(例如，乳腺，肝，脑，前列腺中的肿瘤)，以及心脏梗塞的诊断和血管造影。

尽管磁共振影像在采集速度、空间分辨率和信噪比方面取得了显著的进步，但是依旧会受到噪声的影响。磁共振图像中的噪声有不同的来源，如硬件引起的噪声、几何畸变和磁场强度不均匀性引起的噪声，以及在图像采集过程中病人运动产生的移动伪影。然而，利用药代动力学模型对动态对比度增强磁共振成像进行分析时，噪声可能引起药代动力学模型参数估计产生偏差，这对利用动态对比度增强磁共振成像进行疾病诊断是十分不利的。如公开号为CN109003232A的中国专利公开了基于频域尺度平滑Shearlet的医学MRI图像去噪方法：基于频域尺度平滑Shearlet去噪方法，提出了新的医学MRI图像Racian噪声模型，然后进行频域平滑多尺度分解与多方向分解，对分解后的高频尺度平滑Shearlet系数再进行新的统一阈值处理，然后通过逆Shearlet变换生成去噪后的MRI图像。

因此，如何去除动态对比度增强磁共振成像中的噪声，同时保留药代动力学增强信息和精细的结构细节，并尽量减少引入新的偏差(例如，过度平滑引起的偏差)，是目前本领域技术人员亟需解决的问题。

发明内容

本发明的目的在于提供一种动态对比增强磁共振成像的降噪方法，本发明提供的降噪方法，不仅可以去除动态对比增强磁共振成像中的噪声，同时保留图像的药代动力学信息和精细的结构细节。

一种动态对比增强磁共振成像的降噪方法，所述降噪方法包括以下步骤：

(1)获得添加噪声后的动态对比增强磁共振图像数据I_noise；

(2)动力学诱导双边滤波器KIBF利用I_noise中待滤波像素点与其周围相邻像素点的信号强度差距得到信号强度加权矩阵，利用I_noise中待滤波像素点与其周围相邻像素点的空间距离大小得到空间距离加权矩阵，将信号强度加权矩阵与空间距离加权矩阵相乘再归一化得整体加权矩阵，将整体加权矩阵与I_noise中周围相邻像素点信号对应相乘加和后获得降噪后的动态增强磁共振成像数据；

(3)根据降噪后的动态增强磁共振成像数据与未添加噪声的动态对比增强磁共振图像数据的均方误差确定动力学诱导双边滤波器KIBF的最佳三维空间过滤窗Ω边长或直径L、空间距离加权系数σ_s和信号强度加权系数σ_z，得到优化后的动力学诱导双边滤波器KIBF；

(4)采用动态增强磁共振成像实验数据和步骤(3)中优化后的动力学诱导双边滤波器KIBF分别替换步骤(2)中的I_noise和动力学诱导双边滤波器KIBF，再进行步骤(2)，获得降噪后的动态增强磁共振成像实验数据。

在步骤(1)中，获得添加噪声后的动态对比增强磁共振图像数据I_noise的方法为：

(1-1)对动态对比增强磁共振图像实验数据进行归一化操作；

(1-2)使用药代动力学模型对步骤(1-1)处理后的动态对比增强磁共振图像实验数据进行逐点模型拟合，得到每个像素点每帧的拟合残差；

(1-3)根据步骤(1-2)估计的拟合残差，计算每帧感兴趣区域内拟合残差的标准差，即得到每帧噪音的标准差SD(m)，其中m为动态增强磁共振成像数据的当前帧数；

(1-4)使用药代动力学模型生成包含一定空间信息的动态对比增强磁共振图像模拟数据I_real，按步骤(1-3)中得到的噪音标准差SD(m)对I_real的每帧图像添加Rician噪音，得到添加噪声后的模拟数据I_noise。

优选的，所述归一化的具体方法为：动态对比增强磁共振图像实验数据中每个像素点的时间序列信号除以该像素点未注入造影剂之前的时间序列信号的均值。

优选的，所述药代动力学模型选自Tofts模型、extended-Tofts模型、Patlak模型或shutter-speed模型。

在步骤(2)中，利用I_noise中待滤波像素点与其周围相邻像素点的空间距离大小得到空间距离加权矩阵的方法为：

选择I_noise的像素点，将动力学诱导双边滤波器KIBF过滤窗Ω的中心点设置为该像素点，计算在空间过滤窗内的像素点与中心待过滤像素点的空间距离，得到三维空间距离加权矩阵w_space，得到空间距离加权矩阵公式如下：

其中，i代表空间过滤窗Ω内中心的待过滤像素点空间坐标位置，j代表空间过滤窗Ω内的任一像素点的空间坐标位置。

在步骤(2)中，利用I_noise中待滤波像素点与其周围相邻像素点的信号强度差距得到信号强度加权矩阵的方法为：

计算在空间过滤窗内的任意像素点与中心待过滤像素点的信号强度差距，得到三维信号强度加权矩阵w_signal，得到信号强度加权矩阵公式如下：

其中，i代表空间过滤窗Ω内中心的待过滤像素点空间坐标位置，j代表空间过滤窗Ω内的任一像素点的空间坐标位置，Z(k)代表像素点k的动态增强磁共振成像时序信号，Z(k)＝{x(k,t),t＝1,2,…,T},x为动态增强磁共振成像信号强度，t为动态增强磁共振成像采集的帧数，T代表动态对比增强磁共振成像采集的总帧数，||.||代表欧氏距离计算公式。

在步骤(2)中，将信号强度加权矩阵与空间距离加权矩阵相乘再归一化得整体加权矩阵的方法为：

将空间距离加权矩阵w_space与信号强度加权矩阵w_signal相乘并归一化得到整体加权矩阵w，得到整体加权矩阵的公式如下：

其中，S(i)是空间距离加权矩阵w_space与信号强度加权矩阵w_signal相乘之后对三维矩阵的每个元素进行加和，S(i)公式如下，

在步骤(2)中，将整体加权矩阵与I_noise中周围相邻像素点信号对应相乘加和后获得降噪后的动态增强磁共振成像数据的方法为：

将I_noise中每帧图像对应的空间过滤窗内像素信号强度与整体加权矩阵w对应相乘，得到空间过滤窗中心的待降噪像素点降噪后的动态对比增强磁共振时序信号I_denoised(i)，具体公式如下：

其中，j代表空间过滤窗Ω内某像素点，Z(j)代表空间过滤窗Ω中j像素点降噪前的动态对比增强磁共振时序信号；

移动空间过滤窗到下一个像素点，重复步骤(2)，直到完成动态对比增强磁共振成像图像中所有像素点的降噪，得到整体的降噪图像。

在步骤(3)中，计算降噪后动态增强磁共振成像数据与I_real的均方误差MSE，均方误差的计算公式如下：

其中，M是动态对比增强磁共振图像圈定的感兴趣区域的像素点个数。

在步骤(3)中，确定动力学诱导双边滤波器KIBF的最佳三维空间过滤窗Ω边长(或直径)L，空间距离加权系数σ_s和信号强度加权系数σ_z的方法为：

(3-1)在5-30mm范围内变化三维空间过滤窗Ω的边长或直径L，重复步骤(2)，计算当前L下的均方误差MSE并统计运行时间，以运行时间为横坐标，MSE为纵坐标，则得到的应为L型曲线，选取L型曲线拐点处边长或直径L作为最佳三维空间过滤窗Ω的边长或直径L；

(3-2)，选用步骤(3-1)中确定的最佳三维空间过滤窗Ω的边长或直径L，在0-10范围内变化空间距离加权系数σ_s和信号强度加权系数σ_z，重复步骤6-11，计算当前σ_s和σ_z组合下的均方误差MSE；

(3-3)，根据最小化算法，重复步骤(3-2)，直到找到阈值范围内最小MSE，获得最小MSE下的σ_s和σ_z组合，即为最佳的σ_s和σ_z组合。

本发明提供的动态对比增强磁共振成像的降噪方法不仅可以去除动态对比增强磁共振成像中的噪声，同时保留图像的药代动力学信息和精细的结构细节。

附图说明

图1为本发明方法总流程图；

图2为动力学诱导的双边滤波器(KIBF)降噪方法流程图；

图3为模拟动态对比增强磁共振图像数据中确定最佳三维空间过滤窗Ω过程中的KIBF运行时间-MSE曲线图；

图4为模拟动态对比增强磁共振图像数据中不同σ_s和σ_z组合的MSE图；

图5为真实动态对比增强磁共振图像某像素点的时间-信号曲线降噪效果对比；

图6为真实动态对比增强磁共振图像降噪前后图像利用extended-Tofts模型拟合后的参数结果；

图7为降噪后图像SNR和CNR的统计结果。

具体实施方式

为使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例对本发明进行进一步的详细说明。应当理解，此处所描述的具体实施方式仅仅用以解释本发明，并不限定本发明的保护范围。

作为具体实施范例，将本发明提供的动态对比增强磁共振成像的降噪方法应用在一个胶质瘤患者的动态对比增强磁共振图像上，其总操作流程见图1，具体的KIBF算法流程图见图2，具体包括以下步骤：

步骤1，对动态对比增强磁共振图像实验数据进行归一化操作，归一化的具体方法为，每个像素点的时间序列信号除以该像素点未注入造影剂之前的时间序列信号的均值。

步骤2，使用药代动力学模型extended-Tofts模型对步骤1处理后的动态对比增强磁共振图像实验数据进行逐点模型拟合，得到每个像素点每帧的拟合残差。

步骤3，根据步骤2估计的拟合残差，计算每帧感兴趣区域(可以为整个头部)内拟合残差的标准差，即得到每帧噪音的标准差SD(m)，其中m为动态增强磁共振成像数据的当前帧数。

步骤4，使用药代动力学模型生成包含一定空间信息的动态对比增强磁共振图像模拟数据I_real，按步骤3中得到的噪音标准差SD(m)对I_real的每帧图像添加Rician噪音，得到添加噪声后的模拟数据I_noise。

步骤5，选择动力学诱导双边滤波器(KIBF)的的三维空间过滤窗Ω的大小，三维空间过滤窗Ω应为正方体，边长设为L；设置空间距离加权系数σ_s和信号强度加权系数σ_z。

步骤6，选择I_noise的像素点，将动力学诱导双边滤波器(KIBF)过滤窗Ω的中心点设置为该像素点，计算在空间过滤窗内的像素点与中心待过滤像素点的空间距离，得到三维空间距离加权矩阵w_space，得到空间距离加权矩阵公式如下：

其中，i代表空间过滤窗Ω内中心的待过滤像素点空间坐标位置，j代表空间过滤窗Ω内的任一像素点的空间坐标位置。

步骤7，计算在空间过滤窗内的任意像素点与中心待过滤像素点的信号强度差距，得到三维信号强度加权矩阵w_signal，得到信号强度加权矩阵公式如下：

其中，i代表空间过滤窗Ω内中心的待过滤像素点空间坐标位置，j代表空间过滤窗Ω内的任一像素点的空间坐标位置，Z(k)代表像素点k的动态增强磁共振成像时序信号，Z(k)＝{x(k,t),t＝1,2,…,T},x为动态增强磁共振成像信号强度，t为动态增强磁共振成像采集的帧数，T代表动态对比增强磁共振成像采集的总帧数，||.||代表欧氏距离计算公式。

步骤8，将空间距离加权矩阵w_space与信号强度加权矩阵w_signal相乘并归一化得到整体加权矩阵w，得到整体加权矩阵的公式如下：

其中，S(i)是空间距离加权矩阵w_space与信号强度加权矩阵w_signal相乘之后对三维矩阵的每个元素进行加和，S(i)公式如下，

步骤9，将I_noise中每帧图像对应的空间过滤窗内像素信号强度与整体加权矩阵w对应相乘，得到空间过滤窗中心的待降噪像素点降噪后的动态对比增强磁共振时序信号I_denoised(i)，具体公式如下：

其中，j代表空间过滤窗Ω内某像素点，Z(j)代表空间过滤窗Ω中j像素点降噪前的动态对比增强磁共振时序信号。

步骤10，移动空间过滤窗到下一个像素点，重复步骤6到步骤9，直到完成动态对比增强磁共振成像图像中所有像素点的降噪，得到整体的降噪图像。

步骤11，计算降噪后动态增强磁共振成像数据与I_real的均方误差(MSE)，均方误差的计算公式如下：

其中，M是动态对比增强磁共振图像圈定的感兴趣区域的像素点个数。

步骤12，在合理范围内变化三维空间过滤窗Ω的边长(或直径)L，重复步骤6-11，计算当前L下的均方误差MSE并统计运行时间，以运行时间为横坐标，MSE为纵坐标，得到图3，选取L型曲线拐点处边长作为最佳三维空间过滤窗Ω的边长(或直径)，在本实施例中L为15毫米。

步骤13，选用步骤12中确定的最佳三维空间过滤窗Ω的边长15毫米，在合理范围内变化空间距离加权系数σ_s和信号强度加权系数σ_z，重复步骤6-11，计算当前σ_s和σ_z组合下的均方误差MSE。

步骤14，根据最小化算法，重复步骤13，直到找到阈值范围内最小MSE，获得最小MSE下的σ_s和σ_z组合，即为动态增强磁共振成像实验数据的最佳的σ_s和σ_z组合，如图4所示，本实例最佳σ_s和σ_z组合的分别为6和1.5。

步骤15，根据步骤12中得到的最佳三维空间过滤窗Ω的边长15毫米和步骤14得到σ_s＝6和σ_z＝1.5，用动态增强磁共振成像实验数据替换步骤6-10中的I_noise并重新进行步骤6-10，获得降噪后的动态增强磁共振成像实验数据。

为展示本发明提出的方法在动态对比增强磁共振图像降噪方面的效果，下面将结合附图说明该具体实施范例的实验结果：

图5展示了本实施例中动态对比增强磁共振图像中的某个像素点的降噪前，利用常用的高斯滤波器(Gaussian Filter)降噪(σ＝1.5)以及本发明提供的动力学诱导双边滤波器(KIBF)方法降噪的时间-信号曲线，可以看到本发明提供的方法，相对于原始信号，有很好的降噪效果，同时保留了信号的细节，而常用的高斯滤波器(Gaussian Filter)降噪方法相比于原始信号产生比较大的偏差。

图6展示了利用extended-Tofts模型分别处理原始数据，高斯滤波器(GaussianFilter)降噪后数据和本发明提供的动力学诱导双边滤波器(KIBF)方法降噪的数据后得到的参数图结果，K^trans为造影剂转移常数，v_p为血液组分占比，v_e为组织间隙组分占比，v_e在正常组织区域由造影剂泄露较少难以估计，因此只展示了造影剂泄露较多的肿瘤区域。从图6中可以看出，本发明提供的动力学诱导双边滤波器(KIBF)方法在这三个参数图上参数分布更加平滑，同时相比于高斯滤波器(Gauss Filter)的结果能更多的保留图像的细节。

图7展示了降噪后图像SNR和CNR的统计结果。SNR和CNR的公式如下，可以看到本发明提供的动力学诱导双边滤波器(KIBF)方法在SNR和CNR相比于原始图像和高斯滤波器(Gaussian Filter)过滤后的图像都有较大的提升。

Claims (6)

1.一种动态对比增强磁共振成像的降噪方法，其特征在于，所述降噪方法包括以下步骤：

(1)获得添加噪声后的动态对比增强磁共振图像数据I_noise；

(2)动力学诱导双边滤波器KIBF利用I_noise中待滤波像素点与其周围相邻像素点的信号强度差距得到信号强度加权矩阵，利用I_noise中待滤波像素点与其周围相邻像素点的空间距离大小得到空间距离加权矩阵，将信号强度加权矩阵与空间距离加权矩阵相乘再归一化得整体加权矩阵，将整体加权矩阵与I_noise中周围相邻像素点信号对应相乘加和后获得降噪后的动态增强磁共振成像数据；

(3)根据降噪后的动态增强磁共振成像数据与未添加噪声的动态对比增强磁共振图像数据的均方误差确定动力学诱导双边滤波器KIBF的最佳三维空间过滤窗Ω边长或直径L、空间距离加权系数σ_s和信号强度加权系数σ_z，得到优化后的动力学诱导双边滤波器KIBF；

(4)采用动态增强磁共振成像实验数据和步骤(3)中优化后的动力学诱导双边滤波器KIBF分别替换步骤(2)中的I_noise和动力学诱导双边滤波器KIBF，再进行步骤(2)，获得降噪后的动态增强磁共振成像实验数据；

在步骤(1)中，获得添加噪声后的动态对比增强磁共振图像数据I_noise的方法为：

(1-1)对动态对比增强磁共振图像实验数据进行归一化操作；

(1-2)使用药代动力学模型对步骤(1-1)处理后的动态对比增强磁共振图像实验数据进行逐点模型拟合，得到每个像素点每帧的拟合残差；

(1-3)根据步骤(1-2)估计的拟合残差，计算每帧感兴趣区域内拟合残差的标准差，即得到每帧噪音的标准差SD(m)，其中m为动态增强磁共振成像数据的当前帧数；

(1-4)使用药代动力学模型生成包含一定空间信息的动态对比增强磁共振图像模拟数据I_real，按步骤(1-3)中得到的噪音标准差SD(m)对I_real的每帧图像添加Rician噪音，得到添加噪声后的模拟数据I_noise；

在步骤(2)中，利用I_noise中待滤波像素点与其周围相邻像素点的空间距离大小得到空间距离加权矩阵的方法为：

选择I_noise的像素点，将动力学诱导双边滤波器KIBF过滤窗Ω的中心点设置为该像素点，计算在空间过滤窗内的像素点与中心待过滤像素点的空间距离，得到三维空间距离加权矩阵w_space，得到空间距离加权矩阵公式如下：

其中，i代表空间过滤窗Ω内中心的待过滤像素点空间坐标位置，j代表空间过滤窗Ω内的任一像素点的空间坐标位置；

在步骤(2)中，利用I_noise中待滤波像素点与其周围相邻像素点的信号强度差距得到信号强度加权矩阵的方法为：

计算在空间过滤窗内的任意像素点与中心待过滤像素点的信号强度差距，得到三维信号强度加权矩阵w_signal，得到信号强度加权矩阵公式如下：

其中，i代表空间过滤窗Ω内中心的待过滤像素点空间坐标位置，j代表空间过滤窗Ω内的任一像素点的空间坐标位置，Z(k)代表像素点k的动态增强磁共振成像时序信号，Z(k)＝{x(k,t),t＝1,2,…,T},x为动态增强磁共振成像信号强度，t为动态增强磁共振成像采集的帧数，T代表动态对比增强磁共振成像采集的总帧数，||.||代表欧氏距离计算公式；

在步骤(2)中，将整体加权矩阵与I_noise中周围相邻像素点信号对应相乘加和后获得降噪后的动态增强磁共振成像数据的方法为：

将I_noise中每帧图像对应的空间过滤窗内像素信号强度与整体加权矩阵w对应相乘，得到空间过滤窗中心的待降噪像素点降噪后的动态对比增强磁共振时序信号I_denoised(i)，具体公式如下：

其中，j代表空间过滤窗Ω内某像素点，Z(j)代表空间过滤窗Ω中j像素点降噪前的动态对比增强磁共振时序信号；

移动空间过滤窗到下一个像素点，重复步骤(2)，直到完成动态对比增强磁共振成像图像中所有像素点的降噪，得到整体的降噪图像。

2.根据权利要求1所述的动态对比增强磁共振成像的降噪方法，其特征在于，所述归一化的具体方法为：动态对比增强磁共振图像实验数据中每个像素点的时间序列信号除以该像素点未注入造影剂之前的时间序列信号的均值。

3.根据权利要求1所述的动态对比增强磁共振成像的降噪方法，其特征在于，所述药代动力学模型选自Tofts模型、extended-Tofts模型、Patlak模型或shutter-speed模型。

4.根据权利要求1所述的动态对比增强磁共振成像的降噪方法，其特征在于，在步骤(2)中，将信号强度加权矩阵与空间距离加权矩阵相乘再归一化得整体加权矩阵的方法为：

将空间距离加权矩阵w_space与信号强度加权矩阵w_signal相乘并归一化得到整体加权矩阵w，得到整体加权矩阵的公式如下：

其中，S(i)是空间距离加权矩阵w_space与信号强度加权矩阵w_signal相乘之后对三维矩阵的每个元素进行加和，S(i)公式如下，

5.根据权利要求1所述的动态对比增强磁共振成像的降噪方法，其特征在于，在步骤(3)中，计算降噪后动态增强磁共振成像数据与I_real的均方误差MSE，均方误差的计算公式如下：

其中，M是动态对比增强磁共振图像圈定的感兴趣区域的像素点个数。

6.根据权利要求1所述的动态对比增强磁共振成像的降噪方法，其特征在于，在步骤(3)中，确定动力学诱导双边滤波器KIBF的最佳三维空间过滤窗Ω边长或直径L，空间距离加权系数σ_s和信号强度加权系数σ_z的方法为：

(3-1)在5-30mm范围内变化三维空间过滤窗Ω的边长或直径L，重复步骤(2)，计算当前L下的均方误差MSE并统计运行时间，以运行时间为横坐标，MSE为纵坐标，则得到的应为L型曲线，选取L型曲线拐点处边长或直径L作为最佳三维空间过滤窗Ω的边长或直径L；

(3-2)，选用步骤(3-1)中确定的最佳三维空间过滤窗Ω的边长或直径L，在0-10范围内变化空间距离加权系数σ_s和信号强度加权系数σ_z，重复步骤(2)和步骤(3)中的计算降噪后动态增强磁共振成像数据与I_real的均方误差MSE，计算当前σ_s和σ_z组合下的均方误差MSE；

(3-3)，根据最小化算法，重复步骤(3-2)，直到找到阈值范围内最小MSE，获得最小MSE下的σ_s和σ_z组合，即为最佳的σ_s和σ_z组合。