www.pudn.com > matlab-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. % "max_box_radius" defines the maximum search radius for the balancing of % 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 % Last modified on 18/10/2004 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function branch_cuts=BranchCuts(residue_charge, max_box_radius, IM_mask); [rowdim, coldim]=size(residue_charge); branch_cuts=~IM_mask; %Define initial branch cuts borders as the mask. residue_charge_masked=residue_charge; residue_charge(logical(~IM_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); %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 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) %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 IM_mask(m,n)==0 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? 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 IM_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)]); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % 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; 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; end return; else for i=1:radius 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