Since a Bspline surface of degree (p,q) defined by e+1 rows and f+1 columns of control points has the following equation
it requires two sets of parameters for surface interpolation and approximation. Suppose we have m+1 rows and n+1 columns of data points D_{ij}, where 0 <= i <= m and 0 <= j <= n. Thus, m+1 parameters s_{0}, ..., s_{m} in the udirection (i.e., one for each row of data points) and n+1 parameters t_{0}, ..., t_{n} in the vdirection (i.e., one for each column of data points) are required so that point (s_{c},t_{d}) in the domain of the surface corresponds to point S(s_{c},t_{d}) on the surface, which, in turn, corresponds to the data point D_{cd}. Putting this in an equation form, the point on the surface, evaluated at (s_{c},t_{d}) as shown below,
corresponds to the data point D_{cd}, where s_{c} and t_{d} are parameters in the u and vdirection, respectively.
In the equation of the desired Bspline surface, the udirection corresponds to the index i in N_{i,p}(u) and P_{ij}. Since i runs from 0 to m, N_{0,p}(u), N_{1,p}(u), ..., N_{m,p}(u) are coefficients of the columns of the control points. Therefore, in the udirection, we need m+1 parameters, and using a parameter computation method that takes degree p and the data points on column j, we can compute m+1 parameters u_{0,j}, u_{1,j}, ..., u_{m,j}. This is shown in the figure below. The desired parameters s_{0}, s_{1}, ..., s_{m} are simply the average of each row. More precisely, parameter s_{i} is the average of parameters on row i, s_{i} = (u_{i,0} + u_{i,1} + ... + u_{i,n})/(n+1).
The computation of parameters for the vdirection is similar. Each row of data points has n+1 points and hence requires n+1 parameters. Thus, for data points on row i, we can compute n+1 parameter values v_{i,0}, v_{i,1}, ..., v_{i,n}. Since we have m+1 rows, these values can be organized into a (m+1)×(n+1) matrix, and parameter t_{j} is simply the average of the parameters on column j, t_{j} = (v_{0,j} + v_{1,j} + ... + v_{m,j})/(m+1). In this way, we obtain n+1 parameters for the vdirection.
This procedure is summarized as follows:
Input: (m+1)×(n+1) data points
D_{i,j};
Output: Parameters in the udirection s_{0}, ..., s_{m} and parameters in the vdirection t_{0}, ..., t_{n}; Algorithm:
for j := 0 to n do /* column j */
/* calculate t_{0}, ..., t_{n} */ for i := 0 to m do /* row i */

Then, parameter values s_{0}, s_{1}, ..., s_{m} and degree p jointly define a knot vector U in the udirection. Similarly, parameters t_{0}, t_{1}, ..., t_{n} and degree q jointly define a knot vector V in the vdirection. Click here to review the way of computing a knot vector from a set of parameters.
Note that the above is only a conceptual algorithm and is inefficient. Once you know what is going on, it is quite easy to convert it to a very efficient algorithm. Note also that this algorithm only works for the uniformly spaced, chord length and centripetal methods. For the universal method, because there is no data points involved, we can apply uniform knots to one row and one column of data points for computing the parameters.