Find multiple circles which have common area of overlap in MATLAB

Question

I have two matrices. Each matrix is of dimension 3 * k and represents k circles with each column in the form of [x y r] where (x,y) is the center of the circle and r is it's radius. So one matrix which has 3 circles is in the form of [x1 x2 x3; y1 y2 y3; r1 r2 r3]

I need to find what are the circles which overlap and the overlapped area between the two matrices. Suppose, consider a circle in first matrix. Now consider every other circle in second matrix. Now, I need to find out what circles in the second matrix overlap with the circle considered. I need to do this for every circle in the first matrix. Similarly for every circle in the second matrix with respect to first matrix.

So I need for each circle in the two matrices (There are k1+k2 of them), how much area is overlapping with the other matrix and what are the circles in the other matrix that are overlapping.

Clearly there may be multiple circles which overlap.The dimension k may be different for two matrices. I have 2 extra matrices of each matrix, one sorted by it's x coordinates and other sorted by y coordinates if that helps in computation.The problem is that there are a lot of circles in the matrix and i'm looking for an efficient way to do it. Moreover I would like to extend this this more than two matrices and then an efficient algorithm would greatly improve the execution time.

A preview of the two images (whose corresponding matrices of circles i have with me) is given in the this link: https://www.dropbox.com/s/om5has5uw91dj6p/overlap.jpg

Answer 1

You need two things:

1. Find out what pairs of circles have overlap (distance < sum of radii)
2. compute area of overlap for those pairs

EDITED I modified the code to add some vectorization and trap the case where one circle is entirely inside the other - where the Wolfram formula was breaking down. Demonstrated speed with 1000 x 2000 circles compared.

Here is some code that does these things (with simple example). It seems to be reasonably efficient. Took 1.4 seconds to compare 1000 with 2000 circles (1.66 million overlapping).

% circles [x; y; r]
m1 = rand(3, 1000);
m2 = rand(3, 2000);

tic
% distances between circles:
dx = bsxfun(@minus, m1(1,:), m2(1,:)');
dy = bsxfun(@minus, m1(2,:), m2(2,:)');
sumr = bsxfun(@plus, m1(3,:), m2(3,:)');

% distance between centers:
dist = sqrt(dx.^2 + dy.^2);

% pairs that have overlap:
pairs = find(dist < sumr);

s1 = size(m1,2);
s2 = size(m2,2);

% calculate overlaps
c1 = repmat(1:s1, [s2 1]); % circle from m1 corresponding to value in pairs()
c2 = repmat(1:s2, [s1 1])'; % ditto for m2

o = overlap(m1(3, c1(pairs)), m2(3, c2(pairs)), dist(pairs)');
toc

Put the following in its own file overlap.m in your path: this computes the overlap

function o = overlap(r1, r2, d)
r12 = r1.^2;
r22 = r2.^2;
d2 = d.^2;
a1 = acos((d2 + r12 - r22)./(2.*d.*r1));
a2 = acos((d2 + r22 - r12)./(2.*d.*r2));
o = r12.*a1 + r22.*a2-.5*sqrt((-d+r1+r2)*(d+r1-r2)*(d-r1+r2)*(d+r1+r2));

% fix situation where r1 - r2 > d:
% i.e. one circle fully inside the other
f = find(abs(imag(a1)) + abs(imag(a2)) > 0);
o(f) = pi * min(r1(f),r2(f)).^2;

Answer 2

I thought of using rasterization to approximate the overlap between circles. The idea is create N points to approximate a circular path, then pass those to the poly2mask function. This returns a binary mask image filled where the circle is defined (the algorithm should work with subpixel accuracy). Doing this for all pairs of circles, we can then compute the intersection by just applying logical AND to two masks and count the number of "on" pixels left in the result (or better yet use bwarea for a slightly better estimation). Again this is all an approximation because we are using a discrete grid of pixels...

On the plus side we don't have to worry about mathematical equations. In fact this approach should work for any polygonal shapes where it would be difficult to come up with a closed-form solution for the intersection.

Unfortunately this brute-force method proved to be slower than I expected (no where near as fast as @Floris's code), I am posting my attempt anyway. I borrowed the idea of checking whether the two circles interest from the other code, this should cut down the time significantly if most of the circles you have are far apart.

Just for fun, I added some code to visualize the process!

ANIMATION = true;

% size of image we are working in
%sz = [180 360];
sz = [100 100];

% generated two sets of random circles (specified as columns [x;y;r])
% (I am just trying to create circles wholly visible within the box)
k1 = 100; k2 = 80;
M1 = bsxfun(@times, rand(3,k1), [sz(:)*0.6+sz(:)*0.2;0.1*min(sz)+0.1*min(sz)]);
M2 = bsxfun(@times, rand(3,k2), [sz(:)*0.6+sz(:)*0.2;0.1*min(sz)+0.1*min(sz)]);

% animation
if ANIMATION
    clf
    hImg = imshow(zeros(sz), 'InitialMag','fit');
    hLine(1) = line(NaN, NaN, 'Color','r', 'LineWidth',2);
    hLine(2) = line(NaN, NaN, 'Color','b', 'LineWidth',2);
    axis on
end

% used to approximate circles
num = 50;
t = linspace(0, 2*pi, num);
ct = cos(t); st = sin(t);

% test to find which circles intersect
dist = pdist2(M1(1:2,:)', M2(1:2,:)');
sumr = bsxfun(@plus, M1(3,:)', M2(3,:));
skipIdx = (sumr < dist);

% compute overlap between all pairs
overlap = zeros(k1,k2);
for i=1:k1
    for j=1:k2
        % early skip if circles dont interset
        if skipIdx(i,j), continue, end

        % compute each circle points
        x1 = M1(3,i)*ct + M1(1,i); y1 = M1(3,i)*st + M1(2,i);
        x2 = M2(3,j)*ct + M2(1,j); y2 = M2(3,j)*st + M2(2,j);

        % rasterize circles
        BW1 = poly2mask(x1, y1, sz(1), sz(2));
        BW2 = poly2mask(x2, y2, sz(1), sz(2));

        % compute area of intersection of the two masks
        %overlap(i,j) = sum(BW1(:)&BW2(:));
        overlap(i,j) = bwarea(BW1&BW2);

        if ANIMATION
            % update animation
            set(hImg, 'CData',BW1&BW2)
            set(hLine(1), 'XData',x1, 'YData',y1)
            set(hLine(2), 'XData',x2, 'YData',y2)
            title(sprintf('%g',overlap(i,j)))
            drawnow
        end
    end
end

---

Just to get an idea of how accurate the approximation is, I captured an iteration where one of the circles was completely inside the other. That way we can compare the actual area of intersection against our approximation.

For the image shown below we had:

K>> M1(:,i)
ans =
   58.2106    % x
   24.0996    % y
    8.9387    % r

K>> pi*M1(3,i)^2
ans =
  251.0122

K>> overlap(i,j)
ans =
  252.5000

I think that's close enough :) Using a higher resolution grid of pixels will improve the approximation but will definitely slow it down even more.