Polynomial reconstruction on unstructured grids

For a 1D grid I can calculate a Lagrange polynomial through an arbitrary set of points for the reconstruction of a polynomial function.

In 2D I have an unstructured grid and want to interpolate the value on one node from the values on neighboring nodes. A difference from the 1D case is that a well-defined Lagrange polynomial has fixed stencil sizes, like 3 points for a 1st order reconstruction, 6 for a 2nd order reconstruction and so on.

What is the best approach for a reconstruction from the neighboring points? My application is a reconstruction inside a finite volumes scheme with Voronoi cells, which unlike triangles have an arbitrary (but in many meshes roughly constant) number of neighbors.

It is not polynomial, but you may be interested in the Natural Neighbors interpolation (it fits well with Voronoi diagram). To evaluate the interpolant at a given point, insert the point into the diagram, and compute the volumes of the intersections between the new Vornoi cell (of the added point) and the Voronoi cells in the diagram before inserting the points.This gives the interpolation weights. The interpolant is $C^2$ almost everywhere (on the points it is $C^1$).