CN103279633A - Brain fiber three-dimensional display method based on diffusion-weighted magnetic resonance data - Google Patents

Abstract

The invention discloses a brain fiber three-dimensional display method based on diffusion-weighted magnetic resonance data and belongs to the field of medical image processing in the field of information technology, and brain fiber tracts can be reconstructed fast and accurately with the method. A new spherical surface deconvolution algorithm fit by combination of discrete nonnegative constrains and the Gaussian function is provided to establish a spherical surface deconvolution model of the diffusion-weighted magnetic resonance data, meanwhile, a group fiber tracking algorithm with the fiber tracts as basic tracking units is provided, and the group fiber tracking algorithm can be applied to track the spherical surface deconvolution model to reconstruct brain fibers fast. Data generated by the group fiber tracking algorithm are visualized through a visualization tool kit so that the brain fibers can be displayed accurately. The method simplifies the calculation, promotes tracking speed, and meanwhile has good anti-interference performance.

Description

Brain fibre three-dimensional display packing based on diffusion-weighted MR data

Technical field

The present invention relates to the field of medical image processing under the areas of information technology.

Background technology

Utilizing diffusion-weighted magnetic resonance imaging (DW-MRI) data reconstruction white matter of brain nerve fibre is the important means that present live body shows the white matter of brain fiber orientation, it has opened up brand-new approach for the space microstructure of research white matter of brain, for the research of growth, schizophrenia, congenital and acquired leukoencephalopathy and the dementia etc. of brain provides new method.Its basic thought is the diffusion tensor of at first estimating each voxel of diffusion tensor magnetic resonance image (MRI), obtains the machine direction distribution function with this, utilizes the fiber tracking algorithm to form the fibre space microstructure with anatomy meaning then.Because the intuitive of this method, and be easy to realize the fibre structure comparative analysis of different objects, thereby caused whole world researcher's concern.Use at present wider be based on the diffusion tensor model (Diffusion Tensor Imaging, Streamline fiber tracking display packing DTI), but also have problem at present:

At first, have the diffusion tensor method of estimation now and be difficult to portray complicated brain fiber microstructure.The diversity of the complicacy of human brain structure and fiber connection mode is the main cause that faces difficulty at present.Tensor commonly used at present estimates it is to utilize the diffusion tensor imaging model to set up voxel fiber diffusion tensor, and utilizes least square method to calculate the optimum second-order tensor matrix of acquisition.But an essence limitation of DTI imaging model is that it supposes that each voxel only contains a nerve fibre.In fact, the resolution of DW-MRI data is usually at 3～15mm ³Between, the fiber bundle diameters in the fiber tracking algorithm is the millimeter level, and actual nerve fibre diameter is micron order.Therefore, using the DTI model is difficult to portray such as complicated fibre structure such as intersection, bifurcated, bottle-neck type.

Secondly, use wider fiber restructing algorithm such as Streamline algorithm, Bayesian probability track algorithm at present, its tracing process is based on all that single fiber carries out, and in fact the brain fiber number of one group of common DTI data just reaches the 1000-10000 root, therefore will follow the tracks of out all fibers calculates length consuming time, and counting yield is low.Recently, paid close attention to greatly based on the fibrous bundle track algorithm of technology for global optimization, become the new trend of quick reconfiguration nerve fibre.

And aspect the white matter of brain fibre three-dimensional is visual, the method that generally adopts showed control for using Direct X at present, and is with the modeling of brain fiber reconstruct data, by this model being imported virtual three-dimensional space that it is visual.Programming has higher requirements but Direct X is, needs the professional to set up fiber model, is not suitable for common medical worker and uses.

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:

$G\left(f(x_0\→x_n)\right)={\∫}_{x_0}^{x_n}L(f\left(t\right),n,\mathrm{\α})\mathrm{dt}---\left(1\right)$

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 ₀, x ₁..., 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 ⁱ(t+1)=x ⁱ(t)+χv ⁱ(t),i=1,...,n-1 （2）

Wherein, v ⁱ(t) be position x ⁱOn 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:

$S\left(f(x_0\→x_n)\right)=\mathrm{\χ}\underset{i-1}{\overset{n}{\mathrm{\Σ}}}f\left(x_i\right)---\left(3\right)$

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:

$S=S\left(g\right)/{\mathrm{<figure-callout\; id="0"\; label="S"\; filenames="DEST\_PATH\_HDA00003477761300012.png"\; state="\{\{state\}\}">S</figure-callout>}}_0=R\&CircleTimes;P\left(v\right)---\left(4\right)$

Wherein g is unit diffusion pulsed gradient vector, and S (g) is the measuring-signal at pulsed gradient direction g place, S ₀Measuring-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:

$R(g,v_j)=e^{-l{\left(g^T{v}_j\right)}^2}---\left(5\right)$

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:

$S\left(g\right)/{\mathrm{<figure-callout\; id="0"\; label="S"\; filenames="DEST\_PATH\_HDA00003477761300012.png"\; state="\{\{state\}\}">S</figure-callout>}}_0=\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}{e}^{-l{\left(g^T{v}_j\right)}^2}P\left(v_j\right)---\left(6\right)$

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:

$E_n=\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{(S\left(g_i\right)/S_0-\underset{j=i}{\overset{n}{\mathrm{\&Sigma;}}}{e}^{-l{\left(g_i^T{v}_j\right)}^2}P\left(v_j\right))}^2---\left(7\right)$

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:

$f(p,\mathrm{\&sigma;})={\mathrm{\&omega;}}_j{e}^{(-{\cos }^{-1}\left(\right|(p,v)\left|\right)/{\mathrm{\&sigma;}}^2)}---\left(8\right)$

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:

$E=\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{({\mathrm{\&lambda;}}_i-\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&omega;}}_j{e}^{({-\cos }^{-1}\left(\right|(p_j,v_i)\left|\right)/{\mathrm{\&sigma;}}^2)})}^2---\left(9\right)$

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;

$A_{i,j}=e^{({-\cos }^{-1}\left(\right|(p_j,v_i)\left|\right)/{\mathrm{\&sigma;}}^2)}---\left(11\right)$

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).

$\mathrm{FOD}\left(v\right)=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&omega;}}_j{e}^{({-\cos }^{-1}\left(\right|(p_j,v)\left|\right)/{\mathrm{\&sigma;}}^2)}---\left(12\right)$

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 ⁱ(t), i=1 ..., m;

Calculate diffusion tensor D 2.2 utilize least square method ⁱ(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 ⁱ(t+1)=x ⁱ(t)+χv ⁱ(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;

$S\left(f(x_0\&RightArrow;x_n)\right)=\mathrm{\&chi;}\underset{i-1}{\overset{n}{\mathrm{\&Sigma;}}}f\left(x_i\right)---\left(15\right)$

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 ², 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.

Description of drawings

Accompanying drawing 1 is the SD model of Hough data of the present invention, and (a) (b) is respectively SD and Q-ball regional model figure

(a) in the accompanying drawing 2 (b) is respectively Streamline fiber tracking effect based on DTI and Q-ball model

(c) be colony's fiber tracking effect based on the SD model

Accompanying drawing 3 is based on colony's fiber tracking method for visualizing design sketch of SD model

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:

$S=S\left(g\right)/{\mathrm{<figure-callout\; id="0"\; label="S"\; filenames="DEST\_PATH\_HDA00003477761300012.png"\; state="\{\{state\}\}">S</figure-callout>}}_0=R\&CircleTimes;P\left(v\right)---\left(1\right)$

Wherein g is unit diffusion pulsed gradient vector, and S (g) is the measuring-signal at pulsed gradient direction g place, S ₀Measuring-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:

$R(g,v_j)=e^{-l{\left(g^T{v}_j\right)}^2}---\left(2\right)$

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:

$S\left(g\right)/{\mathrm{<figure-callout\; id="0"\; label="S"\; filenames="DEST\_PATH\_HDA00003477761300012.png"\; state="\{\{state\}\}">S</figure-callout>}}_0=\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}{e}^{-l{\left(g^T{v}_j\right)}^2}P\left(v_j\right)---\left(3\right)$

(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:

$E_n=\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{(S\left(g_i\right)/S_0-\underset{j=i}{\overset{n}{\mathrm{\&Sigma;}}}{e}^{-l{\left(g_i^T{v}_j\right)}^2}P\left(v_j\right))}^2---\left(4\right)$

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:

$f(p,\mathrm{\&sigma;})={\mathrm{\&omega;}}_j{e}^{(-{\cos }^{-1}\left(\right|(p,v)\left|\right)/{\mathrm{\&sigma;}}^2)}---\left(5\right)$

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:

$E=\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{({\mathrm{\&lambda;}}_i-\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&omega;}}_j{e}^{({-\cos }^{-1}\left(\right|(p_j,v_i)\left|\right)/{\mathrm{\&sigma;}}^2)})}^2---\left(6\right)$

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;

$A_{i,j}=e^{({-\cos }^{-1}\left(\right|(p_j,v_i)\left|\right)/{\mathrm{\&sigma;}}^2)}---\left(8\right)$

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).

$\mathrm{FOD}\left(v\right)=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&omega;}}_j{e}^{({-\cos }^{-1}\left(\right|(p_j,v)\left|\right)/{\mathrm{\&sigma;}}^2)}---\left(9\right)$

Its effect as shown in Figure 1, we have set up the SD model of Hough data.

(3) set 50 particle: x ⁱ(t), i=1 ..., 50;

(4) utilize least square method to calculate diffusion tensor D ⁱ(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 ⁱ(t+1)=x ⁱ(t)+χv ⁱ(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;

$S\left(f(x_0\&RightArrow;x_n)\right)=\mathrm{\&chi;}\underset{i-1}{\overset{n}{\mathrm{\&Sigma;}}}f\left(x_i\right)---\left(12\right)$

(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.

Claims (4)

1. based on the brain fibre three-dimensional display packing of diffusion-weighted magnetic resonance (DW-MRI) data, mainly comprise following three steps:

Step 1, adopt discrete fibre direction density function to set up sphere deconvolution model:

A kind of new discrete nonnegativity restrictions is 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;

Step 2, carry out the 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 proposes a kind of parallel colony's fiber tracking algorithm that carries out based on fibrous bundle; In 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;

Step 3, the white matter of brain fibre three-dimensional is visual:

We utilize VTK (Visualization Toolkit) medical tool case and three-dimensional reconstruction technology thereof to make fiber orientation visual; This method for visualizing is at first set up the class of a store data, and the fiber data that generates in the step 2 is mapped to such; Set up an example that shows class then, regulate the position, color, drawing effect of demonstration etc. by the parameter that disposes this demonstration class; The example that at last this is comprised fiber data and display parameter imports the demonstration control of VTK, and then the fibre three-dimensional that reconstruct can be come out shows.

2. brain fibre three-dimensional display packing as claimed in claim 1, it is characterized in that: the sphere deconvolution model based on discrete fibre direction density function in the step 1 may further comprise the steps:

1.1 set up the convolution model:

With observation signal be expressed as wall scroll fiber signals response function R and machine direction density function 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 ₀Measuring-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 (1) (2) formula is:

1.2 find the solution the machine direction density function:

Given m observation signal S (g _i), i=1 ..., m, definition energy function

Wherein, g _iRepresent i diffusion gradient pulse vector.Result of study shows increases the anti-noise ability that non-negativity constraint will significantly improve deconvolution, and 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 battle array of m * n;

Find the solution (8) formula with non-negative least square method and obtain ω, we can obtain the direction density function of fiber by formula (9).

。

3. brain fibre three-dimensional display packing as claimed in claim 2 is characterized in that: the colony's fiber reconstructing method in the step 2, and the process in posterior probability path that obtains overall fibrous bundle is as follows:

2.1 set m particle: x in the first step ⁱ(t), i=1 ..., m;

Calculate diffusion tensor D 2.2 utilize least square method ⁱ(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) （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

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 the appropriate value of fiber path, β the particle of selecting to have optimum value saves as And according to formula (11), suppose that the step-length of institute's directed quantity is identical, i.e. χ _i=χ, i=1 ..., n calculates β particle respectively at the fiber path of l section;

x ⁱ(t+1)=x ⁱ(t)+χv ⁱ(t),i=1,…,β （11）

Wherein, v ⁱ(t) be position x ⁱOn 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:

2.7 match discretize section obtains overall fibrous bundle trend; 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, and then we can be similar to formula (14) and draw fiber path.

G(f(x ₀→x _n))≈S(f(x ₀→x _n)) （14） 。

4. as claim 3 described white matter of brain fibre three-dimensional method for visualizing, it is characterized in that, step 3 specifically:

3.1 show obtaining of data:

The fiber path S that is produced by colony fiber restructing algorithm is by formula (12), 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 is used for the vtkData class of store data in the VTK storehouse for this reason and preserve the fiber data that is generated by colony's fiber restructing algorithm, concrete operation method is: the coordinate of each discrete data point of fiber is added a skew (x ₀, y ₀, z ₀) after deposit in the Data class x wherein in ₀, y ₀, z ₀Be respectively the side-play amount of data in three dimensions, 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.

Images