Plane-Based Sampling for Ray Casting Algorithm in Sequential Medical Images

This paper proposes a plane-based sampling method to improve the traditional Ray Casting Algorithm (RCA) for the fast reconstruction of a three-dimensional biomedical model from sequential images. In the novel method, the optical properties of all sampling points depend on the intersection points when a ray travels through an equidistant parallel plan cluster of the volume dataset. The results show that the method improves the rendering speed at over three times compared with the conventional algorithm and the image quality is well guaranteed.


Introduction
Modeling three-dimensional (3D) volume of biomedical tissues from 2D sequential images is an important technique to highly improve the diagnostic accuracy [1]. Volume rendering refers to the process that maps the 3D discrete digital data into image pixel values [2]. It can be classi�ed into two categories: one is direct volume rendering which generates images by compositing pixel values along rays cast into a 3D image, and the other one is indirect volume rendering which visualizes geometry element graphics extracted from the volume data [3]. e importance of volume rendering is resampling and synthesizing image [4]. Ray casting, splatting, and shear-warp are the three popular volume rendering algorithms now [5].
Ray Casting Algorithm (RCA) is a direct volume rendering algorithm. e traditional RCA is widely used for it can precisely visualize various medical images with details of boundary and internal information from sequential images, while real-time rendering with traditional RCA is still an obstacle due to its huge computation.
In recent years, numerous techniques have been proposed to accelerate the rendering speed. In general, there are three primary aspects, including hardware-based, parallel, and soware-based acceleration algorithms. Liu et al. [6] proposed a method combined that Graphics Processing Unit (GPU) and octree encoding and accelerated RCA at a rate of 85 times. Wei and Feng [7] presented a GPU-based real-time ray casting method for algebraic B-spline surfaces via iterative root-�nding algorithms. �hang et al. [8] accelerated RCA on Compute Uni�ed Device Architecture (CUDA), which can perform more samplings within a ray segment using cubic Bspline.
However, both hardware-based and parallel techniques are inseparable from the development of computer hardware. By comparison, soware-based algorithms can be quickly transplanted among different machines. What is more, they can show �exibility of the procedure and re�ect the thoughts of researchers. Yang et al. [9] sampled points based on all intersection points at which the ray transacts with the voxel. All intersections in a voxel depend on four vertexes on one face. However, the condition whether two intersection points were on adjacent or opposite surface in a voxel was neglected. Ling and Qian [10] used a bounding volume method to avoid casting the viewing rays that do not intersect with the volume. Since such situation can be judged quickly by comparing the world coordinates of sampling point with the volume dataset, it did not obviously speed up the rendering process.
Recently, Qian et al. [11] replaced the sampling points with intersection points when rays travel through three groups of parallel planes along three orthometric axes to reduce the rendering time. However, it cannot guarantee the image density when the distance between adjacent parallel planes far surpasses the sampling interval. is paper proposes an improved RCA to speed the rendering process. e main idea is, when the ray travels through one group of equidistant parallel planes of the volume, intersection points are obtained. en the properties of sampling points between adjacent intersection points can be calculated by the formula of de�nite proportion and separated points. By this method, a small number of intersection points are considered; meanwhile the method does not sacri�ce the sampling density.

Ray Casting Algorithm
Overview. e traditional RCA involves two steps: (1) assign optical properties such as color and opacity to all 3D discrete vertexes according to their gray value, and (2) apply a sampling and composing process. For each output image pixel in sequence, do the following.
(i) Cast the ray through the volume from back to front.
(ii) Sample the color and opacity at each regular sampling point along the ray. (1) e rendering time is mainly comprised of four parts in the above-mentioned rendering process [11]. ey are converting gray value into optical property (about 30%), computing position of sampling points (about 3%), sampling optical properties (about 39%), and compositing properties into output pixel color (about 6%). e time proportion of sampling is the highest. Moreover, the time ratio of four parts is not constant. e greater the sampling data is, the larger the proportion of sampling time is. erefore, sampling has a direct impact on speed of RCA.

Traditional Sampling Method.
Traditionally, the optical property of each sampling point depends on eight vertexes of its voxel by trilinear interpolation [12,13]. In detail, there are four steps for the sampling one point. First, locate its voxel and convert the world coordinates of sampling point into voxel's local coordinates. e following three steps are processes of linear interpolations along three different axes in order. e interpolation diagram of Ray Casting Algorithm is shown in Figure 1.
For example, to sample point ( ) in white circle (Figure 1), �rst obtain the voxel ( ) and local coordinates ( ) of , which are expressed in (2). en the optical property of four points ( 2 3 4 ) on the plane through is deduced according to eight vertexes ( ∼ 8 ) along zaxis. e next property of two points ( 5 6 ) forming the line segment through is computed along -axis. At last is obtained along -axis by de�nite proportional division point formula.
In Figure 1, assume the pixel spacing along -, -, -axes is Δ , Δ , and Δ , respectively, with ( ): where operator [⋅] represents taking the �oor integral. e property of can be calculated by 5 and 6 , which are obtained by 2 3 , and 4 . e relationship between them is shown in = + According to the above equations, 17 additions and 16 multiplications are executed for sampling each point such as (see Figure 1), including 3 additions and 9 multiplications to locate the voxel ( ) and get the local coordinates. In Figure  1, there are 6 sampling points in two voxels, 102 additions, and 96 multiplications performed. To simplify the calculation of sampling process, a new RCA based on plane clusters sampling is proposed.

Proposed Plan-Based Sampling
Method. e basic idea of the plan-based sampling method is to acquire all sampling points based on intersection points when ray travels through a group of parallel planes in the volume data �eld.
e sampling process, speci�cally, consists of three steps. First, intersections and the corresponding plane are obtained based on some necessary initial conditions. en the optical property of all the intersection points is obtained by linear interpolation according to vertexes on plane clusters. e optical property of sampling points between intersection points along the ray is computed by de�nite proportion and separated point formula.
Assuming that the direction vector of ray is and the extent of gridding volume data is , with the spacing Δ Δ Δ along ---axes, respectively, the three plane clusters are as follows: Parallel plane clusters along axis are selected. Let the origin point of ray be . e ray intersects with plane at entry point and belongs to the voxel . e coordinates of and voxel are deduced next. e derivation is shown as follows. Since where means the distance from to along ray, the value of can be obtained from erefore, Δ and of can be expressed as follows: Considering that belongs to voxel , then and are expressed as follows: erefore, when is given, , and can be obtained through the above equations.
From the mathematical derivation, when original position, direction vector, and the extent of volume data are given, all the intersections and associated voxels can be quickly obtained.
In Figure 1, the property of entry point can be computed by the property ( 3 ) of three vertexes on voxel , that is, In the same way, the property of exit point can be obtained. At last the property is expressed as follows: In addition, when one component of the direction vector is zero, a plane cluster along another axis can be chosen. If two components are zero, the plane clusters along the third axis are taken into account.

Comparison of Two Sampling Methods.
In the new RCA sampling process, only intersection points on a plane cluster along one axis need to be considered without converting coordinates. While in the conventional sampling process, the world coordinates of each sampling point are converted into voxel's local coordinates and computed by trilinear interpolation [14,15].
As is shown in Figure 1, there are 6 sampling points between and . 15 additions and 19 multiplications are executed to sample and , and 24 additions and 12 multiplications are run to sample six points based on and . Totally, 39 additions and 31 multiplications are taken compared with 102 additions and 96 multiplications with trilinear interpolation. Furthermore, not all vertexes are referred because some vertexes (such as 4 7 in Figure 1) are not used as reference by the new method. us, in theory, the calculation amount is reduced to less than one third on the whole.  a question whether the image quality can be guaranteed. It can be seen in Figures 2 and 3 that images reconstructed by RCA based on plan cluster sampling method are almost the same as those based on traditional trilinear interpolation in RCA. ey can clearly show the details of the boundary and internal information of the volume with the new sampling method. erefore, the image quality can be well ensured. By comparing the amount of computation (39/102-31/ 96) in the two sampling methods, the new method can reduce the amount of traditional one to about one third. It can be seen that the total rendering time ( Table 1) using new method is less than one third of that using conventional trilinear interpolation. It indicates that the time saved to inquire the property of the vertexes not for reference should not be underestimated.

Experiments and Analysis
Moreover, it is shown that the acceleration rate of the head images is higher than that of the heart images. e main difference between them is that the spacing of head CT sequences is denser than the heart data. erefore, the denser the data is, the more efficient the new method is.

Conclusion
is paper presented a novel RCA based on a parallel plan cluster sampling method. e proposed method can efficiently speed up the sampling process at more than three times and still clearly display the boundary and internal information of the volume; thus the image quality is well guaranteed. In addition, the comparison of acceleration rate indicates that the new method is more effective for dataset with denser spacing. e new method can meet the real-time requirements of interactive rendering.