We present an efficient algorithm to obtain a triangulated graph surface for scattered points (x[?] y[?])[?], i= 1...n, with assoiciated function values fnof;[?]. The maximal distance between the triangulated graph surface and the function values is measured in z-direction (z= fnof; (x, y) and lies within a user-defined tolerance. The number of triangles is minimized locally by adapting their shape to different second-degree least squares approximations of the underlying data. The method consists of three major steps: (i) subdividing the given discrete data set into clusters such that each cluster can be approximated by a quadratic polynomial within a prescribed tolerance; (ii) optimally triangulating the graph surface of each quadratic polynomial; and (iii) ''stitching'' the triangulations of all graph surfaces together. We also discuss modifications of the algorithm that are necessary to deal with discrete data points, without connectivity information, originating from a general two-manifold surface, i.e., a surface in three-dimensional space that is not necessarily a graph surface.