  Given cell-average data values of a piecewise smooth bivariate function $f$
within a domain $\Omega$, we look for a piecewise adaptive approximation to
$f$. We are interested in an explicit and global (smooth) approach. Bivariate
approximation techniques, as trigonometric or splines approximations, achieve
reduced approximation orders near the boundary of the domain and near curves of
jump singularities of the function or its derivatives. Whereas the boundary of
$\Omega$ is assumed to be known, the subdivision of $\Omega$ to subdomains on
which $f$ is smooth is unknown. The first challenge of the proposed
approximation algorithm would be to find a good approximation to the curves
separating the smooth subdomains of $f$. In the second stage, we simultaneously
look for approximations to the different smooth segments of $f$, where on each
segment we approximate the function by a linear combination of basis functions
$\{p_i\}_{i=1}^M$, considering the corresponding cell-averages. A discrete
Laplacian operator applied to the given cell-average data intensifies the
structure of the singularity of the data across the curves separating the
smooth subdomains of $f$. We refer to these derived values as the signature of
the data, and we use it for both approximating the singularity curves
separating the different smooth regions of $f$. The main contributions here are
improved convergence rates to both the approximation of the singularity curves
and the approximation of $f$, an explicit and global formula, and, in
particular, the derivation of a piecewise smooth high order approximation to
the function.
