CN101872385B - Fast-marching fiber tracking method based on topology preservation - Google Patents

Abstract

The invention belongs to the field of resonance diffusion imaging, in particular to a fast-marching fiber tracking method based on topology preservation, comprising the steps of: reading DTI (Diffusion Tensor Imaging) data; manually selecting a seed point and initializing; evoluting from the seed point to peripheral neighbor points by using a fast-marching method and adopting a curvature weighting speed function, recording topology information in an evolution process, i.e. a source node of each evoluting point; calculating all paths by adopting a time gradient descent method; and selecting a real fiber path by a link matrix. The invention has the advantages of more according with the real fiber path by fiber tracking and reducing urious positive branches; and the obtained fiber is smooth and better reflects the fiber direction.

Description

Fast-marching fiber tracking method based on the topology maintenance

Technical field

The invention belongs to magnetic resonance diffusion tensor imaging field, relate to the fast-marching fiber tracking method that a kind of topology keeps.

Background technology

Diffusion tensor imaging (DTI) is that method is sent out in a kind of new non-invasive imaging that grows up on magnetic resonance imaging (MRI) basis, it utilizes the micromechanism of the anisotropic principle detection of the free thermal motion of hydrone tissue in the tissue, reaches the purpose of researching human body function.Diffusion tensor imaging is unique non-invasive formation method that can show the white matter of brain fibrous bundle at live body at present, it makes the institutional framework and the connectedness thereof of assessment white matter of brain fibrous bundle become possibility, have general MRI and check incomparable superiority, and opened up new bright prospects for neuroimaging.

In recent years, the impressive progress of DTI research is be applied to the white matter of brain fiber visual, wherein fiber tracking or be called the imaging of white matter of brain fiber be DTI visual in a focus of research.In the white matter of brain fiber, the hydrone that spreads along its long axis direction spreads owing to resistance is less comparatively fast, and spreads slower owing to resistance is big perpendicular to the hydrone of its long axis direction diffusion.Fiber tracking uses the trend and the distribution of continuous curve representation fiber, can obtain three-dimensional continuous white matter of brain institutional framework by fiber tracking method, can show white matter of brain fiber details.Present reconstruction to the white matter of brain nerve fibre is the research emphasis in the magnetic resonance diffusion imaging field, the white matter of brain fiber orientation is studied and helped to obtain the information that brain structure is connected with function with contact, and can be for because the medical diagnosis on disease that fiber lacks or textural anomaly causes provides effective information.This not only helps the structure of deep understanding human brain fiber, and very big value is arranged clinically.

In DTI, diffusion tensor is the symmetric matrix of a 3*3, and 6 components are wherein arranged.This diagonalization of matrix can be obtained 3 eigenwerts and pairing proper vector thereof.Wherein the direction of the pairing proper vector of eigenvalue of maximum (principal vector) is exactly the main direction of water diffusion, also is considered to the direction of the interior fiber of this tensor place voxel (will represent certain thickness three-dimensional volume element to be called voxel) usually.Present fiber tracking method has a variety of, and they probably can be divided into two kinds: based on the tensor territory with based on the global energy Method for minimization.Fiber tracking algorithm based on the tensor territory mainly is to utilize local tensors information to carry out fiber tracking, and DTI can produce the preferred diffusion direction of each voxel, and the arrangement of each some tensor is called the tensor territory on the space.Carry out fiber tracking at first and be calculating this working direction by certain the some beginning on the fiber, follow the tracks of a segment distance along this vectorial direction after, with new some point to start with on the track, these points are coupled together again, just can show tracked fiber.Fiber tracking algorithm key based on the tensor territory is determining of current some dispersal direction, but there is a shortcoming in it, can not handle fiber bifurcated and intersection, and algorithm often goes to stop at crotch region and intersection.And can handle bifurcated and decussating fibers based on the minimized method of global energy, and this method has increased the consideration of fiber uncertainty and randomness, by minimizing cost function, seeks optimal path, and the path with minimum cost is corresponding with true path.Fast-marching fiber track algorithm based on level set (a kind of numerical technique that is used for interface tracking and shape modeling) is a kind of based on the minimized method of global energy, it comes control curved surface to develop by defining a velocity function, seeks optimal path by detecting all possible paths of ordering with neighbours on every side in the iterative process that at every turn develops.Prop up but introduced too many false sun in the evolutionary process, and describe not obvious for fiber orientation.

Fiber tracking provides a kind of method to study white matter of brain institutional framework and connectivity, but fiber tracking still has certain limitation.Evaluation to live body fiber tracking result still lacks goldstandard.DTI is the unique method that live body shows the nerve fibre beam trajectory, because tissue specimen dissect, in the processing procedures such as freezing, dehydration, fixing, section and dissolving, its micromechanism must change, and then the generation geometry deformation, the application organizes method has great difficulty in external checking live body tracking results.Also do not have at present a kind of algorithm can obtain everybody generally approval, therefore propose a kind of reliably, effectively, the fiber tracking algorithm is a research emphasis in current research field fast.

Summary of the invention

(problem such as Fast Marching, algorithm FM) exist more false sun to prop up, and fiber orientation is not obvious the present invention proposes the fast-marching fiber tracking method that a kind of topology keeps at advancing fast based on global energy is minimized.For this reason, the present invention adopts following technical scheme:

The fast-marching fiber tracking method that a kind of topology keeps is characterized in that, comprises the following steps:

The first step: DTI data processing;

Second step: seed points is chosen and initialization

(4) specify initial moving point, i.e. seed points r ', the reference position that develops as curved surface;

(5) determine velocity function

Wherein,

R ' is a moving point, and n (r ') expression curved surface is to the evolution direction of r ', and n (r) expression curved surface is to the evolution direction of r, e ₁(r) the principal vector direction of expression r, e ₁The principal vector direction of (r ') expression r ', evolution direction n (r) use n (r)=| r-r ' | approximate, fiber curvature n (r ')=| r '-r " |;

(6) determine arrowband point and point of distance, and determine T time of arrival of arrowband point according to following formula, be initially T=TIME_MAXE the time of arrival of point of distance, and set up the minimum ordering heap of arrowband point:

Wherein T (r ') is the time of arrival of r ';

The 3rd step: adopt the velocity function of curvature weighting to carry out the curved surface evolution that topology keeps;

(4) in minimum ordering heap, export the some r ' that arrives the time T minimum, it is labeled as moving point, and deletion from minimum ordering heap;

(5) the abutment points r of expedition point r ' is for each abutment points r, if moving point is then constant; If point of distance then is revised as the arrowband point, calculate its time of arrival of T (r), write down its father node r ', and put it in the minimum ordering heap; If the arrowband point then recomputates its time of arrival of T (r),, adjust minimum ordering heap if T time of arrival (r), then upgrades the time of arrival and the father node r ' of this point less than the time of arrival of last iteration;

(6) surpass volume data volume scope if curved surface develops, then stop iteration; Set a time threshold, when surpass this time threshold the time of arrival of point, stop iteration; Perhaps pre-defined iterations reaches iterations and stops, then loop ends, otherwise forward (1) to, continue to carry out;

The 3rd step: adopt the time gradient descent method to determine all fibres path;

The 4th step: utilize UNICOM's matrix to select the true fiber path.

Fast-marching fiber tracking method of the present invention, in evolutionary process, add topology and keep model, and curvature is incorporated into velocity function, flexional is considered in the global energy scope, thereby better control evolutionary process, avoided original algorithm of advancing fast to use local similar as unique problem of estimating.Method provided by the invention can be good at reacting the bifurcated information of fiber, meets the model of true fiber, has kept good topological structure, in terms of existing technologies, it is accurate to have the fiber tracking result, and false sun props up minimizing, to advantages such as noise robustness are good.

Description of drawings

Fig. 1 is that flow process of the present invention is always schemed;

The process flow diagram of the algorithm of advancing fast that Fig. 2 the present invention adopts;

Fig. 3 is the algorithm synoptic diagram of advancing fast that the present invention adopts;

The evolutionary model synoptic diagram that the topology that Fig. 4 the present invention adopts keeps;

Fig. 5 adopts method of the present invention to handle the computer screen output result that the projection fibre corona radiata obtains;

The partial enlarged drawing of Fig. 5 of Fig. 6 computer screen output.

Embodiment

In view of the problem that fast-marching fiber tracking method exists, the fast-marching fiber tracking method that provides a kind of topology to keep comprises the following aspects referring to Fig. 1 the present invention:

The A.DTI data processing

B. seed points is chosen and initialization

C. adopt the velocity function of curvature weighting to carry out the curved surface evolution that topology keeps

D. adopt the time gradient descent method to determine all paths

E. utilize UNICOM's matrix to select the true fiber path

Referring to Fig. 2, introduce the fast-marching fiber tracking method of the present invention's topology maintenance below in detail and how to utilize it to realize fiber tracking.

The DTI data processing

Read picture and be organized into three-dimensional data, owing to the influence of electronic interferences and external environment, medical image often contains noise in imaging process, therefore uses wave filter to carry out denoising.

Seed points is chosen and initialization

The user specifies seed points, as the reference position (being the reference position of fiber tracking) that curved surface develops, carries out initialization then:

1) moving point: moving point be exactly in the grid time of arrival known point, the seed points of user's appointment just in the time of initial, the evolution curved surface is through the point of time T=0 of this point.

2) arrowband point: all movable neighborhoods of a point, its time of arrival, time T=1/F and put into a minimum ordering heap the inside with the arrowband point, and wherein F is the velocity function of curvature weighting, be the important function of control curved surface evolution direction and speed, by formula in the present invention calculate (2).

3) point of distance: be the point of infinite distance time of arrival, T=MAX_TIME, and except moving point and arrowband point all is point of distance.

Adopt the velocity function of curvature weighting to carry out the curved surface evolution that topology keeps

1) in minimum ordering heap, export the some r ' that arrives the time T minimum, it is labeled as moving point, and deletion from minimum ordering heap.

2) i.e. this neighbours' point spatially of the abutment points of investigating r ', as shown in Figure 3, for each abutment points r, if moving point is then constant; If point of distance then is revised as the arrowband point, and calculate its time of arrival, write down its father node r ', and put it in the minimum ordering heap according to formula (1).If the arrowband point then recomputates its time of arrival according to formula (1), if time of arrival T less than time of arrival of last iteration, then upgrade the time of arrival and the father node r ' of this point, adjust minimum ordering heap.

In iterative process each time, curved surface will develop to r from 1 r ', r is the neighbor node of r ', and in the arrowband, has minimum time of arrival, in like manner, r ' is by r " develop, think in evolutionary model so r ", and → r ' → r is the topology information in the evolutionary process, and r " → r ' → r is on fiber path, as shown in Figure 4.

$T\left(r\right)=T\left(r^{\′}\right)+\frac{|r-r^{\′}|}{F\left(r\right)}...\left(1\right)$

Wherein T (r ') is the time of arrival of r ', and F (r) velocity function has guaranteed time T (r) from seed points, along n (r) direction, advances with speed F (r).

$F\left(r\right)=\frac{C\left(r\right)}{1-\min (\left(\right|e_1\left(r\right)\·n\left(r\right)\left|\right),\left(\right|e_1\left(r^{\′}\right)\·n\left(r\right)\left|\right),\left(\right|e_1\left(r\right)\·e_1\left(r^{\′}\right)\left|\right))}...\left(2\right)$

$C\left(r\right)=\left\{\begin{array}{cc}n\left(r^{\&prime;}\right)\&CenterDot;n\left(r\right)& n\left(r^{\&prime;}\right)\&CenterDot;n\left(r\right)>0\\ 0& n\left(r^{\&prime;}\right)\&CenterDot;n\left(r\right)<0\end{array}\right....\left(3\right)$

Wherein r ' is a moving point, and n (r ') expression curved surface is to the evolution direction of r ', and n (r) expression curved surface is to the evolution direction of r, e ₁(r) the principal vector direction of expression r, e ₁The principal vector direction of (r ') expression r '.Evolution direction n (r) use n (r)=| r-r ' | approximate, the topological n (r ') of reservation=| r '-r " | be fiber curvature.When C (r)＜0, F (r)=0 to guarantee the curved surface expansion always of developing, can not shrink.

The curvature weighting velocity function is based on the consideration of following four aspects:

A. when same moving point arrives different neighborhoods (different arrowbands point) speed and equates, select the direction of curvature minimum to advance.

B. when different moving points arrived same vicinity (identical arrowband point) and equate time of arrival, selecting the less moving point of curvature was father node.

C. when fiber orientation becomes not obvious, when the vector field direction is inconsistent, the orientation preferentially that curvature deflection is less, keep fiber continuously and slickness.

D. guarantee the direction that curved surface develops and extends along fiber, avoid occurring X and T word shape fiber path.

The velocity function of curvature weighting is taken flexional into account, and when the evolution curved surface reaches the fiber border, because the restriction of curvature, irregular vector field will be left out, and quantitative change is big because it causes bending energy.In fibrous bundle inside, the velocity function of curvature weighting is by minimizing flexional simultaneously, and the fiber tracking process that is can be advanced along the direction of tensor field, has reduced false sun and has propped up.

3) surpass volume data volume scope if curved surface develops, then stop iteration; Define a time threshold, when surpass this threshold value the time of arrival of point, think that this point and seed points correlation degree are lower, and do not have the white matter of brain fiber between the seed points, stop iteration; Perhaps pre-defined iterations reaches iterations and stops, then loop ends, otherwise forward 1 to), continue to carry out.

Adopt the time gradient descent method to determine all paths

Suppose to exist a fiber γ, and suppose that length is L, τ is a bit on the γ, and time T is under the effect of speed F

Accumulated result.The algorithm of advancing fast guarantees the evolutionary process of the minimal path generation from seed points A to r, and T (r) is exactly a minimum cost.If there is a paths in A to r then satisfies:

$T\left(r\right)=\underset{\mathrm{\&lambda;}}{\min }{\&Integral;}_A^r|\&dtri;T\left(r\left(\mathrm{\&tau;}\right)\right)|\mathrm{d\&tau;}...\left(4\right)$

After iteration finished, by the time gradient descent method, calculating the shortest path that turns back to seed points was exactly a fiber path from the wavefront position so.

Utilize UNICOM's matrix to select the true fiber path

$\mathrm{\&phi;}\left(\mathrm{\&gamma;}\right)=\min \underset{\mathrm{\&tau;}}{F\left(\mathrm{\&gamma;}\left(\mathrm{\&tau;}\right)\right)}=\underset{\mathrm{\&tau;}}{\min }\frac{1}{|\&dtri;T\left(\mathrm{\&gamma;}\left(\mathrm{\&tau;}\right)\right)|}...\left(5\right)$

It is the point of T from any time of arrival afterwards that the algorithm iteration of advancing fast finishes, can link to each other with seed points by the minimal path algorithm, the judgement of true path will be passed through the matrix φ of UNICOM, set a φ threshold value, by the φ Function Mapping, the most probable path that has than strong connectedness with seed points is selected, and those paths that do not satisfy condition are considered to the impossible fiber that exists, and is then deleted.

The method of advancing fast that the topology that adopts the present invention to propose keeps is rebuild for the white matter of brain fiber, and reconstruction algorithm is from seed points.Data reconstruction adopts the EP/SE sequence from GE MEDICAL SYSTEMS.Seed points is selected the intersection of corpus callosum and sagittal figure, and Fig. 5 is for selecting the projection fibre corona radiata, and the result that the method for advancing fast of using topology of the present invention to keep obtains well distinguishes different fiber orientation among the figure.Fiber path presents tree structure as can be seen from Figure 6, has well reflected the bifurcated information of fiber, meets the model of true fiber, and the result has kept good topological structure.

Claims (1)

1. the fast-marching fiber tracking method that topology keeps is characterized in that, comprises the following steps:

The first step: DTI data processing;

Second step: seed points is chosen and initialization

(1) specifies initial moving point, i.e. seed points r ', the reference position that develops as curved surface;

(2) determine velocity function $F\left(r\right)=\frac{C\left(r\right)}{1-\min (\left(\right|e_1\left(r\right)\&CenterDot;n\left(r\right)\left|\right),\left(\right|e_1\left(r^{\&prime;}\right)\&CenterDot;n\left(r\right)\left|\right),({|e}_1\left(r\right)\&CenterDot;e_1\left(r^{\&prime;}\right)|))},$ Wherein, $C\left(r\right)=\left\{\begin{array}{cc}n\left(r^{\&prime;}\right)\&CenterDot;n\left(r\right)& n\left(r^{\&prime;}\right)\&CenterDot;n\left(r\right)>0\\ 0& n\left(r^{\&prime;}\right)\&CenterDot;n\left(r\right)<0\end{array}\right.,$ R ' is a moving point, and n (r ') expression curved surface is to the evolution direction of r ', and n (r) expression curved surface is to the evolution direction of r, e ₁(r) the principal vector direction of expression r, e ₁The principal vector direction of (r ') expression r ', evolution direction n (r) use n (r)=| r-r ' | approximate, fiber curvature n (r ')=| r '-r " |;

(3) determine arrowband point and point of distance, and determine T time of arrival of arrowband point according to following formula, be initially T=TIME_MAXE the time of arrival of point of distance, and set up the minimum ordering heap of arrowband point:

$T\left(r\right)=T\left(r^{\&prime;}\right)+\frac{|r-r^{\&prime;}|}{F\left(r\right)},$ Wherein T (r ') is the time of arrival of r ';

The 3rd step: adopt the velocity function of curvature weighting to carry out the curved surface evolution that topology keeps

(1) in minimum ordering heap, export the some r ' that arrives the time T minimum, it is labeled as moving point, and deletion from minimum ordering heap;

(2) the abutment points r of expedition point r ' is for each abutment points r, if moving point is then constant; If point of distance then is revised as the arrowband point, calculate its time of arrival of T (r), write down its father node r ', and put it in the minimum ordering heap; If the arrowband point then recomputates its time of arrival of T (r),, adjust minimum ordering heap if T time of arrival (r), then upgrades the time of arrival and the father node r ' of this point less than the time of arrival of last iteration;

(3) surpass volume data volume scope if curved surface develops, then stop iteration; Set a time threshold, when surpass this time threshold the time of arrival of point, stop iteration; Perhaps pre-defined iterations reaches iterations and stops, then loop ends, otherwise forward (1) to, continue to carry out;

The 4th step: adopt the time gradient descent method to determine all fibres path

Suppose to exist a fiber γ, and suppose that length is L, τ is a bit on the γ, and time T is under the effect of speed F

Accumulated result, if exist a paths to satisfy from seed points A to r:

After iteration finished, by the time gradient descent method, calculating the shortest path that turns back to seed points promptly was a fiber path from the wavefront position so;

The 5th step: utilize UNICOM's matrix to select the true fiber path

It is the point of T from any time of arrival afterwards that iteration finishes, all can be by the minimal path algorithm to link to each other with seed points, the judgement of true path is undertaken by the matrix φ of UNICOM, set a φ threshold value, by the φ Function Mapping, the most probable path that has than strong connectedness with seed points is selected, and those paths that do not satisfy condition are considered to the impossible fiber that exists, and is then deleted.

Images