MATLAB: Pseudo-random distribution of points with minimum distance

distancepointspseudo random

Hi all,

I use the following code to generate an exact number of points (650 in this case) within a one hectar area:

numPoints = 650; % number of points to be generated
width = 100; % length in m of one square side
x = 0;
y = 0;
figure('Position', [300 300 900 900])
rectangle('Position', [x, y, width, width],'LineWidth',2,'LineStyle','--');
grid on;
hold on;
xRandom = 50 + (width * rand(1, numPoints) - width / 2);
yRandom = 50 + (width * rand(1, numPoints) - width / 2);
plot(xRandom, yRandom, 'b.', 'MarkerSize', 8);
hold on;
plot(xRandom, yRandom, 'ro', 'MarkerSize', 20);
title(['',num2str(numPoints),' points inside one hectar'], 'Interpreter', 'None');
xlabel('length in meters');
ylabel('length in meters');
axis equal tight;

What I would like to achieve is that the minimum distance between these points is 3 meters. This condition is fulfilled if the red circles (with 3 m diameter – only approximated here using MarkerSize) only touch each other but don't overlap as they currently do (see image below).

Does anybody know how to accomplish this?

PS: at a completely regular spacing of 650 points in one hectar there would be 4 m distance between each point. Maybe a possible solution would be to create a regular spacing first and then add or substract a smaller random increment?

Best Answer

See Bruno Luong's code under

https://de.mathworks.com/matlabcentral/answers/432516-model-of-a-crowd-on-concert-venue-or-how-to-distribute-random-points-according-to-the-2d-window-dist

Related Solutions

MATLAB: Random Distribution Points inside rectangle

Try this:

numPoints = 1000;
width = 2200;
x = -1100
y = -1100;
rectangle('Position', [x, y, width, width],'LineWidth',2,'LineStyle','--');
grid on;
xRandom = width * rand(1, numPoints) - width / 2;
yRandom = width * rand(1, numPoints) - width / 2;
hold on;
plot(xRandom, yRandom, 'r.', 'MarkerSize', 10);
title('Random Points Inside a Rectangle', 'FontSize', fontSize, 'Interpreter', 'None');
xlabel('X', 'FontSize', fontSize);
ylabel('Y', 'FontSize', fontSize);
axis equal;
% Put axis at origin
ax = gca;
ax.XAxisLocation = 'origin';
ax.YAxisLocation = 'origin';
% Set up figure properties:
% Enlarge figure to full screen.
set(gcf, 'Units', 'Normalized', 'OuterPosition', [0 0 1 1]);
% Get rid of tool bar and pulldown menus that are along top of figure.
set(gcf, 'Toolbar', 'none', 'Menu', 'none');
% Give a name to the title bar.
set(gcf, 'Name', 'Demo by ImageAnalyst', 'NumberTitle', 'Off')

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

Best Answer

Related Solutions

MATLAB: Random Distribution Points inside rectangle

MATLAB: Random polygons inside a semi-circle

Related Question