EWA Splatting

In this paper, we present a framework for high quality splatting based on elliptical Gaussian kernels. To avoid aliasing artifacts, we introduce the concept of a resampling ﬁlter combining a reconstruction kernel with a low-pass ﬁlter. Because of the similarity to Heckberts EWA (elliptical weighted average) ﬁlter for texture mapping, we call our technique EWA splatting. Our framework allows us to derive EWA splat primitives for volume data and for point-sampled surface data. It provides high image quality without aliasing artifacts or excessive blurring for volume data, and additionally features anisotropic texture ﬁltering for point-sampled surfaces. It also handles non-spherical volume kernels efﬁciently, hence it is suitable for regular, rectilinear, and irregular volume datasets. Moreover, our framework introduces a novel approach to compute the footprint function, facilitating efﬁcient perspective projection of arbitrary elliptical kernels at very little additional cost. Finally, we show that EWA volume reconstruction kernels can be reduced to surface reconstruction kernels. This makes our splat primitive universal in rendering surface and volume data.


Introduction
Volume rendering is an important technique in visualizing acquired and simulated data sets in scientific and engineering applications.The ideal volume rendering algorithm reconstructs a continuous function in 3D, transforms this 3D function into screen space, and then evaluates opacity integrals along line-of-sights.In 1989, Westover [18,19] introduced splatting for interactive volume rendering, which approximates this procedure.Splatting algorithms interpret volume data as a set of particles that are absorbing and emitting light.Line integrals are precomputed across each particle separately, resulting in footprint functions.Each footprint spreads its contribution in the image plane.These contributions are composited back to front into the final image.
We introduce a new footprint function for volume splatting algorithms integrating an elliptical Gaussian reconstruction kernel and a low-pass filter.Our derivation proceeds along similar lines as Heckbert's elliptical weighted average (EWA) texture filter [4], therefore we call our algorithm EWA volume splatting.
EWA volume rendering is attractive because it prevents aliasing artifacts in the output image while avoiding excessive blurring.Moreover, it works with arbitrary elliptical Gaussian reconstruction kernels and efficiently supports perspective projection.Our method is based on a novel framework to compute the footprint function, which relies on the transformation of the volume data to so-called ray space.This transformation is equivalent to perspective projection.By using its local affine approximation at each voxel, we derive an analytic expression for the EWA footprint in screen space.The rasterization of the footprint is performed using forward differ-encing requiring only one 1D footprint table for all reconstruction kernels and any viewing direction.
Our splat primitive can be integrated easily into conventional splatting algorithms.Because of its flexibility, it can be utilized to render rectilinear, curvilinear, or unstructured volume data sets.By flattening the 3D Gaussian kernel along the volume gradient we will show that EWA volume splats reduce to surface splats that are suitable for high quality iso-surface rendering.
The paper is organized as follows: We discuss previous work in Section 2. Next, we review a typical volume rendering pipeline and the volume rendering equation in Section 3. Specifically, we elaborate how the volume rendering equation is computed by typical splatting algorithms.In Section 4, we present our EWA volume rendering framework.We start by analyzing the aliasing problem due to improper sampling of the output function resulting from volume rendering.In a next step, we introduce the EWA resampling filter, which integrates an arbitrary elliptical Gaussian reconstruction kernel and a Gaussian low-pass filter.Our derivation is based on the local affine transformation of the volume data such that the reconstruction kernels can be integrated analytically.Furthermore, we show how the EWA reconstruction kernels can be continuously adapted from volumes to surfaces in Section 5. Finally, Sections 6 and 7 discuss our implementation and results.

Previous Work
The original work on splatting was presented by Westover [18].Basic splatting algorithms suffer from inaccurate visibility determination when compositing the splats from back to front.This leads to visible artifacts such as color bleeding.Later, Westover [19] solved the problem using an axis-aligned sheet buffer.However, this technique is plagued by disturbing popping artifacts in animations.Recently, Mueller and Crawfis [14] proposed to align the sheet buffers parallel to the image plane instead of parallel to an axis of the volume data.Additionally, they splat several slices of each reconstruction kernel separately.This technique is similar to slice-based volume rendering [2,1] and does not suffer from popping artifacts.Mueller and Yagel [15] combine splatting with ray casting techniques to accelerate rendering with perspective projection.Laur and Hanrahan [7] describe a hierarchical splatting algorithm enabling progressive refinement during rendering.Furthermore, Lippert [9] introduced a splatting algorithm that directly uses a wavelet representation of the volume data.
Westover's original framework does not deal with sampling rate changes due to perspective projections.Aliasing artifacts may occur in areas of the volume where the sampling rate of diverging rays falls below the volume grid sampling rate.Swan et al. [17] use a distance-dependent stretch of the footprints to make them act as low-pass filters.This antialiasing method is closely related to EWA volume splatting, and we will discuss it further in Section 7.
Additional care has to be taken if the 3D kernels are not radially symmetric, as is the case for rectilinear, curvilinear, or irregular grids.In addition, for an arbitrary position in 3D, the contributions from all kernels must sum up to one.Otherwise, artifacts such as splotches occur in the image.For rectilinear grids, Westover [19] proposes using elliptical footprints that are warped back to a circular footprint.To render curvilinear grids, Mao et al. [10] use stochastic Poisson resampling to generate a set of new points whose kernels are spheres or ellipsoids.They compute the elliptical footprints very similar to Westover [19].As pointed out in Section 4, our technique can be used with irregular grids to efficiently and accurately project and rasterize the elliptical splat kernels.
We develop EWA volume splatting along similar lines to the seminal work of Heckbert [4], who introduced EWA filtering to avoid aliasing of surface textures.We recently extended his framework to represent and render texture functions on irregularly pointsampled surfaces [21].Section 5 will show the connection between EWA volume and surface splatting.

The Volume Rendering Pipeline
We distinguish two fundamental approaches to volume rendering: backward mapping algorithms that shoot rays through pixels on the image plane into the volume data, and forward mapping algorithms that map the data onto the image plane.In the following discussion, we will describe a forward mapping technique.Mapping the data onto the image plane involves a sequence of intermediate steps where the data is transformed to different coordinate systems, as in conventional rendering pipelines.We introduce our terminology in Figure 1.Note that the terms space and coordinate system are synonymous.The figure summarizes a forward mapping volume rendering pipeline, where the data flows from the left to the right.As an overview, we briefly describe the coordinate systems and transformations that are relevant for our technique.The volume data is initially given in object coordinates.To render the data from an arbitrary viewpoint, it is first mapped to camera space using the viewing transformation.We deal with the effect of this transformation in Section 4.3.The camera coordinate system is defined such that its origin is at the center of projection.
We further transform the data to ray space, which is introduced in Section 3.2.Ray space is a non-cartesian coordinate system that enables an easy formulation of the volume rendering equation.In ray space, the viewing rays are parallel to a coordinate axis, facilitating analytical integration of the volume function.We present the transformation from camera to ray space in Section 4.4; it is a key element of our technique.Since its purpose is similar to the projective transform used in rendering pipelines such as OpenGL, it is also called the projective mapping.
Evaluating the volume rendering equation results in a 2D image in screen space.In a final step, this image is transformed to viewport coordinates.Focusing on the essential aspects of our technique, we are not covering the viewport transformation in the following explanations.However, it can be easily incorporated in an implementation.Moreover, we do not discuss volume classification and shading in a forward mapping pipeline, but refer to [13] or [20] for a thorough discussion.

Splatting Algorithms
We review the low albedo approximation of the volume rendering equation [5,12] as used for fast, direct volume rendering [19,6,13,8].The left part of Figure 2 illustrates the corresponding situation in 2D.Starting from this form of the rendering equation, we discuss several simplifying assumptions leading to the well known splatting formulation.Because of their efficiency, splatting algorithms [19,13] belong to the most popular forward mapping volume rendering techniques.
We slightly modify the conventional notation, introducing our concept of ray space.We denote a point in ray space by a column vector of three coordinates x = (x0, x1, x2) T .Given a center of projection and a projection plane, these three coordinates are interpreted geometrically as follows: The coordinates x0 and x1 specify a point on the projection plane.The ray intersecting the center of projection and the point (x0, x1) on the projection plane is called a viewing ray.Using the abbreviation x = (x0, x1) T , we refer to the viewing ray passing through (x0, x1) as x.The third coordinate x2 specifies the Euclidean distance from the center of projection to a point on the viewing ray.To simplify the notation, we will use any of the synonyms x, (x, x2) T , or (x0, x1, x2) T to denote a point in ray space.The volume rendering equation describes the light intensity I λ (x) at wavelength λ that reaches the center of projection along the ray x with length L: where g(x) is the extinction function that defines the rate of light occlusion, and c λ (x) is an emission coefficient.The exponential term can be interpreted as an attenuation factor.Finally, the product c λ (x)g(x) is also called the source term [12], describing the light intensity scattered in the direction of the ray x at the point x2.Now we assume that the extinction function is given as a weighted sum of coefficients g k and reconstruction kernels r k (x).This corresponds to a physical model where the volume consists of individual particles that absorb and emit light.Hence the extinction function is: In this mathematical model, the reconstruction kernels r k (x) reflect position and shape of individual particles.The particles can be irregularly spaced and may differ in shape, hence the representation in ( 2) is not restricted to regular data sets.We substitute (2) into (1), yielding: To compute this function numerically, splatting algorithms make several simplifying assumptions, illustrated in the right part of Figure 2. Usually the reconstruction kernels r k (x) have local support.The splatting approach assumes that these local support areas do not overlap along a ray x, and the reconstruction kernels are ordered front to back.We also assume that the emission coefficient is constant in the support of each reconstruction kernel along a ray, hence we use the notation c λk (x) = c λ (x, x2), where (x, x2) is in the support of r k .Moreover, we approximate the exponential function with the first two terms of its Taylor expansion, thus e x ≈ 1−x.Finally, we ignore self-occlusion.Exploiting these assumptions, we rewrite (3), yielding: where q k (x) denotes an integrated reconstruction kernel, hence: Equation ( 4) is the basis for all splatting algorithms.Westover [19] introduced the term footprint function for the integrated reconstruction kernels q k .The footprint function is a 2D function that specifies the contribution of a 3D kernel to each point on the image plane.
Integrating a volume along a viewing ray is analogous to projecting a point on a surface onto the image plane, hence the coordinates x = (x0, x1) T are also called screen coordinates, and we say that I λ (x) and q k (x) are defined in screen space.
Splatting is attractive because of its efficiency, which it derives from the use of pre-integrated reconstruction kernels.Therefore, during volume integration each sample point along a viewing ray is computed using a 2D convolution.In contrast, ray-casting methods require a 3D convolution for each sample point.This provides splatting algorithms with an inherent advantage in rendering efficiency.Moreover, splatting facilitates the use of higher quality kernels with a larger extent than the trilinear kernels typically employed by ray-casting.On the other hand, basic splatting methods are plagued by artifacts because of incorrect visibility determination.This problem is unavoidably introduced by the assumption that the reconstruction kernels do not overlap and are ordered back to front.It has been successfully addressed by several authors as mentioned in Section 2. In contrast, our main contribution is a novel splat primitive that provides high quality antialiasing and efficiently supports elliptical kernels.We believe that our novel primitive is compatible with all state-of-the-art algorithms.

Aliasing in Volume Splatting
Aliasing is a fundamental problem of any rendering algorithm, arising whenever a rendered image or a part of it is sampled to a discrete raster grid, i.e., the pixel grid.Aliasing leads to visual artifacts such as jagged silhouette edges and Moiré patterns in textures.Typically, these problems become most disturbing during animations.From a signal processing point of view, aliasing is well understood: before a continuous function is sampled to a regular sampling grid, it has to be band-limited to respect the Nyquist frequency of the grid.This guarantees that there are no aliasing artifacts in the sampled image.In this section we provide a systematic analysis on how to band-limit the splatting equation.
The splatting equation ( 4) represents the output image as a continuous function I λ (x) in screen space.In order to properly sample this function to a discrete output image without aliasing artifacts, it has to be band-limited to match the Nyquist frequency of the discrete image.In theory, we achieve this band-limitation by convolving I λ (x) with an appropriate low-pass filter h(x), yielding the antialiased splatting equation Although I λ (x) is formulated as a continuous function in (4), in practice this function is evaluated only at discrete positions, i.e., the pixel centers.Therefore we cannot evaluate (6), which requires that I λ (x) is available as a continuous function.However, we make two simplifying assumptions to rearrange the integral in (6).This leads to an approximation that can be evaluated efficiently.First, we assume that the emission coefficient is approximately constant in the support of each footprint function q k , hence c λk (x) ≈ c λk for all x in the support area.Together with the assumption that the emission coefficient is constant in the support of each reconstruction kernel along a viewing ray, this means that the emission coefficient is constant in the complete 3D support of each reconstruction kernel.In other words, we ignore the effect of shading for antialiasing.Note that this is the common approach for antialiasing surface textures as well.
Additionally, we assume that the attenuation factor has an approximately constant value o k in the support of each footprint function.Hence: for all x in the support area.A variation of the attenuation factor indicates that the footprint function is partially covered by a more opaque region in the volume data.Therefore this variation can be interpreted as a "soft" edge.Ignoring such situations means that we cannot prevent edge aliasing.Again, this is similar to rendering surfaces, where edge and texture aliasing are handled by different algorithms as well.
Exploiting these simplifications, we can rewrite (6) to: Following Heckbert's terminology [4], we call: an ideal resampling filter, combining a footprint function q k and a low-pass kernel h.Hence, we can approximate the antialiased splatting equation ( 6) by replacing the footprint function q k in the original splatting equation ( 4) with the resampling filter ρ k .This means that instead of band-limiting the output function I λ (x) directly, we band-limit each footprint function separately.Under the assumptions described above, we get a splatting algorithm that produces a band-limited output function respecting the Nyquist frequency of the raster image, therefore avoiding aliasing artifacts.Remember that the reconstruction kernels are integrated in ray space, resulting in footprint functions that vary significantly in size and shape across the volume.Hence the resampling filter in ( 8) is strongly space variant.Swan et al. presented an antialiasing technique for splatting [17] that is based on a uniform scaling of the reconstruction kernels to band-limit the extinction function.Their technique produces similar results as our method for radially symmetric kernels.However, for more general kernels, e.g., elliptical kernels, uniform scaling is a poor approximation of ideal low-pass filtering.Aliasing artifacts cannot be avoided without introducing additional blurriness.On the other hand, our method provides non-uniform scaling in these cases, leading to superior image quality as illustrated in Section 7.Moreover, our analysis above shows that band-limiting the extinction function does not guarantee aliasing free images.Because of shading and edges, frequencies above the Nyquist limit persist.However, these effects are not discussed in [17].

Elliptical Gaussian Kernels
We choose elliptical Gaussians as reconstruction kernels and lowpass filters, since they provide certain features that are crucial for our technique: Gaussians are closed under affine mappings and convolution, and integrating a 3D Gaussian along one coordinate axis results in a 2D Gaussian.These properties enable us to analytically compute the resampling filter in (8) as a single 2D Gaussian, as will be shown below.In this section, we summarize the mathematical features of the Gaussians that are exploited in our derivation in the following sections.More details on Gaussians can be found in Heckbert's master's thesis [4].
We define an elliptical Gaussian GV(x − p) centered at a point p with a variance matrix V as: where |V| is the determinant of V.In this form, the Gaussian is normalized to a unit integral.In the case of volume reconstruction kernels, GV is a 3D function, hence V is a symmetric 3 × 3 matrix and x and p are column vectors (x0, x1, x2) T and (p0, p1, p2) T , respectively.We can easily apply an arbitrary affine mapping u = Φ(x) to this Gaussian.Let us define the affine mapping as Φ(x) = Mx + c, where M is a 3 × 3 matrix and c is a vector (c0, c1, c2) T .We substitute x = Φ −1 (u) in ( 9), yielding: Moreover, convolving two Gaussians with variance matrices V and Y results in another Gaussian with variance matrix V + Y: Finally, integrating a 3D Gaussian GV along one coordinate axis yields a 2D Gaussian G V, hence: where x = (x0, x1) T and p = (p0, p1) T .The 2 × 2 variance matrix V is easily obtained from the 3 × 3 matrix V by skipping the third row and column: In the following sections, we describe how to map arbitrary elliptical Gaussian reconstruction kernels from object to ray space.Our derivation results in an analytic expression for the kernels in ray space r k (x) as in Equation ( 2).We will then be able to analytically integrate the kernels according to Equation ( 5) and to convolve the footprint function q k with a Gaussian low-pass filter h as in Equation (8), yielding an elliptical Gaussian resampling filter ρ k .

The Viewing Transformation
The reconstruction kernels are initially given in object space, which has coordinates t = (t0, t1, t2) T .Let us denote the Gaussian reconstruction kernels in object space by r k (t) = G V k (t − t k ), where t k are the voxel positions in object space.For regular volume data sets, the variance matrices V k are usually identity matrices.For rectilinear data sets, they are diagonal matrices where the matrix elements contain the squared distances between voxels along each coordinate axis.Curvilinear and irregular grids have to be resampled to a more regular structure in general.For example, Mao et al. [11] describe a stochastic sampling approach with a method to compute the variance matrices for curvilinear volumes.
We denote camera coordinates by a vector u = (u0, u1, u2) T .Object coordinates are transformed to camera coordinates using an affine mapping u = ϕ(t), called viewing transformation.It is defined by a matrix W and a translation vector d as ϕ(t) = Wt + d.We transform the reconstruction kernels G V k (t − t k ) to camera space by substituting t = ϕ −1 (u) and using Equation (10): where u k = ϕ(t k ) is the center of the Gaussian in camera coordinates and r k (u) denotes the reconstruction kernel in camera space.
According to (10), the variance matrix in camera coordinates V k is given by

The Projective Transformation
The projective transformation converts camera coordinates to ray coordinates as illustrated in Figure 3. Camera space is defined such that the origin of the camera coordinate system is at the center of projection and the projection plane is the plane u2 = 1.Camera space and ray space are related by the mapping x = m(u).Using the definition of ray space from Section 3, m(u) and its inverse m −1 (x) are therefore given by: where l = (x0, x1, 1) T .Unfortunately, these mappings are not affine, so we cannot apply Equation (10) directly to transform the reconstruction kernels from camera to ray space.To solve this problem, we introduce the local affine approximation mu k of the projective transformation.It is defined by the first two terms of the Taylor expansion of m at the point u k : where x k = m(u k ) is the center of a Gaussian in ray space.The Jacobian Ju k is given by the partial derivatives of m at the point u k : In the following discussion, we are omitting the subscript u k , hence m(u) denotes the local affine approximation (17).We substitute u = m −1 (x) in ( 14) and apply Equation (10) to map the reconstruction kernels to ray space, yielding the desired expression for r k (x): where V k is the variance matrix in ray coordinates.According to (10), V k is given by: Note that for uniform or rectilinear data sets, V k has to be computed only once per frame, since V k is the same for all voxels and W changes only from frame to frame.However, since the Jacobian is different for each voxel position, V k has to be recalculated for each voxel.In the case of curvilinear or irregular volumes, each reconstruction kernel has an individual variance matrix V k .Our method efficiently handles this situation, requiring only one additional 3 × 3 matrix multiplication.In contrast, previous techniques [19,11] cope with elliptical kernels by computing their projected extents in screen space and then establishing a mapping to a circular footprint table.However, this procedure is computationally expensive.It leads to a bad approximation of the integral of the reconstruction kernel as pointed out in [15,17].
As illustrated in Figure 4, the local affine mapping is exact only for the ray passing through u k or x k , respectively.The figure is exaggerated to show the non-linear effects in the exact mapping.The affine mapping essentially approximates the perspective projection with an oblique orthographic projection.Therefore, parallel lines are preserved, and approximation errors grow with increasing ray divergence.However, the errors do not lead to visual artifacts in general [15], since the fan of rays intersecting a reconstruction kernel has a small opening angle due to the local support of the reconstruction kernels.
A common approach of performing splatting with perspective projection is to map the footprint function onto a footprint polygon in camera space in a first step.In the next step, the footprint polygon is projected to screen space and rasterized, resulting in the so-called footprint image.As mentioned in [15], however, this requires significant computational effort.In contrast, our framework efficiently performs perspective projection by mapping the volume to ray space, which requires only the computation of the Jacobian and two 3 × 3 matrix multiplications.For spherical reconstruction kernels, these matrix operations can be further optimized as shown in Section 6.

Integration and Band-Limiting
We integrate the Gaussian reconstruction kernel of ( 19) according to (5), resulting in a Gaussian footprint function q k : where the 2 × 2 variance matrix Vk of the footprint function is obtained from V k by skipping the last row and column, as shown in (13).
Finally, we choose a Gaussian low-pass filter h = G V h (x), where the variance matrix V h ∈ Ê 2×2 is typically the identity matrix.With (11), we compute the convolution in (8), yielding the EWA volume resampling filter:

Reduction from Volume to Surface Reconstruction Kernels
Since our EWA volume resampling filter can handle arbitrary Gaussian reconstruction kernels, we can represent the structure of a volume data set more accurately by choosing the shape of the reconstruction kernels appropriately.For example, we can improve the precision of isosurface rendering by flattening the reconstruction kernels in the direction of the surface normal.We will show below that an infinitesimally flat Gaussian volume kernel is equivalent to a Gaussian surface texture reconstruction kernel [21].In other words, we can extract and render a surface representation from a volume data set directly by flattening volume reconstruction kernels into surface reconstruction kernels.Our derivation is illustrated in Figure 5.
3D viewing transformation integration 2D to 3D parameterization 3D to 2D projection We construct a flattened Gaussian reconstruction kernel in object space by scaling a spherical Gaussian in one direction by a factor 1/s, hence its variance matrix is: A scaling factor s = 1 corresponds to a spherical 3D kernel.In the limit, if s = ∞, we get a circular 2D kernel.
To render this reconstruction kernel, we first apply a 3D transformation matrix W, which may contain arbitrary modeling transformations concatenated with the viewing transformation.Then we use the local affine approximation of Equation (17) to map the kernel to ray space.The variance matrix V of the reconstruction kernel in ray space is computed as in (20).We introduce the matrix T 3D to denote the concatenated 3D mapping matrix T 3D = JW and write V as: Hence, the elements vij of V are given by: where we denote an element of T 3D by tij.We compute the 2D Gaussian footprint function by integrating the reconstruction kernel.According to (13), its 2D variance matrix is obtained by skipping the third row and column in V.As s approaches infinity, we therefore get the following 2D variance matrix V: Conveniently, the 2D variance matrix can be factored into a 2D mapping T 2D , which is obtained from the 3D mapping matrix by skipping the third row and column: Let us now analyze the 2D mapping matrix T 2D .First, we need an explicit expression for the Jacobian J of the projective mapping.Using (15) and (18), it is given by: where l = (u0, u1, u2) T .With T 3D = JW, we use the first two rows of J and the first two columns of W to factor T 2D into: where wij denotes an element of W. This can be interpreted as a concatenation of a 2D to 3D with a 3D to 2D mapping, resulting in a compound 2D to 2D mapping similar as in conventional texture mapping [3].We illustrate this process schematically in Figure 5 and more intuitively in Figure 6.The first stage is a parameterization of a 3D plane.It maps a circular 2D texture kernel onto a plane defined by the two vectors (w00, w10, w20) T and (w01, w11, w21) T in 3D camera space, resulting in an ellipse.The second stage is an oblique parallel projection with an additional scaling factor 1/u2, which is the local affine approximation of the perspective projection.Finally, combining the projected ellipse with a low-pass filter as in Equation ( 8) yields a texture filter that is equivalent to Heckbert's EWA filter [4].This is the same result as we derive in [21].
We compare splatting with volumetric kernels and splatting with surface kernels in Section 7.

Implementation
We implemented a volume rendering algorithm based on the EWA splatting equation.Our implementation is embedded in the VTK (visualization toolkit) framework [16].We did not optimize our code for rendering speed.We use a sheet buffer to first accumulate splats from planes in the volume that are most parallel to the projection plane [19].In a second step, the final image is computed by compositing the sheets back to front.Shading is performed using the gradient estimation functionality provided by VTK and the Phong illumination model.We summarize the main steps which are required to compute the EWA splat for each voxel: First, we compute the camera coordinates u k of the current voxel k by applying the viewing transformation to the voxel center.Using u k , we then evaluate the Jacobian J as given in Equation (25).In line 4, we transform the Gaussian reconstruction kernel from object to ray space.This transformation is implemented by Equation (20), and it results in the 3 × 3 variance matrix V k of the Gaussian in ray space.Remember that W is the rotational part of the viewing transformation, hence it is typically orthonormal.Moreover, for spherical kernels, V k is the identity matrix.In this case, evaluation of Equation ( 20) can be simplified significantly.Next, we project the voxel center from camera space to the screen by performing a perspective division on u k .This yields the 2D screen coordinates xk .Now we are ready to setup the resampling filter ρ k according to Equation ( 22).Its variance matrix is derived from V k by omitting the third row and column, and adding a 2 × 2 identity matrix for the low-pass filter.Moreover, we compute the determinants 1/|J −1 | and 1/|W| −1 that are used as normalization factors.
Finally, we rasterize the resampling filter in line 7.As can be seen from the definition of the elliptical Gaussian (9), we also need the inverse of the variance matrix, which is called the conic matrix.Let us denote the 2 × 2 conic matrix of the resampling filter by Q.Furthermore, we define the radial index function Note that the contours of the radial index, i.e., r = const.are concentric ellipses.For circular kernels, r is the squared distance to the circle center.The exponential function in (9) can now be written as e − 1 2 r .We store this function in a 1D lookup table .To evaluate the radial index efficiently, we use finite differencing.Since r is biquadratic in x, we need only two additions to update r for each pixel.We rasterize r in a rectangular, axis aligned bounding box centered around xk as illustrated in Figure 7.Typically, we use a threshold c = 4 and evaluate the Gaussian only if r(x) < c.Heckbert provides pseudo-code of the rasterization algorithm in [4].

Results
The EWA resampling filter has a number of useful properties, as illustrated in Figure 8.When the mapping from camera to ray space minifies the volume, size and shape of the resampling filter are dominated by the low-pass filter, as in the left column of Figure 8.In the middle column, the volume is magnified and the resampling filter is dominated by the reconstruction kernel.Since the resampling filter unifies a reconstruction kernel and a low-pass filter, it provides a smooth transition between magnification and minification.Moreover, the reconstruction kernel is scaled anisotropically in situations where the volume is stretched in one direction and shrinked in the other, as shown in the right column.In the bottom row, we show the filter shapes resulting from uniformly scaling the reconstruction kernel to avoid aliasing, as proposed by Swan et al. [17].Essentially, the reconstruction kernel is enlarged such that its minor radius is at least as long as the minor radius of the lowpass filter.For spherical reconstruction kernels, or circular footprint functions, this is basically equivalent to the EWA resampling filter.However, for elliptical footprint functions, uniform scaling leads to overly blurred images in the direction of the major axis of the ellipse.We compare our method to Swan's method in Figure 8 (see colorplate).The images on the left were rendered with EWA volume splats, those on the right with Swan's uniformly scaled kernels.We used a square zebra texture with x and y dimensions of 1024 × 512 in the first row, and 1024 × 256 in the second row.This leads to elliptical reconstruction kernels with a ratio between the length of the major and minor radii of 2 to 1 and 4 to 1, respectively.Clearly, the EWA filter produces a crisper image and at the same time does not exhibit aliasing artifacts.As the ratio between the major and minor radii of the reconstruction kernels increases, the difference to Swan's method becomes more pronounced.For strongly anisotropic kernels, Swan's uniform scaling produces excessively

To appear in the Proceedings of IEEE Visualization 2001
blurred images, as shown on the right in Figure 8 .Each frame took approximately 6 seconds to render on an 866 MHz PIII processor.
In Figure 9 (see colorplate), we compare EWA splatting using volume kernels on the left to surface reconstruction kernels on the right.The texture size is 512 × 512 in x and y direction.Typically, the perspective projection of a spherical kernel is almost a circle.Therefore, rendering with volume kernels does not exhibit anisotropic texture filtering and produces textures that are slightly too blurry, similar to isotropic texture filters such as trilinear mipmapping.On the other hand, splatting surface kernels is equivalent to EWA texture filtering.Circular surface kernels are mapped to ellipses, which results in high image quality because of anisotropic filtering.
In Figure 10 (see colorplate), we show a series of volume renderings of the UNC CT scan of a human head (256 × 256 × 225), the UNC engine (256 × 256 × 110), and the foot of the visible woman dataset (152 × 261 × 220).The texture in the last example is rendered using EWA surface splatting, too.The images illustrate that our algorithm correctly renders semitransparent objects as well.The skull of the UNC head, the bone of the foot, and the iso-surface of the engine were rendered with flattened surface splats oriented perpendicular to the volume gradient.All other voxels were rendered with EWA volume splats.Each frame took approximately 11 seconds to render on an 866 MHz PIII processor.

Conclusion and Future Work
We present a new splat primitive for volume rendering, called the EWA volume resampling filter.Our primitive provides high quality antialiasing for splatting algorithms, combining an elliptical Gaussian reconstruction kernel with a Gaussian low-pass filter.We use a novel approach of computing the footprint function.Exploiting the mathematical features of 2D and 3D Gaussians, our framework efficiently handles arbitrary elliptical reconstruction kernels and perspective projection.Therefore, our primitive is suitable to render regular, rectilinear, curvilinear, and irregular volume data sets.Finally, we derive a formulation of the EWA surface reconstruction kernel, which is equivalent to Heckbert's EWA texture filter.Hence we call our primitive universal, facilitating the reconstruction of surface and volume data.
We have not yet investigated whether other kernels besides elliptical Gaussians may be used with this framework.In principle, a resampling filter could be derived from any function that allows the analytic evaluation of the operations described in Section 4.2 and that is a good approximation of an ideal low-pass filter.
To achieve interactive frame rates, we are currently investigating the use of graphics hardware to rasterize EWA splats as texture mapped polygons.We also plan to use sheet-buffers that are parallel to the image plane to eliminate popping artifacts.To render non-rectilinear datasets we are investigating fast back-to-front sorting algorithms.Furthermore, we want to experiment with our splat primitive in a post-shaded volume rendering pipeline.The derivative of the EWA resampling filter could be used as a bandlimited gradient kernel, hence avoiding aliasing caused by shading for noisy volume data.Finally, we want to exploit the ability of our framework to render surface splats.In conjunction with voxel culling algorithms we believe it is useful for real-time iso-surface rendering.
Swan et al.

Figure 1 :
Figure 1: The forward mapping volume rendering pipeline.

Figure 2 :
Figure 2: Volume rendering.Left: Illustrating the volume rendering equation in 2D.Right: Approximations in typical splatting algorithms.

Figure 3 :
Figure 3: Transforming the volume from camera to ray space.Top: camera space.Bottom: ray space.

Figure 4 :
Figure 4: Mapping a reconstruction kernel from camera to ray space.Top: camera space.Bottom: ray space.Left: local affine mapping.Right: exact mapping.

1 2DFigure 5 :
Figure 5: Reducing a volume reconstruction kernel to a surface reconstruction kernel by flattening the kernel in one dimension.Top: rendering a volume kernel.Bottom: rendering a surface kernel.

Figure 8 :
Figure 8: Properties of the EWA resampling filter

Figure 10 :
Figure 10: Semitransparent objects rendered using EWA volume splatting.The skull of the UNC head, the iso-surface of the engine, and the bone of the foot are rendered with flattened surface splats.