How to create a solid spherical cluster with random distribution of points

Hi all, Kindly help me with this.. I need to create a solid (3D) spherical cluster of n random points with 0 as centroid and radius as 2. Points not just on the surface but inside the spherical cluster too. Somewhat like this http://www.mlahanas.de/MOEA/HDRMOGA/Image105.gif

Moreover i want to save all the distances of these random points from centroid in a variable.

Thanks a ton

[Walter Roberson]

What properties should the random distribution have? Uniform in volume? Equal-spaced concentric shells? Equal angle? Expotential Do the points have an associated volume or can they be placed arbitrarily close to each other?

[Matt Kindig]

Vinita, all of Walter's questions are worth considering. My proposed solution does what you want, but will (by the nature of uniform distributions) tend to make points closer to the centroid really close together, which might need be useful to you.

Answer 1

[James Tursa]

Here is my cut at it.

The asin( ) correction is to take care of the fact that we need progressively fewer samples as we get closer to the poles as a function of phi. The relative proportion is governed by the cos function, integrate that to get a sin distribution, and then invert that to get the actual function we need to use on our uniform random samples.

The ^(1/3) is to take care of the fact that as we move out in radius the distribution needs to progressively get more samples. Volume increases by cube of radius, so invert that to get the (1/3) exponent to use on our uniform random samples.

Seems to get the expected results for various spherical sections and spherical caps.

% Uniformly distributed points inside sphere of radius R
% James Tursa
%/
centroid = [0 0 0]; % Center of sphere
R = 2;   % Radius of sphere
n = 10000000;  % Number of points
phi = asin(2*rand(n,1)-1);   %phi between -pi/2 and pi/2, tapering at poles
theta = 2*pi*rand(n,1);   %theta between 0 and 2pi
radius = R*rand(n,1).^(1/3);  %radius between 0 and R, biased outwards for volume
[x,y,z] = sph2cart(theta, phi, radius);
xyz = [ x+centroid(1), y+centroid(2), z+centroid(3)];
%\
% Visualize
%/
plot3(x,y,z,'.');
grid on
%\
% Check expected number inside smaller volumes
%/
V = (4/3)*pi*R^3;
h   = 0.50*R;
h34 = 0.75*R;
Vhalf = (4/3)*pi*h^3;
V34   = (4/3)*pi*h34^3;
r = sqrt(x.*x+y.*y+z.*z);
a = sqrt(R^2 - (R-h)^2);
Vcap = (pi*h/6)*(3*a^2+h^2);
disp(' ');
fprintf('R = %f , h = %f , g = %f\n',R,h,h34);
disp(' ');
fprintf('Expected number of samples inside radius h = %d\n',n*Vhalf/V);
fprintf('Actual   number of samples inside radius h = %d\n',sum(r<h));
disp(' ');
fprintf('Expected number of samples inside radius g = %d\n',n*V34/V);
fprintf('Actual   number of samples inside radius g = %d\n',sum(r<h34));
disp(' ');
fprintf('Expected number of samples inside spherical cap h     = %d\n',round(n*Vcap/V));
fprintf('Actual   number of samples inside spherical cap x > h = %d\n',sum(x> h));
fprintf('Actual   number of samples inside spherical cap x <-h = %d\n',sum(x<-h));
fprintf('Actual   number of samples inside spherical cap y > h = %d\n',sum(y> h));
fprintf('Actual   number of samples inside spherical cap y <-h = %d\n',sum(y<-h));
fprintf('Actual   number of samples inside spherical cap z > h = %d\n',sum(z> h));
fprintf('Actual   number of samples inside spherical cap z <-h = %d\n',sum(z<-h));
disp(' ');

[Richard Brown]

I can confirm those are correct - you can derive the relationships pretty easily by considering the volume element

dV = r^2 cos(phi) dr dtheta dphi

e.g. if we want the expression for phi in terms of a uniform random variable u, then we require

cos(phi) dphi = a du

for some constant a. Then, integrating gives

sin(phi) = au + b

To find the constants a and b we know phi(0) = -pi/2 => b = -1, and phi(1) = pi/2 => a = 2. So

phi(u) = asin(2u - 1)

etc.

[Vinita]

This is great.

But now the only problem is that I need only the magnitude of the distances of the points from the centroid(only the magnitude of polar coordinates).Look at this figure :-

http://1.bp.blogspot.com/_wrPEFjzD-9o/TTiqlLzm_fI/AAAAAAAAFe8/rY3WVP39nFo/s400/graph+paper+polar+coordinates.PNG

So I need only thr r (magnitude/distance part) from the centroid for each point.

What we have right now is x,z & z coordinates of each point.

Thats it

[James Tursa]

Ummmm ... the r you are talking about sounds like the r that is already in the code. Am I missing something?

Answer 2

[Matt Kindig]

centroid = [0 0 0];
Rad = 2;   %desired radius
n = 100;  %or however many points you want
phi = 2*pi*rand(n,1);   %phi between 0 and 2pi
theta = pi*rand(n,1);   %theta between 0 and pi
radius = Rad*rand(n,1);  %radius between 0 and Rad
[x,y,z] = sph2cart(theta, phi, radius);
xyz = [ x+centroid(1), y+centroid(2), z+centroid(3)];
%'radius' contains the distances to the centroid

[Ashish Uthama]

and then use:

plot3(x,y,z,'k*'); axis equal; rotate3d('on')

to visualize.

[John D'Errico]

These points as generated will NOT be uniformly spaced in the sphere.

[Vinita]

Thanks a ton Matt Kindig !!!!

Could you pls tell me if these points are randomly distributed or not.

Because thats the main concern

[Richard Brown]

They'll be clustered about the centre and also more dense towards the north and south pole

Answer 3

[Richard Brown]

Here's the lazy (no thinking, no maths) way. It will be slow(ish), but serves to illustrate the idea. I'm sure you can think up optimisations if you require them.

n = 1000;
i = 0;
X = zeros(3, n);
while i < n
  x = 4*rand(3, 1) - 2;
  if norm(x) <= 2
    i = i + 1;
    X(:, i) = x;
  end
end

plot3(X(1,:), X(2,:), X(3,:), '.'), axis equal