Fill area with random circles having different diameters
Show older comments
I should fill the area of a 500x500 square with random circles having random diameters between 10 and 50 (without overlap). Then, I need the output file of the generated coordinates. Can someone help me, please?
2 Comments
KSSV
on 12 Jun 2018
How many number of circles? Any further conditions?
Image Analyst
on 12 Jun 2018
Is this homework? Do you want the circles in a digital image (if so, what resolution), or do you want them in an overlay (graphical plane) like you can do with plot(), scatter(), or viscircles()? Be sure to look at Walter's answer below.
Accepted Answer
Anton Semechko
on 14 Jun 2018
Hopefully this iteration of code (see below) is the last one.
% Unconstrained circle packing example. Only centroids constrained to the rectangle.
ab=[500 500]; % rectangle dimensions; [width height]
R_min=5; % minimum circle radius
R_max=25; % maximum circle radius
cnst=false;
[C,R]=random_circle_packing_rectangle(ab,R_min,R_max,cnst);
% Constrained circle packing example
ab=[500 500]; % rectangle dimensions
R_min=5; % minimum circle radius
R_max=25; % maximum circle radius
cnst=true;
[C,R]=random_circle_packing_rectangle(ab,R_min,R_max,cnst);
% Example of how to export data
FileName=sprintf('%s%sRandom circle packing data.csv',cd,filesep);
csvwrite(FileName,[C R]);
fprintf('Data saved to %s\n',FileName)
============================================================================================================
function [C,R]=random_circle_packing_rectangle(ab,R_min,R_max,cnst,vis)
% Random circle packing inside a rectangle.
%
% INPUT:
% - ab : 1-by-2 vector specify rectangle dimensions; ab=[width height].
% ab=[500 500] is the default setting. Coordinates of the
% lower-left and upper-right corners of the rectangle are
% (0,0) and (ab(1),ab(2)), respectively.
% - R_min : minimum circle radius. R_min=min(ab)/100 is the default setting.
% - R_max : maximum circle radius. R_max=min(ab)/20 is the default setting.
% - cnst : set cnst=true to ensure all circles fit into the rectangle.
% cnst=false is the default setting, meaning that only circle
% centroids will be constrained to the boundary and interior of
% the rectangle.
% - vis : set vis=false to suppress visualization. vis=true is the
% default setting.
%
% OUTPUT:
% - C : Q-by-2 array of sphere centroids
% - r : Q-by-1 array of sphere radii
%
% Default settings
if nargin<1 || isempty(ab)
ab=[500 500];
elseif numel(ab)~=2 || ~isnumeric(ab) || ~ismatrix(ab) || sum(ab(:)<1E-6 | isinf(ab(:)))>0
error('Invalid entry for 1st input argument (ab)')
else
ab=ab(:)';
end
if nargin<2 || isempty(R_min)
R_min=min(ab)/100;
elseif numel(R_min)~=1 || ~isnumeric(R_min) || R_min>(min(ab)/4-2E-12) || R_min<min(ab)/1E3
error('Invalid entry for 2nd input argument (D_min)')
end
if nargin<3 || isempty(R_max)
R_max=min(ab)/20;
elseif numel(R_max)~=1 || ~isnumeric(R_max) || R_max<R_min || R_max>(min(ab)/4-2E-12)
error('Invalid entry for 3rd input argument (D_max)')
end
if nargin<4 || isempty(cnst)
cnst=false;
elseif numel(cnst)~=1 || ~islogical(cnst)
error('Invalid entry for 4th input argument (cnst)')
end
if nargin<5 || isempty(vis)
vis=true;
elseif numel(vis)~=1 || ~islogical(vis)
error('Invalid entry for 5th input argument (vis)')
end
% Grid used to keep track of unoccupied space inside the rectangle
dx=max(min(ab)/2E3,R_min/50);
x=0:dx:ab(1);
y=0:dx:ab(2);
[x,y]=meshgrid(x,y);
G=[x(:) y(:)];
clear x y
% Start by placing circles along the edges if cnst=false
dR=R_max-R_min;
Cc=bsxfun(@times,ab,[0 0; 1 0; 1 1; 0 1]); % corner vertices
if ~cnst
[Xa,Xb]=deal(Cc,circshift(Cc,[-1 0]));
Rc=dR*rand(4,1)+R_min;
[Rc_a,Rc_b]=deal(Rc,circshift(Rc,[-1 0]));
[C,R]=deal(cell(4,1));
for i=1:4
[Ci,Ri]=SampleLineSegment(Xa(i,:),Xb(i,:),Rc_a(i),Rc_b(i),R_min,R_max);
Ci(end,:)=[];
C{i}=Ci;
Ri(end)=[];
R{i}=Ri;
end
C=cell2mat(C);
R=cell2mat(R);
% Update grid
for i=1:size(C,1), G=update_grid(C(i,:),R(i),G,R_min); end
else
% Remove all grid points less than R_min units from the boundary
G_max=G+R_min+1E-12;
G_min=G-R_min-1E-12;
chk_in=bsxfun(@le,G_max,ab) & bsxfun(@ge,G_min,[0 0]);
chk_in=sum(chk_in,2)==2;
G(~chk_in,:)=[];
clear G_max G_min chk_in
C=[]; R=[];
end
% Begin visualization
if vis
hf=figure('color','w');
axis equal
set(gca,'XLim',[0 ab(1)]+ab(1)*[-1/20 1/20],'YLim',[0 ab(2)]+ab(2)*[-1/20 1/20],'box','on')
hold on
drawnow
Hg=plot(G(:,1),G(:,2),'.k','MarkerSize',2);
t=linspace(0,2*pi,1E2)';
P=[cos(t) sin(t)]; % unit circle
if ~isempty(C)
for i=1:size(C,1)
Pm=bsxfun(@plus,R(i)*P,C(i,:));
h=fill(Pm(:,1),Pm(:,2),'r');
set(h,'FaceAlpha',0.25)
end
drawnow
end
end
f=5:5:100;
fprintf('Progress : ')
fprintf(' %-3u ',f)
fprintf(' (%%complete)\n')
fprintf(' ')
Ng=size(G,1);
f=size(G,1)-round((f/100)*Ng);
% Use rejection sampling to populate interior of the rectangle
flag=true;
n=0; cnt=0; m=0;
while ~isempty(G) && cnt<1E4 && (~vis || (vis && ishandle(hf)))
n=n+1;
% New circle
if flag && (cnt>500 || size(G,1)<0.95*Ng)
flag=false;
Rg=R_max*ones(size(G,1),1);
end
i=[];
if cnt<=500 && flag
X_new=ab.*rand(1,2); % centroid
else
i=randi(size(G,1));
X_new=G(i,:)+(dx/2)*(2*rand(1,2)-1);
X_new=min(max(X_new,[0 0]),ab);
if cnt>1E3
Rg(:)=max(0.95*Rg,R_min);
end
end
if isempty(i)
R_new=dR*rand(1)+R_min; % radius
else
R_new=(Rg(i)-R_min)*rand(1)+R_min;
end
% Check if the circle fits inside the rectangle when cnst=true
if cnst
X_new_max=X_new+R_new+1E-12;
X_new_min=X_new-R_new-1E-12;
chk_in=X_new_max<=ab & X_new_min>=[0 0];
if sum(chk_in)<2
cnt=cnt+1;
continue
end
end
% Does the new circle intersect with any other circles?
if ~isempty(C)
d_in=sqrt(sum(bsxfun(@minus,C,X_new).^2,2));
id=d_in<(R+R_new);
if sum(id)>0
cnt=cnt+1;
if ~isempty(i)
Rg(i)=min(0.95*Rg(i),min(0.99*(R_new+dx/2),R_max));
Rg(i)=max(Rg(i),R_min);
end
continue
end
end
% Accept new circle
cnt=0;
m=m+1;
C=cat(1,C,X_new);
R=cat(1,R,R_new);
[G,id]=update_grid(X_new,R_new,G,R_min);
if ~flag, Rg(id)=[]; end
% Visualization
if vis && ishandle(hf)
Pm=bsxfun(@plus,R_new*P,X_new);
h=fill(Pm(:,1),Pm(:,2),'r');
set(h,'FaceAlpha',0.25)
if mod(m,50)==0
delete(Hg)
Hg=plot(G(:,1),G(:,2),'.k','MarkerSize',2);
drawnow
end
end
% Progress update
if size(G,1)<=f(1)
f(1)=[];
fprintf('* ')
end
end
fprintf('\n')
% Show rectangle
if vis && ishandle(hf)
Cc=[Cc;Cc(1,:)];
plot(Cc(:,1),Cc(:,2),'--k','LineWidth',2)
delete(Hg)
end
if nargout<1, clear C R; end
function [G,id]=update_grid(X_new,R_new,G,R_min)
% Remove grid points within R_new+R_min units of new circle
D=sum(bsxfun(@minus,G,X_new).^2,2);
id=D<(R_new+R_min+1E-12)^2;
G(id,:)=[];
function [C,R]=SampleLineSegment(Xa,Xb,Ra,Rb,R_min,R_max)
% Place circles along line segment between points Xa and Xb
r=Xb-Xa;
L=norm(r);
r=r/norm(L);
dR=R_max-R_min;
C=Xa; R=Ra;
while true
R_new=dR*rand(1)+R_min;
C_new=C(end,:)+r*(R(end)+R_new+R_max*rand(1));
D=L - norm(C_new + r*(R_new+Rb) - Xa); % will there be enough space left for the end point with radius Rb?
if D<2*(R_min+1E-12)
C=cat(1,C,Xb);
R=cat(1,R,Rb);
break
else
C=cat(1,C,C_new);
R=cat(1,R,R_new);
end
end
10 Comments
mikymike89
on 22 Aug 2018
Edited: mikymike89
on 22 Aug 2018
Walter Roberson
on 22 Aug 2018
Put one more 'end' at the end of the existing function if you are including it in a script.
mikymike89
on 22 Aug 2018
Anton Semechko
on 22 Aug 2018
Code I posted above had some minor bugs. Attached is the final version of the code, free of bugs, and with improved boundary sampling.
mikymike89
on 22 Aug 2018
Walter Roberson
on 22 Aug 2018
If you have the Image Processing Toolbox,
viscircles(C, R)
Walter Roberson
on 22 Sep 2022
Look further down in the same post, the function is defined there.
function [C,R]=SampleLineSegment(Xa,Xb,Ra,Rb,R_min,R_max)
% Place circles along line segment between points Xa and Xb
Philip Steinwascher
on 12 Jan 2023
This code is great. I was looking for something almost like that. How would you need to edit the code for fitting a variable number of different circles for which one sets the exact radii plus a fixed distance between them?
Walter Roberson
on 12 Jan 2023
Requiring a fixed distance between N-dimensional spheres in an N-dimensional space always has a maximum of (N+1) objects placed. After (N+1), you can place additional objects but they will not be the given distance from all other points.
More Answers (5)
Walter Roberson
on 12 Jun 2018
How do you know when to give up on the random placement?
N = 500;
minr = 10; maxr = 50;
maxtries = 500000;
intcoords = false;
sq = zeros(N, N, 'uint8');
[X, Y] = ndgrid(1:N, 1:N);
nc = 0;
iteration = 0;
trynum = 0;
maxgoodtry = 0;
fmt = 'iteration #%d, nc = %d, try #%d';
fig = gcf;
set(fig, 'Units', 'pixels', 'Position', [0 0 N+30, N+30]);
image(sq);
colormap(jet(2))
axis([0 N+1 0 N+1]);
title(sprintf(fmt, iteration, nc, trynum));
drawnow();
rx = []; ry = []; rr = [];
while trynum <= maxtries
iteration = iteration + 1;
if intcoords
r = randi([minr, maxr]);
cxy = randi([r+1, N-r], 1, 2);
else
r = minr + rand(1,1) * (maxr-minr);
cxy = r + 1 + rand(1,2) * (N - 2*r - 1);
end
mask = (X-cxy(1)).^2 + (Y-cxy(2)).^2 <= r^2;
if nnz(sq & mask) > 0
trynum = trynum + 1;
else
sq = sq | mask;
image(sq);
nc = nc + 1;
rx(nc) = cxy(1); ry(nc) = cxy(2); rr(nc) = r;
title(sprintf(fmt, iteration, nc, trynum));
drawnow();
maxgoodtry = max(maxgoodtry, trynum);
trynum = 1;
end
end
fprintf('finished at iteration %d, hardest success took %d tries\n', iteration, maxgoodtry);
fprintf('Number of circles: %d\n', length(rx));
9 Comments
Walter Roberson
on 12 Jun 2018
By the way, there are other algorithms that are much more efficient.
mikymike89
on 12 Jun 2018
Walter Roberson
on 12 Jun 2018
coords = [rx(:), ry(:), rr(:)];
csvwrite('YourOutput.csv', coords);
Walter Roberson
on 12 Jun 2018
In my testing, I found cases where locations could still be found after 350000 rejections.
The number of circles that can be packed with those radii is generally between 170 and 220, partly dependent upon how persistent you are.
I know a much more efficient method than I used above, but there is no point in implementing it until fundamental questions are answered:
- are the coordinates of the circle centers intended to be integers or continuous?
- are the radii intended to be integers or continuous?
- The question states that overlap is not permitted, but is touching of circles permitted?
- Is it okay to approximate whether a square is within a circle according to whether the center of the square is within the circle, or is it necessary to be more precise about it? Like if 1% of a square is theoretically occupied by a circle to its left, does that "consume" the whole circle because a small fraction of it is used, forbidding the case where a circle on its right also uses 1%? In this particular left/right example, the infinitely precise versions of the circles would not overlap, but clearly it is also possible for there to be cases where two circles overlap by just a sliver within a unit square, with neither one of them happening to pass over the center of the circle
- Is it necessary to keep going until every gap is smaller than the minimum radius, to pack as densely as possible, even though a rejection method could take an indefinitely large amount of time to "randomly" locate the last few gaps?
- Are there conditions on the ratio of circle radii? That is, as you place the larger circles, you quickly get to the point where large circles will no longer fit, but you might be able to fit four-ish times as many circles of half the radii. So if you pack as densely as possible, then a histogram of radii would be highly biased towards the smaller diameters. But some circle packing situations ask for fairly uniform distribution of radii, expecting you to stop when you can no longer keep the distribution relatively equal, even though that would leave gaps that could be filled by more smaller circles
- some circle packing situations do not expect you to fill every possible gap, recognizing that random selection of coordinates might not be able to locate the gaps. But they do expect you not to give up immediately. If you are not expecting that every possible gap be filled, then how many tries of being rejected is "good enough" ?
mikymike89
on 12 Jun 2018
Walter Roberson
on 13 Jun 2018
Generalization to rectangle, and also puts a circle in each corner and a circle on each of the four edges.
%M x N, height by width
M = 480;
N = 640;
minr = 10; maxr = 50;
maxtries = 500000;
sq = zeros(M, N, 'uint8');
[Y, X] = ndgrid(1:M, 1:N);
nc = 0;
iteration = 0;
trynum = 0;
maxgoodtry = 0;
fmt = 'iteration #%d, nc = %d, try #%d';
fig = gcf;
set(fig, 'Units', 'pixels', 'Position', [0 0 N+30, M+30]);
image(sq);
colormap(jet(2))
axis equal
axis([0 N+1 0 M+1]);
title(sprintf(fmt, iteration, nc, trynum));
drawnow();
rx = []; ry = []; rr = [];
while trynum <= maxtries
iteration = iteration + 1;
r = randin(minr, maxr);
switch nc
case 0
cx = 1; cy = 1;
case 1
cx = 1; cy = M;
case 2
cx = N; cy = 1;
case 3
cx = N; cy = M;
case 4
cx = randin(r, N-r);
cy = 1;
case 5
cx = randin(r, N-r);
cy = M;
case 6
cx = 1;
cy = randin(r, M-r);
case 7
cx = N;
cy = randin(r, M-r);
otherwise
cx = randin(r, N-r);
cy = randin(r, M-r);
end
mask = (X-cx).^2 + (Y-cy).^2 <= r^2;
if nnz(sq & mask) > 0
trynum = trynum + 1;
else
sq = sq | mask;
image(sq);
nc = nc + 1;
rx(nc) = cx; ry(nc) = cx; rr(nc) = r;
title(sprintf(fmt, iteration, nc, trynum));
drawnow();
maxgoodtry = max(maxgoodtry, trynum);
trynum = 1;
end
end
fprintf('finished at iteration %d, hardest success took %d tries\n', iteration, maxgoodtry);
fprintf('Number of circles: %d\n', length(rx));
And you will need
function r = randin(minval, maxval, varargin)
r = rand(varargin{:}) * (maxval - minval) + minval;
end
To reduce computation time, you can use a lower maxtries limit. In my tests, it sometimes takes over 390000 tries to randomly place a circle in one of the few remaining gaps.
Note: this code does not implement my improved algorithm at all.
Walter Roberson
on 10 Jun 2020
With multiple ranges of permitted radii, based upon the first set of code.
%M x N, height by width
M = 240;
N = 320;
resolution = 100;
minr = ceil([0.14, 0.07] * resolution);
maxr = floor([0.18, 0.1] * resolution);
num_needed = 38;
num_first_type = 28;
circle_type = [1 * ones(1,num_first_type), 2 * ones(1,num_needed - num_first_type)];
circle_type = circle_type(randperm(length(circle_type)));
maxtries = 500000;
sq = zeros(M, N, 'uint8');
[Y, X] = ndgrid(1:M, 1:N);
intcoords = false;
sq = zeros(N, N, 'uint8');
[X, Y] = ndgrid(1:N, 1:N);
nc = 0;
iteration = 0;
trynum = 0;
maxgoodtry = 0;
fmt = 'iteration #%d, nc = %d, try #%d';
fig = gcf;
set(fig, 'Units', 'pixels', 'Position', [0 0 N+30, N+30]);
image(sq);
colormap(jet(2))
axis([0 N+1 0 N+1]);
title(sprintf(fmt, iteration, nc, trynum));
drawnow();
rx = zeros(1,num_needed);
ry = zeros(1,num_needed);
rr = zeros(1,num_needed);
while nc < num_needed && trynum <= maxtries
ct = circle_type(nc+1);
iteration = iteration + 1;
if intcoords
r = randi([minr(ct), maxr(ct)]);
cxy = randi([r+1, N-r], 1, 2);
else
r = minr(ct) + rand(1,1) * (maxr(ct)-minr(ct));
cxy = r + 1 + rand(1,2) * (N - 2*r - 1);
end
mask = (X-cxy(1)).^2 + (Y-cxy(2)).^2 <= r^2;
if nnz(sq & mask) > 0
trynum = trynum + 1;
else
sq = sq | mask;
image(sq);
nc = nc + 1;
rx(nc) = cxy(1); ry(nc) = cxy(2); rr(nc) = r;
title(sprintf(fmt, iteration, nc, trynum));
drawnow();
maxgoodtry = max(maxgoodtry, trynum);
trynum = 1;
end
end
fprintf('finished at iteration %d, hardest success took %d tries\n', iteration, maxgoodtry);
fprintf('Number of circles: %d\n', length(rx));
jaehee jeong
on 11 Jun 2020
@Walter
I really want to make it up to you.
Thank you
davide cademartori
on 19 Jan 2023
Hello, thanks for all your work.
I have a couple of further questions related to the code:
1)Is it possible to implement the partial overlapping of the circles up to r/2 of the new circle introduced at every iteration?
2)Do you have any suggestion on how to implement a criterion which stops the surface filling when there is a continuous path of connected circles from the bottom to the top of rectangle/square?
Thanks again for all the work you have already done.
Anton Semechko
on 12 Jun 2018
Edited: Anton Semechko
on 12 Jun 2018
Here is a demo that uses rejection sampling and Delaunay triangulation
function [C,r]=random_circle_packing_demo
%
% - C : Q-by-2 array of sphere centroids
% - r : Q-by-1 array of sphere radii
%
% Problem settings
a=500; % size of the square
D_min=10; % min diameter of inscribed circle
D_max=50; % max diameter of inscribed circle
% Start by placing randomply spaced points along the boundary
d_min=D_min/a;
d_max=D_max/a;
d=d_max-d_min;
Xb=cell(4,1);
for i=1:4
t=[];
t_net=0;
while t_net<1
t=cat(1,t,d*rand(1)+d_min);
t_net=t_net+t(end);
end
t=t/t_net;
[t_min,id_min]=min(t);
while t_min<d_min
t_new=d*rand(1)+d_min;
if t_new>t_min
dt=t_new-t_min;
if t_new>t_min*(1+dt)
t(id_min)=t_new;
t=t/(1+dt);
[t_min,id_min]=min(t);
end
end
end
t=cat(1,0,t);
t(end)=[];
Xi=ones(numel(t),2);
if i==1
Xi(:,1)=a*cumsum(t);
Xi(:,2)=0;
elseif i==2
Xi(:,1)=a;
Xi(:,2)=a*cumsum(t);
elseif i==3
Xi(:,1)=a*(1-cumsum(t));
Xi(:,2)=a;
else
Xi(:,1)=0;
Xi(:,2)=a*(1-cumsum(t));
end
Xb{i}=Xi;
if i==1
figure('color','w')
axis equal
set(gca,'XLim',[0 a]+a*[-1/20 1/20],'YLim',[0 a]+a*[-1/20 1/20],'box','on')
hold on
drawnow
end
plot(Xi(:,1),Xi(:,2),'.k','MarkerSize',10)
end
Xb=cell2mat(Xb);
% Use rejection sampling to populate interior of the square
X=Xb;
D=D_max-D_min;
n=0; cnt=0;
while cnt<100 % terminate when no additional points can be inserted after 100 iterations
n=n+1;
%disp([n size(X,1) cnt])
% New point
X_new=a*rand(1,2);
% Randomly chosen distance threshold, D_thr, between D_min and D_max. The point
% process being simulated depends on how this parameter is chosen. If
% you want the spheres to be packed closer together, bias D_thr towards
% D_min.
D_thr=D*rand(1)+D_min;
% Reject new point if it is closer to X than D_thr
d_out=sum(bsxfun(@minus,X,X_new).^2,2);
if min(d_out)<D_thr^2
cnt=cnt+1;
continue
end
% Accept new point
cnt=0;
X=cat(1,X,X_new);
plot(X_new(:,1),X_new(:,2),'.k','MarkerSize',10)
if mod(n,10)==0, drawnow; end
end
% Delaunay triangulation
DT=delaunayTriangulation(X);
triplot(DT)
[C,r]=incenter(DT);
id=r<D_min/2 | r>D_max/2;
C(id,:)=[];
r(id)=[];
plot(C(:,1),C(:,2),'.r','MarkerSize',10)
N=size(C,1);
% Plot the spheres
figure('color','w')
axis equal
set(gca,'XLim',[0 a]+a*[-1/20 1/20],'YLim',[0 a]+a*[-1/20 1/20],'box','on')
hold on
t=linspace(0,2*pi,1E2)';
P=[cos(t) sin(t)];
for i=1:N
Pi=bsxfun(@plus,r(i)*P,C(i,:));
plot(Pi(:,1),Pi(:,2),'-k')
if mod(i,10)==0, drawnow; end
end
4 Comments
mikymike89
on 12 Jun 2018
Image Analyst
on 12 Jun 2018
Change 100 to a bigger number to give it more trials to find a diameter that fits:
while cnt<100
Alternatively you could use the Euclidean Distance Transform to immediately locate where the largest circle that can fit would like. However computing the EDT does take some time so you wouldn't want to do it when the field is mostly empty but when you get to where you're rejecting a large number of circles before finding a "good" one, then it may be advantageous to use the EDT.
Anton Semechko
on 12 Jun 2018
To better fill in the area, set 'D_thr' variable inside the while loop to D_min:
D_thr=D_min;
Simplest way to "export" results of the simulation is to call
[C,r]=random_circle_packing_demo;
from command window. Once the simulation is complete, right click on the 'C' variable (centroids) in the workspace, select "Copy", then paste "C" into spreadsheet. Repeat the same for 'r' variable (radii).
Anton Semechko
on 12 Jun 2018
Actually, setting D_thr to D_min will not produce tighter packing of the circles. After sampling is completed I think that the Voronoi cells corresponding to circles with radii below D_min will have to be modified locally using Lloyd's relaxation algorithm to get tighter packing. I don't really have time to experiment with this. Maybe someone else on here can help you.
Anton Semechko
on 12 Jun 2018
Here is another implementation of random circle packing. This version produces much tighter packing than my previous demo thanks to Walter's idea of using a grid to keep track of unoccupied space inside the square.
function [C,R]=random_circle_packing_demo2
%
% - C : Q-by-2 array of sphere centroids
% - r : Q-by-1 array of sphere radii
%
% Problem settings
a=500; % size of the square
D_min=10; % min circle diameter
D_max=50; % max circle diameter
% Unit circle
t=linspace(0,2*pi,1E2)';
P=[cos(t) sin(t)];
% Sampling grid
x=linspace(0,a,2E3);
dx=x(2)-x(1);
[x,y]=meshgrid(x);
G=[x(:) y(:)];
clear x y
% Begin visualization
hf=figure('color','w');
axis equal
set(gca,'XLim',[0 a]+a*[-1/20 1/20],'YLim',[0 a]+a*[-1/20 1/20],'box','on')
hold on
drawnow
f=5:5:100;
fprintf('Progress : ')
fprintf(' %-3u ',f)
fprintf(' (%%complete)\n')
fprintf(' ')
f=size(G,1)-round((f/100)*size(G,1));
% Use rejection sampling to populate interior of the square
C=[]; R=[];
D=D_max-D_min;
flag=true;
n=0; cnt=0; m=0;
while ~isempty(G) && ishandle(hf)
n=n+1;
% New circle
if cnt>500, flag=false; end
if cnt<500 && flag
X_new=a*rand(1,2); % centroid
else
i=randi(size(G,1));
X_new=G(i,:)+(dx/2)*(2*rand(1,2)-1);
end
R_new=(D*rand(1)+D_min)/2; % radius
% Does the new circle intersect with any other circles?
if n>1
d_in=sqrt(sum(bsxfun(@minus,C,X_new).^2,2));
id=d_in<(R+R_new);
if sum(id)>0
cnt=cnt+1;
continue
end
end
% Accept new circle
cnt=0;
m=m+1;
C=cat(1,C,X_new);
R=cat(1,R,R_new);
G=update_grid(X_new,R_new,G,D_min/2);
% Visualization
if ishandle(hf)
Pm=bsxfun(@plus,R_new*P,X_new);
h=fill(Pm(:,1),Pm(:,2),'r');
set(h,'FaceAlpha',0.25)
if mod(m,10)==0, drawnow; end
end
% Progress update
if size(G,1)<=f(1)
f(1)=[];
fprintf('* ')
end
end
fprintf('\n')
if nargout<1
clear C R
end
function G=update_grid(X_new,R_new,G,R_min)
% Remove grid points within R_new+R_min units of new circle
D=sum(bsxfun(@minus,G,X_new).^2,2);
id=D<(R_new+R_min+1E-12)^2;
G(id,:)=[];
18 Comments
mikymike89
on 12 Jun 2018
Anton Semechko
on 13 Jun 2018
Modified code that incorporates requested changes is attached. Data generated by the function will be automatically exported to the 'Random circle packing data.csv' file in your working directory. You can open this file in Excel. First two columns will contain (x,y) co-ordinates of the circle centroids and third column will contain the radii. To change dimensions of the rectangle, modify variable 'ab' at the top of the function. For example, to pack circles inside 500-by-200 rectangle, set ab=[500 200].
function [C,R]=random_circle_packing_demo3
%
% - C : Q-by-2 array of sphere centroids
% - r : Q-by-1 array of sphere radii
%
% Problem settings
ab=[500 200]; % rectangle dimensions; ab=[width height]
D_min=10; % min circle diameter
D_max=50; % max circle diameter
% Unit circle
t=linspace(0,2*pi,1E2)';
P=[cos(t) sin(t)];
% Sampling grid
dx=min(ab)/2E3;
x=0:dx:ab(1);
y=0:dx:ab(2);
[x,y]=meshgrid(x,y);
G=[x(:) y(:)];
clear x y
% Begin visualization
hf=figure('color','w');
axis equal
set(gca,'XLim',[0 ab(1)]+ab(1)*[-1/20 1/20],'YLim',[0 ab(2)]+ab(2)*[-1/20 1/20],'box','on')
hold on
drawnow
f=5:5:100;
fprintf('Progress : ')
fprintf(' %-3u ',f)
fprintf(' (%%complete)\n')
fprintf(' ')
f=size(G,1)-round((f/100)*size(G,1));
% Place four circles at the corners and another four at the edge midpoints
Cc=[0 0; 1 0; 1 1; 0 1]; % corners
Ce=(Cc+circshift(Cc,[-1 0]))/2; % edge midpoints
C=bsxfun(@times,ab,[Cc;Ce]);
D=D_max-D_min;
R=(D*rand(8,1)+D_min)/2;
for i=1:8
G=update_grid(C(i,:),R(i),G,D_min/2);
Pm=bsxfun(@plus,R(i)*P,C(i,:));
h=fill(Pm(:,1),Pm(:,2),'r');
set(h,'FaceAlpha',0.25)
end
drawnow
% Use rejection sampling to populate interior of the square
flag=true;
n=0; cnt=0; m=0;
while ~isempty(G) && ishandle(hf)
n=n+1;
% New circle
if cnt>500, flag=false; end
if cnt<500 && flag
X_new=ab.*rand(1,2); % centroid
else
i=randi(size(G,1));
X_new=G(i,:)+(dx/2)*(2*rand(1,2)-1);
end
R_new=(D*rand(1)+D_min)/2; % radius
% Does the new circle intersect with any other circles?
d_in=sqrt(sum(bsxfun(@minus,C,X_new).^2,2));
id=d_in<(R+R_new);
if sum(id)>0
cnt=cnt+1;
continue
end
% Accept new circle
cnt=0;
m=m+1;
C=cat(1,C,X_new);
R=cat(1,R,R_new);
G=update_grid(X_new,R_new,G,D_min/2);
% Visualization
if ishandle(hf)
Pm=bsxfun(@plus,R_new*P,X_new);
h=fill(Pm(:,1),Pm(:,2),'r');
set(h,'FaceAlpha',0.25)
if mod(m,10)==0, drawnow; end
end
% Progress update
if size(G,1)<=f(1)
f(1)=[];
fprintf('* ')
end
end
fprintf('\n')
% Show rectangle
Cc=[0 0; 1 0; 1 1; 0 1];
Cc=bsxfun(@times,ab,Cc);
Cc=[Cc;Cc(1,:)];
plot(Cc(:,1),Cc(:,2),'--k','LineWidth',2)
% Write to file
FileName=sprintf('%s%sRandom circle packing data.csv',cd,filesep);
csvwrite(FileName,[C R]);
fprintf('Data saved to %s\n',FileName)
if nargout<1
clear C R
end
function G=update_grid(X_new,R_new,G,R_min)
% Remove grid points within R_new+R_min units of new circle
D=sum(bsxfun(@minus,G,X_new).^2,2);
id=D<(R_new+R_min+1E-12)^2;
G(id,:)=[];
mikymike89
on 13 Jun 2018
Anton Semechko
on 13 Jun 2018
Glad to help. It was an interesting problem.
I am not sure what you mean by "maximizing fill on the edge". Do you want to sample the boundary of the rectangle first and then its interior?
mikymike89
on 13 Jun 2018
Image Analyst
on 13 Jun 2018
It may not look very natural (not random sizes) but to maximize fill along the outer edges you could just line up the smallest discs possible shoulder-to-shoulder all the way around the rectangle.
mikymike89
on 13 Jun 2018
Jesus Obregon
on 9 Oct 2018
Someone please could tell me how works the ishandle function, I don't undesrtand the matlab available explanation
Walter Roberson
on 9 Oct 2018
ishandle(hf)
is testing to see whether hf is still a valid handle to something. The case where it would not be a valid handle to something is the case where somehow the figure has stopped existing, such as if the user clicked on the border X to close the figure while the code was running.
jahanzaib ahmad
on 26 Nov 2018
i want to change the rectangle to circle or semi-circle . . how to do that ?
jaehee jeong
on 7 Jun 2020
Hello, I want to limit the number of circles,
How can I limit the number of circles in the rectangle?
the code keep going until whole circles occupy the area.
Thank you
Walter Roberson
on 7 Jun 2020
Right after
n=n+1;
You can insert
if n > TheLimitYouWant; break; end
jaehee jeong
on 8 Jun 2020
@walter ahha! Love you so much!! Thank you
jaehee jeong
on 8 Jun 2020
ah one more question,, sorry for bothering..
Is it possible to fill with other circles that have different range of radius?
I tried but it failed.. I'm new to matlab and having troubleㅠ.ㅠ
so, I would like to have radius between 1. min 0.14~max0.18, and want to add circles :min 0.07~ max 0.1
summary: Circle 1: number of 28, radius= min 0.14~max0.18
circle 2= number of 10, radius= min 0.07~ max 0.1
please help me ㅠㅡㅠ
Walter Roberson
on 8 Jun 2020
Edited: Walter Roberson
on 8 Jun 2020
When you say "number of 28" do you mean 28 of those random circles? So you would be looking for a total of 28+10 = 38 random circles?
This code is Anton's logic and I have not gone through the logic in detail. Unfortunately Anton does not appear to have contributed in the last two years so I do not know if Anton will see your question and respond.
If I were modifying my code from https://www.mathworks.com/matlabcentral/answers/405186-fill-area-with-random-circles-having-different-diameters#answer_324277 then I would probably use something like
resolution = 100;
minr = ceil([0.14, 0.07] * resolution);
maxr = floor([0.18, 0.1] * resolution);
circle_type = [1 * ones(1,28), 2 * ones(1,10)];
circle_type = circle_type(randperm(length(circle_type)));
And then as I looped through for the K'th circle instead of using randi([minr, maxr]) it would be
ct = circle_types(K);
randi([minr(ct), maxr(ct)])
That is, my logic would be to create a fixed pool of circle types, taken in random order, and at each point I would place a circle according to whatever size whose "turn" it is. This easily generalizes to any (finite) number of circle types.
jaehee jeong
on 9 Jun 2020
Edited: jaehee jeong
on 9 Jun 2020
I put your feedback under this 'while' loop, and it appeared error in 'cxy = r + 1 + rand(1,2) * (N - 2*r - 1);' part
Do I need to erase condition above while loop?
The number of matrix dont fit..sorry i'm really new but i need to do this..
N = 1; // because I want the result rectangle 1*1
resolution = 100;
minr = ceil([0.14, 0.07] * resolution);
maxr = floor([0.18, 0.1] * resolution);
circle_type = [1 * ones(1,28), 2 * ones(1,10)];
circle_type = circle_type(randperm(length(circle_type))); (new code)
while trynum <= maxtries
iteration = iteration + 1;
if intcoords
ct = circle_types(K);
randi([minr(ct), maxr(ct)]) (your new code)
cxy = randi([r+1, N-r], 1, 2);
else
r = minr + rand(1,1) * (maxr-minr);
cxy = r + 1 + rand(1,2) * (N - 2*r - 1);
end
mask = (X-cxy(1)).^2 + (Y-cxy(2)).^2 <= r^2;
if nnz(sq & mask) > 0
trynum = trynum + 1;
else
sq = sq | mask;
image(sq);
nc = nc + 1;
rx(nc) = cxy(1); ry(nc) = cxy(2); rr(nc) = r;
title(sprintf(fmt, iteration, nc, trynum));
drawnow();
maxgoodtry = max(maxgoodtry, trynum);
trynum = 1;
end
end
RAJESH
on 19 Jan 2023
jahanzaib ahmad
on 16 Feb 2019
Edited: jahanzaib ahmad
on 16 Feb 2019
0 votes
can some one do such packing with polygons ? without using delauny triagulation
7 Comments
Walter Roberson
on 16 Feb 2019
It is possible, but more complicated.
There is a case that is relative easy to handle. If the polygons are regular, and the minimum distance between them exceeds the difference between furthest and closest point to the centroid, and the packing is not required to be minimal: then in that case you can proceed through placement of circles, with the circles being cicumscriptions of the polygons, and then once the polygon is placed, rorate it around a random angle for its precise placement.
However, if the polygons are permitted to nearly touch, then the exact orientation of each polygon can matter and much more care needs to be taken in placement.
If the poygons are not regular polygons, then unless the minimum distance between them is high and you do not need any kind of minimum packing, then orientation of placement can become a significant problem. Randomly dropping non-overlapping toothpicks until it is no longer possible to fit any more is a tougher problem than randomly dropping squares.
jahanzaib ahmad
on 17 Feb 2019
my problem is polygons are permitted to touch each other and polygons are convex .
..can u please give some more information .what i have read so far those ways are very time consuming .take and place method and other is translation method .
Walter Roberson
on 23 Feb 2019
Would it be acceptable to confine the polygons to a small number of rotational angles?
Does the solution have to apply to arbitrary convex polygons? Regular n-gons are easier to deal with due to the rotational symmetry and the clustering around the centroid.
What is the stopping criteria to be?
jahanzaib ahmad
on 25 Feb 2019
Walter Roberson yes arbitrary convex polygons are fine but with some max and minimum size polygons . the stopping criteria is maximum packing .
if it is possible .oherwise e.g N polygons of convex shape be fixed in square .
but not by delauny triangulation method .
jahanzaib ahmad
on 25 Feb 2019
Edited: jahanzaib ahmad
on 25 Feb 2019
i m trying to fix polygons first i give them random position and then i check intersection if intersect then again i translate them to random position .if i increase the number of polygons the time is increase very much high . on the other hand thanslating them to +rand take them to (1,1) as translate function add to previous co-ordinates . i tried to translate them negative random number .but still i m not happy with this ..so futher sugestion required .plus i m beginner so i make code of everything instead of functions . .i m using polyxpoly function .
i m not trying to implement this paper but just for an idea .
Walter Roberson
on 25 Feb 2019
There is a notable difference between what the paper talks about and what you are talking about. In the paper, they are permitted to rotate the object at need, but your requirements appear to be that the polygons are to be at random orientations. The algorithms in the paper are permitted to use heuristics and computations to figure out where to best put a piece, whereas your requirements appear to be that pieces are added at random locations if they fit.
The paper deals with packing algorithms. Your question appears to deal with a situation similar to if you were to pour a bunch of objects into a box, remove the ones that overlap something beneath it, pour those removed ones on top, remove the ones that overlap, and so on, until finally no more pieces can possibly fit.
jahanzaib ahmad
on 25 Feb 2019
yes.. Walter Roberson i cant see any material like polygon packing or random placement on mathworks .so i m doing that . i have made code just trying to improve it and i will post it ..
i m trying is randomly placing it and then checking overlap and placeing it to new position. but such a method is slow .can u please suggest any good way e.g i m trying now it that if a polygon has occupied some place in square i dont want to produce random number in that place so that new polygons has few chances to goto that place again so nimber of overlaps/intersections will reduce .? what do u say ? or any other method i can search for free space in a box and place new ploygono there in random way ?
Mohammad Faris
on 20 Feb 2019
0 votes
Hey everybody.
i have some issue in 3d dimension
i wanna fill a cylinder with some spherical glasses. can we rework this method for 3d dimension?
7 Comments
jahanzaib ahmad
on 4 Jul 2019
jahanzaib ahmad
on 4 Jul 2019
i made this using random take and place method very conventional method .but computational time is very high .
1- generating
2-placing
3 finding interrsection
4 translating ( or rotating ) and again finding intersection
jahanzaib ahmad
on 4 Jul 2019
Edited: Walter Roberson
on 4 Jul 2019
@ Walter Roberson i m trying this way to make polygons packing .generating a random centroid and making a convex polygon and then using scale function i place that polygon at that centroid position .and then generating till the required area percentage condition fullfilled .please guide me for improvement and i want to avoid polyxpoly function i think it make the code slow. 50 percent area is easy but 70 percent takes alot of time ..plus if there is intersection i m reducing the size of the polygon thats again not a good way because i want polygon between two min and max size .
to measure the size of the polygon i m using minboudrect funciton .(attached )
rA=0.5; % 50 percent area
m=1;
N=9;
x=rand; % random centroid x and y
y=rand;
X=rand(N,2);
XX=X(:,1);
XY=X(:,2);
K=convhull(XX,XY);
AB=[XX(K),XY(K)];
polyin = polyshape(AB);
seivesize=[4.75 9.5]; % reuired size of polygon
D=(seivesize(1)+((seivesize(2)-(seivesize(1))).*rand))/100; % random size of polygon
[rx,ry] = minboundrect(AB(:,1),AB(:,2));
side1=sqrt(((rx(1)-rx(2))^2)+(ry(1)-ry(2))^2);
side2=sqrt(((rx(2)-rx(3))^2)+(ry(2)-ry(3))^2);
sizeofpolygons1=min(side1,side2);
FACTOR=D/sizeofpolygons1; % reducing the size of the polygon by factor
poly1 = scale(polyin,FACTOR,[x,y]);
%plot(poly1);
area2=area(poly1);
AC=[poly1.Vertices(:,1),poly1.Vertices(:,2);poly1.Vertices(1,1),poly1.Vertices(1,2)];
A{m}=AC;
rA=rA-area2; Sarea=0;
%m=m+1;
while rA >= area2
chk3=1;
while (chk3)~=0
x=rand;
y=rand;
P=[x,y];
chk2=[];
for i=1:length(A)
AD=A{i};
AD1=polyshape(AD);
chk1=isinterior(AD1,P);
chk2=[chk2,chk1];
end
chk3=any(chk2);
end
N=9;
X=rand(N,2);
XX=X(:,1);
XY=X(:,2);
K=convhull(XX,XY);
AB=[XX(K),XY(K)];
polyin = polyshape(AB);
seivesize=[4.75 9.5];
D=(seivesize(1)+((seivesize(2)-(seivesize(1))).*rand))/100;
[rx,ry] = minboundrect(AB(:,1),AB(:,2));
side1=sqrt(((rx(1)-rx(2))^2)+(ry(1)-ry(2))^2);
side2=sqrt(((rx(2)-rx(3))^2)+(ry(2)-ry(3))^2);
sizeofpolygons1=min(side1,side2);
FACTOR=D/sizeofpolygons1;
poly1 = scale(polyin,FACTOR,[x,y]);
%area1=area(poly1);
AC=[poly1.Vertices(:,1),poly1.Vertices(:,2);poly1.Vertices(1,1),poly1.Vertices(1,2)];
for j=1:length(A)
AF=A{j};
[Cx,Ca] = polyxpoly(AF(:,1), AF(:,2),AC(:,1),AC(:,2));
while ~isempty(Cx)
polyin1 = polyshape(AC);
FACTOR1=.9;
poly1 = scale(polyin1,FACTOR1,[x,y]);
AC=[poly1.Vertices(:,1),poly1.Vertices(:,2);poly1.Vertices(1,1),poly1.Vertices(1,2)];
[Cx,Ca] = polyxpoly(AF(:,1), AF(:,2),AC(:,1),AC(:,2));
end
end
m=m+1;
A{m}=AC;
polyin1 = polyshape(AC);
area2=area( polyin1);
Sarea=Sarea+area2;
rA=rA- area2;
end
axis equal
axis([-0.1 1.1 -0.1 1.1]);
figure(2);
for i =A
plot(i{:}(:,1), i{:}(:,2),'k');
hold on;
end
jahanzaib ahmad
on 4 Jul 2019
Kaitlyn Casulli
on 3 Apr 2020
Looks like this thread got away from the original question... but I am also interested in doing a 3D arrangement. I was able to find some starter code online, but the time it takes is ridiculously long because it does not have a grid inside and rather tests ALL possible points in the volume. When you're trying to have a closely packed arrangement, this gets difficult to find an empty spot very quickly! I've also noticed some issues where it allows for overlapping, which I do not want.
I am interested in filling a rectangular area with solid spheres (or, ideally, ellipsoids) within certain diameter constraints.
Yi Liu
on 31 Dec 2020
The following function does fill a cylinder space with radom sized ellipsoids without overlapping. The resulting ellipsoids are colour coded by their sizes (small-blue to large-orange). The following plots are generated by calling:
[C,R]=random_ellipsoid_packing_demo(500,600,0.3);
The function is:
function [C,R]=random_ellipsoid_packing_demo(d,h,ffactor)
% Randomly pack different size ellipsoids in a cylinder value with dim d,
% height h, and a filling factor ffactor (0.0 - 1.0).
%
% - C : Q-by-3 array of ellipsiod centroids
% - R : Q-by-3 array of ellipsiod radii
%
% Problem settings
%a=500; % size of the square
tvol = pi*(d/2)^2*h;
D_min=5; % min ellipsiod diameter
D_max=25; % max ellipsiod diameter
% Unit circle
t=linspace(0,2*pi,1E2)';
P=[cos(t) sin(t)];
% Unit ellipsiod (sphere)
[xs,ys,zs]=ellipsoid(0,0,0,1,1,1);
% Sampling grid
x=linspace(0,d,1E3);
y=linspace(0,d,1E3);
z=linspace(0,h,1E3);
dx=x(2)-x(1);
[x,y,z]=ndgrid(x,y,z);
G=[x(:) y(:) z(:)];
clear x y z
cmap=jet;
Ridx = linspace(D_min,D_max,256);
% Begin visualization
hf=figure('color','w');
plot3(d/2*P(:,1)+d/2,d/2*P(:,2)+d/2,h*ones(size(t)),':r')
hold on
plot3(d/2*P(:,1)+d/2,d/2*P(:,2)+d/2,0*ones(size(t)),':r')
axis equal
set(gca,'XLim',[0 d]+d*[-1/20 1/20],'YLim',[0 d]+d*[-1/20 1/20],'ZLim',[0 h]+h*[-1/20 1/20],'box','on')
drawnow
C=zeros(1,3); R=zeros(1,3);
D=D_max-D_min;
flag=true;
n=0; cnt=0; m=0;
fillfact=0;
fillper=0;
while ~isempty(G) && ishandle(hf) && fillfact<ffactor*tvol
n=n+1;
% New ellipsiod
if cnt>500, flag=false; end
if cnt<500 && flag
X_new=[d,d,h].*rand(1,3); % centroid
else
i=randi(size(G,1));
X_new=G(i,:)+[dx/2,dx/2,h/2].*(2*rand(1,3)-1);
end
R_new = D*rand(1,3)+D_min; % radius
% Does the ellipsiod intersect with boundary (cylinder)?
if ((X_new(1)-d/2-R_new(1))^2+(X_new(2)-d/2-R_new(2))^2) > (d/2)^2 ...
|| ((X_new(1)-d/2+R_new(1))^2+(X_new(2)-d/2+R_new(2))^2) > (d/2)^2 ... % partly outside boundary
|| (X_new(3)-R_new(3))<0 || (X_new(3)+R_new(3))>h % top or bottom
continue
end
% Does the new ellipsiod intersect with any other ellipsiods?
if n>1
d_in=sqrt(sum(bsxfun(@minus,C,X_new).^2,2));
id=d_in<(R+R_new);
if (sum(id)>0)
cnt=cnt+1;
continue
end
end
% Accept new ellipsiod
cnt=0;
m=m+1;
C=cat(1,C,X_new);
R=cat(1,R,R_new);
fillfact=fillfact+4/3*pi*R_new(1)*R_new(2)*R_new(3);
idx = find(mean(R_new)<Ridx,1);
ccmap=zs*0;
ccmap(:,:,1)=cmap(idx,1);
ccmap(:,:,2)=cmap(idx,2);
ccmap(:,:,3)=cmap(idx,3);
% Visualization
if ishandle(hf)
hh=surf(xs*R_new(1)+X_new(1),ys*R_new(2)+X_new(2),zs*R_new(3)+X_new(3),ccmap);
set(hh,'EdgeColor','None','FaceAlpha',0.3,'FaceColor',cmap(idx,:))
if mod(m,20)==0, drawnow; end
end
% Progress update
if mod(fillfact/tvol*100,5)<0.05 && ceil(fillfact/tvol*100)>fillper
fillper=ceil(fillfact/tvol*100);
fprintf('Progress (%% filled): %s\n',num2str(fillfact/tvol*100,'%2.1f'))
end
end
title(['Fill Factor = ',num2str(fillfact/tvol*100,'%.1f'),'%'])
figure
[np,bin]=hist(R(:),ceil(D));
yyaxis left
bar(bin,np)
ylabel('Number of Ellipsoids')
xlabel('Radius')
yyaxis right
plot(bin,np.*bin.^3*pi*4/3,'-*')
ylabel('Volume of Ellipsoids at Radius')
title('Statistics')
if nargout<1
clear C R
end
The resulting plots are:
Stats of the filling are:
RAJESH
on 13 Jan 2023
the above code is not running
Categories
Find more on Surface and Mesh Plots in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
