Variational Subdivision
for Laplacian Splines

Technical Report
TR97-291

Abstract

The fundamental problem of geometric design is the representation of curved shapes. Traditionally such shapes are represented by parametric spline curves, e.g. NURBS, which are defined as the minimizers of variational problems. For example, cubic B-splines minimize the bending energy functional.

More recently subdivision curves and surfaces emerged as an alternative means for representing curved shapes. In this framework a curve or surface is defined as the limit of a repeated averaging process of control points. The subdivision scheme is determined by the subdivision mask S which specifies the weights used during the averaging.

A methodology for deriving subdivision schemes from homogeneous variational problems has been outlined in [Warren97]. This report shows the computational details of this approach and presents a method for deriving a subdivision scheme with predefined local support wich produces limit surfaces close to the minimizer of the variational problem. In particular this document contains the detailed derivation of Laplacian Splines discussed in section 4.2 of [Warren97].

The result will be a local stationary subdivision scheme for bounded rectangular grids which produces a limit surface that is close to the minimizer of the original variational problem. By increasing the support of the subdivision basis functions the resulting surface will come arbitrarily close to the real variationally optimal surface.

The Mathematica source code for this report can be found under URL http://www.cs.rice.edu/~jwarren/Laplace-tr.nb and an on-line HTML version under http://www.cs.rice.edu/~jwarren/Laplace-tr .

Introduction

Laplacian Splines are defined as the minimizer to the variational problem derived from Laplace's equation

where the domain Omega is the real plane.The strong form of this equation is exactly Laplace's equation

We want to refer to functions that satisfy this equation as Laplacian Splines.

This report presents a subdivision scheme which produces the minimizer of the variational probelm as stated in (1) for a bounded rectangular grid. This subdivision scheme has global support. We will also show how a local subdivision scheme can be constructed which produces a surface that is close to the minimizer of (1). By increasing the support of the local basis function of this approximating subdivision scheme, the resulting limit surface can lie arbitrarily close to the real Laplacian spline.

Definition of the Energy Matrix

The energy matrix for Laplacian splines is derived from Laplace's equation. We use a generating function representation of the energy mask e(z).

GenRow[i,j,n] now generates the coefficients of the energy matrix for control point (i,j) in an (n+1)x(n+1) grid. You may want to re-evaluate GenRow below for different values of i and j. GenRow2[i,j,n] returns the coefficients of the energy matrix for the subdivided grid (only for control points corresponding to control points on the previous level grid, i.e. before subdivision).

En[n] computes the energy matrix for a (n+1) x (n+1) grid, i.e. a (n+1)^2 x (n+1)^2 matrix computing the energy associated with each of the control points.

UEn[n] computes the energy matrix after upsampling, i.e. a matrix which computes the energy associated with the control points in a (n/2+1) x (n/2+1) grid and then upsamples them to the corresponding (n+1) x (n+1) grid by inserting zeroes as th energy associated with new control points.

Finding The Exact Subdivision Mask

En[n] as defined above is singular, i.e. not invertible. Thus we can not solve En S = UEn by inverting En. To make En invertible we replace one of the (linearly dependent) rows in En by a constraint representing the preservation of moments. Call these modified matricec NewUn and NewUEn respectivey.

Evaluate the following two lines to compute the exact subdivision matrix carrying control points from a 5x5 to a 9x9 control grid and from a 9x9 to a 17x17 grid. The resulting sequence of subdivision matrices will produce a limit surface that exactly solves Laplace's equation.

Unfortunately the resulting subdivision rules have global support (the exact solution must have global support as the variational problem is global). Here are some of the resulting exact subdivision rules rounded to three digits. Note the fast decay of the coefficients.

An internal vertex rule:

An internal horizontal edge rule:

An internal face rule:

Apply the subdivision matrices to some control points

Note that the limit surface will be the exact solution to Laplace's equation.

Finding a Sparse Approximation to the Subdivision Mask

The exact solution to the equation En S = U En as computed above leads to subdivision rules with global (but exponentially decaying) support.

In the sequel we want to find a subdivision scheme with local support that approximates the solution well. Our approach is to set up an unknown subdivision matrix with variable entries at locations that enforce the desired symmetries and sparsity. Call this subdivision matrix Shat. Now we can solve for the actual coefficients in Shat by minimizing the expression
En Shat = U En in the infinity norm. Recall that the infinity norm of a matrix is the maximum of the sums of the absolute values of entries in any row.

Set up a sparse subdivision matrix with the desired symmetry and sparsity

Some helper functions

Definition of the Sparse Matrix

Now Shat[n] computes a sparse (variable) subdivision matrix carrying the control points from a bounded nxn grid to a (2n-1)x(2n-1) grid.

For example Shat[3] subdivides a 3x3 grid into a 5x5 grid. Note that we enforced sparsity and symmetry in the structure of this matrix.

Set up a Linear Programming Problem to Solve for the Actual Coefficients in the Subdivision Matrix

Starting from the equation E S = U E we set up a linear programming problem to solve for S. The linear program will find an S such that the norm || E S - U E || is minimal, where ||.|| denotes the L_infinity norm, i.e. the sum of the absolut values of the elements in any one row should be minimal.

One crucial observation is that the symmetry and sparsity structure of S as defined above laeds only to a finite and fixed number of different rows in E S - U E, i.e. it is sufficient to minimize this small and fixed number of rows once to find a minimal solution for any size of S.

We want to solve here separately for the purely interior subdivisin rules as solving the whole system is beyond the abilities of Mathematica. We solved for the complete system using the CPLEX linear programming package.

Some helper functions to set up the linear program

Get the distinct rows in E S - U E

Evaluate the following line to find all the different rows in E S - U D.

Note that using a larger system leads to the same distinct rows in E S - U D:

Partition these rows of E S - U E into three categories:
purely interior
interior and boundary
interior, boundary and corner

Set up the systems of constraints

Set up constraints to enforce constant and linear precision

Finally solve for the coefficients for the interior

Remember that all variables cann be positive or negative but the solver handles only positive variables. Thus all variables X in the system have to replaced by terms of the form (Xpositive - Xnegative).

Here are the resulting subdivision rules for the interior:

Sparse Solution for the Complete Grid

Mathematica's linear optimizer "ConstrainedMin" simply seems to run forever on the linear program arising from the complete set of constraints. We used an external solver, CPLEX, to solve for these coefficients instead.

Note that by construction these subdivision rules have constant precision (enforced in the constraints).

Here are the results:

Now apply the sparse S to some control points and compare this visually to the exact solution

The sparse approximation based on 5x5 basis functions:

The exact subdivision process with globally supported basis functions:

References

[Warren97] Joe Warren: Subdivision Schemes for Variational Splines. July 1997, submitted for publication to CAGD.

Acknowledgements

This work has been supported in part by NSF grant CCR-9500572 and TATP grant 003604-010.