Given a set S̄ of n points in R3, sampled from an unknown bivariate function f(x, y) (i.e., for each point p ∈ S̄, zp = f(xp, yp)), a piecewise-linear function g(x, y) is called an ε-approximation of f(x, y) if for every p ∈ S̄, |f(x,y)-g(x,y)|≤ε. The problem of computing an ε-approximation with the minimum number of vertices is NP-Hard. We present a randomized algorithm that computes an ε-approximation of size O(c2log2 c) in O(n2+δ+c3log2clog n/c) expected time, where c is the size of the ε-approximation with the minimum number of vertices and δ is any arbitrarily small positive number. Under some reasonable assumptions, the size of the output is close to O(clog c) and the expected running time is O(n2+δ). We have implemented a variant of this algorithm and include some empirical results.