A GPU-accelerated LBM smoke simulation method

Abstract

A GPU-accelerated LBM smoke simulation method includes: 1) constructing organ soft tissue and surgical instrument models; 2) allocating storage space, defining data structures, and defining kernel functions; 3) creating and initializing particle entities for the smoke effect on the CPU, and allocating the necessary space for the buffer used by Compute Shaders on the CPU side; 4) calling the defined computation shader interface and binding the buffer to the kernel functions of the collision, flow, and boundary condition steps of the LBM method, which correspond to collision handling, flow handling, and boundary condition handling, respectively; 5) starting the loop calculation process on the CPU, calling the Compute Shader code, and executing parallel computation on the GPU; and 6) clearing the memory space allocated to the buffer on the CPU. This invention utilizes the parallelism of the LBM method, performs parallel computation using Compute Shaders, and numerically solves the LBM equations by binding kernel functions, significantly shortening the computation time and exhibiting strong real-time performance and realism.

Description

Technical Field

[0001] This invention belongs to the field of medical technology of virtual reality, and specifically relates to a simulation method for smoke simulation in a virtual surgical simulation system. Background Technology

[0002] Electrocoagulation hemostasis is one of the most common surgical training procedures in surgery. Clinically, high-frequency electrocoagulation devices are generally used to burn and stop bleeding in wounds, which inevitably involves the generation of electrocoagulation fumes.

[0003] Smoke simulation can be broadly categorized into image-based methods and physics-based methods. In the early stages of smoke simulation, due to limitations in computer performance, image-based fractal textures and procedural textures were widely used. However, these methods suffered from inflexibility and easily revealed artificial traces in the simulated smoke movement. To address these shortcomings, researchers proposed physics-based discretized smoke simulation methods. By enhancing the depiction of the physical properties of the smoke, these methods more realistically simulate the entire smoke movement process. Physics-based discretized simulation methods can be broadly classified into two categories: mesh methods and particle methods. Mesh methods treat smoke as a discrete mesh structure. By solving partial differential equations and calculating the smoke's motion state using the mesh, parameters such as density and velocity on the mesh are updated to simulate the smoke's movement and diffusion process in real time. Particle methods treat the particles in the smoke as discrete particles, simulating the interactions and movements between these particles to construct the appearance of the smoke.

[0004] Currently, there are still some problems in smoke simulation, mainly 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 entire smoke simulation process has a large amount of computation and low real-time performance. In order to solve the complex computational 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 filament-mesh particle-particle (FMPP) method based on fast boundary wall turbulent smoke simulation 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 high-resolution smoke motion 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 issue of the realism of smoke simulation. While meeting the real-time requirements, the realism of the simulation still needs to be improved, including the study of physical properties such as object collisions and dispersion in the smoke process. To address this, Frank et al. proposed a two-way coupled simulation framework, using the SPH method to simulate diffusion regions such as sprays, achieving good realism in the simulated smoke.

[0005] Since none of the technologies mentioned so far can perfectly meet the requirements of smoke simulation in virtual surgical systems, it is necessary to study a suitable method that can achieve good real-time performance while ensuring realism. Summary of the Invention

[0006] The purpose of this invention is to provide a GPU-accelerated LBM smoke simulation method to improve the stability and realism of smoke effects.

[0007] This invention is achieved through the following technical solution.

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

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

[0010] Step 2: Before using GPU acceleration, make the following preparations: specify the calling interface for the compute shaders, define the data types required during the computation process, pre-allocate the storage space needed for the computation to store global constants and particle data, and define the kernel functions required for the computation. Once these preparations are complete, the Compute Shaders can be automatically imported into the Unity project, and then the computation can be performed using C# scripts.

[0011] Step 3: Create and initialize the particle entities for 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 computation 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 handling, flow handling, and boundary condition handling, 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-th direction, respectively, and Ω(f i ) is the collision operator, generally a single relaxation time collision operator is used, its formula is shown in the figure, τ is the single relaxation time, which satisfies v is the viscosity coefficient in fluid kinematics.

[0016] The formula for flow 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-th direction after the next time step. Here, Δt is the elapsed time between two time steps, s. i It is the direction vector to the adjacent unit.

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

[0020]

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

[0022] As attached Figure 2 The diagram shown illustrates a bending boundary condition. A boundary curve exists between the grid points, where P... f′ For the adjacent nodes of the fluid node, P f Point P is a boundary fluid point. s Point P is a solid point. w This is the physical boundary point between the solid and the fluid. a For fluid point P f Pointing to solid point P s speed, e a- This refers to the velocity in the opposite direction, i.e., the velocity 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 It is known, but due to the solid point P s Flowing towards fluid point P f The distribution function f i (p s (Unknown) Curve boundary conditions utilize the known distribution function f. i (p f To solve for the unknown distribution function f i (p s ).

[0023] Let the fluid point Pf to solid point P s If the distance is one unit (1), then μ = |p f -p w |(0<μ<1) represents the fluid point P f To the physical boundary P w The distance of fluid point P f to solid point P s The proportion of the total length of the connecting lines can be obtained by interpolation:

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

[0025] Simultaneously, the bounce boundary condition is applied to the wall node P. w :

[0026]

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

[0028]

[0029] After applying positional constraints, if the distance between two sliding constraint points is still less than the safety threshold L0, either point will be deleted to avoid redundant sliding constraint points interfering with the stitching interaction simulation process. Where u w For physical boundary point P w The speed at that location.

[0030] Step 5: Initiate the loop calculation process on the CPU, call the Compute Shader code, and utilize its provided functions to perform parallel calculations on the GPU. Then, through the predefined 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 complete, clean up the memory space allocated to the buffer on the CPU to avoid memory leaks. This process is done by using the cleanup interface provided by Compute Shader.

[0032] Compared with the prior art, the beneficial effects of the present invention are:

[0033] This invention provides a smoke simulation method for virtual surgery. Due to the inherent parallelism of the LBM method, it utilizes a Compute Shader for parallel computation and numerically solves the LBM equations by binding kernel functions, significantly reducing the computation time. Furthermore, comparative experimental results with existing methods demonstrate that the proposed smoke simulation method exhibits strong real-time performance and realism. Attached Figure Description

[0034] Figure 1 This is a schematic diagram of the D3Q15 scheme (D represents Dimension and Q represents Quantities of VelocityDirections).

[0035] Figure 2 This is a schematic diagram of bending boundary conditions.

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

[0037] Figure 4 This is a simulation of smoke effects in a virtual surgery. (a) Initiation of hemostasis procedure; (b) Generation of electrocoagulation smoke.

[0038] Figure 5 This is a flowchart of the present invention. Detailed Implementation

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

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

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

[0042] Step 2: Make preparations in advance.

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

[0044] 2) Define the specific data types required during the calculation process so that data calculation and storage can be conveniently performed during function calls.

[0045] 3) Define the amount of storage space required for the calculation, which will be used to store global constants and particle data during the calculation process.

[0046] 4) Define the kernel function required for the computation and assign it the corresponding thread group vector. This kernel function is the corresponding equation in the LBM method.

[0047] Step 3: Create and initialize the particle entities for the smoke effect on the CPU, reserve CPU memory for these particle systems, and allocate the necessary 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 each small grid having a side length of 10, and the initial vortex quantity of the vortex particles being 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 randomly launched at ±15° along the vertical direction at a rate of 500 vortex particles per second.

[0048] Step 4: Call the defined computation 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 handling, flow handling, and boundary condition handling, 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-th direction, respectively, and Ω(f i ) is the collision operator, generally a single relaxation time collision operator is used, its formula is shown in the figure, τ is the single relaxation time, which satisfies ν is the viscosity coefficient of the fluid kinematics.

[0052] The formula for flow 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-th direction after the next time step. Here, Δt is the elapsed time between two time steps, s. i It is the direction vector to the adjacent unit.

[0055] When using the D3Q15 scheme, expanding the macroscopic Navier-Stokes equations using the Chapman–Enskog expansion yields the corresponding momentum equilibrium distribution function of the LBM, i.e.:

[0056]

[0057]

[0058] Where w i In 1992, Qian et al. proposed the weighting factor of the D3Q15 scheme, as shown in the attached diagram. Figure 1 As shown, c s =δx / δt, where δx is the velocity of the lattice, typically taken as a value of δt. ρ is the density at the macroscopic level, and u is the macroscopic flow velocity.

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

[0060]

[0061] In handling boundary conditions, since standard bounce boundary conditions can only achieve first-order accuracy, for interactive smoke simulation, there are curved boundaries that cannot adapt to any unit point or midpoint, and half-bounce boundaries also cannot meet the requirements. Therefore, this paper adopts bounce boundary conditions. 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 a node.

[0062] As attached Figure 2 The diagram shown illustrates a bending boundary condition. A boundary curve exists between the grid points, where P... f′ For the adjacent nodes of the fluid node, P f Point P is a boundary fluid point. s Point P is a solid point. w This is the physical boundary point between the solid and the fluid. a For fluid point P f Pointing to solid point P s speed, e a- This refers to the velocity in the opposite direction, i.e., the velocity 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 It is known, but due to the solid point Ps Flowing towards fluid point P f The distribution function f i (p s (Unknown) Curve boundary conditions utilize the known distribution function f. i (p f To solve for the unknown distribution function f i (p s ).

[0063] Let the fluid point P f to solid point P s If the distance is one unit (1), then μ = |p f -p w |(0<μ<1) represents the fluid point P. f To the physical boundary P w The distance of fluid point P f to solid point P s The proportion of the total length of the connecting lines can be obtained by interpolation:

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

[0065] Simultaneously, the bounce boundary condition is applied to the wall node P. w :

[0066]

[0067] Then perform interpolation calculations on it:

[0068]

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

[0070]

[0071] After applying positional constraints, if the distance between two sliding constraint points is still less than the safety threshold L0, either point will be deleted to avoid redundant sliding constraint points interfering with the stitching interaction simulation process. Where u w For physical boundary point P w The speed at that location.

[0072] Step 5: Initiate the loop calculation process on the CPU, call the Compute Shader code, and utilize its provided functions to perform parallel calculations on the GPU. Then, through the predefined 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 complete, clean up the memory space allocated to the buffer on the CPU to avoid memory leaks. This process is done by using the cleanup interface provided by Compute Shader.

[0074] Table 1 shows the average time taken by the present invention and the traditional LBM algorithm without GPU acceleration under different smoke particle numbers.

[0075] Table 1

[0076]

[0077] As attached Figure 3 This is the effect of the smoke interaction simulation process, in which Figure 3 (a) shows the actual smoke conditions under natural conditions; Figure 3 (b) This illustrates the simulated smoke effect of the present invention. It can be seen that the present invention effectively simulates the details of real smoke, achieving a high degree of realism. (Appendix) Figure 4 For the smoke simulation effect in the virtual surgery system, among which Figure 4 (a) The initial state before starting hemostasis. Figure 4 (b) Effect diagram of surgical electrocoagulation smoke generation.

Claims

1. A GPU-accelerated LBM smoke simulation method, characterized by: Includes the following steps: Step 1: Construct organ soft tissue and surgical instrument models. First, collect and input diagnostic images, including CT and MRI scan images and related video data. Then, use the Blender 3D modeling tool to process them to generate OBJ format stereoscopic images of soft tissue models and surgical instruments. Step 2: Before using GPU acceleration, prepare the following: specify the calling interface for the compute shader, define the data types required during the computation process, pre-allocate the storage space needed for the computation to store global constants and particle data during the computation process, define the kernel function required for the computation, and specify the corresponding thread group vector for it. This kernel function corresponds to the various equations in the LBM method. After the above preparations are completed, ComputeShaders is automatically imported into the Unity project and then the computation is performed using C# scripts. Step 3: Create and initialize the particle entities for the smoke effect on the CPU, reserve CPU memory for these particle systems, and allocate the necessary 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 each small grid having a side length of 10. The initial vortex momentum of the vortex particles is 0.5 s. -1 The buoyancy of the particles is 0.35N; 3000 vortex particles are placed at the center of the coordinate axis and randomly launched along the vertical direction at ±15°, with a launch speed of 500 vortex particles per second. 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 handling, flow handling, and boundary condition handling, respectively. The formula for local collision is: in, and They are respectively The distribution function before and after the collision in the direction. It is a collision operator, and a single relaxation time collision operator is generally used; The formula for flow is: in, This represents the new distribution function in the i-th direction after the next time step; Δt is the elapsed time between two time steps, s. i It is the direction vector to the adjacent unit; Once the equilibrium distribution function is determined, a simple summation is used to calculate the macroscopic variables, where density... And the velocity u is: In handling boundary conditions, a bounce boundary condition is adopted: Let the fluid point P f to solid point P s If the distance is one unit (1), then , , representing fluid point P f To the physical boundary P w The distance of fluid point P f to solid point P s The proportion of the total length of the connecting lines can be obtained by interpolation: Simultaneously, the bounce boundary condition is applied to the wall node P. w : Finally, the boundary conditions can be expressed as: Among them, W i C represents the weighting factor for the D3Q15 scheme. s For lattice velocity, C s = δx / δt, where δx is the lattice step size and δt is the time step size; S i This is the direction vector to the adjacent cell; After applying positional constraints, if the distance between two sliding constraint points is still less than the safety threshold L0, either point will be deleted to prevent redundant sliding constraint points from interfering with the stitching interaction simulation process. Where u w For physical boundary point P w The speed at that location; Step 5: Start the loop calculation process on the CPU, call the Compute Shader code, and use its provided functions to perform parallel calculations on the GPU; then, through the interface of the predefined Compute Shader, update the calculation results on the CPU and send the smoke particle information processed by the GPU back to the CPU memory. Step 6: After all calculations are complete, clear the memory space allocated to the buffer on the CPU.

Images