Generate random points inside a circle

46 views (last 30 days)
Mahesh Ramaraj
Mahesh Ramaraj on 23 Jan 2011
Edited: Torsten on 8 May 2022
I am trying to generate random points inside a circle given the radius and the center. If anyone has a sample code or can help me with this, thanks for the help.

Sign in to answer this question.

Accepted Answer

Paulo Silva
Paulo Silva on 23 Jan 2011
function [x y]=cirrdnPJ(x1,y1,rc)
%the function, must be on a folder in matlab path
a=2*pi*rand
r=sqrt(rand)
x=(rc*r)*cos(a)+x1
y=(rc*r)*sin(a)+y1
end
%to test the function
clf
axis equal
hold on
x1=1
y1=1
rc=1
[x,y,z] = cylinder(rc,200);
plot(x(1,:)+x1,y(1,:)+y1,'r')
for t=1:1000 %loop until doing 1000 points inside the circle
[x y]=cirrdnPJ(x1,y1,rc)
plot(x,y,'x')
%pause(0.01) %if you want to see the point being drawn
end
  7 Comments
Show 5 older comments
Mahesh Ramaraj
Mahesh Ramaraj on 23 Jan 2011
I am. Thanks a ton guys.
Alan Weiss
Alan Weiss on 27 Feb 2012
You might be interested that essentially this same code is in an obscure example in MATLAB documentation:
http://www.mathworks.com/help/toolbox/bioinfo/ug/bs3tbev-1.html#bs3tbev-15

Sign in to comment.

More Answers (6)

Bruno Luong
Bruno Luong on 23 Jan 2011
Careful, mathematically there is always a tiny probability the rejection method (Matt's proposal) returns number of points less than anticipated and it could lead to crash, or run forever if safety looping to fill would be implemented. In practice that might never happens, but if I was a designer of the code for flying an airplane or controlling a robot for brain surgery, then I would never use such code for a peace of mind (mine not the patient).
Here is a direct method
% Data
n = 10000;
radius = rand;
xc = randn;
yc = randn;
% Engine
theta = rand(1,n)*(2*pi);
r = sqrt(rand(1,n))*radius;
x = xc + r.*cos(theta);
y = yc + r.*sin(theta);
% Check
plot(x,y,'.')
% Bruno
Matt Fig
Matt Fig on 23 Jan 2011
Here is a function to do this without skew:
function [X,Y] = rand_circ(N,x,y,r)
% Generates N random points in a circle.
% RAND_CIRC(N) generates N random points in the unit circle at (0,0).
% RAND_CIRC(N,x,y,r) generates N random points in a circle with radius r
% and center at (x,y).
if nargin<2
x = 0;
y = 0;
r = 1;
end
Ns = round(1.28*N + 2.5*sqrt(N) + 100); % 4/pi = 1.2732
X = rand(Ns,1)*(2*r) - r;
Y = rand(Ns,1)*(2*r) - r;
I = find(sqrt(X.^2 + Y.^2)<=r);
X = X(I(1:N)) + x;
Y = Y(I(1:N)) + y;
Now from the command line:
[x,y] = rand_circ(10000,4,5,3);
plot(x,y,'*r')
axis equal
  5 Comments
Show 3 older comments
Amit Goel
Amit Goel on 7 May 2022
why Ns = round(1.28*N + 2.5*sqrt(N) + 100); % 4/pi = 1.2732??
please share the logic behid
Torsten
Torsten on 7 May 2022
Edited: Torsten on 8 May 2022
If you choose 4/pi * N points between -r and r for X and Y, then approximately N points are inside the circle X^2+Y^2 <= r^2 and N*(4/pi-1) points are outside since the area of circle and square are pi*r^2 resp. 4*r^2. The 2.5*sqrt(N) + 100 - term is not clear to me (maybe a confidence level for the uniform distribution).

Sign in to comment.

Bruno Luong
Bruno Luong on 23 Jan 2011
Here is another non-rejection method based on Gaussian RANDN. To goal is to avoid calling sine/cosine. It seems only slightly faster according my test (2011A PreRe)
% Data
n = 1e6;
radius = rand;
xc = randn;
yc = randn;
% Engine 1, rand/cos/sin
tic
theta = rand(1,n)*(2*pi);
r = sqrt(rand(1,n))*radius;
x = xc + r.*cos(theta);
y = yc + r.*sin(theta);
t1=toc;
% Engine 2, randn
tic
x = randn(1,n);
y = randn(1,n);
r2 = x.^2+y.^2;
r = sqrt(rand(1,n)./r2)*radius;
x = xc + r.*x;
y = yc + r.*y;
t2=toc;
fprintf('rand/sin/cos method: %g [s]\n', t1);
fprintf('randn method : %g [s]\n', t2);
% Check
plot(x,y,'.')
% Bruno
Bruno
  1 Comment
Jeff Mason
Jeff Mason on 17 May 2011
How do we know that gaussian x & y give uniform coverage of the circle?

Sign in to comment.

Antti
Antti on 18 May 2011
Here is a multidimensional version of the Brunos code:
%%Modified randn method for hypercircle
% Data
n = 1e6;
n_dim = 3;
pc = randn(1, n_dim);
% Engine 3, multidimensional randn
tic
p = randn(n, n_dim);
r2 = sum(p.^2, 2);
r = sqrt(rand(n, 1) ./ r2)*radius;
p = repmat(pc, n, 1) + repmat(r, 1, n_dim) .* p;
t3 = toc;
fprintf('randn, 3-dim, method : %g [s]\n', t3);
% Check
plot3(p(:,1), p(:,2), p(:,3), '.')
% Antti
  2 Comments
John D'Errico
John D'Errico on 24 Jan 2019
I just saw this answer. It is incorrect, because it creates r incorrectly.
r = sqrt(rand(n, 1) ./ r2)*radius;
Use of sqrt there is a bad idea, since it will not produce uniform sampling in other numbers of dimensions than 2.
Do NOT use the above code in higher dimensions, as it is NOT a uniform sampling.
Bruno Luong
Bruno Luong on 16 Feb 2019
Edited: Bruno Luong on 16 Feb 2019
To generate uniform distribution in spherical shell on higher dimensions
V(a,b) := { X in R^m : a <= |X| <= b }
the correct modification is :
m = 4;
a = 2;
b = 3;
n = 1000; % number of points
s = randn(m,n);
r = (rand(1,n)*(b^m-a^m)+a^m).^(1/m);
c = r./sqrt(sum(s.^2,1));
X = bsxfun(@times, s, c); % (m x n) or simply X = s .* c on MATLAB with auto expansion

Sign in to comment.

James Tursa
James Tursa on 25 Jan 2011
FYI, if you wish to extend methods to higher dimensions, the randn methods are the way to go. Rejection methods in higher dimensions are prone to very small probabilities of being inside the hyper-sphere and thus needing many random numbers generated to get the desired number of samples.
  3 Comments
Show 1 older comment
Aaron Chacon
Aaron Chacon on 16 May 2011
Another problem with rejection methods is that if you reject enough numbers in the pseudo-random number sequence, the next numbers in the sequence can become predictable, and you may not fill the hyper-sphere uniformly.
Walter Roberson
Walter Roberson on 16 May 2011
Is this a serious problem with anything more sophisticated than a linear congruential generator? In particular, with Twister being good up to some 600+ dimensions, is it worth thinking about for Twister?

Sign in to comment.

Greg
Greg on 16 May 2011
Hi guys, I had a bit of a simple question on the above code:
- has anyone tested whether using Sobol/Halton sequences to generate the random numbers improves (if at all possible) the above results?
I've been doing a bit of MC simulations for a project, and some literature suggests that these sequences are better than the 'rand' function (not sure how to measure whether this is the case or not, but graphic display is the only thing that springs to mind). Also unsure what the benefit of improved results v.s. processing time is.
Greg

Categories

Find more on Random Number Generation in Help Center and File Exchange

Tags

Products

Community Treasure Hunt

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

Start Hunting!