Is it possible to determine the Delaunay cell containing a given point efficiently?

Yes, it is possible to find the nearest point set associated with the Delaunay cell containing a given point $c\in R^d$. As we discussed in Section 3.3, the Delaunay complex can be represented by the convex hull of appropriately lifted points in $R^{d+1}$, and the non-vertical facets coincide with the Delaunay cells once they are projected to the original space. Thus the problem of determining the Delaunay cell containing a given point $c$ can be reduced to finding the first facet of a polyhedron ``shoot'' by a ray.

To be more precise, let $f(x) = x_1^2 + \cdots + x_d^2$, and let $\tilde{p}= (p, f(x)) \in R^{d+1}$ for $p \in S$. Then the lower hull $P$ of the lifted points $\tilde{S}:=\{\tilde{p}: p\in S\}$:

$$P = conv(\tilde{S}) + nonneg(e^{d+1})$$

represents the Delaunay complex. Here $e^{d+1}$ is the unit vetor in $R^{d+1}$ whose last component is $1$. For any vector $\tilde{y} \in R^{d+1}$ and $y_0 \in R$, let $\tilde{y}^T x \ge -y_0$ denote a general inequality of a variable vector $x \in R^{d+1}$. For such an inequality to represent a valid inequality of $P$ (see Section 2.2), it must be satisfied by all points in $\tilde{S}$:

$$\tilde{y}^T \tilde{p} \ge -y_0, \forall \tilde{p} \in \tilde{S},$$ (9)

and by any points shifted vertically upwards, i.e.

$$\tilde{y}^T ( \tilde{p} + \alpha e^{d+1}) \ge -y_0, \forall \tilde{p} \in \tilde{S}$$    and any $$\alpha \ge 0.$$

Under the first inequality (9), the last inequality is equivalent to

$$\tilde{y}_{d+1} \ge 0.$$ (10)

Now every Delaunay cell is a projection of a non-vertical facet of $P$. We are thus looking for an inequality $\tilde{y}^T x \ge -y_0$ satisfying (9), (10) and $y_{d+1}\neq 0$. By scaling with $\tilde{y}_{d+1}>0$, we may assume $\tilde{y}_{d+1} = 1$. For a given point $c$, let $\tilde{c} = (c, 0)^T$, and let $L(\lambda)= \tilde{c} + \lambda \; e^{d+1}$, $\lambda \ge 0$. Determining the Delaunay cell containing $c$ is equivalent to finding the last inequality ``hit'' by the halfline $L$. More precisely, it is to find a non-vertical facet inequality such that the intersecion point of the corresponding hyperplane $\{ x : y^T x = -y_0\}$ and the half line $L(\lambda), \lambda \ge 0$ is highest possible.

By substituting $L(\lambda)$ for $x$ in $y^T x = -y_0$ with $\tilde{y}_{d+1} = 1$, we obtain

$$\lambda = -y_0 - y^T c,$$

where $y$ denotes the vector $\tilde{y}$ without the last coordinate $\tilde{y}_{d+1}$. The LP formulation is therefore:

$$z := y_0 + y^T c$$ minimize (11)

$$f(p) + y_0 + y^T p \ge 0 p \in S.$$ subject to

While an optimal solution $(y_0, y)$ to this LP may not determine any facet in general, the simplex method always returns an optimal basic solution which determines a facet inequality in this case. The Delaunay cell containing $c$ is the one determined by the set of points in $S$ whose corresponding inequalities are satisfied by equality at the optimal solution. If the LP solution is not degenerate, the dual variables that are positive at the dual optimal solution coincides with the former set.

It is important to note that the above LP might be unbounded. If it is unbounded, it can be easily shown that $c$ is not in any Delaunay cell, i.e., not in the convex hull of $S$. A certificate of unboundedness actually induces a hyperplane strongly separating $c$ from $S$. (why?)