Volume Visualization: A Technical Overview with a Focus on Medical Applications

Multiplanar Reformation

MPR is an image processing technique, which extracts two-dimensional (2D) slices from a 3D volume using arbitrarily positioned orthogonal or oblique planes.⁸ Although it is still a 2D method, it has the advantages of ease of use, high speed, and no information loss. The observer can display a structure of interest in any desired plane within the data set, and four-dimensional (4D) MPR can be performed in real time using graphics hardware.⁹ In addition, to accurately visualize tubular structures such as blood vessels, curved MPR may be employed to sample a given artery along a predefined curved anatomic plane.¹⁰ Curved MPR is an important visualization method for patients with bypass grafts and tortuous coronary arteries.¹¹ Figure 1 illustrates a typical implementation of plane-based MPR showing three cross planes with arbitrary orientation, displayed within a volume-rendered heart. Figure 2a describes a curved MPR, defined by the green lines within a volume-rendered cardiac image.¹² MPR can also serve as a useful adjunct to the axial images for detecting an inflamed appendix¹³ (Fig. 2b). In addition, curved MPR was combined with DVR to enhance the visualization of tubular structures such as trachea and colon.¹⁴

A major disadvantage of MPR is that it is a 2D display technique and orientation dependent. Therefore, when MPR is used to display coronary arteries, it may introduce false-positive or false-negative stenoses.¹⁵ Furthermore, although algorithms have been developed for automatically defining vessel centerlines, the curved MPR is still operator dependent, which limits the diagnostic accuracy.¹⁶

Surface Rendering

SR¹⁷ is a common method of displaying 3D images. ISR can be considered as surface modeling, while DSR is a special case of DVR. ISR requires that the surfaces of relevant structure boundaries within the volume be identified a priori by segmentation. In recently years, the performance of surface reconstruction and display was improved by employing graphics processing unit (GPU).^18,19

For DSR, the surfaces are rendered directly from the volume without intermediate geometric representations, setting thresholds or using object labels to define a range of voxel intensities to be viewed. Recently, there have been many publications on isosurface extraction and rendering.

Hadwiger et al.²⁰ proposed a real-time pipeline for implicit surface extraction and display, which was based on a two-level hierarchical data representation. Hirano et al.²¹ illustrated a semi-transparent multiple isosurface algorithm by employing the span space to rapidly search the intersections of cells with isosurfaces, while Petrik and Skala²² transformed the data into an alternative one-dimensional (1D) space to speed up isosurface extraction. A fully parallel isosurface extraction algorithm was presented by Zhang et al.²³ In addition, a high degree of realism can be achieved with lighting models that simulate realistic viewing conditions. The first image of Figure 3 illustrates an implementation of ISR.

SR is often applied to contrast-enhanced CT data for displaying skeletal and vascular structures, and is also usually used in describing vascular disease and dislocations.²⁴ In the process of detecting acute ischemic stroke, Schellinger and colleagues²⁵ combined maximum intensity projection (MIP) with shaded SR to visualize proximal middle cerebral artery mainstem occlusion (Fig. 4a). SR has also been exploited to describe polyps within the 3D endoluminal CT colonographic images²⁶ (Fig. 4b).

It is difficult to justify the accuracy and reliability of the images generated with shaded SR, i.e., the shiny surfaces might be misleading, causing the relation between image data and brightness in the resultant image to become more complex, a property which could affect the diagnosis. One example illustrating this phenomenon might be that vascular structures displayed with shaded SR may be easily interpreted as a severe stenosis just because the threshold slightly deviates from an optimal value.^1,5

Volume Rendering

DVR displays the entire 3D dataset as a 2D image, without computing any intermediate geometry representations.^27–29 This algorithm can be further divided into image-space DVR, such as software- ^30,31 and GPU-based raycasting,³² and object-space DVR, such as splatting,^33,34 shell rendering,³⁵ texture mapping (TM),³⁶ and cell projection.³⁷ Shear-warp³⁸ can be considered as a combination of these two categories. In addition, MIP,³⁹ minimum intensity projection (MinIP), and X-ray projection⁴⁰ are also widely used methods for displaying 3D medical images. This paper focuses its attentions on DVR, and the datasets discussed here are assumed to be represented on cubic and uniform rectilinear grids, such as are provided by standard 3D medical imaging modalities. Figure 5 describes the DVR pipeline with corresponding computations. Figure 6 shows an example of our DVR results applied to a CT cardiac volume.

DVR is an efficient method for evaluating tissue structures, showing the vascular distribution and anatomic abnormality.⁴¹ This technique frees the user from making a decision as to whether or not surfaces are presented. When compared with SR, the main advantage of DVR is that interior information is retained, and so provides more information about the spatial relationships of different structures.²⁷ In clinical applications, sensitivity and specificity are the major diagnostic criteria that must be considered. Even though DVR generally has high sensitivity and specificity for diagnosis,⁴² it is computationally intensive, so interactive performance is not always optimum. Another disadvantage is that it may be difficult to interpret the “cloudy” interiors. In addition, for analyzing large amounts of high-resolution medical data, DVR alone is probably not an optimal solution, since so many individual voxels are aggregated in a pixel on the screen during the global rendering process. The use of DVR alone may increase the chance that important information is hidden and overlooked. Therefore, a great deal of effect is currently being executed on merging DVR with newly developed computer-aided detection techniques,⁴³ providing more intelligent visualization for early detection, diagnosis, and follow-up. Moreover, even though DVR is a powerful tool to show entire anatomic structures (for more detailed observation of specific lesions), slab imaging, where thick slices are rendered with DVR or MIP, has generally been used in clinical diagnosis.⁴⁴ In addition, the combination of cross-sectional MPR and DVR can significantly increase the interpretation rate of anatomical structures.⁴⁵

DVR has had a wide range of clinical applications. For example, to perform renal donor evaluation, Fishman et al.⁴⁶ used DVR and MIP to display the renal vein of CT angiography data, as shown in Figure 7; for the DVR-generated image, the left gonadal vein (large arrow) is well defined (a), while the locations of the renal vein and gonadal vein (arrow) are inaccurately depicted in the MIP generated image (b). Gill et al.⁴⁷ used DVR and MinIP to show the central airway and vascular structures of a 60-year-old male who underwent double lung transplantation for idiopathic pulmonary fibrosis, evaluating posttreatment in a noninvasive manner (Fig. 7c, d).

Optical Models

The volume-rendering integral is based on an optical model developed by Blinn,⁴⁸ which describes a statistical simulation of light passing through, and being reflected by, clouds of similar small particles. Such optical models may be based on emission or absorption individually, or both, depending on the applications.⁴⁹ To reduce the computational cost, Blinn assumed that the volume is in a low albedo environment. In this case, a light emission–absorption model is an optimal balance between realism and computational complexity.²⁹ Equation 1 demonstrates this procedure.

1

The solution to this equation is given below, showing the intensity of each pixel.

2

Here describes the volume transparency. The first term I₀ illustrates light coming from the background, and D is the extent of the ray over which light is emitted. The last term demonstrates the behavior of the volume emitting and absorbing incoming light. The source term indicates the color change, and the extinction coefficient τ(λ) defines the occlusion of light per unit length due to light scattering or extinction.

Color and Opacity Mapping

Composition

During the DVR process, a number of composition schemes are commonly employed, including X-ray projection, MIP, MinIP, and alpha blending. In X-ray projection,⁴⁰ the interpolated samples are simply summed, giving rise to an image typical of those obtained in projective diagnostic imaging. It has the advantage of clearly representing vessels or bones at the cost of losing the depth information (Fig. 9a). MIP^39,82 and MinIP are widely used techniques in 3D CT and MR angiography. The image results from the voxels having the MIP or MinIP along the viewing rays traversing through the object. MIP and MinIP images can be generated rapidly and can clearly display vessels, tumor, or bones,²⁴ and the image generation has been accelerated by graphics hardware.⁸³ These two techniques still have some inherent disadvantages because only a small fraction of the data is used, and static MIPs or MinIPs are not particularly useful due to lack of depth information. However, MIP and MinIP can be interactively explored or used as animation sequences, where the limited-depth-perception drawback is less obvious. Local MIP, MinP, or closest vessel projection (CVP) are often used in slab imaging, i.e., thin-slab MIP or MinP, for vascular structure diagnosis.^44,84 For vascular diagnosis, CVP is superior to MIP; however, the user needs to set an appropriate threshold for the local maximum, which is determined by a specific dataset, making the application of CVP more difficult than MIP. Figure 9b shows the result of MIP.

Alpha blending^30,49 is a popular optical blending technique, often implemented by using the Riemann sum to discretize the continuous function Eq. 2, resulting front-to-back and back-to-front alpha blending Eqs. 3 and 4 respectively, depending on the compositing order.

3

4

where now the pre-multiplied α_iC_i is the total light emitted from a point i, which is also called the opacity associated color.⁸⁵T_j(0 ≤ T_j ≤ 1) is the transparency of a point, and α_j is its opacity, which is defined as . Figure 9c describes the DVR result using alpha blending without shading, while (d) shows the result with shading.

Volume Illumination and Illustration

In volume illumination, the normal at every sampling point is calculated by interpolation using the intensity changes across that voxel. These approximated voxel normals are then used in a Phong or Blinn–Phong model for shading computations, with the results being employed in the DVR composition. The shading computations may be accelerated by graphics hardware.⁸⁶ In addition, the shading model can be used with volumetric shadows to capture chromatic attenuation for simulating translucent rendering,⁸⁷ and can be combined with clipped volumes to increase the visual cues.⁸⁸ Figure 10 illustrates a DVR of CT cardiac data sets with and without illumination, demonstrating that images with shading are visually more pleasing. However, clinicians sometimes deliberately disable illumination since it may complicate certain diagnostic tasks.⁸⁴

This technique is helpful to provide insight into medical data, aiding cognition during data exploration. For example, to illustratively define a mapping from volumetric attributes to visual styles, a fuzzy-logic rule-based semantic layer was introduced by Rautek et al.⁶⁸ who extended their work⁸⁹ by introducing an interaction-dependent semantic model. To achieve realistic anatomical volume rendering and illustration, Silverstein et al.⁹⁰ presented a color-mapping algorithm, which combines natural colorization with the perceptual discriminating capacity of grayscale. For the purpose of enhancing anatomical structure explanation, Wang et al.⁹¹ introduced a knowledge-based system to encode a set of illustrative models into the volume-rendered images. Rezk-Salama et al.⁹² generated light field representations from volumetric data, so users could illustratively explore resultant images at run time. Jainek et al.⁹³ presented a hybrid rendering method, combining an illustrative mesh of brain structures with DVR, resulting in a focus-and-context data representation. In addition, to effectively emphasize data in the focus-and-context style, Svakhine et al.⁹⁴ illustrated an illustration-guided DVR algorithm, enhancing anatomical structures and their spatial relationships.

Illumination is an important aspect of volume illustration. Rheingans and Ebert⁹⁵ proposed an illumination approach similar to volume shading, using non-photorealistic rendering to enhance physics-based DVR. Dong and Clapworthy⁹⁶ extended the texture synthesis approach from surface to non-photorealistic DVR, and this volume illumination model was used in conjunction with a TF by Bruckner et al.⁹⁷ to preserve contextual information during the rendering process of interior structures. To enhance material discrimination, Noordmans and colleagues⁹⁸ extended the color- and opacity-based illumination model to wavelength-dependent spectral DVR, and the performance of this spectral DVR was optimized with a linear spectral factor model by Bergner et al.⁹⁹ Recently, for the purpose of generating specified colors or spectral properties in surface or volume rendering, this group used constructed spectra to enhance a palette of given lights and material reflectances.¹⁰⁰ In addition, Strengert et al.¹⁰¹ subsequently employed this spectral model to GPU-based raycasting.

To increase depth perception, Kersten et al.¹⁰² used a purely absorptive lighting model in DVR to enhance the display of translucent volume, while Bruckner and Gröller⁷² employed a style TF to add data-dependent image-based shading to volume-rendered images, and Joshi et al.¹⁰³ exploited tone shading with distance color blending in vascular structure volume visualization. Furthermore, to facilitate tissue illumination, Ropinski and colleagues¹⁰⁴ added a dynamic ambient occlusion simulation model to DVR, and Mora et al.¹⁰⁵ introduced a rendering pipeline for isotropically emissive visualization to facilitate intuitive medical data exploration. To intuitively explore hidden structures and analyze their spatial relations in volumes, Chan et al.¹⁰⁶ introduced a relation-aware visualization pipeline, in which inner structural relationships were defined by a region connection calculus and were represented by a relation graph interface. In addition, Rautek et al.¹⁰⁷ gave a comprehensive review on illustrative visualization and envisioned their potential medical applications. Such illustrative techniques are sometimes integrated in commercial systems; however, there have been no systematic investigations regarding whether any of these techniques actually provides a clinical benefit.¹

Algorithm Overview

Raycasting²⁷ is a popular technique used to display a 3D dataset in two dimensions, whose basic idea is to cast a ray from each pixel in the viewing plane into the volume, sampling the ray with a predetermined step in a front-to-back or back-to-front order by trilinear interpolation. A TF is used to map the scalar value to color and opacity. Finally, the acquired optical values at the sampling points along the casting ray are composited using Eqs. 3 or 4 to approximately compute the DVR integral Eq. 2, obtaining the corresponding final pixel color on the output image. Figure 11 illustrates the pipeline, and Figure 12 shows an implementation on four medical images.

Applications

This technique was a significant 3D display method in some applications during the time when graphics hardware-accelerated DVR is not commonly available. For example, Sakas et al.¹⁰⁸ used raycasting to display 3D ultrasound data of a fetus, and Höhne³¹ and Tiede et al.¹⁰⁹ employed this algorithm to display and investigate the structures of human anatomy.

Improvements

Levoy¹¹⁰ presented an adaptive refinement approach for empty space skipping, and a similar idea was used in a hierarchical spatial enumeration technique,¹¹¹ which first created a binary pyramid, which was used to build a min–max octree to encode the empty and non-empty space, and this min–max octree structure was exploited with a range tree by Lee and Ihm¹¹² to skip coherent space in volume. A similar approach was used by Grimm et al.¹¹³ who proposed an algorithm based on a bricked volume layout.

Algorithm Overview

This technique was developed to improve the calculation speed of DVR at the expense of lower accuracy.^33,34 Splatting calculates the influence of each voxel in the volume upon multiple pixels in the output image. This algorithm represents the volume as an array of overlapping basis functions called reconstruction kernels, which are commonly rotationally symmetric Gaussian functions with amplitudes scaled by the voxel values. These kernels are classified as colors and opacities by a TF, and are projected onto the output image plane in a depth-sorted order. Each projected kernel leaves a footprint (or splat) on the screen, and all of these splats are composed in back-to-front order to yield the final image. This process is described in Figure 13, and Figure 14 presents examples of images rendered with the splatting algorithm.

Applications

Vega-Higuera et al.¹¹⁶ exploited texture-accelerated splatting to visualize the neurovascular structures surrounded by osseous tissue in CT angiography data. Splatting was also employed to rapidly generate digitally rendered radiographs in the iterative registration of medical images.¹¹⁷ Audigier and colleagues¹¹⁸ used splatting with raycasting to guide the interactive 3D medical image segmentation, providing the users with visual feedback at the end of each iterative segmentation step.

Improvements

Since basic splatting algorithms³³ suffer from “color bleeding” artifacts, Westover³⁴ employed an axis-aligned sheet buffer to solve this problem. However, this technique needs to maintain three stacks of sheets and introduces “popping” artifacts, which refers to the abrupt change in image appearance as the view switches from one of the three image stacks to another. To address this issue, Mueller et al.¹¹⁹ aligned the sheet buffers parallel to the image plane, and they later accelerated this image aligned splatting algorithm with GPUs.¹²⁰ Recently, Subr et al.¹²¹ introduced a self-sufficient data structure called FreeVoxels to eliminate the color bleeding errors. In addition, the aliasing problem in splatting was addressed by Zwicker et al.¹²² with the use of an elliptical weighted average filter. To improve splatting performance, Ihm and Lee¹¹⁴ employed an octree to encode the object- and image-space coherence, where only voxels of interest were visited and splatted, while Xue et al.¹²³ increased the speed of splatting through texture mapping.

Algorithm Overview

Shear-warp is a hybrid algorithm that attempts to combine the advantages of the image- or object-order based volume-rendering methods,³⁸ and is based on factorizing the viewing transformation into 3D shear and 2D warp. The shear matrix transforms all the viewing rays so that they are parallel to the principal viewing axis in “sheared object” space, allowing the volume and image to be traversed simultaneously. Figure 15 illustrates the transformation from the original object space to sheared object space for parallel projection.

We can mathematically express the factorization of the viewing matrix M_view, which comprises a permutation matrix P and a view matrix M, such that , and the view matrix M is factorized into a shear S and a warp component . Figure 16 shows the rendering examples.

Applications

Perspective shear-warp was implemented by Schulze et al.¹²⁴ in a virtual reality environment, while Anagnostou et al.¹²⁵ used this algorithm to display time-varying data. Furthermore, Wan et al.¹²⁶ used it to apply 3D warping in an interactive stereoscopic rendering for virtual navigation.

Improvements

This algorithm was thoroughly analyzed by Sweeney and Mueller,¹²⁷ and they proposed several measures to alleviate the short-comings of this algorithm, such as blurry images on zoom and staircase artifacts on flat surfaces. The VolumePro graphics processing board is a hardware implementation of the raycasting using shear-warp factorization for real-time visualization.¹²⁸ One limitation of VolumePro is that it can only provide orthographic projections, which renders it unsuitable for applications such as virtual endoscopy. Therefore, Lim and Shin¹²⁹ employed the subvolume feature in VolumePro to approximate perspective projection, only rendering the subvolumes within the viewing frustum. To alleviate blur in the zoomed image, Schulze et al.¹³⁰ applied a pre-integrated classification to shear-warp. Later, Kye and Oh¹³¹ used a super-sampling to improve the image quality, and a min–max block data structure was designed to maintain the rendering performance.

Algorithm Overview

The shell rendering algorithm is a software-based hybrid of surface and volume rendering, which is based on a compact data structure referred to as a shell, i.e., a set of nontransparent voxels near the extracted object boundary with a number of attributes associated with each related voxel for visualization.³⁵

The shell data structure can store the entire 3D scene or only the hard (binary) boundary. For the hard boundary, the shell is crisp and only contains the voxels on the object surface, and shell rendering degenerates to SR. For a fuzzy boundary, the shell includes some voxels in the vicinity of the extracted surface, and shell rendering is identified as DVR. Figure 17 shows examples of shell SR and DVR, and the detailed descriptions of this algorithm can be found in Udupa and Odhner.³⁵

Applications

Lei et al.¹³⁴ employed this algorithm to render the segmented structures of vessels and arteries of contrast-enhanced magnetic resonance angiography image. However, the explicit surface extraction creates errors. To address the problem, Bullitt et al.¹³⁵ selectively dilated the segmented object boundaries along all axes, and visualized the extracted fuzzy shell with raycasting.

Improvements

To accelerate the shell rendering speed, Falcão and colleagues¹³⁶ added the shear-warp factorization to the shell data structure, and Botha and Post¹³⁷ used a splat-like elliptical Gaussian to compute the voxel contribution energy to the rendered image. Later, Grevera et al.¹³⁸ extended the point-based shell element to a new T-shell element for isosurface rendering.

Algorithm Overview

In 2D texture mapping (2DTM), the volume is decomposed into three stacks of perpendicularly object-aligned polygons. For the current viewing direction, the stack whose slicing direction (normal) is within ±45° of the current viewing direction is chosen for rendering. During rasterization, each of the polygonal slices is textured with the image information obtained from the volume via bilinear interpolation. Finally, the textured slices are alpha-blended in a back-to-front order to produce the final image. Figure 18 describes the working pipeline, while Figure 19 illustrates rendering results with this algorithm.

Three-dimensional texture mapping (3DTM) uploads the volume to the graphics memory as a single 3D texture, and a set of polygons perpendicular to the viewing direction is placed within the volume and textured with the image information by trilinear interpolation. Compared with 2DTM, there are no orientation limitations for the decomposed polygon slices, making the texture access more flexible. The compositing process is similar to 2DTM, in that the set of textured polygon planes are alpha-blended in a back-to-front order. Figure 20 illustrates the algorithm pipeline, and Figure 21 demonstrates an implementation of 3DTM.

Applications

Applications of TM-based medical volume rendering, as reported in the literature, can be divided into three main applications. The first is multimodal medical image rendering and tissue separation. Sato et al.¹⁴¹ designed a multidimensional TF to identify tissue structures in multimodal medical images generated with 3DTM. A two-level rendering technique was integrated with 3DTM by Hauser et al.,¹⁴² allowing different rendering techniques to be selected for different tissues. Later, Hadwiger et al.¹⁴³ improved this two-level DVR algorithm with respect to image quality and performance, minimizing the computational cost of each pass, adding dynamic filtering, and using depth and stencil buffering to achieve correct compositing of objects created with different rendering techniques. Sereda et al.¹⁴⁴ used the LH histogram in a multidimensional TF domain for volume boundary enhancement and tissue detection. In an attempt to improve arterial stenosis diagnosis, Ljung et al.¹⁴⁵ applied a probabilistic TF model to 3DTM, replacing the optical lookup table with a tissue material selection table. Xie et al.¹⁴⁶ used Boolean operations in a 3DTM-based DVR for interactive multi-object medical volume clipping and manipulation, and Kim and colleagues¹⁴⁷ used a dual TF-based algorithm to render fused images of PET and CT in real time. Abellán and Tost⁶⁴ exploited 3DTM for multimodal medical data rendering, in which various fusion techniques were integrated into a common hardware-accelerated DVR framework, which were employed to display registered dual-modal brain datasets.

The second application area is the display of specific tissues and organs for diagnosis and therapy. Holmes et al.¹⁴⁸ used 3DTM to achieve a real-time display of trans-urethral ultrasound for prostate treatment. This technique was also used by Etlik et al.¹⁴⁹ to find bone fractures, and also by Wenger et al.¹⁵⁰ to visualize diffusion tensor MRI data in the brain. Wang et al.¹⁵¹ exploited 3DTM to display 3D ultrasound for cardiac ablation guidance. To provide visual information relating to structures located beneath the skin surface, Sharp et al.¹⁵² employed this technique to visualize the light diffusion in 3D inhomogeneous tissue. This visualization approach was also employed by Lopera et al.¹⁵³ to view a stenosis in arbitrary image planes, and the authors used the results to demonstrate the subsequent arterial bypass. For the purpose of showing early opacification of hepatic veins, Siddiki and colleagues¹⁵⁴ exploited 3DTM to display the arterial phase images. Aspin et al.¹⁵⁵ combined this rendering method with data enhancement techniques to increase the users’ cognitive ability in diagnosing soft tissue pathologies in volume-rendered MR images.

The third area of application is in dynamic imaging and deformation. Levin et al.¹⁵⁶ developed a TM-based platform that was used to interactively render multimodality 4D cardiac datasets. In addition, Yuan et al.¹⁵⁷ designed a TF for nonlinear mapping density values in the dynamic range medical volumes, while Correa and colleagues¹⁵⁸ presented a 3DTM algorithm that sampled the deformed space instead of the regular-grid points, simulating volume deformations caused by clinical manipulations, such as cuts and dissections.

Improvements

Introduction to GPU Architecture

GPUs incorporate multi-pipelined stream processors and high bandwidth memory access, and include their own DRAM (device memory) optimized to process 2D and 3D geometries in parallel.¹⁷⁶ Data are transferred between the CPU main memory and the GPU graphics memory via a dedicated Accelerated Graphics Port or a Peripheral Component Interconnect Express (PCI-Express) bus. Modern GPUs are programmable and employ a SWAR (SIMD, i.e., Single Instruction Stream Multiple Data Stream, Within A Register) execution model. Figure 22 illustrates the architecture and working pipeline of a typical GPU.

As illustrated in Figure 22, the geometry is specified as a set of vertices and a set of facets. The vertex stream is processed by the vertex units and then assembled into primitives, which are grouped together over triangles and are scan-converted by the rasterizer. The outcome of this operation is a large stream of elements (fragments), which are further processed by the fragment units. These fragments are written into a frame buffer and then copied into the video memory to produce an image on the screen. Recent GPUs include a non-visible rendering buffer, called the pixel buffer (P-Buffer), which can be employed if higher precision is required in the intermediate calculations. The nature of the data on which GPUs operate makes them suitable for a stream processing architecture.¹⁷⁷

Recently, NVidia and AMD announced their new software development platforms, i.e., Compute Unified Device Architecture (CUDA)¹⁷⁸ and Close To Metal (CTM),¹⁷⁹ respectively, which are designed to directly access the native instruction set and memory in parallel on GeForce (G8X and later generations) and Radeon (R580 and later generations) graphics cards, respectively. CUDA is a C language development environment, while CTM is a virtual machine running assembler code. The goal of these two new architectures is to provide developers with a streaming-like programming model for compute-intensive parallel computations, allowing low-level, deterministic, and repeatable access to hardware. However, these new generalizations sacrifice some of the efficiencies and graphic specified features.¹⁸⁰ In addition, the bus bandwidth and latency between CPU and GPU still remains a bottleneck for the further performance enhancement.

Algorithm Overview

A general strategy of exploiting the programmable GPU for raycasting was proposed by Krüger and Westermann³² and Purcell et al.,¹¹³ who described a stream model that exploits the intrinsic parallelism and efficient communication on modern GPUs. The main core of this algorithm is to generate a single ray per screen pixel and to trace this ray through the volume on a GPU fragment shader. Assuming that the origin point is O and the ray r(t) direction is v, then r(t) is defined as O + v × t. The algorithm is described by the following items:

• Store the dataset as a 3D texture, whose bounding box is a colored unit cube. The texture coordinates at the eight corners are encoded in the RGB color channels.

• Render the front face of the color cube and store the result into an intermediate texture (Fig. 23a).

• Render the back face (Fig. 23b), subtract the color values of the back face from the front face, normalize the result, and store it in a 2D texture, which represents the direction vectors of the casting ray for each pixel. The origin of the ray is the pixel of the bounding cube front face.

• Cast a ray for each pixel into the 3D texture volume, and uniformly sample the volume along the ray. This is all done in a single fragment program using the GPU Shader Model 3.0 loop capabilities.

• Blend the TF-mapped color and opacity at each sampling point along the ray in a front-to-back order. When the ray leaves the volume or the accumulated opacity exceeds the preset threshold, terminate the ray.

Performance Comparisons

TM-based DVR and GPU-based raycasting are all graphics hardware-accelerated DVR algorithms. The TM uses fixed graphics pipelines, while the raycasting encodes the algorithm on the GPU programmable shaders. For DVR based on TM, the shading calculation is expensive since it needs to apply the computationally expensive shading model to every single fragment, even if some of them do not contribute to the final image. In addition, it is not easy to include some of the traditional acceleration techniques. Finally, for perspective projection, a uniform sampling step cannot be maintained, resulting in striping artifacts. For GPU-based raycasting, since it is a flexible rendering pipeline, the acceleration techniques such as early ray termination, empty space skipping can easily be included. Also, since the shading computation is fragment-based, the shading model need only be applied to the fragments that contribute to the final image. In addition, the sampling step is adaptive, so for perspective projection, uniform sampling steps can be maintained, resulting in improved image quality. Figure 24 illustrates our implementation of GPU-accelerated raycasting of four 3D medical datasets, from which we note that the image quality for raycasting is superior to the images generated with 2DTM (Fig. 19) or 3DTM (Fig. 21).

Applications

As exemplified by Beyer et al.,¹⁸¹ who used GPU-based raycasting to display fused multimodal brain datasets in real time for preoperative planning of neurosurgical procedures, there is a trend toward using TF-based methods for anatomical feature enhancement, tissue separation, and obtaining complementary diagnostic information from different image modalities. In a GPU-based raycasting framework, Weber et al.¹⁸² segmented a volume into regions in the form of branches of a hierarchical contour tree, and created a unique TF for each of these regions. These contour tree-based segmentations were used to optimize the TF design for highlighting the volume-rendered anatomical structures of thoracic CT data, and Termeer et al.¹⁸³ used this rendering technique to build a visualization software platform for the diagnosis of coronary artery disease. The GPU-based raycasting algorithm was also used by Rieder et al.¹⁸⁴ to display multimodal volume data for neurosurgical tumor treatment. They designed individual TFs for each imaging modality, and proposed a distance-based enhancement technique for functional data display and lesion detection, as well as for path to target exploration.

The second main application is volume illustration, applying semantic information to simplify the complexity in the transfer function design and intuitively revealing anatomical structures. Rezk-Salama and colleagues⁶⁷ designed some high-level semantic models for intuitive TF parameter assignments. Rautek et al.⁶⁸ introduced a fuzzy-logic rule based semantic layer to linguistically define a mapping from volumetric attributes to visual styles, and later extended their work⁸⁹ by introducing an interaction-dependent semantic model, and employing graphics hardware to accelerate the fuzzy-logic arithmetic evaluation during the semantic parameter specification process.

The third main application trend is the real-time display of specific tissues or organs. For example, Ljung¹⁸⁵ used this rendering approach to interactively visualize a multi-detector CT-scanned human cadaver to perform a “virtual autopsy.” To accelerate lung cancer screening, Wang and her colleagues¹⁸⁶ developed a software platform based on GPU-accelerated raycasting for displaying volumetric CT images stereoscopically at interactive rates. Krüger et al.¹⁸⁷ used GPU-based rendering technique with MPR to realistically represent volumetric images in a virtual endoscopy application, directly visualizing the individual patient data without time-consuming segmentation. This hardware-accelerated volume visualization technique was also exploited with a vessel detection algorithm by Joshi et al.¹⁰³ to display vascular structures in an angiographic image, while Mayerich et al.¹⁸⁸ designed an anatomical structure encoding algorithm that is used with a GPU-accelerated raycasting to detect and visualize the microvascular network and its anatomical relationships with the surrounding tissues in real time. Recently, Lee and Shin¹⁸⁹ employed the GPU-based volume-rendering algorithm for fast perspective visualization in a virtual colonoscopy application, while Lim et al.¹⁹⁰ used a GPU-based raycasting to visualize 3D ultrasound datasets, during which points within the native pyramidal ultrasound coordinate system were adaptively converted on the fly to a Cartesian grid.

Finally, a DVR system^191,192 was developed in our group, which is based on a GPU-accelerated raycasting. In this system, a novel dynamic texture binding technique was employed to render 4D cardiac datasets with the heart rate being synchronized with the electrocardiogram signal. Subsequently, techniques were developed for dynamic volume rendering, manipulation, and enhancement.^193,194 Furthermore, a new algorithm was presented for rendering segmented datasets with anatomical feature enhancement.^195,196 When compared with previously published results, instead of using different TF for segmented object separation, the new algorithm “paints” various colors acquired from segmentation to the rendered volume, and delivers a higher image quality with smoother subvolume boundaries. Recently, to facilitate image manipulation and image-guided intervention, a new “depth texture indexing” algorithm was developed¹⁹⁷ to allow virtual surgical tools to be included into the 4D medical image visualization environment.

Improvements

The increasing capabilities of modern GPUs provide great potential to enhance the performance of the GPU-based raycasting algorithms. For example, the intrinsically parallel nature of GPUs was exploited by Roettger et al.⁵⁶ to build a real-time raycaster, in which adaptive pre-integrated voxel classification and volume clipping were included. The original GPU-based solutions were based on multiple rendering passes, casting rays through the volume and storing the intermediate results in a temporary framebuffer, which could be accessed in subsequent rendering passes. However, on current graphics hardware architecture, such a multi-pass-based rendering pipeline suffers from low efficiency. To address this issue, Stegmaier and colleagues built a single-pass raycasting framework and extended this framework to spectral DVR for refractive material selective absorption and dispersion using the Kubelka–Munk theory.¹⁰¹

For single-pass raycasting, Hadwiger et al.²⁰ employed this visualization technique to display implicit isosurfaces from a volume defined on a rectilinear 3D grid. A two-level hierarchical grid representation was exploited to allow object- and image-order empty space skipping to be employed to improve the rendering performance, and to circumvent the memory limitations in visualizing datasets that are too large to fit into GPU texture memory. For the similar task of releasing on-board memory size restrictions, Wang and JaJa¹⁹⁸ recently proposed a multicore processor-based raycasting algorithm for isosurface extraction and rendering. In addition, exploded views were introduced into the raycasting algorithm by Bruckner and Gröller¹⁹⁹ to display inner occluded object of interest. Their approach employed a force-based model to partition the volume into several segments for displaying hidden structural details, and the partition manner was controlled by a number of forces and constraints. Similarly, to reveal the interior hidden spatial structures, opacity peeling was employed by Rezk-Salama and Kolb²⁰⁰ in the rendering pipeline. For the purpose of emphasizing objects of interest within the contextual volume, Viola et al.²⁰¹ presented an importance-driven strategy, in which intelligent viewpoint selection and alternative visual representations of the volume were exploited. Furthermore, the importance-driven approach was also used by Wang et al.²⁰² to visualize time-varying data for the purpose of unveiling and enhancing the dynamic features exhibited by these 4D volumetric data.

Compared with multi-pass rendering,⁵⁶ a single-pass solution has the advantages of efficiency, accuracy, and flexibility. Klein et al.²⁰³ proposed techniques to improve the performance and quality of GPU-based volume raycasting. An empty space skipping algorithm that exploits the spatial coherence between consecutively rendered images was employed to estimate the optimal initial ray sampling point for each image pixel. A similar idea to exploit the empty space skipping and frame-to-frame coherence was later used by Lee and Shin¹⁸⁹ to accelerate GPU-based raycasting; however, instead of using super-sampling, they designed a new algorithm to improve image quality by eliminating the artifacts caused by holes or illegal starting positions. Müller et al.²⁰⁴ presented a sort-last parallel architecture to accelerate single-pass volume raycasting on GPU, wherein they partitioned the texture volume into uniform texture bricks. A similar volumetric texture bricking and octree organization technique was later used by Gobbetti et al.²⁰⁵ to improve the performance of GPU raycasting on massive datasets. Furthermore, to render large out-of-core datasets, Moloney et al.²⁰⁶ presented a volume bricking-based sort-first algorithm to balance parallelized data loading and frame-to-frame coherence. The bricked volume technique was also used by Beyer et al.¹⁸¹ to display multimodal volume using GPU raycasting. Similarly, Ljung¹⁸⁵ presented a multiresolution scheme to render large out-of-core CT human cadaver volumes using GPU raycasting, in which techniques such as load-on-demand data management and TF-based data reduction were employed. In addition, to address the difficulties of rendering large datasets on a GPU, a real-time volume compression and decompression technique based on block transformation coding was proposed by Fout and Ma.²⁰⁷ In the GPU-based raycasting pipeline, to render large volumetric data and improve image quality, Wang et al.²⁰⁸ developed a level-of-detail selection algorithm using an image-based quality metric. In addition, to provide novice audiences access to interactive volume renditions, Rezk-Salama et al.⁹² introduced a technique to generate high-quality light field representation and integrated it into the GPU-based raycasting pipeline.

Discussion

DVR is an efficient technique to explore complex anatomical structures within volumetric medical data, with requirement for real-time volumetric and dynamic image visualization; techniques to address these requirements must keep pace with the development of related technologies, and failure to do so will result in visualization being the rate-limiting step for the development of new treatment modalities. Real-time DVR of clinical datasets needs efficient data structures, algorithms, parallelization, and hardware acceleration. The progress of programmable GPUs has dramatically accelerated the performance and medical applications of volume visualization, and opened a new door for future developments of real-time 3D and 4D rendering techniques.

As the capabilities of programmable GPUs continue to increase, the distinction between algorithms that are appropriate for the GPU and those written for the CPU will disappear. The evolution from software-based raycasting to GPU-based raycasting is a good example. For DVR algorithms such as splatting, shell rendering, and shear-warp, hardware accelerations have been proposed and implemented; however, most of these improvements are based on fixed graphics pipelines and do not take full advantages of the programmable features of GPU. As mentioned previously, even if the TM-based DVR algorithms make use of graphics hardware in the volume-rendering pipeline, they differ considerably from the programmable GPU algorithms, such as raycasting, in that all of the volume-rendering computations and corresponding acceleration techniques are implemented on the GPU fragment shaders. The TM-based DVR algorithms only use fixed graphics pipelines that include relative simple texture operations, such as texture extraction, blending, and interpolation.

Since most DVR algorithms are application driven, and medicine is a major application area, in this paper, we described the general techniques of volumetric data visualization with a focus on medical applications. However, for a specific application, different DVR algorithms may generally be used interchangeably when sufficient image quality and rendering performance are available. For the evaluation and comparison of DVR algorithms, Meißner et al.¹¹⁵ quantitatively analyzed four algorithms: software-based raycasting, splatting, shear-warp, and 3DTM. Recently, Boucheny et al.²¹³ presented three perceptual studies to examine the information transformation effectiveness of DVR algorithms, and unsurprisingly, they showed that a static representation is ambiguous. When tuning the rendering parameters, this ambiguity can be counterbalanced by using depth cues (i.e., obtained through motion parallax or stereoscopy). Furthermore, the quantitative evaluation of the four mainstream DVR algorithms was also reported in a publication of our group.¹⁹⁷

Outlook

Modern GPUs support the Shader Model 4.0 features and floating point precision throughout the graphics pipeline. Currently, most of the advanced features of programmable GPUs are accessible through high-level programming languages, such as OpenGL, HLSL, GLSL, and Cg. A promising research direction is in the design of new image rendering and manipulation algorithms, specially for programmable GPU architectures, taking advantage of cache coherence and parallelism to achieve efficiency. With the newly released GPU programmable APIs CUDA¹⁷⁸ and CTM,¹⁷⁹ more complex algorithms can efficiently be implemented on GPUs, and we believe that these two APIs will more efficiently support graphics applications in the future.

In addition, to further enhance the visualization environment, DVR algorithms should be effectively combined with realistic haptic (i.e., force) feedback and crucial diagnostic and functional information extracted from medical data using new algorithms in the research fields such as machine learning and pattern recognition. Furthermore, innovative techniques of robotics, human–computer interaction, and computational vision should be developed to facilitate the medical image exploration, interpretation, processing, and analysis.

Finally, we note that the current volume visualization algorithms should be integrated into standard graphical and imaging processing frameworks such as the Insight Segmentation and Registration Toolkit ²¹⁴ and the Visualization ToolKit,²¹⁵ enabling them to be readily used in clinical medical imaging applications. We also note that even though there have been many efficient and advanced techniques developed for all parts of the volume-rendering pipeline, only a few simple approaches have been developed for clinical visualization. The translation of these advances into clinical practice nevertheless remains a problem. However, we must ensure that such challenges are addressed so as to not introduce roadblocks with respect to related technologies such as image-guided interventions.