In the representation of scattered data by smooth pp (:= piecewise polynomial) functions, perhaps the most important problem is to find an optimal triangulation of the given sample sites (called vertices). Of course, the notion of optimality depends on the desirable properties in the approximation or modeling problems. In this paper, we are concerned with optimal approximation order with respect to the given order r of smoothness and degree k of the polynomial pieces of the smooth pp functions. We will only consider C1 pp approximation with r = 1 and k = 4. The main result in this paper is an efficient method for triangulating any finitely many arbitrarily scattered sample sites, such that these sample sites are the only vertices of the triangulation, and that for any discrete data given at these sample sites, there is a C1 piecewise quartic polynomial on this triangulation that interpolates the given data with the fifth order of approximation.