Given a set (S) over bar of n points in R(3), sampled from an unknown bivariate function f(x, y) (i.e., for each point p is an element of (S) over bar, z(p) = f(x(p),y(p))), a piecewise-linear function g(x, y) is called an epsilon-approximation of f(x, y) if for every p is an element of (S) over bar, \f(x,y) - g(x,y)\ less than or equal to epsilon. The problem of computing an e-approximation with the minimum number of vertices is NP-Hard. We present a randomized algorithm that computes an epsilon-approximation of size O(c(2) log(2) c) in O(n(2+delta) + c(3) log(2) clog n/c) expected time, where c is the size of the a-approximation with the minimum number of vertices and delta is any arbitrarily small positive number. Under some reasonable assumptions, the size of the output is close to O(c log c) and the expected running time is O(n(2+delta)). We have implemented a variant of this algorithm and include some empirical results.