The Delaunay triangulation, in both classic and more generalized sense, is studied in this paper for minimizing the linear interpolation error (measure in $L^p$-norm) for a given function. The classic Delaunay triangulation can then be characterized as an optimal triangulation that minimizes the interpolation error for the isotropic function $||x||^2$ among all the triangulations with a given set of vertices. For a more general function, a function-dependent Delaunay triangulation is then defined to be an optimal triangulation that minimizes the interpolation error for this function and its construction can be obtained by a simple lifting and projection procedure.
The optimal Delaunay triangulation is the one that minimizes the interpolation error among all triangulations with the same number of vertices, i.e. the distribution of vertices is optimized in order to minimize the interpolation error. Such a function-dependent optimal Delaunay triangulation is proved to exist for any given convex continuous function. On an optimal Delaunay triangulation associated with $f$, it is proved that $\nabla f$ at the interior vertices can be exactly recovered by the function values on its neighboring vertices. Since the optimal Delaunay triangulation is difficult to obtain in practice, the concept of nearly optimal triangulation is introduced and two sufficient conditions are presented for a triangulation to be nearly optimal.