Pseudo-random distribution of points with minimum distance

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?

[James Tursa]

If there are only a few points, you could just push the offenders away from each other. But if you have a large number of points, this will turn into more of a packing problem than a random distribution problem. Is that acceptable?

[Florian]

I need to end up with the exact number, so eliminating offenders (as proposed in previous solutions) is not an option. Maybe 'rand' is not the best function to start with...

[James Tursa]

I wasn't suggesting eliminating the offenders, but pushing them away from each other until they are not offending anymore.

Answer 1

[Torsten]

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

[Florian]

Below is what I get using this code. The points are still too close together and I cannot see where to set a minimum distance of 3 m in thsi code. Could you perhaps point me to the right line in the code?

[Bruno Luong]

I cannot see where to set a minimum distance of 3 m in thsi code

dmin = 0.5;

[Torsten]

dmin ?

And in your case, you should set a 2D uniform distribution instead of the 2D epanechnikov-distribution, I guess.

[Florian]

Thanks guys!

[Bruno Luong]

For reference I put here the code with uniform distribution

L = 100; % <-- Choose length of square sides
x0 = 0; y0 = 0; % <-- Choose center of square
n = 300; % <-- Choose number of points
% Generate uniform-distribution on the square
X = rand(n,2)*L - (L/2*[1,1]+[x0,y0]);
XYR = [x0,y0]+[[-1;1;1;-1;-1],[-1;-1;1;1;-1]]*L/2;
XB = interp1((0:4)'*L,XYR,linspace(0,4*L,200));
XB(end,:) = [];
nrepulsion = 500;
% Repulsion of seeds to avoid them to be too close to each other
n = size(X,1);
Xmin = [x0-L/2,y0-L/2];
Xmax = [x0+L/2,y0+L/2];
% Point on boundary
XR = x0+[-1,1,1,-1,-1]*L/2;
YR = y0+[-1,-1,1,1,-1]*L/2;
cla;
hold on
plot(XR,YR,'r-');
h = plot(X(:,1),X(:,2),'b.');
axis equal
dmin = 3; % minima distance between 2 objects
d2min = dmin*dmin;
beta = 0.5;
for k = 1:nrepulsion
    XALL = [X; XB];
    DT = delaunayTriangulation(XALL);
    T = DT.ConnectivityList;
    containX = ismember(T,1:n);
    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)],[n 1],@(x) {x}); 
    maxd2 = 0;
    R = zeros(n,2);
    move = false(n,1);
    for i=1:n
        Ti = T{i};
        P = X(i,:) - XALL(Ti,:);
        nP2 = sum(P.^2,2);
        if any(nP2<4*d2min)
            move(i) = true;
            move(Ti(Ti<=n)) = true;
        end
        maxd2 = maxd2 + mean(nP2);
        b = Ti > n;
        nP2(b) = nP2(b)*5; % reduce repulsion from each point of the border
        R(i,:) = sum(P./max((nP2-d2min),1e-3),1);
    end
    if ~any(move)
        break
    end  
    if k==1
        v0 = (L*5e-3)/sqrt(maxd2/n);
    end
    R = R(move,:);
    v = v0/sqrt(max(sum(R.^2,2)));
    X(move,:) = X(move,:) + v*R;
    
    % Project back if points falling outside the rectangle
    X = min(max(X,Xmin),Xmax);
    
    set(h,'XData',X(:,1),'YData',X(:,2));
    pause(0.01);
end
theta = linspace(0,2*pi,65);
xc = dmin/2*sin(theta);
yc = dmin/2*cos(theta);
% plot circles f diameter dmin around random points
for i=1:n
    plot(X(i,1)+xc,X(i,2)+yc,'k');
end

[Florian]

Perfect. Thanks Bruno!

Answer 2

[Image Analyst]

See my attached demo that I've posted before.

[Florian]

This script doesn't seem to work for me using MATLAB Version: 9.2.0.556344 - R2017a).

I use maxPoints = 1000 and numPoints = 650 with minClosestDistance = 0.03;

The script gives me a message that 650 points have been computed, but when I open x and y there are only NaNs in there.

[Image Analyst]

You should not have set numPoints = 650. Then it will only place an additional 350 points. You should have left it at zero and let it do the count by itself.

When I ran it, I tried a million times to place 1000 points and I was only able to place 789. It looks like Bruno's algorithm works though, at least for the one time I tried it.