Chapter 3
Data Representation
IN Chapter 2, we used a simple application, the visualization of a two-variable function using a height plot, to introduce several main concepts of data visualization, such as datasets, sampling, and rendering. Our discussion revolved around the dataset as a key concept for approximating a continuous signal by means of a discrete representation.
In this chapter, we expand this discussion on discrete data representation and approximation. First, we detail the notions of continuous data sampling and reconstruction. We introduce basis functions, discrete meshes, and cells, as a means of constructing piecewise continuous approximations from sampled data. These notions are illustrated using the simple example of visualizing a two-variable function with a height plot, introduced in Chapter 2. Next, we present the various types of datasets commonly used in the visualization practice and detail their relative advantages, limitations, and constraints. We also discuss several aspects related to efficiently implementing the dataset concept. After completing this chapter, the reader should have a good understanding of the various trade-offs involved in the choice of a dataset model, and be able to decide which is the most appropriate dataset for a given visualization application.
3.1 Continuous Data
3.1.1 What Is Continuous Data?
The main goal of visualization is to produce pictures that enable end users get insight into data that describe various phenomena or processes, including both natural and human-controlled phenomena. Examples are the circulation of air in the atmosphere due to weather conditions; the flow of water in oceans and seas; convection of air in closed spaces such as buildings due to temperature differences caused by heating and cooling; the deformation, electrostatic charging, vibration, and heating of mechanical machine parts subjected to operational stress; and the magnetic and electrostatic fields generated by electrical engines.
Such phenomena are modeled in terms of various physical quantities. These quantities can be either directly measured or computed by software simulations. From the point of view of data representation, these quantities can be classified in two fundamentally different categories: intrinsically continuous and intrinsically discrete ones. Well-known examples of continuous quantities are pressure, temperature, position, speed, density, force, color, light intensity, and electromagnetic radiation. Examples of intrinsically discrete data are the hypermedia, e.g., text and image, contents of web pages; software source code; plain text, as found in documents of various types; and various kinds of database records.
For the sake of brevity, we shall from now on refer to intrinsically continuous data as continuous data and to intrinsically discrete data as discrete data, respectively. Continuous data are usually manipulated by computers in some finite approximative form. To avoid confusion, we shall refer to this as sampled data. Sampled data are also discrete, in the sense that they consist of a finite set of data elements. However, in contrast to what we call intrinsically discrete data, sampled data always originates from, and is intended to approximate, a continuous quantity. As we shall see in Section 3.2, we can always go from sampled data back to a continuous approximation of the original continuous data. In contrast, what we call intrinsically discrete data has no counterpart in the continuous world, as is the case of a page of text, for example. This is a fundamental difference between continuous (or sampled) data and discrete data. As we shall see throughout this book, this difference causes many, often subtle, constraints and choices in designing effective and efficient visualization methods.
Intrinsically continuous data, whether coming in its original continuous representation or represented by some sampled approximation, is the subject of a separate branch of data visualization. This branch is traditionally called scientific visualization, or scivis. The name reflects the fact that continuous data usually measure physical quantities that are studied by various scientific and engineering disciplines, such as physics, chemistry, mechanics, or engineering. Intrinsically discrete data, which has no direct continuous counterpart, is the subject of a relatively younger visualization branch, called information visualization, or infovis. This book mainly focuses on scivis methods and techniques, with the exception of Chapter 11, which performs a short overview of a few main infovis topics. Consequently, the remainder of this chapter will discuss solely the representation of continuous and sampled (scivis) data.
3.1.2 Mathematical Continuity
Mathematically, continuous data can be modeled as a function
where D ⊂ ℝd is the function domain and C ⊂ ℝc is the function codomain, respectively. In a related terminology, f is called a d-dimensional, or d-variate, c-value function. In other words, if we write f (x) = y, where x ∈ D and y ∈ C, this actually means f (x1,...,xd) = (y1,...,yc). In visualization applications, f or its sampled counterpart, introduced in Section 3.2, is sometimes called a field. Mathematically, f is continuous if for every point p ∈ C the following holds:
This continuity criterion is known as the Cauchy criterion, after the French eighteenth-century mathematician who proposed it, or the ϵ − δ criterion, a name that follows from its formulation.
What does continuity mean in the intuitive sense? In plain terms, a function is continuous if the graph of the function is a connected surface without “holes” or “jumps.” Furthermore, we say that a function f is continuous of order k if the function itself and all its derivatives up to and including order k are also continuous in this sense. This is denoted as f ∈ 𝒞k. Figure 3.1 illustrates the continuity concept for the case of a one-dimensional function f : ℝ → ℝ, whose graph is displayed in green. The first image (Figure 3.1(a)) shows a function having a discontinuity at the points x0. Clearly, if we evaluate the jump of f around the point x0, i.e., evaluate f (x0 − δ) and f (x0 + δ), the jump |f (x0 + δ) − f (x0 − δ)| will always have the fixed value b − a, regardless of how small the value δ is, i.e., how close we sample the function to the point x0. The derivative f′(x) is shown in red in the figure. Note that f′(x) is not defined for the discontinuity point x0.1
Figure 3.1. Function continuity. (a) Discontinuous function. (b) First-order 𝒞0 continuous function. (c) High-order 𝒞k continuous function.
Figure 3.1(b) shows a continuous function f ∈ 𝒞0. Here, f consists of three linear components. The graph of f shows no holes or jumps, so clearly Equation (3.1) holds. In contrast, its first derivative f′ is of the same type as the function f shown before in Figure 3.1, i.e., it is discontinuous. However, if we take the three intervals delimited by the two points x1 and x2 and ignore the points x1 and x2 themselves, the derivative f′ is continuous for every interval. Functions f whose derivatives are continuous on compact intervals as in this example are also called piecewise continuously differentiable, piecewise smooth, or piecewise 𝒞∞. Such functions will play an important role in approximating sampled data later in Section 3.2. A final example is shown in Figure 3.1(c). Here, both the function and its derivative are continuous, so f ∈ 𝒞k, where k ≥ 1. Some functions, such as the Gaussian used in Chapter 2, are themselves and all their derivatives continuous. We denote this by f ∈ 𝒞∞.
3.1.3 Dimensions: Geometry, Topology, and Attributes
The triplet 𝒟 = (D,C,f) defines a continuous dataset. In the following, we shall use the notation 𝒟 to refer to a dataset and the notation D to refer to a domain. The dimension d of the space ℝd into which the function domain D is embedded, or contained, is called the geometrical dimension. The dimension s ≤ d of the function domain D itself is called the topological dimension of the dataset. Understanding the difference between geometrical and topological dimensions is easiest by means of an example. If D is a plane or curved surface embedded in the usual Euclidean space ℝ3, then we have s = 2 and d = 3. If D is a line or curve embedded in the Euclidean space ℝ3, then we have s = 1 and d = 3. You can think of the topological dimension as the number of independent variables that we need to represent our domain D. A curved surface in ℝ3 can actually be represented by two independent variables, like latitude and longitude for the surface of the Earth. Indeed, we can describe such a surface by an implicit function f (x, y, z) = 0, which actually says that only two of the three variables x, y, and z vary independently. Similarly, a curve’s domain can be described implicitly by two equations f (x, y, z) = 0 and g(x, y, z) = 0, which means that just one of the variables x, y, and z varies independently. A final concept frequently used in describing functions and spaces is the codimension. Given the previous notation, the codimension of an object of topological dimension s and geometrical dimension d is the difference d − s.
To simplify implementation-related aspects, virtually all data-visualization applications fix the geometrical dimension to d = 3. This makes application code simple, uniform, yet sufficiently generic. Hence, the only dimension that actually varies in visualization datasets is the topological dimension s. In practice, s is one, two, or three, which corresponds to the curve, surface, or volumetric datasets, respectively. Since the topological dimension is the only variable one, in practice it is also called the dataset dimension. In the remainder of this book, when we talk about a dataset’s dimension, we mean its topological dimension, and implicitly assume that the geometric dimension is always three. Topological dimension is going to be important when we talk about sampled datasets and grid cells in Section 3.2.
In practice, dataset dimensions model spatial or temporal components. So, one may ask, what if we want to model a time-dependent volumetric dataset? This would require s = 4 dimensions, which would in turn ask for d = 4 geometric dimensions. Most visualization applications do not support such datasets explicitly, since they are relatively infrequent and supporting them would require working with more than three geometrical dimensions. In practice, time-dependent volumetric datasets are often represented as sequences of three-dimensional datasets where the sequence number represents the time dimension.
The function values are usually called dataset attributes. The dimensionality c of the function codomain C is also called the attribute dimension. The attribute dimension c usually ranges from 1 to 4. Attribute types are discussed later in Section 3.6. For the time being, let us assume for the simplicity of the discussion that c = 1, i.e., our function f is a real-valued function.
3.2 Sampled Data
Scientific visualization aims at displaying various properties of functions, as introduced in the previous section. However, we do not always avail of data in its continuous, functional, representation. Moreover, several operations on the data—including processing, such as filtering, simplification, denoising, or analysis, and rendering—are not easy, nor efficient, to perform on continuous data representations. For these reasons, the overwhelming majority of visualization methods works with sampled data representations, or sampled datasets. Two operations relate sampled data and continuous data:
Sampling: Given a continuous dataset, we have to be able to produce sampled data from it.
Reconstruction: Given a sampled dataset, we have to be able to recover an (approximative) version of the original continuous data.
Sampling and reconstruction are intimately connected operations. The reconstruction operation involves specifying the value of the function between its sample points, using the sample values, using a technique called interpolation. We have seen, in Chapter 2, a simple illustration of sampling and reconstruction for a two-dimensional dataset containing a scalar attribute: the surface representing the graph of a two-variable function. We sampled this dataset using two strategies, i.e., a uniform and a nonuniform sample point density. Next, we reconstructed (and rendered) the surface using a set of quadrilaterals determined by our sample points. This simple example illustrated a correlation between the sample-point density and distribution and the quality of the result of the polygon-based surface rendering. The conclusion was that the reconstruction quality is a function of the amount and distribution of sample points used. In the following, we shall analyze this relationship for the general case of a dataset represented as a d-variate, c-value function.
To be useful in practice, a sampled dataset should comply with several requirements: it should be accurate, minimal, generic, efficient, and simple [Schroeder et al. 06]. By accurate, we mean that one should be able to control the production of a sampled dataset 𝒟s from a continuous one 𝒟c such that 𝒟c can be reconstructed from 𝒟s with an arbitrarily small user-specified error, if desired. By minimal, we mean that 𝒟s contains the least number of sample points needed to ensure a reconstruction with the desired error. As we saw in Chapter 2, minimizing the sample count favors a low memory consumption and high data processing and rendering speed. By generic, we mean that we can easily replace the various data-processing operations we had for the continuous 𝒟c with equivalent counterparts for the sampled 𝒟s. By efficient, we mean that both the reconstruction operation and the data-processing operations we wish to perform on 𝒟s can be done efficiently from an algorithmic point of view. Finally, by simple, we mean that we can design a reasonably simple software implementation of both 𝒟s and the operations we want to perform on it.
Interpolation.
Let us first consider the reconstruction of a continuous approximation from the sampled data. We define reconstruction as follows: given a sampled dataset {pi, fi} consisting of a set of N sample points pi ∈ D and sample values fi ∈ C, we want to produce a continuous function
where ϕi : D → C are called basis functions or interpolation functions. In other words, we have defined the reconstruction operation using a weighted sum of a given set of basis functions ϕi, where the weights are exactly our sample values fi. Since we want that
Equation (3.3) must hold for any function f. Let us consider a function g that is overall zero, except at pj, where g(pj) = 1. Replacing the expression for g in Equation (3.3) we obtain that
Equation (3.4) is sometimes referred to as the orthogonality of basis functions. Let us now consider the constant function g(x) = 1 for any x ∈ D. Replacing g in Equation (3.3), we obtain that
then we can exactly reconstruct g everywhere in D. The property described in Equation (3.5) is called the normality of basis functions or, in some other texts, the partition of unity. Most basis functions used in practice in data approximation are both orthogonal and normal.
Grids and cells.
To reconstruct a sampled function, we can use different basis functions. These offer different trade-offs between the continuity order of the reconstruction, on the one hand, and complexity of the functions, thus the evaluation efficiency of Equation (3.3), on the other hand. To be able to explain the different choices one has for basis functions, we must now introduce the notions of sample grids and grid cells.
A grid, sometimes also called a mesh, is a subdivision of a given domain D ∈ ℝd into a collection of cells, sometimes also called elements, denoted ci. The most commonly used cells are sets of connected lines in ℝ (also called polylines), polygons in ℝ2, polyhedra in ℝ3, and so on. Moreover, the union of the cells completely covers the sample domain, i.e., ∪ici = D, and the cells are nonoverlapping, i.e., ci ∩ cj = 0, ∀i ≠ j. In other words, a grid is a tiling of the domain D with cells. The vertices of these cells are usually the sample points pi ∈ ℝd, though this is not required, as we shall see.
We can now define the simplest set of basis functions, the constant basis functions. For a grid with N cells, we define the basis functions
In Equation (3.6), the superscript 0 denotes that our basis functions are constant, i.e., of zero-order continuity. Clearly, these basis functions are orthogonal and normal. If we now imagine our sample points are inside the grid cells, e.g., at the cell centers, then these basis functions approximate a given function
Figure 3.2. Gaussian function reconstructed with constant basis functions.
Constant basis functions are simple to implement and have virtually no computational cost. Also, they work with any cell shape and in any dimension. However, constant basis functions provide a poor, staircase-like approximation
We can provide a better, more continuous reconstruction of the original function by using higher-order basis functions. The next-simplest basis functions beyond the constant ones are the linear basis functions. To use these, however, we need to make some assumptions about the cell types used in the grid. Let us consider a single quadrilateral cell c having the vertices (v1,v2,v3,v4), where υ1 = (0, 0), υ2 = (1, 0), υ3 = (1, 1), and υ4 = (1, 1), i.e., an axis-aligned square of edge size 1, with the origin as first vertex. We call this the reference cell in ℝ2. In the following, to distinguish coordinates given in the reference cell [0, 1]d from coordinates given in the function domain D ∈ ℝd, we denote the former by r1,...,rd (or r, s, t for d = 3) and the latter by x1,...,xd (or x, y, z for d = 3). Coordinates in the reference cell are also called reference coordinates. For our reference cell, we define now four local basis functions
Figure 3.3. Basis functions, interpolation, and coordinate transformations for the quad cell.
We denoted these local basis functions using capital letters (e.g., Φ) to distinguish them from global, gridwise ones, which are denoted with lowercase letters (e.g., ϕ). These basis functions are indeed orthogonal and normal. For any point (r, s) in the reference cell, we can now use these basis functions to define a function
We would like now to perform exactly the same reconstruction as previously for any quadrilateral cell c = (p1,p2,p3,p4) of some arbitrary grid, such as our height plot. In general, such cells are not axis-aligned unit squares located at the origin, but arbitrarily quadrilaterals having any position and orientation in ℝ3. How can we then use Equation (3.7) on such cells? The answer is relatively simple: For every such arbitrary quadrilateral cell c, we can define a coordinate transformation T : [0, 1]2 → ℝ3 that maps our reference cell to c.
What should such a transformation T look like? First, we want to have a simple expression for T that works for any cell type. Second, we want to map the reference cell vertices vi to the corresponding world cell pi, so T (υi) = pi. Finally, we would like to have, if possible, linear transformation T, for simplicity and computational efficiency. All this suggests that we can design T using our reference basis functions. We do this as follows: Given any cell type having n vertices p1,...,pn in ℝ3, we define the transformation T that maps from a point r, s in the reference cell coordinate system to a point x, y, z in the actual cell to be
In other words, T is a linear combination of the basis functions
To compute the inverse transformation T−1, we must invert the expression given by Equation (3.8). Since this inversion depends on the actual cell type, we shall detail the concrete expressions for T−1 later in Section 3.4.
Putting it all together.
We now have a way to reconstruct a piecewise 𝒞1 function
where cells(pi) denotes the cells that have pi as a vertex.
In other words, Equation (3.10) says that every basis function
The complete sampling and reconstruction process is illustrated schematically in Figure 3.4 for a one-dimensional signal. Sampling the continuous signal f at the points xi produces a set of samples fi. Let us organize the samples along a simple grid consisting of line segments—that is, each cell (xi, xi + 1) holds two consecutive samples (fi,fi + 1). Performing a cell-wise reconstruction, i.e., multiplying the cell samples by the global basis functions ϕi obtained from the reference basis functions
The sampling and reconstruction mechanisms described so far can be applied to more data attributes than surface geometry alone. Let us consider, for example, surface shading. Shading, as introduced in Chapter 2, can be seen as a function s : ℝ3 × ℝ3 → ℝ that yields the light intensity, given a surface point position p ∈ ℝ3 and surface normal n ∈ ℝ3. We have seen how we can approximate the geometry of a surface given a set of surface sample points, using constant or 𝒞1 basis functions, yielding a staircase (e.g., Figure 3.2) or polygonal (e.g., Figure 2.1) representation, respectively. The question is now: Can we apply the same reconstruction mechanisms for the surface shading too? And if so, what are the trade-offs?
Using our basis-function machinery, the answer is now simple. Given a polygonal surface rendering, flat shading, introduced in Chapter 2, assigns the intensity value
Table 3.1. Combinations of geometry and lighting interpolation types.
constant geometry |
bilinear geometry |
|
constant lighting |
staircase shading |
flat shading |
bilinear lighting |
— |
Gouraud shading |
3.3 Discrete Datasets
In the previous sections, we have described how to reconstruct continuous functions from sampled data provided at the vertices of the cells of a given grid. In brief, we can say that, given
a grid in terms of a set of cells defined by a set of sample points,
some sampled values at the cell centers or cell vertices,
a set of basis functions,
we can define a piecewise continuous reconstruction of the sampled signal on this grid, and work with it (e.g., compute its derivatives or draw its graph) in similar ways to what we would have done with the continuous signal the samples came from.
In Section 3.1, we defined a continuous dataset for a function f : D → C as the triplet 𝒟 = (D, C,f). In the discrete case, we replaced the function domain D by the sampling grid ({pi}, {ci}), and the continuous function f by its piecewise k-order continuous reconstruction
In the previous sections, we have shown how to replace a continuous dataset 𝒟c with its discrete counterpart 𝒟s. We have seen that this replacement means, from a mathematical point, working with a piecewise k-order continuous function
3.4 Cell Types
As explained in Section 3.2, a grid is a collection of cells ci, whose vertices are the grid sample points pi. Given some data sampled at the points pi, the cells are used to define supports for the basis functions ϕi used to interpolate the data between the sample points. Hence, the way sample points are connected to form cells is related to the domain D we want to sample, as well as the type of basis functions we want to use for the reconstruction.
The dimensionality d of the cells ci has to be the same as the topological dimension of the sampled domain D, if we want to approximate D by the union of all cells ∪ici. If D is a curve, we will use line cells. If D is a surface, we will use planar cells, such as polygons. If D is a volume (d = 3), we must use volumetric cells, such as tetrahedra. The reference coordinate systems for such cells have d dimensions, which are denoted next by r, s, and t.
In the following, we shall present the most commonly used cell types in data-visualization applications. Figure 3.5 shows, for each cell type, the shape of the cell, its vertices p0,...,pn, and the shape of the reference cell in the reference coordinate system rst. For each cell type, we shall also present the linear basis functions it supports, as well as the coordinate transformation T−1 that maps from locations x, y, z in the actual (world) cell to locations r, s, t in the reference cell.
Figure 3.5. Cell types in world (red) and reference (green) coordinate systems.
3.4.1 Vertex
We start with the simplest cell type of dimension d = 0. This cell is essentially identical to its single vertex, c = {υ1}. Following the normality property (Equation (3.5)), it follows that the vertex has a single, constant basis function:
Essentially, vertex cells are identical to, and provide nothing beyond, the sample points themselves. In practice, one does not make any distinction between sample points and vertex cells, so these cells can be seen more of a modeling abstraction.
3.4.2 Line
The cell type of the next dimension is the line segment, or line. Line cells have dimension d = 1 and two vertices, i.e., c = {υ1,υ2}. Line cells are used to interpolate along any kind of curves embedded in any dimension, e.g., planar or spatial curves. Line cells allow linear interpolation. Given the reference line cell defined by the points υ1 = 0,υ2 = 1, the two linear basis functions are
For a point p contained in a line cell (p1, p2), the transformation T−1 is simply the ratio between the lengths of the segment pp1 and the cell length, that is
3.4.3 Triangle
We move now to the next dimension, d = 2. Here, the simplest cell type is the triangle, i.e., c = {υ1,υ2,υ3}. Triangles can be used to interpolate along any kind of surfaces embedded in any dimension, e.g., planar or curved surfaces. We can do linear interpolation on triangles. Given the reference triangle cell defined by the points υ1 = (0, 0),υ2 = (1, 0),υ3 = (0, 1), the three linear basis functions are
The transformation T−1 for triangular cells is slightly more complicated than the one for line cells:
Figure 3.6 explains this intuitively: Take a triangle (p1,p2,p3), in world coordinates, and denote its area by A123. Take a point p in this triangle, and consider the areas A12 and A13 the areas of the triangles (p, p1,p2) and (p, p1, p3), respectively. These quantities are the so-called barycentric coordinates of p within our triangle. Consider now the area ratios r = A13/A123 and s = A12/A123, which are exactly the expressions in Equation 3.15. As p stays within the initial triangle, both r and s are bounded in [0, 1]. Moreover, when p = p1, we obtain (r, s) = (0, 0); when p = p2, we get (r, s) = (1, 0); and when p = p3, we get (r, s) = (0, 1). The values r and s are precisely our desired reference-cell coordinates.
Figure 3.6. Coordinate transformation Ttri explained using barycentric coordinates.
3.4.4 Quad
Another possibility to interpolate over two-dimensional surfaces is to use quadrilateral cells, or quads. We have seen examples of using quads for approximating curved surfaces since our first visualization example in Chapter 2, the height plot. Quads are defined by four vertices c = (υ0,υ1,υ2,υ3) and have four corresponding basis functions. The reference quad, defined by the points υ1 = (0, 0),υ2 = (1, 0),υ3 = (1, 1),υ4 = (0, 1), is an axis-aligned square of edge size 1. On this reference quad, the basis functions are
A natural question arises whether to use quads or triangles when interpolating on surfaces. The answer is, in most cases, a matter of implementation convenience. Some operations on triangles with linear basis functions, e.g., computing the intersection between the cell and a line or the distance from a cell to a line, are relatively simpler (and possibly faster) to implement than on equivalent quads using bilinear basis functions. However, rendering the same (large) grid using triangles might be slower than when using quads, since roughly twice as many primitives have to be sent to the graphics engine. When designing visualization software, a good practice is to use as few cell types as possible. The same is true for any data-processing operation that has to work on all the cells in a grid, e.g., computing derivatives. A good trade-off between flexibility and simplicity is to support quad cells as input data, but transform them internally into triangle cells, by dividing every quad into two triangles using one of its two diagonals. This simplifies and streamlines the overall software design, as only triangle operations have to be further implemented.
The transformation
If our actual quad cells are rectangular instead of arbitrary quads, like in a uniform or rectilinear grid, we can do better. In that case, the transformation
3.4.5 Tetrahedron
We now move to the next dimension, d = 3. Here, the simplest cell type is the tetrahedron, defined by its four vertices c = (υ1, υ2, υ3, υ4). On the reference tetrahedron defined by the points υ1 = (0, 0, 0), υ2 = (1, 0, 0), υ3 = (0, 1, 0),υ4 = (0, 0, 1), the four linear basis functions are
Given a tetrahedral cell with vertices p1,p2,p3,p4, we can deduce the transformation
Some applications use also pyramid cells and prism cells to discretize volumetric domains (see Figure 3.5). To limit the number of cell types we need to support in a concrete software implementation, pyramid and prism cells can be split into tetrahedral cells, similar to the way we split quad cells into two triangle cells.
3.4.6 Hexahedron
The next d = 3-dimensional cell type is the hexahedron, or hex, defined by its eight vertices c = (υ1,...,υ8). The reference hexahedron is the axis-aligned cube of unit edge length, with υ1 at the origin. On this cell, the eight linear basis functions are
Just as the tetrahedral cell in 3D is analogous to the triangle cell in 2D, hex cells are analogous to quad cells. Hence, the pros and cons of using hex cells instead of tetrahedra are similar to the 2D discussion involving quads and triangles, and so are the solutions. We can split hexahedral cells into six tetrahedra each, and then use only tetrahedra as 3D cell types, simplifying software implementation and maintenance.
Similar to quad cells, the transformation
3.4.7 Other Cell Types
We have presented the most commonly used cell types in d ∈ [0, 3] dimensions. As a closing remark, we mention that various authors and software packages sometimes offer more cell types, such as squares, pixels, triangle strips, polygons in 2D, and cubes and voxels in 3D. Such cells have similar (linear) interpolation functions and properties to the ones described here. Squares and pixels are essentially identical to our rectangle cell, except that they have equal-size edges, and pixels are usually thought to be aligned with the coordinate axes. Triangle strips [Shreiner et al. 03, Schroeder et al. 06] are not a different type of cell from an interpolation point of view, but just a more memory-efficient way to store sequences of triangle cells that share edges, so they can be seen as an implementation-level optimization. Cubes and voxels play the role in 3D that squares and pixels respectively have in 2D.
Finally, let us mention that the piecewise constant and linear basis functions, which we have described so far, are not the only ones used in practice. Some applications use quadratic cells, which are called so because they can support quadratic basis functions.3 Quadratic cells are the equivalents of the cells presented here, but also contain the midpoints of cell edges—shown by the red points in Figure 3.7 for several cell types. Such cells can support quadratic basis functions and provide piecewise quadratic (𝒞2), hence smoother, reconstruction of data, and are often used in numerical applications such as finite element methods [Reddy 93].4 Figure 3.7 (left) shows the quadratic interpolation
Figure 3.7. Converting quadratic cells to linear cells.
However, few visualization applications support such cells natively. In practice, visualization applications usually convert datasets that contain quadratic cells to the cell types we described previously by splitting the quadratic cells using their edge midpoints and/or cell centers into several linear cells (see Figure 3.7 for examples of several 1D and 2D cells), or even simply ignoring the extra midpoints and their associated data. Besides the added complexity, quadratic cells cannot be directly rendered by standard OpenGL, which supports only linear (Gouraud) interpolation of colors over a polygon. Although more recent additions, such as pixel shaders, allow simple implementations to perform quadratic color interpolation at every pixel of a polygon, such extensions are not yet commonly supported by visualization software.
As a general conclusion, when you have to implement and maintain a specific visualization application, it is much simpler, faster, and less error prone to support just the minimal number of cell types that are strictly needed, rather than all possible variants. In general, you should add new cell types to your application data representation only if these allow you to implement some particular visualization or data-processing algorithms much more easily and/or efficiently than cell types your software already supports.
3.5 Grid Types
Now that we have presented the cell types, we can describe the various types of grids we can construct using these cells. Many types of grids exist in the visualization practice. We shall describe the most widely used grid types: uniform grids, rectilinear grids, structured grids, and unstructured grids.
3.5.1 Uniform Grids
We start with the simplest grid type, the uniform grid. In a uniform grid, the domain D is an axis-aligned box, e.g., a line segment for d = 1, rectangle for d = 2, or rectangular solid for d = 3. We can describe this box as a set of d pairs ((m1, M1),..., (md,Md)) where (mi,Mi) ∈ ℝ2,mi < Mi are the coordinates of the box in the ith dimension. On a uniform grid, sample points pi ∈ D ⊂ ℝd are equally spaced along the d axes of the domain D. If we denote by δ1,...,δd ∈ ℝ the spacing, or sampling steps, along the d axes of D then a sample point on a uniform grid can be written as pi = (m1 + n1δ1,...,md + ndδd), where n1,...,nd ∈ ℕ. There are Ni = 1 + (Mi − mi)/δi sample points on every axis. Hence, in a uniform grid, a sample point is described by its d integer coordinates n1,...,nd. These integer coordinates are sometimes also called structured coordinates. A simple example of a uniform grid is a 2D pixel image, where every pixel pi is located by two integer coordinates. The strong regularity of the sample points in a uniform grid makes their implementation simple and economical. We can uniquely order the sample points pi in the increasing order of the indices, starting from n1 to nd, i.e.,
This numbering convention is sometimes called the lexicographic order, since it corresponds to a lexicographic (alphabetic-like) ordering of the strings nd, nd−1, ..., n1 formed by concatenating the index values. If we use this numbering convention, we do not have to store explicit sample point coordinates for uniform grids, as these can be computed from the grid sizes and samples per dimension. Storing a d-dimensional uniform grid amounts, thus, to storing 3d values. For every dimension, we store (mi,Mi,δi), i.e., the spatial extents and sampling step. Other schemes are also possible, such as storing the spatial extents and number of sample points per dimension (mi,Mi,Ni). Moreover, this regular point ordering allows us to define the grid cells implicitly by using the point indexes. For example, let us consider a d-dimensional uniform grid using boxlike cells, i.e., lines in 1D, quads in 2D, and hexahedra in 3D. In such a grid, we can compute the vertex indices v[] of the cell with index c using the function in Listing 3.1. In Listing 3.1, d is a constant indicating the grid dimension, i.e.,
1, 2, or 3. The idea of the algorithm is simple. First, we pass from the cell index c to the cell structured coordinates C[]. Next, we pass from cell structured coordinates to vertex structured coordinates. Finally, we pass from the vertex structured coordinates to the vertex indices. Here, the function lex(n1,...,nd) implements Equation (3.22). Note, however, that cell vertex indices v[] are not returned in the same order as denoted in Figure 3.5, but in order of variation of the vertex structured coordinates. Similar cell index to vertex index translations are also derivable for other cell shapes, such as triangles in 2D and tetrahedra in 3D.
Our first visualization example (see Listing 2.1) used a 2D uniform grid to store the sample values of the two-variable function to be visualized. In the code, the grid was defined by its spatial extents X_min, X_max, Y_min, and Y_max and by the number of sample points per dimension N_x and N_y.
Figure 3.8 shows two examples of uniform grids, one for a 2D rectangular and the other for a 3D box-like domain. The grid edges, that is, the cell edges whose vertices have one of the integer coordinates equal to zero, are drawn in red. The grid corners, i.e., the grid vertices whose integer coordinates are either minimal or maximal, i.e., equal to either 0 or Ni, are drawn in green.
Figure 3.8. Uniform grids. 2D rectangular domain (left) and 3D box domain (right).
The major advantages of uniform grids are their simple implementation and practically zero storage requirements. Regardless of its size, storing a d-dimensional grid itself takes 3d floating-point values, i.e., only 12d bytes of memory.5 Storing the actual sample values at the grid points takes storage proportional to the number of sample points. For example, storing one scalar value for each sample point, as in the case of our height-plot visualization, requires
3.5.2 Rectilinear Grids
However simple and efficient, uniform grids have limited modeling power. Representing non–axis-aligned domain shapes requires framing them in an axis-aligned bounding-box. Since this box is uniformly sampled, we waste memory for the sample values that fall outside the domain itself. A second problem was already illustrated in our first visualization example in Chapter 2. To accurately represent a function with a nonuniform variation rate, such as when drawing the graph of the function e−(x2 + y2), we need either to use a high sampling density on a uniform grid (as in Figure 2.1), or use a grid with nonuniform sample density (as in Figure 2.3(b)).
Rectilinear grids are a first step in this direction. They still keep the axis-aligned, matrix-like point ordering and implicit cell definition used by the uniform grids (Equation (3.22) and Listing 3.1), but they relax the constraint of equal sampling distances for a given axis. Instead, rectilinear grids allow us to define a separate sample step δij for each row of points that shares a coordinate ni in every dimension i ∈ [1,d]. A sample point on a rectilinear grid can be written as pi = (x1,...,xd), where
Figure 3.9. Rectilinear grids. 2D rectangular domain (left) and 3D box domain (right).
Implementing a rectilinear grid implies storing the grid origins (mi,Ni) and sample counts for every dimension d, as for the uniform grid. Additionally, we must store δij,i ∈ [1,d],j ∈ [1,Ni] sample steps. In total, the storage requirements are
Figure 3.9 shows two examples of rectilinear grids. The grid edges and grid corners are colored in red and green, respectively, just as for the uniform grids presented in Section 3.5.1. These grids are similar to the uniform ones shown in Figure 3.8, except that the distances δij between the sample points are now not equal along the grid axes. The 2D rectilinear grid (see Figure 3.9 (left)) is actually a slice extracted from the 3D grid (see Figure 3.9 (right)). Slicing is discussed in detail in Section 8.1.2.
3.5.3 Structured Grids
However useful, rectilinear grids do not remove the two constraints inherent to uniform grids. The sampled domain is still a rectangular box. The sample point density can be changed only one axis at a time. For our height-plot visualization of the exponential function, for example, rectilinear grids do not allow us to place more sample points only in the central peak region, where we need them. To do this, we need to allow a free placement of the sample point locations.
Structured grids serve exactly this purpose. They allow explicit placement of every sample point pi = (xi1,...,xid). The user can freely specify the coordinates xij of all points. At the same time, structured grids preserve the matrix-like ordering of the sample points, which allows an implicit cell construction from the point ordering. Intuitively, a structured grid can be seen as the free deformation of a uniform or rectilinear grid, where the points can take any spatial positions, but the cells, i.e., grid topology, stay the same. Implementing a structured grid implies storing the coordinates of all grid sample points pi and the number of points N1,...,Nd per dimension. This requires a storage space of 3
Structured grids can represent a large set of shapes. Figure 3.10 shows several examples. As shown in Figure 3.10 (left), structured grids can represent domains having a smooth, cornerless border, such as the circular domain in the image. The surface of our familiar function graph introduced in Chapter 2 is also best represented by a structured grid (see Figure 3.10 (middle)). Figure 3.10 (right) shows a 3D structured grid that has hexahedral cells. The grid edges and corners are drawn in red and green, respectively, just as for the uniform and rectilinear grid examples discussed in the previous sections.
Figure 3.10. Structured grids. Circular domain (left), curved surface (middle), and 3D volume (right). Structured grid edges and corners are drawn in red and green, respectively.
3.5.4 Unstructured Grids
Structured grids can be seen as a deformation of uniform grids, where the topological ordering of the points stays the same, but their geometrical position is allowed to vary freely. There are, however, shapes that cannot be efficiently modeled by structured grids. For example, consider a domain consisting of a square with a circular hole in the middle (see Figure 3.11). We cannot cover this domain with a structured grid, since we cannot deform a rectangle to match a rectangle with a hole—the two shapes have different topologies.6 A second limitation of all grid types described so far concerns the specification of grid cells. For uniform, rectilinear, or structured grids, cells are implicitly specified, e.g., by always connecting grid points in the same order. This can be too restrictive in many cases.
Figure 3.11. A domain consisting of a square with a hole in the middle cannot be represented by a structured grid. The domain border, consisting of two separate components, is drawn in red. Unstructured grids can easily model such shapes, whether using (a) a combination of several cell types, such as quads and triangles, or (b) a single cell type, such as, in this case, triangles.
These problems are solved by using unstructured grids. Unstructured grids are the most general and flexible grid type. They allow us to define both their sample points and cells explicitly. An unstructured grid can be modeled as a collection of sample points {pi},i ∈ [0,N] and cells {ci = (υi1,...,υiCi)}. The values υij ∈ [0,N] are called cell vertices, and refer to the sample points pvij used by the cell. A cell is thus an ordered list of sample point indices. This model allows us to define every cell separately and independently of the other cells. Also, cells of different type and even dimensionality can be freely mixed in the same grid, if desired. If cells share the same sample points as their vertices, this can be directly expressed. This last property is useful in several contexts.
First, storing an index, usually represented by an integer, is generally cheaper than storing a d-dimensional coordinate, such as d floating-point values. Second, we can process the grid geometry, i.e., the positions of the sample points pi, independently of the grid topology, i.e., the cell definitions.
Implementing an unstructured grid implies first storing the coordinates of all grid sample points pi, just as for the structured grid, and next storing all vertex indices for all cells. These indices are usually stored as integer offsets in the grid sample point list {pi}. If we use different cell types in the grid then we must either store the cell size, i.e., number of vertices, for every cell, or alternatively organize the cells in separate lists, where each list contains only cells of a given type. In practice, it is preferable to use unstructured grids containing a single cell type, as these are simpler to implement and also can lead to significantly faster application code. The costs of storing an unstructured grid depend on the types of cells used and the actual grid. For example, a grid of C d-dimensional cells with V vertices per cell and N sample points would require dN + CV values.
Figure 3.12 shows several examples of unstructured grids. The first grid (Figure 3.12 left) shows the same circular domain as in Figure 3.10, but now sampled on an unstructured grid. The domain boundary is drawn in red. The second example (Figure 3.12 middle) shows a more complex 2D domain, obtained by slicing a 3D uniform dataset containing an magnetic resonance imaging (MRI) scan with a plane. Slicing is treated in detail in Section 8.1.2. The boundary of the domain is drawn in red, for clarity. It is clearly difficult, if not impossible, to sample domains with such complex and irregular boundaries using structured grids. The unstructured grid has no problem representing such a domain, however. The third example (Figure 3.12 right) shows an unstructured grid representing a 3D surface of a bunny, with the grid edges drawn in red for clarity. As in the previous example, unstructured grids can represent such shapes with no problem.
Figure 3.12. Unstructured grids. Circle (left), head slice (middle), and 3D bunny surface (right).
3.6 Attributes
So far, we have described three of the four ingredients of a discrete dataset
As we said in the beginning of our discussion on data representation, visualization data can be modeled by some continuous or discrete function with values in a domain C ∈ ℝc. Hence, the sample values fi are c-dimensional points. In visualization, the set of sample values of a sampled dataset is usually called attribute data. Attribute data can be characterized by their dimension c, as well as the semantics of the data they represent. This gives rise to several attribute types. These are described in the following sections.
3.6.1 Scalar Attributes
Scalar attributes are c = 1 dimensional. These are represented by plain real numbers. Scalar attributes can encode various physical quantities, such as temperature, concentration, pressure, or density, or geometrical measures, such as length or height. The latter is the case for our function f : ℝ2 → ℝ visualized in our elevation plot example.
3.6.2 Vector Attributes
Vector attributes are usually c = 2 or c = 3 dimensional. Vector attributes can encode position, direction, force, or gradients of scalar functions. Usually, vectors have an orientation and a magnitude, the latter also called length or norm. If only the magnitude is relevant then we actually have a scalar, rather than a vector, attribute. If only the orientation is relevant, we talk about normalized vectors, i.e., vectors of unit length. When such vectors represent the normalized gradient of an implicit function whose graph is a 2D curve or 3D surface, these unit vectors are called also normals, as described in Chapter 2. Some authors regard vectors and normals as two different attribute types. However, in practice there is no structural difference between the two in terms of data representation. In particular, in most visualization software, normals are also stored using the same number of components as vectors, i.e., c = 2 for 2D curve normals and c = 3 for 3D surface normals, as this considerably simplifies the data storage implementation.
3.6.3 Color Attributes
Color attributes are usually c = 3 dimensional and represent the displayable colors on a computer screen. The three components of a color attribute can have different meanings, depending on the color system in use. One of the most-used color systems is the well-known RGB system, where the three color components specify the amount, or intensity, of red, green, and blue that a given color contains. Usually, these components range from 0 (no amount of a given component) to 1 (component is at full brightness). The RGB system is an additive system. That is, every color is represented as a mix of “pure” red, green, and blue colors in different amounts. Equal amounts of the three colors determine gray shades, whereas other combinations determine various hues.
RGB space.
The RGB color space is usually represented as a color cube with one corner at the origin, the diagonally opposite corner at location (1, 1, 1), and the edges aligned with the R, G, and B axes (see Figure 3.13(a)). Every point inside the cube represents one of the displayable colors. The cube’s main diagonal connecting the points (0, 0, 0) and (1, 1, 1) is the locus of all the grayscale values. Brighter colors are located closer to the cube’s outer faces (visible in the figure) whereas darker colors are located closer to the origin (hidden in the figure). RGB color representations are handy from an implementation perspective, as most graphics software, such as OpenGL and image storage formats, manipulate colors in this way.
Figure 3.13. Color-space representations. (a) RGB cube. (b) RGB hexagon. (c) HSV color wheel. (d) HSV color widget (Windows).
However, manipulating a 3D cube representation of a color space can be difficult in practice. Moreover, since we cannot see “inside” the RGB cube, it may be handy to visualize only the outer surface of the cube. If we view the cube along the main diagonal using an orthographic projection, we see a 2D hexagon, as shown in Figure 3.13(b). The cube corners, the primary colors, are shown as small spheres. As we shall see next, this represents the space of all colors having a luminance equal to one. Moreover, all other colors situated inside the RGB cube are just darker versions of these colors.
HSV space.
Another popular color representation system is the HSV system, where the three color components specify the hue, saturation, and value of a given color. The advantage of the HSV system is that it is more intuitive for the human user. Hue distinguishes between different colors of different wavelengths, such as red, yellow, and blue. Saturation represents the color “purity.” Intuitively, this can be seen as how much the hue is diluted with white or how far the color is from a primary color. A saturation of 1 corresponds to the pure, undiluted color, whereas a saturation of 0 corresponds to white. Value represents the brightness, or luminance, or a given color. A value of 0 is always black, whereas a value of 1 is the brightest color of a given hue and saturation that can be represented on a given system. For this reason, the HSV system is sometimes also called HSB, where “B” stands for brightness. The HSV color space is often represented using a color wheel (see Figure 3.13(c)). Every point p inside the color wheel represents an HSV color whose hue is given by the angle (scaled to the [0, 1] range) made by the vector p − o with the horizontal axis, where o is the wheel’s center, and whose saturation is given by the length ‖p − o‖ of the same vector. Value is not explicitly represented by the color wheel. Hence, the color wheel shows all possible hues and saturations for a given value. If value were added, as a dimension orthogonal to the color wheel space, the entire HSV space would thus create a cylinder. The HSV color wheel is quite similar to the RGB color hexagon in terms of color arrangement. Compared to the hexagon, the color wheel has the advantage that the hue parameter is mapped to a visual attribute that varies smoothly, i.e., the wheel angle. For example, a curve of constant saturation is a circle and a hexagon, respectively, in these two representations.
As we shall see next, the value (or luminance) component of an HSV color is equal to the maximum of the R, G, and B components. Hence, all colors shown by an HSV color wheel for a given value V ∈ [0, 1] are equivalent to all points on the outer faces of a cube similar to the whole RGB cube but of edge size V . Such an equal-value surface for V = 0.5 is shown in Figure 3.14 inside the RGB cube. As we shall see in Section 5.3, such a constant-value surface is called an isosurface.
Figure 3.14. Surface in the RGB cube representing colors with the constant value V = 0.5.
In practice, there exist many other visual representations of color spaces. For example, the color selector widget in the Windows operating system uses a different representation, shown in Figure 3.13(d). This widget is similar to the HSV color wheel in the sense that it maps the hue and saturation components to a 2D surface, in this case a square instead of a circle. The value component is explicitly represented as a separate color bar at the right of the colored square, and shows all colors having the H and S components specified by the selected point in the color square and values V ∈ [0, 1]. The color square can be thought of as a “cut-out” of the color wheel along the horizontal positive half-axis followed by a stretch of the resulting shape to a square. Whereas the color wheel maps equal changes of the H and S parameters to unequal changes of position on the wheel by compressing the colors with low S components into a relatively small area around the center, the color square maps the HS space isometrically to the square area. This lets users specify low S colors potentially more easily than with the color wheel.
Converting between RGB and HSV.
Since HSV color specification is often more convenient for the end user while RGB specification is required by various software, it is important to know how to map colors from the HSV to the RGB space. Listing 3.2 provides a simple C++ function that converts from the RGB to the HSV space. In this code, both RGB and HSV colors are represented as arrays of three floating-point values, with the respective color components in the natural array order. The largest RGB component gives the value, or luminance. The saturation represents the ratio between the smallest and largest RGB component (normalized to [0, 1]), i.e., how much the color differs from a grayscale value. Finally, the hue is given as the difference between the medium and smallest RGB components. A special case takes place when the saturation S equals zero. This corresponds to a grayscale, for which the hue H cannot be determined uniquely. In this case, we set H = 0 by convention.
Listing 3.3 shows the implementation of the inverse mapping from HSV to RGB. The code distinguishes six cases, which correspond to sectors of 60 degrees of the color wheel. Within each sector, saturated colors are created as a linear interpolation between two primary colors, for example red and yellow for the first sector (H ∈ [0, 1/6]), as a function of H. The final result is produced by linearly interpolating between the saturated color and white, as a function of S.
Structurally speaking, color attributes can be also seen as vector attributes. However, color is often thought of, and also implemented as, a separate attribute in most data-visualization software systems, for a number of reasons. First, operations on color are quite different than operations on general vectors. For example, it makes sense to decrease the saturation of a color, or compute its complementary color or its grayscale value, but these operations have no natural counterpart on vector data. Conversely, computing the angle between two vectors is a far more common operation on vector attributes, although it could be also performed on colors. Second, the three color components of a color attribute are usually represented as positive normalized real values in [0, 1] or, in case of a discrete color model, as 8-bit integers in [0, 255], whereas vector components usually are arbitrary real numbers. All in all, it follows that it is more natural, less confusing, and also more efficient from an implementation point of view to represent colors and vectors as different attribute types.
Color perception.
The RGB and HSV color spaces are simple and convenient tools for representing colors in a way which efficiently maps to implementation, respectively allows a simple and intuitive way for users to select colors. However, interpreting values such as hue, saturation, and brightness in terms of how humans actually perceive such color properties is a different matter. Consider luminance, for example. The HSV luminance value computed by Listing 3.2 is, indeed, low for colors that we perceive as “dark,” and respectively high for colors that we perceive as “bright.” However, this value is not proportional with what, in general, humans would perceive as being the brightness of a given color. For instance, in the HSV model, both pure blue (R = 0,G = 0,B = 1) and pure yellow (R = 1,G = 1,B = 0) have the same luminance V = 1. However, most humans will argue that they perceive pure yellow as being brighter than pure blue.
To account for such perceptual issues, and also to make the difference of colors (as represented in a color space) be as close as possible to the perceived differences of the same colors (as assessed by human subjects visualizing the respective colors on a given device), many other color representation models, or color spaces, have been designed. For example, in the sRGB color representation [ICC 14], luminance Y of a color is computed as
Intuitively, this models the fact that humans perceive green as contributing most to the brightness of a color, followed by red, and blue as being the least important contributor.
Besides RGB, HSV, and sRGB, many other color spaces, or color representations, exist. The theory and practice of color representation is an extensive subject, which is outside of the scope of this book. For more details, we refer the interested reader to an excellent field guide on color theory, representation, and practice [Stone 03].
3.6.4 Tensor Attributes
Tensor attributes are high-dimensional generalizations of vectors and matrices. Tensor data can most easily be explained by means of an example: Measuring the curvature of geometric objects.
Curvature as a tensor.
Consider first a planar curve C. For every point x0 on C, let us take a local coordinate system xy such that x is tangent to the curve and y is normal to the curve in x0. In this system, the curve can be described as y = f (x) in the neighborhood of x0, where f (x0) = 0. The curvature is then defined as ∂2f/∂x2(x0). The curvature describes how small movements on C change the curve’s normal or, alternatively, how much the curve deviates from the tangent line at a given point. The more the normal changes, the more the curve differs from its tangent line, and the higher its curvature is.
Consider now a surface S (see Figure 3.15). For every point x0 on S, take a local coordinate system xyz such that x and y are tangent to the surface and z coincides with the surface normal n in x0. Around x0, the surface can be described as z = f (x, y), with f (x0) = 0. Similar to the planar case, curvature describes how small movements along S result in changes to the surface normal or, in other words, how the surface deviates, at some given point, from its tangent plane at that point. However, there is a problem. Whereas movements around a point on a curve can happen only in two directions, and so can be described by a single number, movements around a point on a surface can happen in an infinite number of directions, so they cannot be described by a single number. To solve this problem, consider the gradient ∇f of the function f (x, y). The rate of variation of f in some direction s is
Figure 3.15. Curvature of two-dimensional surfaces.
In other words, we can link the directional derivative ∂f/∂s with the gradient by
Now consider expressing the rate of variation of the gradient in a given direction. If we now differentiate Equation 3.24 with respect to x and separately to y, and combine the two resulting expressions using vector notation, we obtain that the rate of variation of ∇f = (∂f /∂x, ∂f/∂y) in some direction s can be expressed as
where H is a 2 × 2 matrix containing the second-order partial derivatives of f (x, y) in the local coordinate system, called the Hessian of f :
From Equations (3.25–3.27), we can easily obtain that
This gives the curvature of our surface at some point x0 in any given direction s. This value is identical to the curvature of the curve 𝒞 determined by the intersection of the surface with a plane that contains the vectors n and s (see Figure 3.15) and is also known as normal curvature. However, to use Equation (3.28), we have to build our local coordinate system every time we want to compute the curvature of some point.
We can simplify this procedure as follows. Consider that our surface 𝒮 is given in global coordinates instead of local ones, for example as an implicit function f (x, y, z) = 0. We can then, after some mathematical manipulations, express the rate of change of the surface normal at the point x0 as
where H is the 3 × 3 Hessian matrix of partial derivatives of f (x, y, z) in the global coordinate system:
Using Equation (3.30), we can compute the curvature of a surface given in global coordinates, without the need of constructing a local coordinate system. For details of the derivation of Equation 3.29, see [Hartmann 99].
Summarizing, we can compute the curvature of a planar curve using its second derivative ∂2f/∂x2, and the curvature of a 3D surface in a given direction using its Hessian matrix H of partial derivatives. Given this property of the Hessian matrix, it is also called the curvature tensor of the given surface.
Tensors, vectors, and scalars.
Besides curvature, tensors can describe other physical quantities that depend on direction, such as water diffusivity or stress and strain in materials. Tensors are characterized by their rank. Scalars are tensors of rank 0. Vectors are tensors of rank 1. The Hessian curvature tensor is a rank 2 symmetric tensor, since it is expressed by a symmetric, rank 2 matrix. In general, tensors of rank k are k-dimensional real-valued arrays.
Table 3.2 summarizes the relationships between the value of a function, its gradient, and its curvature tensor, for a function of two variables f (x, y). Consider the height plot of such a function, like the half-ellipse in Figure 3.16. Given a point (x, y), the function itself, f (x, y), gives us the plot height at that point. Given also a direction s, the gradient ∇f (x, y) and Hessian tensor H(x, y) can be used to retrieve the slope and the curvature of the plot in that direction, respectively. The gradient is also the direction around (x, y) in which the slope is maximal, the maximum value being its norm. The vector ∇f (x, y)⊥, orthogonal to the gradient, is the direction in which the slope is mimimal, the minimum value being zero. Similarly, for the curvature, we can compute the directions e1 and e2 along which the curvature is maximal, respectively minimal, and the maximum and minimum curvature values λ1 and λ2, the latter being the curvatures at (x, y) of the red, respectively yellow, curves in Figure 3.16. The technique to compute these quantities, called principal component analysis, is discussed later on when detailing tensor visualization (Section 7.1).
Table 3.2. Scalars, gradient vectors, and tensors for a function f(x, y).
Quantity |
Data |
Slope |
Curvature |
Value |
f(x, y) |
– |
– |
Value in direction s |
– |
∇f(x, y) · s |
sT H(x, y)s |
Maximum direction |
– |
∇f(x, y) |
e1 |
Maximum value |
– |
‖∇f(x, y)‖ |
λ1 |
Minimum direction |
– |
∇f(x, y)⊥ |
e2 |
Minimum value |
– |
0 |
λ2 |
Figure 3.16. Scalar value, gradient vector, and curvature tensor for a function f(x, y).
In the visualization practice, two-dimensional symmetric tensors are the most common. We shall discuss tensor data in more detail in Chapter 7, when we present tensor visualization methods.
3.6.5 Non-Numerical Attributes
All the attribute types presented so far were real-valued, i.e., points in some space ℝc. The question arises whether attribute types that are not represented by real numbers can be also used. From a purely implementation perspective, the answer is simple: We can use any data type as an attribute type in a discrete dataset, by storing the data values, i.e., instances of the desired data type, at the grid points. Examples of possible non-numerical attribute types are text, images, file names, or even sound samples.
However, when doing this, an essential question immediately arises: What is the meaning of such a dataset? We defined a sampled, or discrete, dataset as the tuple
In many cases, such a relationship is hard to find, or even nonexistent for nonnumerical attributes. Conversely, if the sample values in some discrete dataset have indeed come from the sampling of a continuous function f, it is always possible to define its reconstruction
As mentioned in Section 3.1, scientific visualization, or scivis, is the branch of data visualization that works with sampled datasets, and is the main topic of this book. Purely discrete datasets, which may include points, cells, and attribute values, but no basis functions, are the domain of information visualization, or infovis. Infovis techniques are briefly covered in Chapter 11.
3.6.6 Properties of Attribute Data
The main purpose of attribute data is to allow a reconstruction
Attribute data in sampled datasets have several general properties, as follows.
Completeness.
First, attribute data, i.e., the sample values fi, must be defined for all sample points pi of a dataset 𝒟s. Indeed, if samples are missing at some points pj, it is not possible to apply Equation (3.2) over the cells that share these points pj, i.e., reconstruct
Multivariate data.
A second property of attribute data is that any cell type can contain any number of attributes, of any type, as long as these are defined for all sample points, as discussed previously. Indeed, all cell types presented in Section 3.4 define their own basis functions and thus allow reconstruction following Equation (3.2). Datasets having several attributes are usually called multivariate. If the number of attributes is large, such datasets are sometimes called high variate. The reader may have noticed that, given c real-valued data attributes defined over the same domain D, it is possible to think of these either as a single function f : D → ℝc or as c separate functions fi → ℝ for i ∈ [1..c]. In other words, we can choose whether we want to model our data as a single c-value dataset or as c one-value datasets. Of course, other combinations are possible, too. The answer to this choice is to consider all attributes that have a related meaning as a single higher-dimensional attribute, and to separate attributes with different meanings. For example, consider a two-dimensional (d = 2) image dataset, which has an RGB color c = 3 defined at every point in D ⊂ ℝ2. We could model this either as a single function f : D → ℝ3, where the value of f at some point x yields the color (R, G, B) of that point, f (x) = (R, G, B). Alternatively, we could use three functions r : D → ℝ, g : D → ℝ, and b : D → ℝ, where every function yields one color component, i.e., r(x) = R, g(x) = G, b(x) = B. The first representation is more natural, as the color components have a related meaning. Moreover, operations on color attributes must consider all color components simultaneously.
Node vs. cell attributes.
Some data visualization applications classify attribute data into node or vertex attributes and cell attributes. Node attributes are defined at the vertices of the grid cells, hence they correspond to a sampled dataset that uses linear, or higher-order, basis functions. Cell attributes are defined at the center points of the grid cells, hence they correspond to a sampled dataset that uses constant basis functions. Vertex attributes can be converted to cell attributes and conversely by resampling, as described in Section 3.9.1.
High-variate interpolation.
A final word must be said about interpolation of higher-dimensional attribute data. From the reconstruction formula (see Equation (3.2)) and the definitions of our basis functions Φi : D → ℝ for all our cell types (see Section 3.4), it follows that reconstruction of a c-dimensional function from a c-dimensional attribute data fi ∈ ℝc, where c > 1, is done by reconstructing all its c components separately from the respective components of fi.
Examples of this situation are reconstructing normal, color, and vector attributes (c = 3), and rank 2 or rank 3 tensor attributes (c = 4 and c = 9, respectively).
However, the c attribute components are sometimes related by some constraint, or invariant. We distinguish the following situations:
Normals:
For normal attributes n ∈ ℝ3, the three components nx, ny, nz are constrained to yield unit length normals, i.e.,
Depending on the choice of the basis functions, interpolating the three components nx, ny, nz separately as scalar values may not preserve the unit length property on the interpolated normal n. This is simple to verify when, e.g., using linear basis functions Φ1. There are several solutions to this problem. The easiest is to interpolate the components separately, using whatever basis functions we like, and then enforce the desired constraint on the result by normalizing it, i.e., replacing n with n/ ‖n‖. This gives good results when the sample values do not vary too strongly across a grid cell. However, enforcing the constraint on the interpolated result can fail. For example, imagine a line cell (a, b) having normal attributes na = (−1, 0, 0) and nb = (1, 0, 0) assigned to its end points. In the middle of the cell, the normal interpolated by using, e.g., linear basis functions is the null vector, which cannot be renormalized to unit length. An answer to this problem is to represent the constraint directly in the data attributes, rather than enforcing it after interpolation. For normal attribute types, this means representing 3D normals as two independent orientations, e.g., using two polar coordinates α, β, instead of using three x, y, z coordinates, which are dependent via the unit-length constraint. We can now interpolate the normal orientations α, β using the desired basis functions, to obtain the desired unit-length normal results.
Vectors:
For vector attributes, a similar situation occurs. We can use the “default” solution, that is, interpolate the three components vx, vy, vz of a vector v separately as scalar attributes. Alternatively, we can express the vector using polar coordinates (two angles and one length component), and interpolate these components. The results, as in the case of our normal attribute example described above, will be different. However, whereas for normals we had a clear invariant to enforce, which determined that the polar coordinate solution is better than per-component interpolation, the situation for vector attributes is more subtle. Here, the suitability of one interpolation scheme versus the other one depends on the nature of the continuous data that the vector dataset samples. Knowing how the continuous data actually varies allows selecting the appropriate interpolation scheme.
Colors:
Colors can be interpolated component-wise in the RGB or HSV spaces (or, for that matter, any other space in which we choose to represent them, such as sRGB). Just as for vector data, when interpolating colors it is important to know which invariants (if any) are to be preserved. For instance, linearly interpolating between pure red and pure green in the RGB space will not yield yellow along the interpolation path, while performing the same in the HSV space will do. To visualize this, consider that the interpolation path in RGB space is a straight line in the color cube shown in Figure 3.13(a). In HSV space, the same path becomes a circle arc along the color wheel in Figure 3.13(c). The optimal color space for the desired interpolation, thus, depends on the invariants we wish to maintain along the interpolation path.
Tensors:
Yet a different situation occurs for tensor attributes. Linear componentwise interpolation of the entries of the tensor matrix is the simplest, and often most used, solution. For symmetric tensors, this guarantees that the resulting value will also be a symmetric tensor [Zhukov and Barr 02]. However, in many applications, we are not directly interested in, or using, these tensor components, but derived quantities such as eigenvectors and eigenvalues (described further in Section. 7.1). Linear component-wise interpolation of tensors can create undesired variations in these derived quantities. Since a tensor can be reduced to three eigenvectors and three scalar eigenvalues, directional interpolation, similar to the proposal for vector fields based on polar coordinates outlined above, can be used. However, we now have the problem of matching two sets of three vectors each, in order to determine how the respective vectors will rotate, or vary, along the interpolation path. This matching is not trivial, and incorrect matching decisions can lead to wrong results [Kindlmann et al. 00]. For rank 2 symmetric tensors sampled on triangle grids, an effective interpolation strategy based on a combination of component-wise interpolation and eigenvector interpolation is presented in [Hotz et al. 10].
3.7 Computing Derivatives of Sampled Data
For the same polygonal surface, Gouraud shading usually produces better-looking, easier-to-understand results than flat shading. Since both flat and Gouraud shading are highly optimized nowadays by rendering engines and graphics cards, choosing between the two is not really a matter of optimizing performance. We saw in Chapter 2 that Gouraud shading requires surface normals to be available, in some way, at the polygon vertices. Also, we saw how to compute surface normals using the derivatives of the data, in case these are a continuous function. In this section, we address the problem of computing derivatives when we have a sampled dataset.
One of the requirements mentioned in the previous section for a sampled dataset 𝒟s was that it should be generic, i.e., that we can easily replace the various data-processing operations available for its continuous counterpart 𝒮c with equivalent operations on the 𝒟s. We saw, also, that computing derivatives is such an operation. Derivatives of data are needed for a wide range of processing tasks, ranging from the simple normal computation for Gouraud shading to more complex data filtering and smoothing and feature detection tasks, as described in Chapters 5 and 9.
How should we compute data derivatives in case all we have is a sampled dataset 𝒟s = (pi,ci,fi, Φi)? Since we replaced our continuous signal f with
This is a linear combination of the derivatives of the basis functions in use. However, this expression is not very convenient to use, since the actual basis functions ϕj are quite complex on a general grid (see Equation (3.10)). We can simplify the the computation of the derivatives of
We now use the derivation chain rule and obtain
to obtain
Finally, we can rewrite Equation (3.35) in a convenient and easy-to-remember matrix form, as follows:
The matrix in Equation (3.36), containing the derivatives of the reference cell coordinates ri with respect to the global coordinates xj, is called the inverse Jacobian matrix J−1 = (∂ri/∂xj)ij. As the name suggests, this matrix is indeed the inverse of the Jacobian matrix J = (∂xi/∂rj)ij. In Section 3.2, we introduced a coordinate transform T : [0, 1]d → ℝd that maps points from the reference cell to the actual cells, i.e., T (r1,...,rd) = (x1,...,xd), and its inverse T−1(x1,...,xd) = (r1,...,rd). Section 3.4 presented concrete forms of T−1 for various cell types. Using T−1, we can rewrite the inverse Jacobian as
where
Putting it all together, we get the formula for computing the partial derivatives of a sampled dataset
To use Equation (3.38) in practice, we need to evaluate the derivatives of both the reference basis functions Φk and the coordinate transform T−1 at the desired point x. We have analytic expressions or numeric methods to compute Φk and T−1 for every cell type (see Section 3.4), so the procedure is straightforward. Alternatively, we can evaluate the Jacobian matrix instead of its inverse, using the reference-cell to world-cell coordinate transforms T instead of T−1, then numerically invert J, and finally apply Equation (3.36).
We can make another important observation. For several cells described in Section 3.4, such as lines, triangles, and tetrahedra, the coordinate transformations T−1 are linear functions of the arguments xi, so their derivatives are constant. Hence, the derivatives of
Let us see what this means for a shaded polygonal surface example using, for example, triangular cells. We have seen that this is a piecewise linear (order 1) surface reconstruction. Normals are computed using the surface derivatives (see Equation (2.4)). Hence, surface normals for this polygonal surface are piecewise constant, which is exactly what flat shading does. In other words, flat shading is the mathematically “correct” shading for this polygonal surface.
For uniform and structured datasets, Equation (3.38) for computing partial derivatives is actually quite simple. For example, consider an uniform grid with cell size (δ1,...,δd), as described in Section 3.5.1. For an axis-aligned, box-like cell having the lower-left corner (p1,...,pd) and the upper-right corner (p1 + δ1,...,pd + δd), the coordinate transformation is T−1(x) = ((x1 − p1)/δ1,..., (xd − pd)/δd). Hence, the derivatives of T−1 are
Let us now consider some concrete cell and its basis functions, such as the 2D quad and its bilinear basis functions given in Section 3.4.4. By computing the derivatives ∂Φi/∂r and ∂Φi/∂s for all the four basis functions Φ1,Φ2,Φ3, and Φ4 and substituting these in Equation (3.39), we get the expression of the derivatives
where f1, f2, f3, and f4 are the sample values at the four cell vertices corresponding to basis functions Φ1, Φ2, Φ3, and Φ4. In other words, this means that the partial derivatives of
Computing derivatives of discrete datasets is a delicate process. If the dataset is noisy, the computed derivatives tend to exhibit an even stronger noise than the original data. Care must be used when interpreting information contained in the derivatives. A relatively simple method to limit these problems is to prefilter the input dataset in order to eliminate high-frequency noise, using methods such as the Laplacian smoothing described in Section 8.4. However, one must be aware that smoothing may also eliminate important information from the dataset together with the noise.
3.8 Implementation
In this section, we present an implementation of the dataset concept. The proposed implementation follows several of the requirements for datasets introduced in Section 3.2. First, it allows subsequent application code to treat all dataset types uniformly, by defining a generic Grid interface that is implemented in particular ways by the various concrete grid types. Second, it provides a reasonable balance between implementation simplicity and memory and speed efficiency, making it useful in practical applications. We stress that this implementation is not an optimal one. Higher efficiency can be achieved, both in terms of speed and smart storage schemes that minimize data duplication and copying when working with multiple datasets. For readers interested in studying an efficient and effective implementation of the dataset concept, we strongly recommend a study of [Schroeder et al. 06, Kitware, Inc. 04]. However, describing such a dataset implementation in detail, with all the design issues involved, is beyond the scope of this book.
3.8.1 Grid Implementation
We begin by defining a Grid interface that declares all operations all our grid types should support as shown in Listing 3.4.
The methods numPoints and numCells return the number of grid points and cells, respectively. In a grid, both points and cells are identified by unique integer identifiers (IDs). For uniform, rectilinear, and structured grids, these IDs usually correspond to the lexicographic ordering of the grid elements (see Equation (3.22)). Given an ID, the methods getPoint and getCell return the point coordinates and the cell vertex IDs for the point or cell specified by the passed ID argument. The method findCell returns the cell ID for the cell that contains a given location given in world coordinates, or an invalid cell ID, such as −1, if the location falls outside the grid. Finally, the method world2cell implements the inverse transformation from world coordinates to cell parametric coordinates, i.e., the transform T−1 described in Section 3.4. Specifically, world2cell(c,world,cell) maps from the 3D world coordinates world to the 3D parametric coordinates cell for the cell with the ID c. In the case of 2D or 1D cells, this method will ignore the extra coordinates. This method does not depend on the actual grid type but on the cell type, so it can be implemented in the base class. The world2cell method does not check whether the world point is indeed contained in the given cell—this is the job of the findCell method. If the world point is in the specified cell then the resulting parametric coordinates are in the [0, 1] range. If not, then at least one of the resulting parametric coordinates falls outside this range.
Since we are going to subclass Grid, we implement a dummy virtual destructor. This ensures that application code can safely delete Grid subclasses via references or pointers—the virtual specifier ensures that the appropriate subclass destructor gets invoked [Stroustrup 04]. Next, we implement the different grid types described in the previous sections by subclassing the Grid interface. To simplify the presentation, we shall consider only the case of 2D grids with quad cells. However, the implementation shown here can be easily generalized to other dimensions and cell types. For the same reasons, we skip the implementation of the attribute interpolation using basis functions and the coordinate transformations. These can be programmed by following their description in Sections 3.4 and 3.6.
Uniform grids.
Listing 3.5 shows the implementation of a UniformGrid, which is a straightforward mapping of the structure presented in Section 3.5.1.
For uniform grids, the implementation of findCell is quite simple, as the grid topology is regular. Listing 3.6 exemplifies this for 2D uniform grids. First, we compute the cell coordinates from the point world coordinates. Next, we compute the cell ID from the cell coordinates.
Rectilinear grids.
Let us now consider the rectilinear grid. We can simplify its implementation by subclassing it from UniformGrid. We can inherit most methods except getPoint(), since rectilinear and uniform grids share the same topology but not the same geometry.
In Listing 3.7, we store the 2D grid point x- and y-coordinates in float arrays d[0] and d[1] implemented using the Standard Template Library (STL) class std::vector<float>. Here and in the following, we omit the namespace qualification std:: for brevity. In actual code, either prefix all STL symbols with this namespace (e.g., write std::vector<float>) or specify using namespace std; in every file before making use of these symbols. STL provides a versatile and expressive set of programming building blocks, such as basic data containers and algorithms. Using STL to implement the various aspects of datasets massively simplifies the overall implementation and reduces the code size and the chance for errors such as memory leaks. Moreover, recent STL implementations are carefully optimized for memory and speed, so using STL containers incurs a negligible overhead as compared to classical arrays, for example.
For rectilinear grids, implementing findCell is also quite simple. First, we compute the cell coordinates from the point world coordinates, similar to the uniform grid. However, this process requires two searches (for 2D grids) in order to determine in which intervals [d[0][i], d[0][i+1]](for the x-axis) and [d[1][j], d[1][j+1]](for the y-axis) does the target point fall. These intervals can be found either via linear searches on the two axes, yielding a complexity of max (O(N1),O(N2)). Or, if we explicitly store the vectors of grid coordinates
Structured grids.
Structured grids have the same topology as uniform and rectilinear grids, but a different geometry than rectilinear grids. Hence, we derive our StructuredGrid class from UniformGrid as seen in Listing 3.8.
Since structured grids allow us to set the coordinates of every point independently, we provide a setPoint method to this end. Finally, structured grids do not have a regular geometry that can be exploited to implement a fast and simple version of the findCell method. Since the vertex coordinates can be arbitrary, findCell must perform an actual search that identifies in which cell our target point falls. Iterating over all cells is clearly too slow. An acceptable solution is to use a spatial search structure.
Several types of structures allow fast location of the closest point, or N closest points, of a given point set to a given candidate point. Such structures typically subdivide the convex hull, or the easier to compute bounding-box, of the point set hierarchically in disjoint cells that can be tested quickly. Finding the closest point(s) to the candidate point is done by traversing the spatial subdivision tree starting from the root cell (which contains the complete point set, i.e., the complete grid in our case) and following the path of cells that contain the candidate until a leaf cell is reached. If the leaf cell contains several points, a simple linear search is used to determine the closest one(s). Several spatial search structures exist. A first notable example is octrees, which use axis-aligned boxes as cells and split every cell into four (in 2D) or eight (in 3D) smaller cells of equal size. More efficient examples are kd-trees [Bentley 90, Friedman et al. 77] and box decomposition trees or bd-trees [Arya et al. 98]. These use also axis-aligned cells, but every cell is now split into two subcells. The splitting direction and the splitting line (in 2D) or plane (in 3D) are determined by a number of rules in order to balance the tree optimally. In addition, kd-trees can split a cell into a shrunken version of the cell and the space outside, which creates less-deep trees for highly clustered point sets, hence improving the search time.
A very efficient, simple-to-use open-source software package for spatial search providing kd-trees and bd-trees is the Approximate Near Neighbors (ANN) package written in C++ [Mount 06]. ANN works efficiently even beyond 3D, supports N-nearest-neighbor search with a user-specifiable N, is simple to install and call, and uses preprocessing time and space linear in the point set size and the dimension. Another open-source package that offers efficient spatial search structures such as kd-trees is the Gnu Triangulated Surface (GTS) library [GTS 13]. Figure 3.17 shows the generated kd-tree (left) and bd-tree (right) respectively for a 2D point set.
Figure 3.17. Spatial subdivision of a 2D point cloud using (a) kd-trees and (b) bd-trees.
Spatial search structures such as octrees, kd-trees, and bd-trees work best for point sets. To use them for implementing the findCell operation, we insert the complete set of grid vertices in the data structure and search for the closest vertex to the given point candidate. Once this is found, we test the cells that use this vertex to see which one actually contains the candidate. If none does, which can happen in some special configurations, we search for the second-closest vertex and repeat the procedure until the containing cell is found. An extra speed-up that can be used is based on the heuristic that, in practice, application algorithms search for cells starting from a geometrically coherent set of points. In other words, consecutive calls to findCell will return the same cell, cells that are neighbors, or other cells that are close to each other in the dataset. We can exploit this in the implementation of findCell by caching the previously returned cell and, in subsequent calls, first testing whether the searched location falls in this cell. If so, the function completes without any search. If not, the existing spatial search structure is used.
One limitation of spatial search structures is that they take relatively long to build as compared to the search time. This is not a problem in case of a grid with fixed geometry. If the grid’s geometry changes frequently compared to the number of searches, it can be faster to use a brute-force linear search instead of a spatial search structure.
Unstructured grids.
We now turn our attention to unstructured grids. Unstructured grids have the same geometry as structured grids but have a different topology. Yet, we cannot derive a UnstructuredGrid class from StructuredGrid, since the latter also inherits a regular topology from UniformGrid. Hence, we derive UnstructuredGrid directly from the Grid interface as seen in Listing 3.9.
Unstructured grids allow us to set the vertices that constitute the cells independently. To this end, we provide the setCell method. Finally, unstructured grids can use the same findCell implementation based on spatial search structures as the structured grids.
3.8.2 Attribute Data Implementation
Now that we have the grid functionality, we must implement the data attributes. In Section 3.2, we discussed two interpolation methods for attribute data, constant and linear. Both methods work in a piecewise fashion: You provide the world coordinates xj of a point in a given cell. From these, parametric cell coordinates rj are computed, using the transformation T−1, and these are used to evaluate the cell’s reference basis functions Φi weighted by the cell’s attribute values fi at its sample points (see Equation (3.9)).
An ideal interface for a dataset would let us evaluate its reconstruction function
The data attribute evaluation functionality can be easily added to our dataset class. To do this, we make some design decisions. First, we shall support only piecewise constant and linear interpolations, both for all cell types. This may seem restrictive at first. However, the overwhelming majority of visualization applications use only these interpolation methods, in practice, since they are simple to implement and quick to compute. Second, we shall enforce that attribute values fi are defined for all sample points i. As explained in Section 3.6.6, this is a natural property if we think of reconstructing signals over the whole domain. Third, we must decide which attribute types of those described in Section 3.6 we want to store in our dataset, and how many of each type. In general, we would like to store a user-defined number of each type, and also allow for user-defined attribute types. In the following sample implementation, however, we shall include only a single scalar attribute and leave the extension to several instances of several attribute types as an exercise for the reader. Finally, in the example presented here, we limit ourselves to 2D grids as in the previous code examples shown earlier in this section. Extending these examples to 3D grids is a simple exercise left to the reader.
In practice, we need attribute data for all grid types we have defined. Hence, we shall provide attribute support by adding extra code to the base class Grid, as follows.
Scalar attributes.
Piecewise constant interpolation of scalar attributes is done with the simple code in Listing 3.10.
Given a cell index c, the first method getCOScalar(int) returns the corresponding cell scalar attribute. Cell attributes are stored as a vector of floats. The second method getC0Scalar(int,float,float) performs the piecewise constant interpolation of the scalar data, given a cell index and a location within that cell. This basically amounts to returning the cell scalar value.
For piecewise linear interpolation of scalar attributes, we add the code shown in Listing 3.11 to the Grid class. The first method getC1Scalar(int) returns the sample value stored at a given grid vertex sample point, i.e., implements what is usually called vertex data. As for the cell data, this function returns a reference, so one can both read and write the vertex scalar attributes. The second method getC1Scalar(int,float*) performs the piecewise linear interpolation of the scalar data, given a cell index c and a world coordinate location p (given as a 2-float vector) within that cell, i.e., it implements Equation (3.9). In this function, the MAX_CELL_SIZE constant represents the maximum number of vertices in a cell, which is four for all 2D cells, and eight for all 3D cells discussed earlier. The method world2cell performs the transformation from world coordinates p to parametric coordinates q. Finally, the function Phi models the per-vertex basis functions of the considered cell.
Vector attributes.
Similar methods can be implemented for vector data. We next illustrate this for two-dimensional vector attributes. As for the grid type, extending this to higher-dimensional attributes is a simple exercise. For a dataset having N cells, one typically allocates a vector<float> c0_vectors having 2N float elements such that the vector v = (υx,υy) for cell i is stored at locations 2i, 2i + 1 in the container. In pure object-oriented fashion, it is tempting to consider manipulating vectors via some Vector class and thus storing vector attributes as a vector<Vector> c0_vectors instead of a plain float array. However, most visualization application implementations will not use this model, since it complicates the code significantly. Also, frequent operations on attributes, such as iteration, creating and copying attribute arrays, and reading and writing components of such an array will become slower when using the extra structuring level implied by a Vector class instance.
The implementation of nearest-neighbor interpolation for vector data in Listing 3.12 is basically identical to nearest-neighbor interpolation of scalars (Listing 3.10).
Listing 3.13 shows the implementation of linear interpolation for vector attributes. Just as for cell vector attributes, vertex vector attributes are best stored consecutively component-wise in a float array c1_vectors. The first method getC1Vector(int) returns a pointer to the first component vx of a vertex’s vector attribute. The second method getC1Vector(int,float*,float*) performs the piecewise linear interpolation of the vector data over a cell, similar to the analogous getC1Scalar method for scalar data, and returns the interpolated vector in a user-provided float array v.
Storing and interpolating other attribute types than scalars and vectors, such as normals, colors, and tensors follows the same implementation patterns as illustrated above, so we leave this also as a (slightly more complex) exercise for the interested reader.
Storing several attribute instances.
In actual production code, a dataset needs to support several instances of all the attribute types. This is required when a dataset needs to store a variable number of attributes of a given type, for instance a pressure and a temperature scalar field. This requires an interface by which the user can specify which of the instances of a certain attribute type is to be used. A simple and effective solution is to store the attribute vectors in an associative array that can be addressed by a unique name. With STL, this can be done for scalar attributes by using a map<string, vector<float>*>,for example. Given an attribute name as a string (for example, “temperature”), the map class returns a reference to the corresponding vector<float> that stores the actual data, or null if there is no such attribute. Similar associative arrays can be constructed for every attribute type.
3.9 Advanced Data Representation
In the preceding sections, we have seen how to reconstruct piecewise constant and piecewise linear functions from sampled datasets, represented on grids with various cell types, using constant and linear basis functions. For the vast majority of data visualization applications, these representations are sufficient. However, there are situations when more advanced forms of data manipulation and representation are needed. In this section, first, we will describe the task of data resampling, which is used in the process of converting information between different types of datasets that have different sample points, cells, or basis functions (Section 3.9.1). Then we will describe the process of interpolating data provided as a set of scattered points, in case we do not have cell information (Section 3.9.2).
3.9.1 Data Resampling
In our height-plot visualization example in Chapter 2, we saw several instances of reconstructing functions from a sampled dataset. We used piecewise linear interpolation for the polygonal surface itself. We used piecewise constant interpolation for the flat shading and piecewise linear interpolation for the Gouraud shading, respectively. Finally, we used piecewise constant interpolation for the surface normals. We saw in Chapter 2 that we need normal values at the polygon vertices, the vertex normals, to do the Gouraud shading of the surface. However, piecewise constant normals, i.e., the polygon normals themselves, are discontinuous at the polygon vertices—actually, over the complete polygon edges—so we cannot use them as approximations for the vertex normals. How can we compute vertex normal values from the known polygon normals?
The answer to this question is provided by a more general operation on discrete datasets, called resampling. Consider a source dataset 𝒟 = ({pi}, {ci}, {fi}, {Φi}) and a target dataset
Cell to vertex resampling.
Let us now consider a common resampling operation in data visualization: converting cell attributes (fi) to vertex attributes (
Figure 3.18. Converting cell to vertex attributes. The vertex value
A desirable property of resampling a function
That is, the integrals of the original function
where A(cj) is the area of source grid cell cj and cells(pi) are all cells that have
where A(cj) is the area of source grid cell cj and cells(pi) are all cells that have point pi as vertex. In other words, vertex data is the area-weighted average of the cell data in the cells that use a given vertex. Equation (3.42) is identical to Equation (2.8) used in Section 2.2 to compute vertex normals for Gouraud shading from the polygon normals.
Vertex to cell resampling
Using similar reasoning, we can compute the conversion formula from vertex attributes
where points(ci) denotes all points pj that are vertices of cell ci and C = |points(ci)| is the number of vertices of cell ci. In other words, cell attributes are the average of the cell’s vertex attributes.
In conclusion, we can always convert between cell attributes and vertex attributes. However, this does not mean that a dataset with cell attributes is identical to one with vertex attributes. As explained previously, cell attributes imply using piecewise constant interpolation, whereas vertex attributes imply a higher-order interpolation, such as linear. Resampling data from, e.g., cells to vertices increases the assumed continuity. If our original sampled data were indeed continuous of that order, no problem appears. However, if the original data contained, e.g., zero-order discontinuities, such as jumps or holes, resampling it to a higher-continuity grid also throws away discontinuities which might have been a feature of the data and not a sampling artifact. An example of this delicate problem is given in Section 9.4.7. In contrast, resampling from a higher continuity (e.g., vertex data) to a lower continuity (e.g., cell data) has fewer side effects—overall, the smoothness of the data decreases globally.
Subsampling and supersampling
Data resampling is not limited to converting between cell and vertex samples. Two other frequently used resampling operations are subsampling and supersampling. Subsampling is a resampling operation that reduces the number of sample points. Subsampling is useful when we are interested in optimizing the processing speed and memory demands of a visualization application by working with smaller datasets. In most cases, the subsampled dataset uses the same basis functions as the original dataset and its points are actually a subset of the original dataset points. However, this is not mandatory. After eliminating a certain number of sample points, subsampling operations can choose to redistribute the remaining points in order to obtain a better approximation of the original data in terms of the integral criterion defined earlier. Subsampling implementations can take advantage of the dataset topology. For example, on uniform, rectilinear, and structured grids, an often-used subsampling technique is to keep every kth point along every dimension and discard the remaining ones. This technique, called uniform subsampling, is simple to implement and quite effective when the original dataset is densely sampled with respect to the data variation. However, uniform subsampling makes no assumptions about the spatial data variation, so it might discard important features, such as regions where the data changes rapidly, together with less important, constant regions. Given the way most interpolation functions (such as constant and linear) reconstruct the data, a desirable property of subsampling is to keep most samples in the regions of rapid data variation and cull most samples from the regions of slow data variation. In Section 8.4, we shall present several subsampling methods applicable to both structured and unstructured datasets.
Supersampling (also called refinement) is the inverse of subsampling. Here, more sample points are created from an existing dataset. Similar to subsampling, the supersampled dataset usually includes all the original dataset points and uses the same basis functions, although this is not mandatory. Supersampling is useful in several situations when we need to manipulate or create information on a dataset at a level of detail, or scale, that is below the one captured by the sampling frequency, or density, of that dataset. In that case, we introduce more points by supersampling, and then use these to encode the desired detail information. The counterpart of uniform subsampling is uniform supersampling, which introduces k points into every cell of the original dataset. Similar to sub-sampling, an efficient supersampling implementation usually inserts extra points only in those spatial regions of the dataset where we need to further add extra information. We shall discuss several supersampling methods in Section 8.4.
3.9.2 Scattered Point Interpolation
So far, our definition of a sampled dataset has been based on a grid of cells that represent a tiling of the data domain (see Section 3.2). However, there are situations when we would like to avoid constructing or storing such a grid. Such a case occurs when you have measured some data at some given points that have a complex spatial distribution, and you have no explicit cell information that would connect the points into a tiling of some domain. A classical example is 3D surface scanning, where laser devices are used to measure the position of a set of points on the surface of a 3D object. For complex object surfaces, this process delivers a scattered 3D point set, also called a point cloud, with optional surface normal and surface color per-point attributes. Hence, all the data we have to work with are the points and their corresponding data values {pi, fi}.
For the scanning example, the data values fi are the surface normals and/or colors measured by the scanner device.
Remember, the goal of a dataset is to allow reconstructing some continuous signal from a sampled representation. In our surface scanning example, how can we reconstruct a smooth surface
Constructing a grid from scattered points
First, we can construct a grid from the point set. The way this is done depends on the meaning of the point samples pi. If these points come from the sampling of some supposedly smooth 3D surface, we can construct an unstructured grid with 2D cells, e.g., triangles, which have pi as vertices, and which approximate a smooth surface as much as possible. Several triangulation methods exist to do this. We shall describe such method in more detail in Section 8.3. Once we have this grid, we can use the constant or linear interpolation functions described earlier in this chapter.
Gridless interpolation
A second way is to avoid constructing a grid altogether. This has several advantages. Constructing unstructured grids from large, complex 3D point clouds can be a delicate process. Storing the grid can be a high burden for very large datasets. As described in Section 3.5.4, storing the cell information can double the amount of memory required in the worst case. Moreover, our point set can change in time, because it represents some moving object or because we would like to process it with geometric modeling operations such as editing, filtering, or deformations. Triangulating the point set to compute a grid every time the points change can be a very costly operation. Finally, if we have a large point set representing a complex 3D surface, visualizing it by rendering the corresponding triangulation may not be the fastest option. In such cases, we can do better by using a gridless point representation.
How can we reconstruct a continuous function from a scattered point set without recurring to an explicit grid? What we need is a set of gridless basis functions, which should respect the properties in Equations (3.4) and (3.5) as much as possible. There are several ways to construct such functions. A frequently used choice for gridless basis functions is radial basis functions, or RBFs. These are functions Φ : ℝd → ℝ+, which depend only on the distance
between the current point x = (x1,..., xd) ∈ ℝd and the origin. Moreover, RBFs Φ(r) smoothly drop from one at their origin (r = 0) to a vanishing value for large values of the distance r.
In practice, we would like to limit the effect of a basis function to its immediate neighborhood. This is both computationally efficient and ensures that a data point p contributes only to the interpolation of its immediate neighborhood, thus does not unnecessarily smooth out the information present in the data points which are far away from p. To do this, we specify a radius of influence R, or support radius, beyond which Φ is equal to zero. In this setup, a common RBF is the Gaussian function
The parameter λ ≥ 0 in Equation (3.45) controls the decay speed, or the shape, of the radial basis functions. Setting λ = 0 yields constant cylinder-shaped radial functions, which are equivalent to the constant basis functions we used for grid-based datasets. Higher values of λ yield faster-decaying functions. For d = 2, that is, on a 2D domain, the graph of such a radial basis function looks like our visualization running example introduced in Chapter 2. Besides Gaussian functions, we can use other radial basis functions too. Another popular choice are inverse distance functions defined as
In practice, we would like to make the support radius R variable for each sample point pi. The radius values Ri, i.e., values of R used for all points pi, control the influence of the sample point pi. Higher values of Ri yield smoother reconstructions
Using the reference RBF Φ and the per-point support radius values Ri, we can define the global RBFs
Performance issues
Once we have our global RBFs ϕi, we can reconstruct our surface
Spatial search structures provide efficient retrieval of the nearest neighbors and range neighbors at any given location. A good, scalable implementation of such a search structure is provided by the Approximate Nearest Neighbor (ANN) library [Mount 06]. In some sense, such a data structure plays the role of a grid, i.e., it lets us find which are the sample points that affect a given spatial region. Following the analogy, a grid cell is equivalent to the set of K nearest neighbors of a point. And just as when using grids, we must update our search data structure when our point set changes by reinserting the sample points that have moved. This is an operation that can be done efficiently [Mount 06, Clarenz et al. 04].
Shepard interpolation
Radial basis functions with compact support are an efficient and effective way to smoothly interpolate scattered point data. However, unless extra constraints are set on the support radius values Ri and function shape parameter λ, these functions are not by definition orthogonal and normal, as the constant and linear basis functions defined on a grid were. This means that the reconstruction
A simple solution that offers a second-best alternative to orthogonal and normal basis functions is to constrain the results of the interpolation to lie in the range of the sample positions pi, as follows. Given K nearest neighbors p1,..., pK of some point p, where K is a fixed value for all points, we compute
When using the inverse distance basis functions (see Equation (3.46)), this method is called Shepard’s interpolation method [Shepard 68]. Shepard’s interpolation produces results which are very similar to Gaussian RBF methods.
Figure 3.19 shows an example. The input data consists of a 2D scattered point cloud, with scalar data values fi located at the points pi. Figure 3.19(a) shows a naive visualization where each point is colored by its data value, using a blue-to-green-to-red color mapping. Hence, low value points are rendered with cold colors (blue hues); intermediate values appear as green, and high values are rendered with warm colors (orange, yellow, and red).7 Figure 3.19(a) corresponds to an interpolation
Figure 3.19. Shepard interpolation of a scalar signal from a 2D point cloud.
Figures 3.19(b–d) show the effect of Shepard interpolation computed by
The value λ is set to 5, ensuring that the Gaussian RBF Φ(r) reaches a near-zero value for a distance r = 1. The set NR(p) contains all neighbors pi located within a radius R(p) to the current point p where we want to evaluate
In Figure 3.19(b), ϵ is set to a distance equivalent to a few pixels. The colored area in the figure indicates the spatial extent over which we interpolate
Summarizing, we can represent data purely as scattered point sets {pi} carrying optional data attributes {fi}. Some visualization texts call such representations unstructured point datasets. However, if the function of a dataset is to provide a (piecewise) continuous reconstruction of its data samples using some variant of Equation (3.2), we need to specify also a choice for the basis functions Φi to have a complete dataset (pi, fi, Φi). Several such functions can be used in practice, such as radial basis functions. To efficiently perform the reconstruction, search methods are needed that return the sample points pi located in the neighborhood of a given point p.
Scattered data interpolation and approximation is an extensive subject. More information on this can be found in specialized references such as the recent book by Wendland [Wendland 06]. For a shorter, but still rigorous mathematical treatment of meshless interpolation methods, see [Belytschko et al. 96].
3.10 Conclusion
In this chapter, we have presented the fundamental issues involved in representing data for visualization applications. Visualization data is produced, in many cases, by sampling a continuous signal defined over a compact spatial domain. The combination of the signal sample values, the discrete representation of the signal domain, and the mechanisms used to reconstruct a (piecewise) continuous approximation of the signal is called a dataset. Virtually all data-processing operations in the visualization process involve using the discrete representation provided by the dataset. Hence, having generic, efficient, and accurate ways to represent and manipulate datasets is of utmost importance.
In practice, the signal domain is discretized in a grid that contains a set of cells defined by the sample points. The data samples, also called data attributes, are stored at these points, and can be of several types. The most used are numerical types which permit interpolation: scalar, vector, color, and tensor. Together with these cells, basis functions are provided for signal reconstruction. The most-used basis functions in practice are constant and linear, given the simplicity of implementation and direct support in the graphics hardware.
Several types of grids provide different trade-offs between representation flexibility and storage and computational costs. The most-used types in practice are uniform, rectilinear, structured, and unstructured grids. Finally, gridless interpolation methods exist as well. These provide signal reconstruction from the sample point avoiding the construction and storage of cells altogether. To provide this extra flexibility, special spatial search data structures are required for fast location of neighbor points.
Visualization data can also originate from different sources than sampling a continuous signal. Data attributes such as text, images, or relations are purely discrete, and often not defined on a spatial domain. Such datasets form the target of information visualization applications, and are separately described in Chapter 11.
In the following chapters, we shall see how the various aspects described in this chapter (sampling and reconstruction, interpolation and basis functions, grids and cells) will be reflected in the construction of different types of visualization applications.