9

Possible Duplicate:
Create random number within an annulus

I would like to obtain a uniformly obtained random point within an annulus, that is, the area that lies inside a circle of radius R1, but outside a circle of radius R2, where R1 > R2 and both circles are centered at the same point. I would like to avoid using rejection sampling.

If possible, I would like the solution to be similar to this one —used for calculating random points within a circle— which I find it extremely elegant and intuitive. That is, I would also like to avoid using the square root.

Community
  • 1
  • 1
Ricardo Sanchez-Saez
  • 9,466
  • 8
  • 53
  • 92
  • 1
    Thanks for the edit, didn't know that area was called an annulus. ;-) – Ricardo Sanchez-Saez Oct 25 '12 at 09:33
  • Added clarification noting that I would like to avoid the square root (hence, not exactly a duplicate of the other one). – Ricardo Sanchez-Saez Oct 25 '12 at 10:35
  • I want it to be a case of changing `[0,R]` in your linked question to `[R1,R2]`, but I bet it's not that simple. Or is it? – AakashM Oct 25 '12 at 10:39
  • This question is not an exact duplicate. Unlike the other question, it asks for (and I later provided) a solution that avoids the use of the square root. – Ricardo Sanchez-Saez Oct 25 '12 at 14:46
  • To my mind you're both solving the same problem - it's just you're solving it with an extra constraint on what's allowed. – AakashM Oct 25 '12 at 14:56
  • Yes, but if I really wanted to get a specific solution that complies with the extra constraint, then your definition of "exact duplicate" would not allow me to get one on StackOverflow. Anyway, this is solved already, no actual difference whether it is closed or open. – Ricardo Sanchez-Saez Oct 25 '12 at 17:25

3 Answers3

6

Its very easy. Use polar coordinates, i.e. you generate one random value for the angular value theta, and one for the distance from the origin. As your circles are both at the same origin this gets very easy.

BUT ATTENTION: you can generate the theta value by a uniform random function, that is fine, but for the distance you cant do that, as then the points will cluster around the origin. You have to take into consideration that the perimeter of a circle grows in ^2 (you have to use the inverse which is the square root).

Using a uniform distributed random function rnd (0..1) it would be like this:

theta = 360 * rnd();
dist = sqrt(rnd()*(R1^2-R2^2)+R2^2);

EDIT: For conversion into cartesion coordinates, you just compute:

x =  dist * cos(theta);
y =  dist * sin(theta);
flolo
  • 15,148
  • 4
  • 32
  • 57
  • I doubt that typical `cos` and `sin` functions take degrees, more likely radians. But this is a working method. – Benjamin Bannier Oct 25 '12 at 09:15
  • If possible I would like to avoid sqrt() like in the answer I linked. I don't want to use rejection sampling either. – Ricardo Sanchez-Saez Oct 25 '12 at 09:15
  • Depens on you function. In case you dont use degrees but radians, just replace the 360 with 2*PI – flolo Oct 25 '12 at 09:17
  • 1
    Ah, you didnt mentioned, that you dont want sqrt. (In the link they have as additional cost one more rnd evaluation, that is usally more expensive than sqrt (when it is a good rnd-generator) - its basically the same, the sqrt transform a single uniform distributed variable to the right non-uniform distribution, while the linked solution sums up two, to get the non-uniform distribution). – flolo Oct 25 '12 at 09:32
  • I plotted your solution using Mathematica and it doesn't look uniform. See it here: http://i.imgur.com/8cc0I.png For a solution that it is uniform, see this other answer: http://stackoverflow.com/a/9048443/269753 , and drawn here: http://i.imgur.com/7L04w.png – Ricardo Sanchez-Saez Oct 25 '12 at 10:32
  • you are right, the Radius must be into the sqrt, I corrected it now. – flolo Oct 25 '12 at 10:47
  • Sorry but the new one is no annulus, just regular circle uniform probability: http://i.imgur.com/ljBDd.png For two (hopefully) correct solution, see my answer that I just added. Thanks anyway for trying though! – Ricardo Sanchez-Saez Oct 25 '12 at 11:27
  • @rsanchezsaez: You got the formula wrong in your code, the rnd() must be only multiplied with (R1^2-R2^2) but you multiply also the +R2^2 with it. – flolo Oct 25 '12 at 15:11
  • @flolo: You are right! Sorry about that. Here it is the new draw: http://i.imgur.com/zxsqk.png Also, I think this is a reformulation of this previous answer: http://stackoverflow.com/a/9048443/269753 – Ricardo Sanchez-Saez Oct 25 '12 at 17:18
2

EDIT: Please note that this solution might not be uniform. See comments by Mark Dickinson below.

Ok, I think I figured it out. Note that this solution is heavily inspired in this answer, and that r1 = R1/R1 and r2 = R2/R1.

Pseudo-code:

t = 2*pi*random()
u = random()+random()
r = if u>1 then 2-u else u
r = if r<r2 then r2+r*((R1-R2)/R2) else r
[r*cos(t), r*sin(t)]

Here it is in Mathematica.

f[] := Block[{u, t, r}, u = Random[] + Random[];
r1 = 1; r2 = 0.3;
t = Random[] 2 Pi;
r = If[u > 1, 2 - u, u];
r = If[r < r2, r2 + r*((R1 - R2)/R2), r];
{r Cos[t], r Sin[t]}]

ListPlot[Table[f[], {10000}], AspectRatio -> Automatic]

Distribution of the simple algorithm

What it does is to remap all the numbers falling inside the inner circle into the annulus, spreading them evenly. If somebody finds a problem regarding the uniformity of this solution please comment.

Compare with this other solution found here:

Distribution of the sqrt algorithm

Community
  • 1
  • 1
Ricardo Sanchez-Saez
  • 9,466
  • 8
  • 53
  • 92
  • 5
    This is a neat answer, but unfortunately the result isn't uniform. If you try with e.g., `R2 = 0.95` and `R1 = 1.0`, you'll see that points nearer the outside of the annulus are favoured. A correct way to generate a suitable `r` is to take random variables `u` and `v` with `u` chosen uniformly in `[R2, R1]` and `v` chosen uniformly in `[0, R2 + R1]`, and then take `r = u` if `v < u` and `r = R2 + R1 - u` if `v >= u`. – Mark Dickinson Oct 25 '12 at 17:54
  • Mark, thanks for your comment and for your alternative. Are you sure my solution is non-uniform? I can't notice it even in with R2=0.95 as you suggested: http://i.imgur.com/wmDEo.png – Ricardo Sanchez-Saez Oct 26 '12 at 08:38
  • Pretty sure, yes. :-) Here's a test you can do: take `R2=0.5` and `R1=1.0` and generate 10,000 samples (say). Now divide your annulus into two by adding a circle of radius 0.75, and count how many of the samples you generated lie in the outer annulus (0.75 <= radius <= 1.0), and how many in the inner annulus (0.5 <= radius <= 0.75). If the solution is uniform, you should get a ratio close to 7/5=1.4 (area of outer annulus = 7/16 pi; area of inner ring = 5/16 pi). I haven't tested, but I'm predicting that your solution will give a ratio closer to 5/3 = 1.66666... – Mark Dickinson Oct 26 '12 at 19:32
1

The easiest way to do this is to use rejection sampling. Generate a large number of points uniformly in a square of side length 2*R2, then filter those samples to those inside the outer circle and not in the inner circle.

Not pretty or efficient, but in most cases, sufficient.

Andrew Walker
  • 40,984
  • 8
  • 62
  • 84
  • 3
    He said he didn't want to use rejection sampling. And for rejection sampling one could at least start with points uniform in the larger circle with the method he linked to. – Benjamin Bannier Oct 25 '12 at 09:06