Abstract
Spinal pathology treatment has become an urgent issue to be solved. How to effectively prevent and treat spinal pathology has become a research hotspot in the field of surgery. Aiming at the problem of too long volume rendering time caused by the trilinear interpolation sampling method in the reconstruction and visualization of the vertebra 3D model, an improved ray projection algorithm is proposed to quickly reconstruct a 3D vertebra model from medical CT vertebra images. This method first classifies CT data, assigns corresponding color values and opacity transfer functions to different types of data, and then uses inverse distance-weighted interpolation (IDWI) sampling to replace the trilinear interpolation sampling method for the voxel where the sampling point is located to accelerate the interpolation operation. The color value and opacity of the sampling points are obtained, and finally, the attributes of all the sampling points are synthesized and calculated to obtain the final rendering effect, and the reconstruction of the three-dimensional vertebra model is completed. Experimental results show that the proposed method not only can obtain higher quality rendered images but also has a certain improvement in rendering speed compared with traditional algorithms.
1. Introduction
Nowadays, spine-related diseases are a major problem in modern human society. Modern clinical medicine points out that spinal pathology has become an urgent problem to be solved. Pedicle screw placement is a relatively mature medical method for spinal pathology, and three-dimensional reconstruction of vertebrae is the premise of preoperative simulation and quantification. The three-dimensional vertebral model obtained by medical image visualization can provide more realistic anatomical structure for doctors, so as to help doctors accurately diagnose the disease and design treatment plan, which has important clinical practical value.
Medical image visualization [1] is the human body information obtained by digital imaging technology that is intuitively expressed as three-dimensional effect on the computer, so as to provide structural information that cannot be obtained by traditional means. Volume rendering has become an important method to realize the visualization of medical images because it can show the fine structure and small transformation of objects. As a classical method of volume rendering technology, ray casting algorithm has a wide range of applications in the field of 3D reconstruction visualization of medical images, but its rendering speed is slightly insufficient [2]. At present, the common methods to improve the rendering speed can be roughly divided into three types: hardware-based, software-based, and parallel mode [3]. Among them, the premise of hardware-based and parallel methods is to carry out on specific computer hardware, so its portability is greatly reduced. The acceleration method based on software type does not rely too much on the development of hardware, but improves from the angle of algorithm, which has high flexibility. However, this kind of method has some problems, such as too high complexity to achieve ideal acceleration effect and not taking into account both rendering speed and reconstruction quality.
In order to speed up volume rendering and improve the quality of reconstruction, we proposed a ray casting algorithm based on inverse distance-weighted interpolation sampling to realize the 3D visualization of the vertebral CT image. Firstly, vertebrae and nonvertebrae are distinguished in the 3D volume data field, and the mapping from data to optical features is realized by transfer function. Then, the inverse distance-weighted interpolation method is used to replace the trilinear interpolation in the traditional ray casting algorithm to speed up the resampling process. Finally, the final color of each pixel is obtained by synthesizing the color value and opacity of all the sampling points, and finally, the three-dimensional reconstruction of the vertebral image is realized.
2. Related Work
With the development of computer hardware technology, hardware-based and parallel acceleration methods are first proposed. Cullip and Neumann first proposed a method to accelerate volume rendering using 3D texture hardware [4]. Zhang et al. used GPU-based fast ray casting method for volume rendering in CT 3D reconstruction to shorten the reconstruction time [5]. Ross et al. proposed a CPU-based volume rendering algorithm for 3D ultrasound images, which overcomes the difficulty of low adaptability of GPU-based algorithm [6]. Ma et al. proposed a parallel grid generation algorithm on GPU, which increased the average efficiency by 15 times [7]. Zhou et al. used CUDA to accelerate the implementation of ray casting algorithm, which was 70% faster than GPU [8]. Sans et al. compared the performance of OpenGL, OpenCL, and CUDA for different medical datasets [9].
However, both hardware-based and parallel-based acceleration methods are based on the premise of computer hardware, so they have some limitations. In contrast, the acceleration method based on software is more flexible and convenient, which can be transplanted between different machines quickly and has wider applicability. Mehaboobathunnisa et al. proposed a method of grouping rays projected by similar voxels to reduce the computational complexity of rendering algorithm, but the reconstruction result is not smooth enough due to artifacts [10]. Hadwiger et al. proposed the SparseLeap method which is a novel space hopping method and has been proved that it can avoid the problem of unnecessary space debris [11]. Based on the idea of spatial jump method, Deakin and Knackstedt used Chebyshev distance to guide how to effectively skip the blank area in the ray casting algorithm [12], and then they optimized this method and proposed an effective algorithm to generate anisotropic Chebyshev distance map for accelerating ray casting, but these two methods have some shortcomings in the final reconstruction effect [13]. Liu et al. reduced the amount of computation in the rendering process by adjusting the sampling frequency and adopting different interpolation strategies for the sampling points of different classification groups [14]. Bi et al. applied inverse distance-weighted interpolation algorithm to meteorological data visualization and achieved good rendering effect [15]. In order to improve the speed of ray casting algorithm, this paper proposes a ray casting algorithm based on inverse distance-weighted interpolation sampling, which improves the speed of volume rendering on the premise of meeting the quality requirements of 3D reconstruction.
3. Improved Ray Casting Algorithm
3.1. Traditional Ray Casting Algorithm
The basic idea of the traditional ray casting algorithm can be described as follows: firstly, collecting the sample points of all the voxels in the three-dimensional data field along the ray direction at equal intervals and then obtaining the corresponding color values and opacity values of the sample points, then synthesizing the color values and opacity values of all the sample points on the ray to obtain the final color of the pixel, and finally, calculating each pixel and obtaining the final two-dimensional image. The detailed steps are shown in Table 1.
In traditional ray casting algorithm, step 1 is to initialize the color value C and opacity value α. In step 2, all the voxels in the 3D volume data field are traversed, and the sample points are sampled by trilinear interpolation to obtain the color value and opacity value of the sample points. Finally, the color values and opacity values of all sampling points on the ray are combined in step 3 to get the final projection image.
3.2. Inverse Distance-Weighted Interpolation Method
Inverse distance-weighted interpolation (IDWI) is a kind of interpolation method which takes the distance between the sampling point and adjacent point as weight. Figure 1 shows the schematic diagram of the IDWI method.
The IDWI method considers that each adjacent point will have a certain influence on the sampling point, and the influence is closely related to the distance. The closer the sampling point is, the greater the weight is given to the adjacent point, and the weight contribution is inversely proportional to the distance [16].
Suppose that there are other neighboring points in the neighborhood of sampling point A(x, y, z), denoted as Ai(xi, yi, zi), i = 1, 2, …, n. Let f(A) be the color and opacity value of A, then it can be described as follows:where is the weight of the distance from each adjacent point to the sampling point and can be calculated as follows:where is the Euclidean distance between adjacent points and sampling points and can be calculated as follows:
Obviously, the sum of all weights is 1, that is,
3.3. Ray Casting Algorithm Based on IDWI Sampling
In order to improve the efficiency of the TRC, an improved ray casting algorithm is proposed in this paper. The basic idea is to use inverse distance-weighted interpolation instead of trilinear interpolation in interpolation sampling. The improved ray casting algorithm includes the following three steps:(1)Classification of vertebrae and nonvertebrae: in this step, a three-dimensional volume data field is constructed, and a ray is emitted from the position of each pixel on the screen along the line of sight direction into the data field. According to the different characteristics of the data, the vertebral CT data are classified into different types, that is, vertebral and nonvertebral. In order to show the internal structure of 3D data field more intuitively, it is necessary to assign different color values and opacity to different types of data according to the classification results. The setting of color value and opacity needs to design different transfer functions. The function of transfer function is to complete the mapping from data to optical features (color, opacity, etc.). Formally, the transfer function T can be defined as follows: where is a binary set of color values and opacity, x is the attribute values of volume data, and n is the dimension of x and represents the number of attributes.(2)Achieving the color and opacity values: after separating the region of vertebrae and nonvertebrae, the sample points of all the voxels in the 3D data field which are penetrated by light are collected at equal intervals along the ray direction, and the corresponding color and opacity values of the sample points are calculated according to the attribute values of the eight vertices of the voxels where the sample points are located. The sampling points are above the ray incident from the pixel to the 3D data field, usually in the voxel mesh they pass through, rather than just at the vertex of the voxel mesh. The inverse distance-weighted interpolation method is used to obtain the color value and opacity of the sampling points.(3)Synthesizing the color values: in this step, the final color of the pixels is obtained by compositing the color values and opacity values of all the sampling points on the ray. In the process of synthesis, we use the front-to-back strategy. The formulas can be described as follows: where and are the color values and opacity values after the cast ray passes through the sample point, respectively. and are the color values and opacity values before the incident sampling point i, respectively. and are the color values and opacity values of the current sampling point i, respectively. After the above iteration, each pixel is calculated to get the final three-dimensional image. The detailed steps are shown in Table 2.
3.4. Stability Analysis
In the process of image rendering using the classical ray casting algorithm, there are two kinds of situations in which the projected ray travels in the 3D volume data field: the ray travels in the empty voxel and the ray travels in the nonempty voxel. For analyzing the whole rendering process from the rendering data source, the calculation of empty voxels has no contribution to the final 3D reconstruction effect presented on the screen. However, if the same undifferentiated interpolation is applied to the empty voxels, the rendering time of the ray casting algorithm will be increased, and the real-time performance of the algorithm will be affected. Therefore, we use the spatial jump technique to skip the empty voxels in the direction of the projection ray and only interpolate the nonempty voxels in the 3D volume data field. For each empty voxel, first record the nearest distance from it to the opaque voxel, and then, no matter what the direction of the ray is, as long as the distance is a forward step, it will not intersect with the opaque voxel in this distance range when the ray is casting. Based on this, it can significantly reduce the redundant steps of ray casting, improve the efficiency of the algorithm, and has obvious advantages in computing resources and memory requirements.
4. Experiments and Analysis
4.1. Experimental Setup
The experimental data were obtained from the project cooperation Hospital of Xi’an Zhenwo 3D Technology Co., Ltd., with a total of 77 samples and a resolution of 512 × 512 DICOM format CT image data file, and each CT image interval is 2.5 mm. Hardware environment is Intel (R) core (TM) i5-4590 @ 3.30 GHz, 8 GB, AMD Radon (TM) r5340x, OS is windows10, and all programs are implemented in Python + VTK development environment. In this paper, five vertebrae of lumbar vertebrae are taken as the reconstruction object, and the DICOM format vertebral CT data is used for experiment. The time required by several ray casting algorithms for volume rendering is compared, and the PSNR and SSMI values are used as evaluation indexes to compare the reconstruction quality.
4.2. Reconstruction Results and Analysis
In the same experimental environment, two groups of experiments were carried out: the first group took all slices of lumbar CT data as input, compared with the traditional ray casting (TRC), the ray casting based on bounding box optimization (BRC) [17], the ray casting with viewpoint (VRC) [18], and our method (improved ray casting, IRC), the 3D reconstruction results and execute time are given; in the second group, part of CT slices in the data set was taken to reconstruct a single vertebra, and the time-consuming and reconstruction results of the four volume rendering methods were compared.
Figures 2 and 3 show the comparison of reconstruction effect of lumbar and single vertebral, respectively. In addition, Tables 3 and 4 show the comparison of reconstruction quality indexes to TRC of lumbar and single vertebra, respectively.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
From Figures 2 and 3, we can see that our method IRC algorithm has no obvious difference in reconstruction quality compared with the other three algorithms. However, from the reconstruction quality indicators in Tables 3 and 4, we can see that the PSNR value of IRC is little different from the other algorithms, and the SSIM value is very close. This shows that our improved method not only ensures the quality of 3D reconstruction but also shows the better reconstruction model. Table 5 shows the running time of the four algorithms in two groups of experiments.
From Table 5, we can see that the running time of TRC, BRC, and VRC algorithm for single vertebra is 10.739 s, 8.254 s, and 7.851 s, respectively, which is slower than 6.473 s of our improved ray casting algorithm. Meanwhile, the running time of IRC for whole lumbar spine is 14.346 s, which is 9.979 s, 4.516 s, and 2.162 s faster than TRC (24.325 s), BRC (18.862 s), and VRC (16.508 s), respectively. Therefore, IRC is better than the other three algorithms in running speed for CT data files of the same experimental object. For each interpolation point, the trilinear interpolation method needs 21 times of multiplication and division and 28 times of addition and subtraction, while the inverse distance-weighted interpolation method only has 16 times of multiplication and division and 7 times of addition and subtraction. Because of the reduction of the amount of computation, the IRC shortens the running time.
5. Conclusion
This paper proposed an improved ray casting algorithm to solve the problem of too long rendering time of traditional ray casting algorithm in the process of 3D reconstruction of vertebral CT image. According to the shortcomings of trilinear interpolation method in sampling speed, inverse distance-weighted interpolation sampling is used instead of trilinear interpolation sampling. Compared with the other existing methods, our method can reduce the reconstruction time and improve the rendering speed without reducing the reconstruction quality. How to improve the existing methods so as to obtain more fine reconstruction effect of microstructure and further improve the reconstruction accuracy is the next research direction.
Data Availability
The datasets used in this paper are publicly available.
Conflicts of Interest
The datasets used in this paper are available from the corresponding author.
Acknowledgments
This paper was funded by the National Key Research and Development Program of China (no. 2017YFB1401800), the Philosophy and Social Sciences Research Planning Project of Heilongjiang Province (nos. 20GLB119 and 19GLB327), and the Talents Plan of Harbin University of Science and Technology: Outstanding Youth Project (no. 2019-KYYWF-0216).