Take the 2-minute tour ×
Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

Given a bounded $A \times B$ rectangle with a set of chosen coordinates, generated for example with the command:

A = 1;
B = 1;
randPoints = Table[{RandomReal[{0,A}],RandomReal[{0,B}]},{k,1,10^4}];
ListPlot[randPoints] 

For each point in randPoints, I would like to find the center and radius of the largest circle that contains that point but no other point from the randPoints set, and does not exceed the boundaries of the $A \times B$ box.

How might one do this in Mathematica v9.0? In the example there are 10^4 points, but I'm also envisioning much larger sets, so ideally I'd like a fast routine. This could, however, constitute pushing my luck.

The problem is trivial if one locks the center of the circle on the relevant element in randPoints (just use the Nearest function to find the circle radius), but it seems like an additional trick is required if one only requires this element to be by itself somewhere inside the circle.

share|improve this question
3  
I feel like this question is more about computational geometry than about implementation in Mathematica, and might be more appropriate at the mathematics StackExchange. Anyhow, the desired circle must be centered at a vertex of the Voronoi diagram of all the points in randPoints other than the point of interest; but which vertex, that's the hard part... –  Rahul Sep 4 '13 at 8:50
    
@Kuba: No, if you're on a Voronoi edge then you are equidistant from two nodes, and you can always move along the edge to increase your distance to them; therefore it cannot be a local maximum of the distance-to-nearest-node function, which is the same as the how-big-a-circle-can-I-draw-centered-at-this-point function. Local maxima occur at Voronoi vertices. –  Rahul Sep 4 '13 at 9:12
    
@SeptemberGrass I tend to agree with Rahul on this. Would you like me to migrate this question to Mathematics? Once you have your algorithm you could return and ask how best to implement it in Mathematica. –  Mr.Wizard Sep 4 '13 at 9:37
    
For this type of problems you can use Genetic Algorithms. In Mathematica you can formulate the problem as optimization problem and use Maximization / Minimization with Differential Evolution method. see the example and the C code from ai-junkie.com/ga/intro/gat3.html –  user840 Sep 4 '13 at 9:39
    
Just to clarify, the circles will be overlapping, right? –  gpap Sep 4 '13 at 9:45
show 2 more comments

migrated from mathematica.stackexchange.com Sep 4 '13 at 10:02

This question came from our site for users of Mathematica.

2 Answers

Let the point we want the circle to contain be $\mathbf p_0$, and the other points be $\mathbf p_1, \ldots, \mathbf p_n$. If the circle is centered at $\mathbf q$, the largest radius it can have is $$r(\mathbf q) = \min_{i=1}^n\|\mathbf q-\mathbf p_i\|.$$ For the circle to contain $\mathbf p_0$, we need $\|\mathbf q-\mathbf p_0\|\le r(\mathbf q)$, which if you plug in the definition of $r(\mathbf q)$ just means that $\mathbf q$ is closer to $\mathbf p_0$ than any other $\mathbf p_i$; in other words, it lies in the cell corresponding to $\mathbf p_0$ in the Voronoi diagram of $\{\mathbf p_0, \mathbf p_1, \ldots, \mathbf p_n\}$. So, our problem reduces to finding the maximum of $r(\mathbf q)$ within the Voronoi cell of $\mathbf p_0$, which I'll call $C$.

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 $\mathbf 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$.

So, here is an algorithm:

  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 $\mathbf 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 $\mathbf 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 $\mathbf q$ and radius $r(\mathbf q)$.

Implementation notes: Since you want to do this for each point one by one, you don't have to repeat step 1, which ought to be the most expensive step. Also, in steps 3 and 5, you already have a Voronoi diagram, so you can evaluate $r(\mathbf q)$ without actually iterating over all $n$ points.

share|improve this answer
add comment

enter image description here

enter image description here

This is not efficient but works:

near = {#, Nearest[randPoints, #, 2][[2]]} & /@ randPoints;

gr = {Red, Circle[#[[1]], 0.99 EuclideanDistance[#[[1]], #[[2]]]]} & /@
    near;

Show[lp, Graphics[gr]]

where lp is the ListPlot

share|improve this answer
1  
If I'm not mistaken, all of the circles seem centered on each point? –  SeptemberGrass Sep 4 '13 at 9:32
2  
That isn't what is needed though. Your circles are centered at each point and are not the maximal ones. –  gpap Sep 4 '13 at 9:32
    
The visualization might be more informative using less points. –  Yves Klett Sep 4 '13 at 9:46
    
Less points and AspectRatio->1 or the circles will show up as ovals. –  Pickett Sep 4 '13 at 10:13
    
I am sorry for missing the deeper question, i.e. maximum circle containing point and no other point, not maximum circle around a point. I accept completely the visualization issues viz aspect ratio->1, and the incomprehensive nature of 10000 circles and points. Look forward to the answer. –  ubpdqn Sep 4 '13 at 10:29
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.