枝切法相位解包裹--之放置枝切线

% input:
% residue_charge:残差图，其中正残差为1，负残差为-1，正常点为0
% max_box_radius:平衡残差的最大搜索半径，目前为5，如果这个半径过大，区域将被Branch Cuts所隔离。
% IM_mask:掩膜，有效区域为1，其他为0，是为了枝切线连接而人为设定的边界
% output:
% branch_cuts:包含枝切线的图，正常部分为0，枝切线连接部分为1
function branch_cuts=BranchCuts(residue_charge, max_box_radius, mask)

[rowdim, coldim]=size(residue_charge);

branch_cuts = ~mask;                                %将mask作为枝切线的初始化边界
residue_charge_masked=residue_charge;               %初始化掩膜后的残差电荷图为残差点
residue_charge(logical(~mask))=0;                   %将不在掩膜中的残差点去掉。Remove all residues except those in the mask
cluster_counter=1;                                  %记录每个群中的残差点个数。Keep track of the number of residues in each cluster
satellite_residues=0;                               %记录卫星参差点个数？Keep track of the number of satellite residues accounted for

residue_binary=(residue_charge~=0);                 %逻辑矩阵标注了残差点的位置，即正残差点（+1标记）和负残差点（-1标记）的位置全部用1表明，非残差点不变是0。
residue_balanced=zeros(rowdim, coldim);             %0矩阵，初始时假定所有的残差点都是不平衡的。
[rowres,colres] = find(residue_binary);             %找到残差点的坐标Find the coordinates of the residues
adjacent_residues=zeros(rowdim, coldim);            %定义搜寻窗口中找到的新残差点的位置Defines the positions of additional residues found in the search box
missed_residues=0;                                  %记录有效的残差点的个数？Keep track of the effective number of residues left unbalanced because of

disp('Calculating branch cuts ...');
tic;%开始计时
temp=size(rowres);
for i=1:temp(1)                                     %temp(1)为残差点个数，几行temp(1)或者几列temp(2)都代表有几个
    radius=1;                                       %初始化搜索半径为1。
    r_active=rowres(i);                             %当前活动残差点的行数和列数。Coordinates of the active residue
    c_active=colres(i);
    count_nearby_residues_flag=1;                   %用于说明是否跟踪周围残差点的标记。Flag to indicate whether or not to keep track of the nearby residues
    cluster_counter=1;                              %每次循环重置残差点个数记录Reset the cluster counter
    adjacent_residues = zeros(rowdim, coldim);        %重置邻近的Reset the adjacent residues indicator
    charge_counter=residue_charge_masked(r_active, c_active);         %保存初始残差点电量Store the initial residue charge
    if residue_balanced(r_active, c_active)~=1                        %判断残差点平衡了吗？Has this residue already been balanced?
        while (charge_counter~=0)                                     %一直循环，直到平衡，即Charge_counter=0 Loop until balanced
            
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %This portion of code serves to search the box perimeter，
            %place branch cuts, and keep track of the summed residue charge
            %这段代码用来寻找窗口边缘，布置枝切线，记录累加后的残差点电量。
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            for m=r_active-radius:r_active+radius    %盒子边界像素的坐标Coordinates of the box border pixels  ***I COULD MAKE THIS MORE EFFICIENT***
                for n=c_active-radius:c_active+radius
                    if (abs(m - r_active)==radius || abs(n - c_active)==radius) && charge_counter~=0  %Ensure that only the border pixels are being scrutinised
                        if m<=1 || m>=rowdim || n<=1 || n>=coldim                                      %Is the current pixel on the image border?
                            if m>=rowdim, m=rowdim;
                            end                                              %Make sure that the indices aren't too large for the matrix
                            if n>coldim, n=coldim;
                            end
                            if n<1, n=1;
                            end
                            if m<1, m=1;
                            end
                            branch_cuts=PlaceBranchCutInternal(branch_cuts, r_active, c_active, m, n);%放置枝切线Place a branch cut between the active point and the border
                            cluster_counter=cluster_counter+1;                                      %Keep track of how many residues have been clustered
                            charge_counter=0;                                                       %将电量标记为“平衡”Label the charge as balanced
                            residue_balanced(r_active, c_active)=1;                                 %将中间的残差点标记为平衡。Mark the centre residue as balanced
                        end
                        if mask(m,n)==0  %如果（m,n）该点在掩膜外，即达到图像边缘，连接边缘
                            branch_cuts=PlaceBranchCutInternal(branch_cuts, r_active, c_active, m, n);%Place a branch cut between the active point and the mask border
                            cluster_counter=cluster_counter+1;                                      %Keep track of how many residues have been clustered
                            charge_counter=0;                                                       %Label the charge as balanced
                            residue_balanced(r_active, c_active)=1;                                 %Mark the centre residue as balanced
                        end
                        if residue_binary(m,n)==1                                                   %判断当前点(m,n)是否是残差点 Is the current pixel a residue?
                            if count_nearby_residues_flag==1                                        %记录遇到的残差点？  Are we keeping track of the residues encountered?
                                adjacent_residues(m,n)=1;                                           %标记目前的残差点为近邻Mark the current residue as adjacent
                            end
                            branch_cuts=PlaceBranchCutInternal(branch_cuts, r_active, c_active, m, n);      %放置枝切线 Place a branch cut regardless of the charge_counter value
                            cluster_counter=cluster_counter+1;                                      %记录群中的残差点 Keep track of how many residues have been clustered
                            if residue_balanced(m,n)==0
                               residue_balanced(m,n)=1;                                            %将当前的残差点标记为“平衡” Mark the current residue as balanced
                                charge_counter=charge_counter + residue_charge_masked(m,n);                %更新电量记录 Update the charge counter
                            end
                            if charge_counter==0                                                    %Is the active residue balanced?
                                residue_balanced(r_active, c_active)=1;                             %Mark the centre (active) residue as balanced
                            end
                        end
                    end
                end
            end
            
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %This next portion of code centres the box on the adjacent
            %residues. If the charge is still not balanced after moving
            %through all adjacent residues, increase the box radius and
            %centre the box around the initial
            %residue.如果在窗口移动过所有临近的残差点之后，电量仍然不平衡，那么扩大窗口半径，然后将窗口移到最开始的残差点的位置。
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            if sum(sum(adjacent_residues))==0                           %如果没有邻近的残差点  If there are no adjacent residues:
                radius=radius+1;                                        %扩大窗口半径 Enlarge the box
                r_active=rowres(i);                                     %把窗口放回原来的残差点上 Centre the larger box about the original active residue
                c_active=colres(i);
            else                                                        %如果有临近的残差点  If there are adjacent residues:
                if count_nearby_residues_flag==1                        %Run this bit once per box being searched
                    [r_adjacent,c_adjacent] = find(adjacent_residues);  %找出临近的残差点的坐标  Find the coordinates of the adjacent residues
                    adjacent_size=size(r_adjacent);                     %在周界上有多少个残差点 How many residues are on the perimeter
                    r_active=r_adjacent(1);                             %Centre the search box about a nearby residue
                    c_active=c_adjacent(1);
                    adjacent_residue_count=1;
                    residue_balanced(r_active, c_active)=1;             %Mark the centre (active) residue as balanced before moving on
                    count_nearby_residues_flag=0;                       %不再对周围的残差点计数  Disable further counting of adjacent residues
                else
                    %disp(['Moving block size ', int2str(radius), ' for point ', int2str(i)]);
                    adjacent_residue_count=adjacent_residue_count + 1;  %Move to the next nearby residue
                    if adjacent_residue_count<=adjacent_size(1)
                        r_active=r_adjacent(adjacent_residue_count);    %Centre the search box about the next nearby residue
                        c_active=c_adjacent(adjacent_residue_count);
                    else                                                %Ok, we're done with this box
                        radius=radius+1;                                %Enlarge the box and move on
                        r_active=rowres(i);                             %Centre the larger box about the original active residue
                        c_active=colres(i);
                        adjacent_residues=zeros(rowdim, coldim);        %Reset the adjacent residues indicator
                        count_nearby_residues_flag=1;                   %Enable further counting of adjacent residues
                    end
                end
            end
            
            if radius>max_box_radius                                    %最大半径  Enforce a maximum boundary condition
                %disp('Warning: unreasonably large search area...');
                if cluster_counter~=1                                   %The single satellite residues will still be taken care of
                    missed_residues=missed_residues+1;                  %This effectively means that we have an unbalanced residue
                else
                    satellite_residues=satellite_residues+1;            %This residues is about to be accounted for...
                end
                charge_counter=0;
                
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                %This next portion of code ensures that single satellite
                %residues are not left unaccounted for. The box is simply
                %expanded regardless until the border or ANY other residue
                %is found.
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                while cluster_counter==1                                %If the centre residue is a single satellite, then continue to increase the search radius
                    r_active=rowres(i);                                %Centre the box on the original residue
                    c_active=colres(i);
                    for m=r_active-radius:r_active+radius              %Coordinates of the box border pixels  ***This code works but could definitely be made more efficient***
                        for n=c_active-radius:c_active+radius
                            if (abs(m - r_active)==radius || abs(n - c_active)==radius)                      %Ensure that only the border pixels are being scrutinised
                                if m<=1 || m>=rowdim || n<=1 || n>=coldim                                      %Is the current pixel on the image border?
                                    if m>=rowdim m=rowdim; end                                              %Make sure that the indices aren't too large for the matrix
                                    if n>coldim n=coldim; end
                                    if n<1 n=1; end
                                    if m<1 m=1; end
                                    branch_cuts=PlaceBranchCutInternal(branch_cuts, r_active, c_active, m, n);      %Place a branch cut between the active point and the border
                                    cluster_counter=cluster_counter+1;                                      %Keep track of how many residues have been clustered
                                    residue_balanced(r_active, c_active)=1;                                 %Mark the centre (active) residue as balanced
                                end
                                if mask(m,n)==0                                                          %Does the point fall on the mask
                                    branch_cuts=PlaceBranchCutInternal(branch_cuts, r_active, c_active, m, n);      %Place a branch cut between the active point and the mask border
                                    cluster_counter=cluster_counter+1;                                      %Keep track of how many residues have been clustered
                                    residue_balanced(r_active, c_active)=1;                                 %Mark the centre (active) residue as balanced
                                end
                                if residue_binary(m,n)==1                                                   %Is the current pixel a residue?
                                    branch_cuts=PlaceBranchCutInternal(branch_cuts, r_active, c_active, m, n);      %Place a branch cut regardless of the charge_counter value
                                    cluster_counter=cluster_counter+1;                                      %Keep track of how many residues have been clustered
                                    residue_balanced(r_active, c_active)=1;                                 %Mark the centre (active) residue as balanced
                                end
                            end
                        end
                    end
                    radius=radius+1;                                %Enlarge the box and continue searching
                end
            end
            
        end         % of " while (charge_counter~=0) "
    end             % of " if residue_balanced(r_active, c_active)~=1 "
end

t=toc;
disp(['Branch cut operation completed in ',int2str(t),' seconds.']);
disp(['Unbalanced residues = ', num2str(100*missed_residues/sum(sum(residue_binary))), ' %']);
disp(['Satellite residues accounted for = ', num2str(satellite_residues)]);
% 
% figure;
% mesh(branch_cuts);
% xlabel('X/Pixels','FontSize',14); ylabel('Y/Pixels','FontSize',14);
% title('branch_cuts','FontSize',14);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PlaceBranchCutInternal.m places a branch cut between the points [r1, c1] and
% [r2, c2]. The matrix branch_cuts is binary, with 1's depicting a
% branch cut. 在点[r1, c1]和点[r2, c2]之间放置枝切线.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function branch_cuts=PlaceBranchCutInternal(branch_cuts, r1, c1, r2, c2)

branch_cuts(r1,c1)=1;                       %填充初始点Fill the initial points
branch_cuts(r2,c2)=1;
radius=sqrt((r2-r1)^2 + (c2-c1)^2);         %两点间距离Distance between points
% warning off MATLAB:divideByZero;
m=(c2-c1)/(r2-r1);                          %线的斜率Line gradient
theta=atan(m);                              %线的倾角Line angle

if c1==c2                                   %如果两点列相等Cater for the infinite gradient
    if r2>r1
        for i=1:radius
            r_fill=r1 + i;
            branch_cuts(r_fill, c1)=1; %填充r1至r2之间的所有点
        end
        return;
    else
        for i=1:radius
            r_fill=r2 + i;
            branch_cuts(r_fill, c1)=1; %填充r1至r2之间的所有点
        end
        return;
    end
end

%Use different plot functions for the different quadrants 
%对不同的象限使用不同的绘图函数（有待改进）
if theta<0 && c2<c1                         %
    for i=1:radius                          %Number of points to fill in
        r_fill=r1 + round(i*cos(theta));    %round()四舍五入
        c_fill=c1 + round(i*sin(theta));
        branch_cuts(r_fill, c_fill)=1;
    end
    return;
elseif theta<0 && c2>c1
    for i=1:radius                         %Number of points to fill in
        r_fill=r2 + round(i*cos(theta));
        c_fill=c2 + round(i*sin(theta));
        branch_cuts(r_fill, c_fill)=1;
    end
    return;
end

if theta>0 && c2>c1
    for i=1:radius                          %Number of points to fill in
        r_fill=r1 + round(i*cos(theta));
        c_fill=c1 + round(i*sin(theta));
        branch_cuts(r_fill, c_fill)=1;
    end
    return;
elseif theta>0 && c2<c1
    for i=1:radius                         %Number of points to fill in
        r_fill=r2 + round(i*cos(theta));
        c_fill=c2 + round(i*sin(theta));
        branch_cuts(r_fill, c_fill)=1;
    end
    return;
end