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.

[Jeroen Houwen]

Dear Image Analyst,

I am trying to place the points at a certain distance but the distance value has a normal distribution.

I have adapted your code to try and get a normal distribution in the distance between the points but this does not seem to give a correct normal distribution. Do you know of any previous questions or a solution to my problem? I am trying to use this for a lot of points

distances = sqrt((x-xPossible) .^ 2 + (y - yPossible) .^ 2);

minClosestDistance = normrnd(ClosestDistance,1);

if min(distances) >= minClosestDistance

[Image Analyst]

@Jeroen Houwen, I'm not sure you can enforce a Gaussian distribution. Why do you think that the average spacing should follow a Gaussian distribution? Since there is a minimum spacing, it definitely won't be Gaussian. Even if that min spacing were 0 I don't know that it would be Gaussian. Is there some theory that says it should be Gaussian instead of uniform? Even Uniform wouldn't describe it because the points are not allowed to be outside the rectangle so the distribution tapers down as points get farther away.

Did you construct and look at the Fry plot? The average radial profile looks more log normal.

[Jeroen Houwen]

@Image Analyst there is no theory stating it should be one or the other. I need a Gaussian distribution because I am trying to model metal grains of which the diameter follows a Gaussian distribution. Therefore I am trying to place points at a certain distance which follows the normal distribution and then apply Voronoi.

I will look into Fry plot since I didn't look at that!

Thank you for your fast response!

[Bruno Luong]

"model metal grains of which the diameter follows a Gaussian distribution"

I'm not aware that diameter can be negative.

[Jeroen Houwen]

@Bruno Luong I meant a Gaussian distribution with mean=10. I should have stated my problem a bit better

[Bruno Luong]

Gaussian distribution has non zero probability of negative values. Your request doesn' make any sense, since there is no distance and radius with negative values, so the process can never be a Gaussian distribution.

[Image Analyst]

With a Gaussian distribution, and a finite number of samples drawn from it, it's certainly possible to not encounter any negative values. For example if the mean were 10 and the SD were 0.5 and you drew 10 thousand samples, it's possible that none of those 10,000 numbers would be negative.

% Get 10 thousand numbers from a Gaussian Distribution
r = 10 + 0.5 * randn(10000, 1);
% Find out what the min and max are
fprintf('Min r = %f\nMax r = %f\n', min(r), max(r));
fprintf('Mean r = %f\nStd Dev r = %f\n', mean(r), std(r));
% Fit data, r, to a Gaussian/Normal distribution 
% even though there are no negative numbers.
pd = fitdist(r, 'Normal')

Min r = 7.783199

Max r = 12.011163

Mean r = 10.002350

Std Dev r = 0.497673

pd =

Normal distribution

mu = 10.0011 [9.99139, 10.0108]

sigma = 0.495396 [0.488625, 0.502359]

That said @Jeroen Houwen, looking at my Fry Plot demo below you can see the distribution for distances bounded in a rectangle looks more log normal than normal.

[Bruno Luong]

"With a Gaussian distribution, and a finite number of samples drawn from it, it's certainly possible to not encounter any negative values."

Fine but that has nothing to do my comment, which is one can NEVER generate a Gaussian random distribution with a positive values. I'm not talking about the reverse.

Answer 3

[Image Analyst]

@Jeroen Houwen, see attached demo of Fry Plot.