CN105913480A - Space structure consistency-based brain fiber microstructure reconstructing method - Google Patents

Abstract

A space structure consistency-based brain fiber microstructure reconstructing method is disclosed and comprises the following steps: in a first step, a sphere dictionary deconvolution-based double valve primary function is adopted; in a second step, a space structure consistency model is built, FOD value distribution can be fitted out in a simulated manner via use of mathematic software MATLAB based on a separable space domain W(t+1)c, a main direction of fibers can be obtained by searching an extreme value point among FOD values, and therefore brain fiber microstructure reconstruction can be realized. Via use of the space structure consistency-based brain fiber microstructure reconstructing method which is provided in the invention, calculating efficiency can be improved, and high precision is realized.

Description

A kind of based on space structure conforming brain fiber microstructure reconstructing method

Technical field

The present invention relates to the medical imaging under computer graphics, neuroanatomy field, especially It it is a kind of brain fiber microstructure reconstructing method.

Background technology

Along with the development in epoch, the progress of Medical Imaging Technology, diffusion tensor imaging is god Accounting for increasing power of influence in the research of science, the neuroimaging technology having advanced person is this The individual epoch are indispensable；Diffusion tensor imaging is as emerging a kind of side describing brain structure Method, is also the method for unique a kind of In vivo detection human brain structure, neuromedicine field master simultaneously If the research to cerebral tissue architectural feature；At present, diffusion tensor imaging is just by widely It is applied to the supplementary means of psychiatric condition and diagnosis, it might even be possible to for pre-operative surgical scheme Formulate, it may be said that it has, in the contribution of medical domain, the advantage that can not be substituted；So to based on The algorithm research of diffusion tensor has great meaning for brain science.

Brain white matter integrity directional spreding Model Reconstruction is one of significant process of brain fiber imaging, for Fibre bundle is followed the tracks of provides accurate machine direction to estimate.The constraints of traditional method often relies on In the machine direction information of priori, limit the raising of computational efficiency and precision.Propose new more It is the focus of research for advanced brain white matter integrity directional spreding model.

Summary of the invention

In order to the computational efficiency overcoming existing brain fiber microstructure reconstructing method is relatively low, precision is relatively low Deficiency, the present invention provide a kind of promote computational efficiency, precision higher based on space structure one The brain fiber microstructure reconstructing method of cause property.

The technical solution adopted for the present invention to solve the technical problems is:

A kind of based on space structure conforming brain fiber microstructure reconstructing method, including walking as follows Rapid:

Step one, uses bivalve basic function based on sphere dictionary deconvolution:

At rotating vectorAnd center vectorOn fiber probability-distribution function F (v | u) it is referred to as fODF, wherein m_vRepresent the Spatial Dimension of rotating vector group, m_uRepresent be The Spatial Dimension of center vector group, is described as Fiber morphology structure by sphere the Method of Deconvolution The convolution of kernel function, at diffusion gradient g ∈ S²On measurement signal s (g | u):

$s\left(g\right|u)={\∫}_{S^2}r(g,v)f\left(v\right|u)d\mathrm{\μ}\left(v\right)$

Wherein (g, v) represents kernel function to r, and μ (v) is at S²On Haar estimate；For the sake of Fang Bian, S (g | u) signal can be represented by the discrete sampling μ being evenly distributed in unit sphere:

$s\left(g\right|u)=\underset{i=1}{\overset{m_u}{\Σ}}r(g,v)f\left(v\right|u_i)$

Definition fiber receptance function:

$r(g,v)=e^{-l{\left({\mathrm{gv}}^T\right)}^2}$

Wherein It is the diagonal matrix characterizing diffusion tensor D, λ_maxBeing principal eigenvector in this diagonal matrix, l represents diffusion-sensitive coefficient and anisotropy Interact the influence degree to signal attenuation；The fibre on many shells sampling plan is expanded with this The estimation of dimension receptance function:

$r(g,v)=\{r(g_{b_J},v)=e^{-l_{b_J}{\left(g_{b_J}{v}^T\right)}^2}|J=1,...,q\}$

Wherein J=1 ..., q refers to j-th spherical shell,Correspondence is distributed in b_JGradient vector on shell Collection g, l_bJRepresent the l that the bJ diffusion-sensitive coefficient under j-th spherical shell is referred to；Thus Push away to obtain a kind of new machine direction distribution function form:

$f\left(v\right|u)=\underset{i=1}{\overset{m_u}{\Σ}}{w}_{i}d(v,u_i),u_i\∈u$

d(v,u_i) represent Vector Groups v and u_iOne group of directional spreding base the most complete under individual direction Function, w_iIt is the relative weighting of basic function,It is w_iThe number of middle nonzero element；Apparent diffusion Coefficient can be by approximate evaluation:

g^TDg≈λ_maxg^T(vv^T) g=λ_maxcos²θ(g,v)

(g, v) represents included angle cosine function to θ, and D is to represent by all of d (v, u_i) constitute one group Basic function dictionary；In order to describe the dictionary base distribution under spherical coordinate, original right angle is sat by we D (v, u under mark system_i) be described as under spherical coordinate

It is the minimum angle cosine on discrete set, κ₁Machine direction distribution function is normalized to list Position ball, κ₂Being the parameter for adjusting peak value, τ represents even power；Machine direction distribution function (fODF) estimation of basic function coefficient w in:

$\underset{w}{\mathrm{argmin}}{||\mathrm{\Φ}w-s||}_2^2$

Φ is observing matrix, and s is signal vector；For avoiding pseudo-peak and complicated algorithm and more The calculating of the extensive ill-condition problem of high-order, directly tries to achieve base by non-negative least square method The estimation of function coefficients w:

$w^{*}\←\underset{w\≥0}{\arg }min{||\mathrm{\Φ}w-s||}_2^2$

w^*Represent the optimal solution of w；

Step 2, sets up space structure consistency model:

Utilize Bayesian formula, obtain:

P(x|s)∝P(s|x)P(x)

Posterior probability density P (x | s) is proportional to data likelihood function P (s | x) and priori probability density The product of function P (x)；It is rewritten as afterwards:

$P\left(x\right|s)\∝1/(U^{In}+\underset{i}{\Σ}{\mathrm{\β}}_i{U}_i^{Ex})$

U^InIt is internal energy, U^ExBeing external energy, β is the hyper parameter of prior distribution；Posteriority The maximization of probability is converted into minimizing of total energy function:

$\arg \max P\left(x\right|s)\stackrel{\mathrm{\Δ}}{=}\mathrm{argmin}\{U^{In}+\underset{i}{\Σ}{\mathrm{\β}}_i{U}_i^{Ex}\}$

P (x | s) it is posterior probability density, U^InIt is internal energy, U^ExBeing external energy, β is first Test the hyper parameter of distribution；Wherein internal energy:

$U^{In}\∝{||S-S^{\′}||}_2^2={||S-\mathrm{\Θ}W||}_2^2$

S is to measure signal collection, and what S ' represented is signal to be estimated, and W is w coefficient sets, and Θ is Block diagonal matrix based on observing matrix；External energy:

$U^{Ex}=\underset{c\∈\mathrm{\Ω},\tilde{c}\∈H_c\⋐\mathrm{\Ω}}{\Σ}{||M(w_c-{\stackrel{\‾}{w}}_{\tilde{c}})||}_2$

U^ExRepresent external energy,Representing the arithmetic average of machine direction distribution function, M leads to Cross down-sampled direction v_tGo to train a dictionary base to obtain, w_cRepresent c Ω voxel of ∈ is Number；The estimation of W:

$\begin{array}{cc}\underset{w_c\≥0}{\mathrm{argmin}}\{{||s_c-{\mathrm{\Φw}}_c||}_2^2+{\mathrm{\β}}_c{||M(w_c-{\stackrel{\‾}{w}}_{\tilde{c}})||}_2\}& c\∈\mathrm{\Ω};\tilde{c}\∈H_c\⋐\mathrm{\Ω}\end{array}$

s_cIt is to measure signal coefficient in voxel c, w_cIt it is dictionary coefficient in voxel c；For Obtain the structure of voxel decussating fibers, define global cost function:

$\underset{W\≥0}{\arg \min }\{{||S-\mathrm{\Θ}W||}_2^2+\underset{c\∈\mathrm{\Ω}}{\Σ}(\underset{\tilde{c}\∈H_c\⋐\mathrm{\Ω}}{\Σ}{\mathrm{\σ}}_{c,\tilde{c}}{||w_c-w_{\tilde{c}}||}_2+{\mathrm{\β}}_1||W\left(w_c-{\tilde{w}}_{\tilde{c}}\right)||+J\{w_c|c\∈\mathrm{\Ω}\})\}$

By calculating w_cAnd w_cThe coefficient COS distance of surrounding neighbors obtains, β₁It is artificial fixed One parameter of justice, Q is by training dictionary to obtain on basic function；Space structure consistency model Local linear approximate evaluation:

$\begin{array}{c}(w_c^{(t+1)},{\mathrm{\ξ}}_c^{(t+1)})=\underset{{\mathrm{\ξ}}_c^{(t+1)},w_c^{(t+1)}\≥0}{\arg \min }\{{||s_c-{\mathrm{\Φw}}_c^{\left(t\right)}||}_2^2+\underset{\tilde{c}\∈H_c\⋐\mathrm{\Ω}}{\Σ}{\mathrm{\σ}}_{c,\tilde{c}}{||w_c^{\left(t\right)}-w_{\tilde{c}}||}_2\\ +{\mathrm{\β}}_1{||Q(w_c^{\left(t\right)}-{\tilde{w}}_{\tilde{c}})||}_2+{\mathrm{\β}}_2{||{\mathrm{\ξ}}_c^{\left(t\right)}||}_1+{\mathrm{\δ}}_{\mathrm{\ξ}}{||w_c^{\left(t\right)}-{\mathrm{\ξ}}_c^{\left(t\right)}+p_{c,\mathrm{\ξ}}^{\left(t\right)}||}_2\}\\ p_{c,\mathrm{\ξ}}^{(t+1)}=p_{c,\mathrm{\ξ}}^{\left(t\right)}+w_c^{(t+1)}-{\mathrm{\ξ}}_c^{(t+1)}\end{array}$

T is iteration index, δ_ξIt is predefined aiding constant,It is that the extension of the t time iteration is drawn Ge Lang multiplier vector, above formula is divided into two parts:

$\begin{array}{c}{w}_c^{(t+1)}=\underset{w_c^{(t+1)}\≥0}{\arg \min }\{{||s_c-{\mathrm{\Φw}}_c^{\left(t\right)}||}_2^2+\underset{\tilde{c}\∈H_c\⋐\mathrm{\Ω}}{\Σ}{\mathrm{\σ}}_{c,\tilde{c}}{||w_c^{\left(t\right)}-w_{\tilde{c}}||}_2\\ +{\mathrm{\β}}_1{||Q(w_c^{\left(t\right)}-{\tilde{w}}_{\tilde{c}})||}_2+{\mathrm{\δ}}_{\mathrm{\ξ}}{||w_c^{\left(t\right)}-{\mathrm{\ξ}}_c^{\left(t\right)}+p_{c,\mathrm{\ξ}}^{\left(t\right)}||}_2\}\end{array}$

${\mathrm{\ξ}}_c^{(t+1)}=\underset{{\mathrm{\ξ}}_c^{(t+1)}\≥0}{\mathrm{argmin}}\{{\mathrm{\δ}}_{\mathrm{\ξ}}{||w_c^{\left(t\right)}-{\mathrm{\ξ}}_c^{\left(t\right)}+p_{c,\mathrm{\ξ}}^{\left(t\right)}||}_2+{\mathrm{\β}}_2{||{\mathrm{\ξ}}_c^{\left(t\right)}||}_1\}$

It is optimized for a separable space territory, solves with strengthening Lagrangian method, obtain:

$\begin{array}{c}{w}_c^{(t+1)}={({\mathrm{\Φ}}^T\mathrm{\Φ}+\underset{\tilde{c}\∈H_c\⋐\mathrm{\Ω}}{\Σ}{\mathrm{\σ}}_{c,\tilde{c}}I+{\mathrm{\β}}_c{Q}^{T}Q+{\mathrm{\δ}}_{\mathrm{\ξ}}I)}^{-1}\·\\ {({\mathrm{\Φ}}^T{s}_c+\underset{\tilde{c}\∈H_c\⋐\mathrm{\Ω}}{\Σ}{\mathrm{\σ}}_{c,\tilde{c}}{w}_{\tilde{c}}+{\mathrm{\β}}_1{Q}^{T}Q{\tilde{w}}_{\tilde{c}}+{\mathrm{\δ}}_{\mathrm{\ξ}}({\mathrm{\ξ}}_c^{\left(t\right)}-p_{c,\mathrm{\ξ}}^{\left(t\right)}\left)I\right)}_{w_c^{(t+1)}\≥0}\end{array}$

Wherein I representation unit matrix；

Based on separable space territoryMathematics soft MATLAB emulation is used to simulate FOD value Distribution, obtain the principal direction of fiber by the extreme point in search FOD value.

The technology of the present invention is contemplated that: spherical bivalve base letter apparent, more flexible, more effective Number, and one complete dictionary of mistake of formation characterizes the fibre of many shells basic function weighting on this basis Dimension direction distribution function (is called for short fODF).

Beneficial effects of the present invention is mainly manifested in: lifting computational efficiency, precision are higher.

Detailed description of the invention

The invention will be further described below.

A kind of based on space structure conforming brain fiber microstructure reconstructing method, including walking as follows Rapid:

Step one, uses bivalve basic function based on sphere dictionary deconvolution:

At rotating vectorAnd center vectorOn fiber probability-distribution function F (v | u) it is referred to as fODF, wherein m_vRepresent the Spatial Dimension of rotating vector group, m_uRepresent be The Spatial Dimension of center vector group, is described as Fiber morphology structure by sphere the Method of Deconvolution The convolution of kernel function, at diffusion gradient g ∈ S²On measurement signal s (g | u):

$s\left(g\right|u)={\∫}_{S^2}r(g,v)f\left(v\right|u)d\mathrm{\μ}\left(v\right)$

Wherein (g, v) represents kernel function to r, and μ (v) is at S²On Haar estimate；For the sake of Fang Bian, S (g | u) signal can be represented by the discrete sampling μ being evenly distributed in unit sphere:

$s\left(g\right|u)=\underset{i=1}{\overset{m_u}{\Σ}}r(g,v)f\left(v\right|u_i)$

Definition fiber receptance function:

$r(g,v)=e^{-l{\left({\mathrm{gv}}^T\right)}^2}$

Wherein It is the diagonal matrix characterizing diffusion tensor D, λ_maxBeing principal eigenvector in this diagonal matrix, l represents diffusion-sensitive coefficient and anisotropy Interact the influence degree to signal attenuation；The fibre on many shells sampling plan is expanded with this The estimation of dimension receptance function:

$r(g,v)=\{r(g_{b_J},v)=e^{-l_{b_J}{\left(g_{b_J}{v}^T\right)}^2}|J=1,...,q\}$

Wherein J=1 ..., q refers to j-th spherical shell,Correspondence is distributed in b_JGradient vector on shell Collection g, l_bJRepresent the l that the bJ diffusion-sensitive coefficient under j-th spherical shell is referred to；Thus Push away to obtain a kind of new machine direction distribution function form:

$f\left(v\right|u)=\underset{i=1}{\overset{m_u}{\Σ}}{w}_{i}d(v,u_i),u_i\∈u$

d(v,u_i) represent Vector Groups v and u_iOne group of directional spreding base the most complete under individual direction Function, w_iIt is the relative weighting of basic function,It is w_iThe number of middle nonzero element；Apparent diffusion Coefficient can be by approximate evaluation:

g^TDg≈λ_maxg^T(vv^T) g=λ_maxcos²θ(g,v)

(g, v) represents included angle cosine function to θ, and D is to represent by all of d (v, u_i) constitute one group Basic function dictionary；In order to describe the dictionary base distribution under spherical coordinate, original right angle is sat by we D (v, u under mark system_i) be described as under spherical coordinate

It is the minimum angle cosine on discrete set, κ₁Machine direction distribution function is normalized to list Position ball, κ₂Being the parameter for adjusting peak value, τ represents even power；Machine direction distribution function (fODF) estimation of basic function coefficient w in:

$\underset{w}{\mathrm{argmin}}{||\mathrm{\Φ}w-s||}_2^2$

Φ is observing matrix, and s is signal vector；For avoiding pseudo-peak and complicated algorithm and more The calculating of the extensive ill-condition problem of high-order, directly tries to achieve base by non-negative least square method The estimation of function coefficients w:

$w^{*}\←\underset{w\≥0}{\arg }min{||\mathrm{\Φ}w-s||}_2^2$

w^*Represent the optimal solution of w；

Step 2, sets up space structure consistency model:

Utilize Bayesian formula, obtain:

P(x|s)∝P(s|x)P(x)

Posterior probability density P (x | s) is proportional to data likelihood function P (s | x) and priori probability density The product of function P (x)；It is rewritten as afterwards:

$P\left(x\right|s)\∝1/(U^{In}+\underset{i}{\Σ}{\mathrm{\β}}_i{U}_i^{Ex})$

U^InIt is internal energy, U^ExBeing external energy, β is the hyper parameter of prior distribution；Posteriority The maximization of probability is converted into minimizing of total energy function:

$\arg \max P\left(x\right|s)\stackrel{\mathrm{\Δ}}{=}\mathrm{argmin}\{U^{In}+\underset{i}{\Σ}{\mathrm{\β}}_i{U}_i^{Ex}\}$

P (x | s) it is posterior probability density, U^InIt is internal energy, U^ExBeing external energy, β is first Test the hyper parameter of distribution；Wherein internal energy:

$U^{In}\∝{||S-S^{\′}||}_2^2={||S-\mathrm{\Θ}W||}_2^2$

S is to measure signal collection, and what S ' represented is signal to be estimated, and W is w coefficient sets, and Θ is Block diagonal matrix based on observing matrix；External energy:

$U^{Ex}=\underset{c\∈\mathrm{\Ω},\tilde{c}\∈H_c\⋐\mathrm{\Ω}}{\Σ}{||M(w_c-{\stackrel{\‾}{w}}_{\tilde{c}})||}_2$

U^ExRepresent external energy,Representing the arithmetic average of machine direction distribution function, M leads to Cross down-sampled direction v_tGo to train a dictionary base to obtain, w_cRepresent c Ω voxel of ∈ is Number；The estimation of W:

$\begin{array}{cc}\underset{w_c\≥0}{\mathrm{argmin}}\{{||s_c-{\mathrm{\Φw}}_c||}_2^2+{\mathrm{\β}}_c{||M(w_c-{\stackrel{\‾}{w}}_{\tilde{c}})||}_2\}& c\∈\mathrm{\Ω};\tilde{c}\∈H_c\⋐\mathrm{\Ω}\end{array}$

s_cIt is to measure signal coefficient in voxel c, w_cIt it is dictionary coefficient in voxel c；For Obtain the structure of voxel decussating fibers, define global cost function:

$\underset{W\≥0}{\arg \min }\{{||S-\mathrm{\Θ}W||}_2^2+\underset{c\∈\mathrm{\Ω}}{\Σ}(\underset{\tilde{c}\∈H_c\⋐\mathrm{\Ω}}{\Σ}{\mathrm{\σ}}_{c,\tilde{c}}{||w_c-w_{\tilde{c}}||}_2+{\mathrm{\β}}_1||Q\left(w_c-{\tilde{w}}_{\tilde{c}}\right)||+J\{w_c|c\∈\mathrm{\Ω}\})\}$

By calculating w_cAnd w_cThe coefficient COS distance of surrounding neighbors obtains, β₁It is artificial fixed One parameter of justice, Q is by training dictionary to obtain on basic function；Space structure consistency model Local linear approximate evaluation:

$\begin{array}{c}(w_c^{(t+1)},{\mathrm{\ξ}}_c^{(t+1)})=\underset{{\mathrm{\ξ}}_c^{(t+1)},w_c^{(t+1)}\≥0}{\arg \min }\{{||s_c-{\mathrm{\Φw}}_c^{\left(t\right)}||}_2^2+\underset{\tilde{c}\∈H_c\⋐\mathrm{\Ω}}{\Σ}{\mathrm{\σ}}_{c,\tilde{c}}{||w_c^{\left(t\right)}-w_{\tilde{c}}||}_2\\ +{\mathrm{\β}}_1{||Q(w_c^{\left(t\right)}-{\tilde{w}}_{\tilde{c}})||}_2+{\mathrm{\β}}_2{||{\mathrm{\ξ}}_c^{\left(t\right)}||}_1+{\mathrm{\δ}}_{\mathrm{\ξ}}{||w_c^{\left(t\right)}-{\mathrm{\ξ}}_c^{\left(t\right)}+p_{c,\mathrm{\ξ}}^{\left(t\right)}||}_2\}\\ p_{c,\mathrm{\ξ}}^{(t+1)}=p_{c,\mathrm{\ξ}}^{\left(t\right)}+w_c^{(t+1)}-{\mathrm{\ξ}}_c^{(t+1)}\end{array}$

T is iteration index, δ_ξIt is predefined aiding constant,It is that the extension of the t time iteration is drawn Ge Lang multiplier vector, above formula is divided into two parts:

$\begin{array}{c}{w}_c^{(t+1)}=\underset{w_c^{(t+1)}\≥0}{\arg \min }\{{||s_c-{\mathrm{\Φw}}_c^{\left(t\right)}||}_2^2+\underset{\tilde{c}\∈H_c\⋐\mathrm{\Ω}}{\Σ}{\mathrm{\σ}}_{c,\tilde{c}}{||w_c^{\left(t\right)}-w_{\tilde{c}}||}_2\\ +{\mathrm{\β}}_1{||Q(w_c^{\left(t\right)}-{\tilde{w}}_{\tilde{c}})||}_2+{\mathrm{\δ}}_{\mathrm{\ξ}}{||w_c^{\left(t\right)}-{\mathrm{\ξ}}_c^{\left(t\right)}+p_{c,\mathrm{\ξ}}^{\left(t\right)}||}_2\}\end{array}$

${\mathrm{\ξ}}_c^{(t+1)}=\underset{{\mathrm{\ξ}}_c^{(t+1)}\≥0}{\mathrm{argmin}}\{{\mathrm{\δ}}_{\mathrm{\ξ}}{||w_c^{\left(t\right)}-{\mathrm{\ξ}}_c^{\left(t\right)}+p_{c,\mathrm{\ξ}}^{\left(t\right)}||}_2+{\mathrm{\β}}_2{||{\mathrm{\ξ}}_c^{\left(t\right)}||}_1\}$

It is optimized for a separable space territory, solves with strengthening Lagrangian method, obtain:

$\begin{array}{c}{w}_c^{(t+1)}={({\mathrm{\Φ}}^T\mathrm{\Φ}+\underset{\tilde{c}\∈H_c\⋐\mathrm{\Ω}}{\Σ}{\mathrm{\σ}}_{c,\tilde{c}}I+{\mathrm{\β}}_c{Q}^{T}Q+{\mathrm{\δ}}_{\mathrm{\ξ}}I)}^{-1}\·\\ {({\mathrm{\Φ}}^T{s}_c+\underset{\tilde{c}\∈H_c\⋐\mathrm{\Ω}}{\Σ}{\mathrm{\σ}}_{c,\tilde{c}}{w}_{\tilde{c}}+{\mathrm{\β}}_1{Q}^{T}Q{\tilde{w}}_{\tilde{c}}+{\mathrm{\δ}}_{\mathrm{\ξ}}({\mathrm{\ξ}}_c^{\left(t\right)}-p_{c,\mathrm{\ξ}}^{\left(t\right)}\left)I\right)}_{w_c^{(t+1)}\≥0}\end{array}$

Wherein I representation unit matrix；

Based on separable space territoryMathematics soft MATLAB emulation is used to simulate FOD value Distribution, obtain the principal direction of fiber by the extreme point in search FOD value.

Claims (1)

1. one kind based on space structure conforming brain fiber microstructure reconstructing method, it is characterised in that: this reconstruct side Method comprises the steps:

Step one, uses bivalve basic function based on sphere dictionary deconvolution:

At rotating vectorAnd center vectorOn fiber probability-distribution function f (v | u) be referred to as FODF, wherein m_vRepresent the Spatial Dimension of rotating vector group, m_uRepresent is the Spatial Dimension of center vector group, By sphere the Method of Deconvolution, Fiber morphology structure is described as the convolution of kernel function, at diffusion gradient g ∈ S²On Measurement signal s (g | u):

$s\left(g\right|u)={\∫}_{S^2}r(g,v)f\left(v\right|u)d\mathrm{\μ}\left(v\right)$

Wherein (g, v) represents kernel function to r, and μ (v) is at S²On Haar estimate；For the sake of Fang Bian, and s (g | u) signal Can be represented by the discrete sampling μ being evenly distributed in unit sphere:

$s\left(g\right|u)=\underset{i=1}{\overset{m_u}{\Σ}}r(g,v)f\left(v\right|u_i)$

Definition fiber receptance function:

$r(g,v)=e^{-l{\left({\mathrm{gv}}^T\right)}^2}$

Wherein It is the diagonal matrix characterizing diffusion tensor D, λ_maxIt is that this is right Principal eigenvector in angular moment battle array, l represents that diffusion-sensitive coefficient and anisotropy interact to signal attenuation Influence degree；The estimation of fiber receptance function on many shells sampling plan is expanded with this:

$r(g,v)=\{r(g_{b_J},v)=e^{-l_{b_J}{\left(g_{b_J}{v}^T\right)}^2}|J=1,...,q\}$

Wherein J=1 ..., q refers to j-th spherical shell,Correspondence is distributed in b_JGradient vector collection g, l on shell_bJTable Show the l that the bJ diffusion-sensitive coefficient under j-th spherical shell is referred to；Thus push away to obtain a kind of new machine direction Distribution function form:

$f\left(v\right|u)=\underset{i=1}{\overset{m_u}{\Σ}}{w}_{i}d(v,u_i),u_i\∈u$

d(v,u_i) represent Vector Groups v and u_iOne group of directional spreding basic function the most complete under individual direction, w_iIt is The relative weighting of basic function,It is w_iThe number of middle nonzero element；Apparent diffusion coefficient can be by approximate evaluation:

g^TDg≈λ_maxg^T(vv^T) g=λ_maxcos²θ(g,v)

(g, v) represents included angle cosine function to θ, and D is to represent by all of d (v, u_i) one group of basic function word constituting Allusion quotation；In order to describe the dictionary base distribution under spherical coordinate, we are by d (v, the u under original rectangular coordinate system_i) describe For under spherical coordinate

It is the minimum angle cosine on discrete set, κ₁Machine direction distribution function is normalized to unit ball, κ₂It is For adjusting the parameter of peak value, τ represents even power；Basic function coefficient in machine direction distribution function (fODF) The estimation of w:

$\underset{w}{\arg min}\left|\right|\mathrm{\Φ}w-s|{|}_2^2$

Φ is observing matrix, and s is signal vector；For avoiding pseudo-peak and complicated algorithm and the big rule of higher order The calculating of mould ill-condition problem, the direct estimation being tried to achieve basic function coefficient w by non-negative least square method:

$\underset{w\≥0}{w^{*}\←\arg min}\left|\right|\mathrm{\Φ}w-s|{|}_2^2$

w^*Represent the optimal solution of w；

Step 2, sets up space structure consistency model:

Utilize Bayesian formula, obtain:

P(x|s)∝P(s|x)P(x)

Posterior probability density P (x | s) is proportional to data likelihood function P (s | x) and priori probability density function P (x) Product；It is rewritten as afterwards:

$P\left(x\right|s)\∝1/(U^{In}+\underset{i}{\Σ}{\mathrm{\β}}_i{U}_i^{Ex})$

U^InIt is internal energy, U^ExBeing external energy, β is the hyper parameter of prior distribution；Posterior probability is Bigization is converted into minimizing of total energy function:

$\arg \max P\left(x\right|s)\stackrel{\mathrm{\Δ}}{=}\mathrm{argmin}\{U^{In}+\underset{i}{\Σ}{\mathrm{\β}}_i{U}_i^{Ex}\}$

P (x | s) it is posterior probability density, U^InIt is internal energy, U^ExBeing external energy, β is prior distribution Hyper parameter；Wherein internal energy:

$U^{In}\∝\left|\right|S-S^{\′}|{|}_2^2=||S-\mathrm{\Θ}W|{|}_2^2$

S is to measure signal collection, and what S ' represented is signal to be estimated, and W is w coefficient sets, and Θ is based on observation The block diagonal matrix of matrix；External energy:

$U^{Ex}=\underset{c\∈\mathrm{\Ω},\tilde{c}\∈H_c\⋐\mathrm{\Ω}}{\Σ}\left|\right|M(w_c-{\stackrel{\‾}{w}}_{\tilde{c}})|{|}_2$

U^ExRepresent external energy,Representing the arithmetic average of machine direction distribution function, M is by down-sampled Direction v_tGo to train a dictionary base to obtain, w_cRepresent the coefficient of c Ω voxel of ∈；The estimation of W:

$\underset{w_c\≥0}{\arg min}\left\{\right||s_c-{\mathrm{\Φw}}_c|{|}_2^2+{\mathrm{\β}}_c\left|\right|M(w_c-{\stackrel{\‾}{w}}_{\tilde{c}})\left|{|}_2\right\},c\∈\mathrm{\Ω};\tilde{c}\∈H_c\⋐\mathrm{\Ω}$

s_cIt is to measure signal coefficient in voxel c, w_cIt it is dictionary coefficient in voxel c；In order to obtain body The structure of element decussating fibers, definition global cost function:

$\underset{W\≥0}{\arg \min }\left\{\left|\right|S-\mathrm{\Θ}W|{|}_2^2+\underset{c\∈\mathrm{\Ω}}{\mathrm{\Σ}}\left(\underset{\tilde{c}\∈H_c\⋐\mathrm{\Ω}}{\mathrm{\Σ}}{\mathrm{\σ}}_{c,\tilde{c}}\left|\right|w_c-w_{\tilde{c}}|{|}_2+{\mathrm{\β}}_1|\left|Q\left(w_c-{\tilde{w}}_{\tilde{c}}\right)\right||+J\left\{w_c|c\∈\mathrm{\Ω}\right\}\right)\right\}$

By calculating w_cAnd w_cThe coefficient COS distance of surrounding neighbors obtains, β₁It it is an artificially defined ginseng Number, Q is by training dictionary to obtain on basic function；The local linear approximate evaluation of space structure consistency model:

$\begin{array}{c}\left(w_c^{\left(t+1\right)},{\mathrm{\ξ}}_c^{\left(t+1\right)}\right)=\underset{{\mathrm{\ξ}}_c^{\left(t+1\right)},w_c^{\left(t+1\right)}\≥0}{\arg \min }\left\{\right||s_c-{\mathrm{\Φw}}_c^{\left(t\right)}|{|}_2^2+\underset{\tilde{c}\∈H_c\⋐\mathrm{\Ω}}{\mathrm{\Σ}}{\mathrm{\σ}}_{c,\tilde{c}}\left|\right|w_c^{\left(t\right)}-w_{\tilde{c}}|{|}_2\\ +{\mathrm{\β}}_1\left|\right|Q\left(w_c^{\left(t\right)}-{\tilde{w}}_{\tilde{c}}\right)|{|}_2+{\mathrm{\β}}_2|\left|{\mathrm{\ξ}}_c^{\left(t\right)}\right|{|}_1+{\mathrm{\δ}}_{\mathrm{\ξ}}\left|\right|w_c^{\left(t\right)}-{\mathrm{\ξ}}_c^{\left(t\right)}+p_{c,\mathrm{\ξ}}^{\left(t\right)}\left|{|}_2\right\}\\ p_{c,\mathrm{\ξ}}^{\left(t+1\right)}=p_{c,\mathrm{\ξ}}^{\left(t\right)}+w_c^{\left(t+1\right)}-{\mathrm{\ξ}}_c^{\left(t+1\right)}\end{array}$

T is iteration index, δ_ξIt is predefined aiding constant,It it is the extension Lagrange multiplier of the t time iteration Vector, above formula is divided into two parts:

$\begin{array}{c}{w}_c^{\left(t+1\right)}=\underset{w_c^{\left(t+1\right)}\≥0}{\arg \min }\left\{\right||s_c-{\mathrm{\Φw}}_c^{\left(t\right)}|{|}_2^2+\underset{\tilde{c}\∈H_c\⋐\mathrm{\Ω}}{\mathrm{\Σ}}{\mathrm{\σ}}_{c,\tilde{c}}\left|\right|w_c^{\left(t\right)}-w_{\tilde{c}}|{|}_2\\ +{\mathrm{\β}}_1\left|\right|Q\left(w_c^{\left(t\right)}-{\tilde{w}}_{\tilde{c}}\right)|{|}_2+{\mathrm{\δ}}_{\mathrm{\ξ}}||w_c^{\left(t\right)}-{\mathrm{\ξ}}_c^{\left(t\right)}+p_{c,\mathrm{\ξ}}^{\left(t\right)}|{|}_2\}\end{array}$

${\mathrm{\ξ}}_c^{(t+1)}=\underset{{\mathrm{\ξ}}_c^{\left(t+1\right)}\≥0}{\mathrm{argmin}}\left\{{\mathrm{\δ}}_{\mathrm{\ξ}}\right||w_c^{\left(t\right)}-{\mathrm{\ξ}}_c^{\left(t\right)}+p_{c,\mathrm{\ξ}}^{\left(t\right)}|{|}_2+{\mathrm{\β}}_2\left|\right|{\mathrm{\ξ}}_c^{\left(t\right)}\left|{|}_1\right\}$

It is optimized for a separable space territory, solves with strengthening Lagrangian method, obtain:

$\begin{array}{c}{w}_c^{\left(t+1\right)}={\left({\mathrm{\Φ}}^T\mathrm{\Φ}+\underset{\tilde{c}\∈H_c\⋐\mathrm{\Ω}}{\mathrm{\Σ}}{\mathrm{\σ}}_{c,\tilde{c}}I+{\mathrm{\β}}_c{Q}^{T}Q+{\mathrm{\δ}}_{\mathrm{\ξ}}I\right)}^{-1}\·\\ {\left({\mathrm{\Φ}}^T{s}_c+\underset{\tilde{c}\∈H_c\⋐\mathrm{\Ω}}{\mathrm{\Σ}}{\mathrm{\σ}}_{c,\tilde{c}}{w}_{\tilde{c}}+{\mathrm{\β}}_1{Q}^{T}Q{\tilde{w}}_{\tilde{c}}+{\mathrm{\δ}}_{\mathrm{\ξ}}\left({\mathrm{\ξ}}_c^{\left(t\right)}-p_{c,\mathrm{\ξ}}^{\left(t\right)}\right)I\right)}_{w_c^{\left(t+1\right)}\≥0}\end{array}$

Wherein I representation unit matrix；

Based on separable space territoryMathematics soft MATLAB emulation is used to simulate the distribution of FOD value, The principal direction of fiber is obtained by the extreme point in search FOD value.