2D cubic B-splines

I have been looking at B-splines to interpolate points. Having 1-D B-splines makes perfect sense to me, but haven't been able to find something that explains 2-D B-splines well for me nor provide me an example.

Suppose I wanted to interpolate $$f(x,y) = \exp(xy)$$ on the unit square using two B-splines. The first with control point at (.25, .25) and the second at (.75, .75). As for data points suppose I can sample data from my function.

How can I go about doing this?

What you need is a tensor product b-spline function, i.e.: $$g(x,y) = \sum_{i=1}^{n_x} \sum_{j=1}^{n_y} a_{ij} B_i(x) B_i(y)$$ where $a_{ij}$ are constants to be determined by your data (for either interpolation or approximation). $n_x$ and $n_y$ are the dimensions of $a_{ij}$ matrix in each direction ($x$ or $y$).
You can solve this problem by either interpolation or approximation (recommended) using least squares: $$\chi^2 \sim \sum_{\mu}(D^{\mu} - g(x^{\mu},y^{\mu}))^2$$ for some data $D$ sampled from your function $f$. Although the tensor product b-spline has a 2-dimensional summation, it is easy to transform to 1-dimension and perform standard linear algebra solutions.
PS Tensor product is, in simplifying the definition for your needs, just a fancy way of saying that you construct a matrix $T_{ij}$ from two vectors $u_i$ and $v_j$ by constructing each matrix element $T_{ij}$ from the product: $$T_{ij} = u_i v_j$$ and respecting the order of the indices $i,j$. Sometimes you see this written as $T_{ij} = u_i \otimes v_j$ or equivalently, $$\mathbf{T} = \mathbf{u} \otimes \mathbf{v}.$$