CN105913465A - Overall sparsity regularization model-based fiber reconstructing method - Google Patents

Abstract

An overall sparsity regularization model-based fiber reconstructing method is disclosed and comprises the following steps: 1) a dictionary base reconstruction method-based local sparsity model is built; 2) an overall model is built; 3) a value function of an overall optimization algorithm is obtained. The invention provides the overall sparsity regularization model-based fiber reconstructing method which is high in accuracy.

Description

Fiber reconstruction method based on global sparse regularization model

Technical Field

The invention relates to the field of medical imaging and neuroanatomy under computer graphics, in particular to a fiber reconstruction method.

Background

In order to establish an accurate fiber orientation estimation method, a novel High Angular Resolution Diffusion Imaging (HARDI) technology must be used for solving the area of the neural structure composite fiber; the spherical convolution method is an effective Fiber Orientation (FOD) estimation method, and the model is the convolution of fiber orientation distribution response functions; however, these methods mostly estimate FOD in a voxel-based manner, ignoring spatial prior information with the domain voxels; FOD estimation errors of domain voxels may cause the reconstructed fiber to deviate from the actual fiber due to accumulated errors of fiber tracking; to overcome the limitations, spatial HARDI or global reconstruction methods have been used to reconstruct spatial constraints as parameters of the spatial signal and attempt to reconstruct the fiber bundle while better describing the measurement data by determining their structure. The purpose of global reconstruction is to provide a consistent view of the fiber architecture in the whole; however, these methods can only be constrained by linear smoothness between the connecting fiber segments.

Disclosure of Invention

In order to overcome the defect of low accuracy of the conventional fiber reconstruction method, the invention provides the fiber reconstruction method based on the global sparse regularization model with high accuracy.

The technical scheme adopted by the invention for solving the technical problems is as follows:

a fiber reconstruction method based on a global sparse regularization model 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\left(g\right|u)={\∫}_{S_2}r(g,v)f\left(v\right|u)d\mathrm{\μ}\left(v\right)---\left(1\right)$

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\left(v\right|u)=\underset{i=1}{\overset{m}{\Σ}}{\mathrm{\ω}}_i{d}_i(v,u)---\left(2\right)$

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(\mathrm{\θ},\mathrm{\τ})={\mathrm{\κ}}_1{\[\frac{sin\mathrm{\θ}}{1-{\mathrm{\κ}}_2{\cos }^2}\]}^{\mathrm{\τ}}---\left(3\right)$

wherein the normalization parameter κ₁>0，κ₂∈ (0,1) is a parameter for adjusting the peak value, θ (u, v) ═ v^Tu,v∈S²Is a set of rotation vectors obtained by dividing a unit spherical surface, τ is the number of enhancements of the index d (θ, τ) distributed in the radial direction; 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}{min}\frac{1}{2}\left|\right|S-\mathrm{\Φ}W|{|}_2^{2}s.t.W\≥0---\left(4\right)$

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}{\min }\frac{1}{2}\left|\right|\tilde{S}-\widetilde{\mathrm{\Φ}}W|{|}_2^{2}s.t.W\≥0---\left(5\right)$

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\left(X_e\right|S_e)=P\left(S_e\right|X_e)P\left(X_e\right)=e^{-U^{In}}\underset{i}{\Π}{e}^{-{\mathrm{\β}}_i{U}_i^{Ex}}---\left(6\right)$

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^{In}=\left|\right|S-\mathrm{\Φ}W|{|}_2^2---\left(7\right)$

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}}_{\mathrm{\ρ}}^T\]}^T$

$W={\[W_1^T,W_2^T,...,W_{\mathrm{\ρ}}^T\]}^T---\left(8\right)$

$\mathrm{\Φ}={\left[\begin{array}{cccc}\widetilde{\mathrm{\Φ}}& & & 0\\ & \widetilde{\mathrm{\Φ}}& & \\ & & ...& \\ 0& & & \widetilde{\mathrm{\Φ}}\end{array}\right]}^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^{Ext}=\underset{k=1}{\overset{\mathrm{\ρ}}{\Σ}}\left|\right|M(W_k-{\stackrel{\‾}{W}}_{s_k})|{|}_2+W\≥0---\left(9\right)$

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\left(X_e\right|S_e)\stackrel{\mathrm{\Δ}}{=}argmin\left\{U^{In}+\underset{i}{\mathrm{\Σ}}{\mathrm{\β}}_i{U}_i^{Ext}\right\}---\left(10\right)$

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

$\underset{W\≥0}{min}\left\{\right||S-\mathrm{\Φ}W|{|}_2^2+{\mathrm{\β}}_1\underset{k=1}{\overset{\mathrm{\ρ}}{\Σ}}\left|\right|M(W_k-{\stackrel{\‾}{W}}_{S_k})\left|{|}_2\right\}---\left(11\right)$

defining a local optimization problem:

$\underset{W\≥0}{\min }\left\{\right||\tilde{S}-\widetilde{\mathrm{\Φ}}W|{|}_2^2+{\mathrm{\β}}_1\underset{k=1}{\overset{\mathrm{\ρ}}{\Σ}}\left|\right|M(W-{\stackrel{\‾}{W}}_S)\left|\right|\}---\left(12\right)$

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^{*}={({\widetilde{\mathrm{\Φ}}}^T\widetilde{\mathrm{\Φ}}+{\mathrm{\β}}_1{M}^{T}M)}^{-1}{({\widetilde{\mathrm{\Φ}}}^T\tilde{S}+{\stackrel{\‾}{W}}_S)}_{W^{*}\≥0}---\left(13\right)$

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.

The technical conception of the invention is as follows: the spherical double-leaf basis function with high flexibility can form a complete library to ensure the local sparsity of FOD and establish a global model.

And on the basis of the global model, an optimization cost function is provided, the voxel coefficients are gradually updated, and finally a relatively accurate database is obtained.

The invention has the following beneficial effects: the accuracy is higher.

Detailed Description

The invention is further described below.

A fiber reconstruction method based on a global sparse regularization model comprises the following steps:

1) establishing a local sparse model based on a 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\left(g\right|u)={\∫}_{S_2}r(g,v)f\left(v\right|u)d\mathrm{\μ}\left(v\right)---\left(1\right)$

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 are in phaseThe same diffusivity, so that the final FOD can be described as the sum of the mixing of basis functions, provided that the cross-over configuration describes the fiber with different shape profiles; approximate FOD model is expressed as a function d_iLinear weighted combination of (1):

$f\left(v\right|u)=\underset{i=1}{\overset{m}{\Σ}}{\mathrm{\ω}}_i{d}_i(v,u)---\left(2\right)$

where m is a dictionary of basis functions (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), a basis function called spherical bilobe is proposed herein to represent FOD, as follows:

$d(\mathrm{\θ},\mathrm{\τ})={\mathrm{\κ}}_1{\[\frac{sin\mathrm{\θ}}{1-{\mathrm{\κ}}_2{\cos }^2}\]}^{\mathrm{\τ}}---\left(3\right)$

wherein the normalization parameter κ₁>0，κ₂∈ (0,1) is a parameter for adjusting the peak value, θ (u, v) ═ v^Tu|,v∈S²Is a set of rotation vectors obtained by dividing a unit spherical surface, τ is the number of enhancements of the index d (θ, τ) distributed in the radial direction; 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 (NNLS):

$\underset{W}{min}\frac{1}{2}\left|\right|S-\mathrm{\Φ}W|{|}_2^{2}s.t.W\≥0---\left(4\right)$

the base m can be larger than the sample size n, we default that the estimate of f (v | u) is sparse; we assume that no more than three fiber bundles are within a voxel, so the W coordinate is almost zero or greater; to reconstruct diffuse FOD data from multiple shells, equation (4) can be extended in the q-shell as follows:

$\underset{W}{\min }\frac{1}{2}\left|\right|\tilde{S}-\widetilde{\mathrm{\Φ}}W|{|}_2^{2}s.t.W\≥0---\left(5\right)$

whereinIs 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

Sparse coefficients derived from a non-negative least squares method are susceptible to stability, while the inverse problem is set in individual voxels, and accumulated errors in fiber tracking or loss of unexpected signals may cause deviations of actual fibers and reconstructed fibers, which is the purpose of building a global model; most global reconstruction methods are first implemented by linear smoothing, however this method ignores the spatial consistency of the FOD; the global model herein incorporates spatial information into an estimate of the orientation distribution of each voxel, for the measured signal S_eAnd the predicted signal X of the fiber model_eGlobal excellenceThe chemical process is as follows:

$P\left(X_e\right|S_e)=P\left(S_e\right|X_e)P\left(X_e\right)=e^{-U^{In}}\underset{i}{\Π}{e}^{-{\mathrm{\β}}_i{U}_i^{Ex}}---\left(6\right)$

the aim here is to explore the different states of the model to determine the most suitable data fiber combination; in particular, the model is a combinatorial optimization problem; 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^InControlling the internal structure of each voxel can be expressed as:

$U^{In}=\left|\right|S-\mathrm{\Φ}W|{|}_2^2---\left(7\right)$

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}}_{\mathrm{\ρ}}^T\]}^T$

$W={\[W_1^T,W_2^T,...,W_{\mathrm{\ρ}}^T\]}^T---\left(8\right)$

$\mathrm{\Φ}={\left[\begin{array}{cccc}\widetilde{\mathrm{\Φ}}& & & 0\\ & \widetilde{\mathrm{\Φ}}& & \\ & & ...& \\ 0& & & \widetilde{\mathrm{\Φ}}\end{array}\right]}^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^{Ext}=\underset{k=1}{\overset{\mathrm{\ρ}}{\Σ}}\left|\right|M(W_k-{\stackrel{\‾}{W}}_{s_k})|{|}_2+W\≥0---\left(9\right)$

represents the average diffusion coefficient of coefficient H (where H is 26),is the ambient coefficient; w_kRepresenting the kth coefficient of a coefficient vector WRepresentsRadial sum of white matter region basis functions; w ≧ 0 is for eliminating the negative coefficient.

3) Cost function of global optimization algorithm

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

$\arg \max P\left(X_e\right|S_e)\stackrel{\mathrm{\Δ}}{=}argmin\left\{U^{In}+\underset{i}{\mathrm{\Σ}}{\mathrm{\β}}_i{U}_i^{Ext}\right\}---\left(10\right)$

β therein_i>0, i-1, 2.. is a weighting factor,is the ith external energy, we set to 1, the cost function of the global optimization is expressed as:

$\underset{W\≥0}{min}\left\{\right||S-\mathrm{\Φ}W|{|}_2^2+{\mathrm{\β}}_1\underset{k=1}{\overset{\mathrm{\ρ}}{\Σ}}\left|\right|M(W_k-{\stackrel{\‾}{W}}_{S_k})\left|{|}_2\right\}---\left(11\right)$

defining a local optimization problem:

$\underset{W\≥0}{\min }\left\{\right||\tilde{S}-\widetilde{\mathrm{\Φ}}W|{|}_2^2+{\mathrm{\β}}_1\underset{k=1}{\overset{\mathrm{\ρ}}{\Σ}}\left|\right|M(W-{\stackrel{\‾}{W}}_S)\left|\right|\}---\left(12\right)$

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^{*}={({\widetilde{\mathrm{\Φ}}}^T\widetilde{\mathrm{\Φ}}+{\mathrm{\β}}_1{M}^{T}M)}^{-1}{({\widetilde{\mathrm{\Φ}}}^T\tilde{S}+{\stackrel{\‾}{W}}_S)}_{W^{*}\≥0}---\left(13\right)$

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 value will become more and more accurate.

Claims (1)

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\left(g\right|u)={\∫}_{S_2}r(g,v)f\left(v\right|u)d\mathrm{\μ}\left(v\right)---\left(1\right)$

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\left(v\right|u)=\underset{i=1}{\overset{m}{\Σ}}{\mathrm{\ω}}_i{d}_i(v,u)---\left(2\right)$

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(\mathrm{\θ},\mathrm{\τ})={\mathrm{\κ}}_1{\[\frac{sin\mathrm{\θ}}{1-{\mathrm{\κ}}_2{\cos }^2}\]}^{\mathrm{\τ}}---\left(3\right)$

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}{min}\frac{1}{2}\left|\right|S-\mathrm{\Φ}W|{|}_2^{2}s.t.W\≥0---\left(4\right)$

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}{min}\frac{1}{2}\left|\right|\tilde{S}-\widetilde{\mathrm{\Φ}}W|{|}_2^{2}s.t.W\≥0---\left(5\right)$

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\left(X_e\right|S_e)=P\left(S_e\right|X_e)P\left(X_e\right)=e^{-U^{In}}\underset{i}{\Π}{e}^{-{\mathrm{\β}}_i{U}_i^{Ex}}---\left(6\right)$

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^{In}=\left|\right|S-\mathrm{\Φ}W|{|}_2^2---\left(7\right)$

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}}_{\mathrm{\ρ}}^T\]}^T$

$W={\[W_1^T,W_2^T,...,W_{\mathrm{\ρ}}^T\]}^T---\left(8\right)$

$\mathrm{\Φ}={\left[\begin{array}{cccc}\widetilde{\mathrm{\Φ}}& & & 0\\ & \widetilde{\mathrm{\Φ}}& & \\ & & ...& \\ 0& & & \widetilde{\mathrm{\Φ}}\end{array}\right]}^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^{Ext}=\underset{k=1}{\overset{\mathrm{\ρ}}{\Σ}}\left|\right|M(W_k-{\stackrel{\‾}{W}}_{s_k})|{|}_2+W\≥0---\left(9\right)$

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\left(X_e\right|S_e)\stackrel{\mathrm{\Δ}}{=}\mathrm{argmin}\left\{U^{In}+\underset{i}{\mathrm{\Σ}}{\mathrm{\β}}_i{U}_i^{Ext}\right\}---\left(10\right)$

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\≥0}{min}\left\{\right||S-\mathrm{\Φ}W|{|}_2^2+{\mathrm{\β}}_1\underset{k=1}{\overset{\mathrm{\ρ}}{\Σ}}\left|\right|M(W_k-{\stackrel{\‾}{W}}_{S_k})\left|{|}_2\right\}---\left(11\right)$

defining a local optimization problem:

$\underset{W\≥0}{\min }\left\{\right||\tilde{S}-\widetilde{\mathrm{\Φ}}W|{|}_2^2+{\mathrm{\β}}_1\underset{k=1}{\overset{\mathrm{\ρ}}{\Σ}}\left|\right|M(W-{\stackrel{\‾}{W}}_S)\left|\right|\}---\left(12\right)$

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^{*}={({\widetilde{\mathrm{\Φ}}}^T\widetilde{\mathrm{\Φ}}+{\mathrm{\β}}_1{M}^{T}M)}^{-1}{({\widetilde{\mathrm{\Φ}}}^T\tilde{S}+{\stackrel{\‾}{W}}_S)}_{W^{*}\≥0}---\left(13\right)$

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.