Thread Subject:

Random points in a circular zone

Hi, I want to create a function that give me N random points in a circular zone.

The circular zone can be represented as
A={(x,y)=(R*cos(theta), R*t*sin(theta)), t in [-1, 1], theta in [a,b]}
where 0<a<b<pi.

I have written a function for the circular segment, that is when a=0, but this function doesn't give me a uniform distribution and I don't understand how apply it to the zone.

function [x,y]= RandSeg(N, R, w)
    angle = -w+2*w*rand(N,1); % theta is in [0, w]
    k=linspace(-1,1,N)';
    x=R.*cos(angle);
    y=R.*k.*sin(angle);

Can anyone help me?
Thank you!
Valentina

"Valentina" wrote in message <jpi7su$g55$1@newscl01ah.mathworks.com>...
> Hi, I want to create a function that give me N random points in a circular zone.
>
> The circular zone can be represented as
> A={(x,y)=(R*cos(theta), R*t*sin(theta)), t in [-1, 1], theta in [a,b]}
> where 0<a<b<pi.
>
> I have written a function for the circular segment, that is when a=0, but this function doesn't give me a uniform distribution and I don't understand how apply it to the zone.
>
> function [x,y]= RandSeg(N, R, w)
> angle = -w+2*w*rand(N,1); % theta is in [0, w]
> k=linspace(-1,1,N)';
> x=R.*cos(angle);
> y=R.*k.*sin(angle);
>
> Can anyone help me?
> Thank you!
> Valentina

How about:

d = 2*R;
r2 = R*R;
x = zeros(1,N);
y = x;
count = 1;
while count <= N
  xp = d*(rand-0.5);
  yp = d*(rand-0.5);
  if (xp*xp + yp*yp) < r2
    x(count) = xp;
    y(count) = yp;
    count = count + 1;
  end
end

Not sure it can be vectorised.

"Valentina" wrote in message <jpi7su$g55$1@newscl01ah.mathworks.com>...
> Hi, I want to create a function that give me N random points in a circular zone.
>
> The circular zone can be represented as
> A={(x,y)=(R*cos(theta), R*t*sin(theta)), t in [-1, 1], theta in [a,b]}
> where 0<a<b<pi.
>
> I have written a function for the circular segment, that is when a=0, but this function doesn't give me a uniform distribution and I don't understand how apply it to the zone.
>
> function [x,y]= RandSeg(N, R, w)
> angle = -w+2*w*rand(N,1); % theta is in [0, w]
> k=linspace(-1,1,N)';
> x=R.*cos(angle);
> y=R.*k.*sin(angle);
>
> Can anyone help me?
> Thank you!
> Valentina

Actually, here's a better vectorised method. The square root provides the uniform distribution.

function [x,y]= RandSeg(N, R, w)
  angle = -w+2*w*rand(N,1);
  rad = rand(N,1).^0.5;
  x=R*rad.*cos(angle);
  y=R*rad.*sin(angle);

I though similar question was answered in a previous thread:

R = 7;
a = pi/3;
b = pi;

n=1e4;
w = b-a;
r=R*sqrt(rand(1,n));
theta = a + rand(1,n)*w;
x=r.*cos(theta);
y=r.*sin(theta);

plot(x,y,'.')
axis equal

% Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <jpivel$ra5$1@newscl01ah.mathworks.com>...
> I though similar question was answered in a previous thread:
>
> R = 7;
> a = pi/3;
> b = pi;
>
> n=1e4;
> w = b-a;
> r=R*sqrt(rand(1,n));
> theta = a + rand(1,n)*w;
> x=r.*cos(theta);
> y=r.*sin(theta);
>
> plot(x,y,'.')
> axis equal
>
> % Bruno

Hi Bruno, your function generates n random points in a circular sector, but unfortunately I want a function that generates them in a circular zone and I don't understand how to change it.
Valentina

"Valentina" wrote in message <jpj32t$el4$1@newscl01ah.mathworks.com>...
> Hi Bruno, your function generates n random points in a circular sector, but unfortunately I want a function that generates them in a circular zone and I don't understand how to change it.

How do you define a circular sector (I do not understand you first post)? In the mean time try this

% Radius min and max
Rmin = 2;
Rmax = 5;
% Angle min and max
a = pi/3;
b = pi;

n=1e4;
w = b-a;
r = sqrt(rand(1,n)*(Rmax^2-Rmin^2)+Rmin^2);
theta = a + rand(1,n)*w;
x = r.*cos(theta);
y = r.*sin(theta);

plot(x,y,'.')
axis equal

% Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <jpj4jl$lb3$1@newscl01ah.mathworks.com>...
> How do you define a circular sector (I do not understand you first post)? In the mean time try this
>
> % Radius min and max
> Rmin = 2;
> Rmax = 5;
> % Angle min and max
> a = pi/3;
> b = pi;
>
> n=1e4;
> w = b-a;
> r = sqrt(rand(1,n)*(Rmax^2-Rmin^2)+Rmin^2);
> theta = a + rand(1,n)*w;
> x = r.*cos(theta);
> y = r.*sin(theta);
>
> plot(x,y,'.')
> axis equal
>
> % Bruno

Your last function gives random points in an anular sector, but I want a circular zone. A circular zone is the portion of a disk included between any two parallel chords and their intercepted arcs and I found that it can be represented as I write in my first post.
Valentina

I see.

I have the feeling that there is no close form formula that can generate unform distribution on a sector.

Bruno

On 5/23/2012 1:39 PM, Bruno Luong wrote:
> I see.
>
> I have the feeling that there is no close form formula that can
> generate unform distribution on a sector.
>
> Bruno

An inefficient but theoretically correct scheme is known as rejection.
You generate a bunch of points that are uniformly distributed on a
region that is larger than the region you have in mind, and then throw
away all points that are outside the region you have in mind.

So, for example, you could use Bruno's method of generating points in a
disk, and then remove the points that are not in your zone. The result
will have points uniformly distributed within your zone.

The problem, of course, is you are never sure beforehand how many points
to generate so that you keep the number you are aiming for. In other
words, when you use rejection, you have to generate more points than you
are aiming for, and you don't know how many more.

Alan Weiss
MATLAB mathematical toolbox documentation

FYI: Original thread: http://www.mathworks.com/matlabcentral/newsreader/view_thread/319832#876000

"Valentina" wrote in message <jpi7su$g55$1@newscl01ah.mathworks.com>...
> Hi, I want to create a function that give me N random points in a circular zone.
> The circular zone can be represented as
> A={(x,y)=(R*cos(theta), R*t*sin(theta)), t in [-1, 1], theta in [a,b]}
> where 0<a<b<pi.
- - - - - - - - - -
  The following will give a statistically uniform distribution area-wise in your circular zone without rejection. Let R, a, and b be as defined in your original post. (I have called your t parameter s and your theta parameter t.)

 n = 8000; % Select how many random points you want to generate

 p = (2*a-sin(2*a))+(2*b-sin(2*b)-2*a+sin(2*a))*rand(n,1);
 p = mod(p+pi,2*pi)-pi;
 p = sign(p).*abs(p).^(1/3);
 c1 = 1/(6^(1/3)); % c1 and c2 used in approx. of derivative
 c2 = (c1-2/3/pi^(2/3))/pi^2;
 t = zeros(n,1);
 for k = 1:12 % Twelve trips should suffice for convergence
  p1 = t-sin(t);
  p1 = sign(p1).*abs(p1).^(1/3);
  t = t - (p1-p)./(c1-c2*t.^2); % Pseudo Newton Raphson
 end
 t = mod(t,2*pi)/2;
 s = -1+2*rand(n,1);
 x = R*cos(t); % x and y are cartesian coordinates of the random points
 y = R*s.*sin(t);

Roger Stafford

> The following will give a statistically uniform distribution area-wise in your circular zone without rejection. Let R, a, and b be as defined in your original post. (I have called your t parameter s and your theta parameter t.)
>
> n = 8000; % Select how many random points you want to generate
>
> p = (2*a-sin(2*a))+(2*b-sin(2*b)-2*a+sin(2*a))*rand(n,1);
> p = mod(p+pi,2*pi)-pi;
> p = sign(p).*abs(p).^(1/3);
> c1 = 1/(6^(1/3)); % c1 and c2 used in approx. of derivative
> c2 = (c1-2/3/pi^(2/3))/pi^2;
> t = zeros(n,1);
> for k = 1:12 % Twelve trips should suffice for convergence
> p1 = t-sin(t);
> p1 = sign(p1).*abs(p1).^(1/3);
> t = t - (p1-p)./(c1-c2*t.^2); % Pseudo Newton Raphson
> end
> t = mod(t,2*pi)/2;
> s = -1+2*rand(n,1);
> x = R*cos(t); % x and y are cartesian coordinates of the random points
> y = R*s.*sin(t);
>
> Roger Stafford

Thank you so much Roger, it's exactly what I wanted.
Valentina

n = 8000;
a = pi/6;
b = pi/2+pi/6;
R = 7;
I think rejection method could be the fastest in average.

xmin = R*cos(b);
xmax = R*cos(a);
ymin = -R;
ymax = R;

xy = zeros(0,2);
while true
    M = 2*n - size(xy,1);
    x = xmin + (xmax-xmin)*rand(M,1);
    y = ymin + (ymax-ymin)*rand(M,1);
    in = x.^2 + y.^2 <= R^2;
    xy = cat(1, xy, [x(in) y(in)]);
    if size(xy,1) >= n
        xy(n+1:end,:) = [];
        break
    end
end

plot(xy(:,1),xy(:,2),'.');
axis equal

% Bruno