
function [cents,rs,angles,disk_sequence,arcs_subdomains]=...
    determine_uniondisks_boundary(cents,rs)

%--------------------------------------------------------------------------
% OBJECT.
% This routine given N disks with centers stored in "cents" and radii "rs",
% provides the boundary of the possibly disjoined union of such disks.
% As output of the routine only the disks necessary in the definition of
% the union are given in "cents" and "rs" and later be termed as "active".
% IMPORTANT: The union set must be without "holes".
%--------------------------------------------------------------------------
% INPUT:
% cents: centers of the N disks stored as [x y] in a N x 2 matrix.
% rs: radii of the N disks stored as a N x 1 vector.
%--------------------------------------------------------------------------
% OUTPUT:
% cents: centers of the "active" N disks stored as [x y] in a N x 2 matrix.
% rs: radii of the "active" N disks stored as a N x 1 vector.
% angles: arcs of the disk described in disk_sequence that defines the
%         boundary, as polar coordinated wrt the center in "disk_sequence".
% disk_sequence: what disk defines the i-th arc.
% arcs_subdomains: M x 1 vector such that its i-th component says how many
%         active disks define the i-th connected subregion of the union.
%--------------------------------------------------------------------------
% EXAMPLE:
% >> % OUTPUT OF THE ROUTINE.
%
% centers =
%
%     0.2789    0.9573
%     0.7466    0.6203
%     0.2369    0.6003
%
%
% rs =
%
%     0.1726
%     0.0903
%     0.2553
%
%
% angles =
%
%     1.9119    7.2790
%     5.3082   10.1659
%          0    6.2832
%
%
% disk_sequence =
%
%      3
%      1
%      2
%
%
% arcs_subdomains =
%
%      2
%      1
%      3
%
% >> % the output means that there are 3 disks with centers "cents" and
% radii "rs". The union has arcs that in polar coordinates wrt the centers
% of the pertinent disk as given in "angles". The pertinent disk is given
% in "disk_sequence". The variable "arcs_subdomains" says how many disks in
% each connected region.
%--------------------------------------------------------------------------


[cents,rs]=purge_disk(cents,rs);

[alphas,betas]=determine_intersection_matrices(cents,rs);

[active_arcs,initial_pointer,final_pointer]=determine_active_arcs(alphas,...
    betas);
[arc_sequence,disk_sequence,arcs_sizes]=...
    determine_arc_sequence(active_arcs,initial_pointer,...
    final_pointer,alphas,betas);

[angles,disk_sequence,arcs_subdomains]=reorder_arcs(active_arcs,...
    arc_sequence,disk_sequence,arcs_sizes);









%--------------------------------------------------------------------------
% 1. purge_disk (+ check_indisk, check_ae_disk)
%--------------------------------------------------------------------------

function [centers,rs]=purge_disk(centers,rs)

%--------------------------------------------------------------------------
% Purge some disks that do not contribute to the union of disks being
% strictly inside some other disks.
%--------------------------------------------------------------------------
M=size(centers,1);
purge_list=[];

for ii=1:M
    keep_searching=1;
    center1=centers(ii,:);
    r1=rs(ii);
    for jj=1:M
        center2=centers(jj,:);
        r2=rs(jj);
        if abs(ii-jj) >0 & (keep_searching == 1)
            flag=check_indisk(center1,r1,center2,r2);
            if flag >= 1
                keep_searching=0;
                purge_list=[purge_list; ii];
            end
        end
    end
    
    for jj=ii+1:M
        center2=centers(jj,:);
        r2=rs(jj);
        if abs(ii-jj) >0 & (keep_searching == 1)
            flag=check_ae_disk(center1,r1,center2,r2);
            if flag >= 1
                % warning('\n \t purged almost equal disks.')
                keep_searching=0;
                purge_list=[purge_list; ii];
            end
        end
    end
    
end

ok_centers=setdiff(1:M,purge_list);
centers=centers(ok_centers,:);
rs=rs(ok_centers,:);







%--------------------------------------------------------------------------
% check_indisk
%--------------------------------------------------------------------------

function flag=check_indisk(center1,r1,center2,r2)

%--------------------------------------------------------------------------
% OBJECT
%--------------------------------------------------------------------------
% check if B(center1,r1) is contained in B(center2,r2)
%--------------------------------------------------------------------------

d=norm(center1-center2,2);

if d+r1 <= r2
    flag=1;
else
    flag=0;
end




function flag=check_ae_disk(center1,r1,center2,r2)

% check if B(center1,r1) is contained in B(center2,r2)

d1=norm(center1-center2,2);
d2=abs(r2-r1);

tol=10^(-14);

if d1 <= tol & d2 <= tol
    flag=1;
else
    flag=0;
end








%--------------------------------------------------------------------------
% 2.  determine_intersection_matrices
% 2a. compute_angles
% 2b. decide_alphas_betas
% 2c. decide_arc
%--------------------------------------------------------------------------


%--------------------------------------------------------------------------
% determine_intersection_matrices
%--------------------------------------------------------------------------

function [alphas,betas]=determine_intersection_matrices(cents,rs)

%--------------------------------------------------------------------------
% routine decide_alphas_betas must be fixed.
%--------------------------------------------------------------------------

M=size(cents,1);

alphas=NaN*ones(M,M);
betas=NaN*ones(M,M);

for ii=1:M
    alphas(ii,ii)=0; betas(ii,ii)=2*pi;
end

for ii=1:M
    for jj=ii+1:M
        
        x1=cents(ii,1); y1=cents(ii,2); r1=rs(ii);
        x2=cents(jj,1); y2=cents(jj,2); r2=rs(jj);
        
        C1=[x1 y1]; C2=[x2 y2]; dist_C1_C2=norm(C1-C2);
        tol=10^(-13);
        if abs(dist_C1_C2-(r1+r2)) <= tol
            % warning('Possibly tangent disks (type 1 test).');
            xout =[NaN   NaN]; yout =[NaN   NaN];
        else
            % determine intersection
            [xout,yout] = circcirc(x1,y1,r1,x2,y2,r2);
            
            %         format long e
            %
            %         fprintf('\n \t >> PRE ----------------');
            %         xout
            %         yout
            
            
            % after tests we have noticed that dig=13 may not be enough.
            dig=12;
            spts0=[xout' yout'];
            [spts,ind] = uniquepts(spts0,dig);
            xout=xout(ind); yout=yout(ind);
            
            if norm(size(spts0)-size(spts)) > 0
                %warning('Possibly tangent disks (type 2 test).');
                xout =[NaN   NaN]; yout =[NaN   NaN];
            end
        end
        %         fprintf('\n \t >> POST ----------------');
        %         xout
        %         yout
        %         pause
        points=[xout' yout'];
        
        switch length(xout)
            case 2
                point=points(1,:);
                [angle_loc,radius]=compute_angles([x1 y1],point);
                alphas(ii,jj)=angle_loc;
                
                [angle_loc,radius]=compute_angles([x2 y2],point);
                alphas(jj,ii)=angle_loc;
                
                point=points(2,:);
                [angle_loc,radius]=compute_angles([x1 y1],point);
                betas(ii,jj)=angle_loc;
                
                [angle_loc,radius]=compute_angles([x2 y2],point);
                betas(jj,ii)=angle_loc;
            otherwise
                alphas(ii,jj)=0; betas(ii,jj)=2*pi;
                alphas(jj,ii)=0; betas(jj,ii)=2*pi;
        end
    end
end

[alphas,betas]=decide_alphas_betas(cents,rs,alphas,betas);




%--------------------------------------------------------------------------
% compute_angles
%--------------------------------------------------------------------------

function [angle_pt,radius_pt]=compute_angles(centers,point)

% input:
%
% centers are given by rows: centers=[X Y] with X,Y column vectors.
% point is a row vector.

% angle_pt is in [0,2*pi].

% VECTOR ROUTINE.

% N=size(centers,1);
%
% point_rep=repmat(point,N,1);
%
% ref_point=point_rep-centers;
% X=ref_point(:,1);
% Y=ref_point(:,2);
% [angle_pt,radius_pt] = cart2pol(X,Y);
% angle_pt=rem(angle_pt+2*pi,2*pi);

N=size(centers,1);
ref_point=point-centers;
X=ref_point(:,1);
Y=ref_point(:,2);
[angle_pt,radius_pt] = cart2pol(X,Y);
angle_pt=rem(angle_pt+2*pi,2*pi);







%--------------------------------------------------------------------------
% decide_alphas_betas
%--------------------------------------------------------------------------
function [alphas,betas]=decide_alphas_betas(centers,rs,alphas,betas)

M=size(rs,1);

% plot points
for ii=1:M
    center=centers(ii,:);
    r_ii=rs(ii);
    for jj=1:M
        r_jj=rs(jj);
        if abs(ii-jj) > 0
            angs=[alphas(ii,jj) betas(ii,jj)];
            
            if r_ii == r_jj
                
                % fprintf('\n \t r_ii == r_jj \n')
                
                index=decide_arc(center,r_ii,angs);
                angs=angs(index);
                if angs(1) > angs(2)
                    angs(2)=angs(2)+2*pi;
                end
                alphas(ii,jj)=angs(1);
                betas(ii,jj)=angs(2);
                
            else
                
                % fprintf('\n \t r_ii not = r_jj \n')
                
                if r_ii > r_jj
                    
                    [ang_min,i_min]=min(angs);
                    [ang_max,i_max]=max(angs);
                    width=abs(angs(2)-angs(1));
                    
                    % big disk: from P1 to P2.
                    if width < pi
                        alphas(ii,jj)=ang_max;
                        betas(ii,jj)=ang_min+2*pi;
                        index_P1=i_max;
                        index_P2=i_min;
                    else
                        alphas(ii,jj)=ang_min;
                        betas(ii,jj)=ang_max;
                        index_P1=i_min;
                        index_P2=i_max;
                    end
                    
                    % small disk: from P2 to P1.
                    aa=alphas(jj,ii);
                    bb=betas(jj,ii);
                    
                    if index_P2 == 1
                        ang_P2=aa;
                        ang_P1=bb;
                    else
                        ang_P2=bb;
                        ang_P1=aa;
                    end
                    
                    if ang_P1 < ang_P2
                        ang_P1=ang_P1+2*pi;
                    end
                    
                    alphas(jj,ii)=ang_P2;
                    betas(jj,ii)=ang_P1;
                end
                
            end
        end
        
    end
end




%--------------------------------------------------------------------------
% decide_arc
%--------------------------------------------------------------------------

function index=decide_arc(center,r,angs)

CX=center(1); CY=center(2);
theta1=angs(1); theta2=angs(2);
P1=[CX+r*cos(theta1) CY+r*sin(theta1) 0];
P2=[CX+r*cos(theta2) CY+r*sin(theta2) 0];

A=P2-P1;

B=[CX CY 0]-P1;

C=cross(A,B);
if C(3) > 0
    index=[2 1];
else
    index=[1 2];
end









%--------------------------------------------------------------------------
% 3.  determine_active_arcs
% 3a. determine_active_arcs_row
% 3aa. normalize_arcs
% 3ab. reduce_active_arcs
% 3aba. intersection_1arc_manyarcs
% 3abaa. simple_intersection
%--------------------------------------------------------------------------

%--------------------------------------------------------------------------
% determine_active_arcs
%--------------------------------------------------------------------------
function [active_arcs,initial_pointer,final_pointer]=...
    determine_active_arcs(alphas,betas)

N=size(alphas,2);
active_arcs=[];

final_pointer=0;

for ii=1:N
    
    alpha_row=alphas(ii,:); beta_row=betas(ii,:);
    active_arcs_loc=determine_active_arcs_row(alpha_row,beta_row);
    active_arcs=[active_arcs; active_arcs_loc];
    S=size(active_arcs_loc,1);
    if S > 0
        initial_pointer(ii)=final_pointer(end)+1;
        final_pointer(ii)=size(active_arcs,1);
    else
        initial_pointer(ii)=final_pointer(end);
        final_pointer(ii)=initial_pointer(ii);
    end
    
end






%--------------------------------------------------------------------------
% determine_active_arcs_row
%--------------------------------------------------------------------------
function active_arcs=determine_active_arcs_row(alpha_row,beta_row)

N=length(alpha_row);
active_arcs=[0 2*pi];

for ii=1:N
    alpha_loc=alpha_row(ii); beta_loc=beta_row(ii);
    if (isnan(alpha_loc) == 0) & (isnan(alpha_loc) == 0)
        [alpha_mat,beta_mat]=normalize_arcs(alpha_loc,beta_loc);
        active_arcs=reduce_active_arcs(active_arcs,alpha_mat,beta_mat);
    end
end






%--------------------------------------------------------------------------
% normalize_arcs
%--------------------------------------------------------------------------
function [alpha_mat,beta_mat]=normalize_arcs(alpha,beta)

if (isnan(alpha) == 1) & (isnan(beta) == 1)
    alpha_mat=[];
    beta_mat=[];
end

if (alpha >= 2*pi) & (beta > 2*pi)
    % fprintf('\n \t (alpha >= 2*pi) & (beta > 2*pi)')
    alpha_mat=rem(alpha,2*pi);
    beta_mat=rem(beta,2*pi);
end

if (alpha < 2*pi) & (beta > 2*pi)
    % fprintf('\n \t (alpha < 2*pi) & (beta > 2*pi)')
    alpha_mat=[0 alpha];
    beta_mat=[rem(beta,2*pi); 2*pi];
end

if (alpha < 2*pi) & (beta <= 2*pi)
    alpha_mat=alpha;
    beta_mat=beta;
end






%--------------------------------------------------------------------------
% reduce_active_arcs
%--------------------------------------------------------------------------
function active_arcs_post=reduce_active_arcs(active_arcs,alpha_mat,beta_mat)

L=size(alpha_mat,2);
active_arcs_post=[];

for ii=1:L
    alpha_mat_loc=alpha_mat(ii); beta_mat_loc=beta_mat(ii);
    
    if (isnan(alpha_mat_loc) == 0) & (isnan(beta_mat_loc) == 0)
        
        active_arcs_post_loc=intersection_1arc_manyarcs(active_arcs,...
            alpha_mat_loc,beta_mat_loc);
        
    end
    
    active_arcs_post=[active_arcs_post; active_arcs_post_loc];
end






%--------------------------------------------------------------------------
% intersection_1arc_manyarcs
%--------------------------------------------------------------------------
function active_arcs_post_loc=intersection_1arc_manyarcs(active_arcs,...
    alpha_mat_loc,beta_mat_loc)


M=size(active_arcs,1);

active_arcs_post_loc=[];
for jj=1:M
    alphas=active_arcs(jj,1);
    betas=active_arcs(jj,2);
    if isnan(alphas) == 0
        [alpha_post_loc, beta_post_loc]=simple_intersection(alpha_mat_loc,beta_mat_loc,...
            alphas,betas);
    end
    active_arcs_ij=[alpha_post_loc beta_post_loc];
    active_arcs_post_loc=[active_arcs_post_loc; active_arcs_ij];
end





%--------------------------------------------------------------------------
% simple_intersection
%--------------------------------------------------------------------------
function [a3,b3]=simple_intersection(a1,b1,a2,b2)

a3=max(a1,a2); b3=min(b1,b2);

if a3 > b3
    a3=[]; b3=[];
end









%--------------------------------------------------------------------------
% 4. determine_arc_sequence
% 4a. determine_local_sequence
%--------------------------------------------------------------------------

%--------------------------------------------------------------------------
% determine_arc_sequence
%--------------------------------------------------------------------------

function [arc_sequence,disk_sequence,arcs_sizes]=determine_arc_sequence(active_arcs,initial_pointer,final_pointer,alphas,betas)

L=size(betas,1);
for i=1:L
    for j=1:L
        if betas(i,j) > 2*pi
            betas(i,j)=rem(betas(i,j),2*pi); % normalize betas as in tables
        end
    end
end

T=active_arcs; % table of arcs.
ip=initial_pointer;
fp=final_pointer;
arc_sequence=[];
disk_sequence=[];

N=size(T,1); % number of arcs.
N_done=size(arc_sequence,1);
vN=(1:N)';

arcs_sizes=[];

while N_done < N
    
    arcs_to_do=setdiff(vN,arc_sequence);
    ii=min(arcs_to_do); % initial arc.
    
    [angs_seqs,disk_seqs]=determine_local_sequence(T,ip,fp,ii,alphas,betas);
    arc_sequence=[arc_sequence; angs_seqs];
    disk_sequence=[disk_sequence; disk_seqs];
    
    arcs_sizes=[arcs_sizes; length(angs_seqs)];
    
    N_done=size(arc_sequence,1);
    
end







%--------------------------------------------------------------------------
% determine_local_sequence
%--------------------------------------------------------------------------
function [angs_seqs,disk_seqs]=determine_local_sequence(T,ip,fp,ii,...
    alphas,betas)


angs_seqs=ii; % angs_seqs: sequence of integers, rows of the angles table.

curr_disk_candidates=(find(ii <= fp));
curr_disk=min(curr_disk_candidates);

disk_seqs=curr_disk;

gams=T(ii,2);

gams_hist=[gams];
phis_hist=[];

flag=0;

while flag == 0
    
    betas_ii_row=betas(curr_disk,:);
    
    if gams == 2*pi
        next_disk1=find(betas_ii_row == gams);
        next_disk2=find(betas_ii_row == 0);
        next_disk=[next_disk1; next_disk2];
        if length(next_disk) > 1
            next_disk=setdiff(next_disk,curr_disk);
        else
            gams=0;
        end
    else
        next_disk=find(betas_ii_row == gams);
    end
    
    phis=alphas(next_disk,curr_disk);
    
    
    curr_index_table=(ip(next_disk):fp(next_disk))';
    T_curr=T(curr_index_table,:);
    curr_T_index_loc=find(T_curr(:,1) == phis);
    
    curr_T_index_glob=ip(next_disk)+curr_T_index_loc-1;
    angs_update=curr_T_index_glob;
    gams=T(angs_update,2);
    
    if (curr_T_index_glob == ii)
        flag = 1;
    else
        disk_seqs=[disk_seqs; next_disk];
        angs_seqs=[angs_seqs; angs_update];
        curr_disk=next_disk;
        gams_hist=[gams_hist; gams];
        phis_hist=[phis_hist; phis];
    end
end










%--------------------------------------------------------------------------
% 5. reorder_arcs
%--------------------------------------------------------------------------
function [T_out,disk_sequence_out,arcs_sizes_out]=reorder_arcs(T,arc_sequence,disk_sequence,arcs_sizes)

T=T(arc_sequence,:);

L=length(arcs_sizes);
end_index=0;

T_out=[];
disk_sequence_out=[];
arcs_sizes_out=[];

for ii=1:L
    init_index=end_index+1;
    end_index=sum(arcs_sizes(1:ii));
    
    T_loc=[];
    disk_sequence_loc=[];
    
    jj=init_index;
    
    
    while (jj <= end_index)
        
        if (T(jj,2) == 2*pi) & (jj <= end_index-1)
            if (T(jj,1) == 0)
                % fprintf('\n \t decision 1: regular, unique disk')
                T_loc=[T_loc; T(jj,:)];
                disk_sequence_loc=[disk_sequence_loc; disk_sequence(jj)];
                jj=jj+1;
            else
                if (disk_sequence(jj) == disk_sequence(jj+1))
                    % fprintf('\n \t decision 2: cutting')
                    T_loc=[T_loc; T(jj,1) T(jj+1,2)+2*pi];
                    disk_sequence_loc=[disk_sequence_loc; disk_sequence(jj)];
                    jj=jj+2;
                else
                    % fprintf('\n \t decision 3: regular')
                    T_loc=[T_loc; T(jj,:)];
                    disk_sequence_loc=[disk_sequence_loc; disk_sequence(jj)];
                    jj=jj+1;
                end
            end
        else
            % fprintf('\n \t decision 4: regular, %2.0f',jj)
            T_loc=[T_loc; T(jj,:)];
            disk_sequence_loc=[disk_sequence_loc; disk_sequence(jj)];
            jj=jj+1;
        end
        
    end
    
    if (end_index-init_index > 0)
        % fprintf('\n \t reordering top-bottom see');
        if (T_loc(end,2) == 2*pi) & (T_loc(1,1) == 0) & (disk_sequence_loc(1) == disk_sequence_loc(end))
            % fprintf('\n \t reordering top-bottom');
            aa=T_loc(1,2)+2*pi;
            bb=T_loc(end,1);
            T_loc=[T_loc(2:end-1,:);  bb aa];
            disk_sequence_loc=disk_sequence_loc(2:end);
        end
    end
    
    
    T_out=[T_out; T_loc];
    disk_sequence_out=[disk_sequence_out; disk_sequence_loc ];
    arcs_sizes_out=[arcs_sizes_out; length(disk_sequence_loc)];
    
end



%--------------------------------------------------------------------------
% uniquepts
%--------------------------------------------------------------------------

function [spts,ind] = uniquepts(pts,dig)
 
%--------------------------------------------------------------------------
% Object
%--------------------------------------------------------------------------
% eliminates from an array of points pts those having dig digits
% in common with a subset of other points, keeping only one per subset
%--------------------------------------------------------------------------

atol=100*eps;   
m=10.^(floor(log10(abs(pts))-dig+1));
pts1=round(pts./m).*m;
pts1(find(abs(pts)<atol)) = 0;
[spts,ind]=unique(pts1,'rows');
