Fill area with random circles having different diameters

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?

[KSSV]

How many number of circles? Any further conditions?

[Image Analyst]

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.

Answer 1

[Anton Semechko]

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

[mikymike89]

It's perfect the result from your posted image. Unfortunately, if I use the code, I have the error code "All functions in a script must be closed with an 'end'".

% 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);

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

[Walter Roberson]

Put one more 'end' at the end of the existing function if you are including it in a script.

[mikymike89]

Thanks for your suggestion. But I have already tried to put one more "end" at the end of each function but it does not work

[Anton Semechko]

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]

Thanks, works perfectly. But if I add at the end the following code I write an external file having only the coordinate of the circle with x=0. How can I plot all the generated coordinates?

FileName=sprintf('%s%sRandom circle packing data.csv',cd,filesep);
csvwrite(FileName,[C R]);
fprintf('Data saved to %s\n',FileName)

[Walter Roberson]

If you have the Image Processing Toolbox,

viscircles(C, R)

[Salad Box]

Is 'SampleLineSegment' a pre-defined function? I can't find it anywhere!

[Walter Roberson]

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]

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]

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.

Answer 2

[Walter Roberson]

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

[Walter Roberson]

By the way, there are other algorithms that are much more efficient.

[mikymike89]

Thanks so much for the code. how can I generate an external output file with the circle coordinates?

[Walter Roberson]

coords = [rx(:), ry(:), rr(:)];
csvwrite('YourOutput.csv', coords);

[Walter Roberson]

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]

1st,2nd point- the coordinates and radii can be continuous 3rd point- the touch is permitted 4th,5th,6th,7th point- I don't need a more precise code. It is already excellent for my goal. It doesn't matter if the filling is not so much accurate and remain some space.

I have extra requests if you can do it and I'm sorry if I do it just now.

Before creating circles on the surface, can I generate random circles also in the 4 corners and then on the 4 edges (having the center of the circle along the edge)?

Could it be possible to change the code to generalize it even in the case of a rectangular shape?

[Walter Roberson]

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]

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]

@Walter

I really want to make it up to you.

Thank you

[davide cademartori]

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.

Answer 3

[Anton Semechko]

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

[mikymike89]

Thanks so much. I have a very poor knowledge of Matlab and i'm sorry for the trivial questions, how can I increase the iteration number to better fill the area? Moreover, how can I generate an external output file with the circle coordinates?

[Image Analyst]

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]

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]

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.

Answer 4

[Anton Semechko]

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,:)=[];

[mikymike89]

this last code is very efficient for my case because fill the area well and even really quickly. I did not understand how I can export the coordinates. Is possible to generate automatically an external file with "wrtite" command?

Moreover, I would have the same extra requests made to Walter:

-Before creating circles on the surface, can I generate random circles also in the 4 corners and then on the 4 edges (having the center of the circle along the edge)?

- Could it be possible to change the code to generalize it even in the case of a rectangular shape?

Thanks

[Anton Semechko]

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]

It's perfect. How can I maximize the fill on the edge before filling the area? Thank you so much for your time.

[Anton Semechko]

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]

Yes, exactly. I would prefer to insert the circles before on the edge and only then in the area.

[Image Analyst]

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]

I know that is not very natural but the circles on the edges are placed only in the 4 middle points of it so far. My issue impose that the boundary must be good filled otherwise I'd have a problem for the some following analyses

[Jesus Obregon]

Someone please could tell me how works the ishandle function, I don't undesrtand the matlab available explanation

[Walter Roberson]

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]

i want to change the rectangle to circle or semi-circle . . how to do that ?

[jaehee jeong]

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]

Right after

      n=n+1;

You can insert

      if n > TheLimitYouWant; break; end

[jaehee jeong]

@walter ahha! Love you so much!! Thank you

[jaehee jeong]

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]

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]

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

[Walter Roberson]

See https://www.mathworks.com/matlabcentral/answers/405186-fill-area-with-random-circles-having-different-diameters#comment_892769

[RAJESH]

@Anton Semechko, Could help me with Random Position of short Fiber(RVE) and Generate Fiber.

Answer 5

[jahanzaib ahmad]

can some one do such packing with polygons ? without using delauny triagulation

[Walter Roberson]

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]

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]

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]

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]

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 .

https://www.researchgate.net/publication/224739203_Two-dimensional_packing_for_irregular_shaped_objects

i m not trying to implement this paper but just for an idea .

[Walter Roberson]

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]

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 ?

Answer 6

[Mohammad Faris]

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?

[jahanzaib ahmad]

[jahanzaib ahmad]

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]

@ 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]

[Kaitlyn Casulli]

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]

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]

the above code is not running