Overview

What I learned:

• mathematical foundation and properties of natural cubic splines and bicubic surfaces
• variational formulation of splines and surfaces
• finite element discretization of variational splines and surfaces
• construction of variational splines and surfaces using multi-resolution Gauss-Seidel relaxation
• Poisson interpolation using multi-resolution GSR

Extras:

• extension of conventional variational splines and surfaces to generalized parametric versions
• integration of Poisson interpolation with variational surface construction (with comparison)
• deviation from circle for a generalized beam spline
• deviation from sphere for a generalized thin plate surface with and without freezing
• comparison of generalized thin-plate and pseudo-bicubic surfaces

Natural Cubic Splines and Bicubic Surfaces

Figure 1. Variational spline examples.

Natural cubic splines arise from the historical use of flexible strips of metal, also called splines, to draw smooth shapes. The shape of the physical spline is controlled by fixing its location at a number of points. Mathematically, the spline is modelled as a function f which is defined on an interval [a, b] and characterized by two requirements. First, there is a finite set of points in [a, b] at which the value of f is fixed and which also divide [a, b] into subintervals. Second, f is a cubic in each of the subintervals and its first and second derivative are continuous everywhere in [a, b]. Natural cubic splines have a number of powerful properties [1]:

Uniqueness

The spline f is unique assuming some boundary conditions at a,b are satisfied. In my understanding, non-uniqueness can only be caused by undersampling.

Best Approximation

If g is a known function (with some restrictions on its derivatives), then its spline approximation f minimizes:

Rapid Convergence

If g is a known function with a continuous q-th derivative (q < 5), then the approximating spline f converges to g as quickly as the q-th power of the length of the longest subinterval of [a, b] approaches 0.

Minimum Norm

Under some boundary conditions, the spline function f is the unique admissible function that minimizes:

These properties make natural cubic splines useful in many applications such as curve fitting, numeric integration, and numerically solving differential and integral equations [1]. Furthemore, the definition of natural cubic splines can be generalized in a number of ways and also extended to bicubic surfaces. In many cases, the above properties all have appropriate generalizations [1]. In particular, the Minimum Norm Property for bicubic splines involves the minimization of:

Variational Models of Splines and Surfaces

The variational spline model is motivated by the Minimum Norm Property of natural cubic splines. In this model, the mathematical spline or surface is an unknown function that minimizes an energy functional which corresponds to a physical bending energy. This principle provides for a lot of flexibility in controlling the behaviour of the generated spline. For example, increasing the energy when the spline fails to pass through a given set of control points can force the spline to interpolate those points. Similarly, the smoothness of the spline can be controlled with energy terms involving the spline's derivatives [4]. For an example, Figure 1b shows a mixture of 25% string and 75% beam (described below) behaviour.

Figure 2. Variational surface examples.

In the two-dimensional setting, a variational spline can be represented as the graph of a univariate function f. Then the following functional F forces f to achieve string-like behaviour (Figure 1a):

A spline will behave like a flexible beam (Figure 1c) if its energy is modelled to be proportional to that given by the following functional. These types of splines correspond to natural cubic splines [7].

The principle extends readily to surfaces in three dimensions, with the string and the beam generalizing to become the membrane and the thin plate, respectively (Figure 2). The associated functionals are the following [4]:

Finite Element Discretization

The sought function f corresponding to a variational spline or surface can be discretized on a uniform grid, resulting in a set of samples called nodal variables [4]. One way to compute these samples is to use the finite element method to also discretize the corresponding energy functional and obtain an approximation for the total energy of the spline in terms of the nodal variables. Solving the optimization problem of locating the minimum of that energy equation is equivalent to computing the samples of f. This optimization problem becomes much more tractable if the total energy expression is broken up into local terms that involve only a few neighbouring samples. This is possible because the contribution of each finite element used in the total energy equation depends only on the nodal variables referenced in the appropriate discrete approximations to the first and second derivatives of f.

Figure 3. (a) Nodal variables referenced for a given finite element (shaded);
(b) Terms which locally depend on a given nodal variable (corresponding to the bold point).

There are some additional considerations in setting up the finite element discretization of the energy functional. First, a lot of care has to be taken to avoid redundant calculations and asymmetry along the borders. Second, leveraging the locality of reference described above becomes a major performance improvement for Gauss-Seidel relaxation (the method I have chosen for solving the optimization problem), so it becomes necessary to formulate a discretized local energy and keep the total and local versions compatible. Figure 3 illustrates the relationship between the finite elements and nodal variables for the case of discretizing a variational surface problem that uses a thin plate functional without a mixed term. Figure 4 continues the example and shows that the setup in Figure 3b results in symmetrical calculations along two of the boundaries of the sample grid. The situation is analogous for the remaining two boundaries.

Figure 4. Symmetry in reference patterns along the sample grid's boundary.

Gauss-Seidel Relaxation

GSR is a general iterative method of solving a linear system of equations, such as the one formed by the minimization problem discussed in the previous section. For the two-dimensional case the system can be expressed as Ax = b where x is a column vector of nodal variables. At the k-th iteration, GSR updates each nodal variable to be the following [2]:

The initial values of the nodal variables can be obtained, for example, by linearly interpolating known sample points and boundary conditions. In the setting of solving a variational problem, the matrix A essentially encodes the computation of the derivatives used in the energy functional, and the solution of the resulting system corresponds to having the nodal variables minimize the discretized energy functional. Therefore, in this case it is equivalent to use GSR without specifying A and b explicitly, but instead replacing the update step with an iterative minimum finder that uses the local energy formulation (which is an implicit version of A). The latter approach is straightforward to generalize for surfaces in three dimensions.

For the optimization method, I have chosen to use a steepest descent method with an adaptive step size. It is suitable for this problem, because the local energy function to be minimized is a simple polynomial. Having an adaptive step size which is increased or decreased exponentially means that this method efficiently determines the step size, regardless of the initial guess.

The runtime of GSR depends on the desired precision and number of iterations. I have found that in practice it is best to fix the maximum number of iterations. In this case, assuming the precision is also a constant, the runtime is O(n), where n is the number of nodal variables. Note that using GSR with steepest descent directly on the total energy formulation also works, but results in O(n²) runtime, since it takes O(n) time to recalculate it. One way to improve the runtime of GSR is to apply it in a multi-resolution manner, possibly freezing the nodal variables of the previous levels before proceeding [4]. I have implemented this, but didn't perform a detailed comparison, except the one in the Curve and Surface Fitting section, which does demonstrate that multi-level GSR is much faster.

Figure 5. The effect of freezing the boundary on several levels of a surface.

Figures 1 and 2 are examples of running the 2D and 3D versions of the algorithm, respectively. Additionally, Figure 5 shows several levels of the same surface. The shape of the surface is controlled by specifying an initial set of nodal variables some of which GSR is not permitted to relax. I have found that this kind of point "freezing" is more flexible and convenient than adding extra spring terms to the energy functional to "pull" the spline towards those points, as suggested by Szeliski et al. [4]. (I have used this energy method for 2D splines.) In particular, the subdivision of the multi-resolution approach produces extra nodal variables and if some of them need to be interpolated exactly (this might be desirable on the boundary, for example), it is computationally cheaper to freeze them, than to adjust the energy computation. The issue is illustrated in Figure 5; the bottom row contains three levels of the computed surface with freezing the boundary, and the top row - without.

More General Parametric Splines and Surfaces

The conventional formulation of variational problems discussed above focuses on scalar-valued functions which produce an array of scalar nodal variables when discretized on a uniform grid in their domain. In other words, the discrete versions of variational surfaces, for example, are heightmaps and it is impossible to construct a more general parametric surface. However, the values of a function z=f(x,y) may be viewed as (x,y,f(x,y)) = g(x,y). The finite element discretization of g is essentially the same except that it results in vector-valued nodal variables, but energy functionals relying on the derivatives of f are no longer well-defined.

A look at the discrete approximations to the derivatives of f motivates using simple deltas to define the energy of a vector-valued discretization of g. Since the entire problem is solved using the finite element method, there is no reason to formulate explicit energy functionals for the continuous case. The minima of these generalized energy functions will occur at the same locations in space as those of the original discretizations whose nodal variables are scalar.

Figure 6. Several levels of a spherical surface.

The GSR method described in the previous section can be easily updated to use these new energy definitions. With this change, it can produce general parametric splines and surfaces from a set of initial nodal variable values which are arbitrary positions in space. The only cost is some additional computation for the special case of heightmaps: the gradient used in the steepest descent algorithm will have two extra components which are zero. Figure 6 shows the construction of a spherical surface which can not be represented as a heightmap.

Poisson Interpolation

Perez et al. describe an image interpolation technique [3] which is similar to the problem of constructing surfaces using the variational principle. In addition to images, this method has also been successfully applied to heightmapped terrains [8].

Let f be a known bivariate function defined on a domain D - (G - dG) where G is a closed subset of D and dG is its boundary. One way that a function g which is defined on G can be seen as extending the definition of f to all of D is if g satisfies:

where v is a given vector field. Equivalently, g is the unique solution of the following Poisson equation with Dirichlet boundary condition [3]:

Setting v to be the gradient of another known function h results in "pasting" f and h together. This interpolation procedure can be discretized in a similar way to the already discussed variational surface construction and even combined with it into a single problem which is solvable with multi-resolution GSR.

Figure 7. Poisson editing.

Figure 7 presents the results of a Poisson editing operation. The top row contains a "goal" surface whose gradient is used to obtain v and a "target" surface a part of which is edited to match the goal. The bottom row shows (on the left) the result obtained using a separate stage to construct the target and a separate stage to carry out Poisson interpolation. Bottom right contains the result of a single variational construction stage that constructs part of the target as a thin plate surface and part - as a Poisson interpolation of the goal. The difference in the results is due to the boundary dG of the pasting region being able to "relax" more and sink in the non-integrated version.

Curve and Surface Fitting Experiments

Ahlberg et al. present a discussion of the curve fitting properties of natural cubic splines. In particular, they analyze the total deviation of an approximating periodic spline from the cumulative chord length of a circle [1] (recall that natural cubic splines are conventionally defined in terms of functions). I have carried out a similar experiment using a generalized variational periodic curve with beam behaviour. The results of the experiment are not directly comparable to those of Ahlberg et al., because it involves a curve that approximates a circle and a different construction method (multi-level GSR); additionally the finite elements of a variational curve do not correspond to the intervals of definition of a natural cubic. Therefore, this experiment serves only to demonstrate the curve-fitting properties of generalized beam curves and the convergence behaviour of multi-level GSR.

In the table below, "samples" is the number of grid points at which the spline is discretized and "max deviation" is the maximum deviation from the radius in percent; "total deviation" is also in percent. The points of "level 0" are set to lie on the circle. The GSR algorithm is set to the following parameters: energy precision = 0.01 (for 4 initial samples) or 0.0001 (for 8 initial samples), maximum number of iterations = 99, point precision = 0.001, step increase factor = 1.3, step decrease factor = 0.5, initial step size = 5. As a reference for the magnitude of precision-related settings, the radius of the circle is 175; the energy decreases dramatically with the number of sample points.

A similar fitting experiment can be performed for a thin plate surface approximating a sphere (such as in Figure 6). In this case, it is additionally interesting to determine how accuracy is affected by freezing the points of the preceding levels. The table below presents the results (the settings were the same as above; energy precision = 0.0001). Note that the data does not include strips of points along the boundaries of the mesh, since it is difficult to achieve a periodic surface definition for a sphere without a lot of special cases (such as at the poles, where it is impossible to use a simple local calculation for the numeric derivative approximation).

Thin-Plate vs. Bicubic Surfaces

Although the exact relationship between variational surfaces, such as the ones produced with the thin plate model (Figures 5, 6), and bicubic surfaces is complicated, there is one obvious difference: the thin plate functional is a simpler expression than the one from the Minimum Norm Property of bicubic surfaces. This suggests using the MNP expression as a functional to produce pseudo-bicubic surfaces and compare them with the thin plate ones.

Figure 8. Linearly interpolated, pseudo-bicubic, and thin plate surfaces.

It turns out that pseudo-bicubic surfaces are able to have some extra detail, which thin plate surfaces "smooth away". This is demonstrated in Figure 8. On the left is a version of the input mesh from Figure 5, which has been more finely subdivided to produce some ridge-like structures. In the middle is the resulting pseudo-bicubic surface, which has some extra folds around the places where the original mesh is more angular. On the right is the same thin plate surface from Figure 5 for comparison.

Figure 9. "Facets" emerge at higher resolutions.

However, pseudo-bicubic surfaces have a significant drawback: a naive discretization of the fourth derivative is prone to numeric error, resulting in surfaces that look faceted (especially at higher resolution levels). This is hard to eliminate completely. The issue is shown in Figure 9, which compares a pseudo-cubic and a thin plate surface that both approximate a non-periodic version of the torus (at the same resolution level).