Surface Global Interpolation

Suppose we have m+1 rows and n+1 columns of data points Dij (0 <= i <= m and 0 <= j <= n) and wish to find a B-spline surface of degree (p,q) that contains all of them. Similar to the curve case, we have data points and degrees p and q as input. To define an interpolating a B-spline surface, we need two knot vectors U and V, one for each direction, and a set of control points. Like the curve case, the number of control points and the number of data points are equal (i.e., (m+1)×(n+1) control points in m+1 rows and n+1 columns). With the technique discussed in Parameters and Knot Vectors for Surfaces, we can compute two sets of parameters sc (0 <= c <= m) in the u-direction and td (0 <= d <= n) in the v-direction by setting e and f to m and n, respectively. We also obtain knot vectors U and V for the u- and v-direction, respectively. Therefore, what remains is to find the desired control points!

Global Surface Interpolation

Given a grid of (m+1)×(n+1) data points Dij (0 <= i <= m and 0 <= j <= n) and a degree (p, q), find a B-spline surface of degree (p,q) defined by (m+1)×(n+1) control points that passes all data points in the given order.

Finding a Solution

Suppose the B-spline surface is given as follows:

Since it contains all data points and since parameters sc and td correspond to data point Dcd, plugging u=sc and v=td into the surface equation yields:

Since Ni,p(sc) is independent of the index j, it can be moved out of the summation of j:

If we examine the inner term, we shall see that the index i only appears in Pij. Therefore, we can define the inner expression to be a new term as follows:

More precisely, if i is fixed to the same value, Qid is the point, evaluated at td, on the B-spline curve of degree q defined by n+1 "unknown" control points on row i of the P's (i.e., Pi0, Pi1, ..., Pin). Plugging Qid into the equation of Dcd above yields

Thus, data point Dcd is the point, evaluated at sc, of a B-spline curve of degree p defined by m+1 "unknown" control points on column d of the Q's (i.e., Q0d, Q1d, ..., Qmd). Repeating this for every c (0 <= c <= m), the d-th column of data points (i.e., D0d, D1d, ..., Dmd) is obtained from the d-th column of Q's (i.e., Q0d, Q1d, ..., Qmd) and parameters s0, s1, ..., sm. Since data points D0d, D1d, ..., Dmd, degree p and parameters s0, s1, ..., sm are known, this equation means

Find the d-th column of control points Q0d, Q1d, ..., Qmd, given degree p, parameters s0, s1, ..., sm and the d-th column of data points D0d, D1d, ..., Dmd.

So, it is a curve interpolation problem! More precisely, the curve global interpolation can be applied to each column of the data points to obtain a column of control points Qcd's. Since we have n+1 columns of data points, we shall obtain n+1 columns of Q's.

Now, we can use the same idea but apply to the equation of Qid, which is duplicated below. In this equation, data points on row i of the Q's (i.e., Qi0, Qi1, ..., Qin) are the points on a B-spline curve, evaluated at t0, t1, ..., tn, of degree q defined by n+1 unknown control points Pi0, Pi1, ..., Pin. Therefore, applying curve interpolation with degree q and parameters t0, t1, ..., tn to row i of the Q's gives row i of the desired control points!

Once all rows of control points are found, these control points, along with the two knot vectors and degrees p and q define an interpolating B-spline surface of degree (p,q) of the given data points. Therefore, surface interpolation using B-spline consists of (m+1) + (n+1) curve interpolations! The following summarizes the required steps:

Input: (m+1)×(n+1) data points Dij and degree (p,q);
Output: A B-spline surface of degree (p,q) that contains all data points;
    Compute parameters in the u-direction s0, s1, ..., sm and the knot vector U;
    Compute parameters in the v-direction t0, t1, ..., tn and the knot vector V;
    for d := 0 to n do /* for column d */
      begin /* computing "intermediate data points" Q's */
        Apply curve interpolation to column d of data points (i.e., D0d, D1d, ... Dmd) using
        • degree p
        • parameters s0, s1, ..., sm
        • knot vector U
        The result is column d of the "intermediate data points" Q0d, Q1d, ..., Qmd
    for c := 0 to m do /* for row c */
      begin /* computing the desired control points P's */
        Apply curve interpolation to row c of the Q's (i.e., Qc0, Qc1, ... Qcn) using
        • degree q
        • parameters t0, t1, ..., tn
        • knot vector V
        The result is row c of the desired control points Pc0, Pc1, ..., Pcn
    The computed (m+1)×(n+1) control points Pij, degree (p,q), and knot vectors U and V define an interpolating B-spline surface of degree (p,q) for the given data points.

Note that matrix N used in the first for loop does not change when going from column to column. Therefore, if we naively follow the above algorithm, we would end up solving the system of equation D = N.Q n+1 times. This is, of course, not efficient. To speed up, the LU-decomposition of N should be computed before the first for loop starts, and the interpolation in each iteration is simply a forward substitution followed by a backward substitution. Similarly, a new matrix N will be used for the second for loop. Its LU-decomposition should also be computed before the second for starts. Without doing so, we will perform (m+1)+(n+1)=m+n+2 LU-decompositions when only 2 would be sufficient.

Why Is This Method Global?

Because the interpolating surface is obtained through m+n+2 global curve interpolations, it is obvious that this interpolation technique is global. The following images show the impact of moving a data point. This B-spline surface is obtained by using six rows (m=5) and five columns (n=4) data points with degree (3,3). Images on the top row show the knot curves (i.e., isoparametric curves that correspond to knots), and the data point marked by a yellow circle indicates the point to be moved. The bottom row shows the corresponding rendered surfaces. It is obvious that after changing the position of a data point, the blue knot curve changes its shape drastically, and moves closer to its neighboring data points. The impact propagates to the right as well because the valley in the far right end is a little deeper as can be seen from the red knot curves. The more dramatic change is shown in the rendered surfaces. Although the indicated data point is moved to the right, this move creates a big bulge in the left end. Please also note that the front boundary curve, actually a red knot curve, also changes its shape slightly. Therefore, the impact of changing the position of a single data point can propagate to entire surface, and, hence, this is a global method.

Before Move After Move
Knot Curves