MATLAB: Is it possible to get similar size of polygons or with very less variation

computationalgeometry

https://uk.mathworks.com/matlabcentral/answers/426087-random-polygons-inside-a-semi-circle

Best Answer

yes it is

Related Solutions

MATLAB: Random polygons inside a semi-circle

Edit code:

make polygon more uniform in size at the trade off of randomness remove artifact of the border.
adjust repulsion parameters.
improve sticky-border artifact due to denser points

N = 100; % aproximative number of polygonals to be generated
n = 64; % control size and number of vertexes of polygonal
nrepulsion = 8; % control the size of the polygonal and the randomness of the position
X = randn(N,2);
R = sqrt(rand(N,1));
X = R .* X ./ sqrt(sum(X.^2,2));
X(:,2) = abs(X(:,2));
nb = 100;
theta = linspace(0,pi,nb)';
XC = [cos(theta), sin(theta)];
nb = ceil(nb*2/pi);
XY0 = linspace(-1,1,nb)' .* [1 0];
XY0([1 end],:) = [];
n1 = size(X,1);
n2 = size(XC,1);
n3 = size(XY0,1);
CC = n1+(1:n2-1)' + [0 1];
C0 = (n1+n2)+(1:n3-1)' + [0 1];
C = [CC; C0];
XB = [XC; XY0];
% Repulsion of seeds to avoid them to be too close to each other
for k = 1:nrepulsion-1
    XALL = [X; XB];
    DT = delaunayTriangulation(XALL);
    T = DT.ConnectivityList;
    containX = ismember(T,1:n1);
    b = any(containX,2);
    TX = T(b,:);
    [r,i0] = find(containX(b,:));
    i = mod(i0+(-1:1),3)+1;
    i = TX(r + (i-1)*size(TX,1));
    T = accumarray([i(:,1);i(:,1)],[i(:,2);i(:,3)],[n1 1],@(x) {x}); 
    maxd2 = 0;
    R = zeros(n1,2);
    for i=1:n1
        Ti = T{i};
        P = X(i,:) - XALL(Ti,:);
        nP2 = sum(P.^2,2);
        maxd2 = maxd2 + mean(nP2);
        b = Ti > n1;
        nP2(b) = nP2(b)*3; % reduce repulsion from each point of the border
        R(i,:) = sum(P./nP2,1);
    end
    if k==1
        v0 = 0.005/sqrt(maxd2/n1);
    end
    v = v0/sqrt(max(sum(R.^2,2)));
    X = X + v*R;
    % Project back if points falling outside the half-circle
    X(:,2) = max(X(:,2),0.01);
    r2 = sum(X.^2,2);
    out = r2>1;
    X(out,:) = X(out,:) .* (0.99 ./ sqrt(r2(out)));
end
DT = delaunayTriangulation(X);
[V,P] = voronoiDiagram(DT);
KX = convexHull(DT);
[ib,ik] = ismember(1:N,KX);
r = 2;
r2 = r^2;
warning('off','MATLAB:polyshape:boundary3Points');
warning('off','MATLAB:polyshape:repairedBySimplify');
PXB = polyshape(XC);
for k=1:N
    Pk = V(P{k},:);
    if ib(k) % infinity
        ik0 = ik(k);
          if ik0 == length(KX)
              kp = KX(1);
          else
              kp = KX(ik0+1);
          end
          km = k;
          Pv = Pk(2,:);
          if Pv*Pv' < r2
              t = X(km,:)-X(kp,:);
              nt2 = t*t';
              P0 = t*((Pv*t')/nt2);
              cs2 = P0*P0';
              if cs2 <= r2
                  d = [-t(2),t(1)] / sqrt(nt2);
                  sn = sqrt(r2-cs2);
                  if sn > (Pv-P0)*d'
                      Pk(1,:) = P0 + sn*d;
                  else
                      Pk(1,:) = [];
                  end
              else
                  Pk(1,:) = [];
              end
          else
              Pk(1,:) = [];
          end       
          if ik0 == 1
              km = KX(end);
          else
              km = KX(ik0-1);
          end
          kp = k;
          Pv = Pk(end,:);
          if Pv*Pv' < r2
              t = X(km,:)-X(kp,:);
              nt2 = t*t';
              P0 = t*((Pv*t')/nt2);
              cs2 = P0*P0';
              if cs2 <= r2
                  d = [-t(2),t(1)] / sqrt(nt2);
                  sn = sqrt(r2-cs2);
                  if sn > (Pv-P0)*d'
                      Pk(end+1,:) = P0 + sn*d;
                  end
              end
          end  
      end
      [Pk,sid] = intersect(PXB, polyshape(Pk));
      Pk = Pk.Vertices;
      m = length(Pk);
      if m >= 3
          nb1 = sum(sid == 1);
          m = m - nb1;
          W = rand(n,m-1) .^ (1./(m-1:-1:1));
          W = cumprod([ones(n,1),W],2) .* (1-[W, zeros(n,1)]);
          if nb1>0
              % Consider weight of the borders as weight for two points
              Pk = circshift(Pk,-find(sid==0, 1, 'first'),1);
              nb1 = nb1+2;
              w = linspace(0,2/nb1,nb1);
              Pk = [Pk(1:m-2,:); [w; fliplr(w)]*Pk(m-1:end,:)];
          end
          Pk = W*Pk;
          K = convhull(Pk);
          P{k} = Pk(K,:);
      else
          P{k} = [];
      end
  end
  P(cellfun('isempty',P)) = [];
% Check
fig = figure(1);
clf(fig);
ax = axes('Parent',fig);
hold(ax,'on');
plot(ax, XB([1:end 1],1),XB([1:end 1],2),'k');
for k=1:length(P)
    Pk = P{k};
    fill(ax,Pk(:,1),Pk(:,2),k);
end
axis(ax,'equal');
axis(ax,[-1.1 1.1 -0.1 1.1]);

MATLAB: Fill area with random circles having different diameters

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

Best Answer

Related Solutions

MATLAB: Random polygons inside a semi-circle

MATLAB: Fill area with random circles having different diameters

Related Question