This question is a search for further answers from a question on maths.stackexchange.com.
I've inherited some numerical quadrature code that is designed to integrate sparse 2D data. The quadrature is not being done in any of what I would consider the "normal" ways, but seems to be using an ad-hoc method based on Voronoi cell areas. However the data has some unusual properties that may make this sensible.
The samples fall roughly on a grid such that there is either one or zero sample points inside each grid cell. (But randomly positioned inside the cell, not at the cell center).
The integration then proceeds by adding "ghost" points to the exterior of the region (to bound the Voronoi regions at the non-ghost points), generating a Delaunay triangulation of the domain and then assigning a weight for each sample point equal to the area of the corresponding Voronoi region.
I can see that this scheme tries to minimises the error for C1 functions.
However having gone to the trouble of triangulating the domain, I would have expected the integration to be based upon the neighbouring linear elements, which should give exact answers for linear functions, and thus I'd expect an error estimate (at least for C2 functions) of one order of $h$ higher.
However some numerical experiments I performed suggested that for smooth data both methods had similar convergence rates. (Details can be found in my answer to the original math.SE question)
It was noted in one of the answer to a original question that one possible reason for preferring the Voronoi based scheme is continuity of the weights with changes in the sample locations, (which fails for the Delaunay triangular weights as the mesh can change discontinuously with sample locations)
So my questions are:
Are there known advantages to using the Voronoi method over the Delaunay style method, or vice-versa? Can these be quantified?
Why am I seeing roughly the same convergence rates from the two methods?