Display omitted We propose a novel method to approximate a function on 2D domain by piecewise polynomials. The Voronoi tessellation is used as a partition of the domain, on which the best fitting polynomials in L2 metric are constructed. Our method optimizes the domain partition and the fitting polynomials simultaneously by minimizing an objective function indicating the approximation quality. We also provide the explicit formula of the gradient of the objective function, which makes an efficient gradient-based algorithm workable for the function minimization. We conduct several experiments to demonstrate the efficacy of our new approach for generating piecewise polynomial approximations of analytic functions and color images.