Random non-overlapping circles in MATLAB

For a recent project I had to generate a stimulus set consisting of random non-overlapping circles in MATLAB.  My solution was to keep a list of the outside of a gradually increasing outer ring of circles, and add new ones, expanding the ring.  This requirement is similar to circle-packing, except I needed to have gaps between the circles.

If you use especially large gaps, or the ratio between your maximum and minimum circle sizes is great, I detected the occasional overlap.  But for the parameters I needed, I tested over 10000 tests of 400 circles, with no overlaps.

Here is how it works:

% Parameters for the circles
setup.min_circle_radius = 15;
setup.max_circle_radius = 20;
setup.min_gap = 10;
setup.max_gap = 20;
setup.n_circles = 400;

First I set up some parameters - minimum and maximum circle radiuses and gaps, plus the number of circles I want to generate.

function [ circles ] = CirclePacking( setup )
    % Start with a circle in the middle
    first_circle = generate_random_radius_circle( setup );
    first_circle.position = [ 0, 0 ];

    % Add second circle at random angle
    second_circle = generate_random_radius_circle( setup );

    gap = randi([setup.min_gap, setup.max_gap]);
    total_distance = first_circle.radius + gap + second_circle.radius;
    [x, y] = pol2cart(2 * pi * rand, total_distance);
    second_circle.position = [x, y];

Next I begin with two circles, one in the centre, and another at a random angle from the first.

function [circle] = generate_random_radius_circle( setup )
    circle.radius = randi([setup.min_circle_radius, setup.max_circle_radius]);
end

I represent circles as structs, first defined solely by radius, but then by position. All circles are first generated using this function which makes a random radius.

Back to the main function:

first_circle.before = 2;
first_circle.after = 2;
second_circle.before = 1;
second_circle.after = 1;

% Setup doubly linked list
first_circle.id = 1;
circles(1) = first_circle; % have to set first circle manually
circles = save_circle(second_circle, circles);

I set up a doubly linked list, that is, each circle keeps track of the circle before and after it. With just two, before and after are the same.

function [circles, circle] = save_circle(circle, circles)
    if (~isfield(circle, 'id'))
        circle.id = size(circles,2) + 1;
    end
    circles(circle.id) = circle;
end

'save_circle' adds or saves a circle to the main circles array. I couldn't work out after a google search how to initialise an empty array of structs, so I didn't waste time doing so.

Back to main function:

current_position = 1;
for i = 1:setup.n_circles
    [circles, current_position] = add_new_circle( setup, circles, current_position );
end

To finish the main function, I just call 'add_new_circle', maintaining a variable 'current_position' which points to a particular circle in the circles array.

This is the end of the main function, most of what happens is going on inside 'add_new_circle', but before I go into it, I will show a few more important helper functions:

function circles_overlap = circles_overlap(circle1, circle2, min_gap)
    dist = edist(circle1.position, circle2.position);

    circles_overlap = dist < circle1.radius + circle2.radius + min_gap;
end

function edist = edist(a, b)
    % euclidean distance
    edist = sqrt((a(1)-b(1))^2+(a(2)-b(2))^2);
end

The above functions are self explanatory, just checking whether circles overlap, and a function for euclidean distance.

function next_circle = next_circle(circle, circles)
    next_circle = circles(circle.after);
end

function previous_circle = previous_circle(circle, circles)
    previous_circle = circles(circle.before);
end

The above functions allow me to get the next and previous circles - this is not necessary really, but I find it's easier on the eyes using these functions.

function new_circle = poisiton_new_circle(new_circle, circle1, circle2, gap1, gap2)

    distance1 = new_circle.radius + gap1 + circle1.radius;
    distance2 = new_circle.radius + gap2 + circle2.radius;
    distance_between = edist(circle1.position, circle2.position);

    % Calculate angle to new circle
    angle = get_law_cosines_angle_a(distance1, distance2, distance_between);

    % Calculate angle of first two circles
    vector_first_to_second = circle2.position - circle1.position;
    angle_first_to_second = atan2(vector_first_to_second(2), vector_first_to_second(1));

    angle_of_new_from_first = mod(angle_first_to_second + angle, 2 * pi);
    [x, y] = pol2cart(angle_of_new_from_first, distance1);
    new_circle.position = [x, y] + circle1.position;
end

function angle = get_law_cosines_angle_a(ab,bc,ca)
    num = (ab * ab) + (ca * ca) - (bc * bc);
    den = 2 * ab * ca;

    angle = acos(num / den);
end

The 'position_new_circle' function puts a circle next to circle1 and circle2, taking into account the gaps passed in too. It always results in the anticlockwise order circle1, new_circle, circle2. It is using the law of cosines to do this - I won't explain the maths here.

function [circles, first_circle_index] = add_new_circle( setup, circles, first_circle_index )
    % circles must contain 2 or more circles

    gap_with_first = randi([setup.min_gap, setup.max_gap]);
    gap_with_second = randi([setup.min_gap, setup.max_gap]);

    first_circle = circles(first_circle_index);
    second_circle = next_circle(first_circle, circles);
    new_circle = generate_random_radius_circle(setup);

    % Check backwards for overlaps
    for i = 1:8;

        new_circle = poisiton_new_circle(new_circle, first_circle, second_circle, gap_with_first, gap_with_second);
        before_circle = previous_circle(first_circle, circles);
        number_to_check = min(size(circles,2) ,8);

        overlaps = false;

        for j=1:number_to_check
            % Check if circle overlaps
            if (before_circle.id == second_circle.id)
                break;
            end

            if (circles_overlap(new_circle, before_circle, setup.min_gap))
                first_circle = before_circle;
                first_circle_index = first_circle.id;
                overlaps = true;
                break;
            end

            before_circle = previous_circle(before_circle, circles);
        end

        if (~overlaps)
            break;
        end
    end

    % Update linked list
    circles = update_list(circles, first_circle, second_circle, new_circle);
end

This is the main part of the algorithm that works out where to put a new circle. It starts by setting the gaps, and taking the key circle, and the circle after that. 'first_circle_index' indexes the key circles, all circles are generated around this circle, until overlaps. Often, a new circle will overlap one of the previous circles in the chain. When this occurs I make the overlapping circle the new key circle, and generate circles around that.

I found that checking 8 circles back was more than enough (on the idea that you can fit roughly 6 circles around another circle - if they are similar size) - but this can be modified to check all circles, at a performance hit.

So after positioning a circle between the key circle and the next one, I check circles before the key (making sure I don't check the same circle). If the circles overlap, I set the key circle to the one that overlaps.

If there are no overlaps, then I update the linked list with that circle.

function circles = update_list(circles, circle1, circle2, new_circle)
    new_circle.before = circle1.id;
    new_circle.after = circle2.id;
    [circles, new_circle] = save_circle(new_circle, circles);

    circle1.after = new_circle.id;
    circles = save_circle(circle1, circles);

    circle2.before = new_circle.id;
    circles = save_circle(circle2, circles);
end

Note this could be massively optimised up - but remembering that 'premature optimization is the root of all evil' - I have not yet needed to so.