This work was motivated by the problem of obtaining a smooth density function over a geographical region from data aggregated over irregular subregions. Minimization of a family of roughness criteria given 'volume' data lead to smooth multivariate functions - Laplacian histosplines, having a certain order of the iterated Laplacian of constant value in each of the subregions and satisfying natural boundary conditions on the boundary of the region. For inexact data, e.g., in case of estimating an underlying density given counts of events by subregions, Laplacian smoothing histosplines are constructed, analogous to smoothing splines in the univariate case, and a method for choosing the smoothing parameter is presented. For both cases of exact and inexact data, modified roughness criteria, independent of the region, are discussed, and results known for point-evaluation data are extended to the case of aggregated data. (Author)