Natural neighbor interpolation Recently I am interested in Natural neighbor interpolation, that is :
Given a function $P(x)$ and some interpolation points $\{x_i,P(x_i)\}_{i=1}^N$, we have the interpolation function $$P^*(x)=\sum_{i=1}^N{\omega_i(x)P(x_i)},$$
where $\omega_i(x)=\dfrac{area\{\Omega_x\cap\Omega_{x_i}\}}{area\{\Omega_x\}},\Omega_x=\{y:\|y-x\|_2\leq\|y-x_k\|_2,k=1,\ldots,N\}.$
Now assume that $$L=\max_{j,k}\left\{\dfrac{|P(x_j)-P(x_k)|}{\|x_j-x_k\|}\right\},$$ my question is whether there exists a constant $C$ such that $$\sup_{X,Y}\left\{\dfrac{|P^*(X)-P^*(Y)|}{\|X-Y\|}\right\}\leq CL$$
holds, which means that the Lipschitz constant of the whole domain is bounded by that of the interpolation points? 
Any counterexample or proof is welcome. Thank you.
 A: This estimate can be established under some restrictions how the interpolation points are arranged. But for general arrangement, I doesn't appear that that the estimate can hold. Below is a counterexample, followed by discussion of the literature which ensures the bound in some situations.
First is the counterexample. Let $P(x)| = |x|$. Consider the following 10 interpolation points for a small number $h>0$ and a much smaller number $\epsilon > 0$, 
$$
(-\sqrt{h},0), (-\epsilon,0), (\epsilon, 0), (\sqrt{h},0),\\
(-1,1), (0,1),(1,1),(-1,-1),(0,-1),(1,-1).
$$
The last six points just provide some bounds for the domain so that Voronoi cells in question don't extend to infinity. All the interesting behavior is related to the first four points. 
Now we will focus on interpolating at two points $(0,0)$ and $(0,h)$. Throughout, we emphasize the highest order terms as that is all that is relevant as $h$ and $\epsilon$ approach $0$, but a truly rigorous proof should be a little more careful. Starting with $(0,0)$, the overlap of the Voronoi cells looks like this:

Above, the Voronoi diagram of the known data points is shown with dashed lines while the Voronoi diagram of the augmented point set including the interpolation location $(0,0)$ is the solid lines. The new Voronoi cell of $(0,0)$ mostly overlaps with original Voronoi cells of $(-\epsilon, 0)$ and $(\epsilon,0)$ and slightly overlaps with those of $(0,1)$ and $(0,-1)$. The area of overlap with the latter two is very small ($O(\epsilon^3)$, thus
$$
P^*(0,0) \approx \frac{P(-\epsilon,0) + P(\epsilon,0)}{2} = \epsilon.
$$
Now, we can consider the interpolant at (0,h). The Voronoi diagrams look like this:

Keep in mind that this picture is exaggerated: we are considering the behavior as $h,\epsilon\rightarrow 0$ with $\epsilon \ll h$. 
First consider the overlap with the Voronoi cell of ($\sqrt{h},0)$: the area contains a triangle of height $\frac{1}{2}(1-h)$ and of width $\frac{1}{2}(1-h)\sqrt{h}$. The area is effectively $\frac{\sqrt{h}}{8}$ at the highest order.
Next consider the overlap with the Voronoi cell of $(\epsilon,0)$. That Voronoi cell is essentially rectangular with a width of $\frac{\sqrt{h}}{2}+\epsilon$ and a height of $1/2$. So that cell has an area of effectively $\frac{\sqrt{h}}{4}$.
The area of overlap with the Voronoi cell for $(1,0)$ is $O(h^{\frac{3}{2}})$ since it has a thickness of $O(h)$ vertically and extends a distance $O(h^{\frac{1}{2}})$ horizontally.
Thus the relevant contributions come from $(pm\epsilon,0)$ and $(\pm\sqrt{h},0)$ yielding,
\begin{align}
P^*(0,h) &\approx \frac{P(-\sqrt{h},0) + P(\sqrt{h},0)}{6} + \frac{P(-\epsilon,0) + P(\epsilon,0)}{3}\\
& = \frac{\sqrt{h}}{6} + \frac{2\epsilon}{3} \approx \frac{\sqrt{h}}{6}.
\end{align}
Putting it all together and excluding some constants,
$$
\frac{P^*(0,h) - P^*(0,0)}{h} \approx \frac{\sqrt{h}}{h} \approx \frac{1}{\sqrt{h}}.
$$
So as $h\rightarrow 0$ and more concretely selecting $\epsilon\rightarrow 0$ much faster than $h$ (probably $\epsilon = h^2$ is sufficient), we see that the Lipschitz estimate cannot hold.
Now we return to a positive claim. As long as the interpolation points are well spaced, then bounds on the gradient of the weight functions $\omega_i$ can be established immediately implying the Lipschitz property. The spirit of this estimate can be found in 1. However, the context there is a little strange (interpolation just based on the vertices of a convex polygon), but very similar ideas can be applied to bound the more general interpolant. The requirement is that the points are well-spaced. More specifically, that they generate a good aspect ratio Voronoi diagram in the sense of 2, i.e., there is a bound on the ratio of distance of the nearest and furthest neighbors of any interpolation point in the Voronoi diagram. Note that the long skinny Voronoi cells of $(\pm \epsilon, 0)$ in the counterexample do not satisfy that requirement.
1 Gillette, A.; Rand, A.; Bajaj, C., Error estimates for generalized barycentric interpolation, Adv. Comput. Math. 37, No. 3, 417-439 (2012). ZBL1259.65013 pdf.
2 Hudson, B.; Miller, G.; Phillips, T., Sparse Voronoi refinement. In Proceedings of the 15th International Meshing Roundtable 339-356 (2006). pdf
