A method for simulating the diffusion of contrast agents in blood vessels

Abstract

The present invention proposes a method for simulating the diffusion of contrast agents in blood vessels, which uses an implicit boundary representation method and a two-phase mixed fluid model to simulate the diffusion movement of contrast agents in various complex blood vessels. The method includes three steps: first, a boundary volume map is generated for the blood vessel boundary; then, the simulation parameters are initialized, the properties and positions of the "blood-contrast agent" mixed particles are initialized, and the boundary volume map is loaded; finally, a simulation cycle is performed to calculate the density and force of the "blood-contrast agent" mixed particles at each moment, update their speed and position, and render. This method can realistically and in real time simulate the diffusion effect of contrast agents in various complex blood vessels, solve the problem of realism and immersion in virtual angiography surgery, and has good repeatability and no waste of materials.

Description

Technical Field

[0001] The invention relates to the fields of computer simulation modeling and virtual surgery, and in particular to a method for simulating the diffusion of intravascular contrast agents. Background Art

[0002] Cardiovascular disease is the number one cause of death in the world, accounting for more than 40% of global disease-related deaths and is the disease with the highest number of deaths in the world.

[0003] Vascular interventional surgery is the first choice for the treatment of cardiovascular diseases due to its advantages such as less trauma, quick recovery, and high treatment success rate. Angiography is an important part of vascular interventional surgery. It is used to obtain medical images of blood vessels in real time, guide doctors to send catheters to the opening of coronary arteries, and determine the location of tumors and surrounding vascular stenosis. However, vascular interventional surgery has high requirements for doctors' operations, and virtual cardiovascular surgery can provide doctors with sufficient training opportunities, greatly shorten the surgical training cycle, and has high safety and will not cause any waste of equipment. However, the authenticity and real-time nature of virtual surgery often have a greater impact on the training effect. Among them, the simulation of the angiography process is particularly important. The real and real-time angiography results can help doctors identify and quantify vascular lesions and provide correct guidance for doctors' operations, so that the training can achieve the desired effect.

[0004] Most existing contrast agent diffusion simulation methods use simple convection-diffusion equations to simulate the diffusion of contrast agents in blood flow, ignoring the interactive motion between blood and contrast agents. In addition, existing methods use particle-based methods to sample the boundaries when dealing with vascular boundary problems. This type of method cannot smoothly represent various boundaries, which will introduce unnecessary friction. This type of method requires a lot of time and memory to construct boundary particles before each run, and will also incur huge computational overhead in the process of neighbor search and force solution during simulation. Summary of the invention

[0005] In view of the shortcomings of the existing technology, the purpose of the present invention is to provide a method for simulating the diffusion of contrast agents in blood vessels, which can accurately simulate the process of contrast agent diffusion imaging in various complex blood vessels in real time.

[0006] In order to achieve the above object, the present invention adopts the following technical solutions:

[0007] A method for simulating the diffusion of an intravascular contrast agent comprises the following steps:

[0008] Step 1: Blood vessel boundary condition processing: Use the continuous implicit boundary method to represent the blood vessel wall, and use the cubic polynomial in the sparse grid to discretize the boundary volume contribution of each location in the blood vessel in advance, and store it in the boundary volume map;

[0009] Step 2: Simulation initialization: Initialize simulation parameters, initialize the properties and positions of the "blood-contrast agent" mixed particles, and load the boundary volume map;

[0010] Step 3: Simulation loop; calculate the density and force of each "blood-contrast agent" mixed particle at the current moment, update the speed and position of all particles at the current moment, and render.

[0011] Preferably, the blood vessel boundary condition processing in step 1 includes the following operations:

[0012] 1-1: Use signed distance function to represent blood vessel boundaries;

[0013] 1-2: Use a random cubic polynomial with 32 nodes to discretize the signed distance function on a regular hexahedral grid, and sample the signed distance value at each unit node;

[0014] 1-3: Eliminate the areas far from the boundary and only keep the grid cells within four smooth radii from the boundary, forming a sparse, narrow-band pattern;

[0015] 1-4: Solve the intersection volume of the support domain and the boundary region when the particle is located somewhere in the blood vessel. The intersection volume is defined by the following formula:

[0016]

[0017] Among them, N(x i ) represents particle x i The support domain of x' represents N(x i ) region, φ(x') represents the signed distance value at x', γ(φ(x')) represents a cubic expansion function, whose function variable is φ(x'), and γ is expressed by the following formula:

[0018]

[0019] Where C represents the cubic spline kernel, r is the smooth support radius, and x represents the variable of the function γ, which is the signed distance value φ(x') at x'.

[0020] 1-5: Boundary volume function V at each unit node B Sampling is performed and the 30th-order Gauss-Legendre quadrature method is used to calculate V at each node B value.

[0021] Preferably, the simulation initialization of step 2 includes the following operations:

[0022] 2-1: Initialize the total number of particles, life cycle, smoothing radius, particle size, and simulation step parameters;

[0023] 2-2: Evenly set the position of the "blood-contrast agent" mixed particles so that they fill the blood vessels, thereby simulating the original blood flow in the blood vessels before the contrast agent is injected;

[0024] 2-3: Load the boundary volume map for the boundary model.

[0025] Preferably, the simulation cycle of step 3 includes the following operations:

[0026] 3-1: Use the nearest neighbor search algorithm to find neighbors for particles;

[0027] 3-2: Query the pre-generated boundary volume map according to the current position of the particle to obtain the intersection volume of the particle support radius and the boundary;

[0028] 3-3: Solve the change in the volume fraction of each phase of the "blood-contrast agent" mixed particle and update it. The change in the volume fraction of phase k in the "blood-contrast agent" mixed particle i is defined as the following formula:

[0029]

[0030] Among them, α k is the volume fraction of phase k in the particle, and needs to satisfy u m is the mixed fluid velocity, u mk is the phase k drift velocity, drift velocity u mk Defined as the following formula:

[0031]

[0032] Among them, τ is the parameter that controls the inertia and pressure. The larger the value of τ, the greater the inertia and pressure effect. σ is the parameter that controls the diffusion speed. The larger the value of σ, the faster the diffusion speed. ki is the mass fraction of phase k of particle i, a i is the acceleration of particle i, and the pressure gradient of particle i relative to phase k With volume fraction gradient Respectively expressed as:

[0033]

[0034]

[0035] Among them, p ki is the pressure of mixed particle i phase k, p mi is the pressure of mixed particle i, W ij is the smooth kernel function, m j is the mass of mixed particle j, The mixed particle i is composed of all its neighbor particles j through the smooth kernel function W ij The interpolation density obtained by interpolation, when the mixed particle is close to the boundary, Including boundary density contribution; p ki , p mi , The calculation formula is as follows:

[0036] p ki =α ki p mi

[0037]

[0038]

[0039] Where κ is the gas constant, x* is the distance from the particle x on the boundary. i The nearest point, V B is the boundary density contribution obtained by particle i, that is, the value queried on the boundary volume map;

[0040] 3-4: Solve the acceleration of each "blood-contrast agent" mixed particle, the calculation formula is as follows:

[0041]

[0042] in, is the density of the mixed fluid, ρ k is the static density of phase k, g is the gravitational acceleration, T m is the viscous stress tensor of the mixed fluid, T Dm is the momentum transfer of interphase convection, f r is the friction force of the boundary on the particle;

[0043] 3-5: Update the speed and position of each "blood-contrast agent" mixed particle;

[0044] 3-6: Render all particles.

[0045] Compared with the prior art, the present invention has the following obvious outstanding substantive features and significant advantages:

[0046] 1. The present invention proposes a method for simulating the diffusion of contrast agents in blood vessels. The method uses a two-phase mixed fluid model to simulate the diffusion movement of contrast agents in blood vessels, and processes blood vessel boundary conditions based on an implicit boundary representation method. The method has achieved advanced results, which is of great significance for virtual surgery.

[0047] 2. The present invention uses a high-order quadrature method in advance to calculate the intersection volume of the particle support domain and the boundary area with higher precision. During simulation, only the boundary volume needs to be queried in real time for the calculation of attribute values ​​and forces. This not only reduces the boundary sampling noise, but also improves the efficiency of neighbor search and partial dynamic calculation during simulation. In addition, the method is easy to handle dynamic vascular boundary conditions and has high scalability. BRIEF DESCRIPTION OF THE DRAWINGS

[0048] Figure 1 It is the flowchart of the intravascular contrast agent diffusion simulation method.

[0049] Figure 2 This is an example diagram of the property calculation of fluid particles near the boundary.

[0050] Figure 3 This is an example diagram of the friction calculation when fluid particles approach the boundary.

[0051] Figure 4 A simulated rendering of the diffusion of contrast agent within a blood vessel. DETAILED DESCRIPTION

[0052] The preferred embodiments of the present invention are described in detail below with reference to the accompanying drawings.

[0053] Embodiment 1:

[0054] See also Figure 1 , a method for simulating the diffusion of intravascular contrast agents, comprising the following steps:

[0055] Step 1: Use a continuous implicit boundary method to represent the blood vessel wall, and use a cubic polynomial in a sparse grid in advance to discretely calculate the boundary volume contribution of each location in the blood vessel and store it in a boundary volume map; specifically, it includes:

[0056] 1-1: Use the signed distance function to represent the blood vessel boundary, that is, the signed distance value of any point x in the finite space is equal to the minimum distance from x to the boundary, and it is negative if it is within the boundary and positive if it is in the fluid domain. Therefore, the signed distance function of the blood vessel boundary can be expressed by the following formula:

[0057]

[0058] 1-2: Use a random cubic polynomial with 32 nodes to discretize the signed distance function on a regular hexahedral grid, and sample the signed distance value at each unit node;

[0059] 1-3: Eliminate the area where |φ|>2r, that is, the grid cells with a width of four smooth radii inside and outside the boundary, forming a sparse, narrow-band pattern to reduce the memory usage of the boundary volume map and improve query efficiency. ;

[0060] 1-4: Solve the intersection volume of the support domain and the boundary region when the particle is located somewhere in the blood vessel. The intersection volume can be defined as the following formula:

[0061]

[0062] Among them, N(x i ) represents particle x i The support domain of x' represents N(x i ) region, φ(x') represents the signed distance value at x', γ(φ(x')) represents a cubic expansion function, whose function variable is φ(x'), and γ is expressed by the following formula:

[0063]

[0064] Where C represents the cubic spline kernel, r is the smooth support radius, and x represents the variable of the function γ, which is the signed distance value φ(x') at x'.

[0065] 1-5: Boundary volume function V at each unit node B Sampling is performed and the 30th-order Gauss-Legendre quadrature method is used to calculate V at each node B value.

[0066] Step 2: Simulation initialization. Initialize the properties and positions of the "blood-contrast agent" mixed particles and load the boundary volume map for the boundary model; specifically include:

[0067] 2-1: Set the total number of particles to 114864, the life cycle to 100, the smoothing radius to 0.025, the particle size to 0.00775, and the simulation step to 0.001;

[0068] 2-2: The positions of the mixed particles of "blood-contrast agent" are evenly set so that they fill the blood vessels, and the interval between particles in the blood vessels is 0.03, so as to simulate the original blood flow in the blood vessels before the contrast agent is injected;

[0069] 2-3: Check whether a boundary volume map corresponding to the current blood vessel boundary and simulation parameters exists. If not, generate the boundary volume map using step 1. If it exists, load the corresponding boundary volume map.

[0070] Step 3: Simulation cycle. Calculate the density and force of each "blood-contrast agent" mixed particle at the current moment, update the speed and position of all particles at the current moment, and render; specifically include:

[0071] 3-1: The nearest neighbor search algorithm is used to search for neighbors of particles. The simulation space is divided into uniform grids, and then hash mapping is performed according to the particle position to obtain the corresponding spatial grid index. The neighbor particles of the particle can be determined through the starting and ending positions of the grid where the particle is located and the 26 grids around it.

[0072] 3-2: Query the pre-generated boundary volume map according to the current position of the particle to obtain the intersection volume of the particle support radius and the boundary;

[0073] First, the local coordinates of the particles in the blood vessel are calculated according to the following formula:

[0074] x=R T *(x i -t)

[0075] Among them, R and t are the current rotation and center of mass coordinates of the blood vessel, respectively. i is the coordinate value of the fluid particle in the world coordinate.

[0076] Secondly, the following shape function for interpolation can be generated based on the grid sampling points in the boundary volume map:

[0077]

[0078] Among them, ξ t =(ξ t ,η t ,ζ t ) T represents the node position in the isoparametric space, and (ξ,η,ζ) represents the coordinates in the isoparametric space.

[0079] Finally, the boundary volume is interpolated using the following formula:

[0080]

[0081] T(x)=[diag(ba)] -1 (2x-(b+a))

[0082] Among them, V Bt =V B (x t ) represents the sample value calculated in advance. Formula T(x) is used to map the coordinates in the world space to the coordinates in the isoparametric space of the element, and a and b represent the minimum and maximum limit vectors in the hexahedral unit, respectively.

[0083] 3-3: Solve the change in the volume fraction of each phase of the "blood-contrast agent" mixed particle and update it. The change in the volume fraction of phase k in the "blood-contrast agent" mixed particle i can be defined as the following formula:

[0084]

[0085] Among them, α k is the volume fraction of phase k in the particle, and needs to satisfy u m is the mixed fluid velocity, u mk is the phase k drift velocity, drift velocity u mk It can be defined as the following formula:

[0086]

[0087] Among them, τ is the parameter that controls the inertia and pressure. The larger the value of τ, the greater the inertia and pressure effect. σ is the parameter that controls the diffusion speed. The larger the value of σ, the faster the diffusion speed. ki is the mass fraction of phase k of particle i, a i is the acceleration of particle i, and the pressure gradient of particle i relative to phase k With volume fraction gradient Respectively expressed as:

[0088]

[0089]

[0090] Among them, p ki is the pressure of mixed particle i phase k, p mi is the pressure of mixed particle i, W ij is the smooth kernel function, m j is the mass of mixed particle j, The mixed particle i is composed of all its neighbor particles j through the smooth kernel function W ij The interpolation density obtained by interpolation, when the mixed particle is close to the boundary, Including boundary density contribution. p ki , p mi , The calculation formula is as follows:

[0091] p ki =α ki p mi

[0092]

[0093]

[0094] Where κ is the gas constant, x* is the distance from the particle x on the boundary. i The nearest point, V B is the boundary density contribution obtained by particle i, that is, the value queried on the boundary volume map.

[0095] 3-4: Solve the acceleration of each "blood-contrast agent" mixed particle, the calculation formula is as follows:

[0096]

[0097] in, is the density of the mixed fluid, ρ k is the static density of phase k, g is the gravitational acceleration, T m is the viscous stress tensor of the mixed fluid, T Dm is the momentum transfer of interphase convection, f r is the friction force of the boundary on the particle. Solve the right side of the equation:

[0098]

[0099]

[0100]

[0101] Where x* is the distance from particle x on the boundary i The nearest point, V B is the intersection volume of the support domain of particle i and the boundary region, which can be used to calculate the boundary density contribution, h is the smooth radius domain of the particle, μ i =∑ k α k μ k is the viscosity of each phase in particle i k The polymerization viscosity.

[0102] Friction force f on the fluid particles from the boundary r It is modeled as a viscous force to handle energy dissipation near the boundary. Two orthogonal vectors t1 and t2 are defined on the tangent plane of the boundary surface, which are perpendicular to the normal vector. Then four points x are defined on the plane. 1 / 2 =x*±dt1 and x 3 / 4 =x*±dt2, set the distance d=0.5h, and the boundary friction can be defined as follows:

[0103]

[0104] Among them, s is the spatial dimension, which is set to 3, V j =0.25V B , the speed v at the four points jIt can be obtained by the following formula:

[0105] v j =v rb +ω rb ×(x j -x rb )

[0106] Among them, x rb , v rb and ω rb are the center point position, velocity and angular velocity of the blood vessel boundary respectively.

[0107] 3-5: For each "blood-contrast agent" mixed particle, update its position and velocity using the following formula:

[0108] x next =x cur +vΔt

[0109] v next =v cur +aΔt

[0110] Where a is the particle acceleration calculated in step 3-4.

[0111] 3-6: Use Billboard to render the contrast agent particles. The rendered map is a semi-transparent map that changes from dark to light from the middle to the edge. In addition, the transparency of the map is determined by the calculated volume fraction. The larger the phase volume fraction of the contrast agent in the mixed particles, the darker the rendering color of the map.

[0112] Embodiment 2:

[0113] In step 3, 3-3 needs to calculate the interpolation density of each "blood-contrast agent" mixed particle. When the particle is close to the boundary, the density calculation area of ​​the particle is composed of the fluid domain and the boundary domain, such as Figure 2 As shown. The fluid domain includes particles x i Neighboring fluid particles within the smooth radius domain, these neighboring particles are weighted by a distance-based smooth kernel function to calculate particle x i The boundary domain is the property value of particle x i The intersection area of ​​the smooth radius domain and the boundary, the volume of the intersection area V B The value can be obtained through particle x i The density of the mixed particles can be calculated by the following formula:

[0114]

[0115] Where κ is the gas constant, V B is the boundary density contribution obtained by particle i, and x* is the distance on the boundary from particle xi The closest point is used to represent the bounding volume.

[0116] Embodiment three:

[0117] In the friction calculation of step 3-4 in step 3, the friction force f of the boundary on the fluid particles is r is modeled as a viscous force to handle energy dissipation near the boundary. Figure 3 As shown, two orthogonal vectors t1 and t2 perpendicular to the normal vector are defined on the tangent plane of the boundary surface, and then four points x are defined on the plane. 1 / 2 =x*±dt1 and x 3 / 4 =x*±dt2, set the distance d=0.5h, and the boundary friction can be defined as follows:

[0118]

[0119] Among them, s is the spatial dimension, which is set to 3, V j =0.25V B , the speed v at the four points j It can be obtained by the following formula:

[0120] v j =v rb +ω rb ×(x j -x rb )

[0121] Among them, x rb , v rb and ω rb are the center point position, velocity and angular velocity of the blood vessel boundary respectively.

[0122] Embodiment 4:

[0123] This embodiment uses Billboard to render contrast agent particles. The rendering map used is a semi-transparent map that changes from dark to light from the middle to the edge. The transparency of the map is determined by the calculated volume fraction. The larger the phase volume fraction of the contrast agent in the mixed particles, the darker the rendering color of the map. Figure 4 The rendering results of this method are shown in Figure 2. The results show that the two-phase mixed fluid model and implicit boundary representation method used in this method can realistically simulate the diffusion of contrast agents in blood vessels and greatly improve the efficiency of simulation, which is of great significance for virtual surgery.

[0124] The above describes the embodiments of the present invention in conjunction with the accompanying drawings, but the present invention is not limited to the above embodiments. Various changes can be made according to the purpose of the invention. Any changes, modifications, substitutions, combinations or simplifications made according to the spirit and principle of the technical solution of the present invention should be equivalent replacement methods. As long as they meet the purpose of the invention and do not deviate from the technical principles and inventive concepts of the present invention, they belong to the protection scope of the present invention.

Claims

1. A method for simulating the diffusion of contrast agents in blood vessels, characterized in that: The following steps are included: Step 1: Blood vessel boundary condition processing: Use the continuous implicit boundary method to represent the blood vessel wall, and use the cubic polynomial in the sparse grid to discretize the boundary volume contribution of each location in the blood vessel in advance, and store it in the boundary volume map; Step 2: Simulation initialization: Initialize simulation parameters, initialize the properties and positions of the "blood-contrast agent" mixed particles, and load the boundary volume map; Step 3: Simulate the loop; calculate the density and force of each "blood-contrast agent" mixed particle at the current moment, update the speed and position of all particles at the current moment, and render; The blood vessel boundary condition processing in step 1 includes the following operations: 1-1: Use signed distance function to represent blood vessel boundaries; 1-2: Use a random type cubic polynomial with 32 nodes to discretize the signed distance function on a regular hexahedral grid, and sample the signed distance value at each unit node; 1-3: Eliminate the areas far from the boundary and only keep the grid cells within four smooth radii from the boundary, forming a sparse, narrow-band pattern; 1-4: Solve the intersection volume of the support domain and the boundary region when the particle is located somewhere in the blood vessel. The intersection volume is defined by the following formula: Among them, N(x i ) represents particle x i The support domain of x' represents N(x i ) region, φ(x') represents the signed distance value at x', γ(φ(x')) represents a cubic expansion function, whose function variable is φ(x'), and γ is expressed by the following formula: Where C represents the cubic spline kernel, r is the smooth support radius, and x represents the variable of the function γ, which is the signed distance value φ(x') at x' here; 1-5: Boundary volume function V at each unit node Β Sampling is performed and the 30th-order Gauss-Legendre quadrature method is used to calculate V at each node Β value; The simulation cycle of step 3 includes the following operations: 3-1: Use the nearest neighbor search algorithm to find neighbors for particles; 3-2: Query the pre-generated boundary volume map according to the current position of the particle to obtain the intersection volume of the particle support radius and the boundary; 3-3: Solve the change in the volume fraction of each phase of the "blood-contrast agent" mixed particle and update it. The change in the volume fraction of phase k in the "blood-contrast agent" mixed particle i is defined as the following formula: Among them, α k is the volume fraction of phase k in the particle, and needs to satisfy u m is the mixed fluid velocity, u mk is the phase k drift velocity, drift velocity u mk Defined as the following formula: Among them, τ is the parameter that controls the inertia and pressure. The larger the value of τ, the greater the inertia and pressure effect. σ is the parameter that controls the diffusion speed. The larger the value of σ, the faster the diffusion speed. ki is the mass fraction of phase k of particle i, a i is the acceleration of particle i, and the pressure gradient of particle i relative to phase k With volume fraction gradient Respectively expressed as: Among them, p ki is the pressure of mixed particle i phase k, p mi is the pressure of mixed particle i, W ij is the smooth kernel function, m j is the mass of mixed particle j, The mixed particle i is composed of all its neighbor particles j through the smooth kernel function W ij The interpolation density obtained by interpolation, when the mixed particle is close to the boundary, Including boundary density contribution; p ki , p mi , The calculation formula is as follows: p ki =α ki p mi Where κ is the gas constant, x* is the distance from the particle x on the boundary. i The nearest point, V Β is the boundary density contribution obtained by particle i, that is, the value queried on the boundary volume map; 3-4: Solve the acceleration of each "blood-contrast agent" mixed particle, the calculation formula is as follows: in, is the density of the mixed fluid, ρ k is the static density of phase k, g is the gravitational acceleration, T m is the viscous stress tensor of the mixed fluid, T Dm is the momentum transfer of interphase convection, f r is the friction force of the boundary on the particle; 3-5: Update the speed and position of each "blood-contrast agent" mixed particle; 3-6: Render all particles.

2. The method for simulating intravascular contrast agent diffusion according to claim 1, characterized in that: The simulation initialization of step 2 includes the following operations: 2-1: Initialize the total number of particles, life cycle, smoothing radius, particle size, and simulation step parameters; 2-2: Evenly set the position of the "blood-contrast agent" mixed particles so that they fill the blood vessels, thereby simulating the original blood flow in the blood vessels before the contrast agent is injected; 2-3: Load the boundary volume map for the boundary model.

Images