This question concerns the implementation of an algorithm proposed by Rahul Narain on a former question of mine that was migrated to math.stackexchange: https://math.stackexchange.com/questions/483845/finding-the-largest-circle-that-contains-a-single-point-in-a-set-and-no-other-p/484600?noredirect=1#484600
Here are the details for the algorithm I'd like to implement -
Let $r(\mathbf q) = \min_{i=1}^n\|\mathbf q-\mathbf p_i\|$ be the largest radius of a circle centered at the coordinate $q$ containing some point $p_0$ (thus satisfying the requirement that $\|\mathbf q-\mathbf p_0\|\le r(\mathbf q)$) and containing no other point in the set $(p_0, p_1, ...) \in P$. We wish to find $q$ and $r(\mathbf q)$ for every point $p_i \in P$.
To do so, Rahul Narain proposed the following steps:
Construct the Voronoi diagram $\mathcal V$ of $\{\mathbf p_0,\mathbf p_1,\ldots,\mathbf p_n\}$.
Let $C$ be the cell of $\mathcal V$ corresponding to $p_0$.
Evaluate $r(\mathbf q)$ for all vertices of $C$.
Construct the Voronoi diagram $\mathcal V'$ of $\{\mathbf p_1,\ldots,\mathbf p_n\}$. (There should exist efficient algorithms for removing a single point $p_0$ from a previously computed Voronoi diagram $\mathcal V$.)
Evaluate $r(\mathbf q)$ at all vertices of $\mathcal V'$ that lie inside $C$.
Choose the point from steps 3 and 5 with the largest $r(\mathbf q)$. The desired circle has center $q$ and radius $r(\mathbf q)$.
The idea is that, by the definition of $r(\mathbf q)$ the center of the circle containing some $p_0 \in P$ is going to be closer to $p_0$ than any other point in the set $P$. Thus, we know that $q$ is going to fall in the Voronoi cell (which is a convex polygon) of $p_0$ since, by definition, this cell $C$ contains the set of points closer to $p_0$ than any other point in $P$.
Rahul Narain then argues:
"There are two possibilities: either the maximum lies in the interior of $C$, or it lies on its boundary. In the former case, it must be an unconstrained local maximum of $\min_{i=1}^n\|\mathbf q-\mathbf p_i\|$; any such maximum is a vertex of the Voronoi diagram of $\{\mathbf p_1,\ldots,\mathbf p_n\}$. In the latter case, if $q$ lies on an edge, you can always increase $r(\mathbf q)$ by moving along the edge, so the maximum must occur at a vertex of $C$."
In other words, if $q$ for some point $p_0$ is inside $C$, then it's going to be a point that is as far away as possible from all other points in $P$ excluding $p_0$ (i.e. a vertex of the Voronoi diagram for $P-p_0$ falling inside the convex polygon $C$). Move away from this vertex, and you're strictly moving closer to at least one point in $P-p_0$. If $q$ lies on an edge of $C$, since $C$ is a polygon, any step along an edge will strictly move us closer to $p_0$ or further away, so $q$ must lie on a vertex of $C$.
From this argument, which makes sense to me, we know that $q$ must lie on a vertex $C$, or a vertex inside $C$ in the Voronoi diagram $\mathcal V'$ of the points $\{\mathbf p_1,\ldots,\mathbf p_n\}$ missing the point $p_0$.
Here's where I am thus far with an implementation ---
Let's first generate some random 2D coordinates:
numPoints = 10;
lowerboundPointRange = 0;
upperboundPointRange = 10;
randPoints = RandomReal[{lowerboundPointRange, upperboundPointRange}, {numPoints, 2}];
To plot the Voronoi diagram can write:
Get["ComputationalGeometry`"]
DiagramPlot[randPoints]
To return a set of vertices vdVertices for the Voronoi diagram $\mathcal V$, we can write:
VoronoiDiagram[randPoints]
Which returns something like this (if we set numPoints = 4):
{{{-5.08005, 20.1012}, {4.94272, 2.06684}, {10.0376, 2.1414}, Ray[{4.94272, 2.06684}, {2.70604, -2.84145}], Ray[{-5.08005, 20.1012}, {-9.26944, 26.3737}], Ray[{10.0376, 2.1414}, {15.9904, -0.261357}]}, {{1, {2, 1, 5, 4}}, {2, {3, 2, 4, 6}}, {3, {2, 3, 1}}, {4, {1, 3, 6, 5}}}}
We can understand this as, first, a set of vertices for the Voronoi cell $\mathcal V$, a list of rays defining edges for each Voronoi cell, and finally as a vertex adjacency list. If we just want the vertices, we can select for them with the command:
vdVertices = Select[VoronoiDiagram[randPoints][[1]], Head[#] == List &]
In terms of solving the problem of determining which vertices in $\mathcal V'$ fall in the Voronoi cell $C$ for some point $p_0$, we can use (rm -rf♦)'s answer to the question (How to check if a 2D point is in a polygon?), where he points out that one can use an internal function GraphicsMeshPointWindingNumber to check if a point is in a 2D polygon:
inPolyQ[poly_, pt_] := Graphics`Mesh`PointWindingNumber[poly, pt] =!= 0
(* Examples *)
inPolyQ[{{-1, 0}, {0, 1}, {1, 0}}, {1/3, 1/3}]
(* True *)
inPolyQ[{{-1, 0}, {0, 1}, {1, 0}}, {1, 1}]
(* False *)
For the moment I am stuck at not knowing how to do things like efficiently compute a Voronoi diagram $\mathcal V'$ missing a point from the Voronoi diagram $\mathcal V$, as Rahul Narain suggests, and also how to efficiently find the polytope vertices for $C$ given some $p_0$, and in general how to restrict the search for vertices in $\mathcal V'$ that fall in $C$. Also, there ought to be a better way to compute $r(\mathbf q) = \min_{i=1}^n\|\mathbf q-\mathbf p_i\|$ than scanning through each $p_i \in P-p_0$?
Help on completing the implementation would be much appreciated.
VoronoiDiagram[randPoints]takes a really long time to produce in my machine for even $ 10^3 $ points and you were after something efficient for orders of magnitude of $ 10^4 $ points or more, right? $\endgroup$VoronoiDiagram, though, so I can't tell you how to get the adjacency and correspondence information out of it. $\endgroup$Graphics`Mesh`Voronoiwhich is about 3 orders of magnitude faster thanVoronoiDiagram. It returns only the coordinates of the cell vertices. $\endgroup$