GPU-Accelerated Adaptive Sampling Volume Rendering Method for Non-Uniform Cuboid Meshes

Abstract

The present invention discloses a GPU-accelerated adaptive sampling volume rendering method for non-uniform cuboid meshes. First, based on data bounding box rendering and GPU multi-threaded computing, a set of sampling rays for non-uniform cuboid meshes is generated. Subsequently, the intersection points of each ray with the mesh cell faces are calculated and the coordinate positions of the adaptive sampling points are determined. Then, the non-uniform cuboid mesh data is stored in the 3D texture memory of the GPU in a normalized manner, and texture hardware queries are performed using the weighted and corrected adaptive sampling point coordinates to obtain the trilinear interpolation result scalar values at the sampling points. Finally, the colors and transparencies corresponding to the sampling points are calculated according to the transfer function, and the Under operator in Alpha blending is used to achieve the correct color synthesis of all sampling points on the ray, and finally a visualization result image is formed. Experimental results show that the present invention can effectively improve the computational efficiency of Ray Casting volume rendering for non-uniform cuboid meshes and accelerate the visualization analysis and calculation process of non-uniform cuboid meshes.

Description

Technical Field

[0001] The present invention relates to the technical fields of mechanical finite element calculation and electromagnetic simulation. In particular, in the research on the visualization of three-dimensional scientific data volume rendering, a method for quickly obtaining the scalar field values of ray adaptive sampling points in a non-uniform cuboid grid by using GPU multi-threading and three-dimensional texture hardware interpolation technology is used to achieve the efficient adaptive sampling Ray Casting volume rendering of non-uniform cuboid grid data. Background Art

[0002] The Ray Casting volume rendering method is a three-dimensional visualization technology. This method emits rays directly from the viewpoint to the data volume, and performs data sampling on the ray path to obtain color accumulation to form a two-dimensional image visualization technology. The advantage of this method is that it can directly generate the result image from the original volume data without generating intermediate geometric results, and the visualization efficiency is higher. In addition, by adjusting the semi-transparency parameter, the Ray Casting method can directly reflect the internal structural characteristics of complex data, clearly show the structural details inside the data volume, and the performance is more intuitive. The defect of this method is that the calculation amount is large, and a large amount of computing power will be consumed for sampling and calculating the volume data. Therefore, in practical applications, acceleration hardware such as GPUs is often used to implement this algorithm to improve the calculation efficiency and achieve real-time interactive visualization.

[0003] Non-uniform cuboid grids are usually applied to numerical simulation fields such as three-dimensional finite element calculation and three-dimensional electromagnetic field calculation. This type of grid is composed of cuboid elements and is a three-dimensional extension of two-dimensional non-uniform rectangular grids. The characteristic of this type of grid is that its topological structure is the same as that of three-dimensional regular grids, both are structured right-angled hexahedron grids. The difference is that the distance between adjacent parallel planes of non-uniform cuboid grids changes non-uniformly, and the density of parallel planes can be adjusted according to the regional calculation accuracy. This flexibility makes non-uniform grids more effective than uniform grids when dealing with complex geometric shapes or physical phenomena with drastic changes. Therefore, in terms of the numerical simulation calculation efficiency, using non-uniform cuboid grids often has a higher calculation efficiency than uniform structured grids.

[0004] Although the non-uniform cuboid mesh is similar to the three-dimensional regular mesh in topology, there are significant differences in visualization. At present, there are generally two methods for Ray Casting volume rendering visualization of non-uniform cuboid meshes. The first method is the resampling method, that is, uniformly resampling the non-uniform cuboid mesh, interpolating the original non-uniform mesh to generate a uniform regular mesh, and then loading the uniform mesh into the GPU for visualization. The main problem with this method is that it is difficult to balance accuracy and efficiency. When the uniform resampling interval is relatively sparse, the original non-uniform cuboid data may be skipped, resulting in data feature loss; when the resampling interval is relatively dense, it will cause problems such as too many sampling points and too much calculation, and the visualization efficiency is not high; the second method is the adaptive sampling method. This method adopts an adaptive sampling strategy according to the different densities of the non-uniform grid area, performs encrypted sampling on the dense grid and sparse sampling on the sparse grid area, so as to effectively keep the data features from being lost. Although the second method can effectively capture data features, due to the uneven density of the grid, more query and judgment operations are required to realize the adaptive sampling process compared to the uniform grid. The adaptive calculation process is complex and cannot be unified and accelerated with the help of GPU hardware, making it difficult to give full play to the computational advantages of the simple topological structure of the rectangular grid. Summary of the invention

[0005] The present invention proposes a non-uniform cuboid grid efficient adaptive sampling Ray Casting volume rendering method implemented on a GPU. Firstly, based on data bounding box rendering and GPU multi-threaded calculation, a sampling ray set of the non-uniform cuboid grid is generated and the position of the adaptive sampling point on each ray is calculated. Then, the three-dimensional texture hardware device of the GPU is used to calculate the scalar value of the trilinear interpolation result at the sampling point position after weighted correction in a hardware query manner. Finally, the color and transparency are calculated according to the transfer function, and the Under operator in Alpha fusion is used to realize the correct color synthesis of all sampling points on the ray, so as to finally form a visual result image.

[0006] To achieve the above object, the present invention adopts the following technical solutions:

[0007] (1) Read the non-uniform rectangular grid data, obtain the relevant attribute parameters of the grid data body, establish the world coordinate system with the minimum coordinate point of the grid data as the origin of the world coordinate system, and determine the scene parameters;

[0008] (2) Taking the viewpoint as the starting point of the sampling ray, the parametric equation of the effective sampling ray is determined according to the external contour of the data body and the projection area of ​​the visible surface on the viewing plane;

[0009] (3) Enable multi-threading on the GPU. Each thread is responsible for processing one valid sampling ray, calculating the intersections of the valid sampling ray with all grid surfaces in the X, Y, and Z directions, and obtaining the set E of intersection sequence sets;

[0010] (4) Sort the intersections on each ray in the set E in ascending order of parameter t to obtain the ordered intersection sequence set E';

[0011] (5) Based on the ordered intersection set E', obtain the adaptive sampling point set S;

[0012] (6) Load the non-uniform cuboid grid data into the GPU three-dimensional texture memory, perform coordinate transformation on all sampling points in the set S, and use the transformed coordinates as input parameters for three-dimensional texture query to obtain the set V of trilinear interpolation results;

[0013] (7) Use the Under operator of the Alpha blending algorithm to blend the color values and transparency values sampled from near to far on each ray in the set V to obtain the final color of the pixel point corresponding to the sampling ray, and finally form the output result image.

[0014] In the step (1), the relevant attribute parameters of the non-uniform cuboid grid data include the length dimensions (Lx, Ly, Lz) in the X, Y, and Z directions, and the number of grid points (Nx, Ny, Nz) in the three directions. The relevant scene parameters include the viewpoint position coordinates, the viewing plane, the light source position, the lighting direction, and the result image resolution, etc.;

[0015] In the step (2), first, it is necessary to calculate the projection of the data volume on the viewing plane according to the data volume bounding box to determine the valid pixels on the viewing plane. The ray connecting the viewpoint and the valid pixels is the valid sampling ray; assuming the viewpoint position is P(x p , y p , z p ), then the valid ray passing through the valid pixel point Q(x q , y q , z q ) on the viewing plane (let Rn be the total number of valid pixels, 1 ≤ i ≤ Rn) The parametric equation expression is: (Let Rn be the total number of valid pixels, 1 ≤ i ≤ Rn) The parametric equation expression is:

[0016]

[0017] In the step (3), each GPU thread is responsible for the relevant calculations of one sampling ray The intersections of the ray with all grid surfaces in the X, Y, and Z directions form the set E i , then the set E of intersections formed by all rays intersecting with the data volume is E = {E i|E i is the set of intersection points of the i-th effective ray and all grid planes, where 1 ≤ i ≤ Rn. Taking the two-dimensional case as an example, the ray determined by the viewpoint P and the effective pixel Q passes through the data volume, and the intersection points with all the axis-parallel grid planes in the data volume are marked with the symbol "×", and all the intersection points E ij and the parameters t of the ray parameter equations corresponding to the intersection points are obtained ij to form the intersection point set E, where 1 ≤ i ≤ Rn and 1 ≤ j ≤ Sn i , and Sn i is the number of all sampling points on the ray .

[0018] In the step (4), the intersection points on each ray in the intersection point set E are sorted, that is, the elements in the intersection point set E corresponding to the ray i are sorted. Each intersection point E i in E ij corresponds to a t in the parameter equation ij value. The magnitude of the t ij value indicates the distance relationship between the intersection point and the viewpoint P. Therefore, sorting the elements in the E ij set according to the t i value can obtain the ordered intersection point sequence set E' on the ray i from near to far from the P point.

[0019] In the step (5), the sampling point set S = {S i | S i is all the adaptive sampling points on the ray }, and S i is calculated from E' i . Assuming that there are two adjacent points E on the ray ij (x ij , y ij , z ij ) and E i,j+1 (x i,j+1 , y i,j+1 , z i,j+1 ), then a sampling point S ij is taken between E i,j+1 and E ij . The coordinates of S ij can be expressed as where 1 ≤ j ≤ Sn i , and all the sampling points S on the ray ij constitute the set S i . Because of the adjacent grid plane intersection points E ij and Ei,j+1 Must be the entry point and exit point where the ray passes through a certain grid cell. Therefore, the above sampling point selection strategy can ensure that the ray There is a sampling point in all grid cells passed through by the ray, ensuring that the grid data features are not lost during the sampling process.

[0020] In step (6), the grid data needs to be loaded into the GPU three-dimensional texture memory first to achieve subsequent fast trilinear interpolation. The scalar field value at the sampling point S ij is obtained by trilinear interpolation calculation from the 8 cell vertices of the cell to which the sampling point belongs. Assume that the scalar values at the 8 grid points of the cell where the sampling point is located are v 000 , v 001 , … v 111 , then the scalar value v calculation formula at point S ij is:

[0021] v = v 000 (1 - x d )(1 - y d )(1 - z d ) + v 100 x d (1 - y d )(1 - z d ) + v 001 (1 - x d )(1 - y d )z d + v 101 x d (1 - y d )z d + v 010 x d (1 - y d )(1 - z d ) + v 110 x d y d (1 - z d ) + v 011 (1 - x d )y d z d + v 111 x d y d z d

[0022] where (x0, y0, z0) and (x1, y1, z1) are the grid points V 000 and V 111Coordinate values. According to the hardware characteristics of the GPU three-dimensional texture memory, it can be known that in the case of loading grid data into the three-dimensional texture, the above trilinear interpolation calculation process can be efficiently completed by the hardware. However, this hardware trilinear interpolation is only applicable to normalized regular grid cells. To achieve correct trilinear interpolation of non-uniform cuboid grids on the GPU three-dimensional texture memory, it is necessary to perform a mathematical transformation on the elements in the original sampling point set S. During the transformation process of converting the cuboid grid into a cube grid cell, to ensure that the scalar value at the sampling point remains unchanged after the grid coordinate transformation, it is necessary to ensure that the relative proportional position of the sampling point within the grid cell remains unchanged. Assume that the number of grid points in the three-dimensional directions are Nx, Ny, and Nz respectively, and the lengths of the grid in the three directions are Lx, Ly, and Lz respectively. The sampling point is located in the [k1, k2, k3]th grid cell (k1, k2, and k3 are three-dimensional array subscript indices). After the grid transformation, the corresponding sampling point is Then, before and after the grid coordinate transformation, the following equalities must hold:

[0023]

[0024] Since the three-dimensional texture coordinates in the GPU need to be represented in a normalized form, thus Lx = Ly = Lz = 1.0. Therefore, after substituting into the formula and simplifying, the coordinates of the sampling point S′ ij are: as follows:

[0025]

[0026] Using S′ ij as a parameter to call the three-dimensional texture sampling function in the GPU, the scalar value v′ at the sampling point S′ ij can be quickly obtained. From the previous discussion, it can be known that v = v′. Therefore, v′ can be directly used to query the transfer function to obtain the color value and transparency value at S ij . By executing the above trilinear interpolation steps in a multi-threaded parallel manner, the trilinear interpolation result set V corresponding to the sampling point set S can be obtained.

[0027] In the above step (7), since the sampling points on all rays are sorted from small to large according to the parameter t, the sampling points on the ray can directly use the Under operator in the Alpha blending algorithm for blending to obtain the final color value of the pixel corresponding to the ray. Since the color blending process of each ray is independent of each other, this step is also completed in parallel using GPU multi-threading, and finally the result image is output.

[0028] Compared with the prior art, the present invention has the following two aspects of technical innovations:

[0029] 1. The present invention proposes an adaptive sampling method for Ray Casting volume rendering on a GPU for non-uniform cuboid meshes. By utilizing the multi-thread technology of the GPU, it first calculates the intersection points of the effective projection rays and all mesh cell faces, and then takes the midpoint of adjacent intersection points as the adaptive sampling points for subsequent trilinear interpolation calculations. This adaptive sampling calculation method ensures that there are sampling points in each cell through which the ray passes, guarantees that data features are not lost due to too large a sampling step size, and effectively controls the sampling calculation amount, avoiding repeated sampling of meshes with a relatively large scale.

[0030] 2. The present invention proposes a series of coordinate transformation calculation methods, and realizes fast trilinear interpolation calculation for sampling points of non-uniform cuboid meshes on the GPU three-dimensional texture hardware, effectively improving the execution efficiency of the RayCasting algorithm for non-uniform cuboid mesh data. BRIEF DESCRIPTION OF THE DRAWINGS

[0031] Figure 1 FIG. is a schematic diagram of the principle for determining effective sampling rays by the projection of the data bounding box on the view plane in the present invention;

[0032] Figure 2 FIG. is a schematic diagram of the two-dimensional situation of sampling of non-uniform cuboid meshes by effective rays in the present invention; Ray The intersection situation with non-uniform mesh cell faces (the intersection points are marked with ×), and the adaptive sampling point S ij is located at the midpoint of the connection line of the mesh face intersection points E ij 、E i,j+1 .

[0033] Figure 3 FIG. is a schematic structural diagram of the two-dimensional array Array for storing effective ray intersection points in the present invention E .

[0034] Figure 4 FIG. is the transformation relationship between the cuboid mesh cell structure and the sampling point S ij in the present invention and the transformed cube mesh cell structure and the sampling point S' ij . The length, width and height of the original cuboid mesh cell are L1, L2 and L3 respectively, and the edge length is L after being transformed into a cube mesh cell.

[0035] Figure 5 FIG. is the Ray Casting volume rendering results of 4 test data in the verification experiment. The upper left, upper right, lower left and lower right figures are experimental data 1, data 2, data 3 and data 4 respectively. DETAILED DESCRIPTION OF THE INVENTION

[0036] The following explains the specific implementation manners of the present invention.

[0037] The Ray Casting volume rendering algorithm of the present invention is implemented using the CUDA programming model on the GPU. Assuming that the data type of the non-uniform cuboid grid data B to be processed currently is 32-bit floating-point scalar field data, the lengths of the grid data in the X, Y, and Z directions are (Lx, Ly, Lz) respectively, and the number of grid points in the three directions are (Nx, Ny, Nz) respectively. Then, when performing Ray Casting volume rendering on this data, the specific implementation method for each step corresponding to the content of the present invention is as follows:

[0038] (1) Read the non-uniform cuboid grid data B into the system memory and store it in a three-dimensional array of size [Nx, Ny, Nz]; apply for constant memory space on the GPU to store (Lx, Ly, Lz), (Nx, Ny, Nz) and the screen resolution R N values; take the minimum coordinate point of the data B as the origin, transform the viewpoint position and the viewing plane, and store them in the CUDA global variable memory space applied for. The viewpoint position coordinates P(x p , y p , z p )、the viewing plane normal and the viewing focus are stored in the form of a one-dimensional vector, and the viewport size is made to match the viewing plane resolution;

[0039] (2) According to the viewpoint and viewing plane positions, render the data bounding box to obtain the projection area of the data bounding box on the viewing plane. The connection lines between the viewpoint coordinates and the coordinates of all pixel points within this projection area are all valid sampling rays. The relationship among the viewpoint, the viewing plane, the data bounding box projection, and the valid sampling rays is as Figure 1 shown. Assuming that the coordinates of a valid pixel point within the projection area are Q(x q , y q , z q ), then the parametric equation expression of the valid ray is:

[0040]

[0041] Assume that the total number of pixels within the projection area is Rn, then there are Rn ray parametric equations corresponding to it. Enable a CUDA two-dimensional thread grid on the GPU, and the thread grid scale is that is, specify to start 128 threads within a CUDA thread block, and a total of thread blocks are started to calculate the parametric equations of Rn valid sampling rays, and store the expression form and coefficients of each ray parametric equation in the CUDA global variable array;

[0042] (3) Enable a CUDA two-dimensional thread grid on the GPU (the thread grid scale is ), each thread is assigned to calculate the intersection of a valid ray with all mesh unit faces in the X, Y, and Z directions. For example, we need The grid unit faces for intersection calculation are the axis-parallel faces passing through all grid points in data B. A total of Nx+Ny+Nz faces need to be intersected. Find the intersection point (two-dimensional example Figure 2 shown), then The set of intersections with all unit faces of B is E i , the set of intersection points of all valid rays and grid unit surfaces is E = {E i |E i is the set of intersections between the i-th valid ray and all mesh surfaces, 1≤i≤Rn}. During the implementation, a two-dimensional array Array of size [Rn, Nxyz] is dynamically requested in the GPU global storage. E It is used to store E (it is easy to prove that the number of intersection points between the ray and the grid unit surface does not exceed Nx+Ny+Nz, so the maximum value of the second dimension of the array is set to Nxyz=Nx+Ny+Nz). Each element of the array is the coordinate value of an intersection point and the parameter t of the ray parameter equation corresponding to the intersection point. When the ray passes through the data body, Array E The built-in parameter t for the intersection position of the subsequent vacancies in the array is 2 127 (maximum single-precision floating point value), indicating that the array element is not an intersection point. The two-dimensional array Array corresponding to the intersection set E E Structure Figure 3 As shown, Array E The i-th column set corresponds to the ray The set of all intersection points E i , where E ij is the coordinate of the jth intersection point of the i-th ray, t ij are the parameters of the ray equation at the intersection point (1≤i≤Rn, 1≤j≤Nxyz);

[0043] (4) Sort the intersection points on each ray in the set E from small to large by parameter t, that is, sort the array Array E Each column of is sorted from small to large according to the parameter t. From the meaning of the ray parameter equation, we can know that the ray intersection E ij The corresponding parameter t ij , represents E ij The distance from the ray origin (viewpoint), t ij The larger the value, the ij The farther away from the viewpoint, the intersection point set E corresponding to the i-th ray iSorting according to the parameter t will obtain an ordered intersection point sequence from near to far on the ray, which is convenient for subsequent Alpha blending operations. In the specific implementation process, taking Array E each column as the input, call the Sort function in the Thrust library of CUDA, specify the parameter t as the sorting key value, and implement sorting in ascending order of the t value to obtain the updated array Array E which is the column-ordered intersection point set E';

[0044] (5) Based on the column-ordered intersection point set E', find the adaptive sampling point set S, where S = {S i |S i is all the adaptive sampling points on the ray }. The calculation of S i is calculated from E′ i Assume that there are two adjacent points E on the ray ij (x ij , y ij , z ij ) and E i,j+1 (x i,j+1 , y i,j+1 , z i,j+1 ), then a sampling point S ij is taken between E i,j+1 and E ij with the coordinates where 1 ≤ j ≤ Sn i , and all the sampling points S on the ray ij constitute the set S i . Since the adjacent grid surface intersection points E ij and E i,j+1 must be the entry point and the exit point of the ray passing through a certain grid cell, the above sampling point calculation method can ensure that there is a sampling point in all the grid cells passed by the ray , ensuring that the grid data features are not lost during the sampling process. The selection of S ij and the positional relationship between E ij and E i,j+1 are as shown in Figure 2 .

[0045] When implementing the program, according to the algorithm requirements for calculating the adaptive sampling points S ij , apply for a CUDA two-dimensional global array Array S to store the sampling point coordinate set S, and enable CUDA multi-threading to update the value of the array Array S , and the update method is where 1 ≤ i ≤ Rn and 1 ≤ j ≤ Nxyz - 1. Note that during the calculation of sampling points, it is necessary to first determine whether the t value of the Array E element is 2 127 . If there is a case where the t value is 2 127 , it indicates that the current sampling area has exited the grid area, and the calculation of the current sampling point coordinate Array S [i, j] needs to be terminated;

[0046] (6) After calculating all the valid sampling point coordinate sets S, according to the steps of the Ray Casting volume rendering algorithm, start sampling and calculating the dataset B using the sampling points. The sampling process needs to be completed using trilinear interpolation. Assume that the scalar values at the 8 grid points of the cell where the sampling point S ij is located are v 000 , v 001 …v 111 (as shown in Figure 4 ), then the formula for the scalar value v at point S ij is:

[0047] v = v 000 (1 - x d )(1 - y d )(1 - z d ) + v 100 x d (1 - y d )(1 - z d ) + v 001 (1 - x d )(1 - y d )z d + v 101 x d (1 - y d )z d + v 010 x d (1 - y d )(1 - z d ) + v 110 x d y d (1 - z d ) + v 011 (1 - x d )y d z d + v 111 x d y d z d

[0048] where (x0, y0, z0) and (x1, y1, z1) are the grid points V 000and V 111 coordinate values of

[0049] According to the hardware characteristics of the GPU three-dimensional texture memory, in the case of loading grid data into the three-dimensional texture, the above trilinear interpolation calculation process can be efficiently completed by the hardware. However, this hardware trilinear interpolation is only applicable to normalized regular grid cells. To achieve correct trilinear interpolation of non-uniform cuboid grids on the GPU three-dimensional texture memory, mathematical transformation needs to be performed on the elements in the original sampling point set S. During the transformation process of the cuboid grid into cube grid cells (as shown in Figure 4 ), to ensure that the scalar value at the sampling point remains unchanged after the grid coordinate transformation, it is necessary to ensure that the relative proportional position of the sampling point within the grid cell remains unchanged. Assume that the sampling point is located in the [k1, k2, k3]th grid cell (k1, k2, and k3 are three-dimensional array subscript indices), and the corresponding sampling point after the grid transformation is Then, before and after the grid coordinate transformation, the following equation must hold:

[0050]

[0051] Since the three-dimensional texture coordinates in the GPU need to be represented in a normalized manner, thus Lx = Ly = Lz = 1.0. Substituting into the formula and simplifying, the coordinates of the sampling point S′ ij can be obtained as: as follows:

[0052]

[0053] Using the coordinates of S′ ij as parameters to call the three-dimensional texture sampling function in the GPU, the scalar value v′ at the sampling point S′ ij can be quickly obtained. From the previous discussion, it is known that v = v′. Therefore, v′ can be directly used to query the transfer function to obtain the color value and transparency value at S ij . By performing the above trilinear interpolation steps in parallel using a multi-threaded method, the trilinear interpolation result set V corresponding to the sampling point set S can be obtained.

[0054] When creating a three-dimensional texture on the GPU, a texture object handle of the cudaTextureObject_t type is used, and the CUDA API function cudaCreateTextureObject is called to create a three-dimensional texture object, and the three-dimensional data set B is uploaded to the three-dimensional texture object. The elements in the sampling point array Array S are transformed according to the calculation formula of S′ ij to obtain a new sampling point array Array S′ ; Array S′The sampled point values and the cudaTextureObject_t object in it are used as input parameters to call the CUDA API function tex3D <float>That is, return the color and transparency values at the sampling points, and then apply for an array Array S′ of the same size as the intersection point array Array C for storing the color and transparency array obtained by the query;

[0055] (7) As can be seen from the foregoing steps, each column element in Array C is an adaptive sampling point sequence from near to far on a certain ray. Therefore, by performing Alpha blending on Array C column by column, the final color and transparency values of the effective sampling rays can be obtained. Since the column elements in Array C are sorted from near to far from the viewpoint, the Under operator is used for Alpha blending. Assume that there are N i sampling points on the i-th ray. Array C [i, j].c and Array C [i, j].α are the color value and transparency value of the j-th sampling point on the i-th ray. Then, according to the calculation rule of the Under operator, with j as the loop variable, let j increase gradually from 1 to N i -1, calculate and update Array C as follows:

[0056] Array C [i, j].c = (1 - Array C [i, j - 1].α) · Array C [i, j].c + Array C [i, j - 1].c

[0057] Array C [i, j].α = (1 - Array C [i, j - 1].α) · Array C [i, j].α + Array C [i, j - 1].α

[0058] After updating Array C according to the above rules, blend Array C [i, N i -1].c and Array C [i, N i -1].α with the background color again, and the final result color and transparency value of the i-th pixel can be obtained. The color and transparency values of all pixels form the final result image for output display, and the entire RayCasting volume rendering calculation process ends.

[0059] Experimental verification and effects

[0060] To verify the effectiveness of the present invention, a typical GPU workstation was selected as the experimental hardware platform to verify the computational efficiency of visualizing using the adaptive sampling Ray Casting volume rendering method on non-uniform cuboid grid data.

[0061] The hardware configuration parameters of the GPU workstation selected for the experiment are as follows: The CPU is Xeon(R) E5-2650v2, with a main frequency of 2.60 GHz, a running memory capacity of 64 GB, the GPU model is NVIDIA GeForce RTX 4080, the independent video memory capacity is 16 GB, the number of CUDA cores is 9728, and the operating system is Windows 10 64-bit Professional Edition.

[0062] Four non-uniform cuboid scalar field grid data were selected for the experiment. Each grid point stores a single-precision floating-point scalar value. The grid dimensions and the uniformly resampled grid dimensions under the condition of maintaining the minimum grid features are shown in Table 1. Assuming a screen resolution of 1920×1440, the CPU serial adaptive sampling method, the GPU adaptive sampling method, and the GPU uniform sampling method were selected for Ray Casting computational performance comparison. The average time to draw a frame is shown in Table 1. The rendering results of the test data are as Figure 5 shown.

[0063] Table 1. Computational times of different Ray Casting volume rendering implementation methods

[0064]

[0065] From the statistical data in Table 1, it can be seen that for the 4 test data sets, the method proposed in the present invention shortens the rendering time by approximately 60.79%, 71.4%, 81.24%, and 88.69% respectively compared with the CPU serial adaptive method, and shortens the rendering time by approximately 9.93%, 23.25%, 55.15%, and 57.58% respectively compared with the GPU parallel uniform resampling method. The experimental results show that the technical method of the present invention can effectively improve the computational efficiency of the adaptive sampling Ray Casting volume rendering method.< / float>

Claims

1. A non-uniform cuboid grid adaptive sampling volume rendering method based on GPU acceleration, characterized in that: The following steps are involved: (1) Read the non-uniform rectangular grid data, obtain the relevant attribute parameters of the grid data body, establish the world coordinate system with the minimum coordinate point of the grid data as the origin of the world coordinate system, and determine the scene parameters; (2) Taking the viewpoint as the starting point of the sampling ray, the parametric equation of the effective sampling ray is determined according to the external contour of the data body and the projection area of ​​the visible surface on the viewing plane; (3) Enable multithreading on the GPU. Each thread is responsible for processing a valid sampling ray and calculating the intersection points between the valid sampling ray and all mesh surfaces in the X, Y, and Z directions to obtain the intersection point sequence set E. (4) Sort the intersection points on each ray in the set E from small to large according to the parameter t to obtain an ordered intersection point sequence set E'; (5) Based on the ordered intersection point set E', obtain the adaptive sampling point set S; (6) Load the non-uniform cuboid mesh data into the GPU 3D texture memory, perform coordinate transformation on all sampling points in the set S, use the transformed coordinates as input parameters for 3D texture query, and obtain a trilinear interpolation result set V; (7) Use the Under operator of the Alpha fusion algorithm to fuse the color values ​​and transparency values ​​sampled from near to far on each ray in the set V to obtain the final color of the pixel corresponding to the sampled ray, and finally form the output result image.

2. The GPU-accelerated non-uniform cuboid grid adaptive sampling volume rendering method according to claim 1, characterized in that: In the step (1), the relevant attribute parameters of the non-uniform rectangular grid data include the length dimensions (Lx, Ly, Lz) in the three directions of X, Y, and Z, and the number of grid points in the three directions (Nx, Ny, Nz), and the relevant scene parameters include the viewpoint position coordinates, the viewing plane, the light source position, the lighting direction, and the resulting image resolution.

3. The GPU-accelerated non-uniform cuboid grid adaptive sampling volume rendering method according to claim 1, characterized in that: In step (2), the projection of the data volume on the viewing plane needs to be calculated based on the bounding box of the data volume to determine the valid pixels that need to be calculated on the viewing plane. The ray connecting the viewpoint and the valid pixels is the valid sampling ray. Assume that the viewpoint position is P(x p ,y p ,z p ), then the effective pixel point Q(x q ,y q ,z q )'s effective ray The parametric equation expression is: Where Rn is the total number of valid pixels, 1≤i≤Rn.

4. The GPU-accelerated non-uniform cuboid grid adaptive sampling volume rendering method according to claim 1, characterized in that: In step (3), each GPU thread is responsible for a sampling ray Related calculations, rays The intersection points of all mesh surfaces in the X, Y, and Z directions form a set E i , then the intersection point set formed by all rays and the data volume is E = {E i |E i is the set of intersections between the i-th valid ray and all mesh surfaces, 1≤i≤Rn}.

5. The GPU-accelerated non-uniform cuboid grid adaptive sampling volume rendering method according to claim 1, characterized in that: In step (4), the intersection points on each ray in the intersection point set E are sorted, that is, the rays The corresponding intersection set E i Sort the elements in E i Each intersection point E ij Each corresponds to a t in the parametric equation ij Value, t ij The value indicates the distance between the intersection point and the viewpoint P. ij Value pair E i Sort the elements in the collection and get the ray The ordered intersection sequence set E from point P from near to far i ′.

6. The GPU-accelerated non-uniform cuboid grid adaptive sampling volume rendering method according to claim 1, characterized in that: In step (5), the sampling point set S = {S i |S i For rays All adaptive sampling points on}, S i By E i 'Calculated; Assume that the ray There are two adjacent points E on ij (x ij ,y ij ,z ij ) and E i,j+1 (x i,j+1 ,y i,j+1 ,z i,j+1 ), then E ij With E i,j+1 Take a sampling point S between ij , the coordinates are Where 1≤j≤Sn i ,ray All sampling points S on ij The composition set S i .

7. The GPU-accelerated non-uniform cuboid grid adaptive sampling volume rendering method according to claim 1, characterized in that: In the step (6), the grid data needs to be loaded into the GPU 3D texture memory first to implement the subsequent fast trilinear interpolation; the sampling point S ij The scalar field value at is calculated by trilinear interpolation of the eight unit vertices of the unit to which the sampling point belongs; Assume that the scalar values ​​at the eight grid points of the unit where the sampling point is located are v 000 ,v 001 ,…v 111 , then point S ij The scalar value v is calculated as: v=v 000 (1-x d )(1-y d )(1-z d )+v 100 x d (1-y d )(1-z d )+v 001 (1-x d )(1-y d )z d +v 101 x d (1-y d )z d +v 010 x d (1-y d )(1-z d )+v 110 x d y d (1-z d ) +v 011 (1-x d )y d z d +v 111 x d y d z d in (x0, y0, z0) and (x1, y1, z1) are the grid points V 000 and V 111 The coordinate value of .

8. The GPU-accelerated non-uniform cuboid grid adaptive sampling volume rendering method according to claim 1, characterized in that: In step (6), it is assumed that the number of grid points in the three dimensions is Nx, Ny and Nz respectively, and the length of the grid in the three directions is Lx, Ly and Lz respectively. Located in the [k1, k2, k3]th grid cell (k1, k2 and k3 are the subscript indices of the three-dimensional array), the corresponding sampling point after the grid transformation is Then the following equations must hold before and after the grid coordinate transformation: Since the three-dimensional texture coordinates in the GPU need to be normalized, Lx = Ly = Lz = 1.0, after simplifying the formula, we can get the sampling point S ij ′’s coordinates for: S ij ′ as a parameter to call the 3D texture sampling function in the GPU, and the sampling point S can be quickly obtained. ij ′; use v′ to query the transfer function and get S ij The color value and transparency value at the sampling point are obtained by performing the above trilinear interpolation steps in parallel in a multi-threaded manner to obtain a trilinear interpolation result set V corresponding to the sampling point set S.

9. The GPU-accelerated non-uniform cuboid grid adaptive sampling volume rendering method according to claim 1, characterized in that: In the step (7), the sampling points on all rays are sorted from small to large according to the parameter t, and the sampling points on the rays are directly fused using the Under operator in the Alpha fusion algorithm to obtain the final color value of the pixel corresponding to the ray. The color fusion process of each ray is completed in parallel using GPU multi-threading, and the resulting image is finally output.

Images