You can do this by constructing patches from Catmull-Rom splines. These splines will hit each of the control points and they are continuous in the first derivative (though not the second). I also find them to be extremely easy to work with. The math is straightforward and they behave intuitively with slight changes in the control points.
At the highest level, you'll need 16 points per patch (at the edge of your dataset, you can use the corner and edge points twice in the same spline).
First, you'll need to interpolate across the points p[i][j] in each row in your 4x4 matrix to create a set of four intermediate control points q[i]. Here's a rough ASCII sketch of what I mean.
p00 p01 q0 p02 p03
p10 p11 q1 p12 p13
p20 p21 q2 p22 p23
p30 p31 q3 p32 p33
Now you can interpolate between each of those four intermediate control points to find a final splined point on your surface.
Here is a construction of the Catmull-Rom spline across four points. In this example, you are interpolating between points p[i-1] and p[i], using control points on either side p[i-2] and p[i+1]. u is the interpolation factor ranging from zero to one. τ is defined as the tension on the spline and will affect how tightly your splined surface conforms to your control points.
| 0 1 0 0 | | p[i?2] |
|?τ 0 τ 0 | | p[i?1] |
p(u) = 1 u u2 u3 | 2τ τ?3 3?2τ ?τ | | p[i] |
|?τ 2?τ τ?2 τ | | p[i+1] |
NOTE: it's not immediately obvious how to lay this out in Stackoverflow's gui but u2 and u3 are supposed to represent u squared and u cubed, respectively.
与恶龙缠斗过久,自身亦成为恶龙;凝视深渊过久,深渊将回以凝视…