## Copyright © 2020  Francesco Potortì
##
## This program is free software; you can redistribute and modify it under
## the terms of the GNU General Public License version 3 or later.

## The problem is overdetermined, because 7 independent constraints are
## given in vd, yet only 6 degrees of freedom are available for the
## relative positions and overall size of three circles.  Using ellipses
## would solve the problem, but computation would be either slower or
## more complex.  As a consequence, the result is not exact.  In most
## practical cases however, the error is small enough to be negligible
## from a graphical point of view.

## Author: Francesco Potortì <Potorti@isti.cnr.it>
## $Revision: 2.13 $
## Usage: [CIRCLES, TOTALAREA, AREAS, PARTITIONS, ERR] = venn3circles (VD)
## Description: Draw a 3-circles Venn diagram with areas proportional to data

## Compute and draw three circles that approximate the input data
##
## VD: [A1 A2 A3 A12 A13 A23 A123] the areas of discs and their intersections
##
## CIRCLES: [c1; c2; c3] where cN is [x y r]
## AREAS: the computed areas in the same format as VD
## TOTALAREA: the total area computed of c1 ∪ c1 ∪ c3
## PARTS: partitions of computed areas: c1,c2,c3 only; c12,c13,c23 only; c123
## ERR: maximum relative error on non-null VD areas
##
function [circles totalarea areas parts err] = venn3circles (vd)
  debug = false;

  ## Compute the areas of the seven partitions of the union of the discs
  ## and check that they are non-negative
  vdp = partitions(vd);			# Areas of VD partitions
  darea = sum(vdp);			# Whole area: c1 ∪ c1 ∪ c3
  assert(all(vdp>=0) && darea > 0, "venn3circles:invalid-input",
	 "venn3circles: no three sets exist with the given intersection areas\n");

  ## Now look for three circles whose intersection areas match VD and whose
  ## centers as as close as possible.
  ## This is the initial point; it is made of three partially overlapping
  ## discs whose total area is then scaled to be the same as VD
  c = [[0 0 2] [2 0 2] [1 2 2]];	# each circle is [x y r]
  carea = 27;				# approximate total area of c
  cs = c'*sqrt(darea/carea);		# c scaled to darea

  ## Use qsp to minimise the relative quadratic error between areas
  ##
  ## Center the first circle on the origin and the second one on the x axis
  ## Position the second circle on the right and the third circle on top
  ## Force radii to be non-negative
  ## Do not optimise on err: optimisation would to be less stable
  ##
  cc = [1 2 4 5 7 8];			# indices of centres coordinates
  zda = any(vd == 0);			# zero-divide avoidance
  perw = 1e-6;				# perimeter weight
  vd /= darea; cs /= sqrt(darea);	# scale for areas around 1 (improves sqp convergence)
  phi = @(x) sumsq((intersection3(x)-vd)./(vd+zda)) + perw*sumsq(x(cc)-x(shift(cc,2)));
  g = @(x) x([1 2 5]) - cs([1 2 5]);	# 3 of 9 degrees of freedom are fixed
  h = @(x) x([4 8 3 6 9]);		# probably useless
  [x, obj, info, iter] = sqp(cs, phi, g, h); # find minimum of phi
  vd *= darea; x *= sqrt(darea);	# rescale back

  if (obj > 0.05 || (info != 101 && info != 104)) # sqp did not converge well
    warning("venn3circles:bad-convergence",
            "difficult or impossible circle configuration\n");
  endif

  if debug x, obj, info, iter endif

  draw = (nargout == 0); # if no expected output, then draw circles

  ## Remove circles with vanishing radius; possibly draw remaining circles
  for ci = 1:3				# iterate over the three circles ci
    cbi = 3*(ci-1);			# base index for circle ci
    if (x(cbi+3) < sqrt(darea)/50)	# if the radius is vanishing
      x(cbi+(1:3)) = 0;			# delete the vanishing circle
    elseif draw
      ## Draw circle
      col = 'rgb'(ci);
      drc(x(cbi+(1:3))/darea, col, col, 1);
    endif
  endfor

  ## Return values
  if !draw
    circles = reshape(x, 3, 3)';
    areas = intersection3(x);
    totalarea = sum(areas);
    parts = partitions(areas);
    ref = vd; ref(vd < sqrt(eps)) = darea;
    nz = (vd > 0);			    # non-null areas
    err = max( abs(areas-vd)(nz)./vd(nz) ); # max relative error on non-null areas
  endif

endfunction


## Compute the areas of the seven partitions of the union of the discs
function vdp = partitions(vd)
  A = [ 1  0  0 -1 -1  0 +1 		# Only c1 data
	0 +1  0 -1  0 -1 +1 		# Only c2 data
	0  0 +1  0 -1 -1 +1 		# Only c3 data
	0  0  0 +1  0  0 -1 		# Only data common to c1 & c2
	0  0  0  0 +1  0 -1 		# Only data common to c1 & c3
	0  0  0  0  0 +1 -1 		# Only data common to c2 & c3
	0  0  0  0  0  0 +1 ];		# Data common to c1,c2,c3
  s = size(vd);
  p = A * vd(:);
  vdp = reshape(p, s);
endfunction


## Draw a circle over the current figure
function h = drc (c, f, e, w=1)
  x = c(1); y = c(2); r = c(3);
  d = r*2;
  px = x-r;
  py = y-r;
  if (isnull(f))
    h = rectangle('Position', [px py d d], 'Curvature',[1,1],
		  'EdgeColor', e, 'LineWidth', w);
  else
    h = rectangle('Position', [px py d d], 'Curvature',[1,1],
		  'FaceColor', f, 'EdgeColor', e, 'LineWidth', w, 'FaceAlpha', 0.5);
  endif
  daspect([1,1,1])
endfunction


## Compute the seven areas of intersection of three discs
##
## AC1, AC2, AC3 the three circles, where ACn is [x y r]
##
## INTEAREAS: [A1 A2 A3 A12 A13 A23 A123] areas of discs and their intersections
## CIRCTRI: [x12 y12 x13 y13 x23 y23] circle intersection points which define
##	the inner circular triangle of area A123
##
## Reference:
## From "Area of Common Overlap of Three Circles" by M.P. Fewell, 2006, we use
## the taxonomy of relative circle positions, which distinguishes 9 cases.  In
## all cases but one the areas are simly computed by summing and subtracting the
## pairwise intersection areas of discs computed with intersection2.  The only
## nontrivial case is case 1, where a convex circular triangle exists.  Most of
## this function is concerned with detecting this case and computing the
## circular triangle. The formulas are taken from the same source, with few
## simplifications and a correction.
##
function [intareas circtri] = intersection3 (ac1, ac2, ac3)
  debug = false;

  if (numel(ac1) == 9)
    ac = reshape(ac1',3,3)'; ac1 = ac(1,:); ac2 = ac(2,:); ac3 = ac(3,:);
  endif

  ## Renumber circles so that r1 ≥ r2 ≥ r3
  [c perm] = sortrows([ac1; ac2; ac3],-3);	# circles sorted by decreasing radius

  ## Change coordinates so that x1=0, y1=0, y2=0, y3≥0
  c(:,1:2) -= c(1,1:2);			# center c1 on the origin
  cp = cart2pol(c(:,1:2));		# transform to polar coordinates
  cp(:,1) -= cp(2,1);			# rotate to center C2 on the x axis
  c(:,1:2) = pol2cart(cp);		# transform back to Cartesian coordinates
  if (c(3,2) < 0)			# if third circle is centered below x axis
    c(:,2) *= -1;			# mirror along the x axis
  endif

  ## Normalised circles
  c1 = c(1,:);				# max radius circle
  c2 = c(2,:);				# mid radius circle
  c3 = c(3,:);				# min radius circle
  if debug subplot(1,2,1); drc(ac1,[],'r'); drc(ac2,[],'g'); drc(ac3,[],'b');
      subplot(1,2,2); drc(c1,[],'r'); drc(c2,[],'g'); drc(c3,[],'b'); endif

  [x1 x2 x3] = deal(    0, c2(1), c3(1));	# x coordinates
  [y1 y2 y3] = deal(    0,     0, c3(2));	# x coordinates
  [r1 r2 r3] = deal(c1(3), c2(3), c3(3));	# radii
  [a1 a2 a3] = deal(pi*r1^2, pi*r2^2, pi*r3^2);	# areas

  ## Compute pairwise intersection areas of discs d2, d2, d3
  a12 = intersection2(c1, c2);		# area of d1 ∩ d2
  a13 = intersection2(c1, c3);		# area of d1 ∩ d3
  a23 = intersection2(c2, c3);		# area of d2 ∩ d3

  ## Assume that the circular triangle does not exist
  circtri = [];				# no intersection points

  ## Compute the intersection areas of the three discs
  ## These are cases 1-7 as defined by Fewell
  if (a12*a13*a23 > 0)			# all discs do intersect pairwise

    ## Distances
    d12 = norm((c1-c2)(1:2));		# distance of centres of c1 and c2
    d13 = norm((c1-c3)(1:2));		# distance of centres of c1 and c3
    d23 = norm((c2-c3)(1:2));		# distance of centres of c2 and c3

    ## Check for cases 4-6 as defined by Fewell
    c3inc1 = (d13 + r3 <= r1);		# circle 3 is inside circle 1
    c3inc2 = (d23 + r3 <= r2);		# circle 3 is inside circle 2
    c2inc1 = (d12 + r2 <= r1);		# circle 2 is inside circle 1
    if (c3inc1 || c3inc2 || c2inc1)	# we are in one of cases 4, 5, 6
      ## One disk is entirely contained inside one other disc
      if (c3inc1 && c3inc2 && c2inc1)
	a123 = a3;			# case 6 as defined by Fewell
	if debug disp("case 6") endif
      elseif (c3inc1 && c3inc2 && !c2inc1)
	a123 = a3;			# case 5 as defined by Fewell
	if debug disp("case 5") endif
      elseif (c3inc1)
	a123 = a23;     		# case 4 as defined by Fewell
	if debug disp("case 4") endif
      elseif (!c3inc2)
	a123 = a23;     		# case 4 as defined by Fewell
	if debug disp("case 4") endif
      elseif (!c2inc1)
	a123 = a13;     		# case 4 as defined by Fewell
	if debug disp("case 4") endif
      else				# (!c3inc1 && c3inc2 && c2inc1)
	assert(false, "venn3circles:internal-error", "internal error");
      endif

    else 				# cases 1, 2, 3, 7 as defined by Fewell
      ## No circle is inside any other one

      ## Step 2: relevant intersection point of circles 1 and 2
      assert (d12*d13*d23 > 0,
	      "venn3circles:internal-error", "internal error");
      x12 = (r1^2-r2^2+d12^2)/2/d12;
      y12 = sqrt(2*d12^2*(r1^2+r2^2)-(r1^2-r2^2)^2-d12^4)/2/d12;
      assert(isreal(y12),
	     "venn3circles:internal-error", "internal error");
      if debug drc([x12 y12 0.2], [], 'r'); endif

      ## Step 3: compute angles
      cos1 = (d12^2+d13^2-d23^2)/2/d12/d13; # angle1 is 213 from centres of circles
      sin1 = sqrt(1-cos1^2);
      cos2 = -(d12^2-d13^2+d23^2)/2/d12/d23; # angle2 is pi-123 from centres of circles
      sin2 = sqrt(1-cos2^2);

      ## Step 5: relevant intersection points involving circle 3
      x113 = (r1^2-r3^2+d13^2)/2/d13;
      y113 = -sqrt(2*d13^2*(r1^2+r3^2)-(r1^2-r3^2)^2-d13^4)/2/d13;
      x13 = x113*cos1 - y113*sin1;
      y13 = x113*sin1 + y113*cos1;
      if debug drc([x13 y13 0.2], [], 'b'); endif
      x223 = (r2^2-r3^2+d23^2)/2/d23;
      y223 = sqrt(2*d23^2*(r2^2+r3^2)-(r2^2-r3^2)^2-d23^4)/2/d23;
      x23 = x223*cos2 - y223*sin2 + d12;
      y23 = x223*sin2 + y223*cos2;
      circtri = [x12 y12 x13 y13 x23 y23];	# vertices of circular triangle
      if debug drc([x23 y23 0.2], [], 'g'); endif

      ## Step 4: detect the existence of a circular triangle (case 1)
      if (norm([x12 y12]-[x3 y3]) < r3 && # intersection is inside disc 3
	  norm([x12 -y12]-[x3 y3]) > r3)  # mirrored intersection is out of disc 3
	## Original Fewell's test was:
	## 	(x12-d13*cos1)^2 + (y12-d13*sin1)^2 < r3^2 &&
	## 	(x12-d13*cos1)^2 + (y12+d13*sin1)^2 > r3^2)

	## We are in case 1 as defined by Fewell
	if debug disp("case 1") endif

	## Step 6: chord lengths of the circular triangle
	l1 = norm([x12 y12] - [x13 y13]);
	l2 = norm([x12 y12] - [x23 y23]);
	l3 = norm([x13 y13] - [x23 y23]);

	## Step 7: whether more than half of disc 3 is inside the circular triangle
	in3 = y3 < y13 + (y23-y13)/(x23-x13) * (x3-x13);
	if debug printf("more than half of disc 3 is %s\n",
			[{"not " ""}{in3+1} "inside the circular triangle"]); endif

	## Compute area of circular triangle
	## We deal with a convex circular triangle, made of a central triangle
	## plus three circular segments
	a123 = sqrt((l1+l2+l3)*(-l1+l2+l3)*(l1-l2+l3)*(l1+l2-l3))/4;
	a123 += r1^2 * asin(l1/2/r1) - l1/4 * sqrt(4*r1^2-l1^2);
	a123 += r2^2 * asin(l2/2/r2) - l2/4 * sqrt(4*r2^2-l2^2);
	if !in3
	  ## Use the same formula as the two previous cicular segments
	  a123 += r3^2 * asin(l3/2/r3) - l3/4 * sqrt(4*r3^2-l3^2);
	else
	  ## Correction: Fewell's paper wrong in Step 7
	  ## In fact, we need the complement to the circular segment with
	  ## respect to the whole area of disc 3
	  a123 += r3^2*pi;
	  a123 -= r3^2 * asin(l3/2/r3) - l3/4 * sqrt(4*r3^2-l3^2);
	endif

      else				# cases 2, 3 and 7 as defined by Fewell

	## In these three cases, either both intersection between a pair of circles
	## lie inside the other disc or both lie outside.  The combination of
	## possibilities can discriminate among the three cases 2, 3 and 7.  It is
	## enough to check only one intersection for each pair of circles.
	##
	i12inc3 = norm([x12 y12] - [x3 y3]) < r3; # intersections of c1&c2 in c3
	i13inc2 = norm([x13 y13] - [x2 y2]) < r2; # intersections of c1&c3 in c2
	i23inc1 = norm([x23 y23] - [x1 y1]) < r1; # intersections of c2&c3 in c1

	if (sum([i12inc3 i13inc2 i23inc1]) == 0)
	  ## A concave circular triangle exists outside of the circles
	  a123 = 0;			# case 7 as defined by Fewell
	  if debug disp("case 7") endif

	elseif (sum([i12inc3 i13inc2 i23inc1]) == 1)
	  ## The area common to all circles is the intersection of the other two
	  ## intersection areas
	  a123 = min([a12 a13 a23]);	# case 3 as defined by Fewell
	  if debug disp("case 3") endif

	elseif (sum([i12inc3 i13inc2 i23inc1]) == 2)
	  ## This is a circular quadrilateral
	  assert(i23inc1, "venn3circles:internal-error", "internal error");

	  if (i12inc3)
	    ## The midsize circle c2 contains the circular quadrilateral
	    a123 = a23 - (a2 - a12);
	  else			  # i13inc2
	    ## The smallest circle c3 contains the circular quadrilateral
	    a123 = a23 - (a3 - a13);	# case 2 as defined by Fewell
	  endif;
	  if debug disp("case 2") endif

	else
	  assert(false, "venn3circles:internal-error", "internal error");
	endif				# distinguish cases 4, 5 and 6
      endif				# distinguish cases 2,3,7 from 1
    endif				# distinguish cases 1,2,3,7 from 4-6

  else

    ## Cases 8 and 9: not all discs have a pairwise intersection
    a123 = 0;
    if debug
      if (sum([a1 a2 a3] == 0) == 1)
	disp("case 8")
      else
	disp("case 9")
      endif
    endif
  endif					# distinguish cases 1-7 from 8,9

  ## Results where circles are sorted from largest (c1) to smallest (c3)
  [sa(1,1) sa(2,2) sa(3,3)] = deal(a1, a2, a3);
  [sa(1,2) sa(1,3) sa(2,3)] = deal(a12, a13, a23);
  ## if debug display(sa) endif

  ## Use the sorting permutation matrix perm to restore the original order
  for isa = 1:3
    for jsa = 1:3
      if (isa <= jsa)
	ia = perm(isa);
	ja = perm(jsa);
	a(min(ia,ja),max(ia,ja)) = sa(isa,jsa);
      endif endfor endfor
  ## if debug display(a) endif

  intareas = [a(1,1) a(2,2) a(3,3) a(1,2) a(1,3) a(2,3) a123];
  assert(all(partitions(intareas) > -sqrt(eps)), "venn3circles:internal-error", "internal error");
endfunction


## Compute the area common to two discs
##
## C1, C2 circles to intersect, where Cn is [x y r]
##
## IA the intersection area of discs
## IC lenght of the intersection chord
##
## From: Circle-Circle Intersection by Eric W. Weisstein from Wolfram MathWorld™
##
function [ia ic] = intersection2 (c1, c2)

  R = max(c1(3), c2(3));		# greatest radius
  r = min(c1(3), c2(3));		# smallest radius
  d = norm((c1-c2)(1:2));		# distance of centres

  ## Corner cases
  if (d >= R+r)				# disjoint discs
    ic = 0; ia = 0;
    return;
  elseif (d <= R-r)			# totally overlapping discs
    ic = 0; ia = pi*r^2;
    return;
  endif

  ## General case of partially overlapping discs
  ic = sqrt((-d+r-R)*(-d-r+R)*(-d+r+R)*(d+r+R))/d;
  d1 = (d^2-r^2+R^2)/2/d;
  d2 = d-d1;
  ia = r^2*acos(d2/r) + R^2*acos(d1/R) - ic*d/2;

endfunction


## Debug utility function
##
## Given three circles, draw them on the current figure and compute the data
## ready in the format for input to venn3circles
function vd = v3cdraw (c1, c2, c3)

  if (numel(c1) == 9)
    c = reshape(c1',3,3)'; c1 = c(1,:); c2 = c(2,:); c3 = c(3,:);
  endif

  ## Draw picture
  clf
  drc(c1,[],'r'); drc(c2,[],'g'); drc(c3,[],'b');

  vd = intersection3(c1, c2, c3);
endfunction


## Debug utility function
##
## Same as intersection3 but using a Monte Carlo approach
## Too slow and too noisy to be used with sqp
## Covers all cases, including those not managed by intersection3
## Good for checking the results of intersection3
function vd = intersection3mc (c1, c2, c3)

  if (numel(c1) == 9)
    c = reshape(c1',3,3)';
  else
    c = reshape([c1(:);c2(:);c3(:)],3,3)';
  endif

  nb = max(c(:,2)+c(:,3));		# north window border
  sb = min(c(:,2)-c(:,3));		# south window border
  eb = max(c(:,1)+c(:,3));		# east window border
  wb = min(c(:,1)-c(:,3));		# west window border

  vd = zeros(1,7);
  niter = 100; sn = 1e7;
  for iter = 1:niter
    rs = rand(sn,2);			# samples in [0,1][0,1]
    s = rs.*[eb-wb,nb-sb]+[wb,sb];	# broadcast to get [x,y] samples in the window
    d = s - permute(c(:,1:2),[3 2 1]);	# broadcast to compute coord diffs from centers
    sed = permute(sumsq(d,2),[1 3 2]);	# squared Euclidean distances from centers
    inc = sed < c(:,3)'.^2;		# broadcast to get if sample is in circle
    vds = [inc inc(:,[1 1 2])&inc(:,[2 3 3]) prod(inc,2)];
    vd += sum(vds);
  endfor
  vd *= (eb-wb)*(nb-sb)/niter/sn;

endfunction



%!shared tol
%! tol = -0.1;			# 10% relative tolerance on error after sqp
%!assert(nthargout (5, @venn3circles, [30 20 10 20 10 10 10])           ,4.16e-08, tol)
%!assert(nthargout (5, @venn3circles, [7 13 20 2 0 4 0])		,1.80e-05, tol)
%!assert(nthargout (5, @venn3circles, [0 20 7 0 0 2 0])			,2.30e-05, tol)
%!assert(nthargout (5, @venn3circles, [78 113 100 3 2 2 0])		,3.86e-04, tol)
%!assert(nthargout (5, @venn3circles, [20 11 12 2 3 4 1])		,1.41e-03, tol)
%!assert(nthargout (5, @venn3circles, [313 200 50 20 45 22 17])		,1.23e-06, tol)
%!assert(nthargout (5, @venn3circles, [4468 6738 3934 798 694 1121 229]),1.32e-02, tol)
%!assert(nthargout (5, @venn3circles, [3852 3766 1681 686 174 263 39])	,2.07e-02, tol)
%!assert(nthargout (5, @venn3circles, [385 377 168 69 17 26 5])		,7.16e-02, tol)
%!warning <difficult or impossible> c = venn3circles([385 377 168 69 16 26 16]);
%!error <no three sets exist> venn3circles([385 377 168 69 16 26 17])

%!demo
%! data = [30 20 10 20 10 10 10
%!         7 13 20 2 0 4 0
%!         0 20 7 0 0 2 0			# empty convex circular triangle
%!         78 113 100 3 2 2 0
%!         20 11 12 2 3 4 1
%!         313 200 50 20 45 22 17		# circular quadrangular
%!         4468 6738 3934 798 694 1121 229	# real-life case
%!         3852 3766 1681 686 174 263 39	# real-life case
%!         385 377 168 69 17 26 5];		# 9% error on areas
%! figure('position', [1 1 1000 750]);
%! for idx = 1:rows(data)
%!   vd = data(idx,:);
%!   subplot(3,3,idx);
%!   venn3circles(vd);
%!   title(sprintf("%5s%5s%5s%5s%5s%5s%5s\n%5g%5g%5g%5g%5g%5g%5g", ...
%!                 "A1", "A2", "A3", "A12", "A13", "A23", "A123", vd), ...
%!         'fontname', "Courier");
%!   r = vd(1:3);                          # given areas of the three circles
%!   r(r==0) = [];                         # radii of non -vanishing circles
%!   axis off;
%!   lstr = mat2cell(num2str(r'),ones(size(r)));
%!   lh = legend(lstr, 'location', 'northwest');
%!   legend boxoff
%! endfor
