Chapter 8
Domain-Modeling Techniques
DOMAIN-MODELING techniques form the last class of visualization techniques. By domain-modeling techniques, we mean those operations on datasets that modify the sampling domain representation (e.g., the grid) but not the sampled data. As we shall see, domain-modeling techniques can modify the actual values of the data attributes stored on a given grid, for example in the case of resampling the data on a different grid. However, this modification does not change the reconstructed function, so the meaning of the data attributes stays the same, even though their internal representation may change. In this chapter, we shall present a number of different modeling techniques: cutting (Section 8.1), selection (Section 8.2), constructing grids from scattered points (Section 8.3), and grid-processing techniques (Section 8.4).
8.1 Cutting
Cutting methods are domain-modeling techniques that map the data from a given source domain to a target subdomain. Consider some function f defined on a domain D. Given a domain D′ ∈ 𝔻, how can we compute the restriction of f to D′? Let us now consider that f on D is represented by a sampled “source” dataset 𝒟s = ({pi}, {ci}, {fi}, {Φi}), as described in Section 3.3. Cutting the domain D with the domain D′ means, essentially, resampling f from D to D′. This implies creating a new “target” dataset
The implementation complexity and efficiency of the cutting operation, however, depends strongly on the way we wish to define the cutting dataset. We will now present some of the most widely used variants of the cutting operation: extracting a brick, slicing, cutting with an implicit function, and generalized cutting.
8.1.1 Extracting a Brick
Extracting a brick, also called bricking or extracting a volume of interest (VOI), is a cutting operation that produces a target dataset with the same dimensionality as the source dataset. Moreover, the target grid points are a subset of the source grid points,
Figure 8.1(a) shows a brick extracted from a volumetric magnetic resonance imaging (MRI) scan of a human head. The dataset extent is shown by the wireframe and the brick surface is color-mapped with the scalar values at the respective points of the extracted brick.
Figure 8.1. (a) Brick extraction. (b) Selection of cells with scalar value above 50.
8.1.2 Slicing in Structured Datasets
Slicing is a cutting operation that is very similar to bricking. Given a uniform, rectilinear, or structured grid, we define a slice as all grid points that have one of the structured integer coordinates n1,...,nd equal. Extracting a slice can be seen as a bricking operation where the brick extent (m1,M1),..., (md, Md) is equal to the grid extent for d − 1 of the dimensions, except for the slicing axis s, where ms = Ms. Slicing a d-dimensional dataset generates a d − 1 dimensional dataset. As explained in Chapter 3 where we introduced the geometrical and topological dimensions of a dataset, this means that slicing creates cells of a lower dimension, but whose vertices are points in the same three-dimensional space as the source dataset. The most common type of slicing is extracting a set of planar cells, or a slice, from a volumetric dataset, hence the name “slicing.” Just as bricking, slicing is simple to implement: We iterate over all the sample points in the slice, in order of the structured coordinates ni, i = s that span the slice. Since our source dataset has a uniform, rectilinear, or structured grid, these integer coordinates directly correspond to d − 1-dimensional cells on the slice itself. In the target dataset, we save the sample points and the d − 1-dimensional cells on the slice plane, as well as their data attributes. In case of cell attributes, we must create these, since we also created the cells. We can do this easily by using the method of converting from point to cell attributes presented in Section 3.9.1.
The most common use of slicing is to extract 2D datasets from 3D volumes, and then visualize the extracted slices by one of the 2D visualization methods, such as color mapping or isolines for scalar data or streamlines for vector data. Figure 8.2 shows this technique applied to three slices perpendicular to the x-, y-, and z-axes of the same MRI scan uniform dataset as the one used in Figure 8.1. All three slices are taken at the middle of the respective axes of the dataset. In medical imaging, slices with these orientations are also called sagittal, axial, and coronal slices.
Figure 8.2. Slicing with planes perpendicular to the x-axis (left), y-axis (middle), and z-axis (right).
If desired, slicing can reduce the dimensionality of the dataset by more than one dimension at a time. For example, by extracting all points that share two integer coordinates from a structured dataset, we obtain a curve parallel to one of the integer coordinate axes. Slicing works also in higher dimensions. For example, extracting a time step from a 4D time-dependent volumetric dataset is equivalent to slicing the dataset with respect to the time axis.
8.1.3 Implicit Function Cutting
Slicing is a powerful and widely used visualization technique, especially when one wants to quickly browse through a large, high-dimensional dataset, without having to resort to slower, more complex visualization methods. However, slicing, as described previously, is limited to structured topology datasets. Moreover, even for such datasets, slicing limits the extraction to subsets spanned by the dataset’s integer coordinate axes.
We can generalize the slicing concept that reduces the topological dimensionality of a dataset to different subspaces than those spanned by integer coordinates in structured datasets. One way to do this is to cut an arbitrary dataset with a given lower-dimensional domain. A simple, yet powerful way to specify the cutting domain is to use implicit functions. Given some function ϕ : D → ℝ, where D is the domain of the source dataset, we define the target, or cutting, domain as all points p ∈ D for which ϕ(p) = 0. To cut the source dataset, we now proceed as follows. First, we compute a scalar dataset 𝒟cut that has the same grid as the source dataset and that evaluates ϕ. Second, we compute a contour of 𝒟cut for the value zero. As explained in Section 5.3, this yields an unstructured grid. This grid is our sampled representation of the cutting domain. Finally, we resample the source dataset attributes on this unstructured grid, using one of the available forms of interpolation (e.g., constant or linear), and we obtain the desired result.
Cutting with implicit functions generalizes the axis-aligned slicing presented in Section 8.1.2. Indeed, if we consider the implicit equation of a plane Ax + By + Cz + D = 0 with appropriate coefficients A, B, C, and D, we immediately obtain axis-aligned slice planes, at least for uniform and rectilinear grids. By changing the coefficients, we can obtain slice planes oriented at arbitrary angles with no added difficulty. Cutting through 3D volumes using more complex surfaces, such as spheres, cylinders, quadrics, or ellipsoids, is also trivial. Finally, we can now cut through all types of datasets, including unstructured ones. The added cost of implicit function slicing resides in the computation of the implicit function isosurface, which we did not need for the simpler case of slicing uniform datasets.
8.1.4 Generalized Cutting
We can generalize the cutting technique presented in the previous sections by allowing a target dataset of arbitrary definition. Instead of using a target dataset defined by structured coordinates or an implicit function, we can use an arbitrary grid, as long as its cells are contained in the source dataset. In this case, the cutting operation is identical to the last part of the implicit function cutting procedure described previously. The source attributes are interpolated at the locations of the target grid vertices (if the target dataset uses linear interpolation) or target grid cell centers (if the target dataset uses constant interpolation).
8.2 Selection
In contrast to cutting, which projects the values of a source dataset to a target domain, selection methods extract the data from a source dataset based on data properties. Cutting enforces various geometrical and/or topological properties on the target domain, since the target grid is specified by the user, but cannot explicitly enforce any properties on the data values, as these are fully specified by the source dataset. In contrast, selection explicitly specifies which data values we are interested in, but cannot enforce, in general, a certain topology and/or geometry of the shape or connectivity of the extracted dataset domain.
In the most general case, selection produces just a set of sample points and/or cells from the source dataset for which the data-based selection criterion holds. The simplest variant of selection produces a domain D′ that contains just the sample points whose data values meet the selection criterion
Here, s : D → 𝔹 is a boolean function representing the user-specified selection operation based on the attributes of the point p.
Selecting cells.
If we wish to extract cells, there are several ways to apply the selection criterion on a cell. A cell can meet the selection criterion if one of its vertices, all vertices, or its center point meet the selection criterion as defined for a point. The one-vertex criterion produces more cells, essentially selecting cells that are neighbors of the ones produced by the all-vertex selection criterion. The center point criterion is equivalent to applying the one-point selection criterion on a slightly different sampling grid. If cells are selected in the output dataset, we assume these to have the same interpolation functions as in the input dataset, since we just copied them from the input dataset. If only points are selected in the output dataset, we actually create a scattered dataset (see Section 3.9.2). Finally, the output dataset is assumed to have the same interpolation functions as the input dataset, since it is essentially just a subset of the input points and/or cells.
Since selection generally yields an arbitrary subset of points and cells from the input dataset, its output is an unstructured grid. Implementing selection is relatively simple: Depending on what we want to select (points or cells), we iterate over all the input dataset’s points or cells, apply the selection criterion, and copy the elements that pass the criterion to the output dataset, including their data attributes as well.
Thresholding, segmentation, and contouring.
Many types of selection criteria are used in visualization applications. Selection based on the scalar value matching a given target value s0 produces results related to the contouring operation, as explained later. Selection based on the scalar value being larger or equal (or smaller or equal) than a given threshold s0 produces one (or more, depending on the data monotonicity) compact subsets of the input dataset, also called threshold sets. Such an operation is also known as thresholding or segmentation. A variant of segmentation tests the scalar value against a given value range [smin,smax]. Segmentation is discussed in more detail in Chapter 9 in the context of image-processing algorithms. In addition to using the scalar values themselves, one can use their derivatives, too. In the case of scalar values that represent the luminance, or intensity, of an image, selecting data points based on the derivative values is related to edge-detection methods (see Section 9.3). Similar selection methods can be designed that use vector, color, and tensor data attributes, depending on the data at hand and application type.
Finally, let us mention that selection can also involve other properties than the data attributes of the current point. Selection methods that implement Equation (8.1) are essentially local methods, in the sense that they treat each point or cell of the dataset separately. On one hand, this is advantageous, as it lets us implement such methods simply by designing different types of local selection functions s : D → 𝔹. Moreover, such selection methods can be easily parallelized, as they treat all data points independently. However, in some cases, we are interested in selection criteria that have a quasiglobal or global nature. This means, on one hand, that the selection criterion needs to check more points together to determine whether they pass or fail the test. On the other hand, the selection will output all these points as a set instead of separately.
Such a nonlocal selection function can be described as a function s : 𝔻 → 𝔻, where s(D) = D′ ⊂ D is the result of the selection applied on the domain D. Nonlocal selection operations occur when we must enforce the connectivity of the resulting domain D′. For example, consider the operation “select all connected components
Figure 8.1(b) shows a selection of all cells from our sample MRI dataset whose scalar values are greater than or equal to 50. Data values in such an MRI scan correspond to different types of tissue. In our case, a value of 50 roughly corresponds to skin tissue, while greater values correspond to denser tissues, such as muscles or bone, hence the result of the selection shown in the image.
Selection is related to the contouring operation (see Section 5.3). Indeed, selecting all cells in a dataset whose data values are equal to a given target value τ is conceptually equivalent to producing a piecewise constant approximation of the contour at value τ. In other words, the contour is approximated by a set of cells, which gives it the blocky appearance visible in Figure 8.1(b). The marching squares and marching cubes algorithms discussed in Section 5.3 will compute the same isosurface, but use a piecewise linear approximation. In other words, the contour is approximated by a set of planes (in 3D) or lines (in 2D). Comparing Figures 8.1(b) and 5.17(d), the difference in quality of the two approximations is obvious.
8.3 Grid Construction from Scattered Points
In Section 3.9.2, we described the use of scattered point interpolation as an alternative to grids for reconstructing a piecewise continuous function from sampled data. Gridless methods are attractive when one has to manipulate datasets that contain very large numbers of unstructured point samples, which have a rather high point density. One of the uses of gridless interpolation is to render surfaces represented as 3D dense point clouds. However attractive, gridless methods have also several drawbacks: They trade the grid storage and management for storing and managing some type of spatial search structure for neighboring sample points, and they use radial basis functions that are computationally more expensive compared to piecewise linear basis functions. Moreover, most visualization software packages would require the data to be in a grid-based representation of one of the standard dataset types (see Section 3.3) before it can be processed by the available algorithms. Direct support for processing and visualizing data in gridless representations is not frequent. In such cases, constructing a grid from the scattered point set is a better alternative.
There are several methods that construct grids from scattered points. These differ in the assumptions they make about the original signal the sample points are coming from, the dimension they work in, and the type of cells they produce. We will now present several such methods.
8.3.1 Triangulation Methods
Triangulation methods are probably the most-used class of methods for constructing grids from scattered points. Given a set of points pi (sometimes also called sites), a triangulation method produces a grid (pi, ci) by generating a set of cells ci that have the sample points pi as vertices. The cells ci form a tiling of the convex hull of the point set {pi}. In other words, triangulation methods produce a grid that samples a domain D identical to the convex hull of the triangulated point set.
Delaunay triangulations.
The best-known triangulation method is the Delaunay algorithm [de Berg et al. 00]. This method generates triangular cells ci for a set of 2D points pi ∈ ℝ2 and tetrahedra for a set of 3D points pi ∈ ℝ3. A Delaunay triangulation of a point set consists of a set of triangles that covers the convex hull of the point set. An important property of a Delaunay triangulation is that no point from the input point set {pi} lies in the circumscribed circle of any triangle in the triangulation. Triangulations that obey this property are called conforming Delaunay triangulations. Given a set of scattered points with data values recorded at the point locations, using the Delaunay triangulation is the most “natural” way to create a 𝒞1, piecewise linear, interpolation of the data values over the convex hull of the points. To do this, we define piecewise linear basis functions over the triangles contained in the unstructured grid generated by the Delaunay triangulation, and use these functions to interpolate the vertex data values, as explained in Section 3.3. Figure 8.3(a) shows a Delaunay triangulation of a random point cloud containing 600 points. The point density is higher in the center, which causes the creation of smaller triangles in that area. Another example of Delauney triangulation is shown in Figure 3.12 (middle).
Figure 8.3. (a) Delaunay triangulation and (b) Voronoi diagram of a random point cloud. (c) Angle-constrained and (d) area-constrained Delaunay triangulations.
Voronoi diagrams.
For every Delaunay triangulation, there exists an associated geometric structure called a Voronoi diagram. A Voronoi diagram consists of a set of convex polygonal cells in 2D and polyhedral cells in 3D, respectively. The vertices of the Voronoi cells are the centers of the circumscribed circles of the triangles present in the associated Delaunay triangulation. The edges of the Voronoi cells are line segments contained in the lines perpendicular to, and passing through, the midpoints of the edges of the triangles present in the associated Delaunay triangulation. The centers of the Voronoi cells are the vertices of the Delaunay triangulation, i.e., the given scattered points. Figure 8.3(b) shows the Voronoi diagram of the same point set whose Delaunay triangulation is given in Figure 8.3(a).
Every location x in a Voronoi diagram is included in the Voronoi cell that has as center the closest point p in the input point set {pi}. Hence, Voronoi diagrams can be used to quickly find the closest point p from a given scattered point set to a given test location x. Note that the Voronoi cells corresponding to vertices on the convex hull of the input point set are unbounded, as they contain all points in the 2D plane that are closest to every point on the input’s convex hull. Voronoi diagrams are a natural way to create a 𝒞0, piecewise constant, interpolation of data values sampled at the scattered points, where the supports of the piecewise constant basis functions are the Voronoi cells themselves. However, Voronoi diagrams are not frequently used in practice to produce piecewise constant data approximations, since the Voronoi cells can be n-sided polygons in general, as compared to the simpler triangles of a Delaunay triangulation.
Variation of the basic techniques.
Several variations of the basic Delaunay triangulation idea exist. Angle-constrained triangulations enforce the triangle angles to lie within a given range [αmin,αmax]. For many applications, the approximation quality of a triangle grid is directly related to the triangle shapes and, consequently, to their angles. Triangles with angles close to 60 degrees provide a higher approximation quality, hence the use of angle-constrained triangles. Figure 8.3(c) shows a triangulation of the same point set as in Figure 8.3(a), where all angles lie between 20 and 140 degrees. To satisfy this constraint, 361 extra points, the Steiner points (drawn in yellow), are added to the original 600 points (drawn in red).
Area-constrained triangulations enforce a maximum triangle area and are useful in creating a sampling of a given domain with a user-specified density, which is in turn useful for representing signals with nonuniform variation with a minimal number of sample points, as explained in Chapter 3. Figure 8.3(d) shows the area-constrained triangulation of the same point set as discussed previously. Similar to the angle constraint, the minimal area constraint forces the creation of 1272 extra (Steiner) points. The original point set is colored in red, while the extra points are colored in yellow. In addition to angle and area constraints, geometric constraints can be used too. For example, the triangulation can be forced to cover the inside area of a specified convex or concave polygon whose vertices are part of the input point set, instead of covering the entire convex hull of the input point set. This triangulation variant is useful in automatically creating unstructured grids for domains with complex shapes and boundaries.
In the previous examples, we have used only the Euclidean metric to define the closest-point notion that underlies the construction of Voronoi diagrams and corresponding Delaunay triangulations. Variants of these diagrams can be obtained if we use other metrics. The additively weighted Euclidean metric, where the distance to every site pi is biased by some constant value wi, yields the Johnson-Mehl diagrams whose cell edges are hyperbolic arcs, describing the growth of crystal cells from a set of given seed sites [Okabe et al. 92]. The multiplicatively weighted Euclidean metric, where the distance to every site pi is multiplied by some constant value wi, yields the Apollonius diagrams whose cell edges are circle arcs, which are used to model plant cell growth, the tree coverage of areas in forests, and areas of best reception for radio transmitters [Sakamoto and Takagi 88]. Voronoi diagrams based on the Manhattan distance d(p, q) = ||px − qx + py − qy|| are used to model coverage areas of sites such as fire or police stations in cities where the distances are measured on a Cartesian grid.
Implementation.
Delaunay triangulation and Voronoi diagram generation are involved topics, whose details are beyond the scope of this book. For definitions of and results involving Delaunay triangulations, constrained and conforming versions thereof, and other aspects of triangular mesh generation, see the excellent survey by Bern and Eppstein [Bern and Eppstein 92].
Implementing robust, efficient, and scalable algorithms for these mesh generation methods is a complex task. Fortunately, several high-quality software implementations for Delaunay triangulation and Voronoi diagram computation are available in the open-source arena, such as the Triangle mesh generator [Shewchuk 06, Shewchuk 02]. Triangle provides a rich set of Delaunay triangulation algorithms, including the conforming and area, angle, and geometry constrained variants, as well as the computation of Voronoi diagrams, and is capable of triangulating hundreds of thousands of input points in a few seconds and with high precision on a modern PC. All examples presented in this section are computed with the Triangle software. Another high-quality open-source library providing Delaunay triangulation and Voronoi diagram operations is the Gnu Triangulated Surface Library (GTS) [GTS 13]. The interface of the GTS library is relatively more complex than that of the Triangle library. However, the GTS library offers many extra features, such as set operations on surfaces, multiresolution surface representation capabilities, and kd-trees for fast point location.
8.3.2 Surface Reconstruction and Rendering
A particular use of scattered-point interpolation is to render a 3D surface that is sampled by a point cloud. This task consists of two steps. First, we must specify which is the actual surface that the given point cloud approximates. Second, we must render this surface. We next discuss both steps briefly.
For every point in the point cloud, we assume we have three pieces of information: the point’s location pi, the surface normal ni at that location, and the average distance Ri to the neighboring points on the surface at that location. The question is: how to construct a surface
Using radial basis functions.
We approach this goal by first using 3D radial basis functions (RBFs) to construct a function
where Φ is the reference RBF in 3D and
If our basis functions
We can reconstruct the surface
Using signed distance functions.
A refinement of the previous method is to use a signed distance function. One of the first methods to do this was proposed by Hoppe et al. [Hoppe et al. 92] and works as follows. For every point pi in the point cloud, we compute a tangent plane 𝒯i that approximates our surface 𝒮 in the neighborhood of pi. The plane 𝒯i is defined by its center ci and normal ni. Computing 𝒯i can be done as follows. For every point pi in the point set, we determine a neighbor set Ni = {pj|kRi ≥ pj − pi} that contains all neighbors of pi closer than a fraction k of the support radius Ri of pi. Next, we compute 𝒯i as the plane that minimizes the sum of the squared distances Σp∈Ni d(p, 𝒯i)2 to the points in Ni. It can be shown that the center ci is the centroid of the points in Ni:
where |Ni| denotes the number of points in the neighbor set Ni, and the normal ni of 𝒯i is the eigenvector corresponding to the smallest eigenvalue of the 3 × 3 covariance matrix of the points p ∈ Ni:
with the elements
Here, pj denotes the jth component, or coordinate, of point p. Computing eigenvalues and eigenvectors of matrices was discussed in Chapter 7.
Once we have the tangent planes 𝒯i, we define a function
Finally, our desired surface
Figure 8.4. Scattered point cloud (left) and surface reconstruction with isosurface (right). (Data courtesy of H. Hoppe [Hoppe et al. 92].)
Local triangulations.
Another class of methods constructs an unstructured triangle mesh from a scattered point set by performing local 2D Delaunay triangulations [Linsen and Prautzsch 01, Clarenz et al. 04]. These methods work as follows (see also Figure 8.5). First, we compute a tangent plane 𝒯i for every point pi, using its neighbor set Ni, as described previously. Next, we project the points in Ni on 𝒯i and compute the 2D Delaunay triangulation 𝒯ri of these projections, as described earlier in this section. Next, we add to our mesh those triangles that have pi as a vertex, i.e., the triangle fan around pi. Figure 8.13 shows an application of this method for two different point clouds. Although this method is not guaranteed to produce a consistent triangle mesh, since it treats every point pi separately, it usually produces meshes with no defects such as holes or intersecting triangles. Also, this method is more memory efficient and computationally faster than the isosurface-based method first described, as it does not need to compute a volumetric distance field first. For a point cloud of P points in total and N points in an average neighbor set Ni, we need to perform P 2D Delaunay triangulations of N points each. For most point clouds, values of N ranging between 10 and 50 points give a good compromise between tangent plane stability and geometric noise elimination (which requires larger neighborhoods) and surface feature preservation (which requires smaller neighborhoods).
Figure 8.5. Mesh reconstruction from scattered points with local triangulations.
Multiple local triangulations.
The local triangulation method presented above assumes that the point cloud accurately approximates a single surface without self-intersections. However, this is not always the case. Depending on their creation process, point clouds may contain outlier noise, or points which do not actually belong to the sampling of an actual surface. Where present, outliers will corrupt the estimation of the tangent plane 𝒯i may be inaccurate. Additionally, point clouds may sample a set of intersecting surfaces. In intersection regions, we need to construct several, rather than a single, tangent plane 𝒯i.
Both above problems can be handled if we allow our algorithm to accommodate several surfaces that pass through a point, and in the same time select the most likely such surfaces. For this, we can proceed as follows. Given a point pi, we first compute its neighbor set Ni, as described previously. Next, we construct all triangles T = {ti} that joint points in Ni and have pi as vertex. In contrast to the earlier triangulation method, this is not guaranteed to produce triangles that lie in the same plane. Next, we find the most likely orientations of surface fragments approximated by Ni, by grouping triangles in T that have similar normals. This can be done by a simple bottom-up hierarchical agglomerative clustering of the dataset formed by the triangle normal vectors. The emerging triangle clusters, or patches πi, represent all surfaces that pass through Ni. If, in this process, a point pi yields only small patches πi containing just a few triangles, then pi is classified as outlier and removed. The resulting large patches πi are finally joined into whole triangulated surfaces Si by using a flood-fill process that groups overlapping patches (which share common triangles) in order of their orientation similarity.
Figure 8.6 shows the results of this method. Black dots in Figure 8.6(a) show the input point cloud which contains many noisy outliers. The several surfaces reconstructed by the triangulation method are shown half-transparent. As visible, the method can capture surfaces with complex topologies, such as the outer and inner (cavity) surfaces shown by our model. Figure 8.6(b) shows the reconstructed outer surface for the same input cloud. Each color indicates a different surface component found by the method. As visible, noise points are excluded from the reconstruction. Figure 8.6(c) shows the reconstruction of a medial point cloud, or surface skeleton, of a shoulder-blade bone, also called a scapula (see Section 9.4.8 for a detailed discussion of surface skeletons). As in the first example, noisy outliers are excluded from the reconstruction, and several smooth surface components are found. As such, this local triangulation method can serve not just for surface reconstruction, but also for point cloud segmentation and noisy outlier removal tasks.
Figure 8.6. Segmentation and reconstruction of intersecting surfaces from noisy point clouds.
Alpha shapes.
Another surface reconstruction method for unoriented point sets is provided by alpha shapes [Edelsbrunner et al. 83]. The intuition behind alpha shapes can be best described by Edelsbrunner’s carving analogy: Assume that the 2D or 3D space containing our point set P is filled up with ice cream, and we have a spherical spoon of radius α, we first carve out all ice cream we can without removing, or touching, any point. The resulting shape will have a boundary composed of spherical pieces (in 3D) or circle arcs (in 2D). The alpha shape corresponding to P is obtained by straightening out these curved parts into line segments (in 2D), respectively triangles (in 3D).
Alpha shapes offer a generalization of the notion of hull, or envelope, of a set of points representing a branching structure. Figure 8.7 illustrates this for a 2D point set P.1 Here, points are rendered in yellow, and the circles corresponding to the closest locations of our “carving spoon” to the point set P are rendered in blue. For a value of α = ∞, we obtain the result in Figure 8.7(a). The spoon is now a half-space, thus we cannot carve inside any concavity of P. Consequently, the obtained alpha shape equals the convex hull of P. As we reduce the circle radius α, increasingly less shallow concavities are captured by the alpha shape (Figure 8.7(b). However, as visible, we cannot carve between the closely-located branches in the top of the shape. Reducing α even further yields the alpha shape in Figure 8.7(c). The shape now captures most of the concavities visible in P. However, points on the perceived boundary of the shape implied by P which lie farther apart than α will also get disconnected, as the circle can pass through them to the interior. When α = 0, all space between the points can be carved out, so the corresponding alpha shape is the point set P itself.
Figure 8.7. Alpha shapes of a 2D point cloud.
Figure 8.8 shows two alpha shapes computed for two different α values from a point cloud obtained from a 3D scanner. The left image, obtained for a low α value, gives a good surface reconstruction of the point set. However, as visible in the zoom-in detail, stitches appear close to concave edges, where the α sphere cannot penetrate. The right image, obtained for a higher α value, shows clearly how these stitches fill up a larger part of the point set’s concavities.
Figure 8.8. Alpha shapes of a 3D point cloud.
Alpha shapes can be an attractive and simple-to-use surface reconstruction method for unoriented point clouds which have a regular sampling density. Robust and easy-to-use implementations of 2D and 3D alpha shapes is provided by the C++ Computational Geometry Algorithms Library (CGAL) [CGAL 13]. Note, however, that not all points in the cloud are guaranteed to be contained in the resulting surface (that is, the boundary of the alpha shape). This follows the hull metaphor provided by alpha shapes, which compute an envelope of the input points, and not a surface passing through all points. In particular, sharp convexities or concavities can get smoothed out. Also, for nonuniform point clouds, choosing a suitable α value that yields the desired surface, can be challenging. However, for situations where we are not interested to compute an exact surface approximating a point set, but an envelope thereof which captures concavities better than the convex hull, alpha shapes are a good solution. Such an example is discussed later in Section 11.4.2.
Ball pivoting.
A different technique to reconstruct 3D surfaces from oriented point clouds is ball pivoting [Bernardini et al. 99]. The principle of the method is simple: Given an oriented cloud with points X = {xi} and normals N = {ni}, and a fixed radius value ρ, we first find a seed triangle Tseed having points in X such that a 3D ball Bρ of radius ρ touching these points contains no other data points. The edges of Tseed are added to a so-called expanding front F and Tseed is added to the reconstructed mesh M, respectively. Next, we iteratively pivot a ball Bρ around each edge e = (p1, p2) ∈ F so that Bρ touches p1 and p2. If the pivoting ball hits another point q ∈ X, we distinguish three cases, as follows.
If q ∉ M, we grow F and M to include T′ = (p1, p2, q), i.e., remove e from F and add (p1, q) and (p2, q) to F and T′ to M, respectively.
If q ∈ M \ F, i.e., q is an interior vertex of M, we skip adding T′ to M (otherwise we would create a non-manifold reconstruction at q) and e becomes a boundary edge for M.
If q ∈ F, we perform the operations from case (1), and also check F to remove identical edges with opposite orientation (if any).
If the pivoting ball does not find any point q as above, then e is marked as a boundary edge for M. The process is repeated until F is empty or contains only boundary edges. At this point, we have obtained a connected unstructured mesh M describing a compact surface component. After this, we repeat the process finding another seed triangle Tseed whose vertices are not part of M. This reconstructs subsequent connected mesh components. The entire reconstruction stops when we find no Tseed any longer.
Several technical details are worth mentioning. Key to the efficiency of ball pivoting is finding points within a given radius to a seed point. This can be efficiently done using spatial search structures such as those described in Section 3.8.1. With such techniques, the complexity of ball pivoting is O(N) for N input points, which makes it attractive for reconstructing large point clouds. During the reconstruction, the normal set N is used to orient adjacent triangles consistently in M. In particular, we reject, in case 1, adding triangles T′ = (p1, p2, q) whose normal forms an obtuse angle with the point normal at q. Also, we reject seed triangles Tseed whose normal forms obtuse angles with any of the vertex normals. These constraints help that the reconstructed surface M most likely approximates the smooth surface that (X, N) sample. The radius ρ is typically set to the average interpoint distance in X. However, for highly nonuniform point clouds, reconstruction artifacts can occur (see Figure 8.9): Sparsely sampled areas will create holes in the reconstruction M, whereas concave areas where the sampled surface’s curvature is higher than 1/ρ will not be included in the reconstruction. Note that, apart from the above problems, the reconstruction has no problems in high-curvature convex areas or in nonuniformly sampled areas. These problems can be partially alleviated by running the ball reconstruction several times, at each edge, for increasing values of ρ, starting from a small ρ value.
Figure 8.9. Ball pivoting principle (sketched in 2D). Reconstruction is shown in red.
Figure 8.10(a) shows the surface reconstruction for the point cloud discussed in Figure 8.8 using ball pivoting. Since the input cloud has a rather uniform point density, ball pivoting produces a high-quality surface with no holes or artificial stitches, and which captures the details present in the input cloud. Figure 8.10(b) shows a more challenging case. Here, the point cloud contains the sampling of several intersecting manifolds—specifically, the cloud samples a surface skeleton (Section 9.4.8) of a 3D box, which is shown half-transparent in the figure. The surface skeleton consists of 13 planar manifolds that intersect each other in groups of three along straight lines. Here, ball pivoting creates several “stitches” and small gaps between the mesh reconstructions of these manifolds. This is not unexpected, since ball pivoting is not designed to treat intersecting manifolds. Also, ball pivoting can have problems with the consistent orientation of some faces, in areas that contain noisy point normals. This creates apparent twists in the resulting mesh. Although ball pivoting has such limitations, it is still one of the simplest, fastest, and thus best-known surface reconstruction algorithms in data visualization. A C implementation of ball pivoting is available in the open source mesh processing tool MeshLab [MeshLab 14].
Figure 8.10. Ball pivoting reconstruction of a manifold (a) and non-manifold (b) point cloud.
Poisson reconstruction.
We can extend the basic idea of using distance functions to reconstruct surfaces from point clouds presented earlier in this section in several ways. One such way is to incorporate normal information, as follows. Given a point cloud {xi, ni} of point locations xi and normals ni, the reconstructed surface S can be described implicitly by a function f : ℝ3 → ℝ, so that f is one inside S and zero outside. It follows that the gradient ∇f of f over f equals the inward normal of f. Hence, reconstructing S amounts to finding a scalar function f whose gradient best approximates the point normals ni at the sample points xi, i.e., finding f so that
at xi. If we apply the divergence operator to both sides of Equation 8.7, we obtain
where Δf is the Laplacian of f (see Section 8.4.4). Solving Equation 8.8 for f and a given n is fortunately easy and computationally efficient. Given the solution f, we next extract our desired surface S as the isosurface f = τ, where
Figure 8.11 shows the surface reconstruction of an oriented point cloud containing 128,338 points. The left image shows the actual point cloud. The right image shows the reconstruction result. As visible in the reconstruction image, the Poisson method can faithfully handle highly nonuniform point samples, such as the areas around the model’s eye and ear. However, regions where the cloud contains points sampling intersecting manifolds, such as the hind legs, can be altogether skipped in the reconstruction. Problems arise also in areas where the point cloud contains large holes, such as in the case of the sampling of a surface with boundaries. In such regions, the Poisson method will artificially close the surface with a smooth component. This cannot be avoided, since the reconstructed surfaces are computed via isosurfacing (see Section 5.3.2).
Figure 8.11. Poisson reconstruction of a point cloud.
Surface splatting.
The surface reconstruction methods presented so far produce a surface represented as a triangular mesh, whether via the marching cubes algorithm or directly by triangulating the point set. However, sometimes we need a simple-to-implement and fast, albeit possibly less accurate, method to directly render the surface S from the scattered points, without having to perform any explicit surface reconstruction. To explain how to do this, imagine the restrictions ψi : ℝ2 → ℝ of the 3D radial functions ϕi on our surface 𝒮 (see Figure 8.12). The 3D RBFs ϕi have compact supports on the spheres of radii Ri, so the restrictions ψi also have compact supports, which are the intersections of these spheres with the surface 𝒮. If we assume the surface to be almost flat in a neighborhood of radius Ri around every pi then the restrictions ψi are actually 2D radial basis functions defined on the surface.
Figure 8.12. Radial basis functions for surface reconstruction.
To display the surface 𝒮 described by our point set, we can render the 2D radial basis functions ψi. Just as for the 3D RBFs, the 2D RBFs ψi are transformed versions of a reference 2D radial basis function that we shall call Ψ. If Φ is a 3D Gaussian then Ψ is a 2D Gaussian, whose graph is the familiar shape shown in our elevation plot in Chapter 2 (see Figure 2.1, for example). If Φ is a 3D constant RBF, whose support2 would be a sphere of radius R, then Ψ is a 2D constant RBF, whose support is a disc of radius R. To draw our surface 𝒮, we have to draw the radial domains of ψi, which are nothing but discs of radius Ri, centered at every sample point pi and oriented in the local tangent plane to 𝒮, which is perpendicular to the surface normal ni at pi. We can do this efficiently by taking advantage of the rendering primitives offered by modern graphics hardware. First, we regularly sample the 2D radial basis function Φ on a pixel grid and store it as a 2D transparency texture T, where the value Φ = 1 maps to a totally opaque pixel T = 1 and a value Φ = 0 maps to a totally transparent pixel T = 0. The size of the texture T is taken so that T encloses the compact support of radius R of Φ.
Next, we implement the two-dimensional equivalent of Equation (8.2) by drawing the texture T mapped onto square polygonal supports centered at the sample points pi, rotated to be orthogonal to the surface normals ni and scaled to the radius Ri. For constant RBFs, this actually means drawing an opaque disc of radius Ri at every point pi. For Gaussian RBFs, this draws a set of textures of variable opacity. To sum these up as described in Equation (8.2), we turn on additive alpha blending before rendering the textures. In both cases, we can use any desired lighting model, such as the Phong model described in Chapter 2, to compute the actual color to be used on the rendered elements. We compute the surface lighting at the sample point locations pi only, using the available surface normals ni at those locations. The texture T that encodes a 2D RBF is sometimes called a splat or surfel. Hence, the previous surface reconstruction is also called splatting. We shall encounter splatting in Section 11.4.2 in a different setting when visualizing graph data.
Splats can encode other surface-information data necessary for the rendering besides the normal, color, and radius. All in all, a splat is a rendering element for a 3D surface that is analogous to the pixel, which is the rendering element for a 2D image. If we use constant RBFs, the splats are fully opaque discs stored as 2D textures. The splatting process is simple to implement and actually requires no blending support, but can easily create rendering artifacts. If we use Gaussian RBFs, the splats are variable-transparency 2D textures. However, we must ensure that the sum of the splats’ transparencies at every rendered pixel is exactly one. In case of nonuniform point distributions, this condition does not hold by default. If this condition is not enforced, this leads to pixels on the rendered surface that have transparencies below one, an artifact that is visible as half-transparent “spots” on the surface. More sophisticated definitions of RBFs and blending mechanisms that guarantee the partition of unity property are possible [Ohtake et al. 03].
Figure 8.13 shows two examples that illustrate point-based rendering and surface reconstruction from scattered points. The dragon 3D model, shown in the lower part of the image, is rendered using both disc splats (left) and a triangular mesh that has the same points as vertices (right) constructed by the tangent plane method [Clarenz et al. 04]. Clearly, the two images are very similar. The tangent plane method produces results similar in quality with ball pivoting (compare this image with Figure 8.10(a)). The dinosaur 3D model shown in Figure 8.13 (upper-left) is rendered using a reconstructed triangle mesh from a scattered point set. The quality of the triangle mesh, as well as the point set density, is visible in the detail image (Figure 8.13 (upper-right)).
Figure 8.13. Point-based rendering and surface reconstruction from scattered points.
Sphere splatting.
As mentioned in the previous section, surface splatting requires careful tuning of the shapes of the 2D basis functions ϕi to ensure that they satisfy the partition of unity on the surface 𝒮. In areas where this condition does not hold, gaps will become visible in 𝒮.
An alternative reconstruction method that avoids such problems is sphere splatting. The key idea is to compute a function
To explain the solution to this problem, note that our goal is analogous with “filling” the space between the points pi with spheres Si so that these spheres are as large as possible. Each Si is nothing but the set of points where ϕi = 1. The centers qi of these spheres and their radii Ri give then the translation, respectively the scaling components of our transforms 𝒯i. Such spheres are also called maximally inscribed, since they are contained (inscribed) inside the point cloud but do not go outside the cloud points. Consequently, each Si will be tangent to at least one point pi. For the filling analogy to hold, we must first be able to define what is “inside,” respectively “outside,” a point cloud. For oriented point clouds where each pi also has an associated normal ni, a maximally inscribed sphere Si, which is by definition tangent to pi, is said to be inside the cloud if the vector pi − qi is parallel to the normal ni.
Luckily, such maximally inscribed spheres can be readily computed: Their centers qi form the so-called medial surface, or surface skeleton, of the oriented point cloud (see Equation 9.41). Surface skeleton computation is discussed separately in Section 9.4.8. Together with the sphere centers qi, the medial surface also delivers the radii Ri = ||pi − qi|| of these spheres.
Having the sphere centers qi and radii Ri, we can now render a reconstruction of 𝒮 by rendering triangulated representations of these spheres. However, this method is computationally quite expensive, since we need finely triangulated spheres to obtain a high-quality reconstruction. A more efficient method is to render screen-space representations of such spheres, as follows. For each skeleton point qi with radius Ri, we build a viewplane-aligned quad Qi, also called a billboard. We next texture Qi scaled to (Ri, Ri, 1) with a D × D texture whose texels T(u, v) encode both depth and shading of the sphere. Here, D is set to a large value to accurately capture the rendering of such spheres, e.g., a few hundred pixels. One way to encode depth and shading is to use RGBA textures: The RGB bytes of T(u, v) encode the height at (u/D, v/D) of a half-sphere of radius 1 centered in the texture (see insets in Figure 8.14). The alpha (A) texture byte encodes the sphere color computed, e.g., with Phong shading. Next, the quads Qi are rendered centered at their corresponding skeleton points qi, using hidden surface removal. During rendering, we use as depth for each pixel the RGB color value given by its texturing, to which the depth of qi to the view plane is added. As shading for the current pixel, we use the A (alpha) texture value. The overall effect is as if we rendered half-sphere profiles, or splats, centered at the 3D locations qi. For example, the rendering of the two splats corresponding to the points qi and qj in Figure 8.14 results in the green profile, as if we effectively rendered two spheres of radii Ri and Rj centered at those points.
Figure 8.14. Sphere splatting for surface reconstruction from oriented point clouds.
Figure 8.15 illustrates this method. Figure 8.15(a) shows a polygonal rendering of the scapula bone presented earlier in Figure 8.6(c). For illustration, we consider the oriented point cloud {pi, ni} obtained from this polygonal mesh by removing all triangles, keeping their vertices pi, and computing point normals i by averaging face normals. Figure 8.15(b) shows the medial point cloud {qi} computed from this oriented point cloud. The original polygonal mesh is shown half-transparent around the cloud for illustration purposes. Figure 8.15(c) shows the surface reconstruction using sphere splatting. The result is very similar to the original polygonal mesh (Figure 8.15(a)). Zooming in, however, we can see that the reconstructed surface is nothing but a union of rendered spheres—or, in other words, a piecewise-constant approximation using spherical RBFs.
Figure 8.15. Sphere splatting for surface reconstruction from oriented point clouds.
The method produces results identical to those we would obtain by splatting volumetric spheres in 3D space. However, we work entirely in screen and z buffer space, so the computational cost is significantly smaller—basically rendering N 2D textured screen-aligned quads for a cloud having N points. When the viewpoint is changed, the splatting is performed anew in real time, so we effectively obtain a view-dependent reconstruction of 𝒮. Another attractive aspect of sphere splatting is that it works in a fully unstructured way—the splats can be rendered in any desired order to yield the same reconstruction result. Further implementation details are given in [Jalba et al. 13].
Surface reconstruction from unorganized point clouds is an advanced field, where many new algorithms emerge each year. For both theoretical and practical information on algorithms in this field, including a comprehensive set of software references, we point the reader to the book of Dey [Dey 06]. A recent survey of state-of-the-art surface reconstruction methods for various types of point clouds is given in [Chang et al. 09].
8.4 Grid-Processing Techniques
Grid-processing techniques are methods that change both the grid geometry (locations of grid sample points) and its topology (grid cells). By “grid-processing techniques,” we mean those techniques that manipulate the grid itself, and have no knowledge about data attributes sampled on that grid. Such grid-processing techniques are used in many application domains besides data visualization, such as numerical methods and simulations or computer-graphics applications. There exist a wealth of grid-processing methods. Studying all these methods in depth would be a standalone subject matter by itself. However, since grid processing is important for visualization applications, we shall present a selection of some of the most-used grid-processing methods in data visualization.
8.4.1 Geometric Transformations
Geometric transformations are domain-modeling techniques that change the position of the sample points, or grid points, but do not modify the underlying basis functions, cells, or data attributes. These are probably the simplest grid-processing techniques. Transformations in this class include affine operations such as translation, rotation, and scaling, but also nonaffine operations, such as tapering, twisting, and bending. These transformations are relatively straightforward, so we are not going to detail them further.
A second type of geometric transformation changes the position of the sample points based on the data attributes. We have encountered such techniques in the form of warping, height plots, and displacement plots, for the visualization of scalars (see Section 5.4) and vectors (see Section 6.4).
A third type of geometric transformation changes the position of the sample points based on the characteristics of the grid itself. Grid-smoothing techniques fall in this class. Given the relative importance and complexity of such techniques, we describe them separately in Section 8.4.4.
8.4.2 Grid Simplification
Many visualization applications produce large datasets that take considerable time to manipulate and store. Often, one wants to reduce the size of these datasets, yet keep the data features that are important for the task at hand. In this section, we discuss several methods that allow us to simplify the underlying grid of a dataset. By “grid simplification,” we mean situations involving grids that describe two-dimensional surfaces embedded in 3D, such as isosurfaces or polygonal models. The more general case of simplifying n-dimensional datasets by means of reducing the number of sample points (resampling) was discussed in Section 3.9.1.
We make a second assumption. The simplification criteria we shall look at here are mainly geometric, i.e., based on the shape of the grid. However, such criteria can be adapted to include information about the data attributes stored on the grid itself. An example of this technique is the progressive meshes method [Hoppe 97, Hoppe 98], which is also described later in this section.
Let us first state the generic problem. Consider a surface 𝒮 that is sampled on a geometric grid S. As we saw in Chapters 2 and 3, the quality of the approximation is dependent on the sampling rate, in case we use uniform sampling. However, uniform sampling uses too many grid points in areas of low surface curvature that can be approximated well by larger cells, i.e., fewer sample points. Given a densely sampled grid S, the aim of geometric grid simplification is to produce a surface S′ that contains fewer points than S but still provides a good approximation of S. Usually, the points of S′ are a subset of the points of S obtained by eliminating points from areas that are oversampled with respect to the desired approximation quality. However, depending on the simplification method, the points of the simplified model do not have to be a subset of the original points.
Many grid-simplification algorithms exist. The field is huge, given the applicability in computer graphics, data compression, mesh storage and transmission, shape matching, terrain visualization and rendering, and data visualization. It is, hence, impractical to aim for a comprehensive overview of grid-simplification methods in the current context. For the interested reader, we point to a number of survey papers [Luebke 01, Garland 99, Heckbert and Garland 97]. We next briefly describe a number of the main techniques used in this field, following the classification given by Luebke: triangle mesh decimation, vertex clustering, simplification envelopes, and progressive meshes [Luebke 01].
Triangle mesh decimation.
Decimation algorithms [Schroeder et al. 92] were originally designed to reduce the huge number of triangles produced by the marching-cubes method (see Section 5.3). Given a triangle mesh, the algorithm does multiple passes over the mesh vertices, checking each vertex for removal. A vertex (and its triangle fan) are removed if the removal does not change the mesh topology and if the resulting surface lies within a user-specified distance from the unsimplified surface. The hole left in the mesh is then retriangulated. Decimation continues until a user-specified reduction factor and/or some maximal error criterion are met. The decimated model contains a subset of the original mesh vertices, which is convenient for reusing the vertex information (position, normals, color, and eventual data attributes). However, this constraint can limit the decimation accuracy. A variant of the decimation algorithm can also handle topological changes [Schroeder 97].
Figure 8.16 shows two grid decimation examples produced by the latter algorithm. The bunny geometric model (see Figure 8.16(a)) is simplified to under 10% of the initial 36,000 grid points, yielding the result in Figure 8.16(b). As visible from the grid rendering, the simplification works adaptively. More triangles are kept in the high-curvature area around the bunny’s ear, marked red in Figure 8.16(b) in order to preserve the surface shape. In low-curvature areas, such as the bunny’s back, the simplification reaches its maximum. Figure 8.16(c) shows an isosurface of the skeleton of a human hand from a CT scan. The input data is a 3D uniform grid, so the marching-cubes algorithm produces an unstructured triangle mesh with almost constant point density, since the isosur-face points are always located on the edges of the input grid (see Section 5.3). The isosurface is simplified to 6536 points, i.e., to less than 2% of the original isosurface. After the simplification, surface vertex normals are computed by cell normal averaging (see Chapter 2), which lets us create a rendering of the isosurface (see Figure 8.16(d)) of comparable quality to the original.
Figure 8.16. Decimation of a surface grid. (a) Original grid with 36,000 points and (b) decimated grid with 3510 points. (c) Original isosurface with 373,000 points and (b) decimated version with 6536 points.
Vertex clustering.
Several algorithms simplify meshes by clustering (collapsing) vertices. One such algorithm [Rossignac and Borrel 93] works as follows. Every vertex gets an importance value. Vertices attached to large polygons and vertices of high curvature are more important than vertices attached to small polygons and vertices of low curvature. Next, a 3D grid is overlaid onto the mesh to be simplified. All vertices within a grid cell are collapsed to a new aggregated vertex position within that cell. The polygons whose vertices are collapsed together become degenerate and are removed.
A different clustering algorithm, called floating-cell clustering, works as follows [Low and Tan 97]. Vertices are sorted in importance order, similar to Rossignac and Borrel [Rossignac and Borrel 93]. A cell of user-specified size is centered on the most-important vertex. All vertices within the cell are collapsed to the most-important vertex, which becomes the center of the next cell, and the process repeats. Different error metrics can be used to reposition this most-important vertex in order to ensure a simplified mesh close to the original [Garland and Heckbert 97].
Simplification envelopes.
Given a surface, we construct its simplification envelopes [Cohen et al. 96]. These are two copies of the surface offset at some small distance ϵ from the original surface. Conceptually, the simplification envelopes are identical to the two components of the distance-function isosurface described in Section 8.3.2. However, we compute the envelopes here by displacing each surface vertex in the normal direction with a distance of 0.5ϵ for the outer envelope and −0.5ϵ for the inner envelope, respectively. After building the envelopes, vertices and triangles are iteratively removed from the original surface, and the holes are retriangulated, similar to the decimation method discussed previously. The simplification only occurs if the simplified surface stays between the envelopes, which guarantees that the result never deviates from the original by more than ϵ. Yet, these restrictions can sometimes limit the simplification. The implementation is quite involved, but is fortunately available to the public [Cohen et al. 07].
Progressive meshes.
A progressive mesh consists of a base mesh, created by a sequence of edge collapse operations on a polygonal mesh, and a sequence of vertex split operations [Hoppe 97]. A split is the dual of a collapse, and it replaces a vertex with two edge-connected vertices, creating an extra vertex and two extra triangles. The base mesh can be exactly transformed into the original model via splits, and the model is transformed into the base mesh via collapses. Intermediate versions correspond to progressive simplifications.
The collapses (simplification) are driven by an energy function. Different types of energies can model simplifications driven by mesh geometry, normals, and color, but also additional data attributes [Hoppe 98]. Attributes are classified as discrete, e.g., material and texture IDs, and scalar, e.g., color, normal, and texture coordinates. All edges are put into a priority queue in decreasing order of effect on the energy. The mesh is simplified by collapsing edges in this order, until topological constraints prevent further simplification. The remaining edges and triangles form the base mesh, and the (reversed) sequence of collapse operations become the split operations.
As its name suggests, this method is able to produce a hierarchy of progressively simplified meshes of high quality. Also, from a data-visualization perspective, this method is additionally attractive, as it is able to drive the simplification effort as a function of scalar and vector data attributes, not just geometric grid characteristics [Hoppe 99].
8.4.3 Grid Refinement
Grid refinement is the opposite of grid simplification. Given a coarse grid G that approximates some dataset D, refinement produces a grid G′ that also approximates D but has more sample points than the original grid G. In terms of sampling, the refined grid G′ can be seen as a supersampling of D as compared to the original grid G. Grid refinement has several uses. First, rendering a refined grid can produce a higher-quality image than rendering a coarse grid, for example, by decreasing the banding artifacts caused by the bilinear Gouraud shading. A second use of refinement is as a preprocessing step before applying other grid manipulation operations. For example, smoothing, deforming, or free-editing a refined grid gives better results as compared to performing the same operations on a coarse grid, as the refined grid has more degrees of freedom to accommodate the changes. Also, deformations such as strong stretching change an uniformly sampled grid into a nonuniformly sampled one that exhibits coarsely sampled areas. Grid refinement can be used to bring such grids back to a densely sampled version. Finally, grid refinement can be used to (partially) reverse previous grid simplification operations.
Following this sampling perspective, several grid-refinement methods exist. In case the dataset D represents a surface, we can first reconstruct the surface from the original grid G, using radial basis functions, for example, as explained earlier in this chapter, and then sample the reconstructed surface to a denser level than G to obtain the refined grid G′. Another approach goes directly from G to the refined grid G′, avoiding the cost of explicitly reconstructing the dataset D. In this latter class of methods, several variants exist. For one variant, the points of the original grid G can be a subset of the points of the refined grid G′. In the second variant, the points of G′ are obtained from those of G but do not necessarily need to be a superset of them.
Loop subdivision.
Refining a grid by adding extra sample points to the existing ones is a simple but effective strategy. A simple but effective method in this class is the Loop subdivision scheme [Loop 87, Stam 98]. For surfaces S approximated by unstructured triangular grids G, this procedure works as follows. Each triangle T ∈ G is split into four other triangles, by adding extra vertices at the midpoints of the edges of T. This delivers the topology of the new grid G′. Further on, the coordinates of the initial grid’s vertices, as well as the newly added vertices, are recomputed in order to create a grid G′ that best keeps the features of the original surface S. Here, we distinguish four cases for a vertex v (see also Figure 8.17):
Figure 8.17. Loop subdivision cases.
New internal vertex: For a new internal vertex v′ added half-way an edge e = (v1, v2), where e is not on the boundary of G, v′ is computed by linear interpolation of the coordinates of v1 and v2 and of the two vertices vA and vB of the triangles TA and TB from G which share e, as follows (see also Figure 8.17(a)):
Existing internal vertex: Existing vertices v which are not on G’s boundary (Figure 8.17(b)) are moved to a new position v′ by linearly interpolating the positions of v and of all n vertices vi which are connected to v by triangle edges in G, as follows:
The interpolation factor β is set to
New boundary vertex: New vertices v′ added on boundary edges e = (v1,v2) of G (Figure 8.17(c)) are computed simply by averaging the edge’s endpoints:
Existing boundary vertex: Finally, existing vertices v on boundary edges e = (v1, v2) (Figure 8.17(d)) are moved to new locations v′ via the following linear interpolation:
Figure 8.18 shows a grid-refinement scenario using Loop subdivision. We start with an isosurface of a human colon (Figure 8.18(a)) containing 315,600 points. As visible in the image, the isosurface exhibits strong staircasing artifacts due to some type of nearest-neighbor approximation used at some earlier point on the input data. We can render this isosurface using fewer data points and yet achieve a better visual quality. First, we simplify the surface using the decimation method discussed earlier in this chapter. The simplification yields about 13,000 points but also decreases the visual quality of the surface (see Figure 8.18(b)). If we refine the simplified surface by Loop subdivision, we obtain a surface with around 52,200 points and also a visibly better quality (see Figure 8.18(c)). The final result has a sixth of the points of the initial isosurface and also arguably a better quality and smoother appearance.
Figure 8.18. Refining an isosurface. (a) Original grid. (b) Simplified grid. (c) Refined grid. The zoomed-in insets show the grid quality.
However, Loop subdivision can produce suboptimal results. Consider the triangular grid in Figure 8.19(a), which has 16,031 triangles. This mesh is highly nonuniform, that is, contains large triangles (where the surface curvature is low) and small triangles (where the surface is highly curved, thus has more details). The grid in Figure 8.19(b) shows the result of applying two steps of Loop subdivision. This mesh contains 16,031 ∗ 4 ∗ 4 = 256,496 triangles. Subdivision nicely refines the coarse areas from the initial mesh, such as the rhinoceros’ flanks. However, since every triangle is split, this creates many very small-size triangles in detail areas. These do not bring added value, but make the dataset larger, thus more costly to store and render. This problem did not occur for the colon dataset (Figure 8.18), since the initial mesh there was an isosurface extracted from a uniform grid, hence all its triangles were of relatively similar size.
Figure 8.19. (a) Nonuniform surface mesh and (b) result of two Loop subdivision steps.
Advanced subdivision tools.
Just like grid simplification, refinement can also proceed adaptively by inserting more points where the surface varies more rapidly. As a measure of surface variation, we can use curvature or any scalar signal defined on the surface. Finally, if data attributes are defined on the original surface, these can be transferred to the refined surface by interpolation.
A very good tool for simplification and refinement of triangular and quad grids is the Yams library, implemented in C [Frey 01]. Yams can both simplify and refine surface meshes constrained by either a user-defined target size of the resulting cells, or adapting the size to a user-defined scalar signal defined on the surface. Additionally, Yams can optimize a given mesh, by moving its vertices so that the result stays close to the input mesh but the resulting triangles are close to equilateral triangles.
An alternative is the MeshLab open-source package [MeshLab 14]. MeshLab provides a wide range of surface processing techniques, including but not limited to grid simplification, refinement, remeshing, surface reconstruction from point clouds, and mesh cleaning, e.g., removing duplicate vertices, orienting normals consistently, non-manifold face removal). In contrast to Yams, MeshLab also provides visualization and interactive user interface options for editing meshes, and supports a wider range of mesh formats.
8.4.4 Grid Smoothing
Grid-smoothing methods are a separate class of geometric transformations. Grid-smoothing methods are most often used for grids that represent 2D curved surfaces embedded in ℝ3, although they can be used also for 3D volumetric grids or 1D curve grids. Recall that, for such grids, we can reconstruct a piecewise continuous surface
A simple-to-implement grid smoothing method is the Laplacian smoothing. At the core of this technique is the diffusion equation
If k is constant, this equation can be written as ∂u/∂t = kΔu. Here, Δu is called the Laplacian of u, hence the name of the diffusion process. If u(x, y, z) is a 3D scalar field, the Laplacian of u is defined as the sum of the second partial derivatives of the field with respect to the variables x, y, and z, which is in turn equal to the divergence of the gradient of the field u:
The simplest way to think about diffusion is as smoothing. Given some initial signal u0 : D for t = 0, solving the diffusion equation (Equation (8.13)) describes the smoothing of our signal u0. The time t plays the role of smoothing strength: higher t values correspond to more smoothing. The constant k describes the diffusivity, or diffusion strength, of the process: higher k values yield faster smoothing. The smoothing properties of the diffusion equation are used in many filtering applications for scalar, vector, and tensor datasets. Putting it simply, smoothing removes the high-frequency, small-scale variations in the data, such as small sharp spikes in a mesh or noise grains in an image. Smoothing is useful as a preprocessing operation for data manipulations that are noise sensitive, such as computing derivatives (see Section 3.7).
In particular, we can apply smoothing to the signal that represents the geometric coordinates of the grid. The diffusion process moves the grid vertices in the direction that smooths out the sharp features of the grid. The diffusion time t controls the amount of smoothing: the initial value u(t = 0) is equal to the input grid to be smoothed, whereas high t values produce strongly smoothed grids.
Using the definitions of the derivatives on a discrete grid presented in Chapter 3, it can be shown that the discrete form of Equation (8.13) on a grid containing the sample points {pi} is
Here, (qj)j=1..N denote the neighbors of pi, i.e., those grid points connected to pi by edges. Equation (8.15) is solved iteratively, the superscript n denoting the iteration number. We start with the initial grid points
Figure 8.20. Laplacian smoothing principle for (a) 2D and (b) 3D geometries.
Figure 8.21 shows the Laplacian smoothing in action on an isosurface showing a femur bone. The number of iterations combined with the diffusivity k can be used to control the scale of details one wants to remove from a given dataset. In Figure 8.21(c), small-scale irregularities of the original surface (Figure 8.21(a)) are smoothed out in approximately 100 iterations. If more iterations are applied, larger-scale details would disappear, which is undesirable in this case. Finally, Figure 8.21(d) compares the initial object (green, semitransparent) with the smoothed one (gray, opaque) after 10,000 iterations, by rendering both objects overlapped in the same image. The smoothed object is visibly smaller than the original in the green (convex) regions and larger in the gray (concave) ones.
Figure 8.21. Laplacian smoothing of an isosurface. (a) Original surface. (b) Surface curvature. (c) Smoothed surface. (d) Comparison of original and smoothed surfaces.
Like other low-pass filtering methods, Laplacian smoothing removes both noise and small-scale surface details. An undesired effect is that sharp surface creases, such as the edge of the femur cross section in Figure 8.21, also get smoothed out. Several advanced smoothing methods exist that prevent this. For example, one can use weighted diffusion in Equation (8.13). The diffusion coefficient k is set to be inversely proportional to the surface curvature instead of a using constant value. Weighted diffusion has the effect of strongly flattening small-scale noise but maintaining sharp details such as edges and creases. For this to work, we must use a smoothed (e.g., low-pass filtered) version of the surface curvature so that small-scale noise details are not visible in this signal. Figure 8.21(b) shows the smoothed surface curvature for the bone model color-coded from blue (flat) to red (maximally curved). The edge-preserving effect can be further improved by forcing the diffusion process to be strong in the direction of the surface edges and weak across the edges. Anisotropic techniques that perform this type of smoothing are an advanced subject and are described further in specialized literature [Perona and Malik 90, Weickert 98, Bajaj and Xu 03].
Grid-smoothing and grid-refinement methods can be effectively combined to produce high-quality grids. A simple way to do this is to first refine the grid, in order to create enough sample points so that the smoothing can act on a small spatial scale, and next perform a number of smoothing steps until the desired surface quality is reached. More efficient implementations can alternate smoothing and refinement steps and also introduce grid-decimation steps, all in order to redistribute the points to the optimal locations on the surface in order to achieve the desired quality.
Grid processing is an extensive topic of research. Polygon mesh processing, a key ingredient to all domain processing algorithms presented in this chapter, is covered in great detail by Botsch et al. in [Botsch et al. 10]. A more general resource for the fundamental computational geometry algorithms that underlie all the domain processing algorithms outlined above is the book of Overmars et al. [de Berg et al. 00].
8.5 Conclusion
Grids are a fundamental element of visualization datasets. They provide a discrete representation of a compact spatial domain on which the signals to be visualized are sampled. Grids can vary in several respects: the type of cells, the regularity of the discretization, and the types of basis functions used to perform interpolation on the grid.
In this chapter, we have presented a number of fundamental methods for grid manipulation and processing. Since grids are used in visualization as a representation of the underlying domain of a signal, these methods are referred to as domain-modeling techniques. Cutting techniques extract a lower-dimensional domain from a higher-dimensional one, such as when slicing a volume with a surface, or a subdomain with the same dimensionality, such as in the case of bricking. Selection techniques extract a set of cells or points from a dataset, based on data properties. These techniques are useful when we are interested in extracting a specific subset of interest from a larger dataset.
Grids can be constructed from scattered points using triangulation techniques. An alternative is to use gridless techniques, e.g., using radial basis functions. Triangulation has the advantage of producing a cell-based grid, which can be further used by all grid-based visualization methods. Gridless methods do not require the usually complex triangulation step, require less memory because they do not explicitly store cells, offer fast surface interpolation (rendering) using splatting, but bear additional computational costs. Either of these techniques can be used when we need a continuous domain representation, such as a surface, and all we have is a set of scattered points.
Grids can be processed by a variety of operations, such as simplification, refinement, and smoothing. These operations are especially useful when the grid itself is the visualization target, such as when it represents a surface of interest, e.g., produced by contouring techniques. All in all, domain modeling techniques are an indispensable element of the visualization process, as grids are an indispensable ingredient of datasets.