Suppose we have m+1 rows and n+1 columns of data points D_{ij} (0 <= i <= m and 0 <= j <= n) and wish to find a Bspline 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 Bspline 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 s_{c} (0 <= c <= m) in the udirection and t_{d} (0 <= d <= n) in the vdirection by setting e and f to m and n, respectively. We also obtain knot vectors U and V for the u and vdirection, respectively. Therefore, what remains is to find the desired control points!
Global Surface Interpolation 
Suppose the Bspline surface is given as follows:
Since it contains all data points and since parameters s_{c} and t_{d} correspond to data point D_{cd}, plugging u=s_{c} and v=t_{d} into the surface equation yields:
Since N_{i,p}(s_{c}) 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 P_{ij}. Therefore, we can define the inner expression to be a new term as follows:
More precisely, if i is fixed to the same value, Q_{id} is the point, evaluated at t_{d}, on the Bspline curve of degree q defined by n+1 "unknown" control points on row i of the P's (i.e., P_{i0}, P_{i1}, ..., P_{in}). Plugging Q_{id} into the equation of D_{cd} above yields
Thus, data point D_{cd} is the point, evaluated at s_{c}, of a Bspline curve of degree p defined by m+1 "unknown" control points on column d of the Q's (i.e., Q_{0d}, Q_{1d}, ..., Q_{md}). Repeating this for every c (0 <= c <= m), the dth column of data points (i.e., D_{0d}, D_{1d}, ..., D_{md}) is obtained from the dth column of Q's (i.e., Q_{0d}, Q_{1d}, ..., Q_{md}) and parameters s_{0}, s_{1}, ..., s_{m}. Since data points D_{0d}, D_{1d}, ..., D_{md}, degree p and parameters s_{0}, s_{1}, ..., s_{m} are known, this equation means
Find the dth column of control points Q_{0d}, Q_{1d}, ..., Q_{md}, given degree p, parameters s_{0}, s_{1}, ..., s_{m} and the dth column of data points D_{0d}, D_{1d}, ..., D_{md}. 
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 Q_{cd}'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 Q_{id}, which is duplicated below. In this equation, data points on row i of the Q's (i.e., Q_{i0}, Q_{i1}, ..., Q_{in}) are the points on a Bspline curve, evaluated at t_{0}, t_{1}, ..., t_{n}, of degree q defined by n+1 unknown control points P_{i0}, P_{i1}, ..., P_{in}. Therefore, applying curve interpolation with degree q and parameters t_{0}, t_{1}, ..., t_{n} 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 Bspline surface of degree (p,q) of the given data points. Therefore, surface interpolation using Bspline consists of (m+1) + (n+1) curve interpolations! The following summarizes the required steps:
Input: (m+1)×(n+1) data points
D_{ij} and degree (p,q);
Output: A Bspline surface of degree (p,q) that contains all data points; Algorithm:
Compute parameters in the vdirection t_{0}, t_{1}, ..., t_{n} and the knot vector V; for d := 0 to n do /* for column d */

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 LUdecomposition 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 LUdecomposition should also be computed before the second for starts. Without doing so, we will perform (m+1)+(n+1)=m+n+2 LUdecompositions when only 2 would be sufficient.
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 Bspline 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  
Surfaces 