CN106023138A - Global probability fiber tracking method based on spherical convolution - Google Patents

Abstract

The invention provides a global probability fiber tracking method based on spherical convolution. The method comprises the steps that fiber orientation distribution on a unit sphere is calculated through the convolution of a diffusion signal and a response function; the original fiber orientation distribution fODF is replaced by a spherical Gaussian kernel method to acquire a group of new fODF; global probability fiber tracking based on a group tracking algorithm is carried out to acquire prior probability and observation density, and then posterior distribution is calculated; on the basis of the posterior distribution, the sum and average of global fiber functions are calculated; a tracking start point is established to start tracking; and a relatively accurate fiber structure is acquired by choosing the best direction. The global probability fiber tracking method based on spherical convolution effectively considers local fiber orientation distribution and has high precision.

Description

A kind of global probability fiber tracking method based on spherical convolution

Technical field

The present invention relates to the medical imaging under computer graphics, neuroanatomy field, especially a kind of fiber tracking Method.

Background technology

Fibre straighteness is to obtain neuromechanism and dissect unique work of link information based on DWI molding Tool；It estimates possible fiber path by the diffusion directions following the tracks of local tensors orientation；Fiber path is with a complex web Nerve fibre bundle in the form performance human brain of network chart, and by its investigation nuerological pathology and sick crowd；Actual card Bright, in definitiveness tracking technique, fiber path remains variation, it is impossible to pseudo-at noise, division, head movement and image The uncertainty of track is reconstructed under the influence of shadow；In order to overcome these to limit, fibre straighteness method is developed to quantify Uncertainty with visualization with fiber path；Research worker began to use some fine angular resolution diffusion imaging methods in recent years Realize probabilistic algorithm.

Summary of the invention

In order to overcome the deficiency that local fiber distribution of orientations, precision are relatively low that cannot consider of existing fiber tracking mode, this Based on spherical convolution the global probability fiber tracking that invention provides a kind of effective consideration local fiber distribution of orientations, precision is higher Method.

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

A kind of global probability fiber tracking method based on spherical convolution, comprises the steps:

1) diffusion model is estimated

By solving a linear system and reclaiming sphere gaussian kernel problem estimation spherical convolution SD, S (g)/S₀With N number of ladder Degree direction g approximation is by the convolution of receptance functionAnd the fiber orientation distribution in unit sphereRepresent:

$S\left(g\right)/S_0=F\left(\hat{x}\right)\⊗R(g,\hat{x})---\left(1\right)$

Wherein,Being sample convolution direction, definition receptance function is:

$R(g,\hat{x})=e^{-\mathrm{\δ}{(g\·\hat{x})}^2}---\left(2\right)$

δ is the product of diffusion time and diffusion coefficient, therefore by S (g)/S₀WithObtain fiber orientation distributionEstimate by minimizing following energy:

$E=\underset{i=1}{\overset{N}{\Σ}}{\left(S\right(g)/S_0-\underset{j=i}{\overset{M}{\Σ}}{e}^{-\mathrm{\δ}{(g\·\hat{x})}^2}{F}_j\hat{x}(\hat{x}\left)\right)}^2---\left(3\right)$

It is expressed as a linear system by non-negative least square method；It is jth fiber orientation distribution, passes through ball It is each that face gaussian kernel method substitutes the humorous recovery of ballTo regainIn:

$f\left(\hat{x}\right)=F\left(\hat{x}\right)\·G(g,\hat{x})---\left(4\right)$

Wherein,σ is width parameter；

2) whole world probability fiber tracking, process is as follows:

From x on the c of path₀To x_nGlobal cost function at the effect lower aprons of α be:

$E\left(c\right(x_0\→x_n\left)\right)=\mathrm{\α}\underset{i-1}{\overset{n}{\Σ}}p(v_i,x_i)---\left(5\right)$

Cost function p (v_i,x_i) it is x_iPoint selection direction v_iProbability；For each step, direction v_iIt is from probability density letter Number p (v_i|v_i-1, Φ) and sampling obtains, and it is to be represented by a Bayesian frame:

$p\left(v_i\right|v_{i-1},\mathrm{\Φ})=\frac{p\left(\mathrm{\Φ}\right|v_i)p\left(v_i\right|v_{i-1})}{p\left(\mathrm{\Φ}\right)}---\left(6\right)$

Wherein, Φ is the observed value of one group of three-dimensional diffusion weighted imaging volume, and p (Φ) is a fixed standard of this system The factor, distribution p (v_i|v_i-1) be defined as priori conditions and obtained by a simple component cloth；Observing and nursing p (Φ | v_i) generation For in formula (4)

$p\left(\mathrm{\Φ}\right|v_i)=f\left(\hat{x}\right)---\left(7\right)$

After priori and observation density are all calculated, Posterior distrbutionp p (v_i|v_i-1, Φ) drawn by formula (5)；Finally, Use MCMC technology from p (v_i|v_i-1, Φ) and draw the random sample of machine direction；Definition E^sumAnd E^aveAs from a-quadrant to B Two of region measurements connect, and represent summation and the meansigma methods of whole world fiber function respectively:

$\begin{array}{c}{E}^{sum}\left(c\right(x_0\∈A\→x_n\∈B\left)\right)=\mathrm{\α}\underset{i=0}{\overset{n}{\mathrm{\Σ}}}p\left(v_i\right|v_{i-1},\mathrm{\Φ})\\ E^{ave}\left(c\right(x_0\∈A\→x_n\∈B\left)\right)=\frac{\mathrm{\α}}{n+1}\underset{i=0}{\overset{n}{\mathrm{\Σ}}}p\left(v_i\right|v_{i-1},\mathrm{\Φ})\end{array}---\left(8\right)$

E^aveRepresent the quantity of point, for every fiber average probability, E^aveRepresent the meansigma methods in path, but it is likely to answer For too high E^sumValue and mistake indicate some short or dead paths, and E^sumExplicitly indicating that the summation in path, positive is relevant Property may come from some irrational paths；So needing two global measurings to be used for selecting optimal path；

2) group's track algorithm, process is as follows:

Distribute m particle in starting point on one path and elapse propagation over time；At step number t-1 i.e. v_iAnd x (t-1)_i (t-1), under i=1....m and step-length α, order propagates, at next step, the function that granule is previous step, it may be assumed that

$v_i\left(t\right)=\left\{\begin{array}{cc}{v}_i(t-1)+v_i^L\left(t\right),& ifk=1\\ v_{f_i}\left(t\right),& ifk>1\end{array}\right.---\left(9\right)$

x_i(t+1)=x_i(t)+αv_i(t), i=1 ..., m (10)

Refer to local fiber be oriented in step t from formula (6) andSampled representation whole world fibre orientation；It is Adaptive value according to all particles is selected: for each iteration k, k is parameter, and d top fibers descending is stored asWithD is parameter；According to E^sumAnd E^aveThe two is rung The measurement answered, therefore optimal fiber f_k:

$f_k\left(f\right(1),f(2),...f(p\left)\right)=f_k^{sum}\∩f_k^{ave}---\left(11\right)$

WithRepresenting summation and the meansigma methods of optimal fiber respectively, p is coefficient；All of voxel x_1,2,...,mEvery Individual fiber f_kIt is stored in archive Γ:

$\mathrm{\Γ}=\{x_1(v_1^1,v_1^2,...,v_1^n),x_2(v_2^1,v_2^2,...,v_2^n),...,x_m(v_m^1,v_m^2,...,v_m^n)\}---\left(12\right)$

Wherein, m > 0, n > 0,It it is n the path fibre orientation through voxel x；N=0 represents this voxel in formula (12) In be unnecessary whileIt is suitable；Assume that a particle propagation is by currently putting x_iIt is v with the vector before it_i-1, one Good orientation v_iIt is then and v_i-1Lowest difference fromObtaining, formula is i.e.

The technology of the present invention is contemplated that: first with the fiber in the convolutional calculation unit sphere of diffusion signal and receptance function Distribution of orientations, then replace original fODF by sphere gaussian kernel method, obtain one group of new fODF.

Global probability fiber tracking based on group's track algorithm, after obtaining prior probability and observation density, calculates posteriority Distribution；Summation and the meansigma methods of whole world fiber function it is calculated on the basis of learning Posterior distrbutionp.Set up tracking initiation point Start to follow the tracks of, by selecting optimum orientation to obtain relatively accurate fibre structure.

Beneficial effects of the present invention is mainly manifested in: effectively consider that local fiber distribution of orientations, precision are higher；The side of cooperation Formula obtains optimal fibre structure.

Detailed description of the invention

The invention will be further described below.

A kind of global probability fiber tracking method based on spherical convolution, comprises the steps:

1) diffusion model is estimated

By solving a linear system and reclaiming sphere gaussian kernel problem estimation spherical convolution SD, S (g)/S₀With N number of ladder Degree direction g can approximate by the convolution of receptance functionAnd the fiber orientation distribution fODF in unit sphere,Table Show:

$S\left(g\right)/S_0=F\left(\hat{x}\right)\⊗R(g,\hat{x})---\left(1\right)$

Wherein,Being sample convolution direction, definition receptance function is:

$R(g,\hat{x})=e^{-\mathrm{\δ}{(g\·\hat{x})}^2}---\left(2\right)$

δ is the product of diffusion time and diffusion coefficient, therefore by S (g)/S₀WithObtain fiber orientation distribution FODF,Estimate by minimizing following energy:

$E=\underset{i=1}{\overset{N}{\Σ}}{\left(S\right(g)/S_0-\underset{j=i}{\overset{M}{\Σ}}{e}^{-\mathrm{\δ}{(g\·\hat{x})}^2}{F}_j\hat{x}(\hat{x}\left)\right)}^2---\left(3\right)$

It is expressed as a linear system by non-negative least square method；It is jth fiber orientation distribution fODF, Last in order to prevent the sensitivity to noise, substitute the humorous recovery of ball by sphere gaussian kernel method eachTo regaining 'sIn:

$f\left(\hat{x}\right)=F\left(\hat{x}\right)\·G(g,\hat{x})---\left(4\right)$

Wherein,σ is width parameter, and this numerical value structure is for preventing from merging deficiency Or angular resolution reduces.

2) whole world probability fiber tracking, process is as follows:

From x on the c of path₀To x_nGlobal cost function at the effect lower aprons of α be:

$E\left(c\right(x_0\→x_n\left)\right)=\mathrm{\α}\underset{i=1}{\overset{n}{\Σ}}p(v_i,x_i)---\left(5\right)$

Cost function p (v_i,x_i) it is x_iPoint selection direction v_iProbability；For each step, direction v_iIt is from probability density letter Number p (v_i|v_i-1, Φ) and sampling obtains, and it is to be represented by a Bayesian frame:

$p\left(v_i\right|v_{i-1},\mathrm{\Φ})=\frac{p\left(\mathrm{\Φ}\right|v_i)p\left(v_i\right|v_{i-1})}{p\left(\mathrm{\Φ}\right)}---\left(6\right)$

Wherein, Φ is the observed value of one group of three-dimensional diffusion weighted imaging volume, and p (Φ) is a fixed standard of this system The factor, distribution p (v_i|v_i-1) be defined as priori conditions and obtained by a simple component cloth；Observing and nursing p (Φ | v_i) generation For in formula (4)

$p\left(\mathrm{\Φ}\right|v_i)=f\left(\hat{x}\right)---\left(7\right)$

After priori and observation density are all calculated, Posterior distrbutionp p (v_i|v_i-1, Φ) drawn by formula (5)；Finally, Use MCMC technology from p (v_i|v_i-1, Φ) and draw the random sample of machine direction；Definition E^sumAnd E^aveAs from a-quadrant to B Two of region measurements connect, and represent summation and the meansigma methods of whole world fiber function respectively:

$\begin{array}{c}{E}^{sum}\left(c\right(x_0\∈A\→x_n\∈B\left)\right)=\mathrm{\α}\underset{i=0}{\overset{n}{\mathrm{\Σ}}}p\left(v_i\right|v_{i-1},\mathrm{\Φ})\\ E^{ave}\left(c\right(x_0\∈A\→x_n\∈B\left)\right)=\frac{\mathrm{\α}}{n+1}\underset{i=0}{\overset{n}{\mathrm{\Σ}}}p\left(v_i\right|v_{i-1},\mathrm{\Φ})\end{array}---\left(8\right)$

E^aveRepresent the quantity of point, for every fiber average probability, E^aveRepresent the meansigma methods in path, but it is likely to answer For too high E^sumValue and mistake indicate some short or dead paths, and E^sumExplicitly indicating that the summation in path, positive is relevant Property may come from some irrational paths；So needing two global measurings to be used for selecting optimal path；

2) group's track algorithm, process is as follows:

Distribute m particle in starting point on one path and elapse propagation over time；At step number t-1 i.e. v_iAnd x (t-1)_i (t-1), under i=1....m and step-length α, order propagates, at next step, the function that granule is previous step, it may be assumed that

$v_i\left(t\right)=\left\{\begin{array}{cc}{v}_i(t-1)+v_i^L\left(t\right),& ifk=1\\ v_{f_i}\left(t\right),& ifk>1\end{array}\right.---\left(9\right)$

x_i(t+1)=x_i(t)+αv_i(t), i=1 ..., m (10)

Refer to local fiber be oriented in step t from formula (6) andSampled representation whole world fibre orientation；It is Adaptive value according to all particles is selected: for each iteration k, k is parameter, and d top fibers descending is stored asWithD is parameter；According to E^sumAnd E^aveThe two is rung The measurement answered, therefore optimal fiber f_k:

$f_k\left(f\right(1),f(2),...f(p\left)\right)=f_k^{sum}\∩f_k^{ave}---\left(11\right)$

WithRepresenting summation and the meansigma methods of optimal fiber respectively, p is coefficient；All of voxel x_1,2,...,mEvery Individual fiber f_kIt is stored in archive Γ:

$\mathrm{\Γ}=\{x_1(v_1^1,v_1^2,...,v_1^n),x_2(v_2^1,v_2^2,...,v_2^n),...,x_m(v_m^1,v_m^2,...,v_m^n)\}---\left(12\right)$

Wherein, m > 0, n > 0,It it is n the path fibre orientation through voxel x；N=0 represents this voxel in formula (12) In be unnecessary whileIt is suitable；Assume that a particle propagation is by currently putting x_iIt is v with the vector before it_i-1, one Good orientation v_iIt is then and v_i-1Lowest difference fromObtaining, formula is i.e.

Claims (1)

1. a global probability fiber tracking method based on spherical convolution, it is characterised in that: comprise the steps:

1) diffusion model is estimated

By solving a linear system and reclaiming sphere gaussian kernel problem estimation spherical convolution SD, S (g)/S₀With N number of gradient side To g approximation by the convolution of receptance functionAnd the fiber orientation distribution in unit sphereRepresent:

$S\left(g\right)/S_0=F\left(\hat{x}\right)\⊗R(g,\hat{x})---\left(1\right)$

Wherein,Being sample convolution direction, definition receptance function is:

$R(g,\hat{x})=e^{-\mathrm{\δ}{(g\·\hat{x})}^2}---\left(2\right)$

D is the product of diffusion time and diffusion coefficient, therefore by S (g)/S₀WithObtain fiber orientation distributionLogical Cross and minimize following energy and estimate:

It is expressed as a linear system by non-negative least square method；It is jth fiber orientation distribution, high by sphere It is each that this kernel method substitutes the humorous recovery of ballTo the fODF regainedIn:

$f\left(\hat{x}\right)=F\left(\hat{x}\right)\·G(g,\hat{x})---\left(4\right)$

Wherein,σ is width parameter；

2) whole world probability fiber tracking, process is as follows:

From x on the c of path₀To x_nGlobal cost function at the effect lower aprons of α be:

$E\left(c\right(x_0\→x_n\left)\right)=\mathrm{\α}\underset{i-1}{\overset{n}{\Σ}}p(v_i,x_i)---\left(5\right)$

Cost function p (v_i,x_i) it is x_iPoint selection direction v_iProbability；For each step, direction v_iIt is from probability density function p (v_i|v_i-1, Φ) and sampling obtains, and it is to be represented by a Bayesian frame:

$p\left(v_i\right|v_{i-1},\mathrm{\Φ})=\frac{p\left(\mathrm{\Φ}\right|v_i)p\left(v_i\right|v_{i-1})}{p\left(\mathrm{\Φ}\right)}---\left(6\right)$

Wherein, Φ is the observed value of one group of three-dimensional diffusion weighted imaging volume, p (Φ) be this system a fixed standard because of Son, distribution p (v_i|v_i-1) be defined as priori conditions and obtained by a simple component cloth；Observing and nursing p (Φ | v_i) replace In formula (4)

$p\left(\mathrm{\Φ}\right|v_i)=f\left(\hat{x}\right)---\left(7\right)$

After priori and observation density are all calculated, Posterior distrbutionp p (v_i|v_i-1, Φ) drawn by formula (5)；Finally, use MCMC technology is from p (v_i|v_i-1, Φ) and draw the random sample of machine direction；Definition E^sumAnd E^aveAs from a-quadrant to B region Two measurements connect, represent the whole world summation of fiber function and meansigma methods respectively:

$\begin{array}{c}{E}^{sum}\left(c\left(x_0\∈A\→x_n\∈B\right)\right)=\mathrm{\α}\underset{i=0}{\overset{n}{\mathrm{\Σ}}}p\left(v_i|v_{i-1},\mathrm{\Φ}\right)\\ E^{ave}\left(c\left(x_0\∈A\→x_n\∈B\right)\right)=\frac{\mathrm{\α}}{n+1}\underset{i=0}{\overset{n}{\mathrm{\Σ}}}p\left(v_i|v_{i-1},\mathrm{\Φ}\right)\end{array}---\left(8\right)$

E^aveRepresent the quantity of point, for every fiber average probability, E^aveRepresent the meansigma methods in path, but it is likely to should be High E^sumValue and mistake indicate some short or dead paths, and E^sumExplicitly indicating that the summation in path, positive dependency can Can come from some irrational paths；So needing two global measurings to be used for selecting optimal path；

2) group's track algorithm, process is as follows:

Distribute m particle in starting point on one path and elapse propagation over time；At step number t-1 i.e. v_iAnd x (t-1)_i(t- 1), under i=1....m and step-length α, order propagates, at next step, the function that granule is previous step, it may be assumed that

$v_i\left(t\right)=\left\{\begin{array}{cc}{v}_i\left(t-1\right)+v_i^L\left(t\right),& ifk=1\\ v_{f_i}\left(t\right),& ifk>1\end{array}\right.---\left(9\right)$

x_i(t+1)=x_i(t)+αv_i(t), i=1 ..., m (10)

Refer to local fiber be oriented in step t from formula (6) andSampled representation whole world fibre orientation；It is according to institute The adaptive value having particle is selected: for each iteration k, k is parameter, and d top fibers descending is stored asWithD is parameter；According to E^sumAnd E^aveThe two is rung The measurement answered, therefore optimal fiber f_k:

$f_k\left(f\right(1),f(2),...f(p\left)\right)=f_k^{sum}\∩f_k^{ave}---\left(11\right)$

WithRepresenting summation and the meansigma methods of optimal fiber respectively, p is coefficient；All of voxel x_12...mEach fiber f_kIt is stored in archive Γ:

$\mathrm{\Γ}=\left\{x_1\left(v_1^1,v_1^2,...,v_1^n\right),x_2\left(v_2^1,v_2^2,...,v_2^n\right),...,x_m\left(v_m^1,v_m^2,...,v_m^n\right)\right\}---\left(12\right)$

Wherein, m > 0, n > 0,It it is n the path fibre orientation through voxel x；N=0 represents this voxel While unnecessaryIt is suitable；Assume that a particle propagation is by currently putting x_iIt is v with the vector before it_i-1, one optimal Orientation v_iIt is then and v_i-1Lowest difference fromObtaining, formula is i.e.