www.pudn.com > 2D_phasc_unwrapping_algorithms.zip > BranchCuts.m, change:2008-12-22,size:16911b

```%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% BranchCuts.m generates branch cuts based on the phase residues. This is
% done using the Goldstein method, as described in "Two-dimensional phase
% unwrapping: theory, algorithms and software" by Dennis Ghiglia and
% Mark Pritt.
% "residue_charge" is a matrix wherein positive residues are 1 and
% negative residues are 0.
% residues. If this is too large, areas will be isolated by the branch
% cuts.
% "IM_mask" is a binary matrix. This serves as an artificial border for the
% branch cuts to connect to.
% Created by B.S. Spottiswoode on 15/10/2004
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[rowdim, coldim]=size(residue_charge);

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);                 %Logical matrix indicating the position of the residues
residue_balanced=zeros(rowdim, coldim);             %Initially designate all residues as unbalanced

[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);                                    %Loop through the residues
radius=1;                                       %Set the initial box size
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
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)                                     %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***
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
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                                                   %Is the current pixel a residue?
if count_nearby_residues_flag==1                                        %Are we keeping track of the residues encountered?
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
%centre the box around the initial residue.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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
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)]);
else                                                %Ok, we're done with this box
r_active=rowres(i);                             %Centre the larger box about the original active residue
c_active=colres(i);
count_nearby_residues_flag=1;                   %Enable further counting of adjacent residues
end
end
end

%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***
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
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
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)]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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;
theta=atan(m);                              %Line angle

if c1==c2                                   %Cater for the infinite gradient
if r2>r1
r_fill=r1 + i;
branch_cuts(r_fill, c1)=1;
end
return;
else
r_fill=r2 + i;
branch_cuts(r_fill, c1)=1;
end
return;
end
end

%Use different plot functions for the different quadrants (This is very clumsy,
%I'm sure there is a better way of doing it...)
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

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

```