技术领域
[0001] 本
发明涉及磁共振扩散加权图像平滑方法，属于医学
图像处理领域。
背景技术
[0002] 磁共振扩散成像是目前唯一能够在活体上测量组织内
水分子扩散运动与成像的无创方法，它通过测量和量化组织中水分子的扩散信息来探测组织的微观结构。水分子沿不同方向的扩散信息包含在一组不同扩散加权梯度方向的扩散加权图像(Diffusion Weighted Image,DWI)中，通过对扩散函数进行建模可以解析出每个
体素内的
纤维束结构信息(主要是纤维的走行方向)。根据各体素内纤维束的走行方向，利用纤维束追踪技术可重建出组织纤维束的三维结构，可用于医疗诊断及相关研究等。
[0003] 目前针对复杂组织纤维结构的解析方法主要有：扩散谱、Q-Ball、多张量、高阶张量等方法，这些方法一般要求成像时所采用的扩散加权梯度方向数量M较多(≥100)，然而由于
硬件和采集时间的限制，在临床中磁共振扩散成像
数据采集时一般M≤30，因此这些方法难以应用于临床。
[0004] 而目前能在M≤30的情况下解析纤维结构的有效方法仅有几种，其中较为有效的是使用受限压缩传感(Constrained Compress Sensing，CCS)技术。该方法基于经典的多张量模型，即将磁共振扩散加权
信号建模为如下多个高斯函数加权求和的形式：
[0005]
[0006] 其中，Sk是扩散梯度方向为gk时的磁共振扩散信号，S0是未加扩散梯度场时的磁共振信号，b是扩散加权敏感因子，gk是第k个扩散梯度方向，k＝1,2,...,M，M为磁共振扩散成像时所使用的扩散梯度方向的总数，fj≥0是对应于二阶
扩散张量Dj的加权系数，N是高斯函数的个数，ηk为噪声，每个二阶扩散张量Dj对应于一个主特征方向vjT,该方向即为纤维的走向。首先构造一组可能的纤维走向集合{vj,j＝1,2,...,N}，N＞＞M，根据数据的特性设定扩散张量的特征值，根据特征值和方向即可构造出一组扩散张量{Dj,j＝1,2,...,N}。令：
[0007] y＝[S1,S2,…,SM]T/S0
[0008]T
[0009] f＝[f1,f2,…,fN]T
[0010] η＝[η1,η2,…,ηN]
[0011] 于是根据式(1)有：
[0012] y＝Bf+η
[0013] 然后采用受限压缩传感技术求解如下优化问题：
[0014]
[0015] 得到f的估计值 其中β为正则化系数。 中的每个元素分别表征了其所对应的纤维走行方向在该体素中所占的权重。在实际情况中， 中的大部分元素均为零，表示该体素的纤维结构不存在这些元素所对应的方向。将集合T
中 为0的元素去除后的集合记为{(αi,di),i＝1,2,…,n}，其
中n为 中非零元素的个数，它表征了该体素内的纤维存在多少个不同的走行方向。因此T
集合{(αi,di),i＝1,2,…,n}即可表征一个体素内的纤维结构信息。
[0016] 由于磁共振扩散加权图像噪声较强，而上述方法对噪声比较敏感，尤其是当磁共振扩散成像采集数据所使用的扩散加权梯度方向数量M较少时，计算得到的纤维结构信息准确度较低。
发明内容
[0017] 本发明是为了解决现有方法噪声抑制差，从而使得获得的每个体素的纤维结构信息精确度低的问题，而提供了基于多张量的磁共振扩散加权图像结构自适应平滑方法。
[0018] 基于多张量的磁共振扩散加权图像结构自适应平滑方法，实现该方法的步骤如下：
[0019] 步骤一：选取参数λ、参数a和总
迭代次数kN的值，设定初始邻域半径h(0)，设定当(1) (0)前迭代次数k＝1，计算h ＝ah ，λ的取值范围为0.5～5，a的取值范围为1
kN的取值范围为5～10，h 的取值范围为0.5～1.5；
[0020] 步骤二：计算每个体素p的初始纤维结构信息并设定参数 其中 表示该体素存在
方向为 的纤维，其所占比重为 表示该体素内的纤维存在多少个不同的走行方向，参数 的取值范围为1～2；
[0021] 步骤三：对于每个体素p，记q为以p为中心h(k)为半径的邻域U(p)内的体素，根据每个体素p的纤维结构信息计算其邻域半径内的所有体素对该体素p的权重，并根据获得的所有权重值对磁共振扩散加权图像进行加权平滑处理，然后重新计算每个体素的纤维结构信息；
[0022] 步骤四：判断是否满足迭代终止条件，如果不满足，则扩大邻域半径并继续进行步(k) (k-1)骤三，即如果k≤kN，则另k＝k+1，h ＝ah ，返回步骤三，其中k为迭代次数，kN表示迭代终止次数；否则，完成基于多张量的磁共振扩散加权图像结构自适应平滑。
[0023] 本发明根据邻域体素与中心体素的距离和纤维结构信息的差别，自适应的确定邻域体素的权重，在迭代过程中逐渐扩大邻域半径的同时，对于距离中心体素较近且与中心体素结构信息差别较小的邻域体素赋予较大权重，而对于距离中心体素较远且与中心体素结构信息差别较大的邻域体素赋予较小权重，进而对中心体素进行加权平滑，这样能够很好的对磁共振扩散加权图像进行平滑同时又能够避免图像边缘模糊。因此本发明在每次迭代平滑后，重新计算纤维结构信息，能够达到更精确的从磁共振扩散加权数据中提取纤维结构信息的目的。
[0024] 本发明将自适应平滑方法与纤维结构解析方法(如受限压缩传感方法)结合使用，能够很好的抑制图像噪声，从而提高纤维结构信息的计算
精度，精确度提高了约30％。
附图说明
[0025] 图1为具体实施方式八中仿真磁共振扩散加权图像数据的实际纤维结构示意图；
[0026] 图2为具体实施方式八中数据经第2次迭代平滑后采用受限压缩传感方法计算得到的纤维结构示意图；
[0027] 图3为具体实施方式八中数据经第4次迭代平滑后采用受限压缩传感方法计算得到的纤维结构示意图；
[0028] 图4为具体实施方式八中数据经第6次迭代平滑后采用受限压缩传感方法计算得到的纤维结构示意图；
[0029] 图5为具体实施方式八中数据经第8次迭代平滑后采用受限压缩传感方法计算得到的纤维结构示意图；
[0030] 图6为具体实施方式八中数据未经平滑直接采用受限压缩传感方法计算得到的纤维结构示意图；
[0031] 图7为具体实施方式八中数据经高斯平滑后采用受限压缩传感方法计算得到的纤维结构示意图；
[0032] 图8为具体实施方式九中仿真实体磁共振扩散加权图像数据的实际纤维结构示意图；
[0033] 图9为具体实施方式九中数据经第2次迭代平滑后采用受限压缩传感方法计算得到的纤维结构示意图；
[0034] 图10为具体实施方式九中数据经第4次迭代平滑后采用受限压缩传感方法计算得到的纤维结构示意图；
[0035] 图11为具体实施方式九中数据经第6次迭代平滑后采用受限压缩传感方法计算得到的纤维结构示意图；
[0036] 图12为具体实施方式九中数据未经平滑直接采用受限压缩传感方法计算得到的纤维结构示意图；
[0037] 图13为具体实施方式九中数据经高斯平滑后采用受限压缩传感方法计算得到的纤维结构示意图；
[0038] 图14为具体实施方式一所述的基于多张量的磁共振扩散加权图像结构自适应平滑方法的
流程图。
具体实施方式
[0039] 具体实施方式一：参见图14说明本实施方式，本实施方式所述的基于多张量的磁共振扩散加权图像结构自适应平滑方法，实现该方法的步骤如下：
[0040] 步骤一：选取参数λ、参数a和总迭代次数kN的值，设定初始邻域半径h(0)，设定当(1) (0)前迭代次数k＝1，计算h ＝ah ，λ的取值范围为0.5～5，a的取值范围为1
kN的取值范围为5～10，h 的取值范围为0.5～1.5；
[0041] 步骤二：计算每个体素p的初始纤维结构信息并设定参数 其中 表示该体素存在
方向为 的纤维，其所占比重为 表示该体素内的纤维存在多少个不同的走行方向，参数 的取值范围为1～2；
[0042] 步骤三：对于每个体素p，记q为以p为中心h(k)为半径的邻域U(p)内的体素，根据每个体素p的纤维结构信息计算其邻域半径内的所有体素对该体素p的权重，并根据获得的所有权重值对磁共振扩散加权图像进行加权平滑处理，然后重新计算每个体素的纤维结构信息；
[0043] 步骤四：判断是否满足迭代终止条件，如果不满足，则扩大邻域半径并继续进行步(k) (k-1)骤三，即如果k≤kN，则另k＝k+1，h ＝ah ，返回步骤三，其中k为迭代次数，kN表示迭代终止次数；否则，完成基于多张量的磁共振扩散加权图像结构自适应平滑。
[0044] 具体实施方式二：本实施方式是对具体实施方式一的进一步限定，所述步骤二所述计算每个体素p的初始纤维结构信息 是采用受限压缩传感方法计算得到的。
[0045] 本实施方式所述的每个体素p的初始纤维结构信息 不仅可以采用受限压缩传感方法，也可以采用其它方法。
[0046] 具体实施方式三：本实施方式是对具体实施方式一的进一步限定，所述步骤三所述的根据每个体素p的纤维结构信息计算其邻域半径内的所有体素对该体素的权重，从而对磁共振扩散加权图像进行加权平滑，平滑后重新计算每个体素的纤维结构信息的方法为：
[0047] 计算 其中ρ(p,q)为体素p到体素q的欧氏距离；
[0048] 计算 其中函数 表示体素p与体素q之间纤维结构信息的差别；
[0049] 计算权重 其中核函数Kloc(·)和Kst(·)为两个定义域为[0,∞)的单调递减函数，且满足Kloc(0)＝Kst(0)＝1，当x≥1时Kloc(x)＝Kst(x)＝
0；
[0050] 计算邻域U(p)内所有体素对p的权重之和
[0051] 计算 及 i＝1,2,…,M,其中 表示未加扩散梯度场时体素p的磁共振信号经第k次迭代平滑计算的结果， 表示未加扩散梯度场时体素q的磁共振信号经第k-1次迭代平滑计算的结果， 表示扩散梯度方向为gi时体素p的磁共振信号经第k次迭代平滑计算的结果， 表示扩散梯度方向为gi时体素q的磁共振信号经第k-1次迭代平滑计算的结果，M为所加不同扩散梯度方向的个数。
[0052] 根据迭代平滑计算后的 和 重新计算体素p的纤维结构信息
[0053] 具体实施方式四：本实施方式是对具体实施方式三的进一步限定，所述步骤三所述体素p与体素q之间纤维结构信息的差别 通过如下步骤计算：
[0054] 步骤a、分别对经过第k-1次迭代每个体素p所占比重 和经过第k-1次迭代每个体素q所占比重 进行归一化，即令表示经过第k-1次迭代每个体素p内的纤维存在多少个不同的走行方向、 表示经过第k-1次迭代每个体素q内的纤维存在多少个不同的走行方向；
[0055] 步骤b、通过求解如下线性规划问题得到 的值：
[0056]
[0057]
[0058]
[0059] xij≥0
[0060]
[0061] 式中 xij表示线性规划系数矩阵x中的元素， 表示经过第k-1次迭代的每个体素p的纤维结构信息、 表示经过第k-1次迭代的每个体素q的纤维结构信息、dpi表示体素p的第i个纤维走行方向、dqj表示体素q的第j个纤维走行方向、αqj表示体素q的第j个纤维走行方向的纤维占体素内纤维总量的比重、αpi表示体素p的第i个纤维走行方向的纤维占体素内纤维总量的比重。
[0062] 具体实施方式五：本实施方式是对具体实施方式三的进一步限定，所述步骤三所述根据迭代平滑计算后的 和 重新计算体素p的纤维结构信息是采用受限压缩传感方法。
[0063] 本实施方式所述的迭代平滑计算后的 和 重新计算体素p的纤维结构信息除了采用受限压缩传感方法也可采用其它方法。
[0064] 具体实施方式六：本实施方式是对具体实施方式三的进一步限定，所述步骤三所述核函数Kloc(·)和Kst(·)选取定义域为[0,∞)的单调递减，且满足Kloc(0)＝Kst(0)＝1，当x≥1时Kloc(x)＝Kst(x)＝0的函数。
[0065] 具体实施方式七：本实施方式是对具体实施方式六的进一步限定，所述步骤三所述核函数 Kloc(·)和 Kst(·)分 别 为：
[0066] 具体实施方式八：本实施方式是具体实施方式一所述的基于多张量的磁共振扩散加权图像结构自适应平滑方法的具体
实施例，本实施例以仿真数据为例详述如下：
[0067] 本实施例中仿真数据的生成是采用高斯混合模型法模拟一组
分辨率为16×16的磁共振扩散加权图像数据，该组数据包含30个不同扩散加权梯度方向的磁共振扩散加权图像，扩散敏感因子b值为1000，该数据即包含纤维存在交叉的体素，也包含纤维不存在交叉的体素，其实际纤维结构如图1所示，图中的每个线段表示其所在体素中纤维的一个走行方向。在生成的仿真数据中加入一定强度的莱斯噪声，得到一组
信噪比24的仿真磁共振扩散加权图像数据。使用本发明基于多张量的磁共振扩散加权图像结构自适应平滑方法对该组含噪数据进行处理的步骤如下：
[0068] 步骤一：选取参数λ＝1，参数a＝1.25，总迭代次数kN＝8，初始邻域半径h(0)(1) (0)＝1，设定当前迭代次数k＝1，计算h ＝ah ；
[0069] 步骤二：计算每个体素p的初始纤维结构信息并设定参数 其中 表示该体素存在
方向为 的纤维，其所占比重为 表示该体素内的纤维存在多少个不同的走行方向，参数 且采用受限压缩传感方法计算纤维结构信息；
(k)
[0070] 步骤三：对于每个体素p，记q为以p为中心h 为半径的邻域U(p)内的体素，进行如下计算：
[0071] 计算 其中ρ(p,q)为体素p到体素q的欧氏距离；
[0072] 计算 其中函数 表示体素p与体素q之间纤维结构信息的差别；
[0073] 计算权重 其中核函数Kloc(·)和Kst(·)为两个定义域为[0,∞)的单调递减函数，且满足Kloc(0)＝Kst(0)＝1，当x≥0时Kloc(x)＝Kst(x)＝
0；在本实施例中取
[0074] 计算邻域U(p)内所有体素对p的权重之和
[0075] 计算 及 i＝1,2,…,M,其中 表示未加扩散梯度场时体素p的磁共振信号经第k次迭代平滑计算的结果， 表示扩散梯度方向为gi时体素p的磁共振信号经第k次迭代平滑计算的结果，M为所加不同扩散梯度方向的个数，其中M＝30；
[0076] 采用受限压缩传感方法根据迭代平滑计算后的 和 重新计算体素p的纤维结构信息(k) (k-1)
[0077] 步骤四：如果k≤kN，则另k＝k+1，h ＝ah ，返回步骤三；否则结束。
[0078] 其中步骤三所述体素p与体素q之间纤维结构信息的差别 通过如下步骤计算：
[0079] a、分别对 和 进行归一化，即令
[0080] b、通过求解如下线性规划问题得到 的值
[0081]
[0082]
[0083]
[0084] xij≥0
[0085]
[0086] 式中
[0087] 本实施例中经第2、4、6、8次迭代后重新计算得到的纤维结构信息分别如图2、3、4、5所示。不对数据进行平滑直接采用受限压缩传感
算法得到的纤维结构信息如图6所示。
采用高斯平滑方法对磁共振扩散加权数据平滑后再使用受限压缩传感算法得到的纤维结构信息如图7所示。对比可以看出经本发明方法对数据进行平滑后重新计算得到的纤维结构信息的准确性较高。
[0088] 具体实施方式九：本实施方式是具体实施方式一所述的基于多张量的磁共振扩散加权图像结构自适应平滑方法的具体实施例，本实施例以一组网上公开的仿真实体磁共振扩散加权图像数据中的部分数据为例详述如下：
[0089] 仿真实体磁共振扩散加权图像数据分辨率为128×128，包含64个不同扩散加权梯度方向的磁共振扩散加权图像，扩散敏感因子b值为1500。该组数据的实际纤维结构如图7所示，图中箭头所指的方向为纤维的走行方向。使用本发明基于多张量的磁共振扩散加权图像结构自适应平滑方法对该组仿真实体数据进行处理的具体步骤如下：
[0090] 步骤一：在本实施例中，选取参数λ＝3，参数a＝1.25，总迭代次数kN＝6，初始(0) (1) (0)邻域半径h ＝1，设定当前迭代次数k＝1，计算h ＝ah ，
[0091] 步骤二：计算每个体素p的初始纤维结构信息并设定参数 其中 表示该体素存在
方向为 的纤维，其所占比重为 表示该体素内的纤维存在多少个不同的走行方向，参数 且采用受限压缩传感方法计算纤维结构信息；
[0092] 步骤三：对于每个体素p，记q为以p为中心h(k)为半径的邻域U(p)内的体素，进行如下计算：
[0093] 计算 其中ρ(p,q)为体素p到体素q的欧氏距离；
[0094] 计算 其中函数 表示体素p与体素q之间纤维结构信息的差别；
[0095] 计算权重 其中核函数Kloc(·)和Kst(·)为两个定义域为[0,∞)的单调递减函数，且满足Kloc(0)＝Kst(0)＝1，当x≥0时Kloc(x)＝Kst(x)＝
0；在本实施例中取
[0096] 计算邻域U(p)内所有体素对p的权重之和
[0097] 计算 及 i＝1，2，…，M，其中表示未加扩散梯度场时体素p的磁共振信号经第k次迭代平滑计算的结果， 表示扩散梯度方向为gi时体素p的磁共振信号经第k次迭代平滑计算的结果，M为所加不同扩散梯度方向的个数，其中M＝30；
[0098] 采用受限压缩传感方法根据迭代平滑计算后的 和 重新计算体素p的纤维结构信息
[0099] 步骤四：如果k≤kN，则另k＝k+1，h(k)＝ah(k-1)，返回步骤三；否则结束；
[0100] 其中步骤三所述体素p与体素q之间纤维结构信息的差别 通过如下步骤计算：
[0101] a、分别对 和 进行归一化，即令
[0102] b、通过求解如下线性规划问题得到 的值
[0103]
[0104]
[0105]
[0106] xij≥0
[0107]
[0108] 式中
[0109] 本实施例中经第2、4、6次迭代后重新计算得到的纤维结构信息分别如图9、10、11所示。不对数据进行平滑直接采用受限压缩传感算法得到的纤维结构信息如图12所示。采用高斯平滑方法对磁共振扩散加权数据平滑后再使用受限压缩传感算法得到的纤维结构信息如图13所示。对比可以看出经本发明方法对数据进行平滑后重新计算得到的纤维结构信息的准确性较高。
