Code covered by the BSD License

# venn

by

### Darik (view profile)

30 Nov 2008 (Updated )

Plot 2- or 3- circle area-proportional Venn diagram

venn (varargin)
```function varargout = venn (varargin)

%VENN   Plot 2- or 3- circle area-proportional Venn diagram
%
%  venn(A, I)
%  venn(Z)
%  venn(..., F)
%  venn(..., 'ErrMinMode', MODE)
%  H = venn(...)
%  [H, S] = venn(...)
%  [H, S] = venn(..., 'Plot', 'off')
%  S = venn(..., 'Plot', 'off')
%  [...] = venn(..., P1, V1, P2, V2, ...)
%
%venn(A, I) by itself plots circles with total areas A, and intersection
%area(s) I. For two-circle venn diagrams, A is a two element vector of circle
%areas [c1 c2] and I is a scalar specifying the area of intersection between
%them. For three-circle venn diagrams, A is a three element vector [c1 c2 c3],
%and I is a four element vector [i12 i13 i23 i123], specifiying the
%two-circle intersection areas i12, i13, i23, and the three-circle
%intersection i123.
%
%venn(Z) plots a Venn diagram with zone areas specified by the vector Z.
%For a 2-circle venn diagram, Z is a three element vector [z1 z2 z12]
%For a 3-circle venn, Z is a 7 element vector [z1 z2 z3 z12 z13 z23 z123]
%
%venn(..., F) specifies optional optimization options. VENN uses FMINBND to
%locate optimum pair-wise circle distances, and FMINSEARCH to optimize
%overall three-circle alignment. F is a structure with fields specifying
%optimization options for these functions. F may be a two-element array of
%structures, in which case the first structure is used for FMINBND
%function calls, and the second structure is used for FMINSEARCH function
%calls.
%
%venn(..., 'ErrMinMode', MODE)
%Used for 3-circle venn diagrams only. MODE can be 'TotalError' (default),
%'None', or 'ChowRodgers'. When ErrMinMode is 'None', the positions and
%sizes of the three circles are fixed by their pairwise-intersections,
%which means there may be a large amount of error in the area of the three-
%circle intersection. Specifying ErrMinMode as 'TotalError' attempts to
%minimize the total error in all four intersection zones. The area of the
%three circles are kept constant in proportion to their populations. The
%'ChowRodgers' mode uses the the method proposed by Chow and Rodgers
%[Ref. 1] to draw 'nice' three-circle venn diagrams which appear more
%visually representative of the desired areas, although the actual areas of
%the circles are allowed to deviate from requested values.
%
%H = venn(...) returns a two- or three- element vector to the patches
%representing the circles.
%
%[H, S] = venn(...) returns a structure containing descriptive values
%computed for the requested venn diagram. S is a structure with the
%following fields, where C is the number of circles (N = 2 or 3), Z is
%the number of zones (Z = 3 or 7), and I is the number of intersection
%areas (1 or 4)
%
%
% Position          C*2 array of circle centers
%
% ZoneCentroid      Z*2 array of zone centroids (Can be used for labeling)
%
% CirclePop         C-element vector of supplied circle populations.
%                   (I.e., the 'true' circle areas)
%
% CircleArea        C-element of actual circle areas
%
% CircleAreaError   = (CircleArea-CirclePop)/CirclePop
%
% IntersectPop      I-element vector of supplied intersection populations
%                   (I.e., the 'true' intersection areas)
%
% IntersectArea     I-element vector of actual intersection areas
%
% IntersectError    = (IntersectArea-IntersectPop)/IntersectPop
%
% ZonePop           Z-element vector of supplied zone populations. (I.e.
%                   'true' zone areas
%
% ZoneArea          Z-element vector of actual zone areas.
%
% ZoneAreaError     = (ZoneArea-ZonePop)/ZonePop
%
%
%[H, S] = venn(..., 'Plot', 'off')
%S = venn(..., 'Plot', 'off')
%Returns a structure of computed values, without plotting the diagram. This
%which can be useful when S is used to draw custom venn diagrams or for
%exporting venn diagram data to another application. When Plot is set to off,
%the handles vector H is returned as an empty array. Alternatively, the command
%S = venn(..., 'Plot', 'off) will return only the output structure.
%
%[...] = venn(..., P1, V1, P2, V2, ...)
%Specifies additional patch settings in standard Matlab parameter/value
%pair syntax. Parameters can be any valid patch parameter. Values for patch
%parameters can either be single values, or a cell array of length LENGTH(A),
%in which case each value in the cell array is applied to the corresponding
%circle in A.
%
%Examples
%
%   %Plot a simple 2-circle venn diagram with custom patch properties
%   figure, axis equal, axis off
%   A = [300 200]; I = 150;
%   venn(A,I,'FaceColor',{'r','y'},'FaceAlpha',{1,0.6},'EdgeColor','black')
%
%   %Compare ErrMinModes
%   A = [350 300 275]; I = [100 80 60 40];
%   figure
%   subplot(1,3,1), h1 = venn(A,I,'ErrMinMode','None');
%   axis image,  title ('No 3-Circle Error Minimization')
%   subplot(1,3,2), h2 = venn(A,I,'ErrMinMode','TotalError');
%   axis image,  title ('Total Error Mode')
%   subplot(1,3,3), h3 = venn(A,I,'ErrMinMode','ChowRodgers');
%   axis image, title ('Chow-Rodgers Mode')
%   set([h1 h2], 'FaceAlpha', 0.6)
%
%   %Using the same areas as above, display the error optimization at each
%   iteration. Get the output structure.
%   F = struct('Display', 'iter');
%   [H,S] = venn(A,I,F,'ErrMinMode','ChowRodgers','FaceAlpha', 0.6);
%
%   %Now label each zone
%   for i = 1:7
%       text(S.ZoneCentroid(i,1), S.ZoneCentroid(i,2), ['Zone ' num2str(i)])
%   end
%
%
%Copyright (C) 2008 Darik Gamble, University of Waterloo.
%dgamble@engmail.uwaterloo.ca
%
%References
%1. S Chow and P Rodgers. Extended Abstract: Constructing Area-Proportional
%   Venn and Euler Diagrams with Three Circles. Presented at Euler Diagrams
%   Workshop 2005. Paris. Available online:
%   http://www.cs.kent.ac.uk/pubs/2005/2354/content.pdf
%
%2. S Chow and F Ruskey. Drawing Area-Proportional Venn and Euler Diagrams.
%   Lecture Notes in Computer Science. 2004. 2912: 466-477. Springer-Verlag.
%
%3. MP Fewell. Area of Common Overlap of Three Circles. Australian Government
%   Department of Defence. Defence Technology and Science Organisation. 2006.
%   DSTO-TN-0722. Available online:
%   http://dspace.dsto.defence.gov.au/dspace/bitstream/1947/4551/4/DSTO-TN-0722.PR.pdf

%Variable overview
%   A0, A   Desired and actual circle areas
%               A = [A1 A2] or [A1 A2 A3]
%   I0, I   Desired and actual intersection areas
%               I = I12 or [I12 I13 I23 I123]
%   Z0, Z   Desired and actual zone areas
%               Z = [Z1 Z2 Z12] or [Z1 Z2 Z3 Z12 Z13 Z23 Z123]
%   x, y    Circle centers
%               x = [x1 x2] or [x1 x2 x3]
%               r = [r1 r2] or [r1 r2 r3]
%   d       Pair-wise distances between circles
%               d = d12 or [d12 d13 d23]

%Parse input arguments and preallocate settings
[A0, I0, Z0, nCirc, fminOpts, vennOpts, patchOpts] = parseArgsIn (varargin);
[d, x, y, A, I, Z] = preallocVectors (nCirc);
zoneCentroids = []; %Will only be calculated if needed

r = sqrt(A0/pi);

%Determine distance between first circle pair
d(1) = circPairDist(r(1), r(2), I0(1), fminOpts(1));

%Position of second circle is now known
x(2) = d(1);

%First intersection area
I(1) = areaIntersect2Circ(r(1), r(2), d(1));

if nCirc==3
%Pairwise distances for remaining pairs 1&3 and 2&3
d(2) = circPairDist(r(1), r(3), I0(2), fminOpts(1)); %d13
d(3) = circPairDist(r(2), r(3), I0(3), fminOpts(1)); %d23

%Check triangle inequality
srtD = sort(d);
if ~(srtD(end)<(srtD(1)+srtD(2)))
error('venn:triangleInequality', 'Triangle inequality not satisfied')
end

%Guess the initial position of the third circle using the law of cosines
alpha = acos( (d(1)^2 + d(2)^2 - d(3)^2)  / (2 * d(1) * d(2)) );
x(3) = d(2)*cos(alpha);
y(3) = d(2)*sin(alpha);

%Each pair-wise intersection fixes the distance between each pair
%of circles, so technically there are no degrees of freedom left in
%which to adjust the three-circle intersection. We can either try
%moving the third circle around to minimize the total error, or
%apply Chow-Rodgers

switch vennOpts.ErrMinMode
case 'TotalError'
%Minimize total intersection area error by moving the third circle
pos = fminsearch(@threeCircleAreaError, [x(3) y(3)], fminOpts(2));
x(3) = pos(1);
y(3) = pos(2);
case 'ChowRodgers'
%workspace as a nested fcn
end

%Make sure everything is 'up to date' after optimization
update3CircleData;

end

%Are we supposed to plot?
if vennOpts.Plot
if isempty(vennOpts.Parent)
vennOpts.Parent = gca;
end
hVenn = drawCircles(vennOpts.Parent, x, y, r, patchOpts.Parameters, patchOpts.Values);
else
hVenn = [];
end

%Only determine zone centroids if they're needed
%Needed for output structure
nOut = nargout;
if (nOut==1 && ~vennOpts.Plot) || nOut==2
if nCirc == 2
%Need to calculate new areas
A = A0; %Areas never change for 2-circle venn
Z = calcZoneAreas(2, A, I);
zoneCentroids = zoneCentroids2(d, r, Z);
else
zoneCentroids = zoneCentroids3(x, y, d, r, Z);
end
end

%Figure out output arguments
if nOut==1
if vennOpts.Plot
varargout{1} = hVenn;
else
varargout{1} = getOutputStruct;
end
elseif nOut==2
varargout{1} = hVenn;
varargout{2} = getOutputStruct;
end

function err = threeCircleAreaError (pos)

x3 = pos(1);
y3 = pos(2);

%Calculate distances
d(2) = sqrt(x3^2 + y3^2); %d13
d(3) = sqrt((x3-d(1))^2 + y3^2); %d23

%Calculate intersections
%Note: we're only moving the third circle, so I12 is not changing
I(2:3) = areaIntersect2Circ (r(1:2), r([3 3]), d(2:3)); %I13 and I23
I(4) = areaIntersect3Circ (r, d); %I123

%Replace 0 (no intersection) with infinite error
I(I==0) = Inf;

%Error
err = sum(abs((I-I0)./I0));

end

%Initialize an index matrix to select all 7choose2 zone pairs (21 pairs)
idx = nchoosek(1:7, 2);

%Which zone-zone pairs are considered equal?
%Zones within 10% of each other considered equal
zonePairAreas0 = Z0(idx);

%Percent difference in population between the two members of a pair
ar0 = 2*abs(zonePairAreas0(:,1)-zonePairAreas0(:,2))./sum(zonePairAreas0, 2)*100;
eqPairCutoff = 10;
pairIsEq = ar0<=eqPairCutoff;

%Calculate allowable range for pairs of zones considered unequal
if any(~pairIsEq)
%Sort zone areas
[zUneqAreas0, zUneqAreasSrtIdx] = sort(zonePairAreas0(~pairIsEq,:), 2);

%Make a real index array out of the inconvenient index sort returns
n = sum(~pairIsEq);
zUneqAreasSrtIdx = sub2ind([n,2], [1:n; 1:n]', zUneqAreasSrtIdx);

%rp = (largepopulation/smallpopulation)-1
rp = zUneqAreas0(:,2)./zUneqAreas0(:,1)-1;
rpMin = 1 + 0.3*rp;
rpMax = 1 + 2*rp;
end

%Preallocate zone error vector
zoneErr = zeros(1,21);

%Initialize independent parameters to search over
guessParams = [r(1) x(2) r(2) x(3) y(3) r(3)];

%Search!
pp = fminsearch(@chowRodgersErr, guessParams, fminOpts(2));

[r(1) x(2) r(2) x(3) y(3) r(3)] = deal(pp(1), pp(2), pp(3), pp(4), pp(5), pp(6));

function err = chowRodgersErr (p)

%params = [x2 r2 x3 y3 r3]
[r(1), x(2), r(2), x(3), y(3), r(3)] = deal(p(1), p(2), p(3), p(4), p(5), p(6));

%After changing x2, r2, x3, y3, and r3, update circle areas,
%distances, intersection areas, zone areas
update3CircleData;

if any(pairIsEq)
%For zone pairs considered equal, error is equal to square of the
%distance beyond the cutoff; 0 within cutoff
zAreas = Z(idx(pairIsEq,:));
ar = 2*abs(zAreas(:,1)-zAreas(:,2))./sum(zAreas, 2)*100;
isWithinRange = ar<eqPairCutoff;
ar(isWithinRange) = 0;
ar(~isWithinRange) = ar(~isWithinRange) - eqPairCutoff;

%Amplify error for equal zones with unequal areas
eqZoneUneqAreaErrorGain = 10;
ar(~isWithinRange) = ar(~isWithinRange)*eqZoneUneqAreaErrorGain;

zoneErr(pairIsEq) = ar.^2;
end

if any(~pairIsEq)
%For zone pairs considered unequal, error is equal to square of
%the distance from the allowable range of rp

%rp = (largepopulation/smallpopulation)-1
zUneqPairAreas = Z(idx(~pairIsEq,:));

%Sort based on the population sizes (determined by parent
zUneqPairAreas = zUneqPairAreas(zUneqAreasSrtIdx);
rp = zUneqPairAreas(:,2)./zUneqPairAreas(:,1)-1;

lessThanMin = rp<rpMin;
moreThanMax = rp>rpMax;
rp(~lessThanMin & ~moreThanMax) = 0;

%Determine how far out of range errors are
rp(lessThanMin) = rp(lessThanMin) - rpMin(lessThanMin);
rp(moreThanMax) = rp(moreThanMax) - rpMax(moreThanMax);

%Consider the case where rp < rpMin to be more
%erroneous than the case where rp > rpMax
tooSmallErrorGain = 10;
rp(lessThanMin) = rp(lessThanMin)*tooSmallErrorGain;

zoneErr(~pairIsEq) = rp.^2;
end

%Total error
err = sum(zoneErr);

end %chowRodgersErr

function update3CircleData

%Circle areas
A = pi*r.^2;

%Calculate distances
d(1) = abs(x(2)); %d12
d(2) = sqrt(x(3)^2 + y(3)^2); %d13
d(3) = sqrt((x(3)-d(1))^2 + y(3)^2); %d23

%Calculate actual intersection areas
I(1:3) = areaIntersect2Circ (r([1 1 2]), r([2 3 3]), d); %I12, I13, I23
I(4) = areaIntersect3Circ (r, d); %I123

%Calculate actual zone areas
Z = calcZoneAreas(3, A, I);

end

function S = getOutputStruct

S = struct(...
'Position'              ,[x' y']                ,...
'ZoneCentroid'          ,zoneCentroids          ,...
'CirclePop'             ,A0                     ,...
'CircleArea'            ,A                      ,...
'CircleAreaError'       ,(A-A0)./A0             ,...
'IntersectPop'          ,I0                     ,...
'IntersectArea'         ,I                      ,...
'IntersectError'        ,(I-I0)./I0             ,...
'ZonePop'               ,Z0                     ,...
'ZoneArea'              ,Z                      ,...
'ZoneAreaError'         ,(Z-Z0)./Z0             );
end

end %venn

function D = circPairDist (rA, rB, I, opts)
%Returns an estimate of the distance between two circles with radii rA and
%rB with area of intersection I
%opts is a structure of FMINBND search options
D = fminbnd(@areadiff, 0, rA+rB, opts);
intersectArea = areaIntersect2Circ (rA, rB, d);
dA = abs(I-intersectArea)/I;
end
end

function hCirc = drawCircles(hParent, xc, yc, r, P, V)

hAx = ancestor(hParent, 'axes');
nextplot = get(hAx, 'NextPlot');

%P and V are cell arrays of patch parameter/values
xc = xc(:); yc = yc(:);     %Circle centers
n = length(r);

%Independent parameter
dt = 0.05;
t = 0:dt:2*pi;

%Origin centered circle coordinates
X = r*cos(t);
Y = r*sin(t);

hCirc = zeros(1,n);
c = {'r', 'g', 'b'};                        %default colors
fa = {0.6, 0.6, 0.6};                         %default face alpha
tag = {'Circle1', 'Circle2', 'Circle3'}; 	%default tag

for i = 1:n
xx = X(i,:)+xc(i);
yy = Y(i,:)+yc(i);
hCirc(i) = patch (xx, yy, c{i}, 'FaceAlpha', fa{i}, 'Parent', hParent, 'Tag', tag{i});
if i==1
end
end
set(hAx, 'NextPlot', nextplot);

%Custom patch parameter values
if ~isempty(P)

c = cellfun(@iscell, V);

%Scalar parameter values -- apply to all circles
if any(~c)
set(hCirc, {P{~c}}, {V{~c}});
end

%Parameters values with one value per circle
if any(c)
%Make sure all vals are column cell arrays
V = cellfun(@(val) (val(:)), V(c), 'UniformOutput', false);
set(hCirc, {P{c}}, [V{:}])
end
end

end %plotCircles

function A = areaIntersect2Circ (r1, r2, d)
%Area of Intersection of 2 Circles
%Taken from [2]

alpha = 2*acos( (d.^2 + r1.^2 - r2.^2)./(2*r1.*d) );
beta  = 2*acos( (d.^2 + r2.^2 - r1.^2)./(2*r2.*d) );

A =    0.5 * r1.^2 .* (alpha - sin(alpha)) ...
+ 0.5 * r2.^2 .* (beta - sin(beta));

end

function [A, x, y, c, trngArea] = areaIntersect3Circ (r, d)
%Area of common intersection of three circles
%This algorithm is taken from [3].
%   Symbol    Meaning
%     T         theta
%     p         prime
%     pp        double prime

%[r1 r2 r3] = deal(r(1), r(2), r(3));
%[d12 d13 d23] = deal(d(1), d(2), d(3));

%Intersection points
[x,y,sinTp,cosTp] = intersect3C (r,d);

if any(isnan(x)) || any(isnan(y))
A = 0;
%No three circle intersection
return
end

%Step 6. Use the coordinates of the intersection points to calculate the chord lengths c1,
%c2, c3:
i1 = [1 1 2];
i2 = [2 3 3];
c = sqrt((x(i1)-x(i2)).^2 + (y(i1)-y(i2)).^2)';

%Step 7: Check whether more than half of circle 3 is included in the circular triangle, so
%as to choose the correct expression for the area
lhs = d(2) * sinTp;
rhs = y(2) + (y(3) - y(2))/(x(3) - x(2))*(d(2)*cosTp - x(2));
if lhs < rhs
sign = [-1 -1 1];
else
sign = [-1 -1 -1];
end

%Calculate the area of the three circular segments.
ca = r.^2.*asin(c/2./r) + sign.*c/4.*sqrt(4*r.^2 - c.^2);

trngArea = 1/4 * sqrt( (c(1)+c(2)+c(3))*(c(2)+c(3)-c(1))*(c(1)+c(3)-c(2))*(c(1)+c(2)-c(3)) );
A = trngArea + sum(ca);

end

function [x, y, sinTp, cosTp] = intersect3C (r, d)
%Calculate the points of intersection of three circles

%d = [d12 d13 d23]
%x = [x12; x13; x23]
%y = [y12; y13; y23]

%   Symbol    Meaning
%     T         theta
%     p         prime
%     pp        double prime

x = zeros(3,1);
y = zeros(3,1);

%Step 1. Check whether circles 1 and 2 intersect by testing d(1)
if ~( ((r(1)-r(2))<d(1)) && (d(1)<(r(1)+r(2))) )
%x = NaN; y = NaN;
%bigfix: no returned values for sinTp, cosTp
[x, y, sinTp, cosTp] = deal(NaN);
return
end

%Step 2. Calculate the coordinates of the relevant intersection point of circles 1 and 2:
x(1) = (r(1)^2 - r(2)^2 + d(1)^2)/(2*d(1));
y(1) = 0.5/d(1) * sqrt( 2*d(1)^2*(r(1)^2 + r(2)^2) - (r(1)^2 - r(2)^2)^2 - d(1)^4 );

%Step 3. Calculate the values of the sines and cosines of the angles tp and tpp:
cosTp  =  (d(1)^2 + d(2)^2 - d(3)^2) / (2 * d(1) * d(2));
cosTpp = -(d(1)^2 + d(3)^2 - d(2)^2) / (2 * d(1) * d(3));
sinTp  =  (sqrt(1 - cosTp^2));
sinTpp =  (sqrt(1 - cosTpp^2));

%Step 4. Check that circle 3 is placed so as to form a circular triangle.
cond1 = (x(1) - d(2)*cosTp)^2 + (y(1) - d(2)*sinTp)^2 < r(3)^2;
cond2 = (x(1) - d(2)*cosTp)^2 + (y(1) + d(2)*sinTp)^2 > r(3)^2;
if  ~(cond1 && cond2)
x = NaN; y = NaN;
return
end

%Step 5: Calculate the values of the coordinates of the relevant intersection points involving
%circle 3
xp13  =  (r(1)^2 - r(3)^2 + d(2)^2) / (2 * d(2));
%yp13  = -0.5 / d(2) * sqrt( 2 * d(2)^2 * (r(2)^2 + r(3)^2) - (r(1)^2 - r(3)^2)^2 - d(2)^4 );
yp13  = -0.5 / d(2) * sqrt( 2 * d(2)^2 * (r(1)^2 + r(3)^2) - (r(1)^2 - r(3)^2)^2 - d(2)^4 );

x(2)   =  xp13*cosTp - yp13*sinTp;
y(2)   =  xp13*sinTp + yp13*cosTp;

xpp23 =  (r(2)^2 - r(3)^2 + d(3)^2) / (2 * d(3));
ypp23 =  0.5 / d(3) * sqrt( 2 * d(3)^2 * (r(2)^2 + r(3)^2) - (r(2)^2 - r(3)^2)^2 - d(3)^4 );

x(3) = xpp23*cosTpp - ypp23*sinTpp + d(1);
y(3) = xpp23*sinTpp + ypp23*cosTpp;

end

function z = calcZoneAreas(nCircles, a, i)

%Uses simple set addition and subtraction to calculate the zone areas
%with circle areas a and intersection areas i

if nCircles==2
%a = [A1 A2]
%i = I12
%z = [A1-I12, A2-I12, I12]
z = [a(1)-i, a(2)-i, i];
elseif nCircles==3
%a = [A1  A2  A3]
%i = [I12 I13 I23 I123]
%z = [A1-I12-I13+I123, A2-I12-I23+I123, A3-I13-I23+I123, ...
%     I12-I123, I13-I123, I23-I123, I123];
z = [a(1)-i(1)-i(2)+i(4), a(2)-i(1)-i(3)+i(4), a(3)-i(2)-i(3)+i(4), ...
i(1)-i(4), i(2)-i(4), i(3)-i(4), i(4)];
else
error('')
%This error gets caught earlier in the stack w. better error msgs
end
end

function [Cx, Cy, aiz] = centroid2CI (x, y, r)

%Finds the centroid of the area of intersection of two circles.
%Vectorized to find centroids for multiple circle pairs
%x, y, and r are nCirclePairs*2 arrays
%Cx and Cy are nCirclePairs*1 vectors

%Centroid of the area of intersection of two circles
n = size(x,1);
xic = zeros(n,2);
az = zeros(n,2);

dx = x(:,2)-x(:,1);
dy = y(:,2)-y(:,1);
d = sqrt(dx.^2 + dy.^2);

%Translate the circles so the first is at (0,0) and the second is at (0,d)
%By symmetry, all centroids are located on the x-axis.
%The two circles intersect at (xp, yp) and (xp, -yp)
xp = 0.5*(r(:,1).^2 - r(:,2).^2 + d.^2)./d;

%Split the inner zone in two
%Right side (Area enclosed by circle 1 and the line (xp,yp) (xp,-yp)
%Angle (xp,yp) (X1,Y1) (xp,-yp)
alpha = 2*acos(xp./r(:,1));
%Area and centroid of the right side of the inner zone
[xic(:,1) az(:,1)] = circleChordVals (r(:,1), alpha);
%Angle (xp,yp) (X2,Y2) (xp,-yp)
alpha = 2*acos((d-xp)./r(:,2));
%Area and centroid of the left side of the inner zone
[xic(:,2) az(:,2)] = circleChordVals (r(:,2), alpha);
xic(:,2) = d - xic(:,2);

%Thus the overall centroid  & area of the inner zone
aiz = sum(az,2);
Cx = sum(az.*xic,2)./aiz;

%Now translate the centroid back based on the original positions of the
%circles
theta = atan2(dy, dx);
Cy = Cx.*sin(theta) + y(:,1);
Cx = Cx.*cos(theta) + x(:,1);

end

function centroidPos = zoneCentroids2 (d, r, Z)

centroidPos = zeros(3,2);

%Find the centroids of the three zones in a 2-circle venn diagram
%By symmetry, all centroids are located on the x-axis.
%First, find the x-location of the middle (intersection) zone centroid

%Centroid of the inner zone
centroidPos(3,1) = centroid2CI([0 d], [0 0], r);

%Now, the centroid of the left-most zone is equal to the centroid of
%the first circle (0,0) minus the centroid of the inner zone
centroidPos(1,1) = -centroidPos(3,1)*Z(3)/Z(1);

%Similarly for the right-most zone; the second circle has centroid at x=d
centroidPos(2,1) = (d*(Z(2)+Z(3)) - centroidPos(3,1)*Z(3))/Z(2);

end

function centroidPos = zoneCentroids3 (x0, y0, d, r, Z)

Z = Z(:);

%Get area, points of intersection, and chord lengths
[act, xi, yi, c, atr] = areaIntersect3Circ (r, d);
atr = atr(:);
r = r(:);

%Area and centroid of the triangle within the circular triangle is
xtr = sum(xi/3);
ytr = sum(yi/3);

%Now find the centroids of the three segments surrounding the triangle
i = [1 2; 1 3; 2 3];
xi = xi(i); yi = yi(i);
[xcs, ycs, acs] = circSegProps (r(:), x0(:), y0(:), xi, yi, c(:));

%Overall centroid of the circular triangle
xct = (xtr*atr + sum(xcs.*acs))/act;
yct = (ytr*atr + sum(ycs.*acs))/act;

%Now calculate the centroids of the three two-pair intersection zones
%(Zones 12 13 23)
%Entire zone centroid/areas

%x, y, and r are nCirclePairs*2 arrays
%Cx and Cy are nCirclePairs*1 vectors
i = [1 2; 1 3; 2 3];
[x2c, y2c, a2c] = centroid2CI (x0(i), y0(i), r(i));

%Minus the three-circle intersection zone
xZI2C = (x2c.*a2c - xct*act)./(a2c-act);
yZI2C = (y2c.*a2c - yct*act)./(a2c-act);

x0 = x0(:);
y0 = y0(:);

%Finally, the centroids of the three circles minus the intersection
%areas
i1 = [4 4 5]; i2 = [5 6 6];
j1 = [1 1 2]; j2 = [2 3 3];
x1C = (x0*pi.*r.^2 - xZI2C(j1).*Z(i1) - xZI2C(j2).*Z(i2) - xct*act)./Z(1:3);
y1C = (y0*pi.*r.^2 - yZI2C(j1).*Z(i1) - yZI2C(j2).*Z(i2) - yct*act)./Z(1:3);

%Combine and return
centroidPos = [x1C y1C; xZI2C yZI2C; xct yct];
end

function [x, a] = circleChordVals (r, alpha)
%For a circle centered at (0,0), with angle alpha from the x-axis to the
%intersection of the circle to a vertical chord, find the x-centroid and
%area of the region enclosed between the chord and the edge of the circle
a = r.^2/2.*(alpha-sin(alpha));                         %Area
x = 4.*r/3 .* sin(alpha/2).^3 ./ (alpha-sin(alpha));    %Centroid
end

function [xc, yc, area] = circSegProps (r, x0, y0, x, y, c)

%Translate circle to (0,0)
x = x-[x0 x0];
y = y-[y0 y0];

%Angle subtended by chord
alpha = 2*asin(0.5*c./r);

area = r.^2/2.*(alpha-sin(alpha));                         %Area
d   = 4.*r/3 .* sin(alpha/2).^3 ./ (alpha-sin(alpha));    %Centroid

%Perpindicular bisector of the chord
m = -(x(:,2)-x(:,1))./(y(:,2)-y(:,1));

%angle of bisector
theta = atan(m);

%centroids
xc = d.*cos(theta);
yc = d.*sin(theta);

%Make sure we're on the correct side
%Point of intersection of the perp. bisector and the circle perimeter
xb = (x(:,1)+x(:,2))/2;
xc(xb<0) = xc(xb<0)*-1;
yc(xb<0) = yc(xb<0)*-1;

%Translate back
xc = xc + x0;
yc = yc + y0;
end

function [A0, I0, Z0, nCircles, fminOpts, vennOpts, patchOpts] = parseArgsIn (args)

[A0, I0, Z0] = deal([]);
nIn = length(args);

%Get the easy cases out of the way
if nIn == 0
elseif nIn == 1
%venn(Z)
Z0 = args{1};
nIn = 0;
elseif nIn == 2
if isnumeric(args{2})
%venn (A,I)
[A0, I0] = deal(args{1:2});
nIn = 0;
else
%venn (Z, F)
Z0 = args{1};
args = args(2);
nIn = 1;
end
else
%Find the first non-numeric input arg
i = find(~cellfun(@isnumeric, args), 1);
if i == 2
%venn(Z, ....)
Z0 = args{1};
elseif i == 3
%venn(A, I, ...)
[A0, I0] = deal(args{1:2});
else
end
nIn = nIn - i + 1;
args = args(i:end);
end

error('venn:parseInputArgs:unrecognizedSyntax', 'Unrecogized input syntax')
end
try
[A0, I0, Z0] = parseInputAreas (A0, I0, Z0);
catch
error('venn:parseArgsIn:parseInputAreas', 'Incorrect size(s) for area vector(s)')
end
nCircles = length(A0);
nZones = length(Z0);

%Any arguments left?
if nIn > 0

if isstruct(args{1})
%FMIN search options
f = args{1};

nIn = nIn - 1;
if nIn>0, args = args(2:end); end

if length(f) == 1
%Just double up
fminOpts = [f f];
elseif length(f) == 2
%ok
fminOpts = f;
else
error('venn:parseArgsIn', 'FMINOPTS must be a 1 or 2 element structure array.')
end
else
%Use defaults
fminOpts = [optimset('fminbnd'), optimset('fminsearch')];
end
else
%Use defaults
fminOpts = [optimset('fminbnd'), optimset('fminsearch')];
end

%If there's an even number of args in remaining
if nIn>0
if mod(nIn, 2)==0
%Parameter/Value pairs
p = args(1:2:end);
v = args(2:2:end);
[vennOpts, patchOpts] = parsePVPairs (p, v, nZones);
else
error('venn:parseArgsIn', 'Parameter/Value options must come in pairs')
end
else
vennOpts = defaultVennOptions;
patchOpts = struct('Parameters', [], 'Values', []);
end

end %parseArgsIn

function [vennOpts, patchOpts] = parsePVPairs (p, v, nZones)

p = lower(p);

%Break up P/V list into Venn parameters and patch parameters
vennParamNames = {'plot', 'errminmode', 'parent'};
[isVennParam, idx] = ismember(p, vennParamNames);
idx = idx(isVennParam);
%vennParams = p(isVennParam);
vennVals = v(isVennParam);

%First do Patch options
patchOpts.Parameters = p(~isVennParam);
patchOpts.Values = v(~isVennParam);

%Now do Venn options
vennOpts = defaultVennOptions;

%PLOT
i = find(idx==1, 1);
if i
plot = lower(vennVals{i});
if islogical(plot)
vennOpts.Plot = plot;
else
if ischar(plot) && any(strcmp(plot, {'on', 'off'}))
vennOpts.Plot = strcmp(plot, 'on');
else
error('venn:parsePVPairs', 'Plot must be ''on'', ''off'', or a logical value.')
end
end
end

%ERRMINMODE
i = find(idx==2, 1);
if i
mode = lower(vennVals{i});
okModes = {'None', 'TotalError', 'ChowRodgers'};
[isOkMode, modeIdx] = ismember(mode, lower(okModes));
if isOkMode
vennOpts.ErrMinMode = okModes{modeIdx};
else
error('venn:parsePVPairs', 'ErrMinMode must be None, TotalError, or ChowRodgers')
end
end

%PARENT
i = find(idx==5, 1);
if i
h = v{i};
if length(h)==1 && ishandle(h)
vennOpts.Parent = h;
else
error('venn:parsePVPairs', 'Parent must be a valid scalar handle')
end
end

end %parsePVPairs

function [A0, I0, Z0] = parseInputAreas (A0, I0, Z0)

%Switch to row vectors
A0 = A0(:)';
I0 = I0(:)';
Z0 = Z0(:)';

if isempty(Z0)
%A0 and I0 supplied

Z0 = calcZoneAreas (length(A0), A0, I0);
else
%Z0 supplied
switch length(Z0)
case 3
A0 = Z0(1:2)+Z0(3);
I0 = Z0(3);
case 7
A0 = Z0(1:3)+Z0([4 4 5])+Z0([5 6 6])+Z0(7);
I0 = [Z0(4:6)+Z0(7) Z0(7)];
otherwise
error('')
end
end
end

function vennOpts = defaultVennOptions

vennOpts = struct(...
'Plot'          ,true               ,...
'Labels'        ,[]                 ,...
'PopLabels'     ,false              ,...
'DrawLabels'    ,false              ,...
'Parent'        ,[]                 ,...
'Offset'        ,[0 0]              ,...
'ErrMinMode'    ,'TotalError'       );

end

function [d, x, y, A, I, Z] = preallocVectors (nCirc)

%Initialize position vectors
x = zeros(1, nCirc);
y = zeros(1, nCirc);

if nCirc==2
d = 0;
I = 0;
A = zeros(1,2);
Z = zeros(1,3);

else %nCirc==3
d = zeros(1,3);
I = zeros(1,4);
A = zeros(1,3);
Z = zeros(1,7);
end
end

```