Generate random points inside a circle

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.

[Mahesh Ramaraj]

Thanks Guys, these codes helped me a lot.

Answer 1

[Paulo Silva]

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

[Walter Roberson]

sqrt(rand)is going to concentrate the points non-linearly away from the center of the circle. That might be acceptable, but if the points must be uniformly distributed as if in a square that has had its edges cut off, then a different generation technique would have to be used.

[Paulo Silva]

The question didn't mentioned other conditions beside the points inside the circle, thank you for pointing that out.

[Matt Fig]

I believe Paulo's code does distribute the points randomly as if a unit square had its corners rounded off, Walter. sqrt(2) is the correct transformation for this problem. My code below uses the 'cut the edges' technique and can be used for statistical comparison with Paulo's code, but I think they give statistically equal results.

[Bruno Luong]

I just reread the comments. Paulo's code is exactly like my solution and the SQUAREROOT is correctly used to transform the distribution so as the points are uniform on the disc.

[Paulo Silva]

Mahesh Ramaraj must be happy with all the sample code provided by us :)

[Mahesh Ramaraj]

I am. Thanks a ton guys.

[Alan Weiss]

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

Answer 2

[Bruno Luong]

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

[Mahesh Ramaraj]

I like this...very simple and elegant...

Answer 3

[Matt Fig]

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

[Paulo Silva]

Nice and fast code, can you please explain more about what's that "skew" and/or a reference for more info about it?

[Matt Fig]

Without skew means that there is no preferred direction or area on the disk. Your code does this too, Paulo.

[Paulo Silva]

Thanks Matt, I now know what skew means :)

[Amit Goel]

why Ns = round(1.28*N + 2.5*sqrt(N) + 100); % 4/pi = 1.2732??

please share the logic behid

[Torsten]

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).

Answer 4

[Bruno Luong]

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

[Jeff Mason]

How do we know that gaussian x & y give uniform coverage of the circle?

Answer 5

[Antti]

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

[John D'Errico]

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]

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

Answer 6

[James Tursa]

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.

[Matt Fig]

Ah yes, the curse of dimensionality. Very true. Of course, the rejection method is used all the time for lowly 3-D problems in the industry where I work, despite Bruno's dire warnings.

[Aaron Chacon]

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]

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?

Answer 7

[Greg]

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