Curvature-based transfer functions for direct volume rendering: methods and applications

Direct volume rendering of scalar fields uses a transfer function to map locally measured data properties to opacities and colors. The domain of the transfer function is typically the one-dimensional space of scalar data values. This paper advances the use of curvature information in multi-dimensional transfer functions, with a methodology for computing high-quality curvature measurements. The proposed methodology combines an implicit formulation of curvature with convolution-based reconstruction of the field. We give concrete guidelines for implementing the methodology, and illustrate the importance of choosing accurate filters for computing derivatives with convolution. Curvature-based transfer functions are shown to extend the expressivity and utility of volume rendering through contributions in three different application areas: nonphotorealistic volume rendering, surface smoothing via anisotropic diffusion, and visualization of isosurface uncertainty.

Curvature-Based Transfer Functions for Direct Volume Rendering: Methods and Applications Gordon Kindlmann 1 Ross Whitaker 1 Tolga Tasdizen 1 Torsten Moller2 1 Scientific Computing and Imaging Institute, University of Utah 2 Graphics, Usability, and Visualization (GrUVi) Lab, Simon Fraser University Figure I: Volume renderings of a 64-1 synthetic volume with four different curvature measures. Left to right: first principal curvature ffj, second principal curvature ff2, mean curvature (ffj + ff2)/2, and Gaussian curvature k1k2. Magenta indicates negative curvature, green indicates positive. Iso-curvature contours arc in black, except for zero curvature in blue. Abstract Direct volume rendering of scalar fields uses a transfer function to map locally measured data properties to opacities and colors. The domain of the transfer function is typically the one-dimensional space of scalar data values. This paper advances the use of cur­vature information in multi-dimensional transfer functions, with a methodology for computing high-quality curvature measurements. The proposed methodology combines an implicit formulation of curvature with convolution-based reconstruction of the field. We give concrete guidelines for implementing the methodology, and il­lustrate the importance of choosing accurate filters for computing derivatives with convolution. Curvature-based transfer functions arc shown to extend the expressivity and utility of volume render­ing through contributions in three different application areas: non- photorcalistic volume rendering, surface smoothing via anisotropic diffusion, and visualization of isosurfacc uncertainty. CR Categories: 1.3.7 [Computing Methodologies]: Computer Graphics-Three-Dimensional Graphics and Realism 1.4.7 [Com­puting Methodologies]: Image Processing and Computer Vision- Feature Measurement Keywords: volume rendering, implicit surface curvature, convolution-based differentiation, non-photorcalistic rendering, surface processing, uncertainty visualization, flowlinc curvature 1 {gk|whitaker|tolga}®cs.utah.edu "torsten®cs.sfu.ca IEEE Visualization 2003, October 19-24, 2003, Seattle, Washington, USA 0-7803-8120-3/03/$ 17.00 ©2003 IEEE 1 Introduction Direct volume rendering depicts structure in scalar fields through a simple combination of mappings. At each rendered sample point, locally measured numerical quantities (mainly the data value itself) arc mapped via the transfer function to optical quantities such as opacity and color. Basic computer graphics routines can then shade, composite, and project the samples into a coherent visualization. In nearly all applications, some part of this process involves the first-order differential structure of the field: the gradient. Shad­ing routines generally require the gradient direction as input to a surface-based lighting model. There is also a trend towards using gradient magnitude in the transfer function domain, enabling vi­sualization of material boundaries based on data value and edge strength [Lcvoy 1988], In cither case, first-order differential field structure enhances the clarity and effectiveness of direct volume rendering. This paper further promotes the use of differential struc­ture in direct volume rendering, by demonstrating a family of trans­fer functions based on curvature, using second-order derivatives. The theory behind curvature measurement in volume visualiza­tion draws from two distinct fields of study. The first is differ­ential geometry. We compute isosurfacc curvature directly from the volume data using a simple algebraic framework of differen­tial invariants. Additional insight into the underlying field structure comes from geometry orthogonal to the isosurfaccs, which is also described by an invariant. The second field of study guiding this work is signal processing and filter design. Numerical differentia­tion accentuates noise, and the problem increases with the order of the derivative. Care must bc taken to measure second derivatives ro­bustly and accurately. We describe a convolution-based derivative measurement scheme which allows control of the tradeoff between accuracy, continuity, and filter size. This curvature measurement method obviates intermediate surface construction, surface param­eterization, and explicit finite difference schemes. The ability to compute reliable curvature information at arbitrary points in sam­pled volume data creates new opportunities for direct volume ren­dering with curvature-based transfer functions. This paper selects different components of curvature information to serve as domain variables in multi-dimensional transfer func- 513 Proceedings of the 14th IEEE Visualization Conference (VIS'03) 0-7695-2030-8/03 $ 17.00 © 2003 IEEE tions. We demonstrate that curvature-based transfer functions en­hance the expressive and informative power of direct volume ren­dering with contributions in three different application areas. Non- photorealistic volume rendering is improved in Section 5 with con­tours exhibiting constant thickness in image space, and a flexible method of emphasizing ridge and valley surface creases. Section 6 shows how volume rendering with total curvature provides a tech­nique for studying the behavior of volumetric surface processing algorithms. Section 7 uses flowline curvature as a way to quali­tatively indicate regions of high geometric uncertainty in surface models created by isosurfacing. 2 Previous Work This paper builds on a significant amount of previous research in volumetric curvature measurement and its application to volume vi­sualization. Mongaetal. [1992; 1994] use convolution with Gaus­sian kernels, as well as their first, second, and third derivatives to compute surface curvature characteristics, and to find ridge lines (extrema in normal curvature along the principal curvature direc­tion). Interrante et al. [1995] use ridge and valley lines to sim­plify depiction of the skin surface for surgical planning. Their later work [1997; 1997] systematically explores a variety of methods for using strokes or textures aligned with principal curvature directions to effectively convey the configuration of semi-transparent isosur­faces relative to internal structures. Kinetic visualization [Lum etal. 2003] uses particle motion along the principal curvature direction to clarify depiction of surface orientation. In the work that most di­rectly influences our own, Hladuvka et al. [2000] describe volume renderings with two-dimensional transfer functions using a sim­plified space of principal curvatures, and demonstrate that surface shape is an intuitive space in which to specify color and opacity. Another thread of previous work adapts non-photorealistic ren­dering (NPR) techniques to direct volume rendering [Gooch and Gooch 2001; Ma et al. 2002], Treavett and Chen [2000] create pen-and-ink style visualizations of volume data with a combina­tion of volume-space and image-space measurements, including a simple curvature measure. Rheingans and Ebert [2001] modify the standard volume rendering pipeline by applying a variety of NPR techniques to enhance interesting features and regions, and to provide depth and orientation cues. Csebfalvi et al. [2001] use gradient information to create interactive visualizations emphasiz­ing object contours, quickly revealing internal structures with lit­tle parameter tuning. Lu et al. [2002] mimic stipple drawing by using a variety of feature detection methods to place point primi­tives in the volume with carefully controlled density and location. Lum and Ma [2002] use commodity graphics hardware to create richly detailed renderings of volume datasets using a combination of NPR methods. Nagy et al. [2002] use graphics hardware to en­rich interactive volume renderings with artistic shading methods and curvature-directed strokes and hatches. 3 Curvature Measurement The contribution of this section is a self-contained derivation and intuitive explanation of our curvature measurement method, ending with the three simple steps required to put it into practice. The steps can be implemented without perfect understanding of the underly­ing mathematics. Following Monga etal. [1992], we rely on convo­lution with continuous filters, rather than an explicitly constructed coordinate frame aligned with the gradient direction [Interrante 1997; Nagy et al. 2002], We feel our convolution-based approach is significantly simpler than a previous method based on charac­terizing multiple planar curves passing through the point of inter­est [Hladuvka et al. 2000], Though visualization results may bene­fit from a pre-process smoothing of the data, no pre-computation is required for the curvature measurement itself. The curvature of a surface is defined by the relationship be­tween small positional changes on the surface, and the resulting changes in the surface normal. In sampled volume data, surfaces are implicitly represented as isosurfaces of reconstructed contin­uous data values /(x). Assuming that the values of / increase as we move further inside objects of interest (e.g., a standard CT scan), the surface normal is defined as n = - g/|g|, with the gradi­ent g = Vf = • Curvature information is contained in VnT, a 3x 3 matrix. However, we do not want to evaluate the gradi­ent of a (pre-computed) normalized vector, since this hinders direct convolution-based measurement of the original data. Expanding upon the derivation in [Mitchell and Hanrahan 1992]: Vn1 -V fr Vlsl -±(« VgT gVTlgl |g| Igl2 gVTigTg) = -n[H- |g| gUgrH) 2|g|2 = -Ti-iI-nnT)H. |g| ' I is the identity matrix, and H is the Hessian matrix: gVTigTg) 2lgl(gTg)l/2 H H d2f/d.x2 d2f/dxdv d2f/dxdz d2f/dxdv d2f/dv2 d2f/dvdz d2f/dxdz d2f/dvdz d2f/dz2 (1) The outer product of n with itself, nnT, is a linear operator that projects onto the one-dimensional span of n. I - nnT projects onto the orthogonal complement of the span of n, namely, the tangent plane to the isosurface. Letting P = I - nnT, then t 1 Vn = --PH. (2) Reading Equation 2 from right to left permits some intuitive un­derstanding of VnT. Vector calculus tells us that the Hessian matrix H represents how the gradient g changes as a function of infinites­imal changes of position in K3 [Marsden and Tromba 1996], The changes in g have a component along g (the gradient can change length), and a component within the tangent plane (the gradient can change direction). For the purposes of describing curvature, only the latter component matters. It can be isolated with left- multiplieation by P. Finally, the - l/|g| scaling factor converts in­finitesimal changes of the (un-normalized) gradient g into infinites­imal changes of the unit-length normal n. Both P and H are symmetric matrices, but in general, VnT is not symmetric. However, if v lies in the tangent plane, then Pv = v and v P = vT, so for u and v in the tangent plane. vtPHu vtHu utHv = utPHv . That is, the restriction of VnT = - PH/1g| to the tangent plane is symmetric, and thus there exists an orthonormal basis {pj.p,} for the tangent plane in which VnT is a 2 x 2 diagonal matrix [Hoffman and Kunze 1971], This basis can be easily extended to an orthonor­mal basis for all of K3, {pj.p,.!!}. In this basis, the derivative of the surface normal is Vn1 Kj 0 al 0 KS CT, 0 0 0‘ The bottom row is all zero because no change in position can make the normal n change in length. Motion within the tangent 514 Proceedings of the 14th IEEE Visualization Conference (VIS'03) 0-7695-2030-8/03 $ 17.00 © 2003 IEEE plane, along p, and p,, leads to changes of n along the same di­rections, with proportionalities ifj and k? respectively. There is no cross-term, or twisting of the normal, by the choice of {p, ,p, }. By the definition of surface curvature [Do Carmo 19761, pf and p, are the principal curvature directions, while K-, and K-, are the principal curvatures. As one moves along the normal, off of the surface or deeper into it, the normal tilts according to cr, and cr,. This aspect of implicit surface curvature, termed,fiowtine curvature by ter Haar Romeny et al. [19911, is explored further in Section 7. The above derivation of VnT is not novel. Others have exploited the fact that it has eigenvalues k-,, k-,, and 0 to compute curvature of implicit surfaces [Belyaev et al. 1998; Preusser and Rumpf 20021. Multiplying VnT by P has the effect of isolating k-, and K-, in the {p, ,p,,n} basis: ' 1 0 0 G = VnTP = VnT 0 1 0 = 0 0 0 Our surface curvature measurements are based on G, which we term the geometry tensor. In practice, G will be known only in the (X,Y,Z) basis of the volume axes, and will not have the readily transparent form of Equation 3. Matrix invariants provide the lever­age to extract the desired curvature values k:, and k? from G, re­gardless of the principal curvature direction coordinate frame. The trace of G is k-, + k-,. The Frobenius norm of G, notated |G|F and defined as Vtrace(GGT), is '/icf + &. k*| and K-, are then found with the quadratic formula. To summarize, here are the steps needed to compute curvature at an arbitrary point in a scalar field: 1. Measure the first partial derivatives comprising the gradient g. Compute n = -g/|g|, and P = I - nnT. 2. Measure the second partial derivatives comprising the Hessian H (Equation 1). Compute G = -PHP/|g|. 3. Compute the trace T and Frobenius norm F of G. Then, T + V2F2-T2 T- s/2F2-T2 K> - 2 ' ^2 ~~ 2 ' If the data values inside regions of interest are lower than the background (e.g., inverted from a standard CT scan), the only change in the formulation is the sign of the geometry tensor: G = PHP/|g|. Though not used in this paper, the principal cur­vature directions are easily found as eigenvectors of G. The most important task in the curvature computation is measuring all the necessary partial derivatives (in g and H). It is for this task that we use convolution with continuous filters. 4 Measuring Derivatives with Convolution Image-order volume rendering of discretely sampled data relies on convolution with continuous filters to reconstruct values at arbi­trary locations, such as sample points along a ray. We perform three-dimensional reconstruction as a separable product of one­dimensional convolutions along the X,Y,7 axes. The combination of filters used on each axis determines whether the reconstructed quantity is the interpolated data value, or one of its many partial derivatives. Measuring d'f/dxdy, for example, is done with first- order derivative filters on the X and Y axes, and with a zero-order derivative (interpolation) filter on the Z axis. Accurate curvature measurement depends on a judiciously chosen combination of fil­ters for zero, first, and second derivatives. For reasons of practical­ity and efficiency, we use piecewise polynomial filters with sym­metric support. We apply the filter design framework of Moller et al. [19981 to inform the choice of filters. It is the first time the framework has been used to create second derivative filters. We review the spatial filter design framework by considering a continuous function fit) in one dimension. Sampling fit) at unit intervals produces a sequence of discrete values, notated f[k], where k is an integer. From f[k] we must approximate both the original function fit), as well as its first and second derivatives fit) and fit). Each of the three approximations is done by con­volution with a continuous filter. The result of convolving f[k] with continuous filter vi'(r) is notated fit). By combining the IVth or­der Taylor expansion of f with the convolution definition of f" (r), Moller etal. [19981 derive: fit) = </(r)+«f/'(r)+«f/"(r) + ... + a%flN)(t) +?N <4) The a* are termed error coefficients. From a formulaic standpoint, a* encapsulates a convolution of vr with the fi-th term of the sam­pled Taylor expansion of /. From a practical standpoint, however, the vector of error coeffi­cients is the means by which we "dial up" a filter vt\ according to the desired derivative and accuracy properties. We describe both of these in turn. Looking at Equation 4, we can see that if a filter vr does a perfect job of reconstructing fit), then fit) = fit), so £/q = 1 and all other a)*' (as well as the remainder term rJJ) are zero. Similarly, for ideal first (or second) derivative filters, af (or a*') is unity, and all other coefficients are zero. The measurement behav­ior of a filter vr is governed by its first non-zero error coefficient. To specify that vr should be a rf-order derivative filter, we require a* = 0 for all n < d, and set a^j = 1. In spatial filter design, filter quality is quantified in terms of ac­curacy: the more accurate a filter is, the higher degree polynomial it can reconstruct exactly. A kth degree error filter (&-EF) filter is one which can exactly reconstruct polynomials of degree k- 1 or lower. In the context of Equation 4, we gain accuracy by requir­ing one or more of the error coefficients beyond a^j to be zero. To specify that vr should be a -derivative &-EF filter, we let /V = d + k and require a* = 0 for all n : d < n < N. Note that the number of constraints we set (that is, the number of a)*' we specify) is the sum of the derivative order and the accuracy level. This also turns out to be the lower bound on the filter support (the number of data value samples required in the convolution sum). Continuity is also a desirable filter property, since it avoids visual artifacts in reconstruction, and it helps reduce post­aliasing [Marschner and Lobb 19941. A continuous filter is CM if it has M continuous derivatives. It can be shown that infinitely many piecewise polynomial filters satisfy any given derivative, ac­curacy, and continuity requirements [Moller 19991. The filter wit) is ultimately implemented as a list of coefficients of piecewise poly­nomials in t. The coefficients are found as the solution of a linear system embodying the derivative, accuracy, and continuity proper­ties described above [Moller 19991. Spatial filter design supports the creation and characterization of second derivative filters, but unfortunately curvature is not sim­ply a second derivative. Rather, Section 3 showed how curvature is extracted from both first and second derivatives, through alge­braic transforms. We have conducted an empirical study to better understand how filter quality affects curvature measurements. Our study used measurements of mean curvature (if, + k',)/2 on the Marschner-Lobb [19941 dataset, a 403 sample floating-point vol­ume with isosurfaces resembling a graph of the sine function. Fig­ures 2 shows the analytical form of the data, including colormap- ping with correct mean curvature values. Figure 3 shows results from a variety of filter combinations. Fig­ure 3(a) uses the Catmull-Rom cubic spline and its first and sec­ond derivatives. Because the Catmull-Rom spline is C1 but not C2, 515 Proceedings of the 14th IEEE Visualization Conference (VIS'03) 0-7695-2030-8/03 $ 17.00 © 2003 IEEE 4- (h'| 4- K-,j/2 Figure 2: Marschner-Lobb dataset with ideal reconstruction. In (b), colors indicate mean curvature. its second derivative is discontinuous, leading to markedly discon­tinuous curvature measurements. This example demonstrates how well-regarded reconstruction filters arc not always the basis of use­ful derivative reconstruction filters, especially in the context of sec­ond derivatives. Figure 3(b), on the other hand, better captures the shape and curvature on the thin ridges, but at the expense of using an 83 sample support over which the convolutions arc performed. Large filter supports arc costly, especially in three dimensions, be­cause of poor memory locality. Finding filters which give good results with reasonable supports is an ongoing research problem. Seeking to provide practical advice on a filter set for efficient cur­vature measurement, we performed a search among the filters (cre­ated by the framework above) with 43 sample support. Allowing up to 3rd order continuity (C3) and 4th order accuracy (4-EF), a total of 200 filter combinations were assessed by rendering the Marschncr- Lobb dataset. These renderings (not shown) used an orthonormal projection from directly above the dataset, with shading turned off, to generate grayscale images of mean curvature. Based on the dif­ference between the ideal and rendered curvatures, an "optimal" 43 filter set was found; sec Figure 3(c). An interesting property of this filter set is that for the derivatives, high accuracy is more important than high continuity. For example, there arc 4-support 2nd deriva­tive filters with higher continuity than C°, but they arc 1-EF, not 2-EF. Figure 3(d) uses another set of filters, also created by the frame­work above, but better known through the BC-splincs of Mitchell and Nctravali [1988]. Interpolation is done by the Catmull-Rom fil­ter (B,C = 0,0.5), which is the only 3-EF (the most accurate) inter­polating BC-splinc. The first and second derivatives arc measured with the first and second derivatives of the cubic B-splinc (B,C = 1,0), which is the only C2 (the most smooth) BC-splinc. In our quantitative comparison of 200 filter sets, these did nearly as well as the "optimal" 43 filters. Given that any 43 sample filter set will bc quite limited in the feature detail it can reconstruct and mea­sure, we recommend that when practical considerations dictate a 43 support, the filters listed in Figure 3(d) arc sufficient; only slightly better results will come from filters with higher polynomial degree. When larger filter supports can bc employed, higher accuracy is more important than higher continuity. In most practical appli­cations, the underlying data can bc sufficiently modeled with low degree polynomials. 3-EF or 4-EF filters should bc sufficient; those in Figure 3(b) arc excellent. When measuring derivatives in noisy data, smoothing often improves accuracy. Such smoothing can bc a pre-process. Or, it can be folded into the convolution, by slightly blurring values, instead of interpolating them. For instance, instead of using the Catmull-Rom filter as in Figure 3(d), we can use the B-splinc. We have found that the resulting loss of accuracy on very small features is compensated for by the reduction of spurious cur­vature variations due to noise. For the remainder of this paper, we use the B-splinc and its derivatives for our curvature measurements. (a) / : Catmull-Rom f : I'' deriv of Catmull-Rom /" : 2nd deriv of Catmull-Rom (degree 6) (degree 6) D, CJ, 4-EF (degree 5) (C) / : 0-D. c-. 3-EF (degree 7) f : 1 -D. C2. 2-EF {degree 4) f" : 2D. (''.2 IT (degree]) (d) / : Catmull-Rom f : I'' deriv ofB-spline /" : 2'"l deriv of B-spline Figure 3: Marschner-Lobb dataset with various filtered reconstruc­tions. Below each image is a list of the three filters used to re­construct the data value and the first and second derivatives. The polynomial degrees arc listed for (b) and (c). In (a) and (d) the polynomial degrees arc 3, 2, and 1. 5 Non-Photorealistic Volume Rendering 5.1 Controlling Contour Thickness Contours (sometimes referred to as silhouettes) arc a basic part of all current approaches to non-photorcalistic volume rendering. By emphasizing the transition between front-facing and back-facing surface locations, contours delineate object shape and clarify sites of occlusion. Contours of polygonal surfaces arc drawn on edges between faces with different visibilities, as determined by a sign change in v • n, the dot product between the view vcctor v and the face normal n [Gooch and Gooch 2001]. Lacking polygonal dis­cretization, volume rendered contours arc generally produced with a univariate function of v- n. For instance, the surface color is dark­ened when v- n is within some user-determined range around zero. Figure 4: Contours based solely on v • n 516 Proceedings of the 14th IEEE Visualization Conference (VIS'03) 0-7695-2030-8/03 $ 17.00 © 2003 IEEE An unfortunate consequence of this approach is uncontrolled variation in the apparent contour thickness. As shown in Figure 4, where the surface is nearly flat, a large region of surface normals is nearly perpendicular to the view vector, making the contours too thick. Conversely, in fine structures, where the emphasis provided by contours could be especially helpful, they can be too thin. We propose regulating the thickness of volume rendered con­tours by using curvature. Specifically, the normal curvature along the view direction, kx. restricts the range of surface normal ori­entations rendered as contour. We find ks, in terms of the tangent projection P and the geometry tensor G, exploiting the fact that re­stricted to the tangent plane, G is exactly the Weingarten map [Do Carmo 19761: (Pv)T Pv |Pv| |Pv| 1 Gv Gv IPvl- vTPv Figure 5: Creating contour thickness T based on surface normal n, view vector v, and radius of surface curvature R = 1/kx. Contours with approximate image space thickness T are ren­dered with a two-dimensional transfer function of ks, and v n. Our method is an approximation because ks, is assumed to be constant over the contour extent. Figure 5 illustrates a slice through a sur­face, on the plane spanned by the surface normal and the view vec­tor as it grazes the surface. Within this plane, m is the single surface normal orthogonal to v. A certain range of surface normals around m must be rendered as contour so that, when projected on the image plane, the contour region appears to have thickness T. This range is a function of T and of R, the radius of curvature at this point along the contour. Letting 0T denote the maximum angle between m and a surface normal within the contour: |m-n| > cos(0T) =$■ |v - n| < sin(0r) =Hv-n < s/R2-(R-T)2 _ jlRT-T2 r2 1(2-1 R\ R =$■ | v - n | < s/T ks,(2 - F k'v (5) Figure 6(a) graphs the inequality in Equation 5; it is simply a portion of a circle. When curvature kx is small, a narrow range of surface orientations should contribute to the contour, while a larger range of orientations must be colored to create the contour on high curvature features. The lookup table in Figure 6(b) is one way to implement the contour transfer function. High values of ks, are clamped to 1/T to ensure that very small features are not missed. The transition between dark and light in the contour trans­fer function may be blurred somewhat to soften the contour edges, and to avoid accentuating minor surface variations near the contour. (a) Contour function (b) 2-D lookup table Figure 6: Graph of contour function, and the lookup table used to implement thickness-controlled contours in practice. Figure 7: Thickness-controlled contours, using 7 = 1 (left) and T = 2.5 (right). Figure 7 demonstrates thickness-controlled contours with two different values for the thickness parameter T. The contour thick­ness around the ears and back of the bunny are now equal, and the surface details on the lower body are more clearly depicted, in con­trast to Figure 4. Also, varying T succeeds in adjusting the image space thickness of the contours. Since the dimensions of T are vox­els, setting 7 = 1 means that contours should be about as thick as the size of one voxel projected onto the image plane. As with pre­vious methods for volume rendered contours, it is possible for our thickness-controlled contours to be drawn at locations where there is actually no change in surface visibility. Also, contours may be drawn improperly in regions where curvature is too low to be mea­sured accurately. We have not found either of these problems to be serious impediments to the illustrative purpose of contours. 5.2 Emphasizing Valleys and Ridges One ingredient missing in the current vocabulary of non- photorealistic volume rendering techniques is the use of curvature magnitudes to enhance depiction of surface details. Conveying overall surface shape is effectively done by using textures aligned with curvature directions [Interrante et al. 1997; Interrante 1997; Girshick et al. 2000; Lum et al. 2003], but we wish to avoid the overhead of pre-computing and storing the necessary volumetric texture. Though visually similar to depiction of ridge and valley lines, our approach is simpler than previous volumetric methods in that we do not use third derivatives [Monga et al. 1994], nor do we perform a local search to determine curvature extrema [Interrante et al. 1995], Rather, we work in the two-dimensional space of prin­cipal curvatures ((CpKs). similar to the approach of Hladuvka et al. [2000], with the important difference that the geometry tensor can distinguish between convex and concave surface geometry. Figure 8 demonstrates transfer functions of the (k} , k-,) principal curvatures. Figure 8(a) portrays curvature space with a sequence of quadratic surface patches (each a 203 volume) rendered with the transfer function shown. Figure 8(b) shows how curvature trans­fer functions can emphasize surface variations not otherwise read- 517 Proceedings of the 14th IEEE Visualization Conference (VIS'03) 0-7695-2030-8/03 $ 17.00 © 2003 IEEE 6 Visualizing Surface Smoothing (a) Volume rendered diagram of (v^v,) space. The colors in the (K'pk',) transfer function domain are mapped onto the patches with corresponding surface curvature. (b) Left: Visualization of ear curvature using transfer function from (a); Right: ridge and valley emphasis implemented with inset transfer function, combined with Gooch shading Figure 8: Principal curvature transfer functions ily seen. Following a technical illustration convention of darken­ing valleys and lightening ridges [Gooch et al. 1999], we darken or lighten volume rendered samples according to their location in (, k2 ) space, using the transfer function shown in Figure 8(h). In addition to the curvature-based transfer function, Gooch shading is used to convey overall shape [Gooch et al. 1998]. Figure 9 highlights the complementary nature of our two curvature-based non-photorealistie volume rendering effects. Thickness-controlled contours delineate the shapes of large struc­tures, while ridges and valleys highlight tine surface details. In these renderings, surface color is assigned as a function of three vectors (v, n, 1: view, surface normal, light) and three different kinds of curvature (>'», K'(, k2: view direction, principal curvatures). We follow the example of previous work demonstrating the benefit of expressing multi-dimensional transfer functions as the separable product of simpler functions [Kniss et al. 2002]: RGB(v,n,l, kv, K,, k2) = g(n-l) c(k-v, v-n) e(ic,, k2) . Gooch shading comes from g(n-l); e(KV,n-v) creates thickness- controlled contours; e(fq, k2) provides ridge and valley emphasis. The contour thickness T is the only parameter to tune in c( K'v, v • n). The ridge and valley emphasis does require subjective adjustments, but the intuitive nature of (k'(, k-,) space makes this relatively sim­ple. By far the most challenging parameter setting remains the opacity assignment, for which we use a two-dimensional transfer function of data value and gradient magnitude, as guided by previ­ous work [Kindlmann and Durkin 1998; Kniss et al. 2002], Medical imaging and three-dimensional range sensing technologies produce high resolution descriptions of object surfaces, but they of­ten suffer from noise. Level set methods for denoising the data have become popular because of their inherent ability to handle com­plex topology [Osher and Sethian 1988], Level sets can be evolved by minimizing an energy function based on principal curvatures, the lowest-order differential invariants of surface shape [Tasdizen et al. 2002]. Minimizing the surface integral of total cunvtture kt = \/k\ + k; gives an isotropic surface smoothing, analogous to Gaussian image blurring. A variant of this energy function can be used to preserve prominent features on the surface while elimi­nating noise [Tasdizen et al. 2002], generalizing anisotropic image diffusion [Perona and Malik 1990] to surface processing. We propose that curvature-based direct volume rendering is a natural and effective tool for visualizing curvatures and surfaces as they evolve during a level set solution. Level sets and direct vol­ume rendering share an implicit surface representation on a com­mon underlying grid. Because curvature-based volume rendering requires no pre-processing, it permits direct inspection of interme­diate smoothing results, without an intervening meshing step. To visualize the spatial and temporal structure of total curvature on the evolving surface, we use two one-dimensional transfer functions. One assigns surface color as a function of total curvature; the other assigns opacity as a function of data value. Figure 10: Six snapshots of anisotropic surface smoothing Figure 10 shows volume renderings of six instances during the anisotropic smoothing of an isosurface from MRI data. All ren­dering parameters are fixed; only the volume dataset is changing. Noise in the original surface (top-left picture) is seen as randomly distributed high curvature points on the surface. As the surface evolves, most of the high curvature points fade away, while the rest coalesce into prominent creases. This is predicted by mathe­matical analysis of the anisotropic diffusion equations [Perona and Malik 1990], Our visualization experiment verifies this property, as well as illuminating other features of the process. For instance, as the surface smoothly evolves, the network of creases undergoes discrete changes in connectivity. An interesting research topic sug­gested by these visualizations is steering the energy function in or­der to create or maintain a particular topology in the crease network. 518 Proceedings of the 14th IEEE Visualization Conference (VIS'03) 0-7695-2030-8/03 $ 17.00 © 2003 IEEE Figure 9: Curvature-based non-photorcalistic volume rendering. The two small images depict the contribution of thickness-controlled con­tours (top) and ridge and valley emphasis (bottom). The large images show the combination of the two effects. 7 Visualizing Isosurface Uncertainty If a grayscale image is thought of as a height field on which a fluid is flowing down ward,,/fovt'/iw cunmure is the curvature of the path of the flow, as viewed from above [ter Haar Romeny et al. 19911. In three dimensions, flowlinc curvature Ky characterizes the degree to which an isosurfacc changes its orientation as a function of small changes in isovalue. With zero flowlinc curvature, adjacent isosur- faccs arc parallel. Like the principal curvatures, flowlinc curvature is a differential invariant. However, instead of describing change within a single isosurfacc, it describes change between isosurfaccs. It can be computed from VnT (Equation 2): Ky = |Vn (nn )|F ; K'l 0 al ' 0 0 0 0 K-, CT, 0 0 0 0 0 0 0 0 1 = \/of + C7£ . To our knowledge, flowlinc curvature has not been used in vol­ume visualization. Wc propose that flowlinc curvature can be used to visualize the uncertainty of material boundaries as modeled by isosurfacc extraction. The physical shape of a boundary between two materials is a fixed and intrinsic property of the object being sampled in the volume dataset. At locations where small changes in isovalue produce large changes in isosurfacc orientation, the iso­surfacc is probably a poor model of the material boundary. On the other hand, where the shape of an isosurfacc is robust against small isovalue changes, wc can be more certain of its geometric accu­racy. By colormapping flowlinc curvature onto volume rendered surfaces, wc can qualitatively indicate uncertainty of surface shape. Regions of especially high flowlinc curvature arc visually flagged as suspect. One common material boundary that can be hard to extract via isosurfacing is the surface of bone in CT data. Figure 11 shows six pscudo-isosurfacc volume renderings of a thumb. Each is ren­dered with a univariate opacity function changing from transparent to fully opaque over a narrow range of data values, effectively cre­ating a threshold at the specified isovalue. Flowlinc curvature is eolormapped onto the surfaces through a single transfer function common to all the renderings. The (anatomically erroneous) holes on the bone surface arc consistently lined by high flowlinc curva­ture, and the sites of hole formation can largely be predicted by flowlinc curvature patterns at lower isovalucs. There arc, however, locations with high surface curvature and low flowline curvature. Figure 11: CT isosurfaccs at values 1055, 1095, 1175, 1255, 1335, and 1405, eolormapped by flowlinc curvature. Circle indicates site of high surface curvature and low flowlinc curvature. such as the (circled) lower edge of the top bone in the second im­age. Because there is no single isovalue for which the resulting isosurfacc has consistently low flowlinc curvature, these visualiza­tions illustrate how hard it is to find a CT isovalue which represents the bone surface with high certainty. 8 Discussion, Future Work Wc have demonstrated that curvature plays a valuable role in cre­ating informative volume visualizations. Our method conforms to a standard post-classification volume rendering pipeline: the color and opacity assignments at cach raycast sample arc generated from transfer functions, on the basis of locally computed quantities, such as gradient and curvature. Rhcingans and Ebert [20011 have de­scribed how this approach engenders the difficult task of transfer function creation and modification. However, as noted in Scc­tion 5.2, wc seek to control this complexity by enforcing separabil­ity of the transfer function into simpler, intuitive components. The conceptual orthogonality of the various NPR effects maps to or­thogonality between the corresponding transfer function domains. Biological surfaces tend to be somewhat noisy, and as a result. 519 Proceedings of the 14th IEEE Visualization Conference (VIS'03) 0-7695-2030-8/03 $ 17.00 © 2003 IEEE where the surface curvature is very low. as on the scalp or the top of the skull, it is hard to measure curvature accurately. Adaptive methods which measure derivatives over a larger area in such cir­cumstances could help, as could proper application of scale-space methods [terHaar Romeny etal. 1991]. As mentioned in Section 4. more analysis is needed to determine the precise relationship be­tween the characteristics of the convolution filters, and the quality of curvature measurements produced. Finally, ourproof-of-concept raycaster does not exploit graphics hardware, and is far from inter­active. Previous work on hardware-based convolution [Hadwiger et al. 2001] and hardware-based NPR [Lum and Ma 2002; Nagy et al. 2002], however, suggests the exciting possibility of high- quality curvature measurements at interactive rates. 9 Acknowledgments We thank the anonymous reviewers for constructive comments. The bunny CT scan in Figures 4 and 7 is courtesy of the Stanford Vol­ume Data Archive. The CT data in Figures 8(b). 9. and 11 is cour­tesy of the Visible Human Project. Funding for this research in­cludes NIH grants 5P20HL68566-02 and 2P41RR12553-04. Every rendered image in this paper can be regenerated exactly with open-source software and public datasets; see Chttp: //www. sci .Utah. edu/~gk/vis03>. References Bfj.yaf.v, A., Pasko, A., and Kunii, T. 1998. Ridges and ravines on implicit surfaces. In Computer Graphics International '98, 530-535. Csf.rfalvi, B., Mroz, L., Hauser, H., Konig, A., andGrollf.r, Ii. 2001. Fast visualization of object contours by non-photorealistic volume rendering. Computer Graphics Forum 20, 3, C452-C460. DOCARMO, M. 1976. Differential Geometry of Curves and Surfaces. Prentice-Hall, Inc., Bnglewood Cliffs, New Jersey. Girshick, A., Interrante, V., Hakf.r, S., and Lf.moine, T. 2000. Line direction matters: An argument for the use of principal direc­tions in 3D line drawing. In First International Symposium on Non- Photorealistic Animation and Rendering, 43-52. GOOCH, B., AND Gooch, A. 2001. Non-Photorealistic Rendering. A K Peters, Ltd., Natick, MA. Gooch, A., Gooch, B., Shirley, P., and Cohen, Ii. 1998. A non- photorealistic lighting model for automatic technical illustraction. In Proceedings SIGGRAPH '98, 447-452. Gooch, B., Sloan, P.-P., Gooch, A., Shirley, P., and Riesenfei.d, R. 1999. Interactive technical illustration. In 1999 ACM Symposium on Interactive 3D Graphics, 31-38. Hadwiger, M., Theussl, T., Hauser, H., andGrollf.r, Ii. 2001. Hardware-accelerated hiqh-quality filtering on pc graphics hardware. In Proceedings Vision. Modeling, and Visualization 2001. Hladuvka, J., Konig, A., andGrollf.r, Ii. 2000. Curvature-based transfer functions for direct volume rendering. In Spring Conference on Computer Graphics 2000, vol. 16, 58-65. Hoffman, K., and Kunzf,, R. 1971. I,inear Algebra. Prentice-Hall, Inc., linglewood Cliffs, MJ, ch. 10.2. Interrante, V., Fuchs, H., and Pizf.r, S. 1995. linhancing trans­parent skin surfaces with ridge and valley line. In Proceedings IEEE Visualization 1995, 52-59. Interrante, V., Fuchs, H., and Pizf.r, S. M. 1997. Conveying the 3D shape of smoothly curving transparent surfaces via texture. IEEE Transactions on Visualization and Computer Graphics J, 2 (April-June), 98-117. Interrante, V. 1997. Illustrating surface shape in volume data via prin­cipal direction-driven 3D line integral convolution. In Proceedings SIG­GRAPH '97, vol. 31, 109-116. Kindlmann, G., and Durkin, J. 1998. Semi-automatic generation of transfer functions for direct volume rendering. In IEEE Symposium On Volume Visualization, 79-86. Kniss, J., Kindlmann, G., and Hansen, C. 2002. Multidimensional transfer functions for interactive volume rendering. IEEE Transactions on Visualization and Computer Graphics 8, 3 (July-Sep), 270-285. Levoy, M. 1988. Display of surfaces from volume data. IEEE Computer Graphics & Applications 8, 5, 29-37. Lu, A., Morris, C., Brert, D., Rhf.ingans, P., and Hansen, C. 2002. Non-photorealistic volume rendering using stippling techniques. In Proceedings IEEE Visualization 2002, 211-218. Lum, Ii., AND Ma, K.-L. 2002. Hardware-accelerated parallel non- photorealistic volume rendering. In International Symposium on Non- Photorealistic Animation and Rendering, 67-74. Lum, Ii. B., Stompf.l, A., and Ma, K.-L. 2003. Using motion to illustrate static 3D shape - kinetic visualization. IEEE Transactions on Visualization and Computer Graphics 9, 2 ( April-June), 115-126. Ma, K.-L., Hf.rtzmann, A., Interrante, V., and Lum, Ii. B. 2002. Course 23: Recent advances in non-photorealistic rendering for art and visualization. In Course Notes for SIGGRAPH 2002. Marschnf.r, S., and Lorr, R. 1994. An evaluation of reconstruction filters for volume rendering. In Proceedings of Visualization '94, R. D. Bergeron and A. Ii. Kaufman, lids., 100-107. MARSDEN, J., AND TROMRA, A. 1996. Vector Calculus. W.H. Freeman and Company, New York, ch. 2.6, 4.2. Mitchell, D., and Hanrahan, P. 1992. Illumination from curved reflectors. In Proceedings SIGGRAPH '92,283-291. Mitchell, D., and Nf.travali, A. 1988. Reconstruction filters in com­puter graphics. In Proceedings SIGGRAPH '88, 221-228. Mollf.r, T., Muller, K., Kurzion, Y., Machiraju, R., and Yagf.l, R. 1998. Design of accurate and smooth filters for function and derivative reconstruction. In IEEE Symposium on Volume Visualiza­tion, 143-151. MOLLER, T. 1999. Spatial Domain Filter Design. PhD thesis, Ohio State University. Monga, O., Benayoun, S., and Faugeras, O. 1992. From partial derivatives of 3D density images to ridge lines. In IEEE Conference on Vision and Pattern Recognition (CVPRI. Monga, O., Lf.ngagne, R., and Df.riche, R. 1994. Crest lines ex­traction in volume 3D medical images: a multi-scale approach. Rapport de Recherche 2338, INRIA, July. Nagy, Schneider, J., and Wf.stermann, R. 2002. Interactive volume illustration. In Proceedings Vision. Modeling, and Visualization 2002. Osher, S., and Sethian, J. 1988. Fronts propagating with curvature dependent speed: Algorithms based on hamilton-jacobi formulation. J. Computational Physics 79, 12-49. Pf.rona, P., and Malik, J. 1990. Scale space and edge detection using anisotropic diffusion. IEEE Trans, on Pattern Analysis and Machine Intelligence 12,7 (July), 629-639. Preusser, T., and Rumpf, M. 2002. A level set method for anisotropic geometric diffusion in 3D image processing. SIAM J. Appl. Math. 62, 5, 1772-1793. Rhf.ingans, P., and Brert, D. 2001. Volume illustration: Nonphotore- alistic rendering of volume models. IEEE Transactions on Visualization and Computer Graphics 7, 3 (July-Sep), 253-264. Tasdizen, T., Whitaker, R., Burchard, P., and Oshf.r, S. 2002. Geometric surface smoothing via anisotropic diffusion of normals. In Proceedings IEEE Visualization 2002, 125-132. ter Haar Romeny, B., Florack, L., Koendf.rink, J., and VlF,RGF,VF,R, M. 1991. Lecture Notes in Computer Science, vol. 511. Springer-Verlag, ch. Scale-space: Its natural operators and differential invariants, in Information Processing in Medical Imaging, 239-255. Trf.avett, S., and Chen, M. 2000. Pen-and-ink rendering in volume visualization. In Proceedings IEEE Visualization 2000, 203-210. 520 Proceedings of the 14th IEEE Visualization Conference (VIS'03) 0-7695-2030-8/03 $ 17.00 © 2003 IEEE