Summary of the invention
The tracking accuracy that existing white matter of brain 3 D displaying method exists is low, tracking velocity slow, the problem of operating difficulties in order to overcome, and the present invention proposes the brain fibre three-dimensional display packing that a kind of resolution height, tracking velocity are fast, be easy to medical worker's operation.
The technical solution used in the present invention is:
(1) adopt discrete fibre direction density function set up the sphere deconvolution (spherical deconvolution method, SD) model:
This paper proposes a kind of new discrete nonnegativity restrictions in conjunction with the sphere Deconvolution Algorithm Based on Frequency of Gaussian function match: this algorithm at first adopts discrete fibre direction point to set up the convolution model, isolate the machine direction density function to the influence of Deconvolution Algorithm Based on Frequency with this, thereby obtain high angular resolution identification; Match compensates discretization error by the sphere Gaussian function then, thereby eliminates the discretization error of being brought by the discrete pulse function, obtains the machine direction density function more accurately.
(2) carry out colony's fiber reconstruct based on sphere deconvolution model:
Fiber tracking can be regarded a maximum a posteriori probability path problems of seeking between seed point and target area as; In the track algorithm such as Streamline algorithm and Bayesian probability track algorithm of classics, fiber tracking is based on that single fiber carries out, and we propose a kind of parallel colony's fiber tracking algorithm that carries out based on fibrous bundle; In our track algorithm, set suitable number of particles, particle obtains fiber at random according to fiber probability density direction; Fiber is divided into different sections, by use local function model for the tax of each fiber in each section with different weights; The fiber path of selecting weight maximum in each section is local path; Travel through all fibres section, obtain the local path of each section; Local path of following the tracks of successively in each zone obtains the global path of fibrous bundle.
The overall cost function G that we define given fibrous bundle f is:
In the formula, L (f (t), n are the local cost function of fibrous bundle α), and f (t) is the local cost function of certain root fiber f of partial weight maximum, n be in the fibrous bundle with the vector angle of the fiber f fiber number less than parameter alpha; In fact, the path of fibrous bundle can be approximately the path of the fiber of certain root weight maximum in the fibrous bundle, can greatly promote operation efficiency like this;
Article one, the white matter of brain fiber path can be regarded the path of some discrete points in the image space, i.e. f={x as
0, x
1..., x
n, suppose that the step-length of institute's directed quantity is identical, i.e. χ
i=χ, i=1 ..., n; A recursion can be write as in path in discrete time:
x
i(t+1)=x
i(t)+χv
i(t),i=1,...,n-1 (2)
Wherein, v
i(t) be position x
iOn machine direction, suppose that A is the initiation region, and B is the target area, then the total value S of a fiber path from A to B can be approximately:
We then can draw the trend of fibrous bundle fast with the near-optimization machine direction of all fibres direction in the fibrous bundle.
(3) three-dimensional visualization of data shows, its major technique is characterized as:
We are by calling the reconstruct data that VTK medical image show tools comes visual fiber.VTK has powerful said three-dimensional body drafting, iso-surface patch function.Adopt modular OO technology, basic function is packed with class, and system architecture is simple, and the user only need call the class of corresponding function, just can not realize and do not need to write algorithm from bottom, thereby make the user can be easy to study and grasp VTK in very short time.Simultaneously, show the storehouse than Direct X, its curved surface display capabilities is more powerful, is more suitable for the isostructural demonstration of fiber.
The implementation procedure that concrete fiber shows is as follows:
1. it is as follows to ask for sphere deconvolution model process:
1.1 ask for the convolution model.With observation signal be expressed as wall scroll fiber signals response function R and machine direction density letter P (v) convolution on the whole:
Wherein g is unit diffusion pulsed gradient vector, and S (g) is the measuring-signal at pulsed gradient direction g place, S
0Measuring-signal during for the no pulse gradient, v are the unit direction vector, and P (v) is the weighted array of a series of basis functions on the sphere.
The S discretize is about to discrete n the direction v of unit of being of unit sphere
j(j=1 ..., n), the density function value of this direction is P (v
j); The signal response function adopts as drag:
V wherein
jBe j sample direction vector on the sphere, l is the parameter that sign fiber anisotropy degree and b value (diffusion-sensitive coefficient) influence jointly to signal attenuation.The convolutional that is got discretize by (4) (5) formula is:
1.2 find the solution the machine direction density function based on formula (6).Given m observation signal S (g
i), i=1 ..., m, the definition energy function:
Wherein, g
iRepresent i diffusion gradient pulse vector.Because increase the anti-noise ability that non-negativity constraint will significantly improve deconvolution, this method is found the solution (7) formula with non-negative least square method.In order to reduce the higher continuous machine direction density function of precision, choose the sphere Gaussian function:
P wherein
jBe the average of Gaussian function, σ is the variance of Gaussian function, and v is three-dimensional vector of unit length, definition least square surface fitting function:
Then problem is converted to following system of equations and finds the solution problem:
Aω=λ (10)
Wherein ω is for comprising ω
jN dimension unknown number column vector, λ is for comprising λ
iThe m dimensional vector, A is the matrix of m * n;
Find the solution (11) formula with non-negative least square method and obtain ω, we can obtain the direction density letter (FOD) of fiber by formula (12).
The process in the posterior probability path of the overall fibrous bundle of 2 acquisitions is as follows:
2.1 set m particle: x in the first step
i(t), i=1 ..., m;
Calculate diffusion tensor D 2.2 utilize least square method
i(t);
2.3 the fiber that obtains is divided into l interval, in each interval, studies the local cost function of each bar fiber:
Suppose that we study i interval, are defined as follows fiber weight calculation model:
W
i=W(f
j,θ,s) (13)
Add up for certain root fiber f
j, be in the cylinder space of s at its field radius, vector direction and f
jThe deviation angle of fiber vector direction is less than the fibre number n of θ
i, and with it as fiber f
jWeight factor W
j
2.4 the fiber of repeating step 2.3 in interval i all is traversed, and selects the d root fiber f that weight coefficient surpasses setup parameter λ
i, i=1 ..., d also records corresponding coefficient n
i, i=1 ..., n in d is interval as i with its local cost function f (t)
iThe local cost function L of the fibrous fibrous bundle of root (f (t))
i
2.5 repeating step 2.3,2.4 is until l interval local cost function L (f (t))
iAll obtained, we stipulate this n
iIndividual particle advances side by side along this track;
2.6 calculate the appropriate value of fiber path.β the particle that selection has optimum value saves as
And according to formula (14), calculate β particle respectively at the fiber path of l section.
x
i(t+1)=x
i(t)+χv
i(t),i=1,…,β (14)
3 with reconstruct data data three-dimensional demonstrationization process:
3.1 show obtaining of data.The fiber path S that is produced by colony fiber restructing algorithm is by formula (15), comes match by a series of discrete data point, and at first we need resolve the fibrous bundle data.Set up an example Data who in the VTK storehouse, is used for the vtkData class of store data for this reason and preserve the fiber data that is generated by colony's fiber restructing algorithm, concrete operation method is: deposit in the Data class after the coordinate of each discrete data point of fiber is added a skew,, and it is stored in the internal memory; We can revise the fiber that will show by the data of revising among the Data;
Show experimental data 3.2 set up three dimensions.Data three-dimensional among the Data be shown, we need set up a Virtual Space for it, comprising background data (VtkWindow) and agent (VtkActor) two parts; The VtkWindow class example that we set up a Window by name generates a spatial context, can regulate background color, whether have that coordinate shows, content such as rotatable whether by setup parameter; Simultaneously, we set up the VtkActor class of an Actor by name, and the data among the Data are passed to Actor, and by the parameter of configuration Actor, we can set and show whether thing is stretched, shows the color of thing and the side-play amount of position; The demonstration control that calls VTK then loads the Window example, for Window sets up a Virtual Space, simultaneously the data message among the Actor is shown to the coordinate of each data point correspondence; By this kind method, we can show the brain fibre three-dimensional quickly and accurately.
The more traditional Streamline algorithm based on the Q-ball imaging model of the colony's fiber tracking algorithm based on sphere convolution model that the present invention proposes has significant advantage.
We adopt the Hough data model that contains two bundle decussating fiberss to verify the superiority of this method.Its concrete parameter is as follows: voxel amount 50 * 50 * 1, voxel size 1 * 1 * 1mm, coefficient of diffusion b=1500s/mm
2, 76 diffusion gradient directions.Carry out digital simulation under the same conditions with SD and two kinds of methods of Q-ball respectively, as shown in Figure 1.
We have chosen an intersection region in the model as the research sample, (a) are the regional model figure of SD formation method, (b) are the regional model figure of Q-ball method.As seen, the SD model keeps single characteristic direction at non-infall, approximate with Q-ball model estimated result, and move towards clear at this model of fiber infall, more can reflect interfibrous distributed intelligence exactly, realistic fiber situation than the Q-ball model of pie.
We set up 50 seed points at reference position m place, normal distribution μ=0 is satisfied in generation, r=0, μ=0, r=5%, μ=0, the zero-mean Gaussian noise of r=10% (wherein, r represents noise level) each 50 groups, and it is put on three groups of DW-MRI data, use the Streamline algorithm based on DTI respectively, Streamline algorithm based on Q-ball, based on SD colony track algorithm 450 groups of DW-MRI data are implemented fiber tracking, the record track path is counted out with the seed that arrives n place, target location, tracking results as shown in Figure 2, and the seed that record arrives the target area counts out, it is as shown in table 1 to calculate its ratio.
Table 1 fiber tracking is comparison sheet as a result
We almost can not arrive the target area as scheduled based on the Streamline algorithm of DTI model, so this model can not solve complex fiber structures from experimental data as can be seen; The Streamline tracking effect of Q-ball model has improvement slightly, yet far away from the colony's track algorithm under the SD model.During noiseless, our method fiber is densely distributed, seed point almost can both arrive the target area, when noise was 5%, the fiber of seed zone still can arrive the target area as scheduled, and the Q-ball model then has part fiber sideslip, when noise jamming increases to 10% left and right sides, undesirable by the tracking fiber effects that the Q-ball model obtains, and SD model gained fiber is subjected to the influence of noise little, still keeps higher accuracy.
This shows that the colony's fiber tracking algorithm keeps track based on the SD model that we propose is accurate, simultaneously by good anti-interference effect.
Embodiment
With reference to accompanying drawing 1-3:
(1) asks for the convolution model.With observation signal be expressed as wall scroll fiber signals response function R and machine direction density letter P (v) convolution on the whole:
Wherein g is unit diffusion pulsed gradient vector, and S (g) is the measuring-signal at pulsed gradient direction g place, S
0Measuring-signal during for the no pulse gradient, v are the unit direction vector, and P (v) is the weighted array of a series of basis functions on the sphere.
The S discretize is about to discrete n the direction v of unit of being of unit sphere
j(j=1 ..., n), the density function value P (v of this direction
j); The signal response function adopts as drag:
V wherein
jBe j sample direction vector on the sphere, l is the parameter that sign fiber anisotropy degree and b value (diffusion-sensitive coefficient) influence jointly to signal attenuation.The convolutional that is got discretize by (1) (2) formula is:
(2) find the solution the machine direction density function based on formula (3).Given 450 observation signal S (g
i), i=1 ..., 450, the definition energy function:
Wherein, g
iRepresent i diffusion gradient pulse vector.Because increase the anti-noise ability that non-negativity constraint will significantly improve deconvolution, this method is found the solution (4) formula with non-negative least square method.In order to reduce the higher continuous machine direction density function of precision, choose the sphere Gaussian function:
P wherein
jBe the average of Gaussian function, σ is the variance of Gaussian function, and v is three-dimensional vector of unit length, definition least square surface fitting function:
Then problem is converted to following system of equations and finds the solution problem:
Aω=λ (7)
Wherein ω is for comprising ω
jN dimension unknown number column vector, λ is for comprising λ
iThe m dimensional vector, A is the matrix of m * n;
Find the solution (8) formula with non-negative least square method and obtain ω, we can obtain the direction density letter (FOD) of fiber by formula (9).
Its effect as shown in Figure 1, we have set up the SD model of Hough data.
(3) set 50 particle: x
i(t), i=1 ..., 50;
(4) utilize least square method to calculate diffusion tensor D
i(t);
(5) fiber that obtains is divided into l interval, in each interval, studies the local cost function of each bar fiber:
Suppose that we study i interval, are defined as follows fiber weight calculation model:
W
i=W(f
j,θ,s) (10)
Add up for certain root fiber f
j, be in the cylinder space of s at its field radius, vector direction and f
jThe deviation angle of fiber vector direction is less than the fibre number n of θ
i, and with it as fiber f
jWeight factor W
j
(6) fiber of repeating step (5) in interval i all is traversed, and selects the d root fiber f that weight coefficient surpasses setup parameter λ
i, i=1 ..., d also records corresponding coefficient n
i, i=1 ..., n in d is interval as i with its local cost function f (t)
iThe local cost function L of the fibrous fibrous bundle of root (f (t))
i
(7) repeating step (5), (6) are until l interval local cost function L (f (t))
iAll obtained, we stipulate this n
iIndividual particle advances side by side along this track;
(8) appropriate value of calculating fiber path.β the particle that selection has optimum value saves as
And according to formula (11), calculate β particle respectively at the fiber path of l section.
x
i(t+1)=x
i(t)+χv
i(t),i=1,…,β (11)
Effect is as accompanying drawing 2(c) shown in, 50 seed point parallel arranged of setting are followed the tracks of, and good antijamming capability is arranged.
(9) show obtaining of data.The fiber path S that is produced by colony fiber restructing algorithm is by formula (11), comes match by a series of discrete data point, and at first we need resolve the fibrous bundle data.Set up an example Data who in the VTK storehouse, is used for the vtkData class of store data for this reason and preserve the fiber data that is generated by colony's fiber restructing algorithm, concrete operation method is: deposit in the Data class after the coordinate of each discrete data point of fiber is added a skew,, and it is stored in the internal memory; We can revise the fiber that will show by the data of revising among the Data;
(10) set up three dimensions and show experimental data.Data three-dimensional among the Data be shown, we need set up a Virtual Space for it, comprising background data (VtkWindow) and agent (VtkActor) two parts; The VtkWindow class example that we set up a Window by name generates a spatial context, can regulate background color, whether have that coordinate shows, content such as rotatable whether by setup parameter; Simultaneously, we set up the VtkActor class of an Actor by name, and the data among the Data are passed to Actor, and by the parameter of configuration Actor, we can set and show whether thing is stretched, shows the color of thing and the side-play amount of position; The demonstration control that calls VTK then loads the Window example, for Window sets up a Virtual Space, simultaneously the data message among the Actor is shown to the coordinate of each data point correspondence; By this kind method, we can come out the brain fibre three-dimensional quickly and accurately, and effect as shown in Figure 3.