Abstract
In order to efficiently and effectively reconstruct 3D medical images and clearly display the detailed information of inner structures and the inner hidden interfaces between different media, an Improved Volume Rendering Optical Model (IVROM) for medical translucent volume rendering and its implementation using the preintegrated Shear-Warp Volume Rendering algorithm are proposed in this paper, which can be readily applied on a commodity PC. Based on the classical absorption and emission model, effects of volumetric shadows and direct and indirect scattering are also considered in the proposed model IVROM. Moreover, the implementation of the Improved Translucent Volume Rendering Method (ITVRM) integrating the IVROM model, Shear-Warp and preintegrated volume rendering algorithm is described, in which the aliasing and staircase effects resulting from under-sampling in Shear-Warp, are avoided by the preintegrated volume rendering technique. This study demonstrates the superiority of the proposed method.
1. Introduction
In medical visualization applications, such as computer aided diagnosis and therapy, the relative position between the tumor and its adjacent tissues is often determined by physicians through observing the inner pathologic changes of tissue structures and the detailed information of the tissue surface, before making a reasonable therapy planning. The detailed information of inner structures and inner hidden interfaces between different layers can be clearly displayed by using the translucent volume rendering technique, which has been widely applied in medical visualization.
Lighting Absorption and Emission Model combined with Blinn-Phong surface shading model (LAEM-BP) is one of the most widely used lighting models in direct volume rendering in medical visualization [1, 2]. The optical model was described in papers by Sabella [3], Drebin et al. [4], and Levoy [5], in which the volume shading incorporating diffuse and specular shading by approximating the surface normal with the gradient of the 3D field is described, but scattering is ignored in order to achieve fast approximation used by direct lighting. Moreover, in order to realize interactive volume rendering of 3D scalar fields of medical data, the optical model as well as its volume rendering method is implemented in hardware [6, 7]. The model can approximately render obvious translucent tissue surface where the normal is well defined for regions in the volume that have high gradient magnitudes. However, this method cannot clearly display the detailed information of inner structures and the inner hidden interfaces between different media, especially when the scattering becomes the main factor that affects the rending effects and dominates the visual appearance [2].
The effect of multiple scattering and indirect illumination is important for volume rendering applications [1, 2, 8, 9]. Blinn [10] presents a model for the reflection and transmission of light through thin clouds of particles based on probabilistic arguments and single scattering approximation. Kajiya and Von Herzen describe a model for rendering arbitrary volume densities that includes expensive multiple scattering computation [11]. In their methods, the radiative transport equation [12] is employed, and expensive and sophisticated numerical methods must also be employed to compute the radiance distribution to a desired accuracy.
Furthermore, Max et al. [13] use a one-dimensional scattering equation to calculate the light transport, where forward peaked phase function as the hemisphere is needed to be discretized finely. Spherical harmonics are also used by Kajiya and Von Herzen [11] to calculate anisotropic scattering. Monte Carlo methods are robust and simple techniques for solving light transport equation. Hanrahan and Krueger model scattering in layered surfaces with linear transport theory and derive explicit formulas for backscattering and transmission [14]. Pharr and Hanrahan describe a mathematical framework [15] to solve the scattering equation in context of a variety of rendering problems and also give a numerical Monte Carlo sampling method. The abovementioned models are both powerful and robust but will suffer from problems coming from standard Monte Carlo method such as slow convergence and noise. To solve such problems, Stam and Fiume show that the widely used diffusion approximation can produce good results for scattering in dense media [16]. Jensen et al. introduce a computationally efficient and analytical diffusion approximation method for multiple scattering [9], which is especially preferable for homogeneous materials that exhibit considerable subsurface light transport ability; however, the model does not appear to be easily extended to volumes with arbitrary optical properties [2].
As a matter of fact, simulations about full light transport are computationally expensive and interactive operation is limited; for example, the change of the illumination or transfer function cannot executed on-line. On the other hand, although there exist analytical approximation methods, most of them are derived based on some assumptions, such as homogeneous optical properties, density, simple lighting, and unrealistic boundary conditions [2]. So such analytical approximation models cannot be used for arbitrary volume rendering or scanned data where optical properties such as absorption and scattering coefficients are hard to be obtained. The disadvantages of huge computation cost and large memory requirement hamper its application in practice. Therefore, some optical empirical models that are relatively easy to implement in medical volume rendering are widely used in medical applications. Those optical empirical models can be realized by graphics hardware-based or CPU-based volume rendering methods. For example, based on the method in document [1], graphics hardware-based optical empirical models for direct volume rendering are developed in documents [2, 17–20]. As to CPU-based method, it is well known that the Shear-Warp method [21] is currently the fastest CPU-based volume rendering method [22, 23], of which the Shear-Warp shell rendering method [24] is proven to be extremely fast. Besides that, some other empirical optimal models are also in common use; for example, the basic lighting absorption and emission model is modified by adding a direct scattering factor to enhance the display effect of inner hidden interfaces.
Shear-warp is considered to be the fastest software algorithm, but it achieves this rendering speed only by sacrificing interpolation between the slices of the volume data. The resample method in traditional Shear-Warp algorithm leads to aliasing and staircase effects, which is in general insufficient for medical purposes. Preintegrated volume rendering [9] provides an efficient way to interpolate in-between slices of the volume data with some loss in rendering performance [23]. Document [23] integrates preintegrated volume rendering into the Shear-Warp algorithm to overcome its drawback. But how to integrate the lighting model into preintegrated volume rendering has not been considered elaborately in document [23].
In this study, an improved volume rendering optical model (IVROM) for medical translucent volume rendering and its implementation using the Shear-Warp and preintegrated Volume Rendering algorithm are proposed in this paper, which can be readily applied on a commodity PC. In the proposed model, the lighting absorption and emission model are employed; moreover, other factors are also considered such as the volumetric shadowing, direct, and indirect emission. Furthermore, the display effect of inner hidden interfaces is enhanced by employing rendering technique of feeling of unreality. Finally, the realization of Improved Translucent Volume Rendering Method (ITVRM) combining the IVROM model, Shear-Warp, and preintegrated volume rendering algorithm is described in detail, in which the aliasing and staircase effects resulting from under-sampling in traditional Shear-Warp algorithm are avoided by the preintegrated volume rendering. By using the proposed method, the 3D medical image can be reconstructed efficiently, and the detailed information of inner structures, as well as the inner hidden interfaces between different layers, can be clearly displayed.
2. The Description of LAEM-BP
Important terms used in the paper are found in Table 1.
Classical lighting absorption and emission model is described as follows: where is the 3D coordinate at position along the view direction. is the attenuation from to along the view direction, .
Usually, the scattering lighting intensity at position along direction is where is the lighting intensity at position along the ray direction . is a bidirectional reflection distribution function (BRDF).
In order to realize surface shading, Blinn-Phong model is used, which depicts the part of the direct scattering of . If there is no obstacle along the ray direction , is equal to the intensity of light source . Therefore, combining the Blinn-Phong surface shading model with (1), the LAEM-BP can be written as where is the value of Blinn-Phong surface shading model which is approximated by the normalized gradient at position ; denotes the distance from the view point to some point along the view direction; is the length from the view point to some point along the view direction; denotes light source intensity; subscript denotes ray direction along .
3. The IVROM Optical Model
The LAEM-BP is more suitable for obvious translucent tissue surface; however, it cannot clearly display the detailed information of the inner structures and the inner hidden interfaces between different mediums, especially when scattering is the dominant factor in the model. Also, it cannot clearly display the region with low gradient yet belong to different tissues. Because a ray has different attenuation values for different tissues, the intershading can be used to compensate the shortcoming existent in LAEM-BP.
In reality, there exists attenuation phenomenon when a ray goes through the tissue. Suppose that the lighting intensity of the indirect ray along direction whose source intensity is at position is . Then
So, (3) can be rewritten as where denotes the distance along ray direction, and subscript denotes ray direction along . In a randomly distributed media space , when a ray goes through the space , there exists extinction phenomenon of absorption and scattering along the transmission direction due to the interaction between the light and medium particles. Therefore, at position should be considered to include the collimation radiation intensity and scattering radiation intensity. Furthermore, the latter part can be divided into two parts: the direct scattering and indirect scattering. Assume that is a bidirectional reflection distribution function (BRDF); BRDF of a particle at position is where is the reflectance ratio of a particle; is phase function, which denotes its direction of scattering.
In fact, in (3), the effect of direct scattering has to be considered. Let denote the phase function of direct scattering; then where is the gradient at position ; in other word, it is the normal vectors of the inner hidden interfaces between different mediums. Considering positive and reverse direction of hidden interfaces, the absolute value of is given.
If the gradient value in a region is high, the Blinn-Phong surface shading model is used to calculate shading effect; reversely, if the gradient value of the region is low, the above model is used to realize shading computation. Therefore, (5) can be rewritten as where non-photorealistic rendering is used to enhance the display effect of inner hidden interfaces. For example, a weight can be given for the direct and indirect scattering. Suppose direct scattering , where is the given weight. Obviously, the further the distance to inner hidden interfaces is, the bigger the surface scattering intensity is, and vice versa.
As a fact, the weight function is essentially a boundary detection function; then let weight value be proportional to the gradient modulus value. Equation (8) can be rewritten as
Next, the effect of indirect scattering is considered. Let indirect scattering be written as where is the weight, and is [25] as follows: where is the flux the total lighting intensity carried by some photons that correspond to the indirect illumination. is the volume of the sphere containing these photons.
Equation (10) is where is the indirect lighting attenuation coefficient, which is dependent on propagation medium. denotes the indirect attenuation along a ray from position to , .
4. ITVRM Integrating the IVROM Model, Shear-Warp, and Preintegrated Volume Rendering
4.1. Drawback of Shear-Warp
Due to the in-slice sampling of Shear-Warp, the interslice sampling rate varies depending on the viewing angle [26]. Consider Figure 1 where the 2D case is shown. At 4 the distance between adjacent sampling points is ; on views down a major diagonal it is . Basic rules of sampling theory tell that it needs at least a distance of to combat aliasing artifacts and to remain within the Nyquist limit of the sampled signal.
4.2. Avoiding Drawback of Shear-Warp by Preintegrated Volume Rendering
Preintegrated volume rendering overcomes the necessity for extremely high sampling rates by splitting the numerical evaluation of the volume rendering integral into two integrations: one for the continuous scalar field and one for transfer functions; thus, the problematic product of Nyquist frequencies as said in Section 4.1 is avoided [23].
Specify colors and extinction coefficients for each scalar value of the volume data by transfer functions , , and [19, 23]. Equation (13) can be then rewritten as
Let denote the attenuation along the ray from position to . Let be the opacity of the distance, thus . Similarly, let be the indirect opacity of the distance, thus .
The volume rendering integral (13) can be approximated by a Riemann sum of equal ray segments of length . Thus, the direct opacity of the th segment along the ray is approximated by where , are the scalar value at the start and the end of the segment along the ray, respectively. Thus, is a function of , , and , if the latter is not constant. Let ; where
In the same way, let ; the indirect opacity of the th segment along the ray is approximated by where
The intensity at position along the view direction consists of the colors value of direct lighting and that of indirect lighting. The colors value of direct lighting of the th segment along the ray are approximated correspondingly: Analogously to , , is also a function of , , and . Let , where
4.3. Implementation of the ITVRM Algorithm
The proposed ITVRM for medical translucent volume rendering can be easily realized on a commodity PC. The implementation of the CPU-based ITVRM algorithm is described as follows.
4.3.1. Setting of the Transfer Function and the Corresponding Look-Up Table
In order to speed up the medical translucent volume rendering algorithm, the look-up table of the indirect opacity and direct opacity is preset when . Similarly, the primary color can also be preset when . In the proposed method, the term primary color is borrowed form OpenGL terminology in order to denote the color before shading, which is similar to document [19].
4.3.2. Presetting of the Weight of Direct Scattering
As a matter of fact, the nearer the distance from inner hidden interfaces to somewhere is, the larger its surface scattering intensity is, and vice verse. So the weight function is essentially a boundary detection function; hence the weight of direct scattering is proportional to the gradient modulus value. Furthermore, if the effect of (11) is considered, then the weight of indirect scattering should be proportional to the voxel scalar value.
4.3.3. Computing Look-Up Table of Preintegration
The preintegration table is computed. That is, and are preset according to (17) and (19), respectively. is preset according to (23):
4.3.4. Composition of the Intermediate Image
The proposed ITVRM integrates the IVROM model, Shear-Warp, and preintegrated volume rendering algorithm; so it is similar to the Shear-Warp volume rendering. The implementation of the Shear-Warp volume rendering algorithm is usually divided into two steps: shear transformation of 3D data set and warp of 2D image. The proposed volume rendering method ITVRM renders slabs between adjacent slices instead of individual slices, which is different from the traditional Shear-Warp algorithm, shown in Figure 2. The proposed IVROM will be applied in the composition of the intermediate image in the step of shear transformation of 3D data set. The view direction of 3D discrete data set is usually set arbitrarily by users; so the transformation from object space to image space is also arbitrary. The main idea of Shear-Warp algorithm is that the 3D discrete data set is first transformed into an intermediate coordinate system; then in intermediate coordinate system, the view direction is selected to be parallel to an axis of the coordinate system. Since the direction of light source is arbitrary, without loss of generality, for the convenience of description, the light source is assumed to be in the same side with the view point, which can be judged by whether the included angle is bigger than 9 or not. If the direction of light source is in the other side, then the light composition sequence should be reversed.
The presented ITVRM in the paper applies the IVROM optical model combining Shear-Warp and preintegrated volume rendering algorithm. Similar to the Shear-Warp algorithm, each slab between image slices is processed and composed into intermediate image sequencely from forward to back. For example, when the th slab between the th slice and th slice is processed, the steps are as follows.
Step 1 . Calculate the direct scattering term; it is given as
where denotes the current red color component of the RGB value of the pixel at a position of the intermediate image without consideration of direct scattering, and denotes that with direct scattering; denotes the current opacity value at a position of the intermediate image without consideration of direct scattering; denotes the accumulative lighting intensity at a position of the intermediate image without consideration of direct scattering; denotes resampling opacity of current processed voxel in the current slab. denotes the R component in current resampling voxel in the current slab. The superscript is for direct scattering. The R component before shading in current resampling voxel in the current slab is
Then is computed by (22) and (25). So,
where denotes the opacity value at the current position of the intermediate image after computing direct scattering term.
When the opacity and colors of the current slab are rasampled and composited, the data of back slice data are used as the front slice of the next slab procession which are stored using the intermediate image-sized buffer plane shown as Figure 3.
Since the composition sequence of the light is from back to forward, so
where denotes the accumulative intensity of direct lighting at current position of the intermediate image with the consideration of direct scattering.
Similarly, the direct scattering value for green and blue color can be calculated in the same way.
Step 2 . Calculate the accumulative lighting intensity.
Since the composition sequence of the light is from back to forward, so
where denotes the current accumulative intensity of indirect lighting at a position of the intermediate image without consideration of indirect scattering. denotes that with indirect scattering. denotes resampling indirect opacity of current voxel. The superscript is for indirect scattering.
Step 3 . Calculate the indirect scattering term:
where denotes the R component of RGB value at current position of the intermediate image without consideration of indirect scattering; denotes that with indirect scattering term; denotes the R component value of current resampling voxel when computing the indirect scattering term.
Considering each slice processed in sequence from forward to back in Shear-Warp algorithm, let in (29) be the average of six pixels which include the pixels at position of indirect lighting memory of next or/current slice and their neighbor four pixels. Here, is the approximation value of in (12). The approximation of is different from traditional approximation [27]. In traditional approximation [27] of physically based light transport equation, incoming light, at any sample of general light transport scenario, scattered from all directions over the unit sphere is considered as shown in Figure 4(a) and (12). Instead, the approximation of only considers light scattered in the forward direction within the cone of directions as shown in Figure 4(b), which is similar to document [2]. Reasons for the approximation of are as follows. () Each slice is processed in sequence from forward to back in this preintegrated Shear-Warp algorithm. () IVROM and its implementation ITVRM are hoped to be readily applied on a commodity PC. () According to [2], the question of whether the missing paths involving lateral movements outside the cone or any backscattering create a barrier to achieving important visual effects is an empirical one. Human viewers are believed not highly sensitive to the details of indirect volume lighting; so there is reason to hope that our approximation is useful. Here the minimum cone is used to approximately compute the effect of indirect scattering term instead of cone [2, 25].
| (a) Traditional approximation, at any sample , scattered from all directions over the unit sphere |
(b) Our approximation, only considering light scattered in the forward direction within the cone of directions
The indirect scattering value of G and B components of the intermediate image can be calculated by the same way.
Step 4 . Read all the data of indirect lighting memory of next slice to be processed into the memory of the current slice; then set all the indirect lighting memory of next slice as value 1.
Step 5 . Process next slice.
4.4. Memory Consumption
The proposed ITVRM integrates the IVROM model, Shear-Warp, and preintegrated volume rendering algorithm; so it is similar to the Shear-Warp volume rendering. In the above implementation processes of the CPU-based ITVRM algorithm, the memory consumption of this ITVRM method is increased a little bit in comparison with that of the traditional Shear Warp algorithm [21]. As to this method, extra memories for the storage of the look-up table memory of the indirect opacity and preintegration, the direct lighting and indirect lighting memory of current slice, and indirect lighting memory of the next future slice are needed. That is, if the volume resolution is with an intensity value of bits (now an intensity value of CT slices is usually 12 bits), and the resolution of the largest intermediate image is , extra memories of are needed. Its increasing quantity is little if compared with the necessary memory for the storage of original mass 3D data set. So it can be concluded that the advantage of shear-warp algorithm over the texture-based is not canceled off by integrating lighting models.
5. Experimental Results
The algorithm is programmed in C language. The hardware environment of the host computer is Pentium4 CPU with frequency 3G, 1024 M RAM, and Graphic Card Geforce 6800 GT.
In order to investigate the effect and feasibility of the proposed optical model applied for translucent volume rendering of 3D medical image, the medical CT slices are reconstructed by both methods of ITVRM and traditional Shear-Warp volume rendering, respectively. The proposed ITVRM integrates the optical model IVROM, and the traditional Shear-Warp volume rendering uses traditional LAEM-BP. They are implemented under the same condition, that is, the same transfer function and view angle. For example, if setting the view angle and , that means the view light is rotated by about axis and about axis. The medical CT slices used are data set of with an intensity value of 16 bits.
Results are separately shown in Figures 5–10. Figure 5 is the result of ITVRM by optical model IVROM; Figure 6 is the result of direct scattering of IVROM; Figure 7 is the result of indirect scattering of IVROM; Figure 8 is the result of Shear-Warp volume rendering by traditional LAEM-BP. From the comparison between Figures 5 and 8, it can be obviously seen that the detailed information of inner structures and the inner hidden interfaces between different mediums can be displayed more clearly in Figure 5; thus the proposed method is more suitable for medical application, such as diagnosis.
Figures 9 and 10 show the 3D results of translucent volume rendering for human brain by the proposed ITVRM using optical model IVROM and the traditional Shear-Warp using LAEM-BP, respectively. Note that the soft tissue of the brain and the gradient norm difference of different tissues, such as cerebral white matter and ectocinerea, have nearly no difference between them. So the detailed information of the soft tissue of brain cannot be depicted clearly in Figure 10 which is reconstructed by traditional model LAEM-BP. Compared with Figure 10, it is obvious that the detailed information of inner structures and the inner hidden interfaces between different mediums can be displayed more clearly in Figure 9.
Figure 11 is the Reconstruction effect comparison between ITVRM and traditional Shear-Warp. Compared with Figure 11(a), it is obvious that the aliasing and staircase effects resulting from under-sampling in Shear-Warp are avoided by using the proposed ITVRM algorithm. Figure 12 shows the 3D results of volume rendering for human head by preintegrated volume rendering without shading and combining the proposed IVROM model. Compared with Figure 12, it is obvious that the detailed information of interfaces of the 3D medical image can be displayed more clearly in Figure 11(b).
(a) The result by traditional Shear-Warp
(b) The result by ITVRM
The above experimental results verify that it is very effective to the proposed ITVRM method applying the presented optical model IVROM, compared with the LAEM-BP model; thus it is more valuable for medical application.
The comparison of runtime for 3D reconstruction between the proposed between ITVRM and traditional Shear-Warp is shown in Table 2. In this test, the condition is kept unchanged except the is rotated by the increment of 30 degrees each step, so it has a total of 12 steps in a cycle, and then the running time is averaged. From Table 2, it can be seen that the running time for the proposed ITVRM using optical model IVROM is just a little longer than that for the traditional Shear-Warp using LAEM-BP model. The reasons can be explained that the calculation for direct and indirect lighting cumulation intensity, as well as the calculation of indirect scattering, needs more time than that of the traditional LAEM-BP model.
Besides, in the processing of setting transfer function, the runtime for translucent volume rendering method by the proposed ITVRM combining IVROM model is also a little longer than the traditional Shear-Warp combining LAEM-BP model. The reason is because that the execution of checking the look-up table for the weight of indirect opacity, direct scattering, and indirect scattering needs more time for the proposed optical model IVROM. From the above analysis, it can be drawn that even the running time for the proposed ITVRM combining IVROM model is longer than that of the traditional Shear-Warp combining LAEM-BP model, but it can be still realized on a commodity PC; moreover, the quality of the reconstructed 3D image of ITVRM combining IVROM is greatly improved.
6. Conclusions
In this paper, an improved volume rendering optical model (IVROM) for translucent volume rendering and its implementation using the Shear-Warp and preintegrated Volume Rendering algorithm are proposed in this paper, which can be readily applied on a commodity PC. In the proposed model, the lighting absorption and emission model are employed; moreover, the factors of volumetric shadowing, direct, and indirect emission are also considered. In addition, the realization of translucent volume rendering method ITVRM combining IVROM model, Shear-Warp, and preintegrated volume rendering algorithm is presented in detail. So the 3D medical image could be reconstructed efficiently, and the detailed information of inner structures, as well as the inner hidden interfaces between different mediums, can be display clearly. Experiment results demonstrate the good performance of the proposed method. Therefore, it is very preferable for practical applications.
Acknowledgments
This work is supported by the Natural Science Foundation of Guangdong province, China (no. 8451064101000631), the Ph.D. Programs Foundation of Ministry of Education of China (no. 200805610018), China Postdoctoral Science Foundation-funded project (no. 20090450866), Cooperation project of Industry, Education and Academy, sponsored by Guangdong province government and Education Department of Chinese government (no. 2009B090300057), Special funds to finance operating expenses for basic scientific research of Central Colleges (South China University of Tech., no. 2009ZM0077), and science and technology key project (no. 2009-Z-108-1) of Panyu District, Guangzhou City.