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: 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:

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][], 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
• After some reading I found that, with regard to $\mathcal V'$, for our purposes it is sufficient to compute the Voronoi diagram of only the points neighbouring $\mathbf p_0$ ("neighbouring" being in terms of adjacency of Voronoi cells). In step 3, $r(\mathbf q)$ is just $\|\mathbf q-\mathbf p_0\|$, while in step 5, a Voronoi vertex is adjacent to three Voronoi cells, and $r(\mathbf q)$ is the distance to any of the three corresponding points. I'm not familiar with Mathematica's VoronoiDiagram, though, so I can't tell you how to get the adjacency and correspondence information out of it. – Rahul Sep 5 '13 at 18:09
• I don't know if this helps, but there is a GraphicsMeshVoronoi 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
• If the point $\mathbf p_0$ lies on the convex hull of the point set then there are arbitrarily large circles that contain it. You can detect this by checking if its Voronoi cell $C$ is unbounded. In that case, there is no largest enclosing circle; you could say that its radius is infinite, but you can't really pick a meaningful center except maybe to say that it lies at infinity in the direction that the Voronoi cell extends. – Rahul Sep 5 '13 at 21:26 Here's a direct implementation of Rahul Narain's algorithm using Version 10 functions (no attempt has been made to optimize):

SeedRandom
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[];


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][];
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[])[[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[], vert];
allvert = Join[rq, rq2];
tp = Transpose[{allcent, allvert}];
radius = Max[dist = EuclideanDistance[#1, #2] & @@@ tp];
cen = (tp[[Ordering[dist, -1]]] // Flatten[#, 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][];
vorc = VoronoiMesh[ptscomp];
nf = Nearest[ptscomp];
cell = pointCell[p, MeshPrimitives[vor, 2]];
rm = RegionMember[cell];
rq = (nf /@ cell[])[[All, 1]];
vert = With[{pp = MeshPrimitives[vorc, 0][[All, 1]]}, Pick[pp, rm[pp]]];
rq2 = (nf /@ vert)[[All, 1]];
allcent = Join[cell[], 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 Aug 29 '14 at 18:12
• @RahulNarain Thanks. And good job on coming up with that algorithm, it made my job easier :) – RunnyKine Aug 29 '14 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; 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, #[] &]


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.

Here is an answer utilizing the fact that by definition, Voronoi mesh nodes are equidistant to the input points located in the 3 adjacent meshes. That is, the circle centered at a mesh node, with radius the distance to the nearest input point is the largest circle that contains no input points. Here is a visualization:

SeedRandom;
pts = RandomReal[10, {7, 2}];
vm = VoronoiMesh[pts];
nodes = vm["Coordinates"][[vm["InteriorVertices"]]];
Show[
vm,
] Notice how the green dots (nodes) are the centers of circles that touch 3 red dots (the input points).

Since we actually want the largest circle that contains exactly 1 point, we can start by removing that point, generate the Voronoi mesh without that point, then find the largest circle that contains the removed point. Here is the code:

largestCircle[pts_, pt_] := Module[
{other, vm, nodes, interior, radii, distances, circles, ord},
other = DeleteCases[pts, pt];
vm = VoronoiMesh[other];
(* the following are the candidate circle centers *)
nodes = vm["Coordinates"][[vm["InteriorVertices"]]];

(* The distance from the node to the nearest point is the radius of the circle *)

(* This is the distance from the point pt to each of the nodes *)
distances = DistanceMatrix[nodes, {pt}][[All, 1]];

(* look at only circles whose radius is larger than the distance to the point *)

(* pick the largest circle if there is one
* if there isn't a largest circle, than the point is an exterior point,
* and the circle is unbounded *)
If[Length[circles]>0,
ord = circles[[First @ Ordering[radii[[circles]], -1]]];
{}
]
]


And a comparison with @RunnyKine's answer:

SeedRandom;
pts = RandomReal[10,{100,2}];

r1=largestC[ pts[],pts];//AbsoluteTiming
r1

r2=largestCircle[pts, pts[]];//AbsoluteTiming
r2


{0.255573, Null}

Circle[{5.58479, 2.60399}, 0.996386]

{0.004053, Null}

Circle[{5.58479, 2.60399}, 0.996386]

My answer is significantly faster. Here is another example, where my answer differs from his due to different treatments of how to handle external points.

SeedRandom;
pts = RandomReal[10,{10,2}];

r1=largestC[ pts[],pts]
r2=largestCircle[pts, pts[]]


Circle[{9.19615, 12.2502}, 6.07274]

Circle[{6.57676, 6.96635}, 3.27242]

And here is a visualization comparing them:

Show[
VoronoiMesh[DeleteCases[pts,pts[]]],
Graphics[{r1, r2, Green, Point[pts[]], Red, Point[DeleteCases[pts,pts[]]]}]
]
` My code doesn't return circles that only touch two input points. Also, it seems that the other answer could produce a larger circle that just touches the removed point.