CN105913465A

CN105913465A - Overall sparsity regularization model-based fiber reconstructing method

Info

Publication number: CN105913465A
Application number: CN201610216710.7A
Authority: CN
Inventors: 冯远静; 何建忠; 吴烨; 徐田田; 张军; 黄奕奇; 周思琪
Original assignee: Zhejiang University of Technology ZJUT
Current assignee: Zhejiang University of Technology ZJUT
Priority date: 2016-04-07
Filing date: 2016-04-07
Publication date: 2016-08-31

Abstract

一种基于全局稀疏正则化模型的纤维重构方法，包括如下步骤：1)建立基于字典基重构方法的局部稀疏模型；2)建立全局模型；3)全局优化算法的价值函数。本发明提供一种准确性较高的基于全局稀疏正则化模型的纤维重构方法。A fiber reconstruction method based on a global sparse regularization model, comprising the following steps: 1) establishing a local sparse model based on a dictionary-based reconstruction method; 2) establishing a global model; 3) a value function of a global optimization algorithm. The invention provides a fiber reconstruction method based on a global sparse regularization model with high accuracy.

Description

A Fiber Reconstruction Method Based on Global Sparse Regularization Model

技术领域technical field

本发明涉及计算机图形学下的医学成像、神经解剖学领域，尤其是一种纤维重构方法。The invention relates to the fields of medical imaging and neuroanatomy under computer graphics, especially a fiber reconstruction method.

背景技术Background technique

为了建立准确的纤维取向估算办法，解决神经结构复合纤维的地区，就必须使用新型的高角分辨率扩散成像(HARDI)技术；球面卷积方法是一种有效的纤维取向(FOD)估计方法，该模型为纤维取向分布响应函数的卷积；然而，这些方法大多以基于体素的方式估计FOD，而忽略了与领域体素的空间先验信息；领域体素的FOD估计误差，可能会由于纤维跟踪的累计误差导致重构纤维偏离实际纤维；为了克服限制，空间HARDI或全局重构方法已经被用来将空间约束作为空间信号的参数重构，并尝试通过判定它们的结构能够更好地描述测量数据的同时重构纤维束。全局重构的目的是在整体中提供纤维架构一致性视图；然而，这些方法只能用连接纤维片段之间的线性平滑约束。In order to establish an accurate method for fiber orientation estimation and resolve the region of complex fibers in the neural structure, it is necessary to use the novel high angular resolution diffusion imaging (HARDI) technique; the spherical convolution method is an effective fiber orientation (FOD) estimation method, the The model is the convolution of the fiber orientation distribution response function; however, most of these methods estimate the FOD in a voxel-based manner, ignoring the spatial prior information related to the domain voxel; the FOD estimation error of the domain voxel may be due to the fiber Accumulated errors in tracking cause the reconstructed fibers to deviate from the actual fibers; to overcome the limitations, spatial HARDI or global reconstruction methods have been used to reconstruct the spatial constraints as parameters of the spatial signal and try to better describe the Fiber bundles are reconstructed while measuring data. The goal of global reconstruction is to provide a consistent view of the fiber architecture in the ensemble; however, these methods can only be constrained by linear smoothness between connected fiber segments.

发明内容Contents of the invention

为了克服已有纤维重构方法的准确性较低的不足，本发明提供一种准确性较高的基于全局稀疏正则化模型的纤维重构方法。In order to overcome the disadvantage of low accuracy of existing fiber reconstruction methods, the present invention provides a fiber reconstruction method based on a global sparse regularization model with high accuracy.

本发明解决其技术问题所采用的技术方案是：The technical solution adopted by the present invention to solve its technical problems is:

一种基于全局稀疏正则化模型的纤维重构方法，包括如下步骤：A fiber reconstruction method based on a global sparse regularization model, comprising the following steps:

1)建立基于字典基重构方法的局部稀疏模型1) Establish a local sparse model based on the dictionary base reconstruction method

扩散信号s(g|u)在无扩散加权下在位置v上和梯度方向g的测量标准化，其被表示为一个单一的纤维响应函数r(g,v)和纤维取向分布(FOD)f(v|u)的卷积：Diffusion signal s(g|u) normalized to measurements at position v and gradient direction g without diffusion weighting, which is expressed as a single fiber response function r(g,v) and fiber orientation distribution (FOD) f( Convolution of v|u):

$s the s ((g g | | u u)) = = {&Integral; &Integral;}_{{S S}_{22}} r r ((g g,, v v)) f f ((v v | | u u)) d d μ μ ((v v)) - - - - - - ((11))$

其中，u∈S²是在采样单元半球得到的中心向量组，μ(v)是哈尔测度，定义这些向量组中的一个组为纤维取向分布函数，应用于多壳方法中，单个纤维响应函数被定义为在这里表示表征扩散敏感系数b_i和各向异性相互作用影响程度的信号衰减，g_i表示第i个扩散梯度；球面去卷积方法是假设所有的纤维有相同的扩散性，因此在交叉构型以不同的形状轮廓描述纤维的条件下，最后的FOD描述为基函数混合的总和；近似的FOD模型表示为函数d_i中的线性加权组合：Among them, u∈S ² is the center vector group obtained in the sampling unit hemisphere, μ(v) is the Haar measure, and one of these vector groups is defined as the fiber orientation distribution function, which is applied to the multi-shell method, and the single fiber response function is defined as it's here Represents the signal attenuation that characterizes the diffusion sensitivity coefficient b _i and the degree of anisotropic interaction, g _i represents the i-th diffusion gradient; the spherical deconvolution method assumes that all fibers have the same diffusivity, so in the cross configuration as Under the condition that different shape profiles describe the fibers, the final FOD is described as the sum of the mixture of basis functions; the approximate FOD model is expressed as a linear weighted combination in the function d _i :

$f f ((v v | | u u)) = = {Σ Σ}_{i i = = 11}^{m m} {ω ω}_{i i} {d d}_{i i} ((v v,, u u)) - - - - - - ((22))$

其中，m为基函数字典(d₁,d₂,...,d_m)的基数，W＝[ω₁,ω₂,...,ω_m]^T是系数向量，ω_i是第i个系数i＝1...q，i,q都是系数；正标量W_i表示基函数的分布d_i(v,u)，采用球面双叶的基函数来表示FOD，如下所示：Among them, m is the cardinality of the basis function dictionary (d ₁ ,d ₂ ,...,d _m ), W=[ω ₁ ,ω ₂ ,...,ω _m ] ^T is the coefficient vector, ω _i is the i-th Coefficients i=1...q, i and q are all coefficients; the positive scalar W _i represents the distribution d _i (v, u) of the basis function, and the spherical bilobal basis function is used to represent the FOD, as shown below:

$d d ((θ θ,, τ τ)) = = {κ κ}_{11} {[[\frac{s the s i i n no θ θ}{11 - - {κ κ}_{22} {cos cos}^{22}}]]}^{τ τ} - - - - - - ((33))$

其中，标准化参数κ₁>0，κ₂∈(0,1)是调节峰值的参数，θ(u,v)＝|v^Tu,v∈S²是在单位球面上进行分割获得的旋转向量组，τ是指在径向分布的指数d(θ,τ)的增强数；字典基重构方法旨在恢复系数ω_i，标准化扩散信号S和观测矩阵Φ；似然分布表达的能量是指重构信号和所述FOD正表示数据的L2范数差，并通过非负最小二乘解决：Among them, the normalization parameter κ ₁ >0, κ ₂ ∈ (0,1) is the parameter to adjust the peak value, θ(u,v)=|v ^T u,v∈S ² is the rotation vector obtained by dividing on the unit sphere group, τ refers to the augmented number of exponent d(θ,τ) in the radial distribution; the dictionary basis reconstruction method aims to recover the coefficient ω _i , normalize the diffusion signal S and the observation matrix Φ; the energy expressed by the likelihood distribution is Reconstruct the L2-norm difference between the signal and the FOD positive representation data, and solve by non-negative least squares:

$\underset{W W}{m m i i n no} \frac{11}{22} | | | | S S - - Φ Φ W W | | {| |}_{22}^{22} s the s . . t t . . W W &GreaterEqual; &Greater Equal; 00 - - - - - - ((44))$

基数m大于样本大小n，默认f(v|u)的估计是稀疏的；假定一个体素内不超过三个纤维束；为了从多壳重构扩散FOD数据，公式(4)在q壳中扩展为如下：The base m is larger than the sample size n, and the estimation of f(v|u) is sparse by default; no more than three fiber bundles are assumed in a voxel; in order to reconstruct diffuse FOD data from multi-shell, formula (4) in q-shell expands as follows:

$\underset{W W}{min min} \frac{11}{22} | | | | \overset{~ ~}{S S} - - \overset{~ ~}{Φ Φ} W W | | {| |}_{22}^{22} s the s . . t t . . W W &GreaterEqual; &Greater Equal; 00 - - - - - - ((55))$

其中，是扩散信号强度的一个矩阵，是观测信号的矩阵，S_i,i＝1,2,...,q为扩散信号强度对第i个壳中的q空间测得的向量，Φ_i为第i个观测矩阵，是通过对Q空间第i个壳中测得的扩散信号强度系数；in, is a matrix of diffuse signal strengths, is the matrix of the observation signal, S _i ,i=1,2,...,q is the vector measured by the diffusion signal intensity to the q space in the i-th shell, Φ _i is the i-th observation matrix, which is obtained by The measured diffusion signal intensity coefficient in the ith shell of Q space;

2)建立全局模型2) Build a global model

全局模型合并空间信息转换成每个体素的方位分布估计，对于测量信号S_e和纤维模型的预测信号X_e，全局优化过程如下所示：The global model incorporates spatial information and converts it into an orientation distribution estimate for each voxel. For the measured signal Se and the predicted signal X _e of the fiber model _, the global optimization process is as follows:

$P P (({X x}_{e e} | | {S S}_{e e})) = = P P (({S S}_{e e} | | {X x}_{e e})) P P (({X x}_{e e})) = = {e e}^{- - {U u}^{I I n no}} \underset{i i}{Π Π} {e e}^{- - {β β}_{i i} {U u}_{i i}^{E E. x x}} - - - - - - ((66))$

通过从后验概率候选中检查可能的分布和采样获得的最优解P(X_e|S_e)，对于感兴趣区域ROI，ROI∈Z^ρ，ρ＝N_x×N_y×N_z，Z^ρ表示ROI区域内的体素，而N_x,N_y,N_z表示其x，y，z轴上的范围，内能U^In控制各体素的内部结构，表示为：The optimal solution P(X _e |S _e ) obtained by examining possible distributions and sampling from the posterior probability candidates, for a region of interest ROI, ROI∈Z ^ρ , ρ=N _x ×N _y ×N _z , Z ^ρ represents the voxel in the ROI area, and N _x , N _y , N _z represent the range on the x, y, and z axes, and the internal energy U ^In controls the internal structure of each voxel, expressed as:

${U u}^{I I n no} = = | | | | S S - - Φ Φ W W | | {| |}_{22}^{22} - - - - - - ((77))$

其中，S和Φ是先前描述的多壳参数的积分，表示如下：where S and Φ are integrals of the previously described multi-shell parameters expressed as follows:

$S S = = {[[{\overset{~ ~}{S S}}_{11}^{T T},, {\overset{~ ~}{S S}}_{22}^{T T},, ... ...,, {\overset{~ ~}{S S}}_{ρ ρ}^{T T}]]}^{T T}$

$W W = = {[[{W W}_{11}^{T T},, {W W}_{22}^{T T},, ... ...,, {W W}_{ρ ρ}^{T T}]]}^{T T} - - - - - - ((88))$

$Φ Φ = = {[\begin{matrix} \overset{~ ~}{Φ Φ} & 00 \\ \overset{~ ~}{Φ Φ} \\ ... ... \\ 00 & \overset{~ ~}{Φ Φ} \end{matrix}]}^{T T}$

而外部能量U^Ext表示一个特定信号和领域体素信号之间空间可能性关系，使得FOD的一致性为连接纤维取向的路径：While the external energy U ^Ext represents the spatial possibility relationship between a specific signal and field voxel signals such that the consistency of FOD is the path connecting the fiber orientation:

${U u}^{E E. x x t t} = = {Σ Σ}_{k k = = 11}^{ρ ρ} | | | | M m (({W W}_{k k} - - {\overset{&OverBar; &OverBar;}{W W}}_{{s the s}_{k k}})) | | {| |}_{22} + + W W &GreaterEqual; &Greater Equal; 00 - - - - - - ((99))$

其中，表示系数H的平均扩散系数，是周围系数；W_k代表系数向量W第k个系数代表白质区域基函数的径向和；W≥0是为了消除负的系数；in, Denotes the mean diffusion coefficient with coefficient H, is the surrounding coefficient; W _k represents the kth coefficient of the coefficient vector W Represents the radial sum of the basis functions of the white matter region; W≥0 is to eliminate negative coefficients;

3)全局优化算法的价值函数3) The value function of the global optimization algorithm

公式(6)中得到后验概率的最大化转化为内部和外部约束函数所组成的总能量函数的最小化，表示为：The maximization of the posterior probability obtained in formula (6) is transformed into the minimization of the total energy function composed of internal and external constraint functions, which is expressed as:

$arg arg max max P P (({X x}_{e e} | | {S S}_{e e})) \overset{Δ Δ}{= =} a a r r g g m m i i n no {{{U u}^{I I n no} + + \underset{i i}{Σ Σ} {β β}_{i i} {U u}_{i i}^{E E. x x t t}}} - - - - - - ((1010))$

其中，β_i>0,i＝1,2...是一个权重因子，是第i个外部能量，全局优化的成本函数被表示为：Among them, β _i >0, i=1,2...is a weight factor, is the i-th external energy, and the cost function for global optimization is expressed as:

$\underset{W W &GreaterEqual; &Greater Equal; 00}{m m i i n no} {{| | | | S S - - Φ Φ W W | | {| |}_{22}^{22} + + {β β}_{11} {Σ Σ}_{k k = = 11}^{ρ ρ} | | | | M m (({W W}_{k k} - - {\overset{&OverBar; &OverBar;}{W W}}_{{S S}_{k k}})) | | {| |}_{22}}} - - - - - - ((1111))$

定义一个局部优化问题：Define a local optimization problem:

$\underset{W W &GreaterEqual; &Greater Equal; 00}{min min} {{| | | | \overset{~ ~}{S S} - - \overset{~ ~}{Φ Φ} W W | | {| |}_{22}^{22} + + {β β}_{11} {Σ Σ}_{k k = = 11}^{ρ ρ} | | | | M m ((W W - - {\overset{&OverBar; &OverBar;}{W W}}_{S S})) | | | |}} - - - - - - ((1212))$

其中，是周围系数的平均值，通过使用增广拉格朗日方法解决，每个体素在最小化之后，所有的体素系数逐步更新；最终精确的FOD被表示为新的基函数加权和，如下所示：in, is the average value of the surrounding coefficients, which is solved by using the augmented Lagrangian method. After each voxel is minimized, all voxel coefficients are gradually updated; the final accurate FOD is expressed as a weighted sum of new basis functions, as follows Show:

${W W}^{* *} = = {(({\overset{~ ~}{Φ Φ}}^{T T} \overset{~ ~}{Φ Φ} + + {β β}_{11} {M m}^{T T} M m))}^{- - 11} {(({\overset{~ ~}{Φ Φ}}^{T T} \overset{~ ~}{S S} + + {\overset{&OverBar; &OverBar;}{W W}}_{S S}))}_{{W W}^{* *} &GreaterEqual; &Greater Equal; 00} - - - - - - ((1313))$

其中，W^*代表新的平均扩散系数，公式(5)选定ROI的每个体素和所有取得的W被存储为一个参考库，然后通过公式(13)获得的W^*值逐步替代初始值的参考库W；而库的更新是动态的，所以采集的W^*值会越来越准确。Among them, W ^* represents the new average diffusion coefficient, and each voxel of the selected ROI by formula (5) and all obtained W are stored as a reference library, and then the W ^* value obtained by formula (13) gradually replaces the initial value of Refer to the library W; and the update of the library is dynamic, so the collected W ^* value will become more and more accurate.

本发明的技术构思为：灵活性高的球面双叶基函数，它能够形成一个完备的库，来保证FOD的局部稀疏，同时建立一个全局模型。The technical idea of the present invention is: a highly flexible spherical double-leaf basis function, which can form a complete library to ensure the local sparseness of FOD and establish a global model at the same time.

在全局模型的基础上，提出优化成本函数，将体素系数逐步更新，最后得到一个相对准确的数据库。On the basis of the global model, an optimization cost function is proposed to update the voxel coefficients step by step, and finally a relatively accurate database is obtained.

本发明的有益效果主要表现在：准确性较高。The beneficial effects of the present invention are mainly manifested in: higher accuracy.

具体实施方式detailed description

下面对本发明作进一步描述。The present invention will be further described below.

1)建立基于字典基重构方法的局部稀疏模型：1) Establish a local sparse model based on the dictionary base reconstruction method:

其中u∈S²是在采样单元半球得到的中心向量组，μ(v)是哈尔测度，定义这些向量组中的一个组为纤维取向分布函数，应用于多壳方法中，单个纤维响应函数被定义为在这里表示表征扩散敏感系数b_i和各向异性相互作用影响程度的信号衰减，g_i表示第i个扩散梯度；球面去卷积方法是假设所有的纤维有相同的扩散性，因此在交叉构型以不同的形状轮廓描述纤维的条件下，最后的FOD可描述为基函数混合的总和；近似的FOD模型表示为函数d_i中的线性加权组合：where u∈S ² is the central vector group obtained in the sampling unit hemisphere, μ(v) is the Haar measure, and one of these vector groups is defined as the fiber orientation distribution function, which is applied to the multi-shell method, and the single fiber response function is defined as it's here Represents the signal attenuation that characterizes the diffusion sensitivity coefficient b _i and the degree of anisotropic interaction, g _i represents the i-th diffusion gradient; the spherical deconvolution method assumes that all fibers have the same diffusivity, so in the cross configuration as Under the condition that different shape profiles describe fibers, the final FOD can be described as the sum of the mixture of basis functions; the approximate FOD model is expressed as a linear weighted combination in the function d _i :

其中m为基函数字典(d₁,d₂,...,d_m)的基数，W＝[ω₁,ω₂,...,ω_m]^T是系数向量，ω_i是第i个系数i＝1...q，i,q都是系数；正标量W_i表示基函数的分布d_i(v,u)，本文提出了一种叫做球面双叶的基函数来表示FOD，如下所示：Where m is the cardinality of the basis function dictionary (d ₁ ,d ₂ ,...,d _m ), W=[ω ₁ ,ω ₂ ,...,ω _m ] ^T is the coefficient vector, ω _i is the ith The coefficient i=1...q, i, q are coefficients; the positive scalar W _i represents the distribution d _i (v,u) of the basis function, this paper proposes a basis function called spherical bilobal to represent FOD, as follows Shown:

其中标准化参数κ₁>0，κ₂∈(0,1)是调节峰值的参数，θ(u,v)＝|v^Tu|,v∈S²是在单位球面上进行分割获得的旋转向量组，τ是指在径向分布的指数d(θ,τ)的增强数；字典基重构方法旨在恢复系数ω_i，标准化扩散信号S和观测矩阵Φ；似然分布表达的能量是指重构信号和所述FOD正表示数据的L2范数差，并通过非负最小二乘(NNLS)解决：Among them, the standardized parameter κ ₁ >0, κ ₂ ∈ (0,1) is the parameter for adjusting the peak value, θ(u,v)=|v ^T u|, v∈S ² is the rotation vector obtained by dividing on the unit sphere group, τ refers to the augmented number of exponent d(θ,τ) in the radial distribution; the dictionary basis reconstruction method aims to recover the coefficient ω _i , normalize the diffusion signal S and the observation matrix Φ; the energy expressed by the likelihood distribution is The L2 norm difference between the reconstructed signal and the FOD positive representation data is solved by non-negative least squares (NNLS):

基数m可以大于样本大小n，我们默认f(v|u)的估计是稀疏的；我们假定一个体素内不超过三个纤维束，所以W坐标几乎是大于等于零的；为了从多壳重构扩散FOD数据，公式(4)可以在q壳中扩展为如下：The cardinality m can be larger than the sample size n, and we assume that the estimation of f(v|u) is sparse by default; we assume that there are no more than three fiber bundles in a voxel, so the W coordinate is almost greater than or equal to zero; in order to reconstruct from multi-shell Diffusion FOD data, formula (4) can be expanded in the q-shell as follows:

其中是扩散信号强度的一个矩阵，是观测信号的矩阵，S_i,i＝1,2,...,q为扩散信号强度对第i个壳中的q空间测得的向量，Φ_i为第i个观测矩阵，是通过对Q空间第i个壳中测得的扩散信号强度系数。in is a matrix of diffuse signal strengths, is the matrix of the observation signal, S _i ,i=1,2,...,q is the vector measured by the diffusion signal intensity to the q space in the i-th shell, Φ _i is the i-th observation matrix, which is obtained by The intensity coefficient of the diffusion signal measured in the i-th shell of Q-space.

2)建立全局模型2) Build a global model

从非负最小二乘法衍生的稀疏系数易受到稳定性影响，而逆问题被设置在单独的体素中，纤维跟踪的累积误差或意外信号的丢失可能导致实际纤维和重构纤维的偏差，这就是建立全局模型的目的；首先大多数全局重构方法通过线性平滑来实现，然而这种方法忽略了FOD的空间一致性；本文的全局模型合并空间信息转换成每个体素的方位分布估计，对于测量信号S_e和纤维模型的预测信号X_e，全局优化过程如下所示：Sparse coefficients derived from non-negative least squares are susceptible to stability, and while the inverse problem is set in individual voxels, cumulative errors in fiber tracking or loss of unexpected signals can lead to deviations in actual and reconstructed fibers, which It is the purpose of establishing a global model; first, most global reconstruction methods are implemented by linear smoothing, but this method ignores the spatial consistency of FOD; the global model in this paper combines spatial information and converts it into an orientation distribution estimate for each voxel, for The measurement signal Se and the predicted signal X _e of the fiber model _, the global optimization process is as follows:

本文的目的是探索模型的不同状态来确定最合适的数据纤维组合；具体而言，该模型是一个组合优化问题；通过从后验概率候选中检查可能的分布和采样获得的最优解P(X_e|S_e)，对于感兴趣区域(ROI)，ROI∈Z^ρ，ρ＝N_x×N_y×N_z，Z^ρ表示ROI区域内的体素，而N_x,N_y,N_z表示其x，y，z轴上的范围，内能U^In控制各体素的内部结构，可以表示为：The purpose of this paper is to explore different states of the model to determine the most suitable combination of data fibers; specifically, the model is a combinatorial optimization problem; the optimal solution P obtained by examining possible distributions and sampling from the posterior probability candidates ( X _e |S _e ), for the region of interest (ROI), ROI∈Z ^ρ , ρ=N _x ×N _y ×N _z , Z ^ρ represents the voxel in the ROI area, and N _x , N _y , N _z Indicates the range on the x, y, and z axes, and the internal energy U ^In controls the internal structure of each voxel, which can be expressed as:

其中S和Φ是先前描述的多壳参数的积分，表示如下：where S and Φ are integrals of the previously described multi-shell parameters expressed as follows:

表示系数H(这里H＝26)的平均扩散系数，是周围系数；W_k代表系数向量W第k个系数代表白质区域基函数的径向和；W≥0是为了消除负的系数。 Represents the mean diffusion coefficient of coefficient H (here H=26), is the surrounding coefficient; W _k represents the kth coefficient of the coefficient vector W Represents the radial sum of basis functions of white matter regions; W ≥ 0 is to eliminate negative coefficients.

公式(6)中可以得到后验概率的最大化可转化为内部和外部约束函数所组成的总能量函数的最小化，表示为：In formula (6), it can be obtained that the maximization of the posterior probability can be transformed into the minimization of the total energy function composed of internal and external constraint functions, expressed as:

其中β_i>0,i＝1,2...是一个权重因子，是第i个外部能量，我们设定为1，全局优化的成本函数被表示为：Where β _i >0, i=1,2... is a weighting factor, is the i-th external energy, we set it to 1, and the cost function of global optimization is expressed as:

定义一个局部优化问题：Define a local optimization problem:

是周围系数的平均值，通过使用增广拉格朗日方法解决，每个体素在最小化之后，所有的体素系数逐步更新；最终精确的FOD被表示为新的基函数加权和，如下所示： is the average value of the surrounding coefficients, which is solved by using the augmented Lagrangian method. After each voxel is minimized, all voxel coefficients are gradually updated; the final accurate FOD is expressed as a weighted sum of new basis functions, as follows Show:

W^*代表新的平均扩散系数，公式(5)选定ROI的每个体素和所有取得的W被存储为一个参考库，然后通过公式(13)获得的W^*值逐步替代初始值的参考库W；而库的更新是动态的，所以采集的W^*值将会越来越准确。W ^* represents the new average diffusion coefficient, each voxel of the selected ROI and all obtained W in formula (5) are stored as a reference library, and then the W ^* value obtained by formula (13) gradually replaces the initial value of the reference library W; and the update of the library is dynamic, so the collected W ^* value will become more and more accurate.

Claims

1. A fiber reconstruction method based on a global sparse regularization model is characterized in that: the reconstruction method comprises the following steps:

1) local sparse model established based on dictionary basis reconstruction method

The diffusion signal s (g | u) is normalized to the measurement of the gradient direction g at position v without diffusion weighting, which is expressed as a convolution of a single fiber response function r (g, v) and the fiber orientation distribution FOD f (v | u):

s (g | u) = {&Integral;}_{S_{2}} r (g, v) f (v | u) d μ (v) - - - (1)

wherein u ∈ S²Is a central vector set obtained in a sampling unit hemisphere, mu (v) is a haar measure, one of the vector sets is defined as a fiber orientation distribution function, and in the multi-shell method, a single fiber response function is defined asHerein, theRepresenting the characteristic diffusion sensitivity coefficient b_iSignal attenuation with degree of influence of anisotropic interaction, g_iRepresents the ith diffusion gradient; the spherical deconvolution method assumes that all fibers have the same diffusivity, so under the condition that the cross-over configuration describes the fibers with different shape profiles, the last FOD is described as the sum of the basis function mixtures; approximate FOD model is expressed as a function d_iLinear weighted combination of (1):

f (v | u) = Σ_{i = 1}^{m} ω_{i} d_{i} (v, u) - - - (2)

wherein m is a basis function dictionary (d)₁,d₂,...,d_m) Base of [ W ═ ω [ ]₁,ω₂,...,ω_m]^TIs a coefficient vector, ω_iQ, i, q are coefficients, i ═ 1.. q; positive scalar quantity W_iDistribution d representing basis functions_i(v, u), the FOD is expressed using the basis functions of the spherical double leaf as follows:

d (θ, τ) = κ_{1} {[\frac{s i n θ}{1 - κ_{2} \cos^{2}}]}^{τ} - - - (3)

wherein the normalization parameter κ₁>0，κ₂∈ (0,1) are parameters for adjusting the peak value,is a set of rotation vectors obtained by dividing the unit sphere, and τ is an index distributed in the radial directionThe number of enhancements of (d); dictionary basis weight construction method aiming at restoring coefficient omega_iNormalizing the diffusion signal S and the observation matrix phi; the energy expressed by the likelihood distribution refers to the L2 norm difference between the reconstructed signal and the FOD positive representation data and is solved by non-negative least squares:

\underset{W}{m i n} \frac{1}{2} | | S - Φ W | |_{2}^{2} s . t . W &GreaterEqual; 0 - - - (4)

the base m is larger than the sample size n, the estimate of default f (v | u) is sparse; assuming no more than three fiber bundles within a voxel; to reconstruct diffuse FOD data from multiple shells, equation (4) is extended in the q-shell as follows:

\underset{W}{m i n} \frac{1}{2} | | \tilde{S} - \tilde{Φ} W | |_{2}^{2} s . t . W &GreaterEqual; 0 - - - (5)

wherein,is a matrix of the intensity of the diffuse signal,is a matrix of observed signals, S_iQ is the vector of diffuse signal intensity measured in q space in the ith shell, Φ_iThe ith observation matrix is the diffusion signal intensity coefficient measured in the ith shell of the Q space;

2) building a global model

The global model incorporates the spatial information into an estimate of the orientation distribution of each voxel, for the measurement signal S_eAnd the predicted signal X of the fiber model_eThe global optimization process is as follows:

P (X_{e} | S_{e}) = P (S_{e} | X_{e}) P (X_{e}) = e^{- U^{I n}} \underset{i}{Π} e^{- β_{i} U_{i}^{E x}} - - - (6)

optimal solution P (X) obtained by examining possible distributions and samples from a posteriori probability candidates_e|S_e) For a region of interest ROI, ROI ∈ Z^ρ，ρ＝N_x×N_y×N_z，Z^ρRepresents voxels within the ROI region, and N_x,N_y,N_zRepresents the range on the x, y and z axes and the internal energy U^InThe internal structure of each voxel is controlled, represented as:

U^{I n} = | | S - Φ W | |_{2}^{2} - - - (7)

where S and Φ are the integrals of the previously described multi-shell parameters, expressed as follows:

S = {[{\tilde{S}}_{1}^{T}, {\tilde{S}}_{2}^{T}, ..., {\tilde{S}}_{ρ}^{T}]}^{T}

W = {[W_{1}^{T}, W_{2}^{T}, ..., W_{ρ}^{T}]}^{T} - - - (8)

Φ = {[\begin{matrix} \tilde{Φ} & 0 \\ \tilde{Φ} \\ ... \\ 0 & \tilde{Φ} \end{matrix}]}^{T}

and external energy U^ExtRepresenting the spatial likelihood relationship between one particular signal and the domain voxel signal such that the consistency of the FOD is the path connecting the fiber orientations:

U^{E x t} = Σ_{k = 1}^{ρ} | | M (W_{k} - {\overset{&OverBar;}{W}}_{s_{k}}) | |_{2} + W &GreaterEqual; 0 - - - (9)

wherein,represents the average diffusion coefficient of the coefficient H,is the ambient coefficient; w_kRepresenting the kth coefficient of a coefficient vector WRepresents the radial sum of basis functions of white matter regions; w is more than or equal to 0 to eliminate negative coefficient;

3) cost function of global optimization algorithm

The maximization of the posterior probability obtained in equation (6) translates into a minimization of the total energy function consisting of the internal and external constraint functions, expressed as:

\arg \max P (X_{e} | S_{e}) \overset{Δ}{=} argmin {U^{I n} + \underset{i}{Σ} β_{i} U_{i}^{E x t}} - - - (10)

wherein, β_i>0,iIs a weighting factor of 1,2,is the ith external energy, the cost function of the global optimization is expressed as:

\underset{W &GreaterEqual; 0}{m i n} {| | S - Φ W | |_{2}^{2} + β_{1} Σ_{k = 1}^{ρ} | | M (W_{k} - {\overset{&OverBar;}{W}}_{S_{k}}) | |_{2}} - - - (11)

defining a local optimization problem:

\min_{W &GreaterEqual; 0} {| | \tilde{S} - \tilde{Φ} W | |_{2}^{2} + β_{1} Σ_{k = 1}^{ρ} | | M (W - {\overset{&OverBar;}{W}}_{S}) | |} - - - (12)

wherein,is the average of the surrounding coefficients, solved by using the augmented lagrange method, after each voxel is minimized, all voxel coefficients are updated step by step; the final exact FOD is expressed as a new weighted sum of basis functions, as follows:

W^{*} = {({\tilde{Φ}}^{T} \tilde{Φ} + β_{1} M^{T} M)}^{- 1} {({\tilde{Φ}}^{T} \tilde{S} + {\overset{&OverBar;}{W}}_{S})}_{W^{*} &GreaterEqual; 0} - - - (13)

wherein, W^*Representing the new mean diffusion coefficient, each voxel of the ROI selected by equation (5) and all W taken are stored as a reference library, and then W is obtained by equation (13)^*The values gradually replace the reference library W of the initial values; while the updating of the library is dynamic, so the W collected^*The values will become more and more accurate.