LBM smoke simulation method based on GPU acceleration

Abstract

The invention discloses an LBM smoke simulation method based on GPU acceleration. The LBM smoke simulation method comprises the following steps: 1) constructing an organ soft tissue and a surgical instrument model; 2) carrying out storage space distribution, data structure definition, kernel function definition and the like; 3) creating and initializing a particle entity of a smoke effect in the CPU, and distributing a required space on the CPU side for a buffer used by Compute Shaders; 4) calling a defined calculation shader interface, and binding a buffer to kernel functions of collision, flow and boundary condition steps of the LBM method, which correspond to collision processing, flow processing and boundary condition processing respectively; 5) starting a cyclic calculation process on the CPU, calling a Compute Shader code, and executing parallel calculation on the GPU; and 6) clearing the memory space distributed for the buffer area on the CPU. According to the method, the parallelism of the LBM method is utilized, the Compute Shader is used for carrying out parallel computing, the LBM equation is solved through numerical values in a kernel function binding mode, the operation time is greatly shortened, and the method has high real-time performance and authenticity.

Description

Technical Field

[0001] The invention belongs to the technical field of virtual reality medical treatment, and in particular relates to a simulation method of smoke simulation in a virtual surgery simulation system. Background Art

[0002] Electrocoagulation hemostasis operation is one of the most common surgical training links in surgical operations. In clinical practice, high-frequency electrocoagulation devices are generally used to cauterize wounds to stop bleeding. This process is inevitably accompanied by the generation of electrocoagulation smoke.

[0003] Smoke simulation can be roughly divided into methods based on graphic synthesis and methods based on physics. In the early days of smoke simulation, due to the limitation of computer performance, fractal textures and process textures based on image synthesis were widely used. However, the smoke movement simulated by these methods is not flexible enough and it is easy to see the defects of artificial traces. In response to these shortcomings, researchers have proposed a discretized smoke simulation method based on physical properties, which simulates the entire smoke movement process more realistically by increasing the description of the physical properties of the smoke. Discrete simulation methods based on physics can be roughly divided into two categories: grid method and particle method. The grid method refers to treating smoke as a discrete grid structure, solving partial differential equations and calculating the motion state of smoke by using the grid, thereby updating the density, velocity and other parameters on the grid to simulate the movement and diffusion process of smoke in real time. The particle method treats the particles in the smoke as discrete particles, and constructs the appearance of the smoke by simulating the interaction and movement between particles.

[0004] At present, there are still some problems in smoke simulation, which can be divided into two parts. The first part is the real-time problem of smoke simulation. The real smoke process includes a series of states such as smoke generation, rise, collision with objects, diffusion and disappearance. The whole smoke simulation process has a large amount of calculation and low real-time performance. In view of the complex calculation problem of smoke simulation, Liao et al. (Liao X, Si W, Yuan Z, et al. Animating wall-bounded turbulent smoke via filament-mesh particle-particle method [J]. IEEE Transactions on Visualization and Computer Graphics, 2017, 24 (3): 1260-1273.) proposed a vortex filament-mesh particle-particle (FMPP) method for turbulent smoke simulation based on fast boundary wall to synthesize the lost small-scale vortex structure of smoke and realize the smoke-obstacle interaction and vortex shedding phenomenon at the obstacle boundary in a lightweight way. Wen et al. (Wen J, Ma H. Real-time smoke simulation based on vorticity preserving lattice Boltzmann method [J]. The Visual Computer, 2019, 35: 1279-1292) proposed a vortex particle lattice Boltzmann method based on the lattice Boltzmann method (LBM) to simulate the high-resolution movement of smoke in real time. This method has good parallel computing capabilities to ensure efficiency on parallel computing devices such as GPUs. The second part is the authenticity of smoke simulation. While meeting real-time requirements, the authenticity of the simulation still needs to be improved, including the study of physical properties such as object collision and diffusion in the smoke process. In response to this, Frank et al. proposed a bidirectional coupling simulation framework, using the SPH method to simulate diffusion areas such as spray, and achieved good authenticity in simulating smoke.

[0005] Since the currently mentioned technologies cannot perfectly meet the requirements of smoke simulation in virtual surgery systems, it is necessary to study a suitable method that can have good real-time performance while ensuring authenticity. Summary of the invention

[0006] The purpose of the present invention is to provide a LBM smoke simulation method based on GPU acceleration to improve the stability and authenticity of smoke effects.

[0007] The present invention is achieved through the following technical solutions.

[0008] The present invention provides a GPU-accelerated LBM smoke simulation method, comprising the following steps:

[0009] Step 1: Build organ soft tissue and surgical instrument models. First, collect and input diagnostic images such as CT, MRI scan images and related video data. Then, use the Blender 3D modeling tool to process them to generate OBJ format stereo images of soft tissue models and surgical tools.

[0010] Step 2, before using GPU acceleration, make the following preparations: specify the call interface of the compute shader, define the data type required in the calculation process, pre-allocate the storage space required for the calculation to store the global constants and particle data in the calculation process, and define the kernel function required for the calculation. When the above preparations are completed, the Compute Shaders can be automatically imported into the Unity project used, and then the calculation can be performed using C# scripts.

[0011] Step 3: Create and initialize the particle entities of the smoke effect on the CPU, reserve CPU memory for these particle systems, and allocate the required space on the CPU side for the buffer used by Compute Shaders.

[0012] Step 4, call the defined compute shader interface and bind the buffer to the kernel functions of the collision, flow and boundary condition steps of the LBM method, which correspond to collision processing, flow processing and boundary condition processing respectively.

[0013] The formula for local collision is:

[0014] f i * (x,t)=f i (x,t)+Ω(f i )

[0015] where f i * (x,t) and f i (x, t) are the distribution functions before and after the collision in the i direction, Ω(f i ) is the collision operator, and a single relaxation time collision operator is generally used. Its formula is as shown, τ is the single relaxation time, which satisfies v is the viscosity coefficient of fluid kinematics.

[0016] The flow formula is:

[0017] fi (x+s i Δt,t+Δt)=f i * (x,t)

[0018] where f i (x+s i Δt,t+Δt) represents the new distribution function in the i direction after the next time step. Where Δt is the elapsed time between two time steps, s i is the direction vector to the neighboring cell.

[0019] Once the equilibrium distribution function is determined, macroscopic variables such as density ρ and velocity u can be calculated using simple summations:

[0020]

[0021] In the processing of boundary conditions, since the standard bounce boundary condition can only obtain first-order accuracy, for interactive smoke simulation, there is a curved boundary that cannot adapt to any unit point or midpoint, and the half-bounce boundary cannot adapt to the requirements. Therefore, the present invention adopts the bounce boundary condition. This curved boundary condition can be applied to arbitrarily complex obstacles, even if the boundary is not attached to any node or the midpoint of the node.

[0022] As attached Figure 2 The following is a schematic diagram of the curved boundary condition. There is a boundary curve between the grid points, where P f′ is the adjacent node of the fluid node, P f Point is the boundary fluid point, P s Point is a solid point, P w is the physical boundary intersection point between solid and fluid. a For the fluid point P f Point to solid point P s The speed, e a- It is the reverse speed, that is, the speed from the solid point to the fluid point. At time t, from point P f Flow to point P s The distribution function f i (p f ) is known, but since the solid point P s Flow to fluid point P f The distribution function f i (p s ) is unknown, the boundary condition of the curve is to use the known distribution function f i (p f ) to solve the unknown distribution function f i (p s ).

[0023] Assume that from the fluid point Pf To solid point P s The distance length is a unit 1, then μ=|p f -p w |(0<μ<1) represents the fluid point P f To the physical boundary P w The distance to the fluid point P f To solid point P s The ratio of the total length of the connection line can be obtained by interpolation calculation:

[0024] f i (p w )=μf i (p s )+(1-μ)f i (p f )

[0025] At the same time, a rebound boundary condition is applied to the wall node P w :

[0026]

[0027] Finally, the boundary conditions can be expressed as:

[0028]

[0029] After applying the position constraint, if the distance between two sliding constraint points is still less than the safety threshold L0, any one of them will be deleted to avoid redundant sliding constraint points from interfering with the stitching interaction simulation process. w is the physical boundary point P w The speed at which the

[0030] Step 5: Start the loop calculation process on the CPU, call the Compute Shader code, and use the functions it provides to perform parallel calculations on the GPU. Then, through the defined Compute Shader interface, update the calculation results on the CPU and send the smoke particle information processed by the GPU back to the CPU memory.

[0031] Step 6: After all calculations are completed, clean up the memory space allocated for the buffer on the CPU to avoid memory leaks. This process is completed by referencing the cleanup interface provided by Compute Shader.

[0032] Compared with the prior art, the present invention has the following beneficial effects:

[0033] The present invention provides a smoke simulation method for virtual surgery. Due to the natural parallelism of the LBM method, Compute Shader is used for parallel computing. By binding the kernel function, the LBM equation is numerically solved, which greatly shortens the time required for the calculation. In addition, the experimental results of the comparison with the existing methods show that the proposed smoke simulation method has strong real-time and authenticity. BRIEF DESCRIPTION OF THE DRAWINGS

[0034] Figure 1 This is a schematic diagram of the D3Q15 (D stands for Dimension, Q stands for Quantities of Velocity Directions) solution.

[0035] Figure 2 Schematic diagram of the curved boundary condition.

[0036] Figure 3 This is the smoke simulation effect. (a) is real smoke, and (b) is simulated smoke.

[0037] Figure 4 The smoke simulation effect in virtual surgery, where (a) the hemostasis operation begins and (b) the electrocoagulation smoke is generated.

[0038] Figure 5 It is a flow chart of the present invention. DETAILED DESCRIPTION

[0039] The present invention will be further described through the following embodiments in conjunction with the accompanying drawings so that those skilled in the art can implement the invention with reference to the description.

[0040] The smoke simulation method in a virtual surgery described in this embodiment includes the following steps:

[0041] Step 1: Construct organ soft tissue and surgical instrument models. The virtual surgery simulation system of this embodiment is composed of a liver soft tissue model and surgical instruments. The liver soft tissue model is composed of 4384 triangular basic units, with a length, height and width of 41cm, 34cm, and 34cm respectively. The surgical instrument is composed of 3303 triangular basic units, with a length, height and width of 25cm, 2cm, and 2cm respectively.

[0042] Step 2: Prepare in advance.

[0043] 1) Use the #pragma compilation directive to specify the calling interface of the compute shader.

[0044] 2) Define the specific data types needed in the calculation process so that data can be easily calculated and stored when the function is called.

[0045] 3) Define the storage space required for calculation to store global constants and particle data during the calculation process.

[0046] 4) Define the kernel function required for calculation and assign the corresponding thread group vector to it. The kernel function is the corresponding equations in the LBM method.

[0047] Step 3: Create and initialize the particle entities of the smoke effect on the CPU, reserve CPU memory for these particle systems, and allocate the required space on the CPU side for the buffer used by Compute Shaders. Construct a 100×100×300 grid in the center of the simulation window, with a side length of 10 for each small grid and an initial vortex of 0.5 (s- 1 ), the buoyancy of the particles is 0.35 N. 3000 vortex particles are placed at the center of the coordinate axis, and the vortex particles are randomly emitted along the vertical direction ±15°, with an emission speed of 500 vortex particles per second.

[0048] Step 4, call the defined compute shader interface and bind the buffer to the kernel functions of the collision, flow and boundary condition steps of the LBM method, which correspond to collision processing, flow processing and boundary condition processing respectively.

[0049] The formula for local collision is:

[0050] f i * (x,t)=f i (x,t)+Ω(f i )

[0051] where f i * (x,t) and f i (x, t) are the distribution functions before and after the collision in the i direction, Ω(f i ) is the collision operator, and a single relaxation time collision operator is generally used. Its formula is as shown, τ is the single relaxation time, which satisfies ν is the viscosity coefficient of fluid kinematics.

[0052] The flow formula is:

[0053] f i (x+s i Δt,t+Δt)=f i * (x,t)

[0054] where f i (x+s iΔt,t+Δt) represents the new distribution function in the i direction after the next time step. Where Δt is the elapsed time between two time steps, s i is the direction vector to the neighboring cell.

[0055] When the D3Q15 scheme is used, the macroscopic NS equations are expanded using the Chapman–Enskog expansion to derive the corresponding momentum balance distribution function of the LBM, namely:

[0056]

[0057]

[0058] where w i is the weight factor of the D3Q15 scheme proposed by Qian et al. in 1992. Figure 1 As shown, c s =δx / δt, which is the speed of the grid, and is generally taken as ρ is the macroscopic density and u is the macroscopic flow velocity.

[0059] Once the equilibrium distribution function is obtained, macroscopic variables such as density ρ and velocity u can be calculated using simple summation:

[0060]

[0061] In the processing of boundary conditions, since the standard bounce boundary condition can only obtain first-order accuracy, for interactive smoke simulation, there are curved boundaries that cannot adapt to any unit point or midpoint, and the half-bounce boundary cannot adapt to the requirements. Therefore, this paper adopts the bounce boundary condition. This curved boundary condition can be applied to arbitrarily complex obstacles, even if the boundary is not attached to any node or the midpoint of the node.

[0062] As attached Figure 2 The following is a schematic diagram of the curved boundary condition. There is a boundary curve between the grid points, where P f′ is the adjacent node of the fluid node, P f Point is the boundary fluid point, P s Point is a solid point, P w is the physical boundary intersection point between solid and fluid. a For the fluid point P f Point to solid point P s The speed, e a- It is the reverse speed, that is, the speed from the solid point to the fluid point. At time t, from point P f Flow to point P s The distribution function f i (p f ) is known, but since the solid point Ps Flow to fluid point P f The distribution function f i (p s ) is unknown, the boundary condition of the curve is to use the known distribution function f i (p f ) to solve the unknown distribution function f i (p s ).

[0063] Assume that from the fluid point P f To solid point P s The distance length is a unit 1, then μ=|p f -p w |(0<μ<1) represents the fluid point P f To the physical boundary P w The distance to the fluid point P f To solid point P s The ratio of the total length of the connection line can be obtained by interpolation calculation:

[0064] f i (p w )=μf i (p s )+(1-μ)f i (p f )

[0065] At the same time, a rebound boundary condition is applied to the wall node P w :

[0066]

[0067] Then interpolate it:

[0068]

[0069] Finally, the boundary conditions can be expressed as:

[0070]

[0071] After applying the position constraint, if the distance between two sliding constraint points is still less than the safety threshold L0, any one of them will be deleted to avoid redundant sliding constraint points from interfering with the stitching interaction simulation process. w is the physical boundary point P w The speed at which the

[0072] Step 5: Start the loop calculation process on the CPU, call the Compute Shader code, and use the functions it provides to perform parallel calculations on the GPU. Then, through the defined Compute Shader interface, update the calculation results on the CPU and send the smoke particle information processed by the GPU back to the CPU memory.

[0073] Step 6: After all calculations are completed, clean up the memory space allocated for the buffer on the CPU to avoid memory leaks. This process is completed by referencing the cleanup interface provided by Compute Shader.

[0074] Table 1 shows the average time consumption of the present invention and the traditional LBM algorithm without GPU acceleration under different numbers of smoke particles.

[0075] Table 1

[0076]

[0077] As attached Figure 3 is the effect of smoke interaction simulation process, where Figure 3 (a) is the real smoke state in natural environment; Figure 3 (b) is the simulated smoke effect of the present invention. It can be seen that the present invention can effectively simulate the details of real smoke with a high degree of authenticity. Figure 4 For the smoke simulation effect in the virtual surgery system, Figure 4 (a) Initial state of the hemostasis operation. Figure 4 (b) Effect of electrocoagulation smoke generation during surgery.

Claims

1. A GPU-accelerated LBM smoke simulation method, characterized by: The following steps are involved: Step 1: Build organ soft tissue and surgical instrument models. First, collect and input diagnostic images, including CT, MRI scan images and related video data. Then, use Blender 3D modeling tools to process them to generate OBJ format stereo images of soft tissue models and surgical tools. Step 2: Before using GPU acceleration, make the following preparations: specify the calling interface of the compute shader, define the data type to be used in the calculation process, pre-allocate the storage space required for the calculation to store the global constants and particle data in the calculation process, and define the kernel function required for the calculation; after the above preparations are completed, the compute shaders are automatically imported into the Unity project used, and then the calculation is performed using C# scripts; Step 3: Create and initialize the particle entities of the smoke effect on the CPU, reserve CPU memory for these particle systems, and allocate the required space on the CPU side for the buffer used by Compute Shaders; Step 4: Call the defined compute shader interface and bind the buffer to the kernel functions of the collision, flow and boundary condition steps of the LBM method, which correspond to collision processing, flow processing and boundary condition processing respectively; The formula for local collision is: f i * (x,t)=f i (x,t)+Ω(f i ) Among them, f i * (x,t) and f i (x, t) are the distribution functions before and after the collision in the i direction, Ω(f i ) is the collision operator, and a single relaxation time collision operator is generally used. Its formula is as shown, τ is the single relaxation time, which satisfies ν is the viscosity coefficient of fluid kinematics; The flow formula is: f i (x+s i Δt,t+Δt)=f i * (x,t) Among them, f i (x+s i Δt,t+Δt) represents the new distribution function in the i direction after the next time step; Δt is the elapsed time between two time steps, s i is the direction vector to the adjacent cell; Once the equilibrium distribution function is determined, the macroscopic variables are calculated using a simple summation where the density ρ and velocity u are: In the treatment of boundary conditions, the rebound boundary conditions are adopted: Assume that from the fluid point P f To the solid point P s The distance length is a unit 1, then μ=|p f -p w |,0<μ<1, indicating the fluid point P f To the physical boundary P w The distance to the fluid point P f To the solid point P s The ratio of the total length of the connection line can be obtained by interpolation calculation: f i (p w )=μf i (P s )+(1-μ)f i (p f ) At the same time, a rebound boundary condition is applied to the wall node P w : Finally, the boundary conditions can be expressed as: After applying the position constraint, if the distance between two sliding constraint points is still less than the safety threshold L0, any one of them will be deleted to avoid redundant sliding constraint points from interfering with the stitching interaction simulation process. w is the physical boundary point P w The speed at which Step 5: Start the loop calculation process on the CPU, call the Compute Shader code, and use the functions it provides to perform parallel calculations on the GPU; then update the calculation results on the CPU through the defined Compute Shader interface, and return the smoke particle information processed by the GPU to the CPU memory; Step 6: After all calculations are completed, clean up the memory space allocated for the buffer on the CPU.

Images