The simple solution which I went to is to do the linear interpolation, with the following
approach. First we look at the heightmap "from the top"  just looking at its X and Y coordinates. And logically we split the XY space into a bunch of 1sized squares. So we always
have our x,y point for which we are calculating, inside one of these small tiles. So we can zoom into this tile and take further look:
(x0,y0) (x0+1, y0)
 /
v v
++> X
\ 1 /
 \ / 
4 + 2
 / \ 
/ 3 \
++ < (x0+1, y0+1)

v
Y
in this picture  the corners of the square are the integer X, Y values which have the data.
The point in the center is the "virtual" point which is used to split the square into four triangles  it's x,y coordinates are simply X0+0.5, Y0+0.5.
Then we can pretty simply find which of the triangles holds our "interesting" point  now it's all dead simple  find the equation of the plane that passes through the three points the coordinates of which we know (the Z coordinate of the "center" point is the average of the Z coordinates of all 4 corners)  and then simply substitute the X, Y into the equation, to get Z straight away.
I don't have walking animations yet  for now I simply set the Z coordinate of the avatar, which gets its butt approximately at ground level :) While running around, if one would measure by the proportions to the real world, the staircase effect is reduced to almost unnoticeable ~10cm if we relate to the realworld body proportions. And I hope that once I get the avatar's butt out of the ground, it should be even better visually.
Now, the obvious question  why 4 triangles, instead of 2 ? Because if you take the following layout:
0*
 
*0
and the "*" is "high" an "0" is "low"  depending on the choice to cut this XY box, you get either a pit, or a spike. 4 triangles handle this case noticeably better.