CN111738999B - 基于弥散张量成像的大脑张量模板构建方法 - Google Patents

Abstract

本发明提出了一种基于弥散张量成像的大脑张量模板构建方法，旨在提高大脑张量模板构建的准确度和效率，实现步骤为：对弥散张量图像进行预处理；获取个体平均B0图像集、张量图像集和FA图像集；确定标准模板；对个体FA图像进行迭代配准获取FA模板；获取个体FA图像与FA模板的配准参数；将平均B0图像集和张量图像集空间标准化；获取标准空间下所有被试的平均B0图像和平均张量图像集；获取大脑张量模板。本发明利用FA模板获取个体空间到标准空间的配准参数，作用到张量图像后进行张量重定向，减少了配准信息的丢失，提高了张量模板的准确度，同时所需的配准次数少，提高了模板构建的效率。

Description

基于弥散张量成像的大脑张量模板构建方法

技术领域

本发明属于图像处理技术领域，涉及一种大脑图像模板构建方法，具体涉及一种基于弥散张量成像的大脑张量模板构建方法，可适用于大脑白质的辅助研究。

背景技术

大脑由灰质、白质和脑脊液组成，其中白质构成了中枢神经系统的主要部分。大脑中的白质由数百万条纤维束组成，这些纤维束通过传递动作电位对不同脑区的灰质进行沟通从而实现脑区间的协调运作。因此，当出现白质微结构异常时，可能会导致一些疾病，比如精神分裂症、多发性硬化症等。在对大脑白质疾病进行研究时，通常是对不同个体的弥散张量图像DTI(diffusion tensor imaging)进行组间分析。由于人脑结构具有较大的个体差异性，因此需要一个公共的标准空间，保证不同个体图像的组织扩散信息对应于同一空间位置，再进行后续的组间比较，来检测健康被试与患者之间的大脑白质微小差异。将所有个体的DTI图像进行空间标准化，要求配准的精准度高，要求生成的大脑模板的空间分辨率和对比度高，使之能够分辨出更微小的白质纤维。构建大脑模板的方法很多，具有代表性的是通过结合多个被试个体的大脑图像信息生成平均模板图像，但由于不同被试图像上存在不同程度的伪影和噪声，可能会降低模板的准确性和组平均的有效性，同时配准方法的设计影响着模板生成的效率。

西安电子科技大学在其申请的专利文献“基于白质纤维束的白质微结构特征筛选系统及方法”(申请号：CN201910254191.7，公开号：CN109978872A)中公开了一种大脑张量模板的构建方法，利用图像配准软件，将个体空间下的对照图像组中格式为NIFTI的每张各向异性分数FA(Fractional Anisotropy)图像配准到格式为NIFTI的标准模板，得到配准参数；在标准空间下的所有图像中选取一幅未选图像，将所选图像分别非线性配准到其余未选的每幅图像上，得到所选图像配准到其余的每幅图像上的配准参数，将所有配准参数进行平均，得到所选图像到其余图像上的平均配准参数；将配准参数与平均配准参数融合，得到所选FA图像到标准空间下其余FA图像的变换参数；将所选FA图像的变换参数作用到对应的DTI图像上，得到标准化后的DTI图像；若已选取完标准空间下的所有图像，则已完成所有各项异性图像的标准化，否则，继续在标准空间下的所有图像中选择未选图像；求所有标准化后的DTI图像中对应位置上体素灰度值的平均值；将所有体素灰度的平均值组成一张平均图像，将该平均图像作为构建好的张量模板。该方法存在的不足之处是：对配准参数进行平均、融合等操作得到的变换参数，可能会丢失一部分图像变换信息，导致利用变换参数对相应的DTI图像进行标准化得到的图像准确度不足，而且，对DTI图像进行标准化后，未进行张量方向校正，导致对标准化后的DTI图像平均生成的张量模板准确度不足，同时，为了得到所选图像到其余图像上的所有配准参数，需要进行被试数量×(被试数量-1)次配准操作，工作量大，耗时长，效率低。

发明内容

本发明的目的在于克服上述现有技术存在的缺陷，提出一种基于弥散张量成像的大脑张量模板构建方法，旨在提高大脑张量模板构建的准确度和效率。

为实现上述目的，本发明采用的技术方案包括如下步骤：

(1)对采集的多个被试个体的弥散张量图像进行预处理：

(1a)利用核磁共振采集M个被试个体的维度为4D的大脑矢状位的弥散张量图像，得到弥散张量图像集DTI，DTI＝{DTI₁,DTI₂,...,DTI_i,...,DTI_M}，其中DTI_i表示第i个被试个体的大小为I×J×K×(n+m)的弥散张量图像，包括n个大小为I×J×K维度为3D的未施加敏感梯度方向的图像B0，以及m个大小为I×J×K维度为3D的施加敏感梯度脉冲方向的扩散加权图像DWI，DTI_i＝{B0₁,B0₂,...,B0_l,...,B0_n；DWI₁,DWI₂,...,DWI_k,...,DWI_m}，M≥1，I和J分别表示每幅二维图像的长度和宽度，I≥1,J≥1，K表示采集层数，K≥1，n≥1,m≥6，B0_l表示第l个未施加敏感梯度方向的图像，DWI_k表示第k个施加敏感梯度脉冲方向的扩散加权图像；

(1b)利用图像处理软件对每个DTI_i进行去脑壳处理，并对每个去脑壳处理的图像进行涡流矫正后再进行头动矫正，得到预处理后的弥散张量图像集DTI′，DTI′＝{DTI₁′,DTI₂′,...,DTI_i′,...,DTI_M′}，其中DTI_i′表示DTI_i的预处理弥散张量图像；

(2)获取未施加敏感梯度方向的平均图像集：

计算DTI_i′中n个未施加敏感梯度方向的图像B0对应位置体素灰度值的平均值，得到未施加敏感梯度方向的平均图像集B0′，B0′＝{B0′₁,B0′₂,...,B0′_i,...,B0′_M}，其中B0_i′表示DTI_i′对应的大小为I×J×K的平均图像；

(3)获取张量图像集和FA图像集：

(3a)利用B0_i′和DTI_i′中的m个施加敏感梯度脉冲方向的扩散加权图像DWI，计算DTI_i′中每个体素v的二阶张量D_v，v∈{v₁,v₂,...,v_s,...,v_I×J×K}，其中v_s表示大小为I×J×K的图像中的第s个体素；

(3b)将DTI_i′中所有体素上的二阶张量组合成格式为NIFTI的大小为I×J×K×6的4D张量图像子集DT_i，并将所有张量图像子集组合成DTI′对应的大小为I×J×K×6×M的5D张量图像集DT，DT＝{DT₁,DT₂,...,DT_i,...,DT_M}，其中DT_i表示DTI_i′对应的张量图像子集，DT_i＝{tensor₁,tensor₂,...,tensor_j,...,tensor₆}，tensor_j表示所有体素上的二阶张量中第j个元素组成的大小为I×J×K的张量图像；

(3c)通过每个体素v的二阶张量D_v，计算DTI_i′中每个体素v的各向异性分数FA^*值，并将所有体素的FA^*值组合成格式为NIFTI的大小为I×J×K的FA_i图像，则DTI′的FA图像表示为FA＝{FA₁,FA₂,...,FA_i,...,FA_M}；

(4)确定标准模板：

将脑成像国际联盟ICBM(International Consortium for Brain Mapping)模板作为标准模板，并读取该标准模板的NIFTI格式的大小为I′×J′×K′的FA图像和大小为I′×J′×K′的0、1二值掩模mask图像，其中I′表示标准模板2D图像的长度，I′≥1,J′表示标准模板2D图像的宽度，J′≥1，K′表示标准模板的层数，K′≥1；

(5)对个体FA图像进行迭代配准：

(5a)设迭代次数为t，最大迭代次数为T，T≥1，配准目标模板为I_t，并令t＝1，令I_t为标准模板的NIFTI格式的FA图像；

(5b)将每个FA_i图像与I_t进行线性配准，并将每个线性配准结果与I_t进行非线性配准，得到非线性配准图像集A，A＝{A₁,A₂,...,A_i,...,A_M}，其中A_i表示第i个大小为I′×J′×K′的非线性配准图像；

(5c)计算非线性配准图像集A中所有非线性配准图像对应位置体素的灰度值的平均值，得到大小为I′×J′×K′的平均图像B，并将B中每个体素的灰度值与mask图像对应位置体素的灰度值相乘，得到大小为I′×J′×K′的中间FA模板C；

(5d)判断t＝T是否成立，若是，则C为最终的FA模板，否则，令t＝t+1，I_t＝C，并执行步骤(5b)；

(6)获取每个FA_i图像与FA模板的非线性配准参数：

将每个大小为I×J×K的FA_i图像与大小为I′×J′×K′的FA模板进行线性配准，并将每个线性配准结果与FA模板进行非线性配准，得到非线性配准参数集warp，warp＝{warp₁,warp₂,...,warp_i,...,warp_M}，其中warp_i表示FA_i与FA模板的非线性配准参数；

(7)获取标准空间下的平均图像集B0′和张量图像集DT：

利用warp_i对B0_i′进行转换，得到标准空间下的平均图像集B0″，B0″＝{B0₁″,B0₂″,...,B0_i″,...,B0_M″}，同时利用warp_i对DT_i进行转换，并对转换后的DT_i进行张量重定向，得到标准空间下的张量图像集DT′，DT′＝{DT₁′,DT₂′,...,DT_i′,...,DT_M′}，其中B0_i″表示B0_i′标准空间下的大小为I′×J′×K′的平均图像，DT_i′表示DT_i′标准空间下的大小为I′×J′×K′×6的4D张量图像集；

(8)获取平均图像meanB0和平均张量图像集meanDT：

(8a)将标准空间下的平均图像集B0″中所有图像对应位置体素的灰度值进行平均，得到大小为I′×J′×K′的平均图像meanB0；

(8b)将标准空间下的张量图像集DT′中所有tensor_j图像对应位置体素的灰度值进行平均，得到平均张量图像集meanDT，其中meanDT＝{tensor₁′,tensor₂′,...,tensor_j′,...,tensor₆′}，tensor_j′表示所有DT_i′的第j个张量图像进行平均得到的大小为I′×J′×K′的平均张量图像；

(9)获取大脑张量模板：

利用平均图像meanB0中每个体素v′的灰度值和平均张量图像集meanDT中每个体素v′的灰度值计算扩散加权图像上每个体素v′的灰度值，将所有体素上的灰度值组成扩散加权图像集DWI′，DWI′＝{DWI₁′,DWI₂′,...,DWI_q′,...,DWI_m′}，并对平均图像meanB0和扩散加权图像集DWI′由上到下进行组合，得到大小为I′×J′×K′×(m+1)的4D大脑张量模板DTI′，DTI′＝{meanB0；DWI₁′,DWI₂′,...,DWI_q′,...,DWI_m′}其中DWI_q′表示由meanB0和meanDT计算得到的第q个大小为I′×J′×K′的扩散加权图像，v′∈{v′₁,v′₂,...,v′_p,...,v′_{I′×J′×K′}}，其中v′_p表示大小为I′×J′×K′的图像中的第p个体素。

本发明与现有的技术相比，具有以下优点：

1、本发明通过弥散张量成像，利用个体空间FA图像与迭代配准生成的FA模板进行配准，得到变换参数，再利用变换参数对张量图像进行转换，并对转换后的张量图像进行方向校正，再通过标准空间下的张量图像和平均B0图像重新计算扩散加权图像，最后将扩散加权图像与平均B0图像组合得到张量模板。避免了现有技术图像对中间配准参数进行平均、融合等操作造成的变换信息丢失，以及未进行张量方向校正带来的图像配准误差。与现有技术相比，有效地提高了大脑张量模板的准确度。

2、本发明利用迭代配准的方法生成FA模板，再通过个体空间的FA图像与FA模板进行配准，得到变换参数，避免了现有技术在被试间进行两两配准得到中间配准参数再平均得到平均配准参数，再将平均配准参数和个体空间与标准空间的配准参数融合，得到变换参数。与现有技术相比，减少了配准次数，有效地提高了配准效率。

附图说明

图1是本发明的实现流程图；

图2是本发明所确定的标准模板的格式为NIFTI的FA图像；

图3是本发明与现有技术所构建的大脑张量模板的精度仿真对比图。

具体实施方式

下面将结合附图和具体实施例，对本发明作进一步的详细描述，需要强调的是，本发明不属于疾病的诊断与治疗方法：

参照图1，本发明包括如下步骤：

步骤1)对采集的多个被试个体的弥散张量图像进行预处理：

步骤1a)在弥散张量成像中，通常使用未施加敏感梯度磁场脉冲加权的图像作为基准，然后将未施加敏感梯度脉冲加权的图像按照敏感梯度脉冲方向弥散，从而得到施加敏感梯度脉冲加权的图像；利用核磁共振采集35个被试个体的维度为4D的大脑矢状位的弥散张量图像，得到弥散张量图像集DTI，DTI＝{DTI₁,DTI₂,...,DTI_i,...,DTI₃₅}，其中DTI_i表示第i个被试个体的大小为256×256×35×35的弥散张量图像，包括5个大小为256×256×35维度为3D的未施加敏感梯度方向的图像B0，以及30个大小为256×256×35维度为3D的施加敏感梯度脉冲方向的扩散加权图像DWI，DTI_i＝{B0₁,B0₂,...,B0_l,...,B0₅；DWI₁,DWI₂,...,DWI_k,...,DWI₃₀}，B0_l表示第l个未施加敏感梯度方向的图像，DWI_k表示第k个施加敏感梯度脉冲方向的扩散加权图像；

步骤1b)由于只关注大脑，因此需要去除脑壳部分的图像信息，利用FSL软件对每个DTI_i进行去脑壳处理，并对每个去脑壳处理的图像进行涡流矫正，矫正在磁共振扫描过程中由于人身体切割磁场产生的电流，可以增强图像质量，由于在磁共振扫描过程中，被试的呼吸、血流脉动等生理因素造成的头部运动在所难免，所以要再进行头动矫正，得到预处理后的弥散张量图像集DTI′，DTI′＝{DTI₁′,DTI₂′,...,DTI_i′,...,DTI₃₅′}，其中DTI_i′表示DTI_i的预处理弥散张量图像。

步骤2)获取未施加敏感梯度方向的平均图像集：

通常为了提高图像的信噪比，会采集多个B0图像，但后续的步骤中只需要1个B0图像，所以要计算DTI_i′中5个未施加敏感梯度方向的图像B0对应位置体素的灰度值的平均值，得到未施加敏感梯度方向的平均图像集B0′，B0′＝{B0₁′,B0₂′,...,B0_i′,...,B0₃′₅}，其中B0_i′表示DTI_i′对应的大小为256×256×35的平均图像。

步骤3)获取张量图像集和FA图像集：

步骤3a)分析DTI数据需要通过拟合张量来计算其他的参数图像，利用B0_i′和DTI_i′中的30个施加敏感梯度脉冲方向的扩散加权图像DWI，计算DTI_i′中每个体素v的二阶张量D_v，首先，假设大脑内的水分子扩散运动满足高斯分布，建立描述DTI_i′的体素v中水分子高斯分布的二阶张量D_v：

其中，D_v的不同元素表示体素v内的水分子在不同方向上的扩散率，然后，通过DWI_k图像上体素v的灰度值(S_k)_v和B0′图像上体素v的灰度值(S₀)_v计算D_v的值：

其中，b是扩散敏感系数，与弥散梯度脉冲强度和持续时间有关，在本实施例中，b＝1000，g_k是扩散敏感梯度施加方向的单位化矢量，是磁共振设备的内部参数，g_k＝[x_k,y_k,z_k]^T，[]^T是转置，x_k、y_k和z_k分别是在与DWI对应的x、y和z方向的扩散敏感梯度，在本实施例中，

v∈{v₁,v₂,...,v_s,...,v_256×256×35}，其中v_s表示大小为256×256×35的图像中的第s个体素；

步骤3b)将DTI_i′中所有体素上的二阶张量组合成格式为NIFTI的大小为256×256×35×6的4D张量图像子集DT_i，并将所有张量图像子集组合成DTI′对应的大小为256×256×35×6×35的5D张量图像集DT，DT＝{DT₁,DT₂,...,DT_i,...,DT₃₅}，其中DT_i表示DTI_i′对应的张量图像子集，DT_i＝{tensor₁,tensor₂,...,tensor_j,...,tensor₆}，tensor_j表示所有体素上的二阶张量中第j个元素组成的大小为256×256×35的张量图像；

步骤3c)因为FA图像可以提供较好的灰质和白质对比，而且体素上的值不会随着坐标系统旋转方向的改变而改变，所以后续使用FA图像获取FA模板；通过每个体素v的二阶张量D_v，计算DTI_i′中每个体素v的各向异性分数FA^*值，按照下式对D_v进行特征分解，得到D_v的三个特征值(λ₁)_v、(λ₂)_v和(λ₃)_v：

其中，(e₁)_v、(e₂)_v和(e₃)_v分别是(λ₁)_v、(λ₂)_v和(λ₃)_v对应的特征向量，再利用D_v的三个特征值(λ₁)_v、(λ₂)_v和(λ₃)_v计算体素v的各向异性分数FA^* _v值：

其中，λ_v＝((λ₁)_v+(λ₂)_v+(λ₃)_v)/3，最后将所有体素的FA^*值组合成格式为NIFTI的大小为256×256×35的FA_i图像，则DTI′的FA图像表示为FA＝{FA₁,FA₂,...,FA_i,...,FA₃₅}。

步骤4)确定标准模板：

将脑成像国际联盟ICBM模板作为标准模板，并读取该标准模板的NIFTI格式的大小为91×109×91的FA图像，参照图2，以及大小为91×109×91的0、1二值掩模mask图像。

步骤5)对个体FA图像进行迭代配准：

步骤5a)设迭代次数为t，最大迭代次数可以根据实际情况而定，在本实施例中，已经验证当最大迭代次数为3时得到的大脑张量模板质量最佳，所以设最大迭代次数为3，配准目标模板为I_t，并令t＝1，令I_t为标准模板的NIFTI格式的FA图像；

步骤5b)将每个FA_i图像与I_t进行线性配准，通过计算FA_i图像中所有体素的坐标仿射变换，设FA_i图像中任一体素的坐标为(x,y,z)，计算公式如下：

其中，x′、y′和z′分别是x、y和z变换后的坐标，W是由3个平移参数、3个旋转参数、3个缩放参数和3个剪切参数周期性排布的维数为4×4的矩阵，并将每个线性配准结果与I_t进行非线性配准，通过计算FA_i′图像中所有体素的坐标非线性变换，计算公式如下：

其中，d_x(x,y,z)、d_y(x,y,z)和d_z(x,y,z)分别是x、y和z维度的位移场，得到非线性配准图像集A，A＝{A₁,A₂,...,A_i,...,A₃₅}，其中A_i表示第i个大小为91×109×91的非线性配准图像；

步骤5c)计算非线性配准图像集A中所有非线性配准图像对应位置体素的灰度值的平均值，得到大小为91×109×91的平均图像B，为了减少大脑边缘噪声在迭代配准过程中的积累和传递，将B中每个体素的灰度值与mask图像对应位置体素的灰度值相乘，得到大小为91×109×91的中间FA模板C；

步骤5d)判断t＝3是否成立，若是，则C为最终的FA模板，否则，令t＝t+1，I_t＝C，并执行步骤5b)。

步骤6)获取每个FA_i图像与FA模板的非线性配准参数：

因为FA模板包含了个体信息和标准模板信息，所以认为个体空间的FA_i图像与FA模板的相近程度比个体空间的FA_i图像与标准模板的FA图像的更高；将每个大小为256×256×35的FA_i图像与大小为91×109×91的FA模板进行线性配准，并将每个线性配准结果与FA模板进行非线性配准，得到非线性配准参数集warp，warp＝{warp₁,warp₂,...,warp_i,...,warp₃₅}，其中warp_i表示FA_i与FA模板的非线性配准参数。

步骤7)获取标准空间下的平均图像集B0′和张量图像集DT：

由于张量图像是有方向信息的，而配准过程中只是将体素的解剖位置进行对齐，体素中的矢量信息并没有随着配准而发生改变，因此需要对图像的方向信息进行校正，所以不能直接利用变换参数对个体空间的DTI′进行标准化，而是先对张量图像进行标准化再方向校正；利用warp_i对B0_i′进行转换，得到标准空间下的平均图像集B0″，B0″＝{B0₁″,B0₂″,...,B0_i″,...,B0₃₅″}，同时利用warp_i对DT_i进行转换，得到转换后的DT_i，利用Camino软件对转换后的DT_i进行张量重定向，得到标准空间下的张量图像集DT′，DT′＝{DT₁′,DT₂′,...,DT_i′,...,DT₃₅′}，其中B0_i″表示B0_i′标准空间下的大小为91×109×91的平均图像，DT_i′表示DT_i′标准空间下的大小为91×109×91×6的4D张量图像集。

步骤8)获取平均图像meanB0和平均张量图像集meanDT：

步骤8a)融合所有被试的图像信息，为构建模板做准备，将标准空间下的平均图像集B0″中所有图像对应位置体素的灰度值进行平均，得到大小为91×109×91的平均图像meanB0；

步骤8b)将标准空间下的张量图像集DT′中所有tensor_j图像对应位置体素的灰度值进行平均，得到平均张量图像集meanDT，其中meanDT＝{tensor₁′,tensor₂′,...,tensor_j′,...,tensor₆′}，tensor_j′表示所有DT_i′的第j个张量图像进行平均得到的大小为91×109×91的平均张量图像。

步骤9)获取大脑张量模板：

通过标准空间下的平均图像meanB0和标准空间下已进行方向校正处理的平均张量图像集meanDT，可以计算出标准空间下的融合个体信息和标准模板信息的张量模板；利用平均图像meanB0中每个体素v′的灰度值(S₀′)_v′和平均张量图像集meanDT中每个体素v′的灰度值D_v′计算扩散加权图像上每个体素v′的灰度值(S_q′)_v′：

将所有体素上的灰度值S_q′组成扩散加权图像集DWI′，DWI′＝{DWI₁′,DWI₂′,...,DWI_q′,...,DWI₃₀′}，并对平均图像meanB0和扩散加权图像集DWI′由上到下进行组合，得到大小为91×109×91×31的4D大脑张量模板DTI′，DTI′＝{meanB0；DWI₁′,DWI₂′,...,DWI_q′,...,DWI₃₀′}其中DWI_q′表示由meanB0和meanDT计算得到的第q个大小为91×109×91的扩散加权图像，v′∈{v′₁,v′₂,...,v′_p,...,v′_91×109×91}，其中v′_p表示大小为91×109×91的图像中的第p个体素。

以下结合仿真实验，对本发明的技术效果作一说明。

1.仿真条件和内容：

利用FSL软件对本实施例中的35个被试个体的弥散张量图像集DTI＝{DTI₁,DTI₂,...,DTI_i,...,DTI₃₅}和FA图像集FA＝{FA₁,FA₂,...,FA_i,...,FA₃₅}，利用西安电子科技大学在其申请的专利文献“基于白质纤维束的白质微结构特征筛选系统及方法”(申请号：CN201910254191.7，公开号：CN109978872A)中涉及的大脑张量模板构建方法进行仿真，并选取与本发明相同的标准模板。

首先将35个FA图像分别线性配准再非线性配准到ICBM模板的NIFTI格式的FA图像，得到配准参数a；然后将每个标准空间下FA图像与除自身之外的FA图像进行配准，每个被试得到34个非线性配准参数，将34个非线性配准参数平均，得到的每个被试到除自身之外的其他所有FA图像的平均配准参数b；再将配准参数a与平均配准参数b融合，得到每个被试的个体FA图像到标准空间下其余FA图像的变换参数c；接着将所选FA图像的变换参数c作用到对应的DTI图像上，得到标准化后的DTI图像；再求所有标准化后的DTI图像中对应位置上体素灰度值的平均值；最后将所有体素灰度的平均值组成一张平均图像，得到张量模板。

2.仿真结果分析：

参照图3，图3(a)表示本实施例中使用本发明实现的大脑张量模板，图3(b)表示本实施例中使用现有技术仿真实现的大脑张量模板，从图中可以看出，与现有技术相比，本发明实现的大脑张量模板轮廓和形状与标准模板更接近，如图3(a)和图3(b)中的白色方框所示，在白质与灰质交界部位的清晰度和对比度更高，准确度有所提高，如图3(a)和图3(b)中的白色箭头所示。

Claims (4)

1.一种基于弥散张量成像的大脑张量模板构建方法，其特征在于，包括如下步骤：

(1)对采集的多个被试个体的弥散张量图像进行预处理：

(1a)利用核磁共振采集M个被试个体的维度为4D的大脑矢状位的弥散张量图像，得到弥散张量图像集DTI，DTI＝{DTI₁,DTI₂,...,DTI_i,...,DTI_M}，其中DTI_i表示第i个被试个体的大小为I×J×K×(n+m)的弥散张量图像，包括n个大小为I×J×K维度为3D的未施加敏感梯度方向的图像B0，以及m个大小为I×J×K维度为3D的施加敏感梯度脉冲方向的扩散加权图像DWI，DTI_i＝{B0₁,B0₂,...,B0_l,...,B0_n；DWI₁,DWI₂,...,DWI_k,...,DWI_m}，M≥1，I和J分别表示每幅二维图像的长度和宽度，I≥1,J≥1，K表示采集层数，K≥1，n≥1,m≥6，B0_l表示第l个未施加敏感梯度方向的图像，DWI_k表示第k个施加敏感梯度脉冲方向的扩散加权图像；

(1b)利用图像处理软件对每个DTI_i进行去脑壳处理，并对每个去脑壳处理的图像进行涡流矫正后再进行头动矫正，得到预处理后的弥散张量图像集DTI′，DTI′＝{DTI₁′,DTI₂′,...,DTI_i′,...,DTI_M′}，其中DTI_i′表示DTI_i的预处理弥散张量图像；

(2)获取未施加敏感梯度方向的平均图像集：

计算DTI_i′中n个未施加敏感梯度方向的图像B0对应位置体素灰度值的平均值，得到未施加敏感梯度方向的平均图像集B0′，B0′＝{B0′₁,B0′₂,...,B0′_i,...,B0′_M}，其中B0_i′表示DTI_i′对应的大小为I×J×K的平均图像；

(3)获取张量图像集和FA图像集：

(3a)利用B0_i′和DTI_i′中的m个施加敏感梯度脉冲方向的扩散加权图像DWI，计算DTI_i′中每个体素v的二阶张量D_v，v∈{v₁,v₂,...,v_s,...,v_I×J×K}，其中v_s表示大小为I×J×K的图像中的第s个体素；

(3b)将DTI_i′中所有体素上的二阶张量组合成格式为NIFTI的大小为I×J×K×6的4D张量图像子集DT_i，并将所有张量图像子集组合成DTI′对应的大小为I×J×K×6×M的5D张量图像集DT，DT＝{DT₁,DT₂,...,DT_i,...,DT_M}，其中DT_i表示DTI_i′对应的张量图像子集，DT_i＝{tensor₁,tensor₂,...,tensor_j,...,tensor₆}，tensor_j表示所有体素上的二阶张量中第j个元素组成的大小为I×J×K的张量图像；

(3c)通过每个体素v的二阶张量D_v，计算DTI_i′中每个体素v的各向异性分数FA^*值，并将所有体素的FA^*值组合成格式为NIFTI的大小为I×J×K的FA_i图像，则DTI′的FA图像表示为FA＝{FA₁,FA₂,...,FA_i,...,FA_M}；

(4)确定标准模板：

将脑成像国际联盟ICBM模板作为标准模板，并读取该标准模板的NIFTI格式的大小为I′×J′×K′的FA图像和大小为I′×J′×K′的0、1二值掩模mask图像，其中I′表示标准模板2D图像的长度，I′≥1,J′表示标准模板2D图像的宽度，J′≥1，K′表示标准模板的层数，K′≥1；

(5)对个体FA图像进行迭代配准：

(5a)设迭代次数为t，最大迭代次数为T，T≥1，配准目标模板为I_t，并令t＝1，令I_t为标准模板的NIFTI格式的FA图像；

(5b)将每个FA_i图像与I_t进行线性配准，并将每个线性配准结果与I_t进行非线性配准，得到非线性配准图像集A，A＝{A₁,A₂,...,A_i,...,A_M}，其中A_i表示第i个大小为I′×J′×K′的非线性配准图像；

(5c)计算非线性配准图像集A中所有非线性配准图像对应位置体素的灰度值的平均值，得到大小为I′×J′×K′的平均图像B，并将B中每个体素的灰度值与mask图像对应位置体素的灰度值相乘，得到大小为I′×J′×K′的中间FA模板C；

(5d)判断t＝T是否成立，若是，则C为最终的FA模板，否则，令t＝t+1，I_t＝C，并执行步骤(5b)；

(6)获取每个FA_i图像与FA模板的非线性配准参数：

将每个大小为I×J×K的FA_i图像与大小为I′×J′×K′的FA模板进行线性配准，并将每个线性配准结果与FA模板进行非线性配准，得到非线性配准参数集warp，warp＝{warp₁,warp₂,...,warp_i,...,warp_M}，其中warp_i表示FA_i与FA模板的非线性配准参数；

(7)获取标准空间下的平均图像集B0′和张量图像集DT：

利用warp_i对B0_i′进行转换，得到标准空间下的平均图像集B0″，B0″＝{B0₁″,B0₂″,...,B0_i″,...,B0_M″}，同时利用warp_i对DT_i进行转换，并对转换后的DT_i进行张量重定向，得到标准空间下的张量图像集DT′，DT′＝{DT₁′,DT₂′,...,DT_i′,...,DT_M′}，其中B0_i″表示B0_i′标准空间下的大小为I′×J′×K′的平均图像，DT_i′表示DT_i′标准空间下的大小为I′×J′×K′×6的4D张量图像集；

(8)获取平均图像meanB0和平均张量图像集meanDT：

(8a)将标准空间下的平均图像集B0″中所有图像对应位置体素的灰度值进行平均，得到大小为I′×J′×K′的平均图像meanB0；

(8b)将标准空间下的张量图像集DT′中所有tensor_j图像对应位置体素的灰度值进行平均，得到平均张量图像集meanDT，其中meanDT＝{tensor₁′,tensor₂′,...,tensor_j′,...,tensor₆′}，tensor_j′表示所有DT_i′的第j个张量图像进行平均得到的大小为I′×J′×K′的平均张量图像；

(9)获取大脑张量模板：

利用平均图像meanB0中每个体素v′的灰度值和平均张量图像集meanDT中每个体素v′的灰度值计算扩散加权图像上每个体素v′的灰度值，将所有体素上的灰度值组成扩散加权图像集DWI′，DWI′＝{DWI₁′,DWI₂′,...,DWI_q′,...,DWI_m′}，并对平均图像meanB0和扩散加权图像集DWI′由上到下进行组合，得到大小为I′×J′×K′×(m+1)的4D大脑张量模板DTI′，DTI′＝{meanB0；DWI₁′,DWI₂′,...,DWI_q′,...,DWI_m′}其中DWI_q′表示由meanB0和meanDT计算得到的第q个大小为I′×J′×K′的扩散加权图像，v′∈{v′₁,v′₂,...,v′_p,...,v′_{I′×J′×K′}}，其中v′_p表示大小为I′×J′×K′的图像中的第p个体素。

2.根据权利要求1所述的基于弥散张量成像的大脑张量模板构建方法，其特征在于，步骤(3a)中所述的计算DTI_i′中每个体素v的二阶张量D_v，实现步骤为：

(3a1)通过B0_i′和DTI_i′中的m个施加敏感梯度脉冲方向的扩散加权图像DWI，计算DTI_i′中每个体素v的二阶张量D_v，v∈{v₁,v₂,...,v_s,...,v_I×J×K}，其中v_s表示大小为I×J×K的图像中的第s个体素；

(3a2)假设大脑内的水分子扩散运动满足高斯分布，建立描述DTI_i′中体素v中的水分子高斯分布的二阶张量D_v：

其中，D_v的不同元素表示体素v内的水分子在不同方向上的扩散率；

(3a3)通过DWI_k图像上体素v的灰度值(S_k)_v和B0_i′图像上体素v的灰度值(S₀)_v计算D_v的值：

其中，b是扩散敏感系数，g_k是扩散敏感梯度施加方向的单位化矢量，g_k＝[x_k,y_k,z_k]^T，[]^T是转置，x_k、y_k和z_k分别是在与DWI对应的x、y和z方向的扩散敏感梯度。

3.根据权利要求1所述的基于弥散张量成像的大脑张量模板构建方法，其特征在于，步骤(3c)中所述的计算DTI_i′中每个体素的各向异性分数FA^*值，实现步骤为：

(3c1)按照下式对体素v的D_v进行特征分解，得到D_v的三个特征值(λ₁)_v、(λ₂)_v和(λ₃)_v：

其中，(e₁)_v、(e₂)_v和(e₃)_v分别是(λ₁)_v、(λ₂)_v和(λ₃)_v对应的特征向量；

(3c2)利用D_v的三个特征值(λ₁)_v、(λ₂)_v和(λ₃)_v计算体素v的各向异性分数FA^* _v值：

其中，λ_v＝((λ₁)_v+(λ₂)_v+(λ₃)_v)/3。

4.根据权利要求1所述的基于弥散张量成像的大脑张量模板构建方法，其特征在于，步骤(5b)中所述的将每个FA_i图像与I_t进行线性配准，并将每个线性配准结果与I_t进行非线性配准，实现步骤为：

(5b1)计算FA_i图像中所有体素的坐标仿射变换，设FA_i图像中任一体素的坐标为(x,y,z)，计算公式如下：

其中，x′、y′和z′分别是x、y和z变换后的坐标，W是由3个平移参数、3个旋转参数、3个缩放参数和3个剪切参数周期性排布的维数为4×4的矩阵；

(5b2)计算FA_i′图像中所有体素的坐标非线性变换，计算公式如下：

其中，d_x(x,y,z)、d_y(x,y,z)和d_z(x,y,z)分别是x、y和z维度的位移场。

Images