CN103970929A - High-order diffusion tensor mixture sparse imaging method for alba fiber tracking - Google Patents

Abstract

High-order diffusion tensor for brain white matter integrity tracking mixes sparse imaging method, comprising: reads brain MR data, obtains the magnetic resonance signal for applying gradient direction The magnetic resonance signal of gradient direction is not applied And gradient direction data ROI region needed for choosing calculates diffusive attenuation signal S (g)/S0; Diffusive attenuation signal in each voxel in ROI is modeled as to the ellipsoid distributed model of spreading morphology one by one; The estimation of machine direction is tensor coefficient to be obtained by calculation to obtain spread function, then calculate the spread function value of each sampling to obtain diffusion model.

Description

The high-order diffusion tensor of following the tracks of for brain white matter integrity mixes sparse formation method

(1) technical field

The present invention relates to the fields such as image processing, medical imaging, computing method, mathematics, three-dimensional reconstruction, Nervous System Anatomy, especially this kind of high-order diffusion tensor mixes sparse imaging.

(2) background technology

It is a category information medical technology that obtains cerebral white matter region fiber orientation that brain white matter integrity is followed the tracks of, and it walks always to estimate the possible path of fiber by following the trail of local tensors.For now, brain white matter integrity tracking is unique method that obtains to nothing wound brain white matter integrity trend with kind of energy in live body.First tractography carries out voxel modeling to original DW-MRI data, obtains the fiber orientation in each voxel, forms and has the fibre space microstructure of Anatomical significance, and then utilize fiber tracking algorithm to connect the machine direction of appointed area.Along with the sampling precision of magnetic resonance diffusion signal and the raising to fiber tracking accuracy requirement, machine direction estimation problem to solve scale increasing, make the stable high resolving power fiber identification difficulty that obtains, this has greatly hindered the application of this technology in clinical medicine.

(3) summary of the invention

In order to overcome the deficiencies such as existing method angular resolution and counting yield are low, the present invention proposes a kind of machine direction distribution estimation method taking high order tensor as the low sample number of high angular resolution high-level efficiency of guiding.

The technical solution adopted in the present invention is as follows:

The high-order diffusion tensor of following the tracks of for brain white matter integrity mixes sparse formation method, it is characterized in that: the sparse formation method of described mixing comprises the following steps:

(1) read brain MR data, obtain the magnetic resonance signal that applies gradient direction do not apply the magnetic resonance signal of gradient direction and gradient direction data choose required ROI region, calculate diffusive attenuation signal S (g)/S ₀.

(2) the diffusive attenuation signal in the each voxel in ROI is modeled as one by one to the ellipsoid distributed model that spreads form, its modeling process is as follows:

2.1) voxel microstructure modeling scheme:

By diffusive attenuation signal S (g)/S ₀be assumed to be wall scroll fiber signals response function R (v, g) and spread function D (v) in Spherical Surface S ²on convolution:

$S\left(g\right)/S_0=R(v,g)\⊗D\left(v\right)={\∫}_{s^2}R(v,g)D\left(v\right)\mathrm{dv}$

Wherein response function while only containing a fiber in expression voxel, represent that machine direction distribution function also claims spread function, d _ijrepresent tensor coefficient, f _j(v) represent j monomial.

2.2) structure mathematical model is as follows:

Because data acquisition is not idealized, therefore often with certain noise, in order to overcome as much as possible the impact of noise, conventionally energy function is minimized to reduce error.If diffusion-weighted magnetic resonance signals has m diffusion gradient pulse direction, f (v) tries to achieve by the energy function minimizing below

$E=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(S\left(g_i\right)/S_0-\underset{S_2}{\∫}R(v,g_i)D\left(v\right)\mathrm{dv})}^2=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(S\left(g_i\right)/S_0-\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\λ}}_j\underset{S_2}{\∫}{e}^{-\mathrm{\μ}{\left({g_i}^T{v}_j\right)}^2}{f}_j\left(v\right)\mathrm{dv})}^2$

Obtain tensor coefficient lambda in order to solve energy function _j, be therefore converted into linear problem

y=Af+ξ

Wherein ξ is noise; A is R ^{n × m}matrix, each element in A f=[λ ₁, λ ₂.. λ _m] ^t;

Y is n dimensional vector, each element y _i=S (g _i)/S ₀

(3) estimation of machine direction is to obtain spread function by calculating tensor coefficient, then the spread function value of calculating each sampling obtains diffusion model.The computing method of tensor coefficient comprise the following steps

3.1) 321 discrete points of uniform sampling on unit hemisphere face, obtain this 321 vector of unit length, calculate the value M of wall scroll fiber response function ₁, set a lower exponent number order evaluator matrix M ₂.Calculate M ₁with M ₂convolution as 2.2) in A.

3.2) can, by non-sparse demapping in sparse territory, make its rarefaction representation: f=Ψ x by some sparse conversion

Make Φ=A Ψ, be called sensing matrix.Therefore we are rewritten as the problems referred to above:

y=Af=AΨx=Φx

3.3) mixed weighting Corresponding Sparse Algorithm solves this inverse problem, and its process comprises the following steps:

Step 3.3.1, we solve following constrained optimization problem and come initialization search volume, to obtain the solution of non-negative initial x:

$x^{\left(0\right)\←\arg \min {\left|\right|{\mathrm{\Φ}}^{0}x-y\left|\right|}_2^2,s,t-{\mathrm{\Φ}}^{0}x\≥0\mathrm{andx}\≥0}$

Step 3.3.2, we train and obtain regularization matrix L by low order mixed weighting rarefaction method

$x^{(t+1)}\←\arg \min {\left|\right|{\mathrm{\Φ}}^0{x}^{\left(t\right)}-y\left|\right|}_2^2+\mathrm{\λ}(\mathrm{\α}{\left|\right|{\mathrm{\ω}}^{\left(t\right)}{x}^{\left(t\right)}\left|\right|}_1+(1+\mathrm{\α}){\left|\right|{\mathrm{Lx}}^{\left(t\right)}\left|\right|}_2^2)s,t\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{x}_i^{\left(t\right)}=1\mathrm{andx}\≥0$

Wherein λ is constrained parameters.α is used for controlling two balances between penalty function condition.Wherein l ₁norm is used for promoting approaching of sparse solution, and l ₂norm is used for promoting the flatness of solution.

L is regular matrix, solves by the following method

$L_{m^{,},n^{,}}=\left\{\begin{array}{cc}{D}_{m^{,},n^{,}}& \left|\mathrm{\&mu;}\right|<\mathrm{\&tau;}\stackrel{\&OverBar;}{D}\\ 0& |\mathrm{\&mu;}\&GreaterEqual;\mathrm{\&tau;}\stackrel{\&OverBar;}{D}\end{array}\right.$

Wherein τ is threshold parameter, is set to 0.1, for the mean value of D, μ=Dx;

Step 3.3.3, upgrades weights:

${\mathrm{\&omega;}}_i^{(t+1)}=\frac{\mathrm{sgn}(x^{\left(t\right)}-\mathrm{\&beta;})}{\left|x_i^{\left(t\right)}\right|+\mathrm{\&delta;}}$

When the solution of twice iteration in front and back satisfies condition time termination of iterations, otherwise return to step 2;

Step 3.3.4, the compressed sensing matrix that the L that we obtain by above training and high-order produce solves final x:

x _finally=[(Φ ¹) ^T(Φ ¹)+λ(L ¹） ^T(L ¹）] ^-1[(Φ ¹) ^Ty]

Step 3.3.5, obtains tensor coefficient by the x substitution f=Ψ x trying to achieve

3.4) high order tensor coefficient is used for to matching diffusion profile, obtains fiber diffusion model, search extreme value is calculated machine direction.Its process comprises step:

Step 3.4.1, to dodecahedron carry out four times fractal, obtain 2562 sampled points that approach with unit sphere sampling, obtain 2562 vector of unit length.The tensor coefficient obtaining by previous step just can obtain final direction distribution function

$D\left(v\right)=\underset{j=1}{\overset{m}{\mathrm{\&Sigma;}}}{\mathrm{\&lambda;}}_j{f}_j\left(v\right)$

Step 3.4.2, obtains q extreme point by population local extremum searching method, and each extreme value neighborhood of a point is searched for respectively to sparse field point (t sparse point altogether).

Step 3.4.3, the direction distribution function value of calculating respectively these 2562 sampled points by the direction distribution function obtaining above also claims FOD value, and in fact most of FOD is 0, its actual calculated amount is only the FOD value of t sparse point.The distribution of the FOD value that use MATLAB Software simulation calculation goes out, obtains machine direction by the extreme point in search FOD value.

Further, described step 2.1) in spread function D (v) use higher order polynomial as spread function.

Described step 2.2) in constructed mathematical model use the method for spherical convolution minimization of energy function:

$E=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{(S\left(g_i\right)/S_0-\underset{S_2}{\&Integral;}R(v,g_i)D\left(v\right)\mathrm{dv})}^2=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{(S\left(g_i\right)/S_0-\underset{j=1}{\overset{m}{\mathrm{\&Sigma;}}}{\mathrm{\&lambda;}}_j\underset{S_2}{\&Integral;}{e}^{-\mathrm{\&mu;}{\left({g_i}^T{v}_j\right)}^2}{f}_j\left(v\right)\mathrm{dv})}^2$

$D\left(v\right)=\underset{i=0}{\overset{l}{\mathrm{\&Sigma;}}}\underset{j=0}{\overset{l-i}{\mathrm{\&Sigma;}}}{d}_{\mathrm{ij}}{v}_1^i{v}_2^j{v}_3^{l-i-j}$

Initialization search volume in described step 3.3.1 is by two groups of linear non-negative space-Φ ⁰initialization x in x>=0 and x>=0 scope.

In described step 3.3.2, train regularization matrix to use punishment weighting l ₁and l ₂penalty function set up the metastable convergence of the sparse property of prioritization scheme balance and flatness.

Regular matrix L in described step 3.3.2 adjusts threshold parameter μ by the mean value of each iterative computation tensor D to search for the point that is less than threshold value μ as adjusting matrix:

$L_{m^{,},n^{,}}=\left\{\begin{array}{cc}{D}_{m^{,},n^{,}}& \left|\mathrm{\&mu;}\right|<\mathrm{\&tau;}\stackrel{\&OverBar;}{D}\\ 0& |\mathrm{\&mu;}\&GreaterEqual;\mathrm{\&tau;}\stackrel{\&OverBar;}{D}\end{array}\right.$

Wherein τ is threshold parameter, for the mean value of D, μ=Dx.

Weighted optimization method method in described step 3.3.3 is the mode by adding sparse weights ω:

${\mathrm{\&omega;}}_i^{(t+1)}=\frac{\mathrm{sgn}(x^{\left(t\right)}-\mathrm{\&beta;})}{\left|x_i^{\left(t\right)}\right|+\mathrm{\&delta;}}$

In described step 3.4.2 and 3.4.3, the FOD value method of Local Extremum is to search for q Local Extremum by particle swarm optimization algorithm, near the sparse point of t neighborhood this q Local Extremum of sampling again, obtain the optimal information that concentrates on machine direction, t the sparse value by the optimal information along machine direction is used for calculating FOD value, avoids redundancy phenomenon to improve counting yield.

In described step 3.4.3, machine direction is the distribution that uses the FOD value that MATLAB Software simulation calculation goes out, and utilizes extremum search method to obtain the direction of fiber.

The present invention relates to the theory of compressed sensing and sparse imaging, therefore compare classic method, computing velocity is fast, and imaging angular resolution is high, and calculating sample number can greatly reduce.It seems according to the author's experiment, income effect of the present invention is current this field best effects.

(4) brief description of the drawings

Fig. 1 is of the present invention

Fig. 2 is simulated data result figure of the present invention.Wherein, simulated data is produced by following formula:

$S\left(g\right)=s_0\underset{i=1}{\overset{2}{\mathrm{\&Sigma;}}}{f}_i{S}_0{e}^{-{\mathrm{bg}}^T\mathrm{Dg}}$

Wherein f represents that i is with the shared ratio of fiber, ${\mathrm{\&Sigma;}}_{i=1}^2{f}_i=1,f_1=0.5,f_2=0.5,S_0=1,b=3000s/{\mathrm{mm}}^2,$ The eigenwert of diffusion tensor D is: λ ₁=1.8 × 10 ^-3mm ²s, λ ₂=0.3 × 10 ^-3mm ²s, λ ₃=0.3 × 10 ^-3mm ²s.Order=18 in experiment, 81 equally distributed diffusion-weighted magnetic resonance imaging directions in hemisphere face, hemisphere face sampling number is 321, in figure, the first row represents angle, the second line display machine direction, the third line represents imaging model, the direction (obtaining by calculating diffuse peak) of two fibers of black line signal.

Fig. 3 is actual clinical effect data figure of the present invention, and has drawn the roughly direction of fiber orientation.Real data is from Harvard University's hospital attached to a medical college (Brigham and Women ' s Hospital, Brockton VA Hospital, McLean Hospital), the brain data of utilizing 3-T GE system and double echo plane imaging sequence to extract from true human brain, data acquisition parameters is: TR=17000ms, TE=78ms. voxel amount is 144 × 144 × 85, becoming image field is 85 axial slices that 24cm. is parallel to AC-PC line, every layer thickness 1.7mm. is from 51 different gradient direction scan-datas, diffusion-sensitive coefficient b=900s/mm2, the scan-data of 8 b=0.

(5) concrete implementation step

With reference to accompanying drawing:

The high-order diffusion tensor of following the tracks of for brain white matter integrity of the present invention mixes sparse formation method, it is characterized in that: the sparse formation method of described mixing comprises the following steps:

(1) read brain MR data, obtain the magnetic resonance signal that applies gradient direction do not apply the magnetic resonance signal of gradient direction and gradient direction data choose required ROI region, calculate diffusive attenuation signal S (g)/S ₀.

(2) the diffusive attenuation signal in the each voxel in ROI is modeled as one by one to the ellipsoid distributed model that spreads form, its modeling process is as follows:

2.1) voxel microstructure modeling scheme:

By diffusive attenuation signal S (g)/S ₀be assumed to be wall scroll fiber signals response function R (v, g) and spread function D (v) in Spherical Surface S ²on convolution:

$S\left(g\right)/S_0=R(v,g)\&CircleTimes;D\left(v\right)={\&Integral;}_{s^2}R(v,g)D\left(v\right)\mathrm{dv}$

Wherein response function while only containing a fiber in expression voxel,

represent that machine direction distribution function also claims spread function, d _ijrepresent tensor coefficient, f _j(v) represent j monomial.

2.2) structure mathematical model is as follows:

Because data acquisition is not idealized, therefore often with certain noise, in order to overcome as much as possible the impact of noise, conventionally energy function is minimized to reduce error.If diffusion-weighted magnetic resonance signals has m diffusion gradient pulse direction, f (v) tries to achieve by the energy function minimizing below

$E=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{(S\left(g_i\right)/S_0-\underset{S_2}{\&Integral;}R(v,g_i)D\left(v\right)\mathrm{dv})}^2=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{(S\left(g_i\right)/S_0-\underset{j=1}{\overset{m}{\mathrm{\&Sigma;}}}{\mathrm{\&lambda;}}_j\underset{S_2}{\&Integral;}{e}^{-\mathrm{\&mu;}{\left({g_i}^T{v}_j\right)}^2}{f}_j\left(v\right)\mathrm{dv})}^2$

Obtain tensor coefficient lambda in order to solve energy function _j, be therefore converted into linear problem

y=Af+ξ

Wherein ξ is noise; A is R ^{n × m}matrix, each element in A f=[λ ₁, λ ₂.. λ _m] ^t;

Y is n dimensional vector, each element y _i=S (g _i)/S ₀

(3) estimation of machine direction is to obtain spread function by calculating tensor coefficient, then the spread function value of calculating each sampling obtains diffusion model.The computing method of tensor coefficient comprise the following steps

3.1) 321 discrete points of uniform sampling on unit hemisphere face, obtain this 321 vector of unit length, calculate the value M of wall scroll fiber response function ₁, set a lower exponent number order evaluator matrix M ₂.Calculate M ₁with M ₂convolution as 2.2) in A.

3.2) can, by non-sparse demapping in sparse territory, make its rarefaction representation: f=Ψ x by some sparse conversion

Make Φ=A Ψ, be called sensing matrix.Therefore we are rewritten as the problems referred to above:

y=Af=AΨx=Φx

3.3) mixed weighting Corresponding Sparse Algorithm solves this inverse problem, and its process comprises the following steps:

Step 3.3.1, we solve following constrained optimization problem and come initialization search volume, to obtain the solution of non-negative initial x:

$x^{\left(0\right)\&LeftArrow;\arg \min {\left|\right|{\mathrm{\&Phi;}}^{0}x-y\left|\right|}_2^2,s,t-{\mathrm{\&Phi;}}^{0}x\&GreaterEqual;0\mathrm{andx}\&GreaterEqual;0}$

Step 3.3.2, we train and obtain regularization matrix L by low order mixed weighting rarefaction method

$x^{(t+1)}\&LeftArrow;\arg \min {\left|\right|{\mathrm{\&Phi;}}^0{x}^{\left(t\right)}-y\left|\right|}_2^2+\mathrm{\&lambda;}(\mathrm{\&alpha;}{\left|\right|{\mathrm{\&omega;}}^{\left(t\right)}{x}^{\left(t\right)}\left|\right|}_1+(1+\mathrm{\&alpha;}){\left|\right|{\mathrm{Lx}}^{\left(t\right)}\left|\right|}_2^2)s,t\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{x}_i^{\left(t\right)}=1\mathrm{andx}\&GreaterEqual;0$

Wherein λ is constrained parameters.α is used for controlling two balances between penalty function condition.Wherein l ₁norm is used for promoting approaching of sparse solution, and l ₂norm is used for promoting the flatness of solution.

L is regular matrix, solves by the following method

$L_{m^{,},n^{,}}=\left\{\begin{array}{cc}{D}_{m^{,},n^{,}}& \left|\mathrm{\&mu;}\right|<\mathrm{\&tau;}\stackrel{\&OverBar;}{D}\\ 0& |\mathrm{\&mu;}\&GreaterEqual;\mathrm{\&tau;}\stackrel{\&OverBar;}{D}\end{array}\right.$

Wherein τ is threshold parameter, is set to 0.1, for the mean value of D, μ=Dx;

Step 3.3.3, upgrades weights:

${\mathrm{\&omega;}}_i^{(t+1)}=\frac{\mathrm{sgn}(x^{\left(t\right)}-\mathrm{\&beta;})}{\left|x_i^{\left(t\right)}\right|+\mathrm{\&delta;}}$

When the solution of twice iteration in front and back satisfies condition time termination of iterations, otherwise return to step 2;

Step 3.3.4, the compressed sensing matrix that the L that we obtain by above training and high-order produce solves final x:

x _finally=[(Φ ¹) ^T(Φ ¹)+λ(L ¹） ^T(L ¹）] ^-1[(Φ ¹) ^Ty]

Step 3.3.5, obtains tensor coefficient by the x substitution f=Ψ x trying to achieve

3.4) high order tensor coefficient is used for to matching diffusion profile, obtains fiber diffusion model, search extreme value is calculated machine direction.Its process comprises step:

Step 3.4.1, to dodecahedron carry out four times fractal, obtain 2562 sampled points that approach with unit sphere sampling, obtain 2562 vector of unit length.The tensor coefficient obtaining by previous step just can obtain final direction distribution function

$D\left(v\right)=\underset{j=1}{\overset{m}{\mathrm{\&Sigma;}}}{\mathrm{\&lambda;}}_j{f}_j\left(v\right)$

Step 3.4.2, obtains q extreme point by population local extremum searching method, and each extreme value neighborhood of a point is searched for respectively to sparse field point (t sparse point altogether).

Step 3.4.3, the direction distribution function value of calculating respectively these 2562 sampled points by the direction distribution function obtaining above also claims FOD value, and in fact most of FOD is 0, its actual calculated amount is only the FOD value of t sparse point.The distribution of the FOD value that use MATLAB Software simulation calculation goes out, obtains machine direction by the extreme point in search FOD value.

Further, described step 2.1) in spread function D (v) use higher order polynomial as spread function.

Described step 2.2) in constructed mathematical model use the method for spherical convolution minimization of energy function:

$E=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{(S\left(g_i\right)/S_0-\underset{S_2}{\&Integral;}R(v,g_i)D\left(v\right)\mathrm{dv})}^2=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{(S\left(g_i\right)/S_0-\underset{j=1}{\overset{m}{\mathrm{\&Sigma;}}}{\mathrm{\&lambda;}}_j\underset{S_2}{\&Integral;}{e}^{-\mathrm{\&mu;}{\left({g_i}^T{v}_j\right)}^2}{f}_j\left(v\right)\mathrm{dv})}^2$

$D\left(v\right)=\underset{i=0}{\overset{l}{\mathrm{\&Sigma;}}}\underset{j=0}{\overset{l-i}{\mathrm{\&Sigma;}}}{d}_{\mathrm{ij}}{v}_1^i{v}_2^j{v}_3^{l-i-j}$

Initialization search volume in described step 3.3.1 is by two groups of linear non-negative space-Φ ⁰initialization x in x>=0 and x>=0 scope.

In described step 3.3.2, train regularization matrix to use punishment weighting l ₁and l ₂penalty function set up the metastable convergence of the sparse property of prioritization scheme balance and flatness.

Regular matrix L in described step 3.3.2 adjusts threshold parameter μ by the mean value of each iterative computation tensor D to search for the point that is less than threshold value μ as adjusting matrix:

$L_{m^{,},n^{,}}=\left\{\begin{array}{cc}{D}_{m^{,},n^{,}}& \left|\mathrm{\&mu;}\right|<\mathrm{\&tau;}\stackrel{\&OverBar;}{D}\\ 0& |\mathrm{\&mu;}\&GreaterEqual;\mathrm{\&tau;}\stackrel{\&OverBar;}{D}\end{array}\right.$

Wherein τ is threshold parameter, for the mean value of D, μ=Dx.

Weighted optimization method method in described step 3.3.3 is the mode by adding sparse weights ω:

${\mathrm{\&omega;}}_i^{(t+1)}=\frac{\mathrm{sgn}(x^{\left(t\right)}-\mathrm{\&beta;})}{\left|x_i^{\left(t\right)}\right|+\mathrm{\&delta;}}$

In described step 3.4.2 and 3.4.3, the FOD value method of Local Extremum is to search for q Local Extremum by particle swarm optimization algorithm, near the sparse point of t neighborhood this q Local Extremum of sampling again, obtain the optimal information that concentrates on machine direction, t the sparse value by the optimal information along machine direction is used for calculating FOD value, avoids redundancy phenomenon to improve counting yield.

In described step 3.4.3, machine direction is the distribution that uses the FOD value that MATLAB Software simulation calculation goes out, and utilizes extremum search method to obtain the direction of fiber.

Claims (9)

1. the high-order diffusion tensor of following the tracks of for brain white matter integrity mixes sparse formation method, it is characterized in that: the sparse formation method of described mixing comprises the following steps:

(1) read brain MR data, obtain the magnetic resonance signal that applies gradient direction do not apply the magnetic resonance signal of gradient direction and gradient direction data choose required ROI region, calculate diffusive attenuation signal S (g)/S ₀.

(2) the diffusive attenuation signal in the each voxel in ROI is modeled as one by one to the ellipsoid distributed model that spreads form, its modeling process is as follows:

2.1) voxel microstructure modeling scheme:

By diffusive attenuation signal S (g)/S ₀be assumed to be wall scroll fiber signals response function R (v, g) and spread function D (v) in Spherical Surface S ²on convolution:

Wherein response function while only containing a fiber in expression voxel, represent that machine direction distribution function also claims spread function, d _ijrepresent tensor coefficient, f _j(v) represent j monomial.

2.2) structure mathematical model is as follows:

Because data acquisition is not idealized, therefore often with certain noise, in order to overcome as much as possible the impact of noise, conventionally energy function is minimized to reduce error.If diffusion-weighted magnetic resonance signals has m diffusion gradient pulse direction, f (v) tries to achieve by the energy function minimizing below

Obtain tensor coefficient lambda in order to solve energy function _j, be therefore converted into linear problem

y=Af+ξ

Wherein ξ is noise; A is R ^{n × m}matrix, each element in A f=[λ ₁, λ ₂.. λ _m] ^t; Y is n dimensional vector, each element y _i=S (g _i)/S ₀

(3) estimation of machine direction is to obtain spread function by calculating tensor coefficient, then the spread function value of calculating each sampling obtains diffusion model.The computing method of tensor coefficient comprise the following steps

3.1) 321 discrete points of uniform sampling on unit hemisphere face, obtain this 321 vector of unit length, calculate the value M of wall scroll fiber response function ₁, set a lower exponent number order evaluator matrix M ₂.Calculate M ₁with M ₂convolution as 2.2) in A.

3.2) can, by non-sparse demapping in sparse territory, make its rarefaction representation: f=Ψ x makes Φ=A Ψ, is called sensing matrix by some sparse conversion.Therefore we are rewritten as the problems referred to above:

y=Af=AΨx=Φx

3.3) mixed weighting Corresponding Sparse Algorithm solves this inverse problem, and its process comprises the following steps:

Step 3.3.1, we solve following constrained optimization problem and come initialization search volume, to obtain the solution of non-negative initial x:

Step 3.3.2, we train and obtain regularization matrix L by low order mixed weighting rarefaction method

Wherein λ is constrained parameters.α is used for controlling two balances between penalty function condition.Wherein l ₁norm is used for promoting approaching of sparse solution, and l ₂norm is used for promoting the flatness of solution.

L is regular matrix, solves by the following method

Wherein τ is threshold parameter, is set to 0.1, for the mean value of D, μ=Dx;

Step 3.3.3, upgrades weights:

When the solution of twice iteration in front and back satisfies condition time termination of iterations, otherwise return to step 2;

Step 3.3.4, the compressed sensing matrix that the L that we obtain by above training and high-order produce solves final x:

x _finally=[(Φ ¹) ^T(Φ ¹)+λ(L ¹） ^T(L ¹）] ^-1[(Φ ¹) ^Ty]

Step 3.3.5, obtains tensor coefficient by the x substitution f=Ψ x trying to achieve

3.4) high order tensor coefficient is used for to matching diffusion profile, obtains fiber diffusion model, search extreme value is calculated machine direction.Its process comprises step:

Step 3.4.1, to dodecahedron carry out four times fractal, obtain 2562 sampled points that approach with unit sphere sampling, obtain 2562 vector of unit length.The tensor coefficient obtaining by previous step just can obtain final direction distribution function

Step 3.4.2, obtains q extreme point by population local extremum searching method, and each extreme value neighborhood of a point is searched for respectively to sparse field point (t sparse point altogether).

Step 3.4.3, the direction distribution function value of calculating respectively these 2562 sampled points by the direction distribution function obtaining above also claims FOD value, and in fact most of FOD is 0, its actual calculated amount is only the FOD value of t sparse point.The distribution of the FOD value that use MATLAB Software simulation calculation goes out, obtains machine direction by the extreme point in search FOD value.

2. the method for claim 1, is characterized in that: described step 2.1) in spread function D (v) use higher order polynomial as spread function.

3. the method for claim 1, is characterized in that: described step 2.2) in the mathematical model of constructed machine direction use the method for spherical convolution minimization of energy function:

4. the method for claim 1, is characterized in that: the initialization search volume in described step 3.3.1 is by two groups of linear non-negative space-Φ ⁰initialization x in x>=0 and x>=0 scope.

5. the method for claim 1, is characterized in that: in described step 3.3.2, train regularization matrix to use punishment weighting l ₁and l ₂penalty function set up the metastable convergence of the sparse property of prioritization scheme balance and flatness.

6. the method for claim 1, is characterized in that: the regular matrix L in described step 3.3.2 adjusts threshold parameter μ by the mean value of each iterative computation tensor D to search for the point that is less than threshold value μ as adjusting matrix:

Wherein τ is threshold parameter, for the mean value of D, μ=Dx.

7. the method for claim 1, is characterized in that: the weighted optimization method method in described step 3.3.3 is the mode by adding sparse weights ω:

8. the method for claim 1, it is characterized in that: in described step 3.4.2 and 3.4.3, the FOD value method of Local Extremum is to search for q Local Extremum by particle swarm optimization algorithm, near the sparse point of t neighborhood this q Local Extremum of sampling again, obtain the optimal information that concentrates on machine direction, t the sparse value by the optimal information along machine direction is used for calculating FOD value, avoids redundancy phenomenon to improve counting yield.

9. the method for claim 1, is characterized in that: in described step 3.4.3, machine direction is the distribution that uses the FOD value that MATLAB Software simulation calculation goes out, and utilizes extremum search method to obtain the direction of fiber.