CN101872385A - Fast Marching Fiber Tracking Method Based on Topology Preservation - Google Patents

Abstract

The invention belongs to the field of resonance diffusion imaging, and relates to a topology-preserving fast-traveling fiber tracking method, which includes reading DTI data; manually selecting seed points and initializing them; using the fast-traveling method to adopt a curvature-weighted velocity function from the seed point to the surrounding neighbors Point evolution, record topological information during the evolution process, that is, the source node of each evolution point; use the time gradient descent method to calculate all paths; select the real fiber path through the connectivity matrix. The invention has the advantages that the fiber tracking is more in line with the real fiber path, and the false male branch is reduced, and the obtained fiber is smooth, which well reflects the fiber direction.

Description

Fast Marching Fiber Tracking Method Based on Topology Preservation

technical field

The invention belongs to the field of magnetic resonance diffusion tensor imaging, and relates to a topology-preserving fast-moving fiber tracking method.

Background technique

Diffusion tensor imaging (DTI) is a new non-invasive imaging method developed on the basis of magnetic resonance imaging (MRI). It uses the anisotropy of free thermal motion of water molecules in tissues to detect the microstructure of tissues. To achieve the purpose of studying human body function. At present, diffusion tensor imaging is the only non-invasive imaging method that can display white matter fiber bundles in vivo. It makes it possible to evaluate the tissue structure and connectivity of white matter fiber bundles, and has incomparable advantages over general MRI examinations. Neuroimaging opens up new and vast vistas.

In recent years, an important progress in DTI research is to apply it to the visualization of white matter fibers, and fiber tracking or white matter fiber imaging is a hot spot in DTI visualization. In white matter fibers, water molecules diffusing along its long axis diffuse faster due to less resistance, while water molecules diffusing perpendicular to its long axis diffuse slowly due to greater resistance. Fiber tracking uses continuous curves to represent the direction and distribution of fibers. Through the fiber tracking method, a three-dimensional continuous white matter tissue structure can be obtained, and the details of white matter fibers can be displayed. At present, the reconstruction of white matter nerve fibers is the focus of research in the field of magnetic resonance diffusion imaging. The study of the direction and connection of white matter fibers is helpful to obtain information on brain structure and functional connections, and can provide insight into the pathological changes caused by fiber loss or structural abnormalities. provide effective information for disease diagnosis. This not only contributes to a deep understanding of the structure of human brain fibers, but also has great clinical value.

In DTI, the diffusion tensor is a 3*3 symmetric matrix with 6 components. Diagonalizing this matrix yields three eigenvalues and their corresponding eigenvectors. The direction of the eigenvector (principal vector) corresponding to the largest eigenvalue is the main direction of water molecule diffusion, and is usually considered to be the voxel where the tensor is located (the volume unit representing a three-dimensional space of a certain thickness is called a voxel) The orientation of the inner fibers. At present, there are many fiber tracking methods, which can be roughly divided into two types: based on tensor fields and based on global energy minimization methods. The fiber tracking algorithm based on the tensor field mainly uses local tensor information for fiber tracking. DTI can generate the preferred diffusion direction of each voxel, and the arrangement of each point tensor in space is called the tensor field. The initial fiber tracking starts from a certain point on a fiber, calculates the forward direction of this point, and after a certain distance is tracked along the direction of the vector, then a new point on the trajectory is used as the starting point, and these points are connected to display Fibers being tracked. The key of fiber tracking algorithm based on tensor field is to determine the direction of current point diffusion, but it has a disadvantage that it cannot handle fiber bifurcation and intersection, and the algorithm often stops at bifurcation area and intersection. The method based on global energy minimization can deal with bifurcation and crossing fibers. This method increases the consideration of fiber uncertainty and randomness. By minimizing the cost function, the optimal path is sought, and the path with the minimum cost is consistent with the real path. correspond. The Fast Marching Fiber Tracking algorithm based on level set (a numerical technique for interface tracking and shape modeling) is a method based on global energy minimization, which controls the surface evolution by defining a velocity function, each evolution iteration In the process, the optimal path is sought by detecting all possible paths with surrounding neighbor points. However, too many false male branches were introduced during the evolution process, and the description of the fiber direction was not obvious.

Fiber tracking provides a method to study white matter organization and connectivity, but fiber tracking still has certain limitations. There is no gold standard for evaluating the results of live fiber tracking. DTI is the only way to display the trajectory of nerve fiber bundles in vivo, because the microstructure of tissue samples will inevitably change during the process of anatomy, freezing, dehydration, fixation, slicing and dissolution, resulting in geometric deformation. The application of histological methods in It is very difficult to verify the in vivo tracking results in vitro. At present, there is no algorithm that can be generally recognized by everyone, so proposing a reliable, effective and fast fiber tracking algorithm is a research focus in the current research field.

Contents of the invention

Aiming at the problems that the Fast Marching (FM) algorithm based on global energy minimization has many false positive branches and the fiber direction is not obvious, the present invention proposes a topology-preserving Fast Marching fiber tracking method. For this reason, the present invention adopts following technical scheme:

A topology-preserving fast-moving fiber tracking method, characterized in that it comprises the following steps:

The first step: DTI data processing;

Step 2: Seed point selection and initialization

(4) Designate the initial active point, that is, the seed point r', as the starting position of the surface evolution;

其中，

r′为活动点，n(r′)表示曲面到r′的演化方向，n(r)表示曲面到r的演化方向，e₁(r)表示r的主向量方向，e₁(r′)表示r′的主向量方向，演化方向n(r)使用n(r)＝|r-r′|近似，纤维曲率n(r′)＝|r′-r″|；(5) Determine the speed function

in,

r′ is the active point, n(r′) represents the evolution direction from the surface to r′, n(r) represents the evolution direction from the surface to r, e ₁ (r) represents the principal vector direction of r, e ₁ (r′) Indicates the principal vector direction of r′, the evolution direction n(r) is approximated by n(r)=|rr′|, and the fiber curvature n(r′)=|r′-r″|;

(6) Determine narrowband point and away from point, and determine the time of arrival T of narrowband point according to the following formula, the time of arrival away from point is initially T=TIME_MAXE, and set up the minimum sorting heap of narrowband point:

其中T(r′)为r′的到达时间；

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

The third step: use the curvature-weighted velocity function to perform topology-preserving surface evolution;

(4) Output the point r' with the smallest arrival time T from the minimum sorting heap, mark it as an active point, and delete it from the minimum sorting heap;

(5) Investigate the adjacent point r of the active point r′. For each adjacent point r, if it is an active point, it will remain unchanged; Parent node r', and put it into the minimum sorting heap; if it is a narrowband point, recalculate its arrival time T(r), if the arrival time T(r) is less than the arrival time of the last iteration, update the point Arrival time and parent node r′, adjust the minimum sort heap;

(6) If the surface evolution exceeds the volume range of the volume data, stop the iteration; set a time threshold, when the arrival time of the point exceeds the time threshold, stop the iteration; or pre-define the number of iterations, and stop when the number of iterations is reached, then the loop ends , otherwise go to (1) and continue to execute;

The third step: use the time gradient descent method to determine all fiber paths;

Step 4: Use the connectivity matrix to select real fiber paths.

In the fast-moving fiber tracking method of the present invention, a topology-preserving model is added to the evolution process, and the curvature is introduced into the velocity function, and the bending energy is considered within the global energy range, thereby better controlling the evolution process and avoiding the original rapid The marching algorithm uses local similarity as the only metric for the problem. The method provided by the invention can well reflect the bifurcation information of the fiber, conform to the model of the real fiber, and maintain a good topology structure. Compared with the prior art, it has the advantages of accurate fiber tracking results, reduced false male branches, and is beneficial to It has the advantages of good noise robustness.

Description of drawings

Fig. 1 is a general flow chart of the present invention;

The flowchart of the fast marching algorithm that Fig. 2 the present invention adopts;

Fig. 3 is a schematic diagram of the fast marching algorithm adopted in the present invention;

Fig. 4 is a schematic diagram of the evolution model of topology preservation adopted by the present invention;

Fig. 5 adopts the method of the present invention to process the output result of the computer screen obtained by projecting the fiber crown;

Figure 6 is a partially enlarged view of Figure 5 output from the computer screen.

Detailed ways

Referring to Fig. 1, the present invention provides a topology-preserving fast-moving fiber tracking method in view of the problems existing in the fast-traveling fiber tracking method, including the following aspects:

A. DTI data processing

B. Seed point selection and initialization

C. Topology-preserving surface evolution using a curvature-weighted velocity function

D. Determine all paths using the time gradient descent method

E. Use the connectivity matrix to select real fiber paths

Referring to FIG. 2 , the topology-preserving fast-travel fiber tracking method of the present invention and how to use it to realize fiber tracking will be introduced in detail below.

·DTI data processing

Read pictures and organize them into three-dimensional data. During the imaging process, due to electronic interference and the influence of the external environment, medical images often contain noise, so filters are used for denoising.

·Seed point selection and initialization

The user specifies the seed point as the starting position of surface evolution (that is, the starting position of fiber tracking), and then initializes:

1) Active point: The active point is the point in the grid whose arrival time is known. Initially, it is the seed point specified by the user, and the time T=0 when the evolution surface passes through this point.

2) Narrow-band points: the neighborhood of all active points, whose arrival time T=1/F, and put the narrow-band points into a minimum sorting heap, where F is the curvature-weighted speed function, which controls the evolution direction and speed of the surface The important function of is calculated according to formula (2) in the present invention.

3) Far away point: the point whose arrival time is infinite, T=MAX_TIME, all the far away points except the active point and the narrowband point.

Topology-preserving surface evolution using curvature-weighted velocity functions

1) Output the point r' with the smallest arrival time T from the minimum sorting heap, mark it as an active point, and delete it from the minimum sorting heap.

2) Investigate the adjacent points of r', that is, the neighbor points of this point in space, as shown in Figure 3, for each adjacent point r, if it is an active point, it will remain unchanged; if it is a distant point, it will be modified to a narrowband point, And calculate its arrival time according to formula (1), record its parent node r', and put it into the minimum sorting heap. If it is a narrow-band point, recalculate its arrival time according to formula (1). If the arrival time T is less than the arrival time of the last iteration, update the arrival time of this point and the parent node r′, and adjust the minimum sorting heap.

During each iteration, the surface will evolve from a point r′ to r, r is the neighbor node of r′, and has the minimum arrival time in the narrow band, similarly, r′ is evolved from r″, then In the evolution model, it is considered that r″→r′→r is the topological information in the evolution process, and r″→r′→r is on the fiber path, as shown in Figure 4.

$\mathrm{<span\; class="notranslate">\; T</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; r</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; T</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; r</span>}}^{<span\; class="notranslate">\; \&prime;</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; +</span>\frac{<span\; class="notranslate">\; |</span>\mathrm{<span\; class="notranslate">\; r</span>}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; r</span>}}^{<span\; class="notranslate">\; \&prime;</span>}<span\; class="notranslate">\; |</span>}{\mathrm{<span\; class="notranslate">\; f</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; r</span>}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; .</span><span\; class="notranslate">\; .</span><span\; class="notranslate">\; .</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span>$

Among them, T(r') is the arrival time of r', and F(r) is a speed function, which ensures that the time T(r) advances from the seed point along the n(r) direction at the speed F(r).

$\mathrm{<span\; class="notranslate">\; f</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; r</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>\frac{\mathrm{<span\; class="notranslate">\; C</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; r</span>}<span\; class="notranslate">\; )</span>}{\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; min</span>}<span\; class="notranslate">\; (</span><span\; class="notranslate">\; (</span><span\; class="notranslate">\; |</span>{\mathrm{<span\; class="notranslate">\; e</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; r</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; \&CenterDot;</span>\mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; r</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; |</span><span\; class="notranslate">\; )</span><span\; class="notranslate">\; ,</span><span\; class="notranslate">\; (</span><span\; class="notranslate">\; |</span>{\mathrm{<span\; class="notranslate">\; e</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; r</span>}}^{<span\; class="notranslate">\; \&prime;</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; \&CenterDot;</span>\mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; r</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; |</span><span\; class="notranslate">\; )</span><span\; class="notranslate">\; ,</span><span\; class="notranslate">\; (</span><span\; class="notranslate">\; |</span>{\mathrm{<span\; class="notranslate">\; e</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; r</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; \&CenterDot;</span>{\mathrm{<span\; class="notranslate">\; e</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; r</span>}}^{<span\; class="notranslate">\; \&prime;</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; |</span><span\; class="notranslate">\; )</span><span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; .</span><span\; class="notranslate">\; .</span><span\; class="notranslate">\; .</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 2</span>}<span\; class="notranslate">\; )</span>$

$\mathrm{<span\; class="notranslate">\; C</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; r</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>\left\{\begin{array}{cc}\mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; r</span>}}^{<span\; class="notranslate">\; \&prime;</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; \&CenterDot;</span>\mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; r</span>}<span\; class="notranslate">\; )</span>& \mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; r</span>}}^{<span\; class="notranslate">\; \&prime;</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; \&Center\; Dot;</span>\mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; r</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; ></span>\mathrm{<span\; class="notranslate">\; 0</span>}\\ \mathrm{<span\; class="notranslate">\; 0</span>}& \mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; r</span>}}^{<span\; class="notranslate">\; \&prime;</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; \&Center\; Dot;</span>\mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; r</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; <</span>\mathrm{<span\; class="notranslate">\; 0</span>}\end{array}\right.<span\; class="notranslate">\; .</span><span\; class="notranslate">\; .</span><span\; class="notranslate">\; .</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 3</span>}<span\; class="notranslate">\; )</span>$

where r′ is the active point, n(r′) represents the evolution direction from the surface to r′, n(r) represents the evolution direction from the surface to r, e ₁ (r) represents the principal vector direction of r, e ₁ (r′ ) represents the principal vector direction of r′. The evolution direction n(r) is approximated by n(r)=|rr′|, and the retained topology n(r′)=|r′-r″| is the fiber curvature. When C(r)<0, F(r )=0, to ensure that the surface evolution always expands and never shrinks.

The curvature-weighted velocity function is based on the following four considerations:

a. When starting from the same active point and reaching different neighborhoods (different narrow-band points) at the same speed, choose the direction with the smallest curvature to advance.

b. When starting from different active points and arriving at the same neighborhood (same narrow-band point) with the same arrival time, select the active point with smaller curvature as the parent node.

c. When the direction of the fiber becomes inconspicuous and the direction of the vector field is inconsistent, the direction with smaller curvature deflection is preferred to maintain the continuity and smoothness of the fiber.

d. Ensure that the surface evolves along the direction of fiber extension, avoiding X and T-shaped fiber paths.

The curvature-weighted velocity function takes the bending energy into account, and when the evolved surface reaches the fiber boundary, the irregular vector field will be excluded because of the curvature limit, because it causes the bending energy to become larger. At the same time, inside the fiber bundle, the curvature-weighted velocity function minimizes the bending energy, enabling the fiber tracking process to advance along the direction of the tensor field, reducing false positives.

3) If the surface evolution exceeds the volume range of the volume data, stop the iteration; define a time threshold, when the arrival time of the point exceeds the threshold, it is considered that the point has a low degree of correlation with the seed point, and there is no white matter between the seed point and the seed point Fiber, stop the iteration; or pre-define the number of iterations, reach the number of iterations to stop, then the loop ends, otherwise go to 1) and continue to execute.

Determine all paths using the time gradient descent method

的累计结果。快速行进算法保证从种子点A到r的最小路径发生的演化过程中，并且T(r)就是最小代价。如果A到r存在一条路径则满足：Suppose there is a fiber γ, and suppose the length is L, τ is a point on γ, time T is under the action of velocity F

cumulative results. The fast marching algorithm guarantees the evolution process of the minimum path from the seed point A to r, and T(r) is the minimum cost. If there is a path from A to r then it satisfies:

$\mathrm{<span\; class="notranslate">\; T</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; r</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>\underset{\mathrm{<span\; class="notranslate">\; \&lambda;</span>}}{\mathrm{<span\; class="notranslate">\; min</span>}}{<span\; class="notranslate">\; \&Integral;</span>}_{\mathrm{<span\; class="notranslate">\; A</span>}}^{\mathrm{<span\; class="notranslate">\; r</span>}}<span\; class="notranslate">\; |</span><span\; class="notranslate">\; \&dtri;</span>\mathrm{<span\; class="notranslate">\; T</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; r</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; \&tau;</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; )</span><span\; class="notranslate">\; |</span>\mathrm{<span\; class="notranslate">\; d\&tau;</span>}<span\; class="notranslate">\; .</span><span\; class="notranslate">\; .</span><span\; class="notranslate">\; .</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 4</span>}<span\; class="notranslate">\; )</span>$

Then from the end of the iteration, the shortest path back to the seed point is calculated from the wavefront position through the time gradient descent method, which is a fiber path.

Use the Unicom matrix to select real fiber paths

$\mathrm{<span\; class="notranslate">\; \&phi;</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; \&gamma;</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; min</span>}\underset{\mathrm{<span\; class="notranslate">\; \&tau;</span>}}{\mathrm{<span\; class="notranslate">\; f</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; \&gamma;</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; \&tau;</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; =</span>\underset{\mathrm{<span\; class="notranslate">\; \&tau;</span>}}{\mathrm{<span\; class="notranslate">\; min</span>}}\frac{\mathrm{<span\; class="notranslate">\; 1</span>}}{<span\; class="notranslate">\; |</span><span\; class="notranslate">\; \&dtri;</span>\mathrm{<span\; class="notranslate">\; T</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; \&gamma;</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; \&tau;</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; )</span><span\; class="notranslate">\; |</span>}<span\; class="notranslate">\; .</span><span\; class="notranslate">\; .</span><span\; class="notranslate">\; .</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 5</span>}<span\; class="notranslate">\; )</span>$

After the iteration of the fast marching algorithm is over, any point with an arrival time of T can be connected to the seed point through the minimum path algorithm. The judgment of the real path needs to set a threshold value of φ through the connectivity matrix φ, map through the φ function, and the seed point The most probable paths of points with strong connectivity are selected, while those fibers whose paths do not meet the condition are considered impossible, are deleted.

The white matter fiber is reconstructed by adopting the topology-preserving fast marching method proposed by the present invention, and the reconstruction algorithm starts from the seed point. The reconstructed data come from GE MEDICAL SYSTEMS, using EP/SE sequence. The seed point is selected at the junction of the corpus callosum and the sagittal view. Figure 5 shows the results obtained by selecting the projection fiber corona and using the topology-preserving fast-tracking method of the present invention. In the figure, different fiber orientations are well distinguished. It can be seen from Figure 6 that the fiber path presents a tree structure, which well reflects the bifurcation information of the fiber, conforms to the model of the real fiber, and maintains a good topology.

Claims (1)

1. A fast-moving fiber tracking method of topology preservation, characterized in that, comprising the following steps: The first step: DTI data processing; Step 2: Seed point selection and initialization (1) Designate the initial active point, that is, the seed point r', as the starting position of the surface evolution; (2) Determine the speed function

in,

r′ is the active point, n(r′) represents the evolution direction from the surface to r′, n(r) represents the evolution direction from the surface to r, e ₁ (r) represents the principal vector direction of r, e ₁ (r′) Indicates the principal vector direction of r′, the evolution direction n(r) is approximated by n(r)=|rr′|, and the fiber curvature n(r′)=|r′-r″|; (3) Determine the narrowband point and the far away point, and determine the arrival time T of the narrowband point according to the following formula, the arrival time of the faraway point is initially T=TIME_MAXE, and set up the minimum sorting heap of the narrowband point: $T\left(r\right)=T\left(r^{\&prime;}\right)+\frac{|r-r^{\&prime;}|}{f\left(r\right)},$ where T(r') is the arrival time of r'; The third step: use the curvature-weighted velocity function to perform topology-preserving surface evolution; (1) Output the point r' with the smallest arrival time T from the minimum sorting heap, mark it as an active point, and delete it from the minimum sorting heap; (2) Investigate the adjacent point r of the active point r′. For each adjacent point r, if it is an active point, it will remain unchanged; Parent node r', and put it into the minimum sorting heap; if it is a narrowband point, recalculate its arrival time T(r), if the arrival time T(r) is less than the arrival time of the last iteration, update the point Arrival time and parent node r′, adjust the minimum sort heap; (3) If the surface evolution exceeds the volume range of the volume data, stop the iteration; set a time threshold, when the arrival time of the point exceeds the time threshold, stop the iteration; or pre-define the number of iterations, and stop when the number of iterations is reached, then the loop ends , otherwise go to (1) and continue to execute; The third step: use the time gradient descent method to determine all fiber paths; Step 4: Use the connectivity matrix to select real fiber paths.

Images