Hardware and software improvements of volume splatting

A comparative analysis of hardware and software improvements of volume splatting E. Verg´es , S. Grau and D. Tost May 18, 2006 Abstract This paper compares different hardware-based accelerations of the three classical splatting strategies: composite-every-sample,object-space sheet-buffer and image-space sheet-buffer. Specifically, we analyze the use of point sprites and 2D textures for the splat’s projection, the use of frame-buffer objects for the buffer composition and the use of GPGPU techniques for the transformation of voxels. We compare the efficiency of the hardware-accelerated splatting with 3D-texture mapping. In addition, we compare the rendering speed-up provided by hardware accelerations with two software-based space-leaping techniques: run-length encoding of labeled voxel models and voxel arrays. 1 INTRODUCTION The development of the last generations of graphics hardware has dramatically changed computer graphics research. Scientific visualization have gained many benefits from these advances, because they provide significant accelerations of rendering. Besides, in visualization, the ever-growing size of the datasets yields to a permanent demand of more computing power as well as hardware-based and software-based optimizations. A good example of the influence of graphics hardware technology in visualization is the current popularity of 3D-texture slicing for volume rendering. This technique consists basically of loading the entire volume in video memory and resampling it onto parallel proxy-geometry planes that are composed back-to-front to create the final image. Its speed is due to its use of the embedded trilinear interpolation hardware. Since this idea was first proposed in 1993 [CN93], the 3D-texture slicing technique has evolved and now, it incorporates volume shading [GK96], surface shading [WE98] and interactive classification [MHS99]. This evolution has been done in parallel and thanks to the graphics hardware improvements. The classical volume ray-casting can also be improved by a clever use of hardware capabilities. As an example, Mitsubishi’s VolumePro is a special purpose ray-casting implementation [PHK∗99]. UltraVis [Kni00] is an assembler raycasting implementation for specific types of CPUs. Finally, Kruger and Wester- 1 man [KW03] have proposed a multipass ray-casting with early ray-termination that takes profit from 3D-textures, z-test and fragment shaders. In the last five years, splatting, originally designed to render volume datasets, has gained popularity because its use has been extended to render point-based surface representations [PZBG00]. Splatting considers the volume or the surface as an array of overlapping kernels that are projected onto the screen plane in order to compose the image. Many attempts have been done to accelerate this technique using the GPU capabilities for both surface [RPZ02] and volume splatting [XC04]. Most of the research in this area focuses on the acceleration of the kernels projection using texture maps [CM93], 1D and 2D look-up tables (fast- Splats) [HMSC00] and point sprites [BK03]. The acceleration of other stages of the splatting pipeline such as the use of hardware-assisted opacity convolution [MSHC99], plane composition [NM05] and voxels geometrical transformation have been less addressed and restricted to image-space sheet-buffer splatting. Moreover, several software-based accelerations of splatting have been proposed, such as run-length encoding [KM01] and computation of render lists [MH01]. These techniques are primarily aimed at reducing the number of processed voxels. They have proved to greatly reduce the computational cost of splatting. However, the acceleration provided by these techniques in comparison with that provided by hardware is not very well known. In this paper, we address splatting of volume models. We analyze and compare different hardware-driven splatting accelerations for three different voxel model traversal schemes. We implement them and compare the results with software accelerations and 3D-texture mapping. 2 PREVIOUS WORK The splatting algorithm was proposed by Lee Westover [Wes89]. This algorithm gains its speed by exploiting the similarity of the kernel’s projection. In orthographic views, all the kernels have the same projection or footprint. Thus, the footprint can be computed once, in a pre-process, stored as a look-up-table and used for the projection of all the voxels. In perspective views, the footprints must be distorted according to the distance of the voxels to the observer. In the original approach of the algorithm, all the voxels are splatted directly in the image. This is why the algorithm is known as composite-every-sample. However, this method may cause color bleeding and sparkling artifacts because the visibility ordering of splats is imperfect. To correct this error, Westover [Wes90] proposed the object-space sheet-buffer splatting that splats the voxels slice-by- slice into sheet planes of the voxel model most parallel to the image plane and composites each sheet to the final image. This approach corrects color bleeding but it introduces noticeable popping up artifacts when the camera moves around the volume, because the sheet planes chosen change abruptly. Mueller and Crawfis [MC98] provided a solution to this problem that also enhances the approximation of the light transport inside voxels: the image-space sheet-buffer splatting. In this approach, the sheet-buffers are parallel to the image plane. 2 Therefore, voxels can contribute to more than one sheet. Different footprints corresponding to different intersections of the voxels with the sheet slab must be computed. This algorithm is composed of two main steps. In the first one, the voxel model is traversed, the voxels to be rendered are first transformed according to the viewing matrix, then depth-sorted and finally inserted into buckets such that each bucket covers the distance between suceesive sheet buffers. In the second stage, the sheets are processed in Front-to-Back (FTB) order. For each voxel in a sheet, the proper footprint is chosen according to a fast indexing scheme. In the image-space sheet-buffer splatting [MSHC99], early splat elimination is possible in FTB composition by subdividing the image into small tiles and avoiding to splat voxels that cover tiles that have already reached the maximum opacity. The detection of opaque tiles is efficiently performed using a hardware-assisted opacity convolution filter. One of the major advantages of splatting is that only relevant voxels must be splatted and empty and non-selected voxels can be skipped. This idea was first suggested by Yagel et al. [YESK95] for rendering Computational Fluid Dynamics (CFD). They suggested to construct a fuzzy set composed by an array of planes of the model and, for each plane, a list of voxels with their associated coordinates in the plane and their values. Crawfis [Cra96] introduced the idea of the ListSplat, a list of isosurface voxels that can be splatted directly without depth sorting because they are supposed to all be a homogeneous color. A similar idea is exploited in the RTVR system [MH01] that uses an intermediate array of slice-sorted RenderLists for each structure of the volume that store the voxels of the slices that are relevant for rendering. This structure is also used in the Two-level rendering proposed by Hauser et al. [HMBG01] [HBH03]. Mueller at al. [MSHC99] enhanced the efficiency of the view-aligned sheet-buffer splatting by organizing the selected voxels into buckets, each one corresponding to a sheet-buffer. The selection of the voxels for their insertion in the buckets is fast, based on a binary search in a per-value ordered list of voxels similarly to the work of Ihm at al. [IL95]. More recently, Orchard and M¨oller [OM01], proposed to use a list of adjacency data structure in which each non-empty voxel in a scan list is linked to the next non empty voxel in the scan-line. Finally, Kilthau and M¨oller [KM01] proposed to use run-length encoding (RLE) in order to skip empty voxels. They construct 24 RLE replications of the volume, which allows them to orderly traverse the volume according to any of the 48 orders. The main drawback to these two last approaches is their storage overhead. Many efforts have been done in accelerating splatting using hardware. One of the first proposed methods [LH91] [WG91] consists of approximating the splat by a collection of polygons, thus taking profit of the hardware-supported polygon rendering pipeline. Crawfis and Max [CM93] replaced the polygons by a 2D texture map. These approaches were tested in composite-every-sample traversals and orthographic projections in which only one footprint is necessary. Huang et al. [HMSC00] argued that image-space sheet-buffer splatting requires at least 128 footprint sections, which supposes over than 8MB texture maps storage. For radially symmetric splats, they propose to use a less-memory consuming one-dimensional table that holds the values of the splat along a radial line from 3 the splat center. Moreover, they explore directly copying into the image the block of pixels of a 2D footprint using BitBLT, but conclude that the image quality of this strategy is low. More recently, Xue and Crawfis [XC04] proposed two splatting strategies that work on the GPU. The first strategy consists of using a vertex shader program to generate and render quadrilaterals centered around the voxels center. This strategy works on previous generation hardware. In addition, it requires sorting the voxels along the viewing direction and it has high memory requirements. The second strategy, point-convolution rendering, first projects all the voxels as point primitives into an off-screen PBuffer with additive blending. Next, the GL convolution flag is activated and a texture is copied from the PBuffer using glTexSubImage2D such that each texel is a convolution between the PBuffer pixel and the kernel filter. This strategy is very efficient in terms of computational cost but it only renders x-ray style images for orthographic views, i.e images that do not provide depth clues, because opacity is simply accumulated without alpha-blending. Very recently, Vega-Figueroa et al. [VHFG05] propose to use Point Sprites to render neurovascular data. This reduces to one point per voxel the geometric processing tasks instead of the four-points needed for the quadrilaterals. This idea is also exploited in the GPU-based implementation of the image-space sheet-buffer splatting proposed by Neophitou and Mueller [NM05]. In addition, this paper proposes to use an OpenGL PBuffer object to store the buffers. It first splats onto an auxiliary buffer the density value of all the voxels of a slice using textured point sprites. Then, it classifies and shades all the pixels of the buffer using a fragment shader that computes the gradient vectors at the pixels on the basis of their density central difference. Finally, it composes the buffer into the final image. The authors use the early z-rejection test to eliminate empty-space pixels and those that are already opaque before the fragment processing. 3 HARDWARE ACCELERATIONS We have developed hardware-based accelerations for the three main splatting strategies: Composite-Every-Sample (CES), Object-space Sheet-buffer Splatting (OSS) and Image-space Sheet-buffer Splatting (ISS). Six main operations can be identified in the rendering pipelines of these three strategies: •Access to the voxel Value (V) •Object-space Gradient computation (G) •Shading (S) •Viewing geometrical Transformation of a voxel (T) •Voxel Splatting with or without α-blending (SP) •Bucket insertion (BI) •Buffer Composition (BC) 4 Figure 1: Rendering pipelines of the three splatting strategies from top to down: Composite-Every-Sample (CES), Object-space Sheet-buffer Splatting (OSS) and Image-space sheet-buffer Splatting (ISS). pis the point in world coordinates and p′in viewing coordinates, gis the gradient, vthe voxel value, cthe resulting color and Iithe i-th bucket and IFthe frame-buffer image. 8cm The rendering pipelines that we have implemented for these three strategies are depicted in Figure 1. The ISS pipeline is split into two-steps: the creation of the buckets and the splatting itself. In the second step, the point’s coordinates and color are fetched directly from their associated buckets. Observe that in order to get comparable results, the ISS pipeline that we have implemented is the classical one that shades in object-space and not the one suggested by [NM05], that computes pixel gradients and applies shading in the sheet-buffers. The 3D texture-mapping that we have used in the comparisons uses pre-shaded RGBαtexture loaded into the GPU memory. For splatting, the voxel access and gradient computations are done in the CPU. We have designed and tested hardware-based accelerations of the geometrical transformations (T), shading (S), splatting (SP) and buffer-composition (BC). 3.1 Splats projection Software-based splatting consists of rasterizing one point per voxel and compositing all the pixels of the footprint in the image buffer or the sheet-buffer around the rasterized point. We have analyzed two hardware-based splatting alternatives. First, the splat can be loaded as a 2D texture and a four vertices textured quad can be projected. If the splat is defined in object-space, the texture needs only to be loaded once, but the computation of the quad must be done in object space, which is computationally expensive. Alternatively, the splat can be computed in image-space, for each camera position and, thus, the textured quad can be defined directly in image space setting a 2D ortho 5 camera before its projection. An intermediate possibility between these two is to precompute the splat textures for a set of pre-defined cameras. However, this increases considerably the memory requirements. It is not even feasible in the case of ISS that requires at least eight footprints per camera. The second hardware based strategy is to use the glPointSprite. 3.2 Buffers composition One of the bottlenecks of the OSS and ISS strategies is the sheets composition with α-blending into the image buffer. We have substituted this step by using the Frame Buffer Objects extension (FBO). We splat the voxels using FBOs as offscreen render targets. Next, we activate α-blending; we define a 2D quad of the size of the frame-buffer; we associate to it the texture of the FBO and we render it in the frame-buffer. In this way, the sheet is composed with the frame-buffer using hardware acceleration. 3.3 Shading We have implemented three shading of increasing computational cost: •Per-value shading that simply uses the property value as a gray or RGB color RGB color •Volume shading that gets the emission and opacity from a Look-Up-Table (LUT) implementation of the transfer function •Surface shading that gets a material identifier from LUTs, and then applies the Phong model. Per-value shading is the simplest one. If voxel property values are stored as 1 byte (gray) or 3 bytes (RGB), it requires no computations. In order to accelerate volume shading, we have implemented the LUT as a 1D texture. The voxel property value is used as an index to this texture. In order to combine this shading with the texture-based and point-sprite splat projection, a multitexturing context must be set. The 1D texture is activated first and its corresponding color is blended to all the pixels of the 2D texture. Surface shading can also be handled using a hardware-driven LUT, for a pre-defined set of scalar product between the lighting direction and the gradient vectors. The computational cost of surface shading is, in this case, the same as volume shading without taking into account the table construction. Alternatively, the voxel property value is used as an index to a software-based materials LUT. The optical properties of the surface are set using glMaterial, the gradient vector is used to set glNormal and, if OpenGL lighting is on, Phong’s model is applied. The main drawback of this approach is that, even if only the surface color is set for each voxel, considering constant the specular color and exponent, it requires, for each surface voxel, at least either one glMaterial instruction or one glColorMaterial causing the material color to track the current color. per surface voxel. Other shading 6 models than Phong could be applied using vertex shaders, but we believe that this would not further enhance the speed of shading. 3.4 Viewing transformations In CES, the viewing transformations can be done using OpenGL hardwaresupported matrix operations if voxels are treated as quad or glPointSprite, as described in Section 3.1. In OSS, the viewing transformation is handled during the rendering of the Frame Buffer Objects (FBO) into the frame-buffers. However, in ISS, the viewing transformation is done in the CPU during the creation of the buckets. As an alternative, we have implemented a General Purpose Computation on the GPU (GPGPU) to load the selected voxels into the GPU as a floating point texture, transform them into the viewing coordinate system and generate a transformed list of voxels that can be used to insert each voxel in the proper bucket. 4 SOFTWARE ACCELERATION We have implemented two software accelerations to cull non-selected voxels. In the first strategy, we first traverse the full voxel model and compute a runlength encoding of the classification values of the model, such that each code is composed of a material identifier and the number of consecutive voxels sharing this material. During rendering, users select a set of materials for rendering; the run-length encoding is traversed but only the voxels corresponding to selected codes are processed. In order to allow any camera position, three run-length encoding are computed, one for each of the coordinate planes. Therefore,this model is suitable for OSS splatting or to create the buckets in ISS. It could be extended to CES, but with the drawback that 26 codifications should be computed in this case. The second strategy is used to accelerate ISS. Once the selection is done, the selected voxels are stored in a array that keeps the voxel coordinates and their value. For each camera position, instead of traversing all the voxels, only those of the array are processed. 5 RESULTS We have implemented and tested the different algorithms in our software platform Hipo. Table 1 indicates the abbreviations used for the seven different algorithms analyzed. To test the algorithms we have used six datasets that differ in various aspects: their size, that varies from the smallest 643digitized dinosaur (dino) to the larger micro-CT scanned rabbit femur (femur); their property value type 1 byte o 4 bytes (RGB α); their occupancy ratio measured as the number of 7 Abbreviation Algorithm CES Composite-Every-Sample H-CES Hardware-based Composite-Every-Sample OSS Object-space Sheet-buffer H-OSS Hardware-based Object-space Sheet-buffer ISS Image-space Sheet-buffer H-ISS Hardware-based Image-space Sheet-buffer H-TM Hardware-based 3D Texture Mapping Table 1: Algorithm names and abbreviations. non-empty voxels divided by the total number of voxels; being or not voxelized isosurfaces such as dino and prostate. Table 2 summarizes the characteristics of the different datasets. Figure 2 shows a rendered image of each of the datasets. All the simulation results have been taken on a Pentium IV PC at 3.2 GHz with 3.5 GB of memory and an ATI Radeon X800 Pro graphics card. We have used the OpenGL library on a GNU/Linux operating system. Time is measured in seconds. Dataset Dimensions Size Type Oc.ratio Dino 87*39*62 210366 bool 9.04% Prostate 194*202*256 10032128 float[4] 11.97% Aneurism 256*256*256 16777216 char 0.37% Skull 256*256*256 16777216 char 5.70% Engine 256*256*128 8388608 char 20.17% Femur 485*509*633 156265545 char 32.40% Table 2: Characteristics of the datasets, from left to right: name, size, property value type and occupancy ratio (OR). Table 3 shows the rendering cost of the splat projection of 16 millions of voxels using the CPU-based strategy, the 2D image-space textured quad and the glPointSprite for the different datasets (see Section 3.1). For the two latter strategies, the costs include the texture construction. As it can be observed, the hardware assisted modes are far more efficient than the CPU-based mode whose efficiency slows down as the splat size increases. The textured quad and glPointSprites have similar costs although the textured quad is more efficient. This difference is due to the fact that the implementation of sprites in our Linux operating system is not very efficient on Radeon cards. The same simulations on a NVidia card showed a reduction in the cost of of the glPointSprite to approximatively the same value as the textured quad. Table 4 shows the cost of the CPU-Based buffer composition in comparison to the use of FBOs, for three different image sizes (see Section 3.2). As expected, the larger are the images, the more the FBO accelerates the buffer composition. A drawback of the FBOs is that they require the image dimensions to be powers of 2. The larger image that we have been able to process as a whole in our PC is 1024x512. 8 Figure 2: Rendered image of the datasets. Left column from top to bottom: dino, skull and femur. Right column from top to bottom: prostate, engine and aneurysm. The first row of Table 5 shows the shading cost based on a LUT for volume shading or quantized pre-processed surface shading using software-based LUTs in comparison to 1D-texture LUTs. The LUT have a size of 256. Shading has been applied to 16 millions of voxels. The results are of the same order of magnitude and the CPU-based mode is even slightly better. In the second raw of the table, we show the cost of surface Gouraud shading with pre-computed gradients and one light source, implemented via CPU and with GL lighting. Again, the differences are not very significative, because the geometrical computations (one scalar product and four multiplications) are not very intensive. Finally, Table 6 shows the cost of applying the viewing transformation to 100.000 vertices using the GPGPU program described in Section 3.4. The overall CPU-based process cost is 1.69 10−2seconds. As it can be seen, the geometrical transformation cost is far more efficient in the GPU but the read-pixel operation needed to insert the voxels the buckets are a bottleneck, that globally makes the GPU assisted process slower than the CPU-based one. 9