MATLAB: Random polygons inside a semi-circle

polygonsramdom polygons

how to generate 2D random convex polygons inside a semi-circle 150mm diameter.the size of polygons should have a limit . maximum length of side of polygons is between 5 to 10 mm not smaller nor larger

Best Answer

Edit code:
  1. make polygon more uniform in size at the trade off of randomness remove artifact of the border.
  2. adjust repulsion parameters.
  3. 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]);
Related Question