SQUARE =
There are three simple ways to interpolate inside a 2-d rectangular region. But two of them are analogues of each other.
First, how many points determine a plane? THREE. That is, if you have a surface z(x,y), then EXACTLY 3 points are needed to determine a plane. And if you have 4 points, there will in general be no perfectly planar surface that passes through all 4 points.
For example, consider the unit square.
[x,y] = meshgrid([1 2])
z = [1 2;3 1]
surf(x,y,z)
view(2,12)
I've carefully rotated the plot so you can see what was drawn. We see a surface with a ridge between the points where z==2 and z==3. There were two planar triangles created. Again, NO plane can be constructed that passes through all 4 vertices of the square. If we label the 4 corners of the square as
syms A B C D
SQUARE = [A C; B D]
then the two triangles generated were ABC and BCD. Can you see that? The triangles have shared edge BC.
Now the choice made was an arbitrary one, because as easily, the code could have chosen the alternative. Suppose two planes were drawn using a different pair of triangles? ADC and ADB? Now the triangles would have the shared edge AD. But if the above choice were made (and it is PURELY, ABSOLUTELY ARBITRARY which choice was made) then the suface we would see would show a valley, and not a ridge.
Those are two of the interpolation methods we might employ for a quad gridded surface as a function of two independent variables. The third interpolation method is a little slower, sometimes known as bilinear interpolation.
Inside the unit square, it can be thought of as employing basis functions. We can write Z(x,y) as
Z = @(x,y) x.*y*z(1,1) + (1-x).*y*z(2,1) + x.*(1-y)*z(1,2) + (1-x).*(1-y)*z(2,2)
If you look carefully at that functional form, when you substitute some combination of 0 and 1 in for x and y, we get out only the corresponding vertex of that corner of the unit square.
[Z(0,0),Z(0,1);Z(1,0),Z(1,1)]
So z(x,y) as I wrote it does interpolate the unit square. And it does so in a "smooth" way. Let me see if I can get fsurf to help me out here.
fsurf(Z,[0 1 0 1])
Now, any of those three methods can be used to interpolate. In terms of color interpolation, you just interpolate R(x,y), G(x,y), and B(x,y) separately. The triangle style of interpolation is often used in image processing because it is FASTER. A variation of it is probably employed inside your color printer, when you print out color images. Again, speed is crucial, because we are often talking about millions of pixels. But even for a color image like this, you want speed. Do you really want MATLAB to take several times as long to generate a color image? It is only on very coarse things like this where the artifacts become more obvious, so a 2x2 square, or even the 3x4 grid you showed in your example. Those same artifacts are actually buried in the fzurf image too, but there things are fine enough that you don't see any artifacts.
Interestingly, the bilinear interpolation form also has associated artifacts I can point to. (I once wrote a series of papers on the topic, in great depth. Sadly, they were only internal reports, so I cannot give you a link.) Had they used the bilinear style of interpolation then you, or someone else, would surely have found a reason to complain there too. Oh well. We can never make everyone happy. We can only try to make everyone understand what was done, and why. Mathematics imposes the same limits on all of us.
But do you now see where the "sharp" diagonal edge came from? Visualize each rect inside your picture as broken into two triangles. Now think of how the colors will be interpolated, across each of the pair of triangles that dissect the corresponding rectangle. Now go back to the surf image. Again, TWO triangular planes were generated. Consider that surf and pcolor are actually the same underlying codes!
