# A Novel Ray-Casting Algorithm Using Dynamic Adaptive Sampling

Authors: Huadeng Wang; Guang Xu; Xipeng Pan; Zhenbing Liu; Rushi Lan; Xiaonan Luo
Journal: Wireless Communications and Mobile Computing
 (2020)

---

## Abstract

Ray-casting algorithm is an important volume rendering algorithm, which is widely used in medical image processing. Aiming to address the shortcomings of the current ray-casting algorithms in 3D reconstruction of medical images, such as slow rendering speed and low sampling efficiency, an improved algorithm based on dynamic adaptive sampling is proposed. By using the central difference gradient method, the corresponding sampling interval is obtained dynamically according to the different sampling points. Meanwhile, a new rendering operator is proposed based on the color value and opacity changes before and after the ray enters the volume element, and the resistance luminosity. Compared with the state of other algorithms, experimental results show that the method proposed in this paper has a faster rendering speed while ensuring the quality of the generated image.


---

## Body

## 1. Introduction

The ray-casting algorithm was proposed by Levoy in 1988. As an important volume rendering algorithm, the ray-casting algorithm is widely used in medical images, fluid and quantum mechanics, geographical exploration, and finite element analysis [1, 2]. In the basic principle of the ray-casting algorithm, each pixel in the imaging plane sends out a ray along the line of sight through the volume data, resampling along the ray at the same distance, obtaining the color value and opacity of each resampling point, and then synthesizing the color and opacity of each resampling point on the ray from the front to the back or from the back to the front.Although the traditional ray projection algorithm has a high quality of synthetic image, it has two disadvantages. One is that there are a lot of volume data and a lot of addition and multiplication operations in the interpolation process, which makes the whole rendering process very slow. Second, after the sampling point is determined, the use of cubic linear interpolation will affect the real-time rendering process.Some researchers use the closest interpolation instead of trilinear interpolation to remove some unnecessary sampling points. Some researchers achieve this goal by reducing the intersection of the light field and the data field. For example, Siddon [3] replaced the intersection of ray and voxel with the intersection of ray and plane. Ross et al. [4] employed the Bresenham algorithm to determine the intersection voxels. The data domain octree organization method proposed by Levoy [5] could easily skip empty voxels and reduce the time complexity of the intersection to O(nlog n). Jian [6] and others chose the projection light according to the correlation of the critical voxels, but the sharpness of the image would be worse.There are two ways to improve the ray projection algorithm. One is to simplify the operation of the sampling process, such as using trilinear interpolation instead of complex cubic convolution interpolation [7–10]. The other is to reduce the drawing time by reducing the number of sampling points on each ray. For the existing empty sampling points which are not helpful for image rendering, most of the methods are to eliminate them directly without considering the correlation between them and the surrounding voxels. When users interact with each other, the image will be blurred and the user experience will be affected.Considering that most of the adaptive sampling methods used in most papers are equal-interval sampling, the threshold value of the isosurface is set first, and the isodistance sampling is started from the viewpoint direction along the line of sight direction until the sampling value exceeds the set threshold value. Finally, the linear interpolation is done between the sampling pointsPn and Pn+1 to get the intersection coordinates of the line of sight and the isosurface. The disadvantage of this method is that there are too many sampling points. Therefore, a large number of sampling points are empty voxels, which wastes precious time and is not conducive to user interaction. Another method is based on a given sampling frequency, but the quality of voxel rendering is very poor. Some scholars try to improve the sampling frequency, adjust the sampling frequency by the distance proportion between the tangent point of the surface normal and the viewpoint, and simplify the image synthesis operator. Another is to resample at a lower sampling frequency first. If the data values of two adjacent resample points differ greatly, then increase the sampling frequency in such a high-frequency area. In the place where the voxel value changes slowly and quickly, the sampling interval should be small to avoid missing voxel information. If the voxel value changes slowly, increase the sampling interval appropriately, so that a large number of empty voxels can be skipped. However, if the sampling interval is too large, the rendering quality will be low, resulting in a large “void” effect [11, 12], that is, the phenomenon of the artifact.There are two shortcomings in the current ray casting algorithm for medical image rendering, one is the slow rendering speed, the other is the low efficiency and high time complexity of the whole sampling process. In order to solve the above shortcomings, this paper proposes an algorithm. The algorithm uses the central differential gradient method. By this method, the specific sampling distance can be determined dynamically according to the different distances of sampling points, which can improve the sampling efficiency. On the other hand, according to the change of color value and opacity before and after the ray enters the voxel, we propose a new rendering operator to improve the final image rendering speed.This paper solves the problem of how to find the most suitable sampling interval dynamically. Through the full increment formula of the voxel value difference between two adjacent sampling points, combined with the central difference method, the calculated gradient, and the preset voxel value threshold, the sampling point interval can be calculated dynamically according to the different positions of sampling points. Thus, we can update the sampling interval dynamically on the premise of ensuring the sampling density. Aiming at the problem that the search precision of the intersection point between the light and the actual isosurface will be reduced when the sampling point interval is large, this paper proposes a recursive estimation method based on the intersection point to estimate the intersection point between the sampling point and the isosurface, which is simple and easy to operate. The algorithm we proposed is a novel ray-casting algorithm using dynamic adaptive sampling, which is called the ray-casting DAS algorithm for short. After determining the sampling point, aiming at the general paper’s shortcomings of using the nearest interpolation or trilinear interpolation to synthesize the image according to the opacity and color value of eight data points equidistant from the sampling point, this paper adopts a method of light blocking, which is to change the opacity and color values as well as the photometry before and after the ray enters the volume element and improve the composition operator. The basic algorithm structure and process are listed in Figure1. Experimental results show that the improved method can effectively improve the speed of image rendering and avoid the generation of artifacts. The main innovation of this paper is dynamic adaptive sampling and improving the rendering operator to improve the efficiency and rendering time of the algorithm. Another unique innovation of this paper is the introduction of the distance function to improve the image noise and rendering quality.Figure 1
The brief structure of the ray-casting DAS algorithm.The remainder of this paper is as follows. The second part summarizes the latest research results of some ray projection algorithms and compares their advantages and disadvantages. In the third part, the ray projection algorithm based on adaptive sampling is introduced. The fourth part is the comparison between the ray-casting DAS and other ray projection algorithms in several open-source datasets. The fifth part is the summary of ray-casting DAS and the assumption of improving in the future.


## 2. Related Work

Wu et al. [13] proposed a GPU ray projection method for visualizing 3D pipelines in a virtual sphere. In their method, the pipeline data was initially divided into several data blocks. The pipeline centerline in each block was then segmented and coded, and a thicker pipeline envelope was then constructed using geometric shaders. Finally, elaborate 3D pipes could be rendered using pixel shaders. Compared with the traditional polygon-based pipeline data partition method, this method improves the rendering frame rate under the same pixel-level precision display effect. The visualization of 3D pipe thickness is realized without affecting the rendering efficiency.Luo et al. [14] proposed a rendering algorithm based on the height field [15] and an improved algorithm of finding the intersection of isosurfaces to solve the problem of the slow speed of finding the intersection of light and volume data. However, because the volume data will be imaged many times due to the reflection of the surface of the object, so it needs to do many volume data rendering operations, which hurt the performance of the algorithm.Francisco Sans et al. [16] proposed to combine CPU based ray projection with GPU, and reduce the time complexity of rendering by dynamic gradient interpolation of voxels. It realized single channel GPU rendering and meet the real-time interaction requirements of users.Binyahib et al. [17] proposed a parallel volume rendering algorithm which combines the classic object sequence and image sequence techniques. The algorithm runs on an unstructured grid (and structured grid), so it can deal with block boundary crossing in a complex way. It can also effectively deal with situations that are prone to load imbalance [18]. At the largest scale, it can be expanded to 8192 processors and run on datasets of more than 1 billion cells.Tao et al. [19] proposed a ray projecting algorithm combining the fitness estimation based on ray-casting with the global optimization method of improved adaptive differential evolution. This method eliminated the fine registration step of the famous iterative nearest point (ICP) algorithm [15, 20, 21], and it is the first direct global registration algorithm. Moreover, the calculation of the fitness of the global optimization method was accelerated, and the search space was effectively used to find the optimal transformation solution. The algorithm was successfully implemented in parallel mode on multicore computer processors.


## 3. Our Proposed Approach

### 3.1. Data Preprocessing Based on Fuzzy Enhancement

Because in the ray projection algorithm, in the mapping from three-dimensional to two-dimensional, the image edge, region, texture, and so on are inevitably blurred due to the loss of information. Therefore, the algorithm proposed in this paper first preprocesses the data, uses the method of enhancement operator [22] to calculate the blur rate of the enhanced image, and then selects the best threshold.The purpose of introducing a fuzzy enhancement operator [23] is to enhance the contrast of the image region by different enhancement processing, that is, the black target area is attenuated and the white background area is enhanced. The enhanced image not only retains the information of the original image but also makes the layers of each region clearer, which further reduces the image ambiguity. In this way, the change of the processed data is beneficial to the extraction and selection of sampling points for the subsequent volume rendering algorithm.


### 3.2. Intersection Processes of Light and Plane

First, the parametric equation of light emitted by a pixel in the plane can be expressed as:(1)x=x0+lt,y=y0+mt,z=z0+nt,where l,m,n represents the direction vector of the line in the object space, x0,y0,z0represents a point on the line. x, y, and z denote the three-dimensional coordinates of any point on the line, and t is the parameter of the linear parametric equation, indicating the number of directions of the line.Let the size of the regular data field beL×M×N and logical coordinates x,y,z=i×X,j×Y,k×Z, in which, i=0,1⋯,L−1, j=0,1,⋯,M−1, and k=0,1,⋯,N−1. Here L, M, and N denote the number of X, Y, and Z plane families, respectively, and x,y,z represents the projection coordinates of a specific point in the X,Y,Z plane of the specific plane family. The data distribution of the volume data field is represented by the plane family. The definition of the plane family is as follows:
(2)x=i×X,i=0,1,⋯,L−1,y=j×Y,j=0,1,⋯,M−1,z=k×Z,k=0,1,⋯,N−1.Because the normal vector of lightl,m,n≠0, let us set n≠0. Then, the parameters of the intersection of light and plane z=k×Z,k=0,1,⋯,N−1 can be calculated as follows:
(3)tk=k⋅Z−z0n,k=0,1,⋯,N−1.Further, the intersection of light andz=k×Z can be expressed as:
(4)x0+−z0+kZ⋅ln,y0+−z0+kZ⋅mn,kZ,k=0,1,⋯,N−1.Eq. (4) infers that light and two adjacent planes z=k×Z and z=k+1×Z intersection parameters tk and tk+1 satisfying Eq. (5)
(5)tk+1=tk+Zn.The intersection points of the straight line andz plane family can be obtained by using the above iteration formula:
(6)xk+1=xk+Zn,yk+1=yk+mnZ,k=0,1,⋯,N−2,zk+1=zk+Zn.Similarly, the recurrence formula to the intersection of the line and thex,y plane can be obtained.


### 3.3. Dynamic Adaptive Sampling

The intersection points of the line and plane family obtained by the above operations are not the sampling points we need, but we can get the sampling points on the line based on these intersections. The general method is to obtain the sampling points by equidistant sampling with a fixed distanced. Here, we consider how to dynamically adjust the adaptive sampling interval.Letg be a three-dimensional volume data field. According to the total differential formula of the binary function, the difference between the voxel values of two adjacent sampling points can be obtained as follows:
(7)Δg=gn+1−gn=gx′·△x+gy′·△y+gz′·△z+oρ≈∇g·L·dn,where ▽g is the gradient of the 3D digital theater and dn is the coordinate distance of the sampling point. L is the unit vector of the light direction.Becausegn has been estimated, and gn+1 is unknown. It may be assumed that the voxel value gn+1 of the sampling point Pn+1 is T (T is the set voxel value threshold). According to Eq. (7), the sampling interval can be approximately expressed as:
(8)dn≈T−gn∇g·L.It can be seen from Eq. (8) that if the projection value of the gradient at a point in the direction of light is smaller, and the difference between gn and T is large, then the interval dn will be larger; otherwise, dn will be too small. In order to prevent the sampling interval from being too large or too small, an upper and lower bound can be added to the sampling interval of Eq. (8). It is as follows:
(9)dn=D10≤T−gn∇g·L≤D1,D2T−gn∇g·L≥D2,orT−gn∇g·L<0,T−gn∇g·LD1<T−gn∇g·L<D2.D
1 and D2 are upper and lower bounds of sampling, respectively.In order to simplify the calculation, the center difference method is used to calculate the gradient, and the calculation method is as follows:(10)∇g=12gx+1,y,z−gx−1,y,z,gx,y+1,z−gx,y−1,z,gx,y,z+1−gx,y,z−1.Because the gradient calculation will affect the rendering speed if it is placed in the sampling stage, so we put the gradient calculation in the data preprocessing stage and use a linear table to store the calculated gradient value. When calculating the sampling step length, we look up the gradient value closest to the sampling point in the linear table as the gradient of the sampling point. Assuming the number of sampling points isn, then the time complexity of this step is On.After determining the sampling interval, we need to find the intersection of two adjacent sampling points and the isosurface and then synthesize the image through the color value and opacity of the intersection. In order to improve the imaging quality, this paper presents a recursive method for intersection estimation.Supposing that we need to find the intersection point of the two sampling pointsPn and Pn+1 and the isosurface. First, we take the midpoint coordinate P0 of Pn and Pn+1. If the corresponding voxel value at P0 is greater than the set threshold T, then P0 is taken as the new Pn+1. If the corresponding voxel value at P0 is less than the set threshold T, then P0 is taken as the new Pn, and then continue to take the midpoint of Pn and Pn+1 until the distance between Pn and Pn+1 is very small (less than the given value). Then, the final Pn and Pn+1 are linearly interpolated to obtain the coordinates of P0, which is the coordinates of the intersection.The process of finding the intersection coordinates is shown in Algorithm1. Through the algorithm, we can calculate the intersection point of sampling points Pn and Pn+1 fast. Moreover, we set the threshold of the number of iterations, which improves the convergence speed of the algorithm and indirectly improves the rendering speed of the image.Algorithm 1. solution process of intersection coordinates
Input:Pn, Pn+1 number_Current, Total.
Output:P0
1: Take the coordinateP0 of the geometric midpoint of Pn and Pn+1 to obtain the sampling value g.
2: T is set threshold, and the number current represent the current iterations, which is required to less than the Total.
3: Ifg>T then
4:Pn+1=P0
5: Else
6:Pn=P0
7: WhilePn−Pn+1<ε do
8:   By linear interpolation ofPn and Pn+1, the coordinates of the intersection point of Pn and Pn+1
9: If number_Current>Total then return P0
10: IfP0 is close toPn then
11:Pn=P0
12: Else
13:Pn+1=P0
14: ReturnP0


### 3.4. Improving Rendering Operator

After selecting the sampling point and intersection point in the above steps, the traditional method is to do linear interpolation around the eight closest data points of the intersection point to synthesize the image [24]. For example, the ray projection algorithm based on trilinear interpolation for fiber surface is adopted by Francisco Sans, and the time complexity of the trilinear interpolation is On3.The formula of cubic linear interpolation is as follows:(11)Si,j,k=P000i−1j−1k−1+P100ij−1k−1+P010i−1jk−1+P110ijk−1+P001i−1j−1k+P101ij−1k+P011i−1jk+P111ijk.Among them,Si,j,k represents the voxel value of the corresponding spatial coordinate point Pi,j,k,P000,P001,⋯,P111 represents the voxel value of the eight adjacent data points of Pi,j,k, respectively.In this paper, we consider a method of ray blocking and synthesize a new rendering operator through the change of color and opacity before and after the ray enters and leaves the voxel, as well as the corresponding opacity. The experimental results show that it can effectively reduce the time complexity and improve the rendering speed without affecting the imaging quality.Let the opacity and color values of theith volume element be Pi and Qi, respectively, the opacity and color values before the ray enters the ith volume element P’i and Q’i, respectively, and the opacity and color values after the ray enters are P”i and Q”I, respectively, and let the opacity of the cell be Gi; the composite operator of the point image is
(12)P′′i×Q′′i=Pi×Qi×1−Gi+P′i×Q′i.The relationship between the opacity and opacity of current voxels is(13)P′′i=Pi×1−Gi+P′i.At the same time, we notice that not all the light can finally reach the imaging screen for imaging. The initial point has the smallest opacity. With the increase of data points through which the ray passes, the opacity also increases. Finally, after the ray enters the nth data point, the opacity becomes the maximum 1, so the image synthesis operator can be simplified as:(14)Q′′i=∑i=1nPi×Qi×1−Gi.From the above formula, it can be found that wheng=1, that is the maximum photometry, the right side of the equation is equal to 0, which means that the ray that does not reach the screen does not participate in the calculation of the synthesis operator. At the same time, the time complexity of single Pi and Qi calculation is On, so the time complexity of the final synthesis operator is On2<On3, so the algorithm effectively improves the rendering efficiency. On the other hand, it can not only improve the time of image rendering but also effectively prevent the generation of artifacts and improve the quality of rendering.


## 3.1. Data Preprocessing Based on Fuzzy Enhancement

Because in the ray projection algorithm, in the mapping from three-dimensional to two-dimensional, the image edge, region, texture, and so on are inevitably blurred due to the loss of information. Therefore, the algorithm proposed in this paper first preprocesses the data, uses the method of enhancement operator [22] to calculate the blur rate of the enhanced image, and then selects the best threshold.The purpose of introducing a fuzzy enhancement operator [23] is to enhance the contrast of the image region by different enhancement processing, that is, the black target area is attenuated and the white background area is enhanced. The enhanced image not only retains the information of the original image but also makes the layers of each region clearer, which further reduces the image ambiguity. In this way, the change of the processed data is beneficial to the extraction and selection of sampling points for the subsequent volume rendering algorithm.


## 3.2. Intersection Processes of Light and Plane

First, the parametric equation of light emitted by a pixel in the plane can be expressed as:(1)x=x0+lt,y=y0+mt,z=z0+nt,where l,m,n represents the direction vector of the line in the object space, x0,y0,z0represents a point on the line. x, y, and z denote the three-dimensional coordinates of any point on the line, and t is the parameter of the linear parametric equation, indicating the number of directions of the line.Let the size of the regular data field beL×M×N and logical coordinates x,y,z=i×X,j×Y,k×Z, in which, i=0,1⋯,L−1, j=0,1,⋯,M−1, and k=0,1,⋯,N−1. Here L, M, and N denote the number of X, Y, and Z plane families, respectively, and x,y,z represents the projection coordinates of a specific point in the X,Y,Z plane of the specific plane family. The data distribution of the volume data field is represented by the plane family. The definition of the plane family is as follows:
(2)x=i×X,i=0,1,⋯,L−1,y=j×Y,j=0,1,⋯,M−1,z=k×Z,k=0,1,⋯,N−1.Because the normal vector of lightl,m,n≠0, let us set n≠0. Then, the parameters of the intersection of light and plane z=k×Z,k=0,1,⋯,N−1 can be calculated as follows:
(3)tk=k⋅Z−z0n,k=0,1,⋯,N−1.Further, the intersection of light andz=k×Z can be expressed as:
(4)x0+−z0+kZ⋅ln,y0+−z0+kZ⋅mn,kZ,k=0,1,⋯,N−1.Eq. (4) infers that light and two adjacent planes z=k×Z and z=k+1×Z intersection parameters tk and tk+1 satisfying Eq. (5)
(5)tk+1=tk+Zn.The intersection points of the straight line andz plane family can be obtained by using the above iteration formula:
(6)xk+1=xk+Zn,yk+1=yk+mnZ,k=0,1,⋯,N−2,zk+1=zk+Zn.Similarly, the recurrence formula to the intersection of the line and thex,y plane can be obtained.


## 3.3. Dynamic Adaptive Sampling

The intersection points of the line and plane family obtained by the above operations are not the sampling points we need, but we can get the sampling points on the line based on these intersections. The general method is to obtain the sampling points by equidistant sampling with a fixed distanced. Here, we consider how to dynamically adjust the adaptive sampling interval.Letg be a three-dimensional volume data field. According to the total differential formula of the binary function, the difference between the voxel values of two adjacent sampling points can be obtained as follows:
(7)Δg=gn+1−gn=gx′·△x+gy′·△y+gz′·△z+oρ≈∇g·L·dn,where ▽g is the gradient of the 3D digital theater and dn is the coordinate distance of the sampling point. L is the unit vector of the light direction.Becausegn has been estimated, and gn+1 is unknown. It may be assumed that the voxel value gn+1 of the sampling point Pn+1 is T (T is the set voxel value threshold). According to Eq. (7), the sampling interval can be approximately expressed as:
(8)dn≈T−gn∇g·L.It can be seen from Eq. (8) that if the projection value of the gradient at a point in the direction of light is smaller, and the difference between gn and T is large, then the interval dn will be larger; otherwise, dn will be too small. In order to prevent the sampling interval from being too large or too small, an upper and lower bound can be added to the sampling interval of Eq. (8). It is as follows:
(9)dn=D10≤T−gn∇g·L≤D1,D2T−gn∇g·L≥D2,orT−gn∇g·L<0,T−gn∇g·LD1<T−gn∇g·L<D2.D
1 and D2 are upper and lower bounds of sampling, respectively.In order to simplify the calculation, the center difference method is used to calculate the gradient, and the calculation method is as follows:(10)∇g=12gx+1,y,z−gx−1,y,z,gx,y+1,z−gx,y−1,z,gx,y,z+1−gx,y,z−1.Because the gradient calculation will affect the rendering speed if it is placed in the sampling stage, so we put the gradient calculation in the data preprocessing stage and use a linear table to store the calculated gradient value. When calculating the sampling step length, we look up the gradient value closest to the sampling point in the linear table as the gradient of the sampling point. Assuming the number of sampling points isn, then the time complexity of this step is On.After determining the sampling interval, we need to find the intersection of two adjacent sampling points and the isosurface and then synthesize the image through the color value and opacity of the intersection. In order to improve the imaging quality, this paper presents a recursive method for intersection estimation.Supposing that we need to find the intersection point of the two sampling pointsPn and Pn+1 and the isosurface. First, we take the midpoint coordinate P0 of Pn and Pn+1. If the corresponding voxel value at P0 is greater than the set threshold T, then P0 is taken as the new Pn+1. If the corresponding voxel value at P0 is less than the set threshold T, then P0 is taken as the new Pn, and then continue to take the midpoint of Pn and Pn+1 until the distance between Pn and Pn+1 is very small (less than the given value). Then, the final Pn and Pn+1 are linearly interpolated to obtain the coordinates of P0, which is the coordinates of the intersection.The process of finding the intersection coordinates is shown in Algorithm1. Through the algorithm, we can calculate the intersection point of sampling points Pn and Pn+1 fast. Moreover, we set the threshold of the number of iterations, which improves the convergence speed of the algorithm and indirectly improves the rendering speed of the image.Algorithm 1. solution process of intersection coordinates
Input:Pn, Pn+1 number_Current, Total.
Output:P0
1: Take the coordinateP0 of the geometric midpoint of Pn and Pn+1 to obtain the sampling value g.
2: T is set threshold, and the number current represent the current iterations, which is required to less than the Total.
3: Ifg>T then
4:Pn+1=P0
5: Else
6:Pn=P0
7: WhilePn−Pn+1<ε do
8:   By linear interpolation ofPn and Pn+1, the coordinates of the intersection point of Pn and Pn+1
9: If number_Current>Total then return P0
10: IfP0 is close toPn then
11:Pn=P0
12: Else
13:Pn+1=P0
14: ReturnP0


## 3.4. Improving Rendering Operator

After selecting the sampling point and intersection point in the above steps, the traditional method is to do linear interpolation around the eight closest data points of the intersection point to synthesize the image [24]. For example, the ray projection algorithm based on trilinear interpolation for fiber surface is adopted by Francisco Sans, and the time complexity of the trilinear interpolation is On3.The formula of cubic linear interpolation is as follows:(11)Si,j,k=P000i−1j−1k−1+P100ij−1k−1+P010i−1jk−1+P110ijk−1+P001i−1j−1k+P101ij−1k+P011i−1jk+P111ijk.Among them,Si,j,k represents the voxel value of the corresponding spatial coordinate point Pi,j,k,P000,P001,⋯,P111 represents the voxel value of the eight adjacent data points of Pi,j,k, respectively.In this paper, we consider a method of ray blocking and synthesize a new rendering operator through the change of color and opacity before and after the ray enters and leaves the voxel, as well as the corresponding opacity. The experimental results show that it can effectively reduce the time complexity and improve the rendering speed without affecting the imaging quality.Let the opacity and color values of theith volume element be Pi and Qi, respectively, the opacity and color values before the ray enters the ith volume element P’i and Q’i, respectively, and the opacity and color values after the ray enters are P”i and Q”I, respectively, and let the opacity of the cell be Gi; the composite operator of the point image is
(12)P′′i×Q′′i=Pi×Qi×1−Gi+P′i×Q′i.The relationship between the opacity and opacity of current voxels is(13)P′′i=Pi×1−Gi+P′i.At the same time, we notice that not all the light can finally reach the imaging screen for imaging. The initial point has the smallest opacity. With the increase of data points through which the ray passes, the opacity also increases. Finally, after the ray enters the nth data point, the opacity becomes the maximum 1, so the image synthesis operator can be simplified as:(14)Q′′i=∑i=1nPi×Qi×1−Gi.From the above formula, it can be found that wheng=1, that is the maximum photometry, the right side of the equation is equal to 0, which means that the ray that does not reach the screen does not participate in the calculation of the synthesis operator. At the same time, the time complexity of single Pi and Qi calculation is On, so the time complexity of the final synthesis operator is On2<On3, so the algorithm effectively improves the rendering efficiency. On the other hand, it can not only improve the time of image rendering but also effectively prevent the generation of artifacts and improve the quality of rendering.


## 4. Experimental Results

In this section, several experiments and comparison results are reported to evaluate the performance of the proposed algorithm ray-casting DAS. We compare the ray-casting DAS with other algorithms in several aspects, such as time complexity, image rendering speed, and image quality.


### 4.1. Time Consumption Comparison

The following is a comparison of several algorithms’ test data on the open dataset ISIC2018 [25]. Table 1 compares the drawing time of different algorithms for drawing the same picture of the human leg. VS3D, Dual-FHR, RIFT2, SHSR-UV, and ray-casting DAS in Figure 2 represent the experimental results of Wu’s algorithm [13], Luo’s algorithm [14], Francisco Sans’ algorithm [16], Roba Binyahib’s algorithm [17], and ray-casting DAS algorithm, respectively.Table 1
Comparison of the time consumption of each stage of the above algorithms (unit: seconds).
Algorithms
Data preprocessing time
Sampling time
Drawing time
Total time
VS3D
23
32
33
88
Dual-FHR
17
28
31
76
RIFT2
16
25
31
72
SHSR-UV
17
23
26
66
Ray-casting DAS
24
16
21
61Figure 2
The experimental results of the above algorithms on dataset ISIC2018. The experimental results show that compared with other algorithms, this algorithm can significantly reduce the sampling time required for image rendering, improve the rendering speed, and also improve the overall time consumption.
(a)
VS3D
(b)
Dual-FHR
(c)
RIFT2
(d)
SHSR-UV
(e)
Ray-casting DASAccording to the experimental results in Table1, compared with other algorithms, the preprocessing time of the algorithm in this paper is increased by about 6 to 7 seconds compared with other algorithms. This is because we add gradient calculations in the preprocessing stage, which consumes more time. However, due to the use of dynamic adaptive sampling, the acquisition of sampling points is more efficient, and the interference of unnecessary sampling points is eliminated. Compared with other algorithms, the time consumption of this algorithm in the sampling phase is greatly reduced, reducing 7-16 seconds. On the other hand, due to the improvement of the rendering operator, the idea of early termination is added in the rendering process, which greatly reduces the rendering time by 5-10 seconds. Finally, compared with other better algorithms, the total time is improved by 5-20 seconds, and the speed is significantly improved.Through the results of Figure2, we can find that the algorithm in this paper can better draw the details of the leg, including some texture information. These details are very valuable for medical diagnosis. Compared with other algorithms, the overall occlusion of the scene is captured coherently and close to the real value. The shadow area of the leg is less, and more texture details are preserved with ray-casting DAS.


### 4.2. Rendering Speed and Quality Comparison on CT Dataset

As shown in Figure3, RIFT2, SHSR-UV, ICVC-GPU, DOS-Ray, IMPA, and ray-casting DAS represent the experimental results of Francisco Sans’ algorithm [16], Roba Binyahibs’ algorithm [17], Feiniu Yuan’s algorithm [26], Leonardo’s algorithm [27], Yan Zhang’s algorithm [28], and ray-casting DAS algorithm, respectively. By comparing the results of the different algorithms on the dataset of ISBI2018 [29], the ray-casting DAS algorithm can reduce and avoid the phenomenon of artifact in CT image rendering and improve the visibility of the image.Figure 3
Comparison of the phenomenon of artifact and respective visibility of several algorithms.Table2 compares the rendering time of the ray-casting DAS algorithm with the other five algorithms. The experimental results show that ray-casting DAS has better performance in rendering time and better image signal-to-noise ratio than the previous algorithms.Table 2
Comparison of the rendering speed of various algorithms (drawing times/s, signal-to-noise ratio/dB).
Algorithms
Iteration
Drawing time
Signal-to-ratio
RIFT2 algorithm [16]
46
24
38
SHSR-UV algorithm [17]
53
29
42
ICVC-GPU algorithm [26]
34
21
33
DOS-ray algorithm [27]
49
25
34
IMPA algorithm [28]
45
22
36
Ray-casting DAS algorithm
32
15
60Similarly, from the experimental results in Figure3, it can be seen that the rendering algorithm in this paper can effectively solve the problem of artifacts in CT image rendering. From (a) to (e), we can clearly see the existence of artifacts in CT images, they obviously cause interference to the analysis image, and the image (f) drawn by ray-casting DAS algorithm can effectively eliminate the existence of these artifacts. It can be clearly seen from the graph that in the experimental results of this paper, each part of the image has no artifacts interference, and the image is clear.


### 4.3. Influence of Different Distance Functions on Image Rendering

The distance function expresses the relationship between distance and distance intensity value. The distance determines the distance intensity value. Through the distance function, the distance and depth are reflected in the final drawing process. This paper proposes an interactive algorithm. Interaction is based on distance first. The interaction process specifies the distance value, and the distance value is used in the distance function. Then, the final transparent value is determined by combining the distance function in the rendering process. The corresponding distance function is a piecewise function:(15)fx=high,x∈0,k,low,x∈k,max.The valuek is the distance value obtained by interaction. When the distance is less than k, it is a high-intensity value, and when the depth is greater than k, it is a low-intensity value.In order to improve the quality of rendering and reduce the interference of noise and other factors, we introduce the concept of distance function before color and opacity composition. The four contrast distance functions used in this experiment are as follows:(16)f1x=maxixi,1≤i≤p,xi≤52,minf1=0,(17)f2x=∑i=1pixi4+random0,1,1≤i≤p,xi≤52,minf2=0,(18)f3x=∑i=1pxi2−10cos2πxi,1≤i≤p,xi≤52,minf3=0,(19)f4x=−20∑i=1pexp−0.2130xi2−130exp∑i=1pcos2πxi+20+e,1≤i≤p,xi≤52,minf4=0.Here,xi is a three-dimensional column vector, representing the spatial coordinates of the ith voxel, and p is the total number of all effective voxels.As is shown in Table3, LCC [31] is short for linear correlation coefficient, and SROCC [32] is short for Spearman’s rank-order correlation coefficient, which are both evaluation indexes to measure the similarity between drawn image and real image. The structural similarity theory is simplified as SSIM [33], whose value is between 0 and 1. The larger LCC, SROCC, and SSIM value indicate the higher similarity between the drawn image and the real value of the original image. The peak signal-to-noise ratio (PSNR) [32] measures the ratio of useful information and noise in an image. For the same image, the larger the signal-to-noise ratio is, the smaller the noise is, and the higher the image quality is.Table 3
Comparison of several distance functions.
Algorithms
LCC
SPOCC
SSIM
MSSIM
PSNR
Without distance function
0.69
0.71
0.792
0.813
29 db
Ray-casting DAS with f1 function
0.72
0.76
0.814
0.839
33 db
Ray-casting DAS with f2 function
0.76
0.79
0.857
0.874
36 db
Ray-casting DAS with f3 function
0.75
0.78
0.849
0.867
35 db
Ray-casting DAS with f4 function
0.86
0.85
0.896
0.918
43 dbIn the experiment, four different distance functions are used, and their performances are compared. The experimental results in Table3 show that, compared with the method without distance function, adding distance function can significantly increase the similarity between the drawn image and the real image; at the same time, it has a higher signal-to-noise ratio and relatively less noise content. On the other hand, Figure 4 shows that different distance functions have advantages and disadvantages in image rendering. The distance function f4 is the best in the quality of completion. It can eliminate the interference of noise to the greatest extent. The image drawn is also the clearest and contains the least noise interference.Figure 4
Comparison of algorithms using different distance functions in the MICCAI [30] dataset.
(a)
Eithout distance function
(b)
f 1 function
(c)
f 2 function
(d)
f 3 function
(e)
f 4 function


### 4.4. Image Rendering Quality Comparison on Other Dataset

We use the open-source dataset VisMale [34] for experiments. The simulation environment is VS2016 and OpenGL. Table 4 shows the comparison of the six algorithms for human head image rendering, and ICVC-GPU, MVRC2, SigIg+FMI, SHSR-UV, DBRay, and Ray-casting DAS in Figure 5 represent the experimental results of Feiniu Yuan’s algorithm [26], Mohammandmehdi Bozorgi’s algorithm [32], R. Mehaboobathunnisa’s algorithm [35], Roba Binyahibs’ algorithm [17], Alec G.Moore’s algorithm [36], and ray-casting DAS algorithm proposed in this paper, respectively.Table 4
Comparison of completion time of several algorithms (unit: seconds).
Algorithms
Data preprocessing time
Sampling time
Drawing time
Total time consumption
ICVC-GPU
15
25
26
66
MVRC2
14
22
19
56
SigIg+FMI
13
24
32
69
SHSR-UV
15
24
33
72
DARay
17
26
23
66
Ray-casting DAS
22
19
15
57Figure 5
Comparison of the visual effects of the different algorithms in image rendering.
(a)
ICVC-GPU
(b)
MVRC2
(c)
SigIg+FMI
(d)
SHSR-UV
(e)
DBRay
(f)
Ray-casting DASIn Figure6, information entropy, contrast, sharpness, SNR [31, 33], and PSNR [32] are three objective indicators reflecting image quality. Their definitions are as follows:
(1)
Information entropy represents the amount of information contained in the image, which is proportional to the amount of information contained in the image [33]. The formula is(20)Ent=‐∑i=0L−1∑L−1pi,jlogΛpi,jΛ,where i represents the gray value of the pixel 0≤i≤255, j represents the mean gray value of the domain 0≤j≤255. And the formula pi,j=fi,j/N2 reflects the comprehensive characteristics of the gray value of a pixel position and the gray distribution of its surrounding pixels, where fi,j is the frequency of the feature binary i,j, and N is the scale of the image.
(2)
Contrast reflects the influence of image on visual effect [32]. Contrast is proportional to the clarity of the image. The formula is as follows:(21)con=L1MN∑x=0L−1∑y=0L−1Ix,y−1MN∑x=0L−1∑y=0L−1Ix,y2.Figure 6
Comparison of image quality generated by several algorithms.Among them,L denotes the average gray level of the pixels, M and N denote the width and height of the image, and Ix,y denotes the gray level of the pixels x,y.
(3)
Sharpness represents the clarity of the image, and the calculation formula is as follows(22)Def=∑x=0M−2∑y=0N−2Gx,y.(4)
Firstly, the mean square error (MSE) [37] is defined. The mean square error is used to calculate the mean square value of the pixel difference between the original image and the distorted image. The formula is as follows:(23)MSE=∑m=1M∑n=1NRm,n−Im,n2M×N.Peak signal-to-noise ratio (PSNR) [32] is the ratio of maximum signal quantity to noise intensity. Since digital images are represented by discrete numbers, generally L=255.
(24)PSNR=10lnl2MSE=10ln2552MSE.(5)
Signal to noise ratio (SNR)’s definition is as follows [31, 33](25)SNR=10lg∑m=1M∑n=1NRm,n2∑m=1M∑n=1NRm,n−Im,n2.The higher the SNR is, the less noise the image owns.Through the experimental results in Figure3, we can find that ray-casting DAS is slightly worse than algorithms E in image rendering quality, because the dynamic adaptive sampling used in this paper will inevitably skip some small voxels containing useful information due to the setting of the threshold value and other problems, leading to slightly poor effective information of the image, but it is close to or much better than A, B, C, and D in image rendering quality. Because the distance function is introduced into the process of rendering operator composition and the most appropriate distance function is selected, the image drawn by this algorithm has a high signal-to-noise ratio and contains less noise and other interference factors, compared with other algorithms. This is a significant progress. Simultaneously, the sampling and drawing time of ray-casting DAS, as shown in Table 3, is much better than all the other five algorithms, which indicates that our method is effective and saving time consumption.


## 4.1. Time Consumption Comparison

The following is a comparison of several algorithms’ test data on the open dataset ISIC2018 [25]. Table 1 compares the drawing time of different algorithms for drawing the same picture of the human leg. VS3D, Dual-FHR, RIFT2, SHSR-UV, and ray-casting DAS in Figure 2 represent the experimental results of Wu’s algorithm [13], Luo’s algorithm [14], Francisco Sans’ algorithm [16], Roba Binyahib’s algorithm [17], and ray-casting DAS algorithm, respectively.Table 1
Comparison of the time consumption of each stage of the above algorithms (unit: seconds).
Algorithms
Data preprocessing time
Sampling time
Drawing time
Total time
VS3D
23
32
33
88
Dual-FHR
17
28
31
76
RIFT2
16
25
31
72
SHSR-UV
17
23
26
66
Ray-casting DAS
24
16
21
61Figure 2
The experimental results of the above algorithms on dataset ISIC2018. The experimental results show that compared with other algorithms, this algorithm can significantly reduce the sampling time required for image rendering, improve the rendering speed, and also improve the overall time consumption.
(a)
VS3D
(b)
Dual-FHR
(c)
RIFT2
(d)
SHSR-UV
(e)
Ray-casting DASAccording to the experimental results in Table1, compared with other algorithms, the preprocessing time of the algorithm in this paper is increased by about 6 to 7 seconds compared with other algorithms. This is because we add gradient calculations in the preprocessing stage, which consumes more time. However, due to the use of dynamic adaptive sampling, the acquisition of sampling points is more efficient, and the interference of unnecessary sampling points is eliminated. Compared with other algorithms, the time consumption of this algorithm in the sampling phase is greatly reduced, reducing 7-16 seconds. On the other hand, due to the improvement of the rendering operator, the idea of early termination is added in the rendering process, which greatly reduces the rendering time by 5-10 seconds. Finally, compared with other better algorithms, the total time is improved by 5-20 seconds, and the speed is significantly improved.Through the results of Figure2, we can find that the algorithm in this paper can better draw the details of the leg, including some texture information. These details are very valuable for medical diagnosis. Compared with other algorithms, the overall occlusion of the scene is captured coherently and close to the real value. The shadow area of the leg is less, and more texture details are preserved with ray-casting DAS.


## 4.2. Rendering Speed and Quality Comparison on CT Dataset

As shown in Figure3, RIFT2, SHSR-UV, ICVC-GPU, DOS-Ray, IMPA, and ray-casting DAS represent the experimental results of Francisco Sans’ algorithm [16], Roba Binyahibs’ algorithm [17], Feiniu Yuan’s algorithm [26], Leonardo’s algorithm [27], Yan Zhang’s algorithm [28], and ray-casting DAS algorithm, respectively. By comparing the results of the different algorithms on the dataset of ISBI2018 [29], the ray-casting DAS algorithm can reduce and avoid the phenomenon of artifact in CT image rendering and improve the visibility of the image.Figure 3
Comparison of the phenomenon of artifact and respective visibility of several algorithms.Table2 compares the rendering time of the ray-casting DAS algorithm with the other five algorithms. The experimental results show that ray-casting DAS has better performance in rendering time and better image signal-to-noise ratio than the previous algorithms.Table 2
Comparison of the rendering speed of various algorithms (drawing times/s, signal-to-noise ratio/dB).
Algorithms
Iteration
Drawing time
Signal-to-ratio
RIFT2 algorithm [16]
46
24
38
SHSR-UV algorithm [17]
53
29
42
ICVC-GPU algorithm [26]
34
21
33
DOS-ray algorithm [27]
49
25
34
IMPA algorithm [28]
45
22
36
Ray-casting DAS algorithm
32
15
60Similarly, from the experimental results in Figure3, it can be seen that the rendering algorithm in this paper can effectively solve the problem of artifacts in CT image rendering. From (a) to (e), we can clearly see the existence of artifacts in CT images, they obviously cause interference to the analysis image, and the image (f) drawn by ray-casting DAS algorithm can effectively eliminate the existence of these artifacts. It can be clearly seen from the graph that in the experimental results of this paper, each part of the image has no artifacts interference, and the image is clear.


## 4.3. Influence of Different Distance Functions on Image Rendering

The distance function expresses the relationship between distance and distance intensity value. The distance determines the distance intensity value. Through the distance function, the distance and depth are reflected in the final drawing process. This paper proposes an interactive algorithm. Interaction is based on distance first. The interaction process specifies the distance value, and the distance value is used in the distance function. Then, the final transparent value is determined by combining the distance function in the rendering process. The corresponding distance function is a piecewise function:(15)fx=high,x∈0,k,low,x∈k,max.The valuek is the distance value obtained by interaction. When the distance is less than k, it is a high-intensity value, and when the depth is greater than k, it is a low-intensity value.In order to improve the quality of rendering and reduce the interference of noise and other factors, we introduce the concept of distance function before color and opacity composition. The four contrast distance functions used in this experiment are as follows:(16)f1x=maxixi,1≤i≤p,xi≤52,minf1=0,(17)f2x=∑i=1pixi4+random0,1,1≤i≤p,xi≤52,minf2=0,(18)f3x=∑i=1pxi2−10cos2πxi,1≤i≤p,xi≤52,minf3=0,(19)f4x=−20∑i=1pexp−0.2130xi2−130exp∑i=1pcos2πxi+20+e,1≤i≤p,xi≤52,minf4=0.Here,xi is a three-dimensional column vector, representing the spatial coordinates of the ith voxel, and p is the total number of all effective voxels.As is shown in Table3, LCC [31] is short for linear correlation coefficient, and SROCC [32] is short for Spearman’s rank-order correlation coefficient, which are both evaluation indexes to measure the similarity between drawn image and real image. The structural similarity theory is simplified as SSIM [33], whose value is between 0 and 1. The larger LCC, SROCC, and SSIM value indicate the higher similarity between the drawn image and the real value of the original image. The peak signal-to-noise ratio (PSNR) [32] measures the ratio of useful information and noise in an image. For the same image, the larger the signal-to-noise ratio is, the smaller the noise is, and the higher the image quality is.Table 3
Comparison of several distance functions.
Algorithms
LCC
SPOCC
SSIM
MSSIM
PSNR
Without distance function
0.69
0.71
0.792
0.813
29 db
Ray-casting DAS with f1 function
0.72
0.76
0.814
0.839
33 db
Ray-casting DAS with f2 function
0.76
0.79
0.857
0.874
36 db
Ray-casting DAS with f3 function
0.75
0.78
0.849
0.867
35 db
Ray-casting DAS with f4 function
0.86
0.85
0.896
0.918
43 dbIn the experiment, four different distance functions are used, and their performances are compared. The experimental results in Table3 show that, compared with the method without distance function, adding distance function can significantly increase the similarity between the drawn image and the real image; at the same time, it has a higher signal-to-noise ratio and relatively less noise content. On the other hand, Figure 4 shows that different distance functions have advantages and disadvantages in image rendering. The distance function f4 is the best in the quality of completion. It can eliminate the interference of noise to the greatest extent. The image drawn is also the clearest and contains the least noise interference.Figure 4
Comparison of algorithms using different distance functions in the MICCAI [30] dataset.
(a)
Eithout distance function
(b)
f 1 function
(c)
f 2 function
(d)
f 3 function
(e)
f 4 function


## 4.4. Image Rendering Quality Comparison on Other Dataset

We use the open-source dataset VisMale [34] for experiments. The simulation environment is VS2016 and OpenGL. Table 4 shows the comparison of the six algorithms for human head image rendering, and ICVC-GPU, MVRC2, SigIg+FMI, SHSR-UV, DBRay, and Ray-casting DAS in Figure 5 represent the experimental results of Feiniu Yuan’s algorithm [26], Mohammandmehdi Bozorgi’s algorithm [32], R. Mehaboobathunnisa’s algorithm [35], Roba Binyahibs’ algorithm [17], Alec G.Moore’s algorithm [36], and ray-casting DAS algorithm proposed in this paper, respectively.Table 4
Comparison of completion time of several algorithms (unit: seconds).
Algorithms
Data preprocessing time
Sampling time
Drawing time
Total time consumption
ICVC-GPU
15
25
26
66
MVRC2
14
22
19
56
SigIg+FMI
13
24
32
69
SHSR-UV
15
24
33
72
DARay
17
26
23
66
Ray-casting DAS
22
19
15
57Figure 5
Comparison of the visual effects of the different algorithms in image rendering.
(a)
ICVC-GPU
(b)
MVRC2
(c)
SigIg+FMI
(d)
SHSR-UV
(e)
DBRay
(f)
Ray-casting DASIn Figure6, information entropy, contrast, sharpness, SNR [31, 33], and PSNR [32] are three objective indicators reflecting image quality. Their definitions are as follows:
(1)
Information entropy represents the amount of information contained in the image, which is proportional to the amount of information contained in the image [33]. The formula is(20)Ent=‐∑i=0L−1∑L−1pi,jlogΛpi,jΛ,where i represents the gray value of the pixel 0≤i≤255, j represents the mean gray value of the domain 0≤j≤255. And the formula pi,j=fi,j/N2 reflects the comprehensive characteristics of the gray value of a pixel position and the gray distribution of its surrounding pixels, where fi,j is the frequency of the feature binary i,j, and N is the scale of the image.
(2)
Contrast reflects the influence of image on visual effect [32]. Contrast is proportional to the clarity of the image. The formula is as follows:(21)con=L1MN∑x=0L−1∑y=0L−1Ix,y−1MN∑x=0L−1∑y=0L−1Ix,y2.Figure 6
Comparison of image quality generated by several algorithms.Among them,L denotes the average gray level of the pixels, M and N denote the width and height of the image, and Ix,y denotes the gray level of the pixels x,y.
(3)
Sharpness represents the clarity of the image, and the calculation formula is as follows(22)Def=∑x=0M−2∑y=0N−2Gx,y.(4)
Firstly, the mean square error (MSE) [37] is defined. The mean square error is used to calculate the mean square value of the pixel difference between the original image and the distorted image. The formula is as follows:(23)MSE=∑m=1M∑n=1NRm,n−Im,n2M×N.Peak signal-to-noise ratio (PSNR) [32] is the ratio of maximum signal quantity to noise intensity. Since digital images are represented by discrete numbers, generally L=255.
(24)PSNR=10lnl2MSE=10ln2552MSE.(5)
Signal to noise ratio (SNR)’s definition is as follows [31, 33](25)SNR=10lg∑m=1M∑n=1NRm,n2∑m=1M∑n=1NRm,n−Im,n2.The higher the SNR is, the less noise the image owns.Through the experimental results in Figure3, we can find that ray-casting DAS is slightly worse than algorithms E in image rendering quality, because the dynamic adaptive sampling used in this paper will inevitably skip some small voxels containing useful information due to the setting of the threshold value and other problems, leading to slightly poor effective information of the image, but it is close to or much better than A, B, C, and D in image rendering quality. Because the distance function is introduced into the process of rendering operator composition and the most appropriate distance function is selected, the image drawn by this algorithm has a high signal-to-noise ratio and contains less noise and other interference factors, compared with other algorithms. This is a significant progress. Simultaneously, the sampling and drawing time of ray-casting DAS, as shown in Table 3, is much better than all the other five algorithms, which indicates that our method is effective and saving time consumption.


## 5. Conclusion

Based on the low efficiency of the sampling point selection and slow rendering speed of the existing ray-casting algorithm, a new method based on dynamic sampling and improved rendering operator is proposed in this paper. The experimental results show that this method can effectively shorten the time of sampling and rendering and improve the efficiency of the algorithm on the premise of ensuring the high quality of graphics rendering.At the same time, this paper introduces the concept of distance function in the process of image rendering, which makes the image contrast high and contains less noise. However, there are other shortcomings in this algorithm, such as the need for gradient calculation, and the gradient calculation is placed in the preprocessing stage, which aggravates the calculation amount in the preprocessing stage. It results in more time in the pretreatment stage. The efficiency of this algorithm can be further improved by a more efficient gradient calculation method.



---

*Source: 8822624-2020-10-14.xml*