Mixed partial derivatives for a surface from surface normals for bicubic interpo

So I’m trying to construct continuous terrain from discretely sampled height values and accompanying surface normal using bicubic interpolation.

Brief description of bicubic interpolation for the maths minded who haven’t encountered it before:
Given any function f(x,y), bicubic interpolation constructs a bicubic polynomial which can be used to interpolate in the range (0,0) - (1,1) from the values of f, f x , fy and fxy (aka f, df/dx, df/dy and d2f/dxdy) at the points (0, 0), (1, 0), (0, 1) and (1, 1). Ie given the values and first derivatives at the corners of the unit square you can interpolate the value (or indeed derivative) at any point inside that square. The nice thing about cubic interpolation is that it produces continuous derivatives.

So lets say I’ve got a surface described by z = f(x, y). Now I’ve related the surface normal to the partial derivatives in x and y by taking the cross product of the vectorized gradients in x and y. Ie Surface Normal, N (not normalised).

[quote]N = [0, 1, fy] x [1, 0, fx] = [-fx, - fy, 1]
[/quote]
But really I’d like to be able to calculate all the partial derivatives (ie fxy as well) from the surface normal and the height and I haven’t been able to figure out a way to do this beyond numerical differentiation which ideally I’d like to avoid.

Anyone have any thoughts?

Cheers