Spectral Volume Rendering using GPU-based ... - Semantic Scholar

preprint, original article available at http://dx.doi.org/10.1007/s00371-006-0028-0 The Visual Computer 22(8)

Magnus Strengert · Thomas Klein · Ralf Botchen · Simon Stegmaier · Min Chen · Thomas Ertl

Spectral Volume Rendering using GPU-based Raycasting

Abstract Traditional volume rendering does not incorporate a number of optical properties that are typically observed for semi-transparent materials, such as glass or water, in the real world. Therefore, we have extended GPUbased raycasting to spectral volume rendering based on the Kubelka-Munk theory for light propagation in parallel colorant layers of a turbid medium. This allows us to demonstrate the effects of selective absorption and dispersion in refractive materials, by generating volume renderings using real physical optical properties. We show that this extended volume rendering technique can be easily incorporated into a flexible framework for GPU-based volume raycasting. Our implementation shows a promising performance for a number of real data sets. Particularly, we obtain up to 100 times the performance of a comparable CPU implementation.

increasingly large data sets; but this desire also refers to new methodologies, i.e. we expect an optimal tool to be very flexible and to allow for fast and painless integration of innovative new volume rendering techniques. At the moment, actual volume rendering implementations are almost exclusively based on slice-based methods where axis-aligned or viewport-aligned textured slices are blended together to approximate the volume rendering integral. Without doubt, these slice-based approaches have their benefits. However, slice-based implementations are rasterization-limited and can hardly be optimized from an algorithmic point of view. Furthermore, when applying a perspective projection the integration step size will vary along viewing rays when using planar proxy geometries, leading to visible artifacts. Finally, they do not easily allow for an implementation of techniques with viewing rays changing direction Keywords GPU-Raycasting · Spectral Volume Rendering · as it occurs in refracting volumes. Slice-based techniques, Kubelka-Munk Theory therefore, do not fulfill all of the criteria defined above. On the other hand, the recent advancements in highly programmable graphics processing units provide an ideal platform for efficiently mapping raycasting-based volume 1 Introduction rendering to graphics hardware. A fragment program based Scientific volume visualization plays an important role in raycasting does not suffer from any flexibility issues and, both academia and industry. Depending on the application, therefore, fulfills the second criterion. Furthermore, graphics various visualization techniques are being employed. In or- hardware manufacturers urge us to use any novel fragment der to obtain a productive system that is universally appli- program feature, promising that this functionality will becable — in medicine, geosciences, engineering, etc. — it come very fast in future generations of their hardware. It can is necessary to integrate these techniques into a single tool. thus be assumed that GPU-based raycasting is also futureFurthermore, a volume rendering tool is supposed to be future- proof from a technical point of view and, accordingly, meets proof. Of course, this refers to technical advances, like e.g. in all criteria defined. Although, the traditional volume rendering integral algraphics hardware, to guarantee best performance also for lows one to create renderings that are suitable for the analMagnus Strengert, Thomas Klein, Ralf Botchen, Simon Stegmaier, ysis of typical volume data sets in a large number of appliThomas Ertl cations it does not incorporate important optical properties VIS Group, University of Stuttgart, Germany that are typically observed for real semi-transparent materiUniversitätsstraße 38, 70569 Stuttgart als, such as glass or water. This includes the reflection, reE-mail: {strengert, klein, botchen, stegmaier, ertl} fraction, selective absorption, scattering, or polarization of @vis.informatik.uni-stuttgart.de light. Therefore, many extensions to the basic volume renMin Chen dering equation have been proposed in order to allow for Department of Computer Science, University of Wales Swansea Singelton Park, Swansea, SA2 8PP, United Kingdom more sophisticated global illumination effects, for example, E-mail: [email protected] subsurface scattering, translucency, or volumetric shadows.

2

However, a number of optical characteristics like dispersion or selective absorption are based on the spectral properties of light and the wavelength-dependent characteristics when interacting with different media. While these effects have gained much attention in photorealistic image synthesis, in the field of volume rendering few publications on this subject can be found in the literature. This paper is an extension to our previous work [40] that describes our experiences in developing a framework for GPU-based raycasting and a series of single-pass volume shaders for both standard and non-standard volume rendering techniques. Specifically we have included the following topics. First, we have extended the framework to support multi-pass rendering in off-screen buffers to account for complex volume rendering algorithms that cannot be realized within a single fragment shader. Second, we describe the implementation of spectral volume rendering including wavelength-dependent refraction and absorption based on measured material properties in the case of volume raycasting. Last, we compare our GPU implementation to a CPUbased solution realized in the same framework.

2 Related Work Most of the work in direct volume visualization in recent years has been focused on texture-based approaches. First introduced by Cullip and Neumann [7] and Cabral et. al [6], the basic raycasting concept is realized by sampling the volume data using a stack of, typically, planar slices as a proxy geometry and approximating the evaluation of the volume rendering integral by blending the textured slices in frontto-back or back-to-front order in the framebuffer. This pixelparallel processing of the viewing rays during rasterization of the proxy geometry exploits the unmatched bi- and trilinear interpolation capabilities of modern graphics hardware and is the primary reason for the unsurpassed speed and success of this method. Many enhancements to this simple approach have been proposed that exploit more advanced texture mapping capabilities of today’s graphics hardware to increase the interactivity and applicability of the method, e.g. [32, 42]. The introduction of multidimensional transfer functions [20] and pre-integrated transfer functions [10] has significantly improved the quality of renderings that can be achieved. Also, some acceleration techniques proposed for the original raycasting approach, such as early ray termination and empty space skipping [26], or hierarchical acceleration structures [5, 13, 25] have been successfully adopted to texture-based direct volume rendering. Nevertheless, it is still much harder and requires considerably more effort to integrate such techniques into a slice-based volume renderer compared to the implementation in an obviously much more flexible software-based raycasting code. Although, as we have seen, the use of graphics hardware has a long standing tradition in volume rendering, only recent advances in the programming model and the incredible pace at which the performance of graphics processors has

Magnus Strengert et al.

increased over the last years, provide the means of carrying the benefits of raycasting to hardware-accelerated volume rendering. GPU-based implementations of the raycasting algorithm that overcome some of the shortcomings of the traditional volume slicing approach have been presented for both structured and unstructured grids. These solutions can be divided in two classes. The first [22, 35, 41] performs multiple rendering passes — similar to traditional slice-based volume rendering — in order to traverse the volume and stores intermediate results computed on a per-fragment basis in temporary buffers that are accessed in subsequent rendering passes. Unfortunately, these approaches suffer from the lack of high accuracy blending support in contemporary graphics hardware or have to resort to ping-pong rendering to overcome this problem. The second group executes the whole raycasting algorithm in a single rendering pass, exploiting the functionality of dynamic looping available on recent graphics processors. This was first shown as a simple technology demo by NVIDIA [30] that demonstrates the use of advanced fragment shader functionality to implement a basic single-pass raycasting algorithm for regular volume data in a single fragment program. In [40] we presented a flexible framework for single pass GPU-raycasting that takes advantage of the easily extensible raycasting approach to demonstrate a number of non-standard volume rendering techniques, including translucent material and self-shadowing isosurfaces, as well as an acceleration technique based on exploiting inter-frame coherence. Hadwiger et al [15] presented a GPU-raycasting system for isosurface rendering. They employ a two-level hierarchical representation of the volume data set for efficient empty space skipping and use an adaptive sampling approach for iterative refinement of the isosurface intersections. Spectral color models have been deployed in computer graphics since the 1980s. For example, the reconstruction of the spectral power distribution (SPD), from RGB representations were studied by [11, 17, 28]. Peercy et al. [31] considered a hardware implementation of full-spectral rendering on a Silicon Graphics RealityEngineT M . Johnson and Fairchild developed an OpenGL extension for full-spectral local illumination [18] For a more extensive coverage of spectral modeling and rendering, readers may refer to [12, 16, 36, 38] Noordmans et al. [29] were the first to introduce the spectral color model to volume visualization, and demonstrated that a spectral volume renderer produced physically more realistic visualization in comparison with the traditional RGBα model. Bergner et al. presented a spectral volume rendering system, which facilitates post-illumination and interactive data exploration [3]. They recently examined the incorporation of spectral color representation in several rendering algorithms, including raycasting, volume splatting and Fourier volume rendering [4]. In the field of computer graphics, one particular optical model, namely the Kubelka-Munk model (KM model) [24, 39], was often used in conjunction with spectral color com-


putation. The model, which theorizes the passage of light among multiple inhomogeneous layers of translucent materials, is often referred to as a one-dimensional radiosity model. Haase and Meyer [14] employed the Kubelka-Munk model for pigment mixing in order to improve graphical realism. Dorsey and Hanrahan [9] used the model for simulating metallic patinas, a type of weathering effect of metallic objects. Others deployed the model for non-photorealistic rendering through the physically-based simulation of different paints, such as watercolor [8], wax crayon [37] and oil paints [2]. Recently, Abdul-Rahman and Chen [1] presented a spectral volume rendering integral based on the KM model. They highlighted some relative merits of the Kubelka-Munk theory over the Lambert-Bouguer law in terms of its built-in distance attenuation, and its capability of determining the opacity and transparency optically according to the absorption and scattering properties, in addition to the advantages of more accurate simulation of selective absorption, the use of measured spectral data and post-illumination. Apart from the early work by Peercy [31] and the spectral volume splatting by Bergner et al. [4], most implementation of spectral rendering have been carried out primarily in software.

3

The sums of the two infinite series of reflectance and transmittance define the composited reflectance R and transmittance T. Let Ri and Ti be the accumulated sums of compositing reflectance and transmittance from layer 1 to layer i. Let, Ri+1 and Ti+1 be the sampled reflectance and transmittance of layer i + 1. The compositing reflectance and transmittance from layer 1 to layer i + 1 is as follows: T2i Ri+1 1 − Ri Ri+1 Ti Ti+1 Ti+1 = 1 − Ri Ri+1

Ri+1 = Ri +

(1)

This facilitates a discrete volume rendering integral with built-in one-dimensional radiosity. As sampled the reflectance R and transmittance T values depend on the thickness of the layer (or the sampling distance) δ , it is common to calculate R and T from the sampled absorption attribute K and scattering attribute S as: sinh(bSδ ) a sinh(bSδ ) + b cosh(bSδ ) b T= a sinh(bSδ ) + b cosh(bSδ ) R=

(2)

where 3 Kubelka-Munk Theory in a Nutshell The Kubelka-Munk model is basically a two-flux approach to the general radiation transport theory, which considers two radiation fluxes, namely reflectance (R) and transmittance (T ), that pass through a continuous medium in two opposite directions. While an early version of the model assumed that the medium is completely homogeneous [24], Kubelka later extended the KM model to inhomogeneous layers of paint-like turbid substances, and each elementary layer is assumed to be isotropic and completely diffuse [23]. Consider a volume data set to be evaluated using raycasting in direct volume rendering. Given a series of samples, s1 , s2 , . . . , sn , along a ray cast from the viewing position at a regular interval δ , we approximate the intersection volume between the ray and the data set as a series of homogeneous elementary layers of thickness δ . The KM model can thereby be deployed, and we have two fluxes I f and Ib flow in opposite directions, that is, I f flows forwards following the ray direction, and Ib flows backwards. Note that, these two fluxes were originally called downward and upward fluxes respectively. Let the reflectance and transmittance for each layer be Ri and Ti , i = 1, 2, . . . , n, a light flux from layer i − 1 passes through layer i with the portion Ti , then the next layer with the portion Ti Ti+1 , and so on. Meanwhile, a portion of flux of the amount Ti Ri+1 Ti will be reflected from layer i + 1 and passes backwards though layer i again. The infinite process of interaction between layers can be considered as onedimensional radiosity [9] — a limited form of global illumination.

a=

S+K , S

b=

p

a2 − 1.

Note that the absorption and scattering coefficients K and S are functions of the wavelength and derived from material properties [14, 2]. 4 A Framework for GPU-Raycasting The system for volume raycasting proposed in this paper consists of two major parts: a framework written in plain C and based on OpenGL and GLUT, and a set of shaders that are flexibly loaded and modifiable at runtime. Some examples are given in Secs. 6 and 8. The system is portable between MS Windows and Linux. Independent of the selected volume shader, the framework generates the primary rays, initializes the textures, and provides various functionality like event handling and transfer functions. Besides setting up the rays, this is all standard functionality also found in any slice-based volume rendering application. Furthermore, the framework supports the use of offscreen buffers and multiple render targets for complex shaders, that cannot be carried out in a single-pass. This functionality is built on top of the OpenGL EXT FRAMEBUFFER OBJECT and ARB DRAW BUFFERS extensions. As an example for the use of this functionality we have implemented a completely GPU-based empty space skipping algorithm for GPU-raycasting [19] and the spectral volume rendering approach described in Sec. 6. The volume data set is expected in unsigned char or unsigned short format and stored in a 3D texture. Since

4


6.8/9.2/12.9

7.2/7.4/19.3

3.3/3.4/8.0

6.8/8.5/15.9

Figure 1 Examples for different volume shaders rendering a CT data set (160×430×183 voxels). From left to right: opaque isosurface with self-shadowing, pre-integrated volume rendering, isosurface combined with pre-integrated volume rendering, and shaded volume rendering. The minimum, average, and maximum framerates are given in the lower left corner.

many volume rendering techniques, e.g. lighting or refraction, require gradient information, these can be precomputed using central differences or a 3×3×3 Sobel operator and optionally stored together with the scalar data. Furthermore, an additional filtering step can be applied to achieve smooth gradients. As shown in [34] filtering can be significant for the rendering quality, especially for refractive volume rendering methods or data sets that have high noise. All precomputed gradient components are quantized to 8 bits and stored in the RGB components of the RGBA texture already holding the scalar volume data in the alpha channel. Since the quantization poses a loss of accuracy, it is possible to compute gradients on the fly in the fragment shader if higher quality gradients are required. In this case the single component volume texture is sufficient and larger data sets can be loaded. For the description of the rendering process we first define the following volume properties. Let Sc denote the number of slices in c-direction and Dc the slice distance in cdirection with c ∈ {x, y, z}. Then ec = Ec = Cc =

Sc Dc , (Sc Dc )/ max({ex , ey , ez }) , and Ec /2

(3)

define the volume extents ec , normalized volume extents Ec , and the volume center Cc . For displaying the volume, we bind a user-selected volume shader and as proxy geometry render an axis-aligned box translated by (−Cx , −Cy , −Cz )T with backface-culling enabled. One corner of the bounding box is located at the origin, the opposite one at (Ex , Ey , Ez )T . The box defines the normalized bounding box of the volume considering both the number of slices and their distances. For each bounding box vertex we define texture coordinates identical to the vertex positions. Each fragment can thus access the respective ray’s entry point by just looking at its interpolated texture coordinates and — in combination with the camera position — easily compute the parametric ray of sight. The fragment attribute fragment.position stores the (x, y) window coordinates of the fragment center relative to the lower left corner of the window and fragment’s z window coordinate and not the fragment’s coordinates in object

space; accessing the object space coordinates is, therefore, not possible without the apparently redundant setting of texture coordinates equal to vertex coordinates. A significant benefit from this approach is the inherent minimization of the number of rays that are generated. 5 Raycasting on the GPU In the last years programmability of workstation and consumer level graphics hardware has evolved at an increasing pace. Driven by the steadily growing demands of the game industry, performance of modern graphics processors has exceeded the computational power of CPUs both in raw numbers and in their extraordinary rate of growth. Abandoning the simple fixed-function pipeline which was the characteristic feature of graphics processors only a couple of years ago, today’s GPUs have evolved into very sophisticated, highly programmable SIMD processing units. With the advent of graphics processors supporting the new features of DirectX Pixel Shader 3.0 and equivalent OpenGL extensions, namely dynamic looping and branching, GPUs are becoming more and more “general purpose processing units” comparable to the CPU. The basic raycasting approach fits very well into the intrinsically parallel stream processing semantics of these new fragment processors. For each pixel of the final image a single ray is traced independently through the volume. Therefore, a fragment program implementation of this simple algorithm working on the fragments generated by rasterizing a polygon covering the screen space area of the volume’s projected bounding box is sufficient to compute the correct result. The volume rendering integral for each pixel is then approximately evaluated by sampling the ray at a finite number of positions inside the volume. The contributions of those samples along the ray are accumulated to the overall chromaticity and opacity. Applying an appropriate optical model every desired kind of interaction between light and the volumetric object can be realized. Unfortunately, due to limited capabilities of graphics hardware most of the elaborate optical models already proposed a decade ago [27] have not or only with considerable effort and overhead been integrated


into slice-based volume renderers. In contrast, it is often very easy to include them into a raycasting system [21, 33]. Fig. 1 shows a number of images that have been rendered with the presented GPU-based volume raycasting framework using different types of volume shaders. Performance data for each shader was measured for a 512×512 viewport when rotating the volume once around the y-axis. The sampling step size was chosen to be approximately half the smallest voxel distance. All performance measurements were conducted on a standard PC equipped with a 3.4GHz Pentium 4 processor and a NVIDIA GeForce 7800 GTX based graphics accelerator card. Note that the performance drops significantly in the case of mixed volume and isosurface rendering due to the high overhead of conditional branching. First we start with a brief overview of the GPU-based raycasting. For the primary rays the direction for the respective fragment/pixel together with the entry point, i.e the first intersection of the eye ray with the volume’s bounding-box, has to be determined. Then in each step of the loop the actual interpolated data value for the current sample point is fetched from a 3D texture map and the already accumulated color and opacity values for the fragment are updated according to the chosen optical model or rendering style. In the simplest case, this could be a dependent texture look-up combined with alpha-blending for a basic volume rendering transfer function. Then, the sampling position is advanced along the ray by a specific step size. The loop is terminated either if the ray has left the volume or if some other criteria — depending on the chosen optical model — is met, e.g. early ray termination due to high accumulated opacity or the first hit semantic if an opaque isosurface is encountered. The complete code for a simple fragment program based volume raycaster is shown in Fig. 2. As mentioned in section 4, the bounding box with attached texture coordinates is rendered as proxy geometry to generate all necessary fragments. Each of these fragments corresponds to a single ray with the starting position given by the interpolated texture coordinates of its fragment. The direction for the initial rays is given by subtracting the position of the camera from the ray starting position. Not shown in the code is that the camera position can be easily derived by reversing the translation dictated by the modelview matrix. Since the camera is initially located at the origin, this translates to a single instruction in the fragment program. Using the ray’s entry point and the computed ray direction, it is now possible to sample the volume with a set of two nested REP loops, allowing for an overall number of 65,025 iterations. Note that current GPUs only support a maximum of 255 iterations per loop, which would not be sufficient for traversing reasonably sized data sets. The actual volume sampling is then straightforward since we basically can now follow the software approach without sacrificing the parallel-processing power of modern GPUs. Inside the loop the first step is to sample the volume data set at the current position. Depending on the optical model and the

5

chosen type of classification, the obtained scalar value and optional gradient can be processed further. In this example, we use standard front-to-back alphablending to accumulate the opacity in combination with a pre-integrated transfer function. Note that in the implementation presented, the integration does not have to resort to low-accuracy framebuffer blending operations and can take advantage of highly accurate 32 bit floating point computations. For the special case of pre-integration, we benefit of the flexibility of the GPU-raycasting approach. This tech-

# Get start position MOV position, startPosition; # Determine ray direction SUB direction, position, camera; NRM direction, direction; MUL direction, direction, scaleFactors; # Lookup scalar value in 3D volume texture TEX scalar.g, position, texture[0], 3D; # Move one step forward along ray direction MAD position.xyz, direction, step, position; # Start ray traversal MOV dst, 0.0; REP 255; REP 255; # Lookup scalar value in 3D volume texture TEX scalar.r, position, texture[0], 3D; # Lookup in pre-integrated transfer function TEX src, scalar, texture[1], 2D; # Perform front-to-back blending SUB texblen.r, 1.0, dst.a; MAD_SAT dst, src, texblen.r, dst; # Move one step forward along ray direction MAD position.xyz, direction, step, position; # Test if outside volume and exit loop SGE temp1.xyz, position, volExtentMin; SLE temp2.xyz, position, volExtentMax; DP3 inside.x, temp1, temp2; SEQC inside.x, inside.x, 3.0; BRK (EQ.x); # Save current scalar value for pre-integration MOV scalar.g, scalar.r; ENDREP; # Exit loop if volume left BRK (EQ.x); ENDREP; # Write output color to framebuffer MOV result.color, dst; Figure 2 Basic fragment program code of a volume raycaster using pre-integrated lookup tables. Implemented using the GL NV fragment program2 extension.

6

nique uses a lookup-table dependent on two consecutive volume samples in order to increase the rendering quality [10]. The ray traversal in one single fragment program allows us to store intermediate values and reuse them in the following iterations. In the example this refers to the MOV instruction at the end of the loop. Thus, pre-integration can be implemented with only one additional texture lookup. On the opposite, slice-based approaches cannot easily communicate information between the slices and therefore a common solution would need to sample the texture twice for each slab before performing the lookup. The remaining lines of code inside the loop perform the actual ray traversal and the test whether the ray has already left the volume. The test compares the current position to the minimal and maximal extent of the volume and uses the dynamic flow control of Shader Model 3.0 to exit the loop. However, there are still two subtle issues that have to be addressed. First, when leaving the volume, the sample lying outside the volume must be projected back along the ray direction onto the respective bounding box back side and the volume rendering integral must be evaluated using the reduced distance. However, a memory-efficient 2D pre-integration table is only valid for a single fixed interval length and thus a 3D pre-integration table would be required, consuming an immense amount of memory. We therefore accept the error introduced by terminating the sampling just before the volume is left and neglect a correct boundary handling as is commonly done also when using slice-based volume rendering. Second, as seen in Sec. 4, the normalized bounding box extents are given by the values Ec . Those values correspond to the geometric extent of the volume in each axis, but they cannot directly be mapped to the maximum extent of the volume in texture space. This originates from the fixed interval [0..1] used for texture coordinates in all dimensions independently of the number of voxels provided on each axis. Thus, before accessing the volume data the coordinates need to be scaled by a factor of E1c that is passed to the shaders as program parameters. The scaling can then be either accomplished within the inner loop — which is very costly — or outside the loop by scaling the entry point and the direction vector. It should be noted that the latter must be done after normalization. The resulting raycasting then takes place in texture space and neither the ray direction nor the sampling position can be used any more for geometric computations in object space, which are required for rendering the background. To optimize the rendering performance we also implemented early ray termination, which can be easily added to the code as a single conditional BRK instruction. Unfortunately, this did not improve the overall rendering performance which can likely be explained with the high overhead for conditional branching in the fragment shader as well as the fact that fragments are processed in large groups and therefore an equally large number of coherent fragments is necessary to gain any positive effects. Although the performance benefit that can be achieved by early ray termination


(a)

(b)

Figure 3 Illustration of sample distances in direct volume rendering approaches. In contrast to raycasting (b), the sampling distance is not constant in slice-based volume rendering (a).

is hard to quantify — it strongly depends on the structure of the data set, the chosen transfer function, the current view direction, and the sampling rate — we found that in most cases very little is gained. For typical transfer functions used for pre-integrated volume rendering we have measured everything from a 20% performance loss up to a 60% performance increase. The fact that the principal rasterization complexity is the same for both raycasting and conventional slice-based rendering poses the question about why a single pass volume raycasting should be superior to 3D texture-based slicing. First, it obviously eliminates the necessity for intermediate buffer reads and writes. Second, since basically only a single polygon has to be rendered in order to generate the necessary fragments, raycasting exhibits a very low geometry processing and fragment generation overhead. And, third, raycasting allows for adaptive step sizes (including early ray termination and empty space skipping), and by definition samples the volume at equal distances (Fig. 3), thereby avoiding artifacts [25]. Other optimizations based on the history along the sampling ray also benefit from the flexibility of our raycasting approach as will be demonstrated in Secs. 6 and 7. Moreover, raycasting has a much higher accuracy than slice-based rendering since the entire algorithm is performed in full 32 bit floating point precision. In contrast, slice-based rendering suffers from the lack of high accuracy floating point blending in today’s graphics hardware. Although, it is possible to avoid framebuffer quantization and to emulate full precision floating point blending by ping-pong rendering to an off-screen floating point texture target instead of the framebuffer, this approach further increases the cost associated with buffer accesses and reduces the amount of video memory that can be used to store volume data. 6 Spectral Volume Raycasting In this section we extent the so far presented GPU-raycasting approach by the spectral volume rendering method based on the Kubelka-Munk theory described in Sec. 3. Standard volume rendering using alpha compositing does not account for selective absorption of light, hence the color accumulation only depends on the assigned opacity of the material and all intensity components of the light are attenuated equally. In contrast, the amount of light absorption by real media varies


over the spectrum. That means with an RGBα model, it is not possible to combine different colored semi-transparent layers in an optically correct manner. An example for this is given in Fig. 4a. The semi-transparent homogeneous blue cube absorbs most wavelengths of the red spectrum, thus the yellow squares of the background appear to be green when seen through the cube and the red parts of the background completely loose their color information. Neither effects are reproduced using α -compositing, as shown in Fig. 4b. The RGB transfer function used in (b) was computed from the spectral properties of the cube’s material and the α transfer function adapted to match the appearance of the cube shown in (a). To include those spectral effects into the volume raycasting framework described in Secs. 4 and 5 we need to take the wavelength dependent behavior of the media into account. Therefore, we use a equidistant sampling of the spectrum of visible light in the range from 400nm to 700nm at 10nm intervals, resulting in a total of 31 distinct wavelengths λ that have to be considered independently. Similar to the transfer function used in traditional volume rendering for each scalar value of the volume, a corresponding set of spectral absorption and scattering coefficients Sλ and Kλ has to be defined that maps data values to material properties. The volume sampling then proceeds as usual, with the difference that instead of accumulating color and opacity values, for each wavelength a reflectance value according to the reflectance and transmittance encountered along the sampled path in the volume is computed. In other words, the classification by the transfer function lookup and the volume integration based on alpha blending is replaced by a lookup of the Kubelka-Munk coefficients Sλ and Kλ for the sampled material or scalar value, the evaluation of the Eqs. (2) for the reflectance and transmittance, and the Kubelka-Munk composition as defined by Eq. (1). Although, a straight forward implementation of Eqs. (2) can be done very easily in the fragment program, this is not very efficient since we would need to use a lot of costly shader instructions like the square root and the hyperbolic

7

sine and cosine. Fortunately, when the thickness δ of the layer is known, which is the case for raycasting with a constant step length, R and T can be pre-computed for a finite set of scalar values and wavelengths and stored in a twodimensional texture map. Since both functions are bounded between 0 and 1 it is possible to store these values with high precision in a 16 bit floating point texture. When the ray exits the volume bounding box we have to add the contribution of the background. We assume the background to be a completely opaque layer represented by a reflectance map Rb (λ ). This map is obtained from an RGB image using an algorithm by Glassner [11] that computes the spectral representation of an RGB value under controlled lighting conditions. Note that since there are an unlimited number of possibilities to define a spectrum for an RGB color triplet and the limited number of frequencies that we sample in the spectrum, RGB values are not accurately reproduced in this approach. However, the approximation is sufficiently close for our application. The final spectral reflectance map RI (λ ) of the image can then be computed according to Eq. (1) by compositing Rb (λ ) and the accumulated reflectance of the volume Rn (λ ) as RI (λ ) = Rn (λ ) +

T2n (λ )Rb (λ ) . 1 − Rn (λ )Rb (λ )

(4)

As a last step, the computed reflectance map has to be lit with an appropriate light source to produce the final color values of the image. This is done by applying the spectrum L(λ ) of the CIE standard illuminant D65 representation of natural daylight [43]. Which gives us a spectral representation of the image: Iλ = RI (λ )L(λ ) .

(5)

Finally, the spectral image has to be converted to RGB color space for displaying the result. This is accomplished by converting the color spectrum Iλ to the XYZ color space using the tristimulus color matching functions, x( ¯ λ ), y( ¯ λ ), z¯(λ ), defined for the CIE 1964 standard colorimetric observer [43] by IXYZ = ∑ Iλ (x( ¯ λ ), y( ¯ λ ), z¯(λ ))

(6)

λ

and transforming these values to RGB space IRGB = MXYZ→RGB IXYZ

(a) Kubelka-Munk

(b) RGBα

Figure 4 Comparison of our spectral volume raycasting based on the Kubelka-Munk model (a) to traditional RGBα blending (b) to show the effect of selective absorption. Note how the appearance of color changes when seen through the blue colored cube.

(7)

using the color conversion matrix MXYZ→RGB given in [11]. The images that can be rendered using the algorithm described above in general look flat and unrealistic, especially for nearly opaque materials. Like in RGBα volume rendering we can improve this by adding a shading term. Since the Kubelka-Munk model assumes isotropic light scattering that obeys Lambert’s cosine law we can multiply the local reflectance coefficient by a diffuse shading term without violating the model.

8


Although it is theoretically possible to do the spectral raycasting for all wavelengths in parallel in a single fragment program, there are several drawbacks that prevent the use of this approach. First, the number of temporary registers that would be necessary to store all intermediate values would be very high. Second, since there is no possibility for indirect addressing of registers in fragment programs the number of lines of code necessary to implement this program would be very high resulting in code that is error prone and unmaintainable. And third, and most importantly, this would even for fairly large sampling steps exceed the maximum instruction limit of 65k executed fragment shader instructions per fragment. Therefore, we treat each wavelength separately using a multi-pass solution. Since the summation of color matching functions and the following color transformations are linear operations, additive blending can be used for accumulating the contributions of the 31 rendering passes in a floating point render target. Unfortunately, blending in floating point framebuffer objects is only supported for 16bit precision. But, as only color values are blended the precision should be sufficient. Of course, another possibility would be to use a pair of 32bit floating point textures and do the addition in a ping-pong fashion, but with a much higher overhead. Ignoring the refractive effects, as described in the following section, it would be also possible to treat four frequency bands in parallel, to account for the vector processing capabilities of the GPU.

7 Refractive Volumes Besides the selective absorption property that can be modeled based on the Kubelka-Munk theory there are other important effects that can be observed for transparent media. Fig. 5a shows a photograph of a cube made of plexiglass in a box with differently colored side walls. There are two important phenomena that can be observed in this image. First, the parts of the blue-white colored back plane that can be seen through the cube are distorted due to the refraction of light and the side walls of the box are mirrored in the left and right backwards facing sides of the cube due to total reflection that occurs in certain cases on the transition from an optically denser layer to an optically lighter material. Another effect that can be often observed is dispersion. This is a phe-

nomenon that also arises from the variation of the refraction coefficient of the material over the spectrum of light. Since different wavelengths are refracted to different degrees, an achromatic light ray is split into rays of chromatic light that follow different paths through the media. Typically this phenomenon is associated with prisms. Although the KM-model assumes that the entire volumetric domain of colorant has the same refractive index, it should not prevent the introduction of refraction in such a volume rendering integral in an approximative manner. Firstly, having refraction in the KM model is not in any way less accurate than having refraction in the traditional volume rendering integral, since the traditional volume rendering integral is based on the Lambert-Bouguer law [43], which also assumes that the entire volumetric domain is of the same refractive index. Secondly, if there were refraction between different elementary layers, changing the direction of a flux according to the relative refractive index between the two layers is more accurate than no directional change of the flux. The reason of this is obvious as the dominant direction of transmittance in the phase function of a volume is largely determined by the field of refractive index. Implementing refractive volume rendering is straight forward with the presented raycasting approach, requiring only few modifications. The first is an additional refraction coefficient that has to be stored in the material texture for each of the 31 wavelengths of our spectral rendering model. The second is the update of the current ray traversal direction. Similarly to storing the scalar value of the previous sampling in the case of pre-integrated volume rendering (cf. Sec. 5) we keep the refraction coefficient ηi−1 (λ ) obtained for the last sample. Then in the current integration step the new ray direction di is given by q di = η di−1 − η (di−1 · ni ) + 1 − η 2 (1 − (di−1 · ni )2 ) ni , (8)

where η = ηi−1 (λ )/ηi (λ ) and ni is the volume gradient at the current sampling position, if necessary re-oriented to point in the direction opposite to di . Both di and ni are assumed to be normalized. If the angle between the volume gradient and the incoming direction of the ray exceeds the critical angle ψt = arcsin η −1 the incoming ray is reflected from the material interface, rather than being refracted. This total reflection is indicated by a negative value below the square root in Eq. (8). In this case di is set to di = di−1 − 2(di−1 · ni )ni .

(a)

(b)

Figure 5 Refractive volumes. (a) Photograph of a cube made of plexiglass in a box with differently colored side walls; (b) Volume rendering of a cube data set using spectral refractive volume rendering.

(9)

Unfortunately there are some problems with total reflections. It is possible and in the case of real data sets not unlikely that very long light paths can occur when a ray is reflected multiple times inside the volume. Such problems are typically caused by noise and inaccurate gradients or numerical inaccuracies in the computation. This not only makes the computation very slow, especially as quite large groups of pixels seem to be processed in ”lock-step” mode on current


9

Figure 6 Comparison of RGBα volume raycasting using a pre-integrated transfer function (left), spectral volume raycasting based on the Kubelka-Munk model (middle), and spectral volume raycasting including refractive effects like dispersion (right).

GPUs, but also as the maximum number of fragment program instructions can be easily exceeded in that case. As a solution we have limited the maximal number of total reflections that are allowed when traversing a ray through the volume to a total of five. This effectively contains the number of iterations in the sampling loop and in fact does not introduce significant error, since multiple reflections are not easily understood by the observer and in the case of volume rendering typically lead to noisy images. In this section it becomes also clear why we use a rather costly test for checking whether a ray has left the volume. It would be much cheaper to just render the front and back facing polygons of the volume’s bounding box separately as it is done, for example, in [22, 15] and from that information compute the number of steps that are necessary to sample the volume. This obviously only works for linear rays, where the exit point and direction is known in advance, which is not the case for refractive volumes. 8 Results and Discussion In this section we present results of our spectral refractive volume raycasting algorithm, as well as a comparison to standard volume rendering with pre-integrated transfer functions and a CPU-based implementation. Again all presented performance numbers were measured on a standard desktop PC equipped with a 3.4GHz Pentium 4 processor and a NVIDIA GeForce 7800 GTX based graphics accelerator card. Fig. 6 shows three renderings of a 2562 ×225 CT data set using the same basic rendering parameters, but three different raycasting techniques. The leftmost image shows the result of standard pre-integrated volume rendering. The transfer function in this case was derived from the spectral material properties chosen for generating the other two images. Note that only the color components of the RGBα transfer function were automatically calculated, while the opacity values were chosen to match the appearance of the spectral

renderings as good as possible. The other two images were rendered based on the Kubelka-Munk model using 31 samples of the spectrum of visible light. In contrast to the rightmost image the middle image does not include refraction and its resultant effects. But distortion of the background image as well as the inner structures of the volume and dispersion of the light are clearly visible in the last image. Both effects significantly help to depict the three-dimensional nature of the data set. Parts of the volume that were hardly visible in the first two images due to their high transparency can now be clearly identified. On the other hand, artifacts due to noise in the data set are much more apparent while rendering with refraction as it is easy to notice in the left part of the image. What can also be noticed is that the selective absorption acts as an additional depth cue. This is even more obvious in the images shown in Fig. 7. Note how the red color component of the yellow structures is filtered out by the transparent blue medium and their apparent color changes over green to blue depending on the depth of the structure. In Tab. 1 performance numbers for both data sets and all rendering techniques are presented. Since the data set has to be sampled 31 times for the spectral raycasting, the rendering times are significantly larger than for the RGBα method. Given the higher complexity of the Kubelka-Munk compositing as well as the overhead introduced by combining the 31 intermediate results, the performance obtained for the spectral rendering seems to be reasonable and even better than we anticipated. Additionally, including refraction poses a high computational effort for each sampling step, therefore the rendering times for these effects are significantly higher. Nevertheless, the performance that can be achieved using the GPU is up to two orders of magnitude higher than what can be obtained by a software solution. Since comparing a GPU-based with a CPU-based algorithm is always an intricate task we tried to include the same optimizations in both solutions, e.g. a pre-computed look-up table for reflectance and transmittance values.

10


Table 1 Performance results for rendering the volume data sets shown in Figs. 6 and 7. Timings are given in seconds per frame.

Head Convection

GPU-based RGBα , w/o refraction

GPU-based spectral, w/o refraction

GPU-based spectral w/ refraction

CPU-based spectral w/ refraction

0.3 0.1

11.8 5.5

18.0 9.8

2109.7 960.6

9 Conclusion In this paper we have presented an extension to our framework [40] for hardware-accelerated visualization of volumetric data sets. Our system is based on a single-pass raycasting approach taking advantage of recent advances in programmable graphics hardware. Compared to standard slicebased volume rendering techniques, the system exhibits very high flexibility and allows for an easy integration of nontrivial volume rendering techniques which has been demonstrated for a spectral volume rendering approach based on the Kubelka-Munk model. At the moment spectral volume raycasting is far from being interactive, but allows to create images of higher optical correctness than it is possible using the traditional RGBα volume rendering model. Neverthe-

Figure 7 Spectral volume renderings including refractive effects for two different sets of material properties. The data set is a 2562 ×64 simulation of the convection flow in the earth’s crust.

speedup 117.3 98.5

less, for typical volume data sets the presented GPU-based raycasting implementation of the spectral volume raycasting achieves up to two orders of magnitude higher performance compared to a equivalent software implementation. Acknowledgements This work was partly supported by the German Research Council (DFG) as part of SFB 382. The data sets used are courtesy of UNC at Chapel Hill (Head) and the Department of Radiology, University of Iowa (Foot).

References 1. Abdul-Rahman, A., Chen, M.: Spectral Volume Rendering Based on the Kubelka-Munk Theory. Computer Graphics Forum 24(3), 413–422 (2005) 2. Baxter, W., Wendt, J., Lin, M.C.: IMPaSTo: A Realistic, Interactive Model for Paint. In: Proceedings of the 3rd International Symposium on Non-photorealistic Animation and Rendering, pp. 45–148 (2004) 3. Bergner, S., Möller, T., Drew, M.S., Finlayson, G.D.: Interactive spectral volume rendering. In: Proceedings of IEEE Visualization, pp. 101–108 (2002) 4. Bergner, S., Möller, T., Tory, M., Drew, M.: A Practical Approach to Spectral Volume Rendering. IEEE Transactions on Visualization and Computer Graphics 11(2), 207–216 (2005) 5. Boada, I., Navazo, I., Scopigno, R.: Multiresolution Volume Visualization with a Texture-based Octree. The Visual Computer 17(3), 185–197 (2001) 6. Cabral, B., Cam, N., Foran, J.: Accelerated Volume Rendering and Tomographic Reconstruction using Texture Mapping Hardware. In: Proceedings of the 1994 Symposium on Volume Visualization, pp. 91–98 (1994) 7. Cullip, T.J., Neumann, U.: Accelerating Volume Reconstruction With 3D Texture Hardware. Tech. Rep. TR93-027, University of North Carolina at Chapel Hill (1993) 8. Curtis, C.J., Anderson, S.E., Seims, J.E., Fleischer, K.W., Salesin, D.H.: Computer-generated Watercolor. In: Proceedings of the 24th Annual ACM SIGGRAPH Conference, pp. 421–430 (1997) 9. Dorsey, J., Hanrahan, P.: Modeling and Rendering of Metallic Patinas. In: Proceedings of the 23rd Annual ACM SIGGRAPH Conference, pp. 387–396 (1996) 10. Engel, K., Kraus, M., Ertl, T.: High-Quality Pre-Integrated Volume Rendering Using Hardware-Accelerated Pixel Shading. In: Eurographics / SIGGRAPH Workshop on Graphics Hardware ’01, pp. 9–16 (2001) 11. Glassner, A.S.: How to Derive a Spectrum from an RGB Triplet. IEEE Computer Graphics and Applications pp. 95–99 (1989) 12. Glassner, A.S.: Principles of Digital Image Synthesis. Morgan Kaufmann, San Francisco (1995) 13. Guthe, S., Wand, M., Gonser, J., Straßer, W.: Interactive Rendering of Large Volume Data Sets. In: Proceedings of IEEE Visualization ’02, pp. 53–60. IEEE Computer Society (2002) 14. Haase, C., Meyer, G.: Modeling Pigmented Materials for Realistic Image Synthesis. ACM Transactions on Graphics 11(4), 305–335 (1992) 15. Hadwiger, M., Sigg, C., Scharsach, H., Bhler, K., Gross, M.: RealTime Ray-Casting and Advanced Shading of Discrete Isosurfaces. In: Proceedings of Eurographics ’05, pp. 303–312 (2005)


16. Hall, R.: Illumination and Color in Computer Generated Imagery. Springer-Verlag (1989) 17. Hall, R.A., Greenburg, D.P.: A Testbed for Realistic Image Synthesis. IEEE Computer Graphics and Applications 3(8), 10–20 (1983) 18. Johnson, G.M., Fairchild, M.D.: Full-Spectral Color Calculations in Realistic Image Synthesis. IEEE Computer Graphics and Applications 19(4), 47–53 (1999) 19. Klein, T., Strengert, M., Stegmaier, S., Ertl, T.: Exploiting Frameto-Frame Coherence for Accelerating High-Quality Volume Raycasting on Graphics Hardware. In: Proceedings of IEEE Visualization ’05, pp. 223–230. IEEE (2005) 20. Kniss, J., Kindlmann, G., Hansen, C.: Multidimensional Transfer Functions for Interactive Volume Rendering. IEEE Transactions on Visualization and Computer Graphics 8(3), 270–285 (2002) 21. Kniss, J., Premoze, S., Hansen, C., Ebert, D.: Interactive Translucent Volume Rendering and Procedural Modeling. In: Proceedings of IEEE Visualization ’02, pp. 109–116 (2002) 22. Krüger, J., Westermann, R.: Acceleration Techniques for GPUbased Volume Rendering. In: Proceedings of IEEE Visualization ’03, pp. 287–292 (2003) 23. Kubelka, P.: New contributions to the optics of intensely lightscattering materials, Part II. Nonhomogeneous layers. Journal of the Optical Society of America 44, 330–355 (1954) 24. Kubelka, P., Munk, F.: Ein Beitrag zur Optik der Farbanstriche. Zeitschrift für Technische Physik 12, 593–601 (1931) 25. LaMar, E., Hamann, B., Joy, K.I.: Multiresolution Techniques for interactive Texture-based Volume Visualization. In: Proceedings of IEEE Visualization ’99, pp. 355–361 (1999) 26. Li, W., Mueller, K., Kaufman, A.: Empty Space Skipping and Occlusion Clipping for Texture-based Volume Rendering. In: Proceedings of IEEE Visualization ’03, pp. 317–324 (2003) 27. Max, N.: Optical Models for Direct Volume Rendering. IEEE Transactions on Visualization and Computer Graphics 1(2), 99– 108 (1995) 28. Meyer, G.W.: Wavelength Selection for Synthetic Image Generation. Computer Vision, Graphics and Image Processing 41, 57–79 (1988) 29. Noordmans, H., van der Voort, H., Smeulders, A.: Spectral Volume Rendering. IEEE Transactions on Visualization and Computer Graphics 6(3), 196–207 (2000) 30. NVIDIA Corporation: NVIDIA SDK 8.0. http://developer.nvidia.com/object/sdk home.html (2004) 31. Peercy, M.S.: Linear Color Representations for Full Spectral Rendering. ACM SIGGRAPH Computer Graphics 27(3), 191–198 (1993) 32. Rezk-Salama, C., Engel, K., Bauer, M., Greiner, G., Ertl, T.: Interactive Volume Rendering on Standard PC Graphics Hardware using Multi-textures and Multi-stage Rasterization. In: Proceedings of the ACM SIGGRAPH/EUROGRAPHICS Workshop on Graphics Hardware, pp. 109–118 (2000) 33. Rodgman, D., Chen, M.: Refraction in Discrete Raytracing. In: Proceedings of Volume Graphics 2001, pp. 3–17. Springer, New York (2001) 34. Rodgman, D., Chen, M.: On Volume Denoising Filters for Visualizing Refraction. In: Proceedings of Vision, Modeling, and Visualization ’03, pp. 355–363 (2003) 35. Röttger, S., Guthe, S., Weiskopf, D., Ertl, T.: Smart HardwareAccelerated Volume Rendering. In: Proceedings of EG/IEEE TCVG Symposium on Visualization VisSym ’03, pp. 231–238 (2003) 36. Rougeron, G., Péroche, B.: Color Fidelity in Computer Graphics: A Survey. Computer Graphics Forum 17(1), 3–15 (1998) 37. Rudolf, D., Mould, D., Neufeld, E.: Simulating Wax Crayons. In: Proceedings of the 11th Pacific Conference on Computer Graphics and Applications, pp. 163–172 (2003) 38. Schlick, C.: A Survey of Shading and Reflectance Models. Computer Graphics Forum 13(2), 121–131 (1994) 39. Schuster, A.: Radiation Through a Foggy Atmosphere. Journal of Astrophysics 21, 1–22 (1905)

11

40. Stegmaier, S., Strengert, M., Klein, T., Ertl, T.: A Simple and Flexible Volume Rendering Framework for Graphics-Hardware–based Raycasting. In: Proceedings of the International Workshop on Volume Graphics ’05, pp. 187–195 (2005) 41. Weiler, M., Kraus, M., Merz, M., Ertl, T.: Hardware-Based Ray Casting for Tetrahedral Meshes. In: Proceedings of IEEE Visualization ’03, pp. 333–340 (2003) 42. Westermann, R., Ertl, T.: Efficiently using Graphics Hardware in Volume Rendering Applications. In: SIGGRAPH ’98: Proceedings of the 25th annual Conference on Computer Graphics and interactive Techniques, pp. 169–177 (1998) 43. Wyszecki, G., Stiles, W.S.: Color Science Concepts and Methods, Quantitative Data and Formulae. Wiley-Interscience Publication (1982)