A method for cerebral vascular tracking based on streamline differential equations

Abstract

A cerebral vascular tracking method based on streamline differential equations first introduces voxel direction information into the segmentation of cerebral blood vessels, adds a columnar model with direction information to the 4D curve model, and estimates the tensor by measuring the intensity difference of multiple sampling directions; after obtaining the vascular direction information, the blood flow movement is described as fluid movement in combination with the relevant theory of fluid motion, and the vascular trajectory is described by flow field and streamlines, and the streamline differential equation is established to reconstruct the flow field distribution of the vascular direction. Constraints are imposed on the solution function through the twenty-six neighborhood voxel information and the equality of the overall blood flow inflow and outflow, and a more accurate flow field distribution model is obtained, and the directional distribution of the blood vessel is represented by the flow field distribution; finally, the starting seed point is set, and the spatial coordinates of each point in the streamline are calculated by the numerical integration method to obtain the trajectory distribution of the blood vessel. The present invention greatly reduces the workload of expert labeling and recognition, and realizes a fully automatic, accurate and fast vascular tracking method.

Description

Technical Field

[0001] The invention relates to the fields of medical imaging, image segmentation and tracking under computer graphics, and is a cerebral blood vessel tracking method based on streamline differential equations. Background Art

[0002] Accurate vascular extraction is a prerequisite for 3D visualization and quantification, which helps to better distinguish blood vessels from surrounding tissues, observe blood vessels more clearly for more accurate diagnosis, and provide more precise information for surgical planning.

[0003] Compared with the traditional geometric model method, the deep learning method has shown a strong automatic feature learning ability and is suitable for accurate vascular segmentation. The neural network classifier has shown superior classification and segmentation performance. The trained neural network classifier is equivalent to a nonlinear decision tree of a generalized matched filter. With the outstanding performance of the U-net network structure in medical image segmentation in recent years, more and more studies have focused on using the U-net network structure to achieve cerebrovascular segmentation. Many research works focus on improving the network structure and trying to mine more information from the input image. The structural design of the network focuses on improving the resolution of the network and the comprehensive utilization of existing information. Inspired by U-net, the researchers proposed a multi-scale dense CNN segmentation framework with an encoder-decoder structure for automatic segmentation of cerebrovascular vessels in DSA images. The network introduced an improved dense block to improve the extraction ability of high-level features. Based on two-dimensional images, CE-net uses dense dilated convolution and residual multi-kernel pooling to retain spatial information. In addition to improving the information extraction ability of the network, the combined use of two-dimensional and three-dimensional information is also a research hotspot. For example, D-UNet extracts two-dimensional and three-dimensional information from two downsampling paths for fusion. Researchers also designed an end-to-end DeepVessel-Net network architecture. Although it performed well on synthetic test data, there were still problems with over-segmentation and under-segmentation of blood vessels in the segmentation results of clinical TOF-MRA data. Later, researchers designed a multi-scale convolutional model to enhance the feature extraction capability by introducing dilated convolution. Although it combines two-dimensional and three-dimensional information, it cannot guarantee the spatial consistency of the final segmentation results. In view of the specific anatomical structure of blood vessels, researchers proposed a cerebrovascular segmentation method based on a generative adversarial network, which uses an attention mechanism to enhance the extraction capability of blood vessels and uses the minimum radius and connectivity of blood vessels as indicators to measure the segmentation results. In these works, most CNNs are end-to-end network structures, divided into two parts: encoder and decoder.

[0004] However, the end-to-end network structure requires a large amount of fully spaced densely labeled data for model training. Since the topological structure of cerebrovascular segmentation is complex and changeable, distributed throughout the skull, and has a small intracranial volume ratio, manually marking the complete cerebrovascular system from the scanned data slices is a time-consuming and error-prone task. Therefore, it is very difficult to obtain a set of high-precision labeled data, making deep learning in the field of cerebrovascular segmentation a great challenge. Summary of the invention

[0005] In order to overcome the biggest bottleneck of using deep learning for cerebrovascular segmentation tasks, which is the inability to obtain a large amount of full-space labeled data, the present invention provides a cerebrovascular tracking method based on streamline differential equations. It only needs to set the starting seed point to obtain the trajectory distribution of the blood vessels, which greatly reduces the workload of expert labeling and recognition, and realizes a fully automatic, accurate and fast blood vessel tracking method.

[0006] The technical solution adopted by the present invention to solve its technical problem is:

[0007] A cerebral blood vessel tracking method based on streamline differential equations comprises the following steps:

[0008] Step 1: Data preprocessing

[0009] In order to eliminate the interference of irrelevant background areas (non-vascular areas), the BET tool in FSL was used to remove the skull, and the data was standardized by zero mean and unit variance. Then, the voxel intensity of vascular calcification was set to be equal to the intensity of brain tissue to ensure the accuracy of the vascular direction obtained in step 2.

[0010] Step 2: Obtaining blood vessel direction information

[0011] A cylindrical model with directional information is added to the 4D curve model, n directions are randomly and uniformly sampled on the unit sphere, and the tensor is estimated by measuring the intensity difference of multiple sampling directions. A tensor model is established according to the anisotropy of the blood vessels, and then the least squares method is used to solve it to obtain the directional information of the blood vessels.

[0012] Step 3: Establish a vascular flow field distribution model

[0013] After obtaining the vascular direction information, combined with the relevant theories of fluid motion, the blood flow motion is described as fluid motion, the vascular trajectory is described by flow field and streamlines, and the streamline differential equation is established to reconstruct the flow field distribution of the vascular direction. Constraints are imposed on the solution function through the consistency information function within the voxel and between the twenty-six neighboring voxels, as well as the equality of the overall blood flow inflow and outflow, to obtain a more accurate flow field distribution model, and to represent the directional distribution of the blood vessels with the flow field distribution;

[0014] Step 4: Vessel Tracking

[0015] Vessel tracking can be regarded as the fitting process of particle motion trajectory. Since the particle motion streamline trajectory is composed of discrete spatial points, the fitting process uses the numerical integration method to calculate the spatial coordinates of each point in the streamline. Given a seed point position P0, the motion direction of this point is Get the following

[0016]

[0017] in P1, P2, and P3 are three transition points starting from the starting point. is the direction of the corresponding spatial point, is the direction field solved in step 3, P i+1 P i The next spatial point position is obtained by numerical integration of the starting point, and then the new spatial point is calculated with this point as the new starting point until the termination condition is met.

[0018] Furthermore, in step 2, the specific method of obtaining the blood vessel direction information is as follows:

[0019] First, the blood vessel model is constructed as follows:

[0020]

[0021] Where c(u) represents the position of the centerline of the blood vessel in three-dimensional space, r(u) represents the radius of the sphere with u as the center, and L represents the length of the blood vessel. A cylinder is established along each sampling direction, which is represented by a stack of disks as follows:

[0022]

[0023] Randomly and uniformly sample n directions on the unit sphere, use the piecewise constant intensity model, calculate the average intensity difference between the sphere Ω1 and the region Ω2 in each direction, and the region Ω2 is expressed as Ω2=cyl-Ω1. The measurement value M is calculated by the following formula i :

[0024] M i =(μΩ1-μΩ2) 2

[0025] where μΩ1 and μΩ2 represent the average intensity of regions Ω1 and Ω2, respectively.

[0026] Then construct a 3×3 symmetric positive definite matrix Measurement value M i Modeled as make

[0027]

[0028] Since D is a symmetric matrix, the six components are superimposed into a vector d

[0029] d=[d 11 d 22 d 33 d 12 d 13 d 23 ] T

[0030] So we can get M i =H i d. Finally, the least squares method is used to solve the above equation, which usually gives good results for noisy data sets because more measurements are used in the estimation process than necessary.

[0031]

[0032] The above formula can be used to calculate the eigenvalues ​​λ1, λ2, and λ3 of D and the corresponding eigenvectors e1, e2, and e3. The anisotropy of the blood vessel can intuitively describe the spatial geometry of the tensor. When λ1=λ2=λ3, the tensor is a sphere, indicating that the anisotropy of this region is the same in all directions, that is, the background region; when λ1≈λ2>>λ3, the tensor is a column, indicating that the directionality of this region is regular, and the eigenvector e3 corresponding to the smallest eigenvalue λ3 is the direction of the blood vessel.

[0033] Furthermore, in step 3, the specific method of obtaining and establishing the vascular flow field distribution model is as follows:

[0034] The blood flow in the blood vessels is regarded as a problem of fluid mechanics. According to Riemann geometry, the streamline distribution can be composed of a series of streamline sets S = {S i , i = 1, ..., n}, in the streamline set S, a blood vessel dS connects the two end sections dM and dN, M and N represent the two end areas of the blood vessel, when dS takes the limit (dM → 0, dN → 0), a single blood vessel S is obtained i .

[0035] Assuming that the streamline tangents at any point (x, y, z) in space can be represented by the flow field, the directional distribution of the blood vessels can be represented by the flow field:

[0036]

[0037] Where v(x,y,z) represents the flow field distribution, v x ,v y ,v z They represent the components of the x, y, and z axes in three-dimensional space, the vascular flow field distribution v(x, y, z) and the vascular S i The following relations are satisfied:

[0038]

[0039] The streamline differential equation is defined as follows:

[0040]

[0041] The directional information of blood vessels in each imaging voxel is represented by the local flow field distribution, which is a continuous function about (x, y, z) in space. It can represent the directional information of any point in the voxel space. A set of ternary quadratic polynomials is used to represent the flow field distribution of the voxel:

[0042] υ(x,y,z)=AC(x,y,z)

[0043] The coefficient matrix A is defined as follows:

[0044]

[0045] C(x,y,z)=[x 2 ,y 2 ,z 2 ,xy,xz,yz,x,y,z,1] T , represented by C. a 19 、a 29 and a 39 They are respectively the component values ​​of the blood vessel direction e3 obtained in the first step on the three coordinate axes.

[0046] Combined with the fact that the blood flow inflow and outflow of the whole cerebral blood vessels are equal, the continuous incompressible fluid theory is applied to describe the spatial consistency of the blood vessel trajectory by introducing the concept of fluid divergence in the tensor vector field:

[0047]

[0048] According to the theory that blood flow in the same voxel satisfies the continuous incompressible fluid, that is, divΩ=0; combining the above equations, we get:

[0049]

[0050] The above formula ensures the consistency of blood flow diffusion within the voxel. In fact, blood vessels not only maintain consistency in diffusion within the voxel, but also maintain consistency between voxels. The following uses the consistency function between vascular voxels to represent the distribution characteristics of blood vessels. Assuming that the voxel is a unit cube, N c =(c 0 ,c 1 ,…,c 5 ) represents the six voxels adjacent to the central voxel c, with A c =(A 0 ,A 1 ,…,A5 ) represents the coefficient matrix of the neighboring voxel c. For the neighboring voxels c and c that the same flow beam passes through i The streamline trajectories are S and S i , respectively, by the flow field distribution v(x,y,z)=AC(x,y,z) and v i (x,y,z)=A i C(x,y,z) represents the voxel interface through which the beam passes as Γ, and the consistency function divΓ(x,y,z) of any point on the intersection is defined as follows:

[0051]

[0052] Continuity function on the entire surface Γ That is the surface integral of divΓ(x,y,z):

[0053]

[0054] Each voxel has six neighboring voxels, so the consistency function between voxels is the sum of the continuity functions of the six interfaces:

[0055]

[0056] Combining the above equations, we can obtain Ψ.

[0057] According to the consistency function within and between voxels, combined with the anisotropic information of the blood vessels in the 26 neighborhoods of the voxel, the flow field distribution of the blood vessels is estimated by minimizing the fitting. This method of using the 26 neighborhood information to calculate the flow field distribution of the central voxel avoids the influence of the uneven intensity of a single voxel and the error of the vascular boundary. When selecting the 26 neighborhood information, the angle threshold is used to limit the smoothness of the local voxel direction information, and the voxel information with large direction changes is deleted to ensure the consistency of the local direction information. The flow field coefficient matrix A is calculated by solving the minimum cost function:

[0058]

[0059]

[0060] In order to simplify the calculation, the present invention intends to use the calculated blood vessel direction information of the neighborhood of the central voxel 26 as an approximation of Φ(x, y, z), and finally solve the above formula using the enhanced Lagrangian method to obtain the flow field coefficient matrix A and obtain the flow field distribution model of the blood flow.

[0061] The solution of the present invention first introduces voxel direction information into the segmentation of cerebral blood vessels, adds a columnar model with direction information to the 4D curve model, and estimates the tensor by measuring the intensity difference of multiple sampling directions; after obtaining the blood vessel direction information, combined with the relevant theory of fluid motion, the blood flow motion is described as fluid motion, the blood vessel trajectory is described by flow field and streamlines, and the streamline differential equation is established to reconstruct the flow field distribution in the blood vessel direction. Constraints are imposed on the solution function through the twenty-six neighborhood voxel information and the equality of the overall blood flow inflow and outflow, so as to obtain a more accurate flow field distribution model, and the direction distribution of the blood vessel is represented by the flow field distribution; finally, the starting seed point is set, and the spatial coordinates of each point in the streamline are calculated by the numerical integration method to obtain the trajectory distribution of the blood vessel.

[0062] The beneficial effects of the present invention are as follows: based on the intensity difference of blood flow in the blood vessels, a seed point is given, and a flow field distribution model is established, so that the distribution trajectory of the blood vessels can be automatically obtained, which greatly reduces the workload of expert labeling and identification, and realizes a fully automatic, accurate and fast blood vessel tracking method. BRIEF DESCRIPTION OF THE DRAWINGS

[0063] Figure 1 It is a flow chart of the implementation steps of the present invention.

[0064] Figure 2 It is a model diagram of blood vessel direction. DETAILED DESCRIPTION

[0065] The present invention will be further described below.

[0066] Reference Figure 1 and Figure 2 , a cerebral blood vessel tracking method based on streamline differential equations, comprising the following steps:

[0067] Step 1: Data preprocessing

[0068] In order to eliminate the interference of irrelevant background areas (non-vascular areas), the BET tool in FSL was used to remove the skull, and the data was standardized by zero mean and unit variance. Then, the voxel intensity of vascular calcification was set to be equal to the intensity of brain tissue to ensure the accuracy of the vascular direction obtained in step 2.

[0069] Step 2: Obtaining blood vessel direction information

[0070] A cylindrical model with directional information is added to the 4D curve model, n directions are randomly and uniformly sampled on the unit sphere, and the tensor is estimated by measuring the intensity difference of multiple sampling directions. A tensor model is established according to the anisotropy of the blood vessels, and then the least squares method is used to solve it to obtain the directional information of the blood vessels.

[0071] Furthermore, in step 2, the specific method of obtaining the blood vessel direction information is as follows:

[0072] First, the blood vessel model is constructed as follows:

[0073]

[0074] Where c(u) represents the position of the centerline of the blood vessel in three-dimensional space, r(u) represents the radius of the sphere centered at u, and L represents the length of the blood vessel. A cylinder is established along each sampling direction, represented by a stack of disks as follows

[0075]

[0076] Randomly and uniformly sample n directions on the unit sphere, use the piecewise constant intensity model, calculate the average intensity difference between the sphere Ω1 and the region Ω2 in each direction, and the region Ω2 is expressed as Ω2=cyl-Ω1. The measurement value M is calculated by the following formula i :

[0077] M i =(μΩ i μΩ2) 2

[0078] where μΩ1 and μΩ2 represent the average intensity of regions Ω1 and Ω2, respectively.

[0079] Then construct a 3×3 symmetric positive definite matrix Measurement value M i Modeled as make

[0080]

[0081] Since D is a symmetric matrix, the six components are superimposed into a vector d

[0082] d=[d 11 d 22 d 33 d 12 d 13 d 23 ] T

[0083] So we get M i =H i d. Finally, the least squares method is used to solve the above equation, which usually gives good results for noisy data sets because more measurements are used in the estimation process than necessary;

[0084]

[0085] The above formula can be used to calculate the eigenvalues ​​λ1, λ2, and λ3 of D and the corresponding eigenvectors e1, e2, and e3. The anisotropy of blood vessels can intuitively describe the spatial geometry of the tensor. When λ1=λ2=λ3, the tensor is a sphere, indicating that the anisotropy of all directions in this area is the same, that is, the background area; when λ1≈λ2>>λ3, the tensor is a column, indicating that the directionality of this area is regular, and the eigenvector e3 corresponding to the smallest eigenvalue λ3 is the direction of the blood vessel;

[0086] Step 3: Establish a vascular flow field distribution model

[0087] After obtaining the vascular direction information, combined with the relevant theories of fluid motion, the blood flow motion is described as fluid motion, the vascular trajectory is described by flow field and streamlines, and the streamline differential equation is established to reconstruct the flow field distribution of the vascular direction. Constraints are imposed on the solution function through the consistency information function within the voxel and between the twenty-six neighboring voxels, as well as the equality of the overall blood flow inflow and outflow, to obtain a more accurate flow field distribution model, and to represent the directional distribution of the blood vessels with the flow field distribution;

[0088] In step 3, the specific method of obtaining and establishing the vascular flow field distribution model is as follows:

[0089] The blood flow in the blood vessels is regarded as a problem of fluid mechanics. According to Riemann geometry, the streamline distribution can be composed of a series of streamline sets S = {S i , i = 1, ..., n}, in the streamline set S, a blood vessel dS connects the two end sections dM and dN, M and N represent the two end areas of the blood vessel, when dS takes the limit (dM → 0, dN → 0), a single blood vessel S is obtained i .

[0090] Assuming that the streamline tangents at any point (x, y, z) in space can be represented by the flow field, the directional distribution of the blood vessels can be represented by the flow field:

[0091]

[0092] Where v(x,y,z) represents the flow field distribution, v x ,v y ,v z They represent the components of the x, y, and z axes in three-dimensional space, the vascular flow field distribution v(x, y, z) and the vascular S i The following relations are satisfied:

[0093]

[0094] The streamline differential equation is defined as follows:

[0095]

[0096] The directional information of blood vessels in each imaging voxel is represented by the local flow field distribution, which is a continuous function about (x, y, z) in space. It can represent the directional information of any point in the voxel space. A set of ternary quadratic polynomials is used to represent the flow field distribution of the voxel:

[0097] υ(x ， y,z)=AC(x,y,z)

[0098] The coefficient matrix A is defined as follows:

[0099]

[0100] C(x,y,z)=[x 2 ,y 2 ,z 2 ,xy,xz,yz,x,y,z,1] T , represented by C. a 19 、a 29 and a 39 They are respectively the component values ​​of the blood vessel direction e3 obtained in the first step on the three coordinate axes.

[0101] Combined with the fact that the blood flow inflow and outflow of the whole cerebral blood vessels are equal, the continuous incompressible fluid theory is applied to describe the spatial consistency of the blood vessel trajectory by introducing the concept of fluid divergence in the tensor vector field:

[0102]

[0103] According to the theory that blood flow in the same voxel satisfies the continuous incompressible fluid, that is, divΩ=0; combining the above equations, we get:

[0104]

[0105] The above formula ensures the consistency of blood flow diffusion within the voxel. In fact, blood vessels not only maintain consistency in diffusion within the voxel, but also maintain consistency between voxels. The following uses the consistency function between vascular voxels to represent the distribution characteristics of blood vessels. Assuming that the voxel is a unit cube, N c =(c 0 ,c 1 ,…,c 5 ) represents the six voxels adjacent to the central voxel c, with A c =(A 0 ,A 1 ,…,A 5 ) represents the coefficient matrix of the neighboring voxel c. For the neighboring voxels c and c that the same flow beam passes through i The streamline trajectories are S and S i , respectively, by the flow field distribution υ(x,y,z)=AC(x,y,z) and v i(x,y,z)=A i C(x,y,z) represents the voxel interface through which the beam passes as Γ, and the consistency function divΓ(x,y,z) of any point on the intersection is defined as follows:

[0106]

[0107] Continuity function on the entire surface Γ That is the surface integral of divΓ(x,y,z):

[0108]

[0109] Each voxel has six neighboring voxels, so the consistency function between voxels is the sum of the continuity functions of the six interfaces:

[0110]

[0111] Combining the above equations, we can obtain Ψ.

[0112] According to the consistency function within and between voxels, combined with the anisotropic information of the blood vessels in the 26 neighborhoods of the voxel, the flow field distribution of the blood vessels is estimated by minimizing the fitting. By using the 26 neighborhood information to calculate the flow field distribution of the central voxel, the influence of the uneven intensity of a single voxel and the error of the vascular boundary is avoided. When selecting the 26 neighborhood information, the angle threshold is used to limit the smoothness of the local voxel direction information, and the voxel information with large direction changes is deleted to ensure the consistency of the local direction information. The flow field coefficient matrix A is calculated by solving the minimum cost function:

[0113]

[0114]

[0115] In order to simplify the calculation, the present invention intends to use the calculated blood vessel direction information of the neighborhood of the central voxel 26 as an approximation of Φ(x, y, z), and finally solve the above formula using the enhanced Lagrangian method to obtain the flow field coefficient matrix A and obtain the flow field distribution model of the blood flow.

[0116] Step 4: Vessel Tracking

[0117] Blood vessel tracking can be regarded as the fitting process of particle motion trajectory. Since the particle motion streamline trajectory is composed of discrete spatial points, the fitting process of the present invention uses the numerical integration method to calculate the spatial coordinates of each point in the streamline. Assume that the seed point position P0 is given and the motion direction of the point is Get the following

[0118]

[0119] in P1, P2, and P3 are three transition points starting from the starting point. is the direction of the corresponding spatial point, is the direction field solved in step 3. i+1 P i The next spatial point position is obtained by numerical integration of the starting point, and then the new spatial point is calculated with this point as the new starting point until the termination condition is met.

[0120] The specific implementation described above is only an optimal implementation method of the present invention and is not intended to limit the patent scope of the present invention. Any equivalent structure or equivalent process transformation made by utilizing the spirit and principles of the present invention and the contents of the drawings should be included in the patent protection scope of the present invention.

Claims

1. A cerebral blood vessel tracking method based on streamline differential equations, characterized in that: The method comprises the following steps: Step 1: Data preprocessing The skull was removed using the BET tool in FSL, and the data was standardized using zero mean and unit variance. The voxel intensity of vascular calcification was then set to be equal to that of brain tissue to ensure the accuracy of the vascular direction obtained in step 2. Step 2: Obtaining blood vessel direction information A cylindrical model with directional information is added to the 4D curve model, n directions are randomly and uniformly sampled on the unit sphere, and the tensor is estimated by measuring the intensity difference of multiple sampling directions. A tensor model is established according to the anisotropy of the blood vessels, and then the least squares method is used to solve it to obtain the directional information of the blood vessels. Step 3: Establish a vascular flow field distribution model After obtaining the vascular direction information, combined with the relevant theories of fluid motion, the blood flow motion is described as fluid motion, the vascular trajectory is described by flow field and streamlines, and the streamline differential equation is established to reconstruct the flow field distribution of the vascular direction. Constraints are imposed on the solution function through the consistency information function within the voxel and between the twenty-six neighboring voxels, as well as the equality of the overall blood flow inflow and outflow, to obtain a more accurate flow field distribution model, and to represent the directional distribution of the blood vessels with the flow field distribution; Step 4: Vessel Tracking Blood vessel tracking is regarded as the fitting process of particle motion trajectory. Since the particle motion streamline trajectory is composed of discrete space points, the fitting process uses the numerical integration method to calculate the spatial coordinates of each point in the streamline. Given a seed point position P0, the motion direction of this point is Get the following in P1, P2, and P3 are three transition points starting from the starting point. is the direction of the corresponding spatial point, is the direction field solved in step 3, P i+1 P i The next spatial point position is obtained by numerical integration of the starting point, and then the new spatial point is calculated with this point as the new starting point until the termination condition is met; In the step 3, establishing the blood vessel flow field distribution model comprises the following steps: The blood flow in the blood vessels is regarded as a problem of fluid mechanics. According to Riemann geometry, the streamline distribution is composed of a series of streamline sets S = {S i , i=1,…,n}, in the streamline set S, one blood vessel dS connects the two end sections dM and dN, M and N represent the two end areas of the blood vessel, when dS takes the limit, that is, dM→0, dN→0, a single blood vessel S is obtained i ; Assume that the streamline tangents at any point (x, y, z) in space are represented by the flow field, and the direction distribution of the blood vessels is represented by the flow field: Where v(x,y,z) represents the flow field distribution, v x ,v y ,v z They represent the components of the x, y, and z axes in three-dimensional space, the vascular flow field distribution v(x, y, z) and the vascular S i The following relations are satisfied: The streamline differential equation is defined as follows: The directional information of the blood vessels in each imaging voxel is represented by the local flow field distribution, which is a continuous function about (x, y, z) in space, representing the directional information of any point in the voxel space. A set of ternary quadratic polynomials is used to represent the flow field distribution of the voxel: υ(x,y,z)=AC(x,y,z) The coefficient matrix A is defined as follows: C(x,y,z)=[x 2 ,y 2 ,z 2 ,xy,xz,yz,x,y,z,1] T , represented by C, a 19 、a 29 and a 39 are the component values ​​of the blood vessel direction e3 obtained in the first step under the three coordinate axes; Combined with the fact that the blood flow inflow and outflow of the whole cerebral blood vessels are equal, the continuous incompressible fluid theory is applied to describe the spatial consistency of the blood vessel trajectory by introducing the concept of fluid divergence in the tensor vector field: According to the theory that blood flow in the same voxel satisfies the continuous incompressible fluid, that is, divΩ=0; combining the above equations, we get: The above formula ensures the consistency of blood flow diffusion within the voxel. In fact, the diffusion of blood vessels is not only consistent within the voxel, but also consistent between voxels. The following uses the consistency function between vascular voxels to represent the distribution characteristics of blood vessels; assuming that the voxel is a unit cube, N c =(c 0 ,c 1 ,…,c 5 ) represents the six voxels adjacent to the central voxel c, with A c =(a 0 ,A 1 ,…,A 5 ) represents the coefficient matrix of the neighborhood voxel c. For the neighboring voxels c and c that the same flow beam passes through, i The streamline trajectories are S and S i , respectively, by the flow field distribution υ(x,y,z)=AC(x,y,z) and v i (x,y,z)=A i C(x,y,z) represents the voxel interface through which the beam passes as Γ, and the consistency function divΓ(x,y,z) of any point on the intersection is defined as follows: Continuity function on the entire surface Γ That is the surface integral of divΓ(x,y,z): Each voxel has six neighboring voxels, so the consistency function between voxels is the sum of the continuity functions of the six interfaces: Combining the above equations, we can obtain ψ; According to the consistency function within and between voxels, combined with the anisotropic information of the blood vessels in the 26 neighborhoods of the voxel, the flow field distribution of the blood vessels is estimated by minimizing the fitting. By using the 26 neighborhood information to calculate the flow field distribution of the central voxel, the influence of the uneven intensity of a single voxel and the error of the vascular boundary is avoided. When selecting the 26 neighborhood information, the smoothness of the local voxel direction information is limited by the angle threshold, and the voxel information with large direction changes is deleted to ensure the consistency of the local direction information. The flow field coefficient matrix A is calculated by solving the minimum cost function: The calculated blood vessel direction information of the neighborhood of the central voxel 26 is used as an approximation of Φ(x, y, z). Finally, the enhanced Lagrangian method is used to solve the above equation to obtain the flow field coefficient matrix A and the flow field distribution model of the blood flow.

2. A cerebral blood vessel tracking method based on streamline differential equations as claimed in claim 1, characterized in that: In the step 2, the blood vessel direction information includes the following steps: First, the blood vessel model is constructed as follows: Where c(u) represents the position of the centerline of the blood vessel in three-dimensional space, r(u) represents the radius of the sphere with u as the center, and L represents the length of the blood vessel. A cylinder is established along each sampling direction, which is represented by a stack of disks as follows: Randomly and uniformly sample n directions on the unit sphere, use the piecewise constant intensity model, calculate the average intensity difference between the sphere Ω1 and the region Ω2 in each direction, and the region Ω2 is expressed as Ω2=cyl-Ω1. The measurement value M is calculated by the following formula i : M i =(μΩ1-μΩ2) 2 where μΩ1 and μΩ2 represent the average intensity of regions Ω1 and Ω2, respectively. Then construct a 3×3 symmetric positive definite matrix Measurement value M i can be modeled as make Since D is a symmetric matrix, the six components are superimposed into a vector d d=[d 11 d 22 d 33 d 12 d 13 d 23 ] T So we get M i =Hd, and finally the least square method was used to solve the above equation. The above formula calculates the eigenvalues ​​λ1, λ2, and λ3 of D and the corresponding eigenvectors e1, e2, and e3. The anisotropy of the blood vessel can intuitively describe the spatial geometry of the tensor. When λ1=λ2=λ3, the tensor is a sphere, indicating that the anisotropy of this region is the same in all directions, that is, the background region; when λ1≈λ2>>λ3, the tensor is a column, indicating that the directionality of this region is regular, and the eigenvector e3 corresponding to the smallest eigenvalue λ3 is the direction of the blood vessel.

Images