I am currently working on 2D interpolation and I am using cubic splines. In 1D given N points one would use N+2 splines and would have to apply additional constraints in order to achieve a square system. For example requiring that the interpolant, S, second derivative is equal to zero ,Sxx = 0, at the end points (Natural Spline conditions) I was wondering how one would extend this to 2D? I believe it should be something like Sxx = Syy = 0 at the edge of your domain but I am not sure (if so I believe we would still be missing some equations). Any help would be appreciated.I will provide an example in 2D: Consider the points 1:1:10 in the x component and 1:1:10 in the y component, thus N points in each direction. For cubic splines one need N+2 splines in each direction.  In 1D interpolation we can add the constraint that the second derivative should be zero to arrive at (N+2) equations for (N+2) unknowns.Thus to extend to 2D I will take the tensor product of these points I will have a 10x10 grid (100 pairs of points) and the tensor product of 12 splines in each direction (144 splines). Since I am using these for interpolation then I can retrieve 100 equations by setting up an equation for each point in my grid. I can retrieve an additional 40 equations by asking that S_{xx} = S_{yy} = 0 at the boundaries. 10 at $(x_0,y)$, $(x_N,y)$, $(x,y_0)$, $(x,y_N)$ Thus I would be missing 4 equations. This is where my confusion lies. 
Natural spline conditions for 2D cubic splines
image processing