function M = area_intersect_circle_analytical(varargin)
% Compute the overlap area between 2 circles defined in an array.
% Computation is vectorized, and intersection area are computed an
% analytical way.
%
% Input: Circles data presented in an array G of three columns.
% G contains parameters of the n circles
% . G(1:n,1) - x-coordinate of the center of circles,
% . G(1:n,2) - y-coordinate of the center of circles,
% . G(1:n,3) - radii of the circles
% Each row of the array contains the information for one circle.
%
% Input can also be provided in three different vectors. These
% vectors can be row or column vectors. The 1st one corresponds to
% x-coordinate of the center of circles, the 2nd one to the
% y-cooridnate and the 3rd one to the radii of the circles.
%
%
% Output: Square matrix M(n,n) containing intersection areas between
% circles
% M(i,j) contains the intersection area between circles i & j
% By definition, M(i,i) corresponds to the area of circle i.
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% For 2 circles i & j, three cases are possible depending on the distance
% d(i,j) of the centers of the circles i and j.
% Case 1: Circles i & j do not overlap, there is no overlap area M(i,j)=0;
% Condition: d(i,j)>= ri+rj
% M(i,j) = 0;
% Case 2: Circles i & j fully overlap, the overlap area has to be computed.
% Condition: d(i,j)<= abs(ri-rj)
% M(i,j) = pi*min(ri,rj).^2
% Case 3: Circles i & j partially overlap, the overlap area has to be computed
% decomposing the overlap area.
% Condition: (d(i,j)> abs(ri-rj)) & (d(i,j)<(ri+rj))
% M(i,j) = f(xi,yi,ri,xj,yj,rj)
% = ri^2*arctan2(yk,xk)+ ...
% rj^2*arctan2(yk,d(i,j)-xk)-d(i,j)*yk
% where xk = (ri^2-rj^2+d(i,j)^2)/(2*d(i,j))
% yk = sqrt(ri^2-xk^2)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Guillaume JACQUENOT
% guillaume.jacquenot@gmail.com
% 2007_06_08
% Some explanations, comments were added. Input arguments are detailed.
% Verifications of input arguments are performed.
% The number of calculations has been divided by two.
if nargin==0
% Performs an example
% Creation of 5 circles
x = [0,1,5,3,-5];
y = [0,4,3,7,0];
r = [1,5,3,2,2];
elseif nargin==1
temp = varargin{1};
x = temp(:,1);
y = temp(:,2);
r = temp(:,3);
elseif nargin==3
x = varargin{1};
y = varargin{2};
r = varargin{3};
else
error('The number of arguments must 0, 1 or 3')
end
% Inputs are reshaped in
size_x = numel(x);
size_y = numel(y);
size_r = numel(r);
x = reshape(x,size_x,1);
y = reshape(y,size_y,1);
r = reshape(r,size_r,1);
% Checking if the three input vectors have the same length
if (size_x~=size_y)||(size_x~=size_r)
error('Input of function must be the same length')
end
% Checking if there is any negative or null radius
if any(r<=0)
disp('Circles with null or negative radius won''t be taken into account in the computation.')
temp = (r>0);
x = x(temp);
y = y(temp);
r = r(temp);
end
% Checking the size of the input argument
if size_x==1
M = pi*r.^2;
return
end
% Computation of distance between all circles, which will allow to
% determine which cases to use.
[X,Y] = meshgrid(x,y);
D = sqrt((X-X').^2+(Y-Y').^2);
% Since the resulting matrix M is symmetric M(i,j)=M(j,i), computations are
% performed only on the upper part of the matrix
D = triu(D,1);
[R1,R2] = meshgrid(r);
sumR = triu(R1+R2,1);
difR = triu(abs(R1-R2),1);
% Creating the resulting vector
M = zeros(size_x*size_x,1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Case 2: Circles i & j fully overlap
% One of the circles is inside the other one.
C1 = triu(D<=difR);
M(C1) = pi*min(R1(C1),R2(C1)).^2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Case 3: Circles i & j partially overlap
% Partial intersection between circles i & j
C2 = (D>difR)&(D<sumR);
% Computation of the coordinates of one of the intersection points of the
% circles i & j
Xi(C2,1) = (R1(C2).^2-R2(C2).^2+D(C2).^2)./(2*D(C2));
Yi(C2,1) = sqrt(R1(C2).^2-Xi(C2).^2);
% Computation of the partial intersection area between circles i & j
M(C2,1) = R1(C2).^2.*atan2(Yi(C2,1),Xi(C2,1))+...
R2(C2).^2.*atan2(Yi(C2,1),(D(C2)-Xi(C2,1)))-...
D(C2).*Yi(C2,1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute the area of each circle. Assign the results to the diagonal of M
M(1:size_x+1:size_x*size_x) = pi.*r.^2;
% Conversion of vector M to matrix M
M = reshape(M,size_x,size_x);
% Creating the lower part of the matrix
M = M + tril(M',-1);
% Display results if no argument is provided
if nargin==0
f = figure;
hAxs = axes('Parent',f);
hold on, box on, axis equal
xlabel('x')
ylabel('y','Rotation',0)
title('Compute the intersection area between 2 circles a vectorized way for several cicles')
text(-2,-2,'Numbers on lines represent the intersection')
text(-2,-3,'area between 2 circles')
axis([-8 9 -4 10])
colour = rand(size_x,3);
for t = 1: size_x
plot(x(t)+r(t).*cos(0:2*pi/100:2*pi),...
y(t)+r(t).*sin(0:2*pi/100:2*pi),'color',colour(t,:),'parent',hAxs)
plot(x(t),y(t),'+','color',colour(t,:),'parent',hAxs)
for u = t+1:size_x
if M(t,u)~=0
plot([x(t),x(u)],[y(t),y(u)],'k-','parent',hAxs)
text((x(t)+x(u))/2,(y(t)+y(u))/2,num2str(M(t,u)),...
'HorizontalAlignment','center',...
'BackgroundColor',[1 1 1],'parent',hAxs)
end
end
end
end
% function A=cercle_inter(G)
% % Guillaume JACQUENOT
% % guillaume.jacquenot@gmail.com
% % 2007_05_20
% % % Input: G - array, which contain parameters of the circles
% % G(n,1) - x-coordinate,
% % G(n,2) - y-coordinate,
% % G(n,3) - radius of circle number n=1,2
% % Output: A - area of the overlapping circles
% d=sqrt((G(1,1)-G(2,1))^2+(G(1,2)-G(2,2))^2);
%
% if d>=(G(1,3)+G(2,3))
% % No contact between circles.
% A=0;
% elseif d<=abs(G(1,3)-G(2,3))
% % The smaller circle is inside the big one.
% A=pi*(min(G(1,3),G(2,3)))^2;
% else
% xi=(G(1,3)^2-G(2,3)^2+d^2)/(2*d);
% yi=sqrt(G(1,3)^2-xi^2);
% A=G(1,3)^2*atan2(yi, xi )+...
% G(2,3)^2*atan2(yi,(d-xi))-d*yi;
% end
