Interactive high-quality visualization of color volume datasets using GPU-based refinements of segmentation data

Abstract

Data sets containing colored anatomical images of the human body, such as Visible Human or Visible Korean, show realistic internal organ structures. However, imperfect segmentations of these color images, which are typically generated manually or semi-automatically, produces poor-quality rendering results. We propose an interactive high-quality visualization method using GPU-based refinements to aid in the study of anatomical structures. In order to represent the boundaries of a region-of-interest (ROI) smoothly, we apply Gaussian filtering to the opacity values of the color volume. Morphological grayscale erosion operations are performed to reduce the region size, which is expanded by Gaussian filtering. Pseudo-coloring and color blending are also applied to the color volume in order to give more informative rendering results. We implement these operations on GPUs to speed up the refinements. As a result, our method delivered high-quality result images with smooth boundaries and provided considerably faster refinements. The speed of these refinements is sufficient to be used with interactive renderings as the ROI changes, especially compared to CPU-based methods. Moreover, the pseudo-coloring methods used presented anatomical structures clearly.

Keywords

Color volume rendering color human body Gaussian filtering morphology operations pseudo-coloring GPU-based refinement

1 Introduction

In order to visualize medical images, direct volume rendering (DVR) is frequently used. Computerized tomography, magnetic resonance images, and ultrasound datasets are used in medical imaging areas [1]. Since these images do not contain optical information, DVR normally classifies voxels using a transfer function to map voxel values to colors. There is an increasing need to visualize the original color of human organs. Color volume datasets such as Visible Human (VH) [2], Chinese Visible Human (CVH) [3], and Visible Korean (VK) [4] provide anatomical structure information. However, the detection of an object’s boundary is very difficult for color volume data in comparison to gray-scale volume data, which has a scalar value. Hence, it is difficult to render color volume data using a general DVR method [5]. Color values are generally transformed to a monochrome space using various color spaces when the images are segmented [6 –8]. OTF can be applied to the CIE L*U*V color space, transformed from the RGB color space [9]. However, these methods are effective only for specific organs having a similar color distribution.

Segmentation data can be used as additional volume data containing ROI information. With this data, it is possible to render color volume datasets [10, 11]. However, segmentation is a time-consuming task because it is performed manually and these color volume datasets consist of thousands of slices. Furthermore, as illustrated in Fig. 1, the result may not be precise due to the complicated human anatomical structures. Since segmentation data has a discrete volumetric form, it can cause the aliasing artifacts while rendering.

There are two methods that can be used to render anatomical structures from segmentation data: extracting surfaces [12 –17] and reconstructing a volumetric model [11 , 19]. The use of a surface model to visualize the data by exploiting the existing graphics pipeline is straightforward and delivers smooth results and fast renderings. However, it is difficult to observe the interior of objects, a process that is necessary for observing anatomical structures. DVR can be used to represent the interior of objects with a volumetric model. However, the rendering quality is dependent on the resolution of the volume.

We propose an interactive high-quality color volume rendering method with a volumetric model using fast refinements of the segmentation data. The volume data opacity values are filled using the segmentation data and the user-defined indexed opacity table. To remove artifacts on the boundaries caused by discrete segmentation, we apply Gaussian filtering to the opacity values of the volume data. When the slice interval is greater than the pixel spacing, we use a filter having different sized radii along the axes. A morphological grayscale erosion operation is performed to shrink boundaries that are over-expanded due to Gaussian filtering. This process can also be effective in removing unnecessary objects on boundaries, such as resin that contacts the skin. Since these operations can be parallelized, we implement the operations on graphical processing units (GPUs) in order to process the refinements quickly. Our method provides not only realistic renderings using a real color volume, but also object distinction using pseudo-coloring and blending using user-defined color values.

The remainder of this paper is organized as follows: In the next section, we give the problem background. The proposed method is described in Section 3. In Section 4 we present the results and discussion. Section 5 concludes this paper.

2 Background

2.1 Color volume rendering

To perform DVR, ray integration processes are done on each ray by assuming a volume dataset as an optical medium [20]. In addition, the emission and absorption of the volume are evaluated as the rays pass through the volume. The ray integration process is formulated as follows: $I = \int_{0}^{D} \tilde{c} (s (x (λ))) e^{- \int_{0}^{λ} τ (s (x (λ^{'}))) d λ^{'}} d λ$ (1) where I is the integrated value from the observer to a maximum distance D assuming that there is no medium after D, x(λ) is the position along the ray of distance λ from the observer, s(x) is the scalar value at the position x. The color density $\tilde{c}$ (s) and extinction density τ(s) are evaluated via the transfer function.

Color volume rendering is a type of pre-classification method because emission and absorption properties are stored in the volume. Thus, a transfer function cannot be applied. Ray integration of color volume rendering is described as follows: $I = \int_{0}^{D} \tilde{c} (x (λ)) e^{- \int_{0}^{λ} τ (x (λ^{'})) d λ^{'}} d λ$ (2)

This equation is similar to Equation (1) except that the s term is absent, because $\tilde{c}$ and τ are directly determined by x. This integral should be approximated discretely in order to obtain the color and opacity throughout the volume. Dividing the ray into a constant step size, we get an approximation of the color and opacity using a Riemann sum as follows: $\begin{matrix} I \approx \sum_{i = 0}^{n} {\tilde{C}}_{i} \prod_{j = 0}^{i - 1} (1 - α_{j}), \\ \partial_{i} \approx 1 - e^{- τ (x (id)) d}, \\ {\tilde{C}}_{i} \approx \tilde{c} (x (id)) d, \end{matrix}$ (3) where d = D/n.

2.2 Gaussian filtering and morphological grayscale erosion operation

Gaussian filtering [21] is widely used to blur images. Three-dimensional (3D) Gaussian filtering, which is used in this paper, is described by the following equation: $\begin{matrix} f (x_{0}, y_{0}, z_{0}) = \int_{- \infty}^{\infty} \int_{- \infty}^{\infty} \int_{- \infty}^{\infty} g (x, y, z) I (x_{0} - x, y_{0} - y, z_{0} - z) dxdydz, \\ g (x, y, z) = \frac{1}{(2 π)^{3 / 2} σ_{x} σ_{y} σ_{z}} e^{- \frac{x^{2}}{2 σ_{x}^{2}} - \frac{y^{2}}{2 σ_{y}^{2}} - \frac{z^{2}}{2 σ_{z}^{2}}}, \end{matrix}$ (4) where (x₀, y₀, z₀) is the position to be filtered, g(x, y, z) is the Gaussian weight, and σ _x , σ _y , and σ _z are the standard deviations of each axis. I(x, y, z) is the input value at (x, y, z). To be used in practice, this equation is approximated into a discrete form with a finite range. $f (x_{0}, y_{0}, z_{0}) \approx \frac{1}{W (r)} \sum_{- r}^{r} \sum_{- r}^{r} \sum_{- r}^{r} g (x, y, z) I (x_{0} - x, y_{0} - y, z_{0} - z),$ (5) where r is the radius and W(r) is the normalization factor corresponding to r.

Morphological grayscale erosion [21] is formulated as follows: $(f Θ b) (s) = min {f (s + x) - b (x) | (s + x) \in D_{f}, x \in D_{b}},$ (6) which indicates that (fθb) (s) is the minimum value among the values satisfying the restriction that (s + x) must be in the domain of f, and D_f, and x must be in the domain of b, and D_b, respectively. This generally results in a darker image than the original. In our method, erosion operations shrink the ROI that is expanded through Gaussian filtering. The erosion operation uses a smaller radius than the one used in Gaussian filtering. Thus, fine structures remain after the erosions.

3 Description of method

3.1 Boundary refinements

We refine the color volume data opacity value boundaries. Gaussian filtering is applied to the volume data opacity values. Then, morphological grayscale erosion is performed in order to shrink the boundaries that are over-expanded by the Gaussian filtering.

Generally, color human body data contains only RGB values. Opacity values are assigned to the alpha channel of the volume data using the segmentation data and user-defined indexed opacity table [11]. These values are formed discretely because segmentation tasks are performed on a discrete image space. This process leaves some artifacts during rendering. We apply Gaussian filtering to the volume data opacity values to smooth the object boundaries. Our method applies the Gaussian filter with different r and σ for each axis to reduce discontinuities in datasets that have a larger spacing between the slices than the pixel spacing. This is described as follows: $f (x_{0}, y_{0}, z_{0}) \approx \frac{1}{W (r_{x}, r_{y}, r_{z})} \sum_{- r_{z}}^{r_{z}} \sum_{- r_{y}}^{r_{y}} \sum_{- r_{x}}^{r_{x}} g (x, y, z) I (x_{0} - x, y_{0} - y, z_{0} - z),$ (7) where r_x, r_y, and r_z are the radii of the x, y, and z-axes, respectively. W (r_x, r_y, r_z) is the newly calculated normalization factor under the radii r_x, r_y and r_z. We use a greater r_z than r_x or r_y for the case mentioned previously. This is also beneficial in improving filtering speed for the case when small radii are effective for axial slice filtering.

After Gaussian filtering, the ROI of the original object is expanded to an extent that depends on the radius of the filter. This may cause the undesirable visualization of objects adjacent to the desired objects. Furthermore, there are unwanted mixtures, such as skin boundaries facing resin, which must be peeled out. To remedy this problem, we apply morphological grayscale erosion operations on the opacity values for each axial slice after Gaussian filtering. A disk shaped structuring element, meaning that D_b would be circle centered on s in Equation (6), is used in this paper [21]. The radius of the structuring element, which determines the extent of D_b, is dependent on the pixel pitch of the volume data and the radius of the Gaussian filter. Figure 2 shows an example of our method.

3.2 Fast GPU-based refinements

Applying Gaussian filtering to the opacity values of the volume data requires can have significant computational cost, as indicated in Equation (7). This could prevent our method from providing interactive renderings when the ROI changes, because opacity values of the new ROI must be filtered. To achieve interactive filtering as the ROI changes, we implement GPU-based Gaussian filtering. Since Gaussian filtering can be performed in parallel, GPUs can accelerate the process significantly [22]. For further speed improvement, we use separable filters that have time complexity of O(nr) where n and r are the volume data size and the radius of the filter, respectively [23]. In comparison, the complexity of the brute-force method is O(nr³).

We further optimize filtering by exploiting single instruction multiple data (SIMD) operations that GPUs natively support. GPUs can process a vector containing four scalar values, such as an RGBA channel [24]. This speeds up our process because the number of instructions and threads is reduced. Our method filters four grayscale slices simultaneously by packing them into one RGBA slice.

We also minimize the memory overhead by using the manner of cyclic render-to-texture [25]. Figure 3 presents the overall procedure for filtering the opacity values. We use only two additional slices, which have sizes of L_x×L_y and L_x×L_z, where L_x, L_y, and L_z are the sizes of each volume data axis. During the x-direction filtering stage, each filtering result is stored in the upper slice of the current slice (Fig. 3a). The result of the top slice is stored in an L_x×L_y-sized additional slice as an exception. The storing direction is the opposite of the previous stage during the y-direction filtering stage (Fig. 3b). The z-direction filtering stage is performed with the other additional slice, similar to the x-direction stage (Fig. 3c). The filtering is completed after the slice-shift-back step and stores the result into the original volume space (Fig. 3d). Minimizing the memory overhead is advantageous because there can be a memory space limitation on GPUs, especially for a large dataset.

The erosion operation can be performed on the GPU intuitively [26]. We further optimize this method by exploiting SIMD operations.

3.3 Distinction by pseudo-coloring

Although color volume data shows realistic images, non-photorealistic rendering sometimes provides more informative results for the study of anatomical structures. Our method utilizes pseudo-coloring to visualize and distinguish each object using only a volumetric model. Pseudo-colors can be assigned to each voxel of the color volume data using the segmentation data and a user-defined look-up table. Furthermore, our method can blend real color and pseudo-color using a user-defined weight, which can provide a more plausible and useful result than using only real color. Figure 4 demonstrates an example of our pseudo-coloring method.

4 Results and discussion

Our method is implemented in Direct3D 11 and HLSL on an NVIDIA GTX 780 Ti with 3.0 GB of memory. All experiments were performed using an Intel i7 CPU with 16 GB RAM and the GPU described above. Early ray-termination and empty space leaping were applied to optimize rendering speed [27]. The rendered image’s resolution was 800×800 in all cases. The test datasets are shown in Table 1. All datasets were reconstructed from VK data.

We conducted experiments using the torso dataset, which contains many organs, to demonstrate the effectiveness of our method. Figure 5 shows the rendering results with various refinements of the segmentation data. Figure 5a is the result without any refinement. As shown, the rendered surface was rough because the segmentation data is represented as a discrete form. We applied a 3D Gaussian filter with r = 3 and σ= 1 (Fig. 5c). Consequently, the rendered surface was noticeably smoother than the image shown in Fig. 5a. However, some artifacts were observed at the boundaries where multiple objects met each other (yellow circle in Fig. 5e). The rightmost column of Fig. 5 shows the resultant images that were eroded after Gaussian filtering. As shown, the artifacts between the ribs and lungs were removed (yellow circle in Fig. 5f). This shows that our refinements are an effective way to render smooth boundaries.

Figure 6 shows the results using the liver dataset to compare the effectiveness based on r. The disk shaped structuring element was used with a radius of r = 2. Without any refinement, there were a considerable number of artifacts in the rendering result, especially between adjacent slices (Fig. 6a). This is because the space between the slices (1 mm) is greater than the pixel spacing (0.4 mm) and segmentation tasks were performed on the axial slices. After applying Gaussian filtering with r = 4 and σ= 6, the rendering result was further enhanced by smoothing the boundary (Fig. 6b). However, staircase artifacts remained (red rectangle on Fig. 6b). A greater radius (r = 6) delivered a smoother rendering result (Fig. 6c); however, this produced blurring, resulting in a larger green area than that of Fig. 6b (yellow rectangle on Fig. 6c). Our final method, which used a greater r_z than r_x and r_y, provided an adequately smoothed result while reducing the staircase artifacts (Fig. 6d). This shows that a filter with different sized radii can provide superior results compared to those using a filter of the same radius.

Table 2 shows the refinement times of Fig. 6, which are averaged from ten experimental executions. All experiments were performed using the fully optimized GPU-based method. As r increased, the Gaussian filtering time and the erosion time increased. Our method (third and fourth rows of Table 2) required a small amount of additional time for refinement because the radius of erosion was not influenced by r_z.

Figure 7 shows the rendering results using pseudo-coloring and blending with the torso dataset. Gaussian filtering and erosion operations were applied. Figure 7a is the product of using pseudo-coloring without blending with the original color data. As shown, our method generated a rendering result that clearly distinguished each organ. Blending the original colors with pseudo-colors not only provided a distinction between the organs, but also delivered a sensation of real texture. By tuning the blending weights, we can emphasize specific organs (Fig. 7c).

Our pseudo-coloring method can also be used to illustrate intra-brain structures. Figure 8 presents the results using the brain dataset. Using pseudo-coloring, our method rendered each part of brain, distinguished by an anatomist. The frontal lobe, parietal lobe, occipital lobe, and temporal lobe were all observed (Fig. 8b). Our method also was used to visualize tissues inside the brain, such as the thalamus, by rendering the right brain only (Fig. 8d).

Figure 9 shows the results using the head dataset. Without any refinement, the rendering result was too poor to observe (Fig. 9a). After erosion refinements and Gaussian filtering, blue resins near the skin were removed and the skin boundary was smoothed (Fig. 9b). Finally, our method removed skin pigmentation by blending it with a pre-defined skin color (Fig. 9c).

Table 3 shows the refinement times for the brain dataset. The results were averaged from ten tests. As shown, our basic GPU-based method was approximately 46 times faster than the CPU-based method and has O(nr) time complexity. Our GPU-based optimization using SIMD was 86 times faster in comparison with the CPU-based method. As a result, interactive renderings while changes to the ROI are made are feasible.

5 Conclusion

We proposed an interactive high-quality color volume rendering method that employs Gaussian filtering, morphological grayscale erosion, and blending of original color with user-defined colors. Our method applied Gaussian filtering to the opacity values of the color volume data resulting in smooth boundaries. The morphological grayscale erosion operation effectively contracted the ROI that was expanded by the Gaussian filtering. Moreover, erosion removed artifacts caused by multiple objects that were near or facing each other, such as resin on skin. Pseudo-coloring helped us deliver a clear distinction between each organ; blending real color with user-defined colors gave the images the appearance of real texture. Our GPU-based refinement was approximately 86 times faster than the CPU-based method. This allowed interactive rendering whenever the ROI changed. By implementing changes a cyclic manner, our approach minimized the memory overhead due to refinement.

Footnotes

Acknowledgments

This research was supported by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIP) (No. 2015R1A2A2A01008248). This research was supported by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education (2015R1D1A4A01020277).

References

Levoy

Display of surfaces from volume data, IEEE Computer Graphics and Applications 83 (1988), 29–37.

Ackerman

M.J.

The visible human project, Proceedings of the IEEE 863 (1998), 504–11.

Zhang

S.-X.

, Heng

P.-A.

, Liu

Z.-J.

, Tan

L.-W.

, Qiu

M.-G.

, Li

Q.-Y.

, Liao

R.-X.

, Li

, Cui

G.-Y.

, Guo

Y.-L.

, Yang

X.-P.

, Liu

G.-J.

, Shan

J.-L.

, Liu

J.-J.

, Zhang

W.-G.

, Chen

X.-H.

, Chen

J.-H.

, Wang

, Chen

, Lu

, You

, Pang

X.-L.

, Xiao

, Xie

Y.-M.

and Cheng

J.C.-Y.

The chinese visible human (CVH) datasets incorporate technical and imaging advances on earlier digital humans, Journal of Anatomy 2043 (2004), 165–73.

Park

J.S.

, Chung

M.S.

, Hwang

S.B.

, Lee

Y.S.

, Har

D.-H.

and Park

H.S.

Visible Korean Human: Improved serially sectioned images of the entire body, IEEE Transactions on Medical Imaging 243 (2005), 352–60.

Liu

, Chen

J.X.

and Yang

Real-time photorealistic virtual human anatomy, Computing in Science and Engineering 102 (2008), 41–47.

Pal

N.R.

and Pal

S.K.

A review on image segmentation techniques, Pattern Recognition 269 (1993), 1277–94.

Spirkovska

A summary of image segmentation techniques, National Aeronautics and Space Administration, Ames Research Center, 1993.

Gauch

J.M.

, Hsia

C.W.

Applications in optical science and engineering, International Society for Optics and Photonics, Comparison of three-color image segmentation algorithms in four color spaces (1992), 1168–81.

Ebert

D.S.

, Morris

C.J.

, Rheingans

, Yoo

T.S.

Designing effective transfer functions for volume rendering from photographic volumes, IEEE Transactions on Visualization and Computer Graphics 82 (2002), 183–97.

10.

Park

J.S.

, Chung

M.S.

, Hwang

S.B.

, Lee

Y.S.

, Har

D.-H.

Technical report on semiautomatic segmentation using the adobe photoshop, Journal of Digital Imaging 184 (2005), 333–43.

11.

Kwon

K.-J.

, Shin

B.-S.

Visualization of segmented color volume data using gpu, in: Advances in Artificial Reality and Tele-Existence, Springer, (2006), 1062–69.

12.

Shin

D.S.

, Park

J.S.

, Park

H.S.

, Hwang

S.B.

, Chung

M.S.

Outlining of the detailed structures in sectioned images from visible korean, Surgical and Radiologic Anatomy 343 (2012), 235–47.

13.

Shin

D.S.

, Chung

M.S.

, Lee

J.W.

, Park

J.S.

, Chung

, Lee

S.-B.

, Lee

S.-H.

Advanced surface reconstruction technique to build detailed surface models of the liver and neighboring structures from the visible korean human, Journal of Korean Medical Science 243 (2009), 375–83.

14.

Shin

D.S.

, Park

J.S.

, Lee

S.-B.

, Lee

S.-H.

, Chung

M.S.

Surface model of the gastrointestinal tract constructed from the visible korean, Clinical Anatomy 225 (2009), 601–09.

15.

Jang

H.G.

, Chung

M.S.

, Shin

D.S.

, Park

S.K.

, Cheon

K.S.

, Park

H.S.

, Park

J.S.

Segmentation and surface reconstruction of the detailed ear structures, identified in sectioned images, The Anatomical Record 2944 (2011), 559–64.

16.

Shin

D.S.

, Chung

M.S.

, Park

J.S.

, Park

H.S.

, Lee

S.-B.

, Lee

S.-H.

, Choi

H.-N.

, Riemer

, Handels

, Lee

J.E.

et al., Three-dimensional surface models of detailed lumbosacral structures reconstructed from the visible korean, Annals of Anatomy-Anatomischer Anzeiger 1931 (2011), 64–70.

17.

Shin

D.S.

, Park

J.S.

, Shin

B.-S.

, Chung

M.S.

Surface models of the male urogenital organs built from the visible korean using popular software, Anatomy and cell biology 442 (2011), 151–59.

18.

Uhl

J.-F.

, Park

J.S.

, Chung

M.S.

, Delmas

Three-dimensional reconstruction of urogenital tract from visible korean human, The Anatomical Record Part A: Discoveries in Molecular, Cellular, and Evolutionary Biology 2888 (2006), 893–99.

19.

Huang

Y.-X.

, Jin

L.-Z.

, Lowe

J.A.

, Wang

X.-Y.

, Xu

H.-Z.

, Teng

Y.-J.

, Zhang

H.-Z.

, Chi

Y.-L.

Three dimensional reconstruction of the superior mediastinum from chinese visible human female, Surgical and Radiologic Anatomy 327 (2010), 693–98.

20.

Max

Optical models for direct volume rendering, IEEE Transactions on Visualization and Computer Graphics 12 (1995), 99–108.

21.

Gonzalez

R.C.

, Woods

R.E.

Digital image processing, 2002.

22.

Podlozhnyuk

Image convolution with cuda, NVIDIA Corporation white paper, June 2007 (3).

23.

Gaussian blur (2015). URL http://en:wikipedia:org/wiki/Gaussian blur

24.

Mueller

, Xu

, Neophytou

Why do gpus work so well for acceleration of CT?, SPIE Electronic Imaging’07, (2007).

25.

Yan

, Tian

, Zhu

, Dai

, Qin

Fast cone-beam CT image reconstruction using gpu hardware, Journal of X-ray Science and Technology 164 (2008), 225–34.

26.

Rane

M.A.

Fast morphological image processing on gpu using CUDA, 2013.

27.

Levoy

Efficient ray tracing of volume data, ACM Transactions on Graphics 93 (1990), 245–61.