Obtaining Interpolated Values Of Z Based On X And Y

I'm working on a program, which gives me interpolated values of z from a digitized graph based on the cursor position or input values of x and y. Thing that bothers me for quite some time is how to get properly interpolated z value. I was using SADA (spatial analysis software), where I was creating a 100x100 matrix, but the interpolation methods there are not as accurate as I need.

Basically, from digitizing a graph I get a matrix:

VB:
x1, y1, z1;
x2, y2, z1;
x3, y3, z2;
x4, y4, z2;
.
.
.
x3742, y3742, z24

Now, inputting x and y that are not in the matrix should give me interpolated value of z, which gets even more tricky when the axes are in log scale. The only idea that I have in my mind and that I think would work is this:

for each missing z, from zmin=0,625 to zmax=3 , every 0,001, I would make 1000 (enough accuracy in this case) xs and ys. How to automate the process.