# Implementing an algorithm for finding the largest circle that contains a single point in a set (and no other point)

This question concerns the implementation of an algorithm proposed by Rahul Narain on a former question of mine that was migrated to math.stackexchange: http://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:

1. Construct the Voronoi diagram $\mathcal V$ of $\{\mathbf p_0,\mathbf p_1,\ldots,\mathbf p_n\}$.

2. Let $C$ be the cell of $\mathcal V$ corresponding to $p_0$.

3. Evaluate $r(\mathbf q)$ for all vertices of $C$.

4. 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$.)

5. Evaluate $r(\mathbf q)$ at all vertices of $\mathcal V'$ that lie inside $C$.

6. 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_] := GraphicsMeshPointWindingNumber[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.

-
@cormullion Re: spelling mistake in the title, oops, thanks. –  SeptemberGrass Sep 5 '13 at 8:39
I tried playing with this yesterday (for the record I fully agree with the method proposed) but 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? –  gpap Sep 5 '13 at 10:38
@gpap Right, I had the same experience. My actual point set is something like 10^6, but I can of course select only points within a certain threshold distance of some $p_0$ point. Still, VoronoiDiagram is quite slow. –  SeptemberGrass Sep 5 '13 at 11:13
@gpap It's a neat technique - I was guessing that one would have to settle for some approximation scheme due to pathological cases with multiple local maxima. Or that I would at least have to run a numerical optimization with a guarantee of convergence. Looks like that guess was wrong. –  SeptemberGrass Sep 5 '13 at 11:26
I don't know if this helps, but there is a GraphicsMesh`Voronoi which is about 3 orders of magnitude faster than VoronoiDiagram. It returns only the coordinates of the cell vertices. –  Simon Woods Sep 5 '13 at 19:33

Here's a direct implementation of Rahul Narain's algorithm using Version 10 functions (no attempt has been made to optimize):

SeedRandom[80]
pts = RandomReal[5, {50, 2}];
vor = VoronoiMesh[pts, {{0, 5}, {0, 5}}];

Let's pick a point to draw it's largest circle:

poi = pts[[9]];

The point in red is our point of interest:

Graphics[{GraphicsComplex[MeshCoordinates[vor], {Thick, Blue, MeshCells[vor, 1],
Opacity[0.2], Yellow, MeshCells[vor, 2]}], PointSize[0.013],
Point[pts], Red, Point[poi]}]

Here we go:

Let's find the Voronoi cell enclosing our point of interest (poi)

pointCell[pt_, reg_] := Pick[reg, RegionMember[#, pt] & /@ reg][[1]];
cell = pointCell[poi, MeshPrimitives[vor, 2]];

We create a NearestFunction for all point set excluding our poi

nf = Nearest[Complement[pts, {poi}]];

We apply this NearestFunction to all vertices of the enclosing Voronoi cell

rq = (nf /@ cell[[1]])[[All, 1]];

Next, we compute the Voronoi diagram of the point set excluding our poi

vorc = VoronoiMesh[Complement[pts, {poi}], {{0, 5}, {0, 5}}];

The following picks out the vertices of vorc that is within the Voronoi cell enclosing our poi

rm = RegionMember[cell];
vert = With[{pp = MeshPrimitives[vorc, 0][[All, 1]]}, Pick[pp, rm[pp]]];

Next, we apply nf to these vertices

rq2 = (nf /@ vert)[[All, 1]];

The following lines of code computes the radius and center of the desired circle:

allcent = Join[cell[[1]], vert];
allvert = Join[rq, rq2];
tp = Transpose[{allcent, allvert}];
radius = Max[dist = EuclideanDistance[#1, #2] & @@@ tp];
cen = (tp[[Ordering[dist, -1]]] // Flatten[#, 1] &)[[1]];

Visualize:

Graphics[{GraphicsComplex[MeshCoordinates[vor], {Thick, Blue, MeshCells[vor, 1],
PointSize[0.02], Opacity[0.2], Yellow, MeshCells[vor, 2]}], PointSize[0.013],
Point[pts], Red, Point@poi, Thick, Purple, Circle[cen, radius]}]

Here's the complete code (some optimizations are possible e.g. computing the Voronoi diagram once outside this function):

largestC[p_, points_] :=
Module[{vor = VoronoiMesh[points], ptscomp = Complement[points, {p}],
vorc, nf, cell, rm, vert, nearp, radius, rq, rq2, allvert,
allcent, tp, cen, dist, pointCell},
pointCell[pt_, reg_] := Pick[reg, RegionMember[#, pt] & /@ reg][[1]];
vorc = VoronoiMesh[ptscomp];
nf = Nearest[ptscomp];
cell = pointCell[p, MeshPrimitives[vor, 2]];
rm = RegionMember[cell];
rq = (nf /@ cell[[1]])[[All, 1]];
vert = With[{pp = MeshPrimitives[vorc, 0][[All, 1]]}, Pick[pp, rm[pp]]];
rq2 = (nf /@ vert)[[All, 1]];
allcent = Join[cell[[1]], vert];
allvert = Join[rq, rq2];
tp = Transpose[{allcent, allvert}];
radius = Max[dist = EuclideanDistance[#1, #2] & @@@ tp];
cen = First[tp[[Ordering[dist, -1]]] // Flatten[#, 1] &];
]

Here's the animation code:

Manipulate[
Graphics[{GraphicsComplex[MeshCoordinates[vor], {Thick, Blue, MeshCells[vor, 1],
PointSize[0.02], Opacity[0.2], Yellow, MeshCells[vor, 2]}],
PointSize[0.005], Point[pts], PointSize[0.015], Red, Point@pts[[k]],
Thin, Purple, largestC[pts[[k]], pts]}], {k, 1, Length @ pts, 1}]
-
+1 It's a bit of a relief to see that my algorithm actually works! "I have only proved it correct, not tried it." --Don Knuth –  Rahul Narain Aug 29 at 18:12
@RahulNarain Thanks. And good job on coming up with that algorithm, it made my job easier :) –  RunnyKine Aug 29 at 18:16

Problems like this may benefit from getting rid of the "non-linearity", where the circle centres are fixed and the radii are to be determined. If this formulation is fast enough (is it ever?), then it may be run several times, refining the grid of fixed circle centres until the radii do not change significantly.

The code below is a first pass, for proof of concept. If the algorithm is feasible for the full problem the OP has in mind, it will certainly benefit from optimization by more skilled users.

The algorithm is as follows. Generate a set of random points $p$.

SeedRandom[1]; p = RandomReal[{0, 1}, {10, 2}];

Generate a grid (or even a random selection) of fixed circle centres $c$.

c = Flatten[Table[{i, j}, {i, 0., 1., 0.1}, {j, 0., 1., 0.1}], 1]

Form the nearest function $f$ of the points $p$.

f = Nearest[p]

For each circle centre $c_i$, find the nearest two points $p_1$ and $p_2$ to $c_i$, then form the circle with centre $c_i$ and radius just less than the distance from $c_i$ to the second nearest point to $c_i$. This circle contains the nearest point $p_1$ but just misses containing the second nearest point $p_2$. Label this circle with the point $p_1$ it contains.

u = Map[({p1,p2}=f[#,2]; {p1, Circle[#, 0.999*Norm[p2-#]]})&, c]

There are multiple circles containing any particular point $p_j$. Gather those circles containing each $p_j$.

g = GatherBy[u, #[[1]] &]

Pick from each gather the circle with the largest radius.

v = Map[Sort[#, #1[[2, 2]]<#2[[2, 2]]&]&, g][[All,-1]]

Displaying Sort[v] shows the contained point and the circle surrounding it. The Sort function is used to place the random points $p$ in a standard order for comparison of circle radii from refined grids of fixed circle centres.

Ten thousand random points $p$ and 40000 circle centres $c$ took about 1 second on my machine.

-