Motivated by applications in computer graphics, visualization, and scientific computation, we study the computational complexity of the following problem: given a set S of n points sampled from a bivariate function f(x, y) and an input parameter epsilon > 0, compute a piecewise-linear function Sigma(x, y) of minimum complexity (that is, an xy-monotone polyhedral surface, with a minimum number of vertices, edges, or faces) such that \Sigma(x(p), y(p)) - z(p)\ less than or equal to epsilon for all (x(p), y(p), z(p)) is an element of S. We give hardness evidence for this problem, by showing that a closely related problem is NP-hard. The main result of our paper is a polynomial-time approximation algorithm that computes a piecewise-linear surface of size O(K-o log K-o), where K-o is the complexity of an optimal surface satisfying the constraints of the problem. The technique developed in our paper is more general and applies to several other problems that deal with partitioning of points (or other objects) subject to certain geometric constraints. For instance, we get the same approximation bound for the following problem arising in machine learning: given n "red" and m "blue" points in the plane, find a minimum number of pairwise disjoint triangles such that each blue point is covered by some triangle and no red point lies in any of the triangles.