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.