Create random points in a rectangular domain, but with minimum separation distance

Hello,

I am not a regular user of this program. Hence this query.

I am interested to create random points in a 2D rectangular space, to start with, a uniform distribution. I know how to do this. However, I want to ensure that each of these points are separated by a certain minimum distance.

Can anyone help me with the script please.

Thank you.

Suresh

Answer 1

[Image Analyst]

You basically have to try and see. If it's too close, reject that point and get a new one.

x = rand(1, 10000);
y = rand(1, 10000);
minAllowableDistance = 0.05;
numberOfPoints = 300;
% Initialize first point.
keeperX = x(1);
keeperY = y(1);
% Try dropping down more points.
counter = 2;
for k = 2 : numberOfPoints
  % Get a trial point.
  thisX = x(k);
  thisY = y(k);
  % See how far is is away from existing keeper points.
  distances = sqrt((thisX-keeperX).^2 + (thisY - keeperY).^2);
  minDistance = min(distances);
  if minDistance >= minAllowableDistance
    keeperX(counter) = thisX;
    keeperY(counter) = thisY;
    counter = counter + 1;
  end
end
plot(keeperX, keeperY, 'b*');
grid on;

[SS]

Hi,

Thanks a lot. This is perfect. You have saved a lot of my time. I see the trick was in the variable "distances", where you have performed vector operation. This is something I did not know and I was doing it point by point, which was ridiculous.

If you don't mind extending your help again, what I was finally interested is to draw randomly (e.g. uniform distribution of angle theta, randomly varying between 0 to 180 degrees) inclined line segments (of fixed length and also if possible variable length, in which case the variable length follows some distribution) through these random points, these points acting as the center of the line segments.

Suresh

[Image Analyst]

See this code:

% Creates a random pattern of points with no point closer to another than some specified distance.
% Also puts a ray emanating from each point at random angle and random length.
% Initialize some things.
x = rand(1, 10000);
y = rand(1, 10000);
minAllowableDistance = 0.05;
numberOfPoints = 300;
rMax = 0.25;
% Initialize first point.
keeperX = x(1);
keeperY = y(1);
% Try dropping down more points.
counter = 2;
for k = 2 : numberOfPoints
  % Get a trial point.
  thisX = x(k);
  thisY = y(k);
  % See how far is is away from existing keeper points.
  distances = sqrt((thisX-keeperX).^2 + (thisY - keeperY).^2);
  minDistance = min(distances);
  if minDistance >= minAllowableDistance
    keeperX(counter) = thisX;
    keeperY(counter) = thisY;
    % Draw a line from it.
    % First get a random angle
    angle = 2*pi*rand(1);
    % Next get a random length
    r = rMax * rand(1);
    % Now get endpoint of the ray.
    x2 = thisX + r * cos(angle);
    y2 = thisY + r * sin(angle);
    % Plot the point
    plot(keeperX(counter), keeperY(counter), 'bo', 'LineWidth', 3);
    hold on;
    % Draw the line.
    line([thisX x2], [thisY, y2], 'Color', 'r', 'LineWidth', 3);
    counter = counter + 1;
  end
end
grid on;
% Enlarge figure to full screen.
set(gcf, 'Units', 'Normalized', 'OuterPosition', [0 0 1 1]);

[SS]

Hi,

Once again, perfect and exactly what I wanted. In the meantime, I achieved the same as above but through another approach, wherein I solved an equation of a line and a circle with the random points acting as the centre of a circle and the line passing through the centre. But I like your approach the most.

Thank you for taking time to help me. Much appreciated.

S

[Prashant Katiyar]

Is there anyway we can replace points to small circles of variable diameters? Also if we can variate the distance between each of these circles? Some thing like below?

[William Harvie]

ngl this is probably the least effecient way to do this but it works

[Image Analyst]

@William Harvie, please post a more efficient way to do it if you know of one.

@Prashant Katiyar, try this:

%============================================================================================================================================
% Initialization Steps.
clc;    % Clear the command window.
close all;  % Close all figures (except those of imtool.)
clear;  % Erase all existing variables. Or clearvars if you want.
workspace;  % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 25;
% Get a list of trial (x,y) circle centers.
x = 100 * rand(1, 1000000);
y = 100 * rand(1, 1000000);
% Specify how many circles are desired.
numberOfCircles = 300;
minRadius = 1;
maxRadius = 5;
% Specify how close the centers may be to each other.
% Should be at least twice the max radius if they are not to touch.
minAllowableDistance = max([11, 2 * maxRadius]);
% Initialize first point.
keeperX = x(1);
keeperY = y(1);
radii = minRadius;
% Try dropping down more points.
counter = 2;
for k = 2 : length(x)
    % Get a trial point.
    thisX = x(k);
    thisY = y(k);
    % See how far is is away from existing keeper points.
    distances = sqrt((thisX-keeperX).^2 + (thisY - keeperY).^2);
    minDistance = min(distances);
    if minDistance >= minAllowableDistance
  	  % This trial location works.  Save it.
      keeperX(counter) = thisX;
      keeperY(counter) = thisY;
  	% Get a random radius.
    radii(counter) = minRadius + (maxRadius - minRadius) * rand;
    % Quit if we have enough or ran out of circles to try
    if counter >= numberOfCircles
        break;
    end
    counter = counter + 1;
    end
end
% Plot a dot at the centers
plot(keeperX, keeperY, 'b+', 'MarkerSize', 9);
viscircles([keeperX(:), keeperY(:)], radii, 'Color', 'b');
grid on;
axis equal
numCirclesPlaced = length(keeperX);
caption = sprintf('Could place %d circles', numCirclesPlaced);
title(caption, 'fontSize', fontSize)
g = gcf;
g.WindowState = 'maximized'

[Joan]

@Image Analyst, thank you very much for your code. It is showing me some direction. However, I am a bit stuck. Please add some code that would make the circles enclosed inside a polyshape of coordinates say,

xv=[0 0 20 35 66 85 105 105]; % X-Coordinates

yv=[0 20 20 40 40 20 16 0]; % Y-Coordinates

[tyler]

@Image Analyst

>please post a more efficient way to do it if you know of one

Hello, I know this thread is many years old but I spent several weeks trying to figure this out and came across this resource:

https://research.wmz.ninja/articles/2018/06/matlab-plots-can-be-beautiful-iv.html

It does work quite well, however it does not behave when expanding to 3D; it can literally paint itself into a corner (and in fact does seem to do this with more regularity than not). I'm pretty sure this is because the scatter grows like a mold: there's a max distance a new point can be from any other point, so there'll always be a bias towards where the seed point is. Anyways, I just wanted to share this because it's been giving me a headache and I did find what appears to be a reasonably efficient 2D solution. I'm also copying the code here for posterity, in case anything happens to the link.

function points = poisson_disc_2d(w, h, n, r, k)
%POISSON_DISC_2D Generate 2D points using the Bridson's Poisson-disc
%sampling algorithm.
%Inputs:
%   w - Width of the region.
%   h - Height of the region.
%   n - Maximum number of points.
%   r - Radius of the disc.
%   k - (Optional) Maximum number of trials. Default value is 30.
%Output:
%   points - n x 2 matrix where each row represents a point. n may be less
%            than the specified n.
if nargin < 5
    k = 30;
end
% Background grid for fast occupancy search.
grid_size = r / sqrt(2);
n_w = ceil(w / grid_size);
n_h = ceil(h / grid_size);
grid = zeros(n_w, n_h);
% Generate the initial point.
points = zeros(n, 2);
points(1, :) = [unifrnd(0, w) unifrnd(0, h)];
active_points = [1];
n_points_generated = 1;
grid(floor(points(1,1) / grid_size) + 1, floor(points(1,2) / grid_size) + 1) = 1;
% Generate the remaining points
while ~isempty(active_points) && n_points_generated < n
    % Pick a random point from the active set.
    idx = randi(size(active_points, 1));
    p = points(active_points(idx), :);
    % Try generate a new point from this point within the annulus [r, 2r].
    success = false;
    for kk = 1:k
        % Uniformly samples a point within an annulus.
        new_r = sqrt(r * r + 3 * r * r * rand());
        new_theta = unifrnd(0, pi*2);
        new_p = [new_r * cos(new_theta) + p(1)...
                 new_r * sin(new_theta) + p(2)];
        % Check validity.
        if new_p(1) < 0 || new_p(2) < 0 || new_p(1) > w || new_p(2) > h
            continue;
        end
        grid_x = floor(new_p(1) / grid_size) + 1;
        grid_y = floor(new_p(2) / grid_size) + 1;
        if grid_x <= 0 || grid_y <= 0 || grid_x > n_w || grid_y > n_h
            continue;
        end
        if grid(grid_x, grid_y) > 0
            continue;
        end
        % Check eight adjacent grids.
        dirs = [-1 -1;-1 0;-1 1;0 -1;0 1;1 -1;1 0;1 1];
        valid = true;
        for ii = 1:8
            new_grid_x = grid_x + dirs(ii, 1);
            new_grid_y = grid_y + dirs(ii, 2);
            if new_grid_x == 0 || new_grid_y == 0 || new_grid_x > n_w || new_grid_y > n_h
                continue;
            end
            if grid(new_grid_x, new_grid_y) == 0
                continue;
            end
            p = points(grid(new_grid_x, new_grid_y), :);
            if sqrt((new_p(1) - p(1)).^2 + (new_p(2) - p(2)).^2) < r
                valid = false;
                break;
            end
        end
        if ~valid
            continue;
        end
        % Add this point.
        n_points_generated = n_points_generated + 1;
        grid(grid_x, grid_y) = n_points_generated;
        active_points = [active_points; n_points_generated];
        points(n_points_generated, :) = new_p;
        success = true;
        break;
    end
    % Remove this point if the generation fails.
    if ~success
        active_points(idx,:) = [];
    end
end
if n_points_generated < n
    points = points(1:n_points_generated, :);
end
end

[tyler]

https://www.mathworks.com/matlabcentral/answers/453169-pseudo-random-distribution-of-points-with-minimum-distance#comment_687604

There is also what appears to be Voronoi iteration algorithm; I am currently trying (and failing) to generalize this to 3D. Not sure which I prefer. Again, just sharing to spread other people's knowledge. The more we know, the more we can help others.

Answer 2

[SS]

Hi,

If I may bother you once again please. I am looking to generate random ellipsoids in a cube, given the volume fraction of ellipsoids (say 40%), non-overlapping (can touch), and they must follow certain particle size distribution. One instance could be, 60% of 2mm, 20% of 4 mm, 15% of 8 mm and 5% of size 16 mm. These sizes (2mm, 4 mm...) are the lengths of the largest semi-principal axes.

Thank you for your time.

S

[Image Analyst]

I suggest you post this in a brand new question. But try the ellipsoid function first using the "try and keep/reject" concept. If you can't get it working, post your code in a new question since it's a different topic than this one.

[SS]

Hi,

Thank you, I have posted it as a separate question. Thank you for the hint, in the meantime, I am also trying.

S

[Image Analyst]

Yeah but I see you didn't take my advice, so expect a lot of questions like "Do the axes of the ellipsoids have to align with the xyz axes?", "What have you tried?", "Have you seen the ellipsoid function?", "Is this a homework problem?", and so on.

[SS]

Hi,

I take your point. Once I formulate the code, I will put it up in the discussions. The idea for posting in advance was to know if someone has already done it.

There should be six degrees of freedom (three rotations and three axes). Indeed I was looking at ellipsoid function, where in I use random function for all the six variables. I am figuring out the best way to determine minimum distance, placing the particles, etc. not yet there. It is not a homework problem, but needed for a small element of my own research work.

Thanks again. S