# How can I create randomly scattered points between two circles?

13 views (last 30 days)
Sim on 28 Dec 2023
Edited: David Goodmanson on 17 Jan 2024
How can I create randomly scattered points between two circles?
I would like to have 3 inputs (assuming that the two circles are concentric, with center at (0,0)):
1. the number of scattered points
2. the radius of the inner circle
3. the radius of the outer circle
% create two concentric circles with radii r=1 and r=2
figure('Color',[1 1 1]);
for r = 1:2
th = 0:pi/50:2*pi;
xc = r * cos(th);
yc = r * sin(th);
hold on
plot(xc,yc,'color','r','LineWidth',2)
end
axis equal
set(gca,'Visible','off')
Sim on 28 Dec 2023
Edited: Sim on 28 Dec 2023
I am not sure if this is correct:
% create the concentric circles
figure('Color',[1 1 1]);
for r = 1:2
th = 0:pi/50:2*pi;
xc = r * cos(th);
yc = r * sin(th);
hold on
plot(xc,yc,'color','r','LineWidth',1)
end
% generate N random numbers in the interval (a,b) with the formula r = a + (b-a).*rand(N,1)
N=500;
t = 2 * pi * rand(N,1);
g = 1 + (2-1).*rand(N,1);
x = g.*cos(t);
y = g.*sin(t);
plot(x,y,'o','MarkerFaceColor','b','MarkerEdgeColor','b')
axis equal
set(gca,'Visible','off')

Adam Danz on 28 Dec 2023
Edited: Adam Danz on 29 Dec 2023
a = 4;
b = 10;
N = 1000;
r = a + (b-a) * rand(N,1);
th = 2*pi * rand(N,1);
x = r.*cos(th);
y = r.*sin(th);
where a<=b.
Another way to think about it is in polar coordinates
[x, y] = pol2cart(th, r);
Update
As @John D'Errico points out below, the solution above will not produce a uniform distribution of coordinates in 2D space. This is because points will be compressed as the radii decreases within the disc. To compensate for this, you could use inverse transform sampling to generate random r values from a nonuniform distribution that accounts for the increasing circumference size.
r2 = a + (b-a) * rand(N,1).^.5;
x2 = r2.*cos(th);
y2 = r2.*sin(th);
Compare these results
The first method clusters points near the inner radius. The second method is closer to uniform in 2D space but has a negative bias arond the inner radius. See comments below for more on this topic.
figure()
tiledlayout(2,2)
plot(nexttile(),x,y,'b.')
title('Uniform th and r sampling')
axis equal tight
grid on
rectangle('Position',[-a -a 2*a 2*a],'Curvature',1,'EdgeColor','r')
rectangle('Position',[-b -b 2*b 2*b],'Curvature',1,'EdgeColor','r')
plot(nexttile(),x2,y2,'b.')
title('Inverse transform sampling')
axis equal tight
grid on
rectangle('Position',[-a -a 2*a 2*a],'Curvature',1,'EdgeColor','r')
rectangle('Position',[-b -b 2*b 2*b],'Curvature',1,'EdgeColor','r')
N = 100000;
r = a + (b-a) * rand(N,1);
th = 2*pi * rand(N,1);
x = r.*cos(th);
y = r.*sin(th);
r2 = a + (b-a) * rand(N,1).^.5;
x2 = r2.*cos(th);
y2 = r2.*sin(th);
histogram2(nexttile(),x,y,20,'DisplayStyle','tile','EdgeColor','none')
colorbar
axis equal
histogram2(nexttile(),x2,y2,20,'DisplayStyle','tile','EdgeColor','none')
colorbar
axis equal
Adam Danz on 30 Dec 2023
Edited: Adam Danz on 30 Dec 2023
Great discusion here. @Paul and @David Goodmanson, please consider moving your solutions to a new answer to increase visibility. I can't un-accept my answer at the moment and there are much bettern approaches in this thread.
Sim on 30 Dec 2023
Edited: Sim on 30 Dec 2023
Many many thanks to everyone!
I did not expect that such a question could generate so many interesting comments/discussions/answers!
In these days I am not able to read carefully all your answers, but I will do it as soon as possible!
Meanwhile: have a great start into the New Year!

David Goodmanson on 30 Dec 2023
Edited: David Goodmanson on 17 Jan 2024
What tells the tale is the differential area element in the (r, theta) plane, which is r dr dtheta.
Picking theta from a uniform random variable on 0<theta<2pi and letting x = r*cos(theta), y = r*sin(theta) takes care of the angle part in the usual way. That leaves picking r.
To fill the plane uniformly in r dr with draws from a uniformly distributed variable rho, then those two quantities must be proportonal, so
r dr = c drho ***
for some constant c. Everything else follows from that expression.
Now
r dr = d(r^2/2) = c drho and consequently
r^2/2 = c rho + d (4)
for a couple of constants c and d. For an annulus of radii a and b,
r = a @ rho = 0 and
r = b @ rho = 1
which fixes c and d, and
r^2 = a^2 + (b^2-a^2) rho so
r = sqrt(a^2 + (b^2-a^2) rho) (3)
Then for the choosing of N random points,
r = sqrt(a^2 + (b^2-a^2)*rand(1,N)) ***
which is the correct result and is one of the previously suggested solutions. What does not work are two other previous suggestions,
r = a + (b-a) sqrt(rho) % Incorrect (2)
r^2/2 = (1/2) ( a^2 + 2a(b-a)sqrt(rho) + (b-a)^2 rho )
r = a + (b-a) rho % Incorrect (1)
r^2/2 = (1/2) (a^2 + 2a(b-a) rho + (b-a)^2 rho^2 )
neither of which is of the required form (4).
It's interesting to take a look at the point densities that are produced by these three cases. The density is, for N points (assuming angular symmetry),
sigma = dn/dA = N drho / (2pi r dr) = N / ( 2pi r dr/drho )
Ignoring the N/2pi factor for comparison purposes,then
rho = 0:.001:1;
a = 1;
b = 2;
% case 1
r1 = a+(b-a)*rho;
drdrho1 = (b-a);
sigma1 = 1./(r1.*drdrho1);
% case 2
r2 = a+(b-a)*sqrt(rho);
drdrho2 = (1/2)*(b-a)./sqrt(rho);
sigma2 = 1./(r2.*drdrho2);
% case 3
r3 = sqrt(a^2+(b^2-a^2)*rho);
drdrho3 = (1/2)./sqrt(a^2+(b^2-a^2)*rho) * (b^2-a^2);
sigma3 = 1./(r3.*drdrho3);
plot(r1,sigma1,r2,sigma2,r3,sigma3)
xlabel('r')
ylabel('density' )
legend('[1] a+(b-a)*rho','[2] a+(b-a)*sqrt(rho)', ...
'[3] sqrt(a^2+(b^2-a^2)*rho)','location','south')
As has been previously mentioned, case 1 throws too many points to the inside, and case2 throws too many points to the outside. Case 3 is just right.
oo As a side note, N points go into a specified area. Replacing the draws from the two random uniform variables with integrals over rho and theta, then for case (3), and using 2pi for the dtheta integral due to angular symmetry
r = sqrt(a^2 + (b^2-a^2) rho)
dr = (1/2) (b^2-a^2) 1/sqrt(a^2 + (b^2-a^2) rho) drho
rdr = (1/2) (b^2-a^2) drho
Integral (rdr dtheta) = (1/2) (b^2-a^2) (2pi) = pi (b^2-a^2) % annulus
which is the correct area. That's not a ringing endorsement of (3), though, since *any* well behaved function r = f(rho) such that f(0) = a and f(1) = b (including (1) and (2) ), will produce the same result, albeit with an incorrect density of points.
John D'Errico on 30 Dec 2023
Edited: John D'Errico on 30 Dec 2023
I think the incorrect form stems from an intuitive leap, because when a==0, so the circle will be completely filled, then this
r = a + (b-a)*sqrt(rho)
reduces to
r = b*sqrt(rho)
And that form should fill the circle in a uniform way.
As well, the generation of random points across a triangle uses a similar form with the sqrt. So this is a not uncommon idea. It just does not work for an annulus with a>0.
Sim on 16 Jan 2024
Edited: Sim on 16 Jan 2024
@Adam Danz, @John D'Errico, @Paul, @Torsten many many thanks to all of you for all your answers and comments, that I found extremely interesting and insightful! I have really appreciated them!
.....and now, which one to accept ? This is quite embarassing, since all of your answers are top-notch solutions, and they all deserve to be accepted and used by the community! :-) :-) :-)

John D'Errico on 30 Dec 2023
Edited: John D'Errico on 30 Dec 2023
Let me explain it differently than the rest, with a simple, intuiitive derivation.
Again, the problem is perfectly symmetric in polar angle, so that part is trivial. All we need to understand is the distribution in the radius, r.
Assume you wish to see a uniform sampling on the annulus from a to b in radius. Now consider any fixed radius r, so now just a circle of radius r. We want to see proportionally 2*pi*r*d points around that infinitessimally narrow circular annulus. Essentially, this gives us the radial distribution. The PDF is simply:
syms r d a b
PDF = 2*pi*r*d;
Here d is the packing density, and r will vary from a to b. We need to normalize the pdf so it has unit area on that domain.
PDF = PDF/int(PDF,r,[a,b])
PDF =
I hope that makes sense, in a very intuitive way. Now that we have the PDF, we can compute the CDF of the radial distribution, as simply the integral over r.
syms R
CDF = int(PDF,r,[a,R])
CDF =
Again, r will live on the interval [a,b]. And also see that I used the definite integral from a to the radius R. This is important!
For example, with a=1, b=10, we would have:
fplot(subs(PDF,[a,b],[1,10]),[1,10])
fplot(subs(CDF,[a,b],[1,10]),[1,10])
As you can see, the PDF is linear, and the CDF is a simple quadratic polynomial.
How do we sample from a distribution when the CDF is known, and it has a simple inverse? This itself is simple enough. We use the inverse of the CDF. That is, we choose a UNIFORM random sampling on the interval [0,1]. rand does this perfectly. Then we solve for R.
syms u
CDFinv = solve(CDF == u,R)
CDFinv =
We clearly want the second solution there, since it is the positive one.
CDFinvfun = @(u,a,b) sqrt((b^2-a^2)*u + a^2)
CDFinvfun = function_handle with value:
@(u,a,b)sqrt((b^2-a^2)*u+a^2)
I just wrote it as a function, since matlabFunction did something screwy to me, and it is too early in the AM right now for me to think.
N = 1000000;
Usample = rand(N,1);
Rsample = CDFinvfun(Usample,1,10);
histogram(Rsample,100,'norm','pdf')
And you see a perfectly linear sample PDF drop out ion the interval [a,b]. Now lets try it out.
N = 10000;
Usample = rand(N,1);
Rsample = CDFinvfun(Usample,1,10);
theta = rand(N,1)*2*pi;
x = Rsample.*cos(theta);
y = Rsample.*sin(theta);
plot(x,y,'.')
axis equal
And that is exactly what we want to see. A nice uniform sampling on the annulus.
If what I did does not make sense, then go through it step by step. I tried to be as complete as possible in my Answer, but I am sure I went too fast at some point.
Addendum: Some additional points of interest I may have glossed over above,
1. The entire solution stems directly from the assumption of a uniform scattering on the annulus. That directly implies the PDF, and then the rest follows.
2. The use of a definite integral to compute the CDF is important, since if you just compute an indefinite integral with no other flags to int, the constant of integration will be lost. The definite integral implicitly builds the (correct) constant of integration into the CDF. I think this is an easily missed point.
3. The above analysis can also trivially be extended to 3 or more dimensions, so with annular spheres. There the transformation will be different, but as easily derived based on the resulting polynomial PDF on the finite annular domain.
4. If you follow the solution I've posed, you can see that whan a==0, and we want to simply sample uniformly over the entire circle of radius b, then the intuitive square root transformation that Adam proposed now appears. But that clearly fails when a>0. When a==0, the inverse CDF reduces to simply b*sqrt(u),
Paul on 30 Dec 2023
Hi John,
Just came to point out that finverse is another option to get the inverse of the CDF
syms R a b u
CDF = -(R^2-a^2)/(a^2-b^2)
CDF =
CDFinv = subs(finverse(CDF,R),R,u)
iCDF =