CN103985099B - 一种弥散张量磁共振图像张量域非局部均值去噪方法 - Google Patents

Abstract

本发明公开一种弥散张量磁共振图像张量域非局部均值去噪方法，属于数字图像处理和应用数学交叉学科技术领域，解决弥散张量磁共振图像易于被噪声影响的问题。首先，依次遍历该弥散张量磁共振图像的体素，并以遍历到的每个体素为中心设定相应搜索区域；然后将搜索区域内的所有体素逐个与中心体素进行张量矩阵相似性比较；最后，按照张量矩阵相似性高低赋予搜索区域内的体素不同权值，计算加权均值张量矩阵，获得该中心体素的去噪结果。解决了弥散张量磁共振图像易于被噪声影响的问题。

Description

一种弥散张量磁共振图像张量域非局部均值去噪方法

技术领域

一种弥散张量磁共振图像张量域非局部均值去噪方法，用于弥散张量磁共振图像的去噪，属于数字图像处理和应用数学交叉学科技术领域。

背景技术

弥散磁共振成像是对活体脑部组织的结构和生理机能进行无创研究的成像工具，可用于直接评价活体脑部纤维组织的生理机能，并使用纤维成像技术间接估计和重建脑部神经纤维的三维结构。弥散磁共振成像在中枢神经系统研究中得到深入应用，并可推广到其他人体纤维状组织，在心理、认知、临床医学等领域均具有巨大潜力。弥散张量磁共振图像是应用最广泛的弥散磁共振图像，其利用脑白质中水分子弥散效应对脑白质三维结构的每一体素使用一个3×3对称正定矩阵进行三维建模,该矩阵即为体素弥散张量。该张量矩阵通过对每一体素采集六个非共面的扩散敏感梯度磁场方向的回波衰减信号测量值构成的弥散加权图像和一个不施加扩散敏感梯度磁场的磁共振信号参考测量值解得。在实际扫描过程中，一般采集多于6个的扩散敏感方向信号以及在每个方向上采用更多的编码幅度，以此来抑制成像噪声。得到多于6个的扩散编码方向的测量结果后，可采用最小二乘拟合方法求解弥散矩阵，获得弥散张量磁共振图像。

由于弥散磁共振成像采用的快速扫描平面回波序列易受噪声干扰，弥散张量磁共振图像应用的首要预处理工作就是去噪。现有弥散磁共振去噪方法一般分为三类：第一类是对直接采集的弥散加权磁共振图像，使用传统时域或者频域手段进行去噪，该方法需要依次对多个不同方向的弥散加权磁共振图像处理，计算量过大；第二类是在估计弥散张量磁共振图像的张量矩阵过程中，对该过程进行正则化去噪，由于估计张量矩阵是一个非线性过程，对该过程的正则化容易产生过平滑结果；第三类是对弥散张量磁共振图像在张量域直接进行去噪，但是现有方法一般使用张量矩阵的部分信息，准确度较低。张量域去噪方法效果直观，便于应用，但是不同于普通MRI，由于弥散张量磁共振图像每个体素使用张量模型进行描述，在数学上表示为一个3×3对称正定矩阵，同时包含三维空间的弥散大小和方向信息，普通MRI或灰度图像使用的标量算子无法直接使用，因此弥散张量磁共振图像张量域去噪需要高阶算子实现张量之间的矩阵比较。

非局部均值去噪是近年得到有效应用的图像去噪算法，其使用像素或像素块的灰度相似性替代常用均值去噪的空间位置相似性，通过对高相似性像素(像素块)赋予高权值并通过加权平均的方法获取估计像素的最优值。但在弥散张量磁共振图像中，由于每个体素均为3×3对称正定矩阵描述的张量模型，现有非局部均值去噪使用的基于标量的灰度相似性无法直接应用。

发明内容

本发明针对现有技术的不足之处提供了一种弥散张量磁共振图像张量域非局部均值去噪方法，在弥散张量磁共振图像的张量域，根据非局部均值去噪原理，通过高阶算子对三维张量模型的几何结构及空间方向相似性的综合比较，对高相似性体素赋予高权值并进行加权均值方法获得去噪后张量数据，解决了弥散张量磁共振图像易于被噪声影响的问题。

为实现上述目的，本发明采用的技术方案为：

一种弥散张量磁共振图像张量域非局部均值去噪方法，其特征在于，如下步骤：

(1)在背景区域，计算背景区域的灰度直方图，并使用高斯函数拟合灰度直方图，根据高斯函数的方差确定噪声标准差，并设定衰减系数h和搜索区域半径Ω；

(2)根据噪声标准差、衰减系数h和搜索区域半径Ω依次遍历弥散张量磁共振图像中的所有体素，并以遍历到的每个体素为中心体素，设定半径为Ω的正方搜索区域Q；

(3)将正方搜索区域Q中的所有体素依次与中心体素进行张量矩阵相似性比较，获得描述张量矩阵相似性的相似性距离d；

(4)根据衰减系数h、正方搜索区域Q和张量矩阵相似性的相似性距离d获得归一化参数，计算公式如下：

$Z\left(p\right)=\underset{q\∈Q}{\mathrm{\Σ}}\exp (-\frac{d^2(p,q)}{h^2}),$

其中，p为中心体素，q为比较体素，d(p,q)为中心体素与比较体素的张量矩阵相似性的相似性距离；

(5)根据归一化参数Z(p)、衰减系数h和张量矩阵相似性的相似性距离d计算各体素的权值，计算公式如下：

$w(p,q)=\frac{1}{Z\left(p\right)}\exp (-\frac{d^2(p,q)}{h^2});$

(6)根据各体素权值w(p,q)，计算被遍历到体素的加权均值张量矩阵I，获得该中心体素的去噪结果，计算公式如下：

$\left(p\right)=\exp \left(\underset{q\∈Q}{\mathrm{\Σ}}w(p,q)\log \left(V\left(q\right)\right)\right),$

其中，I(p)为中心体素去噪后张量矩阵，V(q)为比较体素的张量矩阵；

(7)重复上述步骤，直到遍历完该弥散张量磁共振图像的所有体素。

作为优选，所述步骤(1)中，衰减系数h为估计噪声标准差的1-1.5倍。

作为优选，所述步骤(1)中，搜索区域半径Ω为整幅图像长度与宽度相比的较大值的10％-20％。

作为优选，所述步骤(3)中，使用黎曼距离描述张量矩阵相似性的相似性距离d，其公式为：

d(p,q)＝||log(V(p)^-1V(q)||；

其中，V(q)、V(p)分别为比较体素和中心体素的张量矩阵。

作为优选，所述步骤(3)中，使用对数欧式距离描述张量矩阵相似性的相似性距离d，其公式为：

$d(p,q)=\sqrt{\mathrm{Trace}\left[{(l\mathrm{og}\left(V\left(p\right)\right)-\log \left(V\left(q\right)\right))}^2\right]};$

其中，V(q)、V(p)分别为比较体素和中心体素的张量矩阵。

作为优选，所述步骤(3)中，使用偏张量商距离描述张量矩阵相似性的相似性距离d，其公式为：

$d(p,q)=\underset{i,j=1}{\overset{3}{\mathrm{\Σ}}}{V}_{\mathrm{ij}}\left(p\right)\·V_{\mathrm{ij}}\left(q\right)-(1/3)\mathrm{Trace}\left(V\left(p\right)\right)\·\mathrm{Trace}\left(V\left(q\right)\right);$

其中，V(q)、V(p)分别为比较体素和中心体素的张量矩阵，i、j分别为张量矩阵的行标和列标，V_ij(p)、V_ij(q)分别为体素p和q张量矩阵中第i行第j列元素。

作为优选，所述步骤(3)中，使用张量欧式距离描述张量矩阵相似性的相似性距离d，其公式为：

$d(p,q)=\sqrt{\mathrm{Trace}\left[{(V\left(p\right)-V\left(q\right))}^2\right]};$

其中，V(q)、V(p)分别为比较体素和中心体素的张量矩阵。

与现有技术相比，本发明的优点在于：

一、在弥散张量磁共振图像的张量域，通过对三维张量模型的几何结构及空间方向相似性的综合比较，对高相似性体素赋予高权值并进行加权均值方法获得去噪后张量数据，可在有效去噪的基础上，有效保持图像的几何形态结构。

二、衰减系数h为估计噪声标准差的1-1.5倍，该系数与指数曲线相关，控制均值滤波的平滑程度，防止了过小时，无法有效滤除噪声，过大时，模糊图像的问题。

三、搜索区域半径Ω为长度与宽度相比的较大值的10％-20％，在提高计算效率的同时考虑到了成像对象特征。

附图说明

图1为本发明的原理示意图；

图2为受噪声干扰图像示意图；

图3为本发明中使用黎曼距离本文方法去噪结果后的示意图；

图4为本发明中使用对数欧式距离本文方法去噪结果后的示意图；

图5为本发明中使用偏张量商距离本文方法去噪结果后的示意图；

图6为本发明中使用张量欧式距离本文方法去噪结果后的示意图。

具体实施方式

下面结合附图和实施例对本发明做进一步具体说明。

弥散张量磁共振图像采用标准扫描参数，使用3T的Philips InteraAchieva MRI扫描仪，数据集分辨率为128行、128列和53层，每个体素大小为2mm×2mm×2mm，共采集32个不同梯度方向(弥散加权值1000s/mm²)和1个未加权标准数据(弥散加权值0s/mm²)数据集，并使用最小二乘法计算每个体素的弥散张量矩阵V(p)和V(q)，设定半径为Ω的正方搜索区域Q，设定正方搜索区域Q为13行、13列和13层的正方体。

实施例1：

步骤1：在背景区域，计算背景区域的灰度直方图，并使用高斯函数拟合灰度直方图，根据高斯函数的方差确定噪声标准差，并设定衰减系数h为噪声标准差的1.2倍，搜索区域半径Ω为整幅图像长度与宽度相比的较大值的10％。

步骤2：根据噪声标准差、衰减系数h和搜索区域半径Ω依次遍历弥散张量磁共振图像数据集分辨率为128行、128列和53层的所有体素，并以遍历到的每个体素为中心，设定半径为Ω的正方搜索区域Q，正方搜索区域Q为13行、13列和13层的正方体。

步骤3：将正方搜索区域Q中的所有体素依次与中心体素进行张量矩阵相似性比较，使用黎曼距离描述张量矩阵相似性的相似性距离d，其公式为：

d(p,q)＝||log(V(p)^-1V(q)||，

其中，V(q)、V(p)分别为比较体素和中心体素的张量矩阵。

步骤4：根据衰减系数h和张量矩阵相似性的相似性距离d获得归一化参数，计算公式如下：

$Z\left(p\right)=\underset{q\∈Q}{\mathrm{\Σ}}\exp (-\frac{d^2(p,q)}{h^2}),$

其中，p为中心体素，q为比较体素，d(p,q)为中心体素与比较体素的张量矩阵相似性的相似性距离。

步骤5：根据归一化参数Z(p)、衰减系数h和张量矩阵相似性的相似性距离d计算该体素的权值，计算公式如下：

$w(p,q)=\frac{1}{Z\left(p\right)}\exp (-\frac{d^2(p,q)}{h^2}).$

步骤6：根据各体素权值w(p,q)，计算被遍历到体素的加权均值张量矩阵I，获得该中心体素的去噪结果，计算公式如下：

$\left(p\right)=\exp \left(\underset{q\∈Q}{\mathrm{\Σ}}w(p,q)\log \left(V\left(q\right)\right)\right),$

其中，I(p)为中心体素去噪后张量矩阵，V(q)为比较体素的张量矩阵。

步骤7：重复步骤1-6，直到遍历完该弥散张量磁共振图像的所有体素。

实施例2：

步骤1：在背景区域，计算背景区域的灰度直方图，并使用高斯函数拟合灰度直方图，根据高斯函数的方差确定噪声标准差，并设定衰减系数h为噪声标准差的1.2倍，搜索区域半径Ω为整幅图像长度与宽度相比的较大值的10％。

步骤2：根据噪声标准差、衰减系数h和搜索区域半径Ω依次遍历弥散张量磁共振图像数据集分辨率为128行、128列和53层的所有体素，并以遍历到的每个体素为中心，设定半径为Ω的正方搜索区域Q，正方搜索区域Q为13行、13列和13层的正方体。

步骤3：将正方搜索区域Q中的所有体素依次与中心体素进行张量矩阵相似性比较，使用黎曼距离描述张量矩阵相似性的相似性距离d，其公式为：

$d(p,q)=\sqrt{\mathrm{Trace}\left[{(l\mathrm{og}\left(V\left(p\right)\right)-\log \left(V\left(q\right)\right))}^2\right]};$

其中，V(q)、V(p)分别为比较体素和中心体素的张量矩阵。

步骤4：根据衰减系数h和张量矩阵相似性的相似性距离d获得归一化参数，计算公式如下：

$Z\left(p\right)=\underset{q\∈Q}{\mathrm{\Σ}}\exp (-\frac{d^2(p,q)}{h^2}),$

其中，p为中心体素，q为比较体素，d(p,q)为中心体素与比较体素的张量矩阵相似性的相似性距离。

步骤5：根据归一化参数Z(p)、衰减系数h和张量矩阵相似性的相似性距离d计算该体素的权值，计算公式如下：

$w(p,q)=\frac{1}{Z\left(p\right)}\exp (-\frac{d^2(p,q)}{h^2}).$

步骤6：根据各体素权值w(p,q)，计算被遍历到体素的加权均值张量矩阵I，获得该中心体素的去噪结果，计算公式如下：

$\left(p\right)=\exp \left(\underset{q\∈Q}{\mathrm{\Σ}}w(p,q)\log \left(V\left(q\right)\right)\right),$

步骤7：重复步骤1-6，直到遍历完该弥散张量磁共振图像的所有体素。

实施例3：

步骤1：在背景区域，计算背景区域的灰度直方图，并使用高斯函数拟合灰度直方图，根据高斯函数的方差确定噪声标准差，并设定衰减系数h为噪声标准差的1.2倍，搜索区域半径Ω为整幅图像长度与宽度相比的较大值的10％。

步骤2：根据噪声标准差、衰减系数h和搜索区域半径Ω依次遍历弥散张量磁共振图像数据集分辨率为128行、128列和53层的所有体素，并以遍历到的每个体素为中心，设定半径为Ω的正方搜索区域Q，正方搜索区域Q为13行、13列和13层的正方体。

步骤3：将正方搜索区域Q中的所有体素依次与中心体素进行张量矩阵相似性比较，使用黎曼距离描述张量矩阵相似性的相似性距离d，其公式为：

$d(p,q)=\underset{i,j=1}{\overset{3}{\mathrm{\Σ}}}{V}_{\mathrm{ij}}\left(p\right)\·V_{\mathrm{ij}}\left(q\right)-(1/3)\mathrm{Trace}\left(V\left(p\right)\right)\·\mathrm{Trace}\left(V\left(q\right)\right);$

其中，V(q)、V(p)分别为比较体素和中心体素的张量矩阵，i、j分别为张量矩阵的行标和列标，Vij(p)、Vij(q)分别为体素p和q张量矩阵中第i行第j列元素。

步骤4：根据衰减系数h和张量矩阵相似性的相似性距离d获得归一化参数，计算公式如下：

$Z\left(p\right)=\underset{q\∈Q}{\mathrm{\Σ}}\exp (-\frac{d^2(p,q)}{h^2}),$

其中，p为中心体素，q为比较体素，d(p,q)为中心体素与比较体素的张量矩阵相似性的相似性距离。

步骤5：根据归一化参数Z(p)、衰减系数h和张量矩阵相似性的相似性距离d计算该体素的权值，计算公式如下：

$w(p,q)=\frac{1}{Z\left(p\right)}\exp (-\frac{d^2(p,q)}{h^2}),$

其中，p为中心体素，q为比较体素。

步骤6：根据各体素权值w(p,q),计算被遍历到体素的加权均值张量矩阵I，获得该中心体素的去噪结果，计算公式如下：

$\left(p\right)=\exp \left(\underset{q\∈Q}{\mathrm{\Σ}}w(p,q)\log \left(V\left(q\right)\right)\right).$

步骤7：重复步骤1-6，直到遍历完该弥散张量磁共振图像的所有体素。

实施例4：

步骤1：在背景区域，计算背景区域的灰度直方图，并使用高斯函数拟合灰度直方图，根据高斯函数的方差确定噪声标准差，并设定衰减系数h为噪声标准差的1.2倍，搜索区域半径Ω为整幅图像长度与宽度相比的较大值的10％。

步骤2：根据噪声标准差、衰减系数h和搜索区域半径Ω依次遍历弥散张量磁共振图像数据集分辨率为128行、128列和53层的所有体素，并以遍历到的每个体素为中心，设定半径为Ω的正方搜索区域Q，正方搜索区域Q为13行、13列和13层的正方体。

步骤3：将正方搜索区域Q中的所有体素依次与中心体素进行张量矩阵相似性比较，使用黎曼距离描述张量矩阵相似性的相似性距离d，其公式为：

$d(p,q)=\sqrt{\mathrm{Trace}\left[{(l\mathrm{og}\left(V\left(p\right)\right)-\log \left(V\left(q\right)\right))}^2\right]};$

其中，V(q)和V(p)分别为比较体素和中心体素的张量矩阵。

步骤4：根据衰减系数h和张量矩阵相似性的相似性距离d获得归一化参数，计算公式如下：

$Z\left(p\right)=\underset{q\∈Q}{\mathrm{\Σ}}\exp (-\frac{d^2(p,q)}{h^2}),$

其中，p为中心体素，q为比较体素，d(p,q)为中心体素与比较体素的张量矩阵相似性的相似性距离。

步骤5：根据归一化参数Z(p)、衰减系数h和张量矩阵相似性的相似性距离d计算该体素的权值，计算公式如下：

$w(p,q)=\frac{1}{Z\left(p\right)}\exp (-\frac{d^2(p,q)}{h^2}),$

其中，p为中心体素，q为比较体素。

步骤6：根据各体素权值w(p,q)，计算被遍历到体素的加权均值张量矩阵I，获得该中心体素的去噪结果，计算公式如下：

$\left(p\right)=\exp \left(\underset{q\∈Q}{\mathrm{\Σ}}w(p,q)\log \left(V\left(q\right)\right)\right).$

步骤7：重复步骤1-6，直到遍历完该弥散张量磁共振图像的所有体素。

本发明已经通过上述实施例进行了说明，但应当理解的是，上述实施例只是用于举例和说明的目的，而非意在将本发明限制于所描述的实施例范围内。此外本领域技术人员可以理解的是，本发明并不局限于上述实施例，根据本发明的教导还可以做出更多种的变型和修改，这些变型和修改均落在本发明所要求保护的范围以内。本发明的保护范围由附属的权利要求书及其等效范围所界定。

Claims (7)

1.一种弥散张量磁共振图像张量域非局部均值去噪方法，其特征在于，如下步骤：

(1)在背景区域，计算背景区域的灰度直方图，并使用高斯函数拟合灰度直方图，根据高斯函数的方差确定噪声标准差，并设定衰减系数h和搜索区域半径Ω；

(2)根据噪声标准差、衰减系数h和搜索区域半径Ω依次遍历弥散张量磁共振图像中的所有体素，并以遍历到的每个体素为中心体素，设定半径为Ω的正方搜索区域Q；

(3)将正方搜索区域Q中的所有体素依次与中心体素进行张量矩阵相似性比较，获得描述张量矩阵相似性的相似性距离d；

(4)根据衰减系数h、正方搜索区域Q和张量矩阵相似性的相似性距离d获得归一化参数，计算公式如下：

$Z\left(p\right)=\underset{q\∈Q}{\mathrm{\Σ}}\exp (-\frac{d^2(p,q)}{h^2}),$

其中，p为中心体素，q为比较体素，d(p,q)为中心体素与比较体素的张量矩阵相似性的相似性距离；

(5)根据归一化参数Z(p)、衰减系数h和张量矩阵相似性的相似性距离d计算各体素的权值，计算公式如下：

$w(p,q)=\frac{1}{Z\left(p\right)}\exp (-\frac{d^2(p,q)}{h^2});$

(6)根据各体素权值w(p,q)，计算被遍历到体素的加权均值张量矩阵I，获得该中心体素的去噪结果，计算公式如下：

$I\left(p\right)=\exp \left(\underset{q\∈Q}{\mathrm{\Σ}}w\right(p,q\left)log\right(V\left(q\right)\left)\right),$

其中，I(p)为中心体素去噪后张量矩阵，V(q)为比较体素的张量矩阵；

(7)重复上述步骤，直到遍历完该弥散张量磁共振图像的所有体素。

2.根据权利要求1所述的一种弥散张量磁共振图像张量域非局部均值去噪方法，其特征在于：所述步骤(1)中，衰减系数h为估计噪声标准差的1-1.5倍。

3.根据权利要求1所述的一种弥散张量磁共振图像张量域非局部均值去噪方法，其特征在于：所述步骤(1)中，搜索区域半径Ω为整幅图像长度与宽度相比的较大值的10％-20％。

4.根据权利要求1所述的一种弥散张量磁共振图像张量域非局部均值去噪方法，其特征在于：所述步骤(3)中，使用黎曼距离描述张量矩阵相似性的相似性距离d，其公式为：

d(p,q)＝||log(V(p)^-1V(q))||；

其中，V(q)、V(p)分别为比较体素和中心体素的张量矩阵。

5.根据权利要求1所述的一种弥散张量磁共振图像张量域非局部均值去噪方法，其特征在于：所述步骤(3)中，使用对数欧式距离描述张量矩阵相似性的相似性距离d，其公式为：

$d(p,q)=\sqrt{Trace\[{\left(log\right(V\left(p\right))-log(V\left(q\right)\left)\right)}^2\]};$

其中，V(q)、V(p)分别为比较体素和中心体素的张量矩阵。

6.根据权利要求1所述的一种弥散张量磁共振图像张量域非局部均值去噪方法，其特征在于：所述步骤(3)中，使用偏张量商距离描述张量矩阵相似性的相似性距离d，其公式为：

$d(p,q)=\underset{i,j=1}{\overset{3}{\mathrm{\Σ}}}{V}_{ij}\left(p\right)\·V_{ij}\left(q\right)-(1/3)Trace\left(V\right(p\left)\right)\·Trace\left(V\right(q\left)\right);$

其中，V(q)、V(p)分别为比较体素和中心体素的张量矩阵，i、j分别为张量矩阵的行标和列标，V_ij(p)、V_ij(q)分别为体素p和q张量矩阵中第i行第j列元素。

7.根据权利要求1所述的一种弥散张量磁共振图像张量域非局部均值去噪方法，其特征在于：所述步骤(3)中，使用张量欧式距离描述张量矩阵相似性的相似性距离d，其公式为：

$d(p,q)=\sqrt{Trace\[{\left(V\right(p)-V(q\left)\right)}^2\]};$

其中，V(q)、V(p)分别为比较体素和中心体素的张量矩阵。