CN105303577A - Cerebral white matter fiber imaging method - Google Patents

Abstract

The invention discloses a method for eliminating the partial volume effect of cerebral white matter fiber imaging by a mixed response kernel function. The method comprises the following steps: reading cerebral magnetic resonance data, obtaining a magnetic resonance signal S(g) of which the gradient applying direction is g, a magnetic resonance signal S0 without a gradient direction and gradient direction data, selecting a required interested area, and calculating a diffusion attenuation signal S(g)/S0 of the area; utilizing a Richardson-Lucy iterative algorithm to carrying out modeling on the diffusion attenuation signal S(g)/S0 of each voxel in the interested area into models with ellipsoidal distribution one by one, adding total variation-based regularization to eliminate noise influence, obtaining a diffusion direction distribution function f(v) through the Richardson-Lucy iterative algorithm, and separating an isotropous diffusion value f(m+1) and an anisotropic diffusion fiber distribution function fod from the f(v). The invention relates to a probability theory and a total variation-based regularization concept. Compared with the traditional method, the method disclosed by the invention has the advantages of being high in calculation speed, high in imaging angle resolution and favorable in experiment effect and can distinguish ectocinerea from alba.

Description

A kind of method of brain white matter integrity imaging

(1) technical field

The present invention relates to a kind of method of brain white matter integrity imaging, especially based on the partial volume effect in the hybrid response kernel function solution fiber imaging of RL algorithm, utilize diffusion-weighted magnetic resonance imaging (DiffusionWeightedMagneticResonanceImaging, DW-MRI) data and multiple response Kernel-Based Methods carry out the matching of machine direction distribution function in conjunction with RL algorithm, thus obtain the direction of fiber, and obtain machine direction more accurately in grey matter and the obvious intersection of white matter partial volume effect.The machine direction obtained is made to be more conducive to the tracking of fiber.Belong to medical imaging, Nervous System Anatomy field.

(2) background technology

Nuclear magnetic resonance (MRI) be a kind of be widely used in medical imaging without diffusivity method, as unique live body non-invasive method, it obtains in clinical neuromechanism information and the function understood between cerebral cortex region and contact etc. help people and has played huge effect.The trend of brain white matter integrity and mental disorder and cranial surgery medical conditions exist to be maintained close ties with, and these information provide new application prospect for the research of the growth of brain, schizophrenia, congenital and acquired leukoencephalopathy and dementia etc.Fiber imaging algorithm based on diffusion-weighted magnetic resonance imaging (DW-MRI) can obtain machine direction information from DW-MRI data, for clinical diagnose provides foundation, for brain scientific research provides new method.

In all kinds of MRI method, diffusion tensor imaging (DiffusionTensorImaging, DTI) is the one of outbalance, and for multiple brain diseases clinical diagnosis known at present, DTI technology has all played irreplaceable effect.But only containing a fiber in traditional DTI method hypothesis voxel, therefore the fibre structure that such as intersection, bottleneck, dispersion etc. are complicated cannot be differentiated, and the nerve fibre of human brain often exists the complex situations of intersection, branch or fusion, the machine direction that DTI is reconstructed becomes uncertain.

In order to overcome the inherent limitation of DTI, high angular resolution diffusion magnetic resonance imaging (HARDI) technology is arisen at the historic moment.On basis based on HARDI technology, propose the method for multiple fiber reconstruct, such as: Q-ball, diffusion spectrum imaging (DiffusionSpectrumImaging, DSI), sphere deconvolution (Sphericaldeconvolution, SD) etc.At present, although often kind of method all well solves the imaging problem of complicated white matter few fibers, most of HARDI method does not explain that the partial volume effect of non-white matter (grey matter and cerebrospinal fluid) is on the impact of fiber imaging.Not yet there is a kind of mathematical model truly solving non-white confrontation brain fiber Imaging in this field at present.

(3) summary of the invention

In order to overcome weak point above-mentioned in prior art, the present invention proposes a kind of multiple response kernel function that utilizes based on RL algorithm and processes the formation method of the partial volume effect of non-white matter, thus make the machine direction reconstructing high resolution and low error of the fiber of white matter and grey matter intersection, and distinguish white matter of brain and ectocinerea region further.

(1) read brain MR data, obtaining and applying gradient direction is magnetic resonance diffusion signal S (g) of g, does not apply the magnetic resonance diffusion signal S of gradient direction ₀and gradient direction data, pre-service is carried out to the data gathered, comprising: High frequency filter, spatial noise reduction, remove linear drift etc.; Choose required area-of-interest, and calculate diffusive attenuation signal S (the g)/S in this region ₀;

(2) utilize Richardson-Lucy iterative algorithm by diffusive attenuation signal S (the g)/S of each voxel in area-of-interest ₀be modeled as the model with elliposoidal distribution one by one, and increase total variation regularization removal noise effect, modeling process is as follows:

2.1) voxel model micro-structure scheme: by diffusive attenuation signal S (g)/S ₀be assumed to be along rebuilding signal response kernel function R (v, g) and dispersal direction distribution function f (v) of vector v in Spherical Surface S ²on convolution:

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

R (v, g) represents hybrid response kernel function, and it utilizes the isotropy response kernel function in the ultimate fibre response kernel function of white matter of brain and ectocinerea to form, and effectively can solve the partial volume effect in brain white matter integrity imaging, g={g _i∈ R ^{1 × 3}| i=1 ..., n} is Diffusion direction, v={v _j∈ R ^{1 × 3}| j=1 ..., m} is for rebuilding vector; Its mathematical model is:

R(v,g)＝[R _aniR _isot]

Wherein, represent anisotropy response kernel function and isotropy response kernel function in voxel respectively, b is diffusion-sensitive coefficient; Anisotropy response kernel function R _anibe made up of the response core rebuilding direction v along m, each response core is identical round pie, and just their distribution arrangement is different; And isotropy response kernel function is only made up of an independent response core, its shape is spherical shape; $D_{ani}=\left[\begin{array}{ccc}\mathrm{\α}& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right]$ Represent that diffusion is carried out along a principal direction, $D_{iso}=\left[\begin{array}{ccc}\mathrm{\β}& 0& 0\\ 0& \mathrm{\β}& 0\\ 0& 0& \mathrm{\β}\end{array}\right]$ Consistent in its diffusion of all directions, wherein α, β represent fiber diffusion, and its numerical value can be consulted and obtain in pertinent literature;

2.2) a kind of mathematical model based on Richardson-Lucy iterative algorithm is as follows:

Because the data acquisition of magnetic resonance signal is unsatisfactory, often with certain noise, in order to overcome the impact that noise is rebuild machine direction as much as possible, adopt the Richardson-Lucy iterative algorithm in classical image recovery method, the effective nonnegativity ensureing the machine direction distribution function of reconstruct, reduces the quantity at the pseudo-peak of reconstruct; This mathematical model is calculated by the method for maximal possibility estimation to obtain; If diffusion-weighted magnetic resonance signals has n Diffusion direction, and rebuild along m reconstruction vector, then its mathematical model is:

$f{{\left(v\right)}_i}^{(k+1)}=f{{\left(v\right)}_i}^{\left(k\right)}\frac{{\[R{(v,g)}^{T}S\]}_i}{{\[R{(v,g)}^{T}R(v,g)f{\left(v\right)}^{\left(k\right)}\]}_i}$

Wherein k is iterations, and i represents the i-th row of vector, f ^(k)be the length that kth time iteration obtains be the column vector of m, represent the dispersal direction distribution function be evenly distributed on along m reconstruction direction on sphere, S is the length comprising HARDI signal is the column vector of n;

2.3) the TV regularization model based on Richardson-Lucy iterative algorithm is as follows:

Richardson-Lucy algorithm itself has limitation, although this algorithm can the impact on imaging of restraint speckle to a certain extent; But along with the increase of iterations, the optimum solution calculated by Richardson-Lucy algorithm can control by noise; In order to overcome the harmful effect that iterations increase brings, on the basis of this algorithm, add total variation regularization term carry out restraint speckle, its mathematical model is:

$f{{\left(v\right)}_i}^{(k+1)}=f{{\left(v\right)}_i}^{\left(k\right)}\frac{{\[R{(v,g)}^{T}S\]}_i}{{\[R{(v,g)}^{T}R(v,g)f{\left(v\right)}^{\left(k\right)}\]}_i}\×{\mathrm{TV}}_i^{\left(k\right)}$

Wherein,

${\mathrm{TV}}_i^{\left(k\right)}=\frac{1}{1-\mathrm{\λ}div\left(\frac{\▿{\[f^{\left(k\right)}\]}_i}{\[|\▿f^{\left(k\right)}{|}_i\]}\right)}$

TV _i ^(k)it is the total variation regularization term of kth time iteration; the diffuse images direction distribution function f of kth time iteration ^(k)gradient, div represents divergence, and λ is regularization parameter;

(3) obtain dispersal direction distribution function f (v) by iterative computation, the computing method of dispersal direction distribution function f (v) comprise the following steps:

3.1) the individual discrete point of uniform sampling m in unit sphere, is that initial point obtains this m reconstruction vector v with the centre of sphere, calculates the value of fiber response kernel function R (v, g), obtain the circulant matrix of n × (m+1);

3.2) utilize simulated data analog simulation, iterative initial value is set; Make f ⁽⁰⁾for the dispersal direction distribution function of isotropic, its amplitude is set to 1; Because f ⁽⁰⁾initial value is positive, naturally meets so the nonnegativity of algorithm limits; λ value can be selected by experiment, come equilibrium criterion item and regular terms to the impact of iterative algorithm;

3.3) the RL algorithm without regular terms is utilized to carry out pre-service to the voxel of region of interest; Obtain the dispersal direction distribution function f of each voxel, as the initial propagations direction distribution function value of regularization RL algorithm;

3.4) stopping criterion for iteration is set: one is iterations; One is iteration error, makes iteration error be:

$d=\frac{\left|\right|f^{(k+1)}-f^{\left(k\right)}\left|\right|}{\left|\right|f^{\left(k\right)}\left|\right|}$

So iterations is greater than best iterations or iteration error d< ε, and (general ε gets 10 ^-3) as stopping criterion for iteration; (for simulated data, the dispersal direction distribution function solution f of known preferred, so f in above formula iteration error ^(k)=f; )

3.5) by iteration, obtain optimum dispersal direction distribution function f, it is the column vector of m+1 row, wherein the value f of last row _(m+1)be the relative scale size in each voxel shared by isotropic diffusion; And namely front m row form anisotropic fiber direction distribution function fod in voxel, namely machine direction distribution function fod is the front m-1 item of dispersal direction distribution function f; Can know by experiment: work as f _(m+1)during > θ, the diffusion in voxel is isotropy; It is not carried out to the reconstruct of machine direction, and for f _(m+1)the voxel of < θ, utilizes MATLAB to emulate the direction of matching machine direction distribution function fod;

3.6) in perceptive construction on mathematics, three-dimensional imaging is carried out to machine direction distribution function fod, and obtain the principal direction of fiber by the extreme point in search machine direction distribution function value.

The present invention utilizes the RL algorithm of multiple response function to solve non-white matter part to the impact of brain fiber imaging just.The present invention relates to the theory of maximal possibility estimation and Medical Image Processing;

Advantage of the present invention is: computing velocity is fast, and imaging angular resolution is high, and it is high to calculate robustness.It seems according to the author's experience, income effect of the present invention is this field best effects current.

(4) accompanying drawing explanation

Fig. 1 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}{\&Sigma;}}{f}_i{S}_0{e}^{-{\mathrm{bg}}^{T}Dg}$

Wherein f represents that i-th with the ratio shared by fiber, f ₁=0.5, f ₂=0.5, S ₀=1, b=3000s/mm ², the eigenwert of diffusion tensor D is: λ ₁=1.8 × 10 ^-3mm ²/ s, λ ₂=0.3 × 10 ^-3mm ²/ s, λ ₃=0.3 × 10 ^-3mm ²/ s.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, second row represents machine direction, the third line represents imaging model, black line illustrate the direction of two fibers (by calculate diffuse peak obtain).

Fig. 2 is actual clinical effect data figure of the present invention.Real data from Harvard University's hospital attached to a medical college (BrighamandWomen ' sHospital, BrocktonVAHospital, McLeanHospital), utilize the brain data that 3-TGE system and double echo plane imaging sequence extract from true human brain, data acquisition parameters is: TR=17000ms, TE=78ms.Voxel amount is 144 × 144 × 85, become image field be 85 axial slices that 24cm. is parallel to AC-PC line, every layer thickness 1.7mm. from 51 different gradient direction scan-datas, diffusion-sensitive coefficient b=900s/mm2, the scan-data of 8 b=0.

(5) concrete implementation step

The present invention is further illustrated below in conjunction with accompanying drawing.

The method of a kind of brain white matter integrity imaging of the present invention, comprises the steps:

(1) read brain MR data, obtaining and applying gradient direction is magnetic resonance diffusion signal S (g) of g, does not apply the magnetic resonance diffusion signal S of gradient direction ₀and gradient direction data, pre-service is carried out to the data gathered, comprising: High frequency filter, spatial noise reduction, remove linear drift etc.; Choose required area-of-interest, and calculate diffusive attenuation signal S (the g)/S in this region ₀;

(2) utilize Richardson-Lucy iterative algorithm by diffusive attenuation signal S (the g)/S of each voxel in area-of-interest ₀be modeled as the model with elliposoidal distribution one by one, and increase total variation regularization removal noise effect, modeling process is as follows:

2.1) voxel model micro-structure scheme: by diffusive attenuation signal S (g)/S ₀be assumed to be along rebuilding signal response kernel function R (v, g) and dispersal direction distribution function f (v) of vector v in Spherical Surface S ²on convolution:

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

R (v, g) represents hybrid response kernel function, and it utilizes the isotropy response kernel function in the ultimate fibre response kernel function of white matter of brain and ectocinerea to form, and effectively can solve the partial volume effect in brain white matter integrity imaging, g={g _i∈ R ^{1 × 3}| i=1 ..., n} is Diffusion direction, v={v _j∈ R ^{1 × 3}| j=1 ..., m} is for rebuilding vector; Its mathematical model is:

R(v,g)＝[R _aniR _isot]

Wherein, represent anisotropy response kernel function and isotropy response kernel function in voxel respectively, b is diffusion-sensitive coefficient; Anisotropy response kernel function R _anibe made up of the response core rebuilding direction v along m, each response core is identical round pie, and just their distribution arrangement is different; And isotropy response kernel function is only made up of an independent response core, its shape is spherical shape; $D_{ani}=\left[\begin{array}{ccc}\mathrm{\&alpha;}& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right]$ Represent that diffusion is carried out along a principal direction, $D_{iso}=\left[\begin{array}{ccc}\mathrm{\&beta;}& 0& 0\\ 0& \mathrm{\&beta;}& 0\\ 0& 0& \mathrm{\&beta;}\end{array}\right]$ Consistent in its diffusion of all directions, wherein α, β represent fiber diffusion, and its numerical value can be consulted and obtain in pertinent literature;

2.2) a kind of mathematical model based on Richardson-Lucy iterative algorithm is as follows:

Because the data acquisition of magnetic resonance signal is unsatisfactory, often with certain noise, in order to overcome the impact that noise is rebuild machine direction as much as possible, adopt the Richardson-Lucy iterative algorithm in classical image recovery method, the effective nonnegativity ensureing the machine direction distribution function of reconstruct, reduces the quantity at the pseudo-peak of reconstruct; This mathematical model is calculated by the method for maximal possibility estimation to obtain; If diffusion-weighted magnetic resonance signals has n Diffusion direction, and rebuild along m reconstruction vector, then its mathematical model is:

$f{{\left(v\right)}_i}^{(k+1)}=f{{\left(v\right)}_i}^{\left(k\right)}\frac{{\&lsqb;R{(v,g)}^{T}S\&rsqb;}_i}{{\&lsqb;R{(v,g)}^{T}R(v,g)f{\left(v\right)}^{\left(k\right)}\&rsqb;}_i}$

Wherein k is iterations, and i represents the i-th row of vector, f ^(k)be the length that kth time iteration obtains be the column vector of m, represent the dispersal direction distribution function be evenly distributed on along m reconstruction direction on sphere, S is the length comprising HARDI signal is the column vector of n;

2.3) the TV regularization model based on Richardson-Lucy iterative algorithm is as follows:

Richardson-Lucy algorithm itself has limitation, although this algorithm can the impact on imaging of restraint speckle to a certain extent; But along with the increase of iterations, the optimum solution calculated by Richardson-Lucy algorithm can control by noise; In order to overcome the harmful effect that iterations increase brings, on the basis of this algorithm, add total variation regularization term carry out restraint speckle, its mathematical model is:

$f{{\left(v\right)}_i}^{(k+1)}=f{{\left(v\right)}_i}^{\left(k\right)}\frac{{\&lsqb;R{(v,g)}^{T}S\&rsqb;}_i}{{\&lsqb;R{(v,g)}^{T}R(v,g)f{\left(v\right)}^{\left(k\right)}\&rsqb;}_i}\&times;{\mathrm{TV}}_i^{\left(K\right)}$

Wherein,

${\mathrm{TV}}_i^{\left(k\right)}=\frac{1}{1-\mathrm{\&lambda;}div\left(\frac{\&dtri;{\&lsqb;f^{\left(k\right)}\&rsqb;}_i}{\&lsqb;|\&dtri;f^{\left(k\right)}{|}_i\&rsqb;}\right)}$

TV _i ^(k)it is the total variation regularization term of kth time iteration; the diffuse images direction distribution function f of kth time iteration ^(k)gradient, div represents divergence, and λ is regularization parameter;

(3) obtain dispersal direction distribution function f (v) by iterative computation, the computing method of dispersal direction distribution function f (v) comprise the following steps:

3.1) the individual discrete point of uniform sampling m in unit sphere, is that initial point obtains this m reconstruction vector v with the centre of sphere, calculates the value of fiber response kernel function R (v, g), obtain the circulant matrix of n × (m+1);

3.2) utilize simulated data analog simulation, iterative initial value is set; Make f ⁽⁰⁾for the dispersal direction distribution function of isotropic, its amplitude is set to 1; Because f ⁽⁰⁾initial value is positive, naturally meets so the nonnegativity of algorithm limits; λ value can be selected by experiment, come equilibrium criterion item and regular terms to the impact of iterative algorithm;

3.3) the RL algorithm without regular terms is utilized to carry out pre-service to the voxel of region of interest; Obtain the dispersal direction distribution function f of each voxel, as the initial propagations direction distribution function value of regularization RL algorithm;

3.4) stopping criterion for iteration is set: one is iterations; One is iteration error, makes iteration error be:

$d=\frac{\left|\right|f^{(k+1)}-f^{\left(k\right)}\left|\right|}{\left|\right|f^{\left(k\right)}\left|\right|}$

So iterations is greater than best iterations or iteration error d< ε, and (general ε gets 10 ^-3) as stopping criterion for iteration; (for simulated data, the dispersal direction distribution function solution f of known preferred, so f in above formula iteration error ^(k)=f; )

3.5) by iteration, obtain optimum dispersal direction distribution function f, it is the column vector of m+1 row, wherein the value f of last row _(m+1)be the relative scale size in each voxel shared by isotropic diffusion; And namely front m row form anisotropic fiber direction distribution function fod in voxel, namely machine direction distribution function fod is the front m-1 item of dispersal direction distribution function f; Can know by experiment: work as f _(m+1)during > θ, the diffusion in voxel is isotropy; It is not carried out to the reconstruct of machine direction, and for f _(m+1)the voxel of < θ, utilizes MATLAB to emulate the direction of matching machine direction distribution function fod;

3.6) in perceptive construction on mathematics, three-dimensional imaging is carried out to machine direction distribution function fod, and obtain the principal direction of fiber by the extreme point in search machine direction distribution function value.

Claims (3)

1. a method for brain white matter integrity imaging, its formation method comprises the following steps:

(1) read brain MR data, obtaining and applying gradient direction is magnetic resonance diffusion signal S (g) of g, does not apply the magnetic resonance diffusion signal S of gradient direction ₀and gradient direction data, pre-service is carried out to the data gathered, comprising: High frequency filter, spatial noise reduction, remove linear drift etc.; Choose required area-of-interest, and calculate diffusive attenuation signal S (the g)/S in this region ₀;

(2) utilize Richardson-Lucy iterative algorithm by diffusive attenuation signal S (the g)/S of each voxel in area-of-interest ₀be modeled as the model with elliposoidal distribution one by one, and increase total variation regularization removal noise effect, modeling process is as follows:

2.1) voxel model micro-structure scheme: by diffusive attenuation signal S (g)/S ₀be assumed to be along rebuilding signal response kernel function R (v, g) and dispersal direction distribution function f (v) of vector v in Spherical Surface S ²on convolution:

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

R (v, g) represents hybrid response kernel function, and it utilizes the isotropy response kernel function in the ultimate fibre response kernel function of white matter of brain and ectocinerea to form, and effectively can solve the partial volume effect in brain white matter integrity imaging, g={g _i∈ R ^{1 × 3}| i=1 ..., n} is Diffusion direction, v={v _j∈ R ^{1 × 3}| j=1 ..., m} is for rebuilding vector; Its mathematical model is:

R(v,g)＝[R _aniR _isot]

Wherein, represent anisotropy response kernel function and isotropy response kernel function in voxel respectively, b is diffusion-sensitive coefficient; Anisotropy response kernel function R _anibe made up of the response core rebuilding direction v along m, each response core is identical round pie, and just their distribution arrangement is different; And isotropy response kernel function is only made up of an independent response core, its shape is spherical shape; $D_{ani}=\left[\begin{array}{ccc}\mathrm{\&alpha;}& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right]$ Represent that diffusion is carried out along a principal direction, $D_{iso}=\left[\begin{array}{ccc}\mathrm{\&beta;}& 0& 0\\ 0& \mathrm{\&beta;}& 0\\ 0& 0& \mathrm{\&beta;}\end{array}\right]$ Consistent in its diffusion of all directions, wherein α, β represent fiber diffusion, and its numerical value can be consulted and obtain in pertinent literature;

2.2) a kind of mathematical model based on Richardson-Lucy iterative algorithm is as follows:

Because the data acquisition of magnetic resonance signal is unsatisfactory, often with certain noise, in order to overcome the impact that noise is rebuild machine direction as much as possible, adopt the Richardson-Lucy iterative algorithm in classical image recovery method, the effective nonnegativity ensureing the machine direction distribution function of reconstruct, reduces the quantity at the pseudo-peak of reconstruct; This mathematical model is calculated by the method for maximal possibility estimation to obtain; If diffusion-weighted magnetic resonance signals has n Diffusion direction, and rebuild along m reconstruction vector, then its mathematical model is:

$f{{\left(v\right)}_i}^{\left(k+1\right)}=f{{\left(v\right)}_i}^{\left(k\right)}\frac{{\&lsqb;R{(v,g)}^{T}S\&rsqb;}_i}{{\&lsqb;R{(v,g)}^{T}R(v,g)f{\left(v\right)}^{\left(k\right)}\&rsqb;}_i}$

Wherein k is iterations, and i represents the i-th row of vector, f ^(k)be the length that kth time iteration obtains be the column vector of m, represent the dispersal direction distribution function be evenly distributed on along m reconstruction direction on sphere, S is the length comprising HARDI signal is the column vector of n;

2.3) the TV regularization model based on Richardson-Lucy iterative algorithm is as follows:

Richardson-Lucy algorithm itself has limitation, although this algorithm can the impact on imaging of restraint speckle to a certain extent; But along with the increase of iterations, the optimum solution calculated by Richardson-Lucy algorithm can control by noise; In order to overcome the harmful effect that iterations increase brings, on the basis of this algorithm, add total variation regularization term carry out restraint speckle, its mathematical model is:

$f{{\left(v\right)}_i}^{\left(k+1\right)}=f{{\left(v\right)}_i}^{\left(k\right)}\frac{{\&lsqb;R{(v,g)}^{T}S\&rsqb;}_i}{{\&lsqb;R{(v,g)}^{T}R(v,g)f{\left(v\right)}^{\left(k\right)}\&rsqb;}_i}\&times;{{\mathrm{TV}}_i}^{\left(k\right)}$

Wherein,

${{\mathrm{TV}}_i}^{\left(k\right)}=\frac{1}{1-\mathrm{\&lambda;}div\left(\frac{\&dtri;{\&lsqb;f^{\left(k\right)}\&rsqb;}_i}{\&lsqb;|\&dtri;f^{\left(k\right)}{|}_i\&rsqb;}\right)}$

TV _i ^(k)it is the total variation regularization term of kth time iteration; the diffuse images direction distribution function f of kth time iteration ^(k)gradient, div represents divergence, and λ is regularization parameter;

(3) obtain dispersal direction distribution function f (v) by iterative computation, the computing method of dispersal direction distribution function f (v) comprise the following steps:

3.1) the individual discrete point of uniform sampling m in unit sphere, is that initial point obtains this m reconstruction vector v with the centre of sphere, calculates the value of fiber response kernel function R (v, g), obtain the circulant matrix of n × (m+1);

3.2) utilize simulated data analog simulation, iterative initial value is set; Make f ⁽⁰⁾for the dispersal direction distribution function of isotropic, its amplitude is set to 1; Because f ⁽⁰⁾initial value is positive, naturally meets so the nonnegativity of algorithm limits; λ value can be selected by experiment, come equilibrium criterion item and regular terms to the impact of iterative algorithm;

3.3) the RL algorithm without regular terms is utilized to carry out pre-service to the voxel of region of interest; Obtain the dispersal direction distribution function f of each voxel, as the initial propagations direction distribution function value of regularization RL algorithm;

3.4) stopping criterion for iteration is set: one is iterations; One is iteration error, makes iteration error be:

$d=\frac{\left|\right|f^{(k+1)}-f^{\left(k\right)}\left|\right|}{\left|\right|f^{\left(k\right)}\left|\right|}$

So iterations is greater than best iterations or iteration error d< ε, and (general ε gets 10 ^-3) as stopping criterion for iteration; (for simulated data, the dispersal direction distribution function solution f of known preferred, so f in above formula iteration error ^(k)=f; )

3.5) by iteration, obtain optimum dispersal direction distribution function f, it is the column vector of m+1 row, wherein the value f of last row _(m+1)be the relative scale size in each voxel shared by isotropic diffusion; And namely front m row form anisotropic fiber direction distribution function fod in voxel, namely machine direction distribution function fod is the front m-1 item of dispersal direction distribution function f; Can know by experiment: work as f _(m+1)during > θ, the diffusion in voxel is isotropy; It is not carried out to the reconstruct of machine direction, and for f _(m+1)the voxel of < θ, utilizes MATLAB to emulate the direction of matching machine direction distribution function fod;

3.6) in perceptive construction on mathematics, three-dimensional imaging is carried out to machine direction distribution function fod, and obtain the principal direction of fiber by the extreme point in search machine direction distribution function value.

2. the method for claim 1, it is characterized in that: in each region of described 2.1) midbrain, there is isotropy and anisotropy parameter all simultaneously, hybrid response kernel function has fully utilized isotropy and anisotropy response kernel function, isotropic diffusion part in white matter of brain can be removed, make the fiber of imaging have better anisotropy, have better imaging effect; More can distinguish white matter of brain and ectocinerea region further.

3. the method for claim 1, is characterized in that: in described step 2 regularization parameter initial value can not too little can not be too large: if λ is too little, RL algorithmic procedure is mainly controlled by data item model, and regularization Optimization Progress is slower; If λ is too large, RL algorithmic procedure then mainly control by regular terms, acquired results can depart from true solution; The method that can set different parameters value by many experiments selects preferably regularization parameter.