randomly generated spheres in a volume in Matlab

Hey,

I want to generate n spheres of defined radius with in a defined volume at random locations. Following the condition of not intersecting each other.

Can someone help me with this.

Regards

Om

[darova]

Do you have any ideas?

Answer 1

[Jan]

function P = GetRandomSpheres(nWant, Width, Radius)
% INPUT:
%   nWant:  Number of spheres
%   Width:  Dimension of 3d box as [1 x 3] double vector
%   Radius: Radius of spheres
% OUTPUT:
%   P:      [nWant x 3] matrix, centers

P = zeros(nWant, 3);

R2     = (2 * Radius) ^ 2;   % Squared once instead of SQRT each time
W      = Width - 2 * Radius; % Avoid interesction with borders
iLoop  = 1;                  % Security break to avoid infinite loop
nValid = 0;
while nValid < nWant && iLoop < 1e6
  newP = rand(1, 3) .* W + Radius;
  % Auto-expanding, need Matlab >= R2016b. For earlier versions:
  % Dist2 = sum(bsxfun(@minus, P(1:nValid, :), newP) .^ 2, 2);
  Dist2 = sum((P(1:nValid, :) - newP) .^ 2, 2);
  if all(Dist2 > R2)
    % Success: The new point does not touch existing sheres:
    nValid       = nValid + 1;  % Append this point
    P(nValid, :) = newP;
  end
  iLoop = iLoop + 1;
end
% Stop if too few values have been found:
if nValid < nWant
  error('Cannot find wanted number of points in %d iterations.', iLoop)
end
end

And a demo code:

n = 10;
R = 10;
P = GetRandomSpheres(n, [100, 100, 100], R);
figure
axes('NextPlot', 'add', ...
    'XLim', [0, 100], 'YLim', [0, 100], 'ZLim', [0, 100]);
view(3);
[X, Y, Z] = sphere();
for k = 1:n
    surf(X * R + P(k, 1), Y * R + P(k, 2), Z * R + P(k, 3));
end

[Bahareh Vossoughi]

Hi,

I've been trying to create spheres with different radius randomly in a volume without intersecting. How can I rewrite your code?

Here is what I've written but the spheres intersect with each other. I would appreciate your help.

Regards

BV

function [S1, S2] = RandomSpheres (Num1, Num2, RVE, R1, R2)

S1 = zeros (Num1, 3);
S2 = zeros (Num2, 3);
A1 = RVE - 2*R1;  % No intersections with the borders
A2 = RVE - 2*R2;
iloop = 1;
NValid = 0;
while NValid < Num1 && NValid < Num2 && iloop < 1e6
    newS1 = rand(1,3).*A1 + R1; %The centre of spheres
    newS2 = rand(1,3).*A2 + R2;
    Dist1 = sum ((S1(1:NValid, :)- newS1).^2, 2);
    Dist2 = sum ((S2(1:NValid, :)- newS2).^2, 2);
    Dist3 = sum ((S1(1:NValid, :)- newS2).^2, 2);
      if all(Dist1 > (R2+R1)^2) && all(Dist2 > (R2+R1)^2) && all(Dist3 > (R2+R1)^2) && all(Dist1 > (2*R2)^2) && all(Dist2 > (2*R2)^2) && all(Dist3 > (2*R2)^2) && all(Dist1 > (2*R1)^2) && all(Dist2 > (2*R1)^2) && all(Dist3 > (2*R1)^2)                      
           NValid = NValid + 1;
           S1(NValid, :) = newS1;
           S2(NValid, :) = newS2;
      end
          iloop = iloop + 1;
end
end

[Om Jagtap]

You can just do this and use the function shared by Jan in previous answers.

n = 100;
R = 3.5;
P = GetRandomSpheres(n, [100, 100, 100], R);
figure
axes('NextPlot', 'add', ...
    'XLim', [0, 100], 'YLim', [0, 100], 'ZLim', [0, 100]);
view(3);
[X, Y, Z] = sphere();
for k = 1:n
    o=rand(1)
    surf(X * o*R + P(k, 1), Y * o*R + P(k, 2), Z *o* R + P(k, 3));
end

[Bahareh Vossoughi]

I don't want random radius actually. I have two specific radius, and I want to spread multiple numbers of these 2 spheres, randomly in the cube without intersecting.

[Om Jagtap]

Hey,

I am not a very good coder but hope this can help you!

n = 100;
Ra = 3.5;
Rb = 5;
if Ra>Rb
    P = GetRandomSpheres(n, [100, 100, 100], Ra);
elseif Rb>Ra
    P = GetRandomSpheres(n, [100, 100, 100], Rb)
end

figure
axes('NextPlot', 'add', ...
    'XLim', [0, 100], 'YLim', [0, 100], 'ZLim', [0, 100]);
view(3);
[X, Y, Z] = sphere();
for k = 1:n
    o=randi([1,2]);
    if o==1
        surf(X * Ra + P(k, 1), Y * Ra + P(k, 2), Z * Ra + P(k, 3));
    elseif o==2
        surf(X *Rb + P(k, 1), Y * Rb + P(k, 2), Z *Rb + P(k, 3));
    end 
end

[Bahareh Vossoughi]

Hey,

Thanks for your time:) I tried two different while loops and it worked! But I have to add the if loop as the one you've written to compare diameters first.

BV

[Jan]

Now a code, which let you choose a specific radius for each sphere:

function P = GetRandomSpheres(nWant, Width, Radius)
% INPUT:
%   nWant:  Number of spheres
%   Width:  Dimension of 3d box as [1 x 3] double vector
%   Radius: [nWant x 1] vector, radii of spheres
% OUTPUT:
%   P:      [nWant x 3] matrix, centers
P      = zeros(nWant, 3);
R2     = (2 * Radius(:)) .^ 2;  % Squared to avoid expensive SQRT
iLoop  = 1;                     % Security break to avoid infinite loop
index  = 1;
while index <= nWant && iLoop < 1e6
  newP = rand(1, 3) .* (Width - 2 * Radius(index)) + Radius(index);
  % Auto-expanding, need Matlab >= R2016b. For earlier versions:
  % Dist2 = sum(bsxfun(@minus, P(1:nValid, :), newP) .^ 2, 2);
  Dist2 = sum((P(1:index-1, :) - newP) .^ 2, 2);
  if all(Dist2 > R2(1:index-1) + R2(index))
    % Success: The new point does not touch existing sheres:
    P(index, :) = newP;
    index       = index + 1; 
  end
  iLoop = iLoop + 1;
end

% Error if too few values have been found:
if index < nWant
  error('Cannot find wanted number of points in %d iterations.', iLoop)
end
end

And almost the same demo code:

n = 30;
R = randi([2,10], n, 1);  % Or however your radii are defined
% Perhaps:
% R = repelem([10, 20], n / 2);  % 15 with R=10 and 15 with R=20
P = GetRandomSpheres(n, [100, 100, 100], R);
figure
axes('NextPlot', 'add', ...
    'XLim', [0, 100], 'YLim', [0, 100], 'ZLim', [0, 100]);
view(3);
[X, Y, Z] = sphere();
for k = 1:n
    surf(X * R(k) + P(k, 1), Y * R(k) + P(k, 2), Z * R(k) + P(k, 3));
    %         ^^^                 ^^^                 ^^^
    % Bugfix based on Jared Cole's comment. Thanks!
end

[Jared Cole]

There is a minor error in the plotting part of the demo, you need to select each radii when plotting.

surf(X * R(k) + P(k, 1), Y * R(k) + P(k, 2), Z * R(k) + P(k, 3));

[Jan]

@Jared Cole: Thanks! Bug fixed.

[PB]

Dear @Jan,

I want to use your code for similar problem, however, the differences are as follows:

1. I want to generate spheres with radius from R1 to R2 in the specified volume. It can be a cylinder or a cube (in the code, which I have passed below it is cube, but I do not know how to change it to cylinder).
2. The spheres inside the specified volume should have a specified volume fraction, i.e. 0.4 of V.
3. I want to have csv with the center location and the radius of each sphere (it can be obtain from P and R matrixes).

The problem is that I have used your modified code, however, changing the volume fraction does nothing. It is always a similar number spheres despite a different volume fraction.

Code:

function P = GetRandomSpheres2(Fract, Width, Radius)
% INPUT:
%   nWant:  Number of spheres
%   Width:  Dimension of 3d box as [1 x 3] double vector
%   Radius: [nWant x 1] vector, radii of spheres
% OUTPUT:
%   P:      [nWant x 3] matrix, centers
P      = zeros(3000, 3);
R2     = (2 * Radius(:)) .^ 2;  % Squared to avoid expensive SQRT
iLoop  = 1;                     % Security break to avoid infinite loop
index  = 1;
V = 0;
Vmax = Width(1)*Width(2)*Width(3);
Vfract = Fract * Vmax;
nWant = Vfract / (4/3*pi*min(Radius)^3);
while V <= Vfract && iLoop < 1e6
  newP = rand(1, 3) .* (Width - 2 * Radius(index)) + Radius(index);
  % Auto-expanding, need Matlab >= R2016b. For earlier versions:
  % Dist2 = sum(bsxfun(@minus, P(1:nValid, :), newP) .^ 2, 2);
  Dist2 = sum((P(1:index-1, :) - newP) .^ 2, 2);
  if all(Dist2 > R2(1:index-1) + R2(index))
    % Success: The new point does not touch existing sheres:
    P(index, :) = newP;
    V = V + (4/3*pi*Radius(index)^3);
    index = index + 1; 
  end
  iLoop = iLoop + 1;
end
end

Demo file:
Width = [50, 50, 50];
Fract = 0.7;
Radius = [2,4];
%
R = randi(Radius, nWant, 1);  % Or however your radii are defined
P = GetRandomSpheres2(Fract, Width, R);
figure
axes('NextPlot', 'add', ...
    'XLim', [0, 100], 'YLim', [0, 100], 'ZLim', [0, 100]);
view(3);
[X, Y, Z] = sphere();
for k = 1:nWant
    surf(X * R(k) + P(k, 1), Y * R(k) + P(k, 2), Z * R(k) + P(k, 3));
end

Can you help me with this?

Bet regards

[Ismail Can Karaoglu]

Hello @Jan,

I am currently dealing with the packing of equal-sized spheres into another sphere without overlapping. I am trying to adapt your code into my system but I have to create a sphere instead of the box. When I did that, I basically faced with a matrix dimension error because of the fact that MATLAB creates spheres with faces that come out n number of coordinates for each dimension. How can I overcome this situation? If you have any suggestion, I will appreciate it.

Thanks in advance.

[Partha Chakraborty]

Hi @Jan,

I am currently working on randomly packed spheres of equal sizes in a cylindrical beaker. Can you please help me with the code modification?

Thanks,

Partha

[Jan]

@Partha Chakraborty: What does "randomly packed" mean? Randomly means, that the package is not dense, but as far as I understand "packed" means dense.

The problem of a dense packing is very complicated and you can publish books about it. So do you have a certain approach or goal? I think opening a new question is the best approach, because setting some spheres randomly in a box is a different problem than producing a dense package.

[Partha Chakraborty]

Hi @Jan,

Sorry if I could not clarify earlier. It is more suited to randomly distributed. So, I want to model random distribution of n number of spheres in a cylindrical beaker. That's it.

[Jan]

@Partha Chakraborty: In the code of my answer, you can insert a test after the line:

newP = rand(1, 3) .* W + Radius;

Now reject a point, if it is outside the cylinder.

There are cheaper methods that a rejection to created normally distributed points inside a cylinder, but creating new random points is cheap.

[Partha Chakraborty]

@Jan how can I define the dimension of the cylinder ?

[Walter Roberson]

If you center the cylinder on (0,0) and it is radius R and height H, then

newP = (rand(1, 3) * 2 - 1) .* [R, R, H/2];
if newP(1)^2 + newP(2)^2 > R^2; continue; end  %reject if outside the cylinder
x = newP(1); y = newP(2); z = newP(3);

Or you could code

radius = sqrt(rand) * R;
angle = rand() * 2*pi;
vertical = rand() * H;
[x, y, z] = pol2cart(angle, radius, vertical);
z = z - H/2;   %to center around 0

[Walter Roberson]

Testing to ensure the sqrt() is appropriate:

R = 4; H = 3;
N = 1000;
r2 = rand(1,N);
angle = rand(1,N) * 2*pi;
vertical = rand(1,N) * H;

radius = sqrt(r2) * R;

[x, y, z] = pol2cart(angle, radius, vertical);
z = z - H/2;   %to center around 0

scatter3(x, y, z)

Compared to the same random data without the sqrt()

radiusB = r2 * R;

[xB, yB, zB] = pol2cart(angle, radiusB, vertical);
zB = zB - H/2;   %to center around 0

scatter3(xB, yB, zB)

You can see that the version without the sqrt() is denser in the core instead of being equally distributed through space.

[Ju Lia]

Hey!

I used your code with just a few adjustments and it works very well for me and helped me a lot. Now I wanted to see the distribution in a cut through the 3D Plot. slice or contourslice are not working because of the missing volumeinformation. Do you have an idea for a workaround?

Thanks!

[Simple Life]

Hey!

I used your code with some cases and it worked well.

However when particle number is huge since volume fraction is still less than 1, it failed.

Do your have any idea to deal with it?

Thanks for any help.

[Walter Roberson]

I am not sure who you are addressing? The post with volume fraction was

https://www.mathworks.com/matlabcentral/answers/461635-randomly-generated-spheres-in-a-volume-in-matlab#comment_798412

and that code gives up after 1e6 tries.

[Image Analyst]

@Simple Life obviously for a high enough number there will be no room left to place any more spheres of a certain radius or larger. The only way to fit in more spheres and increase the volume fraction of occupied space is to make the radius smaller.

[Jan]

@Simple Life: Yes, such rejection methods tends to fail for denser packages.

Which problem do you want to solve? How dense and how random should your result be?

[Simple Life]

Code used by me is:

clc;clear all;close all;

% step 1. input all parameter
BoxSize = [0.4, 0.8, 0.4];
grianDiameter =  0.028 ; % particle size in mm, diameter, 0-1

wRatioGrain = 0.4;   %weight ratio of grain, abrasive = grain + media

dGrain = 3.2 ;  %density of grain
dMedia = 0.9 ;  %density of media


% step 3. cal total number of grain
vRatioGrain = dMedia / ( dMedia + (1-wRatioGrain) / wRatioGrain *  dGrain );    % to get volumn ratio of grain
vBox =  prod(BoxSize);   % to get total volumn of abrasive
vGrain =  vRatioGrain * vBox;   % to get total volumn of abrasive
nGrain = vGrain / (pi * grianDiameter.^3/6);    % to get total numberof grain
nGrain = round(nGrain);    % to make nGrain is not a decimal

disp(['diameter of particle is(mm): ', num2str(nGrain, '%2d')])
disp(['weight fraction of particle is(wt%):',num2str(wRatioGrain*100,'%.2f')])
disp(['box size is(mm): ',num2str(BoxSize,'%.0d  ')])   
disp([ 'volumn fraction of particle is(vt%): ',num2str(vRatioGrain*100,'%.2f'),'%, '])
disp(['Grain-No:', num2str(nGrain, '%2d')])


%-----------------------------------------------------------------
R = grianDiameter/2;

P = GetRandomSpheres(nGrain, BoxSize, R);

figure
axes('NextPlot', 'add', ...
    'XLim', [0, BoxSize(1)], 'YLim', [0, BoxSize(2)], 'ZLim', [0, BoxSize(3)]);
view(3);
[X, Y, Z] = sphere();
for k = 1:nGrain
    surf(X * R + P(k, 1), Y * R + P(k, 2), Z * R + P(k, 3));
end
box on
ax = gca;
ax.BoxStyle = 'full';
axis equal

when wRatioGrain is higher, say 0.7, the function you provided won't work propertly. Still, you gave us a fabulous function.

Thanks for any help. May you have a good day.

[Jan]

This is a statistical effect: For a random distribution with the density of > 0.7 there is a small chance only to find a remaining free space of the wanted size.

[Walter Roberson]

Wikipedia: "For equal spheres in three dimensions, the densest packing uses approximately 74% of the volume. A random packing of equal spheres generally has a density around 63.5%"

The code was never intended to be a dense-packing algorithm, which requires much different approaches.

[Simple Life]

@Walter Roberson Right answer. Thanks a lot. @Jan and @Walter Roberson

wiki reference : Sphere_packing

[Russell Jones]

Hello @PB,

I am looking into a similar problem where my volume fraction is 0.59, and from testing @Jan's modified code, it unfortunately seems to only work for volume fractions of 0.31 or lower. I have increased the iterations to 1e9 and I still can't get above 0.31. From what I can tell, it would require a different approach. Did you have any success with a different method / does anyone know what method I should attempt for a packing ratio of 0.59?

[Walter Roberson]

@Russell Jones

@Jan wrote code for the requirement that the spheres be non-intersecting. To get to the volume fractions you are talking about, you would need intersecting (touching) spheres, and you will need to use some sort of compaction algorithm. As the theoretical limit is only 63.4% you will need to do a fair amount of compaction to get 0.59. Or you could deliberately generate one of the regular packings.

https://en.wikipedia.org/wiki/Sphere_packing#Irregular_packing

Answer 2

[Mehdi Chouchane]

Here is a code I'm using which is a custom version of a code taken from Matlab Community.

It generates the desired volume fraction of non-overlapping spheres with radii varying between rmin and rmax.

function [centers, rads] = sampleSpheres( Fract,dims,rmin,rmax,verbosity)
% main function which is to be called for adding multiple spheres in a cube
% dims is assumed to be a row vector of size 1-by-ndim
% For demo take dims = [ 2 2 2 ] and n = 50

% preallocate
V = Fract * dims(1) * dims(2) * dims(3) ; 
v = 0 ; 
ndim = numel(dims);
n = round(1.5*V / ((4*pi*((rmax+rmin)/2)^3)/3)) ;
centers = zeros( n, ndim );
rads = zeros( n, 1 );
ii = 1;
n = 0 ; 
while v < V
    [centers(ii,:), rads(ii) ] = randomSphere2( dims,rmin,rmax);
    if nonOver2( centers(1:ii,:), rads(1:ii))
        n = n + 1 ; 
        v = v + (4*pi*rads(ii)^3)/(3) ;
        ii = ii + 1; % accept and move on
        if verbosity == 1
            100*v/V
        end
     end
end
centers = centers(~all(centers == 0, 2),:);
rads = rads(~all(rads == 0, 2),:);
end

function [ c, r ] = randomSphere2( dims,rmin,rmax)
% creating one sphere at random inside [0..dims(1)]x[0..dims(2)]x...
% In general cube dimension would be 10mm*10mm*10mm
% radius and center coordinates are sampled from a uniform distribution
% over the relevant domain.
% output: c - center of sphere (vector cx, cy,... )
%         r - radius of sphere (scalar)
r = rmin + ( rmax - rmin) .* rand(1);
c = bsxfun(@times,(dims - r) , rand(1,3)) + r; % to make sure sphere is placed inside the cube
end

function not_ovlp = nonOver2( centers, rads) % to make sure spheres do not overlap
if numel( rads ) == 1
    not_ovlp = true;
    return; % nothing to check for a single or the first sphere
end
center_dist = sqrt(sum(bsxfun(@minus,centers(1:end-1,:),centers(end,:)).^2,2));
radsum = rads(end) + rads(1:end-1);
not_ovlp = all(center_dist >= radsum);
return;
end

[PB]

Dear Mehdi Chouchane,

Thank you!! It's working like a charm:)

Allll theee beeest regards,

Paweł

[Sonti Srinivasa Murali Kartheek]

Can you please share the demo file to run the function?

[Bhishm Dewangan]

Please share demo file

[Walter Roberson]

sampleSpheres( .1, [2 2 2 2], 0.3, 0.5, 0)

[Russell Jones]

Hello @Mehdi Chouchane,

I tried using your code and it worked better to achieve a higher volumetric fraction. However, I am unable to get the code to work for volumetric fractions greater than 0.5. As I need 0.59, I am struggling to figure out how to add to the code to get to this fraction. Do you have any tips?

[Walter Roberson]

@Russell Jones

As I already discussed with you, this kind of code is for creating non-overlapping spheres that might have any amount of distance between them (though the larger the gap the more likely that another sphere will be randomly placed in the gap volume.)

In order to get the kinds of packings you need, you need code that actively compacts.

For example you might want to see if you can find the discussion from 5-ish years ago from the person who needed to model accreation, which was handled by "dropping" particles from the top and letting them "fall" until they were in a stable situation. That kind of code might still not get you close enough for your purposes, but it would be better than the random-placement-in-a-volume approaches.

Answer 3

[sadedbouta]

Hey, are you help me for this lines in script matlab with Comsol:

n = 100;
Vsum = 0;
% Generate position vectors
Pos = zeros(n,3) ; % XYZ
R = zeros(n,3);
idx = 1; % index for sphere 
flag = 0;
while (Vsum < Vsq * vf)
    r = abs( normrnd(miu,sigma) ); 
    % generates a random number from the normal distribution with mean parameter mu and 
    % standard deviation parameter sigma.
      
    pos = [blk_size * rand(1,1) blk_size * rand(1,1) blk_size * rand(1,1)];
    
       
    for k = 1:idx %Check the distance between the randomly generated sphere and all existing spheres.
        Distance = sqrt((pos(1)-Pos(k,1))^2+(pos(2)-Pos(k,2))^2+(pos(3)-Pos(k,3))^2);
        
        rsum = r   
        %rsum = r+R(k);
        if Distance < 2*rsum     
            flag = 1;
            break;
        end
    end

Can you idea me corrected the error in this script?

[Robert U]

When "idx" exceeds "n" the matrix dimension of Pos is not sufficient to be addressed. Your variable "n" does not have any influence on the max. number of while-loop iterations.

See here as well.

Kind regards,

Robert

Answer 4

[Serigne Saliou Sene]

@JanHello Jan I hope you are well, I would like you to improve my script I have already randomly generated several ellipsoides in the cube I want a plot (volumeCube +aggregate+fiber) Thanks in advance

clear ; clc; close all;

%input
%------------------
H=150;                   %Cube Height (mm)
W=150;                   %Cube Width  (mm)
L=150;                   %Cube Length (mm) 

x=[0 L L 0];                                                               %Specimen dimensions 
y=[0 W];                                                                   %Specimen dimensions
z=[0 0 H H];                                                               %Specimen dimensions

Classes_diameters=[19 12.70 9.5 4.75 2.36];           %Particles classes diameters  (descendingly)
Alpha=0.45;                                         %Fuller's curve exponent
m=3;                                                %Particles shape distribution factor
Particle_ratio=0.50;                                %Particles ratio of the volume
er=0.05;                                            %Spacing factor between particles

dist=W/2;                                           %Cutting distance for ellipses
r_min=2.36;                                            %Minimum ellipse radius involved   

L=3;                                                 %Length of fibers
N=1200;                                               %Number of fibers
DFiber=3;                                             %Diameter of fibers
Orientation=[];                                       %Orientation: Orientation of fibers as [l; m; n] column vector where l, m, and n are direction cosines of orientation in x, y, and z directions, respectively. 
Ndiv=1;                                               %Number of fiber mesh divisions

%------------------------

Classes=Particles_Generation(x,y,z,Classes_diameters,Alpha,m,Particle_ratio);

Plot_Sieve(Classes,x,y,z,Classes_diameters,Alpha,Particle_ratio);

Ellipsoids=Particles_Distribution(Classes,x,y,z,er);

Plot_Ellipsoids(Ellipsoids,x,y,z);

Ellipses=Ellipsoids_to_Ellipses(Ellipsoids,dist,r_min);
 
Plot_Ellipses(Ellipses,x,z);

[Nodes_Fibers, Fibers]=Generate_Fiber(x,y,z,L,N,DFiber,Orientation,Ndiv,Ellipsoids);

Plot_Fiber(x,y,z,Nodes_Fibers,Fibers,DFiber);

Plot_Ellipsoids_Fiber(Ellipsoids,x,y,z,Nodes_Fibers,Fibers,DFiber);

[Walter Roberson]

Are you using the functions from https://www.mathworks.com/matlabcentral/fileexchange/87644-aggregate-packing-generator-with-thick-fibers ?

What is the difference between the output you are getting now and the output you want to get?

[Serigne Saliou Sene]

Yes this is the function I used.

I must add a 6th plot in which I would have the aggregates and the fibers in the packaging (Volume) of the cube clearly visible in order to form the ITZ, as in the photo. Thanks in advance

[Yahriel Salinas-Reyes]

Hello,

I am running into the similar problem and trying to plot the Interfacial Transition Zone (ITZ) all in one 3d plot. Where you able to find a solution or have any help to offer?