Find multiple circles which have common area of overlap in MATLAB

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

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;