Azzera filtri
Azzera filtri

How can I create randomly scattered points between two circles?

16 visualizzazioni (ultimi 30 giorni)
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
axis equal
  1 Commento
Sim il 28 Dic 2023
Modificato: Sim il 28 Dic 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
% generate N random numbers in the interval (a,b) with the formula r = a + (b-a).*rand(N,1)
t = 2 * pi * rand(N,1);
g = 1 + (2-1).*rand(N,1);
x = g.*cos(t);
y = g.*sin(t);
axis equal

Accedi per commentare.

Risposte (3)

Adam Danz
Adam Danz il 28 Dic 2023
Modificato: Adam Danz il 29 Dic 2023
Your solution is fine.
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);
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.
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')
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);
axis equal
axis equal
  10 Commenti
Adam Danz
Adam Danz il 30 Dic 2023
Modificato: Adam Danz il 30 Dic 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 il 30 Dic 2023
Modificato: Sim il 30 Dic 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!
P.S. I thought to temporarily unaccept the @Adam Danz's answer (and please apologise Adam!! Your answer is really great!) just to encourage other replies and give more visibility to everyone...

Accedi per commentare.

David Goodmanson
David Goodmanson il 30 Dic 2023
Modificato: David Goodmanson il 17 Gen 2024
(At Adam's suggestion I moved my previous comment to this answer)
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.
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);
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.
  2 Commenti
John D'Errico
John D'Errico il 30 Dic 2023
Modificato: John D'Errico il 30 Dic 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 il 16 Gen 2024
Modificato: Sim il 16 Gen 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! :-) :-) :-)

Accedi per commentare.

John D'Errico
John D'Errico il 30 Dic 2023
Modificato: John D'Errico il 30 Dic 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:
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:
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);
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);
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),
  1 Commento
Paul il 30 Dic 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 = 

Accedi per commentare.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by