CN105913458A - Alba fiber imaging method based on colony tracking - Google Patents

Abstract

The present invention provides an alba fiber imaging method based on colony tracking. The method comprises the following steps: 1) the estimation of the fiber individual tracking direction: observing the voxel level direction of diffusion signal, fully employing the structuring direction information existing in the neighborhood voxel, and observing the direction of the current tracking points according to the interpolation of the neighborhood voxel main diffusion direction on the fiber individual historical tracking direction; 2) the estimation of the fiber colony tracking direction: employing a group tracking method based on the random matrix for the fiber colony tracking direction, describing the fiber colony tracking direction by using the state iteration update process based on the Bayes framework, and realizing the non-determinacy filtering caused by measurement of noise and motion model noise, where the model is similar to the simple mode of the Kalman filter; and performing the fiber bundle group object state tracking through the model. The alba fiber imaging method based on colony tracking is able to reduce the sensitivity for the signal noise and the model error and the dependence on the local fiber direction reconstruction module and is high in stability.

Description

A kind of brain white matter integrity formation method followed the tracks of based on colony

Technical field

The present invention relates to the medical imaging under computer graphics, Nervous System Anatomy field, especially a kind of white matter of brain is fine Dimension formation method.

Background technology

Brain is to control the mankind to carry out the synthesizer that the sophisticated functions such as logical thinking, learning and memory, motion and emotion are movable Official, human brain working mechanism is probed into be contemporary scientific research forward position and focus.Brain white matter integrity reconfiguration technique will be by having The brain fibre space Microstructure Information having Anatomical significance carries out imaging, it has also become brain Mechanism of Cognition is probed into, neural class disease The important technical of the sick brain science research such as pathological analysis and brain surgery navigation, thus the development of this technology by information, The common concern of Neuscience related researcher.Diffusion-Weighted MR Imaging (Diffusion Weighted Imaging, DWI) with Based on Magnetic resonance imaging (Magnetic Resonance Imaging, MRI) technology, execute in a plurality of directions by measuring Add the Strength Changes that tissue signal produces before and after diffusion gradient field, it is achieved to the observation of water diffusion motion state in tissue. And the track reconstructing of nerve fibre based on water diffusion anisotropic directional spreding model, passes through mould in portraying voxel Intend the fiber tracking process that in signal transmission passage, hydrone is propagated to describe.Nerve fibre reconstruct based on DWI scanning signal Technology has become unique non-invasive methods of fibr tissue image microstructures in live body brain.

Existing fiber formation method is based on filamentary tracking scheme, and the quality of Fiber track result is completely dependent on In the quality of monomer element local fiber directional spreding reconstruct, owing to the existence of model error makes along with the increasing of fiber path length Length causes error accumulation and uncertain increase, and the least deviation in local fiber direction may greatly change final simultaneously Fiber path is distributed, thus causes the unstability of reconstruction result.Current fiber imaging algorithm is often through improving local side To modeling method reduction model error to realizing the raising of tracking accuracy and stability, but the impact of signal noise is difficult to lead to Cross improvement local fiber direction model to be completely eliminated, thus unavoidably produce model error.

Summary of the invention

In order to the sensitiveness to signal noise and model error overcoming existing brain white matter integrity formation method is higher, play a game The dependence of portion's machine direction reconstruction model is higher, the deficiency of less stable, and the present invention provides a kind of and reduces signal noise And the sensitiveness of model error, reduce to the dependence of local machine direction reconstruction model, stability higher based on colony with The brain white matter integrity formation method of track.

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

A kind of brain white matter integrity formation method followed the tracks of based on colony, comprises the steps:

1) individual fibers tracking direction is estimated

For multifilament local direction distributed model, the machine direction of each voxel is distributed during extreme value vector set closes and history The immediate vector in the direction of propagation is as diffusion principal directionCorresponding weight value ω of each voxel diffusion principal direction is respectively by distance Weights ω₁With synteny weights ω₂Together decide on:

$\mathrm{\ω}(u_{k-1}^i,N_{k-1,m}^i,v_{k-1}^i)={\mathrm{\ω}}_1(u_{k-1}^i,N_{k-1,m}^i){\mathrm{\ω}}_2(u_{k-1}^i,N_{k-1,m}^i,v_{k-1}^i)$

${\mathrm{\ω}}_1(u_{k-1}^i,N_{k-1,m}^i)=\exp (-\frac{1}{2{\mathrm{\σ}}^2}{||u_{k-1}^i-N_{k-1,m}^i||}^2)$

${\mathrm{\ω}}_2(u_{k-1}^i,N_{k-1,m}^i,v_{k-1}^i)={\left(\frac{\left|v_{k-1}^i(u_{k-1}^i-N_{k-1,m}^i)\right|}{||u_{k-1}^i-N_{k-1,m}||}\right)}^z$

Wherein,For eight diagonal angle neighboring voxels coordinate points set of current tracking point, N is the seat of field voxel Mark [x, y, z], k is to follow the tracks of the moment, and m is diagonal angle voxel (8 altogether), and i is i-th individual fibers particle；Distance weights ω₁Profit It is weighted with the Euclidean distance of trace point with neighboring voxels coordinate points, weighs, with its distance, the side that neighboring voxels is provided To information weight；Synteny weights ω₂With the history direction of propagation and neighboring voxels position between the orientation of current tracking point Synteny size, as weight metric, shows that current fibre individuality propagates trend and the neighboring voxels tool on historical trace direction There is higher correlation；

Historical trace directionAs weighted term to portray the inertial force that fiber is propagated, it is based ultimately upon individual fibers The multi-voxel proton neighbor interpolation direction of historical trace stateObtained by following formula:

$d_k^i={\mathrm{\ϵ}}_1(0.3v_{k-1}^i+{\mathrm{\ϵ}}_2\underset{m=1}{\overset{8}{\Σ}}{\mathrm{\ω}}_{k,m}^i{a}_{k,m}^i)$

Wherein ε₁, ε₂For vector normalization coefficient；

2) groups of fibers volume tracing direction estimation

By two multi-agent synergy action directionsRealize individual fibers is followed the tracks of the adjustment of trend: take turns based on fiber beam bundle Wide and the constraint of black and white matter frontier distance group revises directionIn group, individual fibers moves towards similarity constraint direction

2.1) group's profile status tracking model is set up

In the k moment by observing along sub-voxel level directionPropagated forward, obtain comprising n_kIndividual Point Target position quantity The fiber beam bundle observation set surveyedAndRepresent that all fibres multiple targets to the k moment is seen Order-checking row；Multiple targets state joint probability density function p (x under utilizing Bayes principle will to observe at individuality_k,X_k|y_k) with high The product of this density function and inverse Wishart density function represents:

p(x_k,X_k|y_k)=p (Y_k,n_k|x_k,X_k)p(x_k,X_k|y_k-1)

Based on measuring individual number n_kEqual-probability distribution it is assumed that k moment measuring assembly Y_kLikelihood function can be simplified to Equation below:

$p\left(Y_k\right|n_k,x_k,X_k)=\underset{j=1}{\overset{n_k}{\Π}}N(y_k^i;{\mathrm{Hx}}_k,X_k)\∝N({\stackrel{\‾}{y}}_k;{\mathrm{Hx}}_k,\frac{X_k}{n_k})\×W({\stackrel{\‾}{Y}}_k;n_k-1,X_k)$

Wherein, N (g) is Gaussian Profile, and W (g) is Wishart distribution, and H is observing matrix, It is respectively the k moment Measure mean vector and measure stroll matrix；

By prior probability function factorization it is:

$\begin{array}{c}p(x_k,X_k|y_{k-1})=p\left(x_k\right|X_k,y_{k-1})p\left(X_k\right|y_{k-1})\\ =N(x_k;x_{k|k-1})IW(X_k;v_{k|k-1},X_{k|k-1})\end{array}$

Wherein, IW (g) is inverse Wishart distribution, x_k|k-1, X_k|k-1It is respectively scalar parameter and the matrix ginseng of one-step prediction Number.N(x_k；x_k|k-1,P_k|k-1) represent that obeying vector is desired for x_k|k-1, covariance matrix is P_k|k-1Gaussian Profile；

Multiple targets state joint posterior probability density function is multiplied with prior probability function by likelihood function and obtains:

p(x_k,X_k|y_k)∝N(x_k；x_k|k,P_k|k)IW(X_k；v_k|k,X_k|k)

Based on Bayesian model, fiber beam bundle tracking mode include center of mass motion state and contour feature respectively with state with Machine variable x_kAnd X_kPrince's iteration more newly obtained；

2.2) the black and white matter profile constraint to group's direction of motion

Using fiber beam bundle profile front and back the moment distance black and white matter border ratio the most constant as weigh fibre bundle walk To being parallel to the foundation on black and white matter border, utilize frontier distance function to barycenter direction of primary motionThe heart is uniformly distributed and angle Np the direction less than predetermined angleIt is evaluated determining the optimal fibre bundle direction of propagation；

With in three pairwise orthogonal axis of orientations of fibre bundle ellipsoid profile Xk withThe axis of orientation of angle minimum is as tracking Axle, remaining two orthogonal direction axle collectively forms an elliptic plane, tracking axis rotates to candidate direction of primary motion w_tOn Form new elliptic plane；The n that adjacent angle is equal is up-sampled at this ellipsoid_sTo direction vector, and measure on these directions The distance on elliptic contour point and black and white matter borderWherein d_t,2j, d_t,sj+1It is respectively in contrary sample direction Frontier distance；All candidate's direction of primary motions are calculated corresponding frontier distance deviation D (w_t), take the minimum frontier distance of generation Candidate's direction of primary motion of deviation most preferably moves towards as colony, thus obtains group orientation

Wherein, δ is the corresponding sample direction frontier distance that previous moment is most preferably moved towards by colony to determine.

2.3) in group, individual fibers moves towards similarity constraint

Space vector according to k moment difference Fiber track pointDistance sizeUtilize parameter r_min, r_maxDivide three different types of fiber track zones of action: (1) distributing areaFibre in this region Show as sucking action between dimension individuality, promote the bunching effect of similar track；(2) repulsive areaHave less Repelling effect is produced to ensure the diversity of fiber population at individual between the trace point of space length；(3) territory, neutral zone, shows Independence in evolutionary process between the fiber path of larger space distance；

$c_{k,2}^i=-\mathrm{\&tau;}\underset{j\&NotEqual;i}{\overset{n_k}{\&Sigma;}}\left\{\begin{array}{cc}{l}_k^{i,j}/(||l_k^{i,j}||r_{\min }),& ||l_k^{i,j}||<r_{\min }\\ l_k^{i,j}/{||l_k^{i,j}||}^3,& r_{\min }<||l_k^{i,j}||<r_{\max }\\ 0,& ||l_k^{i,j}||>r_{\max }\end{array}\right.$

Wherein, τ is the parameter of control action force intensity；

2.4) group based on k-means separates

Separating judgement by realizing group based on the k-means cluster of individual fibers tracking direction in group, multiple subgroups are divided Realize from the cluster by repeatedly carrying out having Liang Gelei center, cluster centre number is set to 2, by two cluster centres Angle exceed given threshold value and each class comprise individual amount account for group sum 1/3 time as occur groups of fibers separate judgement Condition.

The technology of the present invention is contemplated that: for observing the sub-voxel level direction of diffusion signal, make full use of neighborhood body simultaneously The structuring directional information that element exists, by interpolation to neighboring voxels diffusion directions on individual fibers historical trace direction Obtain the observation of the direction to current tracking point.

Groups of fibers volume tracing direction estimation uses multiple targets tracking based on random matrix, with based on Bayesian frame State iteration renewal process describe, it is achieved the uncertain filtering that measurement noise and motion model noise are brought, this mould Type is similar to the simple form of Kalman filter；By carrying out fibre bundle multiple targets status tracking based on this model.

Beneficial effects of the present invention is mainly manifested in: reduces signal noise and the sensitiveness of model error, reduce and play a game The dependence of portion's machine direction reconstruction model, stability are higher.

Below detailed description of the invention, the invention will be further described.

A kind of brain white matter integrity formation method followed the tracks of based on colony, comprises the steps:

1) individual fibers tracking direction is estimated

The extremal vector that tradition tracking process is estimated using individual fibers current tracking point place voxel directional spreding is as biography Broadcast direction, but monomer cellulose fiber directional spreding model can only obtain the overall estimation of diffusion profile in whole voxel, and not Having the directional resolution of sub-voxel level, the tracking point change in location in the most same voxel does not affect the distribution of tracking direction；

Observe for realizing the sub-voxel level direction to diffusion signal, make full use of the structuring side that neighboring voxels exists simultaneously To information, can by interpolation to neighboring voxels diffusion directions on individual fibers historical trace direction obtain to currently with The direction observation of track point；

For multifilament local direction distributed model, the machine direction of each voxel can be distributed during extreme value vector set closes with The immediate vector in the history direction of propagation is as diffusion principal directionEach voxel diffusion principal direction corresponding weight value ω respectively by Distance weights ω₁With synteny weights ω₂Together decide on:

$\mathrm{\&omega;}(u_{k-1}^i,N_{k-1,m}^i,v_{k-1}^i)={\mathrm{\&omega;}}_1(u_{k-1}^i,N_{k-1,m}^i){\mathrm{\&omega;}}_2(u_{k-1}^i,N_{k-1,m}^i,v_{k-1}^i)$

${\mathrm{\&omega;}}_1(u_{k-1}^i,N_{k-1,m}^i)=\exp (-\frac{1}{2{\mathrm{\&sigma;}}^2}{||u_{k-1}^i-N_{k-1,m}^i||}^2)$

${\mathrm{\&omega;}}_2(u_{k-1}^i,N_{k-1,m}^i,v_{k-1}^i)={\left(\frac{\left|v_{k-1}^i(u_{k-1}^i-N_{k-1,m}^i)\right|}{||u_{k-1}^i-N_{k-1,m}||}\right)}^z$

Wherein,For eight diagonal angle neighboring voxels coordinate points set of current tracking point, N is the seat of field voxel Mark [x, y, z], k is to follow the tracks of the moment, and m is diagonal angle voxel (8 altogether), is i-th individual fibers particle；Distance weights ω₁Profit It is weighted with the Euclidean distance of trace point with neighboring voxels coordinate points, weighs, with its distance, the side that neighboring voxels is provided To information weight；Synteny weights ω₂With the history direction of propagation and neighboring voxels position between the orientation of current tracking point Synteny size, as weight metric, shows that current fibre individuality propagates trend and the neighboring voxels tool on historical trace direction There is higher correlation.

For promoting the flatness of fiber propagated forward, the abnormal conditions preventing wide-angle from shifting occur, historical trace directionAlso serve as weighted term with portray fiber propagate inertial force, be based ultimately upon many bodies of individual fibers historical trace state Element neighbor interpolation directionObtained by following formula:

$d_k^i={\mathrm{\&epsiv;}}_1(0.3v_{k-1}^i+{\mathrm{\&epsiv;}}_2\underset{m=1}{\overset{8}{\&Sigma;}}{\mathrm{\&omega;}}_{k,m}^i{a}_{k,m}^i)$

Wherein ε₁, ε₂For vector normalization coefficient；

2) groups of fibers volume tracing direction estimation

The restriction relation between fibre bundle population mass motion trend is propagated, by two for describing individual fibers direction Multi-agent synergy action directionRealize individual fibers is followed the tracks of the adjustment of trend: based on fiber beam bundle profile and black and white matter border away from Group from constraint revises directionIn group, individual fibers moves towards similarity constraint direction

2.1) group's profile status tracking model is set up

Describe with state iteration renewal process based on Bayesian frame, it is achieved that measurement noise and motion model are made an uproar The uncertain filtering that vocal cords come, and there is the simple form being similar to Kalman filter；Enter based on utilizing this model Row fibre bundle multiple targets status tracking；

In the k moment by observing along sub-voxel level directionPropagated forward, can obtain comprising n_kIndividual Point Target position Put the fiber beam bundle observation set of measurementAndRepresent all fibres group to the k moment Target observation sequence；Based on random matrix multiple targets follow the tracks of basic thought utilize Bayes principle will individuality observe under Multiple targets state joint probability density function p (x_k,X_k|y_k) come with the product of Gaussian density function and inverse Wishart density function Represent:

p(x_k,X_k|y_k)=p (Y_k,n_k|x_k,X_k)p(x_k,X_k|y_k-1)

Based on measuring individual number n_kEqual-probability distribution it is assumed that the Section 1 i.e. k moment measuring assembly Y of above formula_kSeemingly So function can be simplified to equation below:

$p\left(Y_k\right|n_k,x_k,X_k)=\underset{j=1}{\overset{n_k}{\&Pi;}}N(y_k^i;{\mathrm{Hx}}_k,X_k)\&Proportional;N({\stackrel{\&OverBar;}{y}}_k;{\mathrm{Hx}}_k,\frac{X_k}{n_k})\&times;W({\stackrel{\&OverBar;}{Y}}_k;n_k-1,X_k)$

Wherein, N (g) is Gaussian Profile, and W (g) is Wishart distribution, and H is observing matrix, It is respectively the k moment Measure mean vector and measure stroll matrix.

For convenience of the calculating of undated parameter and explanation, construct one suitable conjugate gradient descent method so that Posterior distrbutionp with Prior distribution has identical form, by prior probability function factorization is:

p(x_k,X_k|y_k-1)=p (x_k|X_k,y_k-1)p(X_k|y_k-1)

=N (x_k；x_k|k-1)IW(X_k；v_k|k-1,X_k|k-1)

Wherein, IW (g) is inverse Wishart distribution, x_k|k-1, X_k|k-1It is respectively scalar parameter and the matrix ginseng of one-step prediction Number.N(x_k；x_k|k-1,P_k|k-1) represent that obeying vector is desired for x_k|k-1, covariance matrix is P_k|k-1Gaussian Profile；

Multiple targets state joint posterior probability density function is multiplied with prior probability function by likelihood function and obtains:

p(x_k,X_k|y_k)∝N(x_k；x_k|k,P_k|k)IW(X_k；v_k|k,X_k|k)

Based on above-mentioned Bayesian model, fiber beam bundle tracking mode includes that center of mass motion state and contour feature can be distinguished By state stochastic variable x_kAnd X_kPrince's iteration more newly obtained；

2.2) the black and white matter profile constraint to group's direction of motion,

The research of white matter fiber tracts is shown, and its track trend is often parallel to black and white matter border, thus black and white matter border Trend for fiber provides important information, by the fiber beam bundle profile ratio on front and back's moment distance black and white matter border the most not Become the foundation being parallel to black and white matter border as measurement fibre bundle trend, utilize frontier distance function to barycenter direction of primary motionThe heart is uniformly distributed and the angle n less than 5 °_pIndividual directionIt is evaluated determining the optimal fibre bundle direction of propagation.With fibre Dimension bundle ellipsoid profile Xk three pairwise orthogonal axis of orientations inThe minimum axis of orientation of angle as tracking axis, remaining two Orthogonal direction axle collectively forms an elliptic plane, and tracking axis rotates to candidate direction of primary motion w_tThe new ellipse of upper formation is put down Face.The n that adjacent angle is equal is up-sampled at this ellipsoid_sTo direction vector, and measure elliptic contour point on these directions with The distance on black and white matter borderWherein d_t,2j, d_t,sj+1It is respectively the frontier distance in contrary sample direction.To institute Candidate's direction of primary motion is had to calculate corresponding frontier distance deviation D (w_t), take the main fortune of candidate producing minimum border range deviation Dynamic direction is most preferably moved towards as colony, thus obtains group orientation

$D\left(w_t\right)=\sqrt{\underset{j=1}{\overset{n_s}{\&Sigma;}}{\left(\right(\frac{{\mathrm{\&delta;}}_{2j}}{{\mathrm{\&delta;}}_{2j}+{\mathrm{\&delta;}}_{2j+1}})-(\frac{d_{2j}}{d_{2j}+d_{1j+1}}\left)\right)}^2}$

Wherein, δ is the corresponding sample direction frontier distance that previous moment is most preferably moved towards by colony to determine.

2.3) in group, individual fibers moves towards similarity constraint

Based on there is the fact that Similar Track is distributed between adjacent fiber individuality, the interaction that individual fibers in group is followed the tracks of Firmly repulsion model is attracted to portray, according to the space vector of k moment difference Fiber track point with interparticle's Distance sizeUtilize parameter r_min, r_maxDivide three different types of fiber track zones of action: 1. distributing areaIn this region, between individual fibers, show as sucking action, promote the bunching effect of similar track；2. Repulsive areaHave between the trace point of less space length and produce repelling effect to ensure fiber population at individual Diversity；3. territory, neutral zone, shows the independence in evolutionary process between the fiber path of larger space distance；

$c_{k,2}^i=-\mathrm{\&tau;}\underset{j\&NotEqual;i}{\overset{n_k}{\&Sigma;}}\left\{\begin{array}{cc}{l}_k^{i,j}/(||l_k^{i,j}||r_{\min }),& ||l_k^{i,j}||<r_{\min }\\ l_k^{i,j}/{||l_k^{i,j}||}^3,& r_{\min }<||l_k^{i,j}||<r_{\max }\\ 0,& ||l_k^{i,j}||>r_{\max }\end{array}\right.$

Wherein τ is the parameter of control action force intensity.

2.4) group based on k-means separates

Being affected owing to being broken up by structure of fibrous tissue, have bifurcated, the region complex fiber structures such as sprawling is fine The individual direction of propagation of dimension can separate, and further results in the division of fiber tracking group's profile.Drilling for portraying group direction of motion Many principal direction characteristic that these labyrinth regions occur is entered, it is to avoid the single distributional pattern of group's profile has caused during change The space topological continuity of limit, needs to judge whether fibre bundle colony occurs to separate and to the multiple targets occurring that the direction of propagation breaks up Divide.

Judgement is separated by realizing group based on the k-means cluster of individual fibers tracking direction in group.Due to many height Group is separated and can be realized by the cluster repeatedly carrying out having Liang Gelei center, therefore cluster centre number is set to 2. Two cluster centre angles are exceeded in given threshold value and each class comprise individual amount account for group sum 1/3 time as occur fiber The Rule of judgment that group separates.

Claims (1)

1. the brain white matter integrity formation method followed the tracks of based on colony, it is characterised in that: comprise the steps:

1) individual fibers tracking direction is estimated

For multifilament local direction distributed model, the machine direction of each voxel is distributed during extreme value vector set closes and propagates with history The immediate vector in direction is as diffusion principal directionCorresponding weight value ω of each voxel diffusion principal direction is respectively by distance weights ω₁With synteny weights ω₂Together decide on:

$\mathrm{\&omega;}(u_{k-1}^i,N_{k-1,m}^i,v_{k-1}^i)={\mathrm{\&omega;}}_1(u_{k-1}^i,N_{k-1,m}^i){\mathrm{\&omega;}}_2(u_{k-1}^i,N_{k-1,m}^i,v_{k-1}^i)$

${\mathrm{\&omega;}}_1(u_{k-1}^i,N_{k-1,m}^i)=\exp (-\frac{1}{2{\mathrm{\&sigma;}}^2}||u_{k-1}^i-N_{k-1,m}^i|{|}^2)$

${\mathrm{\&omega;}}_2(u_{k-1}^i,N_{k-1,m}^i,v_{k-1}^i)={\left(\frac{\left|v_{k-1}^i(u_{k-1}^i-N_{k-1,m}^i)\right|}{\left|\right|u_{k-1}^i-N_{k-1,m}\left|\right|}\right)}^z$

Wherein,For eight diagonal angle neighboring voxels coordinate points set of current tracking point, N be field voxel coordinate [x, Y, z], k is to follow the tracks of the moment, and m is diagonal angle voxel, 8 altogether, and i is i-th individual fibers particle；Distance weights ω₁Utilize and follow the tracks of Point is weighted with the Euclidean distance of neighboring voxels coordinate points, weighs, with its distance, the directional information that neighboring voxels is provided Weight；Synteny weights ω₂The synteny between the orientation of current tracking point with the history direction of propagation and neighboring voxels position Size, as weight metric, shows that current fibre individuality is propagated trend and had higher with the neighboring voxels on historical trace direction Correlation；

Historical trace directionAs weighted term with portray fiber propagate inertial force, be based ultimately upon individual fibers history with The multi-voxel proton neighbor interpolation direction of track stateObtained by following formula:

$d_k^i={\mathrm{\&epsiv;}}_1(0.3v_{k-1}^i+{\mathrm{\&epsiv;}}_2\underset{m=1}{\overset{8}{\&Sigma;}}{\mathrm{\&omega;}}_{k,m}^i{a}_{k,m}^i)$

Wherein ε₁, ε₂For vector normalization coefficient；

2) groups of fibers volume tracing direction estimation

By two multi-agent synergy action directionsRealize to individual fibers follow the tracks of trend adjustment: based on fiber beam bundle profile with The group of black and white matter frontier distance constraint revises directionIn group, individual fibers moves towards similarity constraint direction

2.1) group's profile status tracking model is set up

In the k moment by observing along sub-voxel level directionPropagated forward, obtain comprising n_kThe fibre that individual Point Target position measures Dimension beam bundle observation setAndThe all fibres multiple targets observation sequence represented to the k moment Row；Multiple targets state joint probability density function p (x under utilizing Bayes principle will to observe at individuality_k,X_k|y_k) close by Gauss The product of degree function and inverse Wishart density function represents:

p(x_k,X_k|y_k)=p (Y_k,n_k|x_k,X_k)p(x_k,X_k|y_k-1)

Based on measuring individual number n_kEqual-probability distribution it is assumed that k moment measuring assembly Y_kLikelihood function can be simplified to as follows Formula:

$p\left(Y_k\right|n_k,x_k,X_k)=\underset{j=1}{\overset{n_k}{\&Pi;}}N(y_k^i;{\mathrm{Hx}}_k,X_k)\&Proportional;N({\stackrel{\&OverBar;}{y}}_k;{\mathrm{Hx}}_k,\frac{X_k}{n_k})\&times;W({\stackrel{\&OverBar;}{Y}}_k;n_k-1,X_k)$

Wherein, N (g) is Gaussian Profile, and W (g) is Wishart distribution, and H is observing matrix, The measurement being respectively the k moment is equal Value vector sum measures stroll matrix；

By prior probability function factorization it is:

$\begin{array}{c}p(x_k,X_k|y_{k-1})=p\left(x_k\right|X_k,y_{k-1})p\left(X_k\right|y_{k-1})\\ =N(x_k;x_{k|k-1})IW(X_k;v_{k|k-1},X_{k|k-1})\end{array}$

Wherein, IW (g) is inverse Wishart distribution, x_k|k-1, X_k|k-1It is respectively scalar parameter and matrix parameter, the N of one-step prediction (x_k；x_k|k_-1,P_k|k-1) represent that obeying vector is desired for x_k|k-1, covariance matrix is P_k|k-1Gaussian Profile；

Multiple targets state joint posterior probability density function is multiplied with prior probability function by likelihood function and obtains:

p(x_k,X_k|y_k)∝N(x_k；x_k|k,P_k|k)IW(X_k；v_k|k,X_k|k)

Based on Bayesian model, fiber beam bundle tracking mode includes that center of mass motion state and contour feature become at random by state respectively Amount x_kAnd X_kPrince's iteration more newly obtained；

2.2) the black and white matter profile constraint to group's direction of motion

Using the most constant for the fiber beam bundle profile ratio on front and back's moment distance black and white matter border flat as weighing fibre bundle trend Row, in the foundation on black and white matter border, utilizes frontier distance function to barycenter direction of primary motionThe heart is uniformly distributed and angle is less than The n of predetermined angle_pIndividual directionIt is evaluated determining the optimal fibre bundle direction of propagation；

With fibre bundle ellipsoid profile X_kThree pairwise orthogonal axis of orientations inThe minimum axis of orientation of angle as tracking axis, its Two remaining orthogonal direction axles collectively form an elliptic plane, and tracking axis rotates to candidate direction of primary motion w_tUpper formation is new Elliptic plane；The n that adjacent angle is equal is up-sampled at this ellipsoid_sTo direction vector, and measure the ellipse on these directions Profile point and the distance on black and white matter borderWherein d_t,2j, d_t,sj+1The border being respectively in contrary sample direction away from From；All candidate's direction of primary motions are calculated corresponding frontier distance deviation D (w_t), take the time producing minimum border range deviation Select direction of primary motion most preferably to move towards as colony, thus obtain group orientation

Wherein, δ is the corresponding sample direction frontier distance that previous moment is most preferably moved towards by colony to determine；

2.3) in group, individual fibers moves towards similarity constraint

Space vector according to k moment difference Fiber track pointDistance sizeUtilize parameter r_min, r_maxDraw Point three different types of fiber track zones of action: (1) distributing areaIndividual fibers in this region Between show as sucking action, promote the bunching effect of similar track；(2) repulsive areaHave less space away from From trace point between produce repelling effect with ensure fiber population at individual diversity；(3) territory, neutral zone, shows bigger sky Independence in evolutionary process between the fiber path of spacing；

$c_{k,2}^i=-\mathrm{\&tau;}\underset{j\&NotEqual;i}{\overset{n_k}{\&Sigma;}}\left\{\begin{array}{cc}{l}_k^{i,j}/\left(\right|\left|l_k^{i,j}\right|\left|r_{\min }\right),& \left|\right|l_k^{i,j}\left|\right|<r_{\min }\\ l_k^{i,j}/\left|\right|l_k^{i,j}|{|}^3,& r_{\min }<\left|\right|l_k^{i,j}\left|\right|<r_{\max }\\ 0,& \left|\right|l_k^{i,j}\left|\right|>r_{\max }\end{array}\right.$

Wherein, τ is the parameter of control action force intensity；

2.4) group based on k-means separates

Separating judgement by realizing group based on the k-means cluster of individual fibers tracking direction in group, multiple subgroups separate logical Cross and repeatedly carry out having the cluster at Liang Gelei center and realize, cluster centre number is set to 2, by two cluster centre angles Exceed given threshold value and each class comprise individual amount account for group sum 1/3 time as occur groups of fibers separate judgement bar Part.