`decimal=2` rang alarm bells...

I'm not a fan of comparing the bilinear interpolation to the non-linear reference surface. It's comparing the results to the wrong thing. We don't even want the results to be equal to the reference surface, we want them to be equal to the bilinear interpolation.

So I'd suggest:
 - Use a (piecewise) linear reference surface when doing these bulk tests so there's no need for the sloppy numerical tolerance.
 - Use a non-linear reference surface for one or two non-trivial interpolation points (e.g. using target_x[1], target_y[0]) and compare (with a tight tolerance) with an explicitly calculated bilinear interpolation value. In other words, include the bilinear calculation of the expected value in the testing code.