I'm trying to come up with a simple and efficient way to create a smooth surface which intersects a number of given "sample" points.
For any X,Y point on the surface, I identify up to 4 sample points in each of the 4 directions (the next higher and lower points on the X, and then the Y axes). Given this point, I want a way to compute a Z value that interpolates between the 4 sample points.
Of course the algorithm, given the X, Y position of any of the 4 sample points, should output the Z value for that point. Note also that there may be less than 4 sample points.
I guess some function of the Z values for the 4 sample points, somehow inversely biased by the distance to the sample point, but I can't figure out how to do this.
Anyone got any ideas about a simple way to do this?
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.