MATALB相位解包裹代码

% BranchCuts.m % % Aaron James Lemmer % November 12, 2013 % % Modifications: % November 29, 2013: Logic was added to check if a residue found in the % search of the box perimeter has already been connected to the image % edge. This prevents additional cuts to the edge being placed that % may isolate regions of the image unnecessarily. % % December 12, 2013: Added capability to account for the image % border/mask. % % April 25, 2014: branch_cuts is now initialized in a way to accomodate % dipole pre-processing. function [branch_cuts] = BranchCuts(branch_cuts, residues, num_residues, border) branch_cuts = branch_cuts + border; %matrix storing the branch cuts, initially contains only dipole branch cuts from pre-processing step [row_coord_res, col_coord_res] = find(residues ~= 0); %find the residue coordinates %This step essentially covers the role of the outer two for loops (in the x %and y coordinates of the image) in Goldstein's C implementation. The for %loop method covers every pixel of the image, alternating between locating %and then balancing each residue. Here, all residues are located first, %then all are balanced sequentially. if (isempty(row_coord_res)) %if there are no residues, display message %change this to fprintf return; %invoke an early return; no residues, so no branch cuts required end %if (isempty(row_coord_res)) if (num_residues ~= length(row_coord_res)) disp(['BranchCuts.m: Residue counts contradict!; num_residues =', int2str(num_residues),... '; length(row_coord_res) =', int2str(length(row_coord_res))]) error('BranchCuts:ResCount','Residue counts contradict.'); %invoke an early return; error in counting residues end %if (num_residues ~= length(row_coord_res)) %allocate memory [num_row, num_col] = size(residues); %same as the size of wrapped_phase residues_balanced = zeros(num_row, num_col); %initially, all residues unbalanced (false = 0) residues_active = zeros(num_row, num_col); %initially, all residues inactive (false = 0) connected2edge = zeros(num_row, num_col); %initially, no residues connected to border (false = 0) max_searchbox_rad = floor(length(residues)/2); for current_residue = 1:num_residues %loop through all the (initially) unbalanced residues % disp(['current_residue = ', int2str(current_residue)]) row_current = row_coord_res(current_residue); %row coordinate of the current residue col_current = col_coord_res(current_residue); %column coordinate of the current residue %disp(['(row_current, col_current) = (', int2str(row_current), ',', int2str(col_current), ');']) %evaluate the box center if (residues_balanced(row_current, col_current) > 0) %if current residue is already balanced (possibly by a previous residue) continue; %pass to next iteration of for loop; i.e., go to next unbalanced residue end %if active residue is already balanced if (row_current <= 1 || row_current >= num_row || col_current <= 1 || col_current >= num_col || border(row_current, col_current) == 1) %if active residue is on image edge/border %If the current residue is on the image edge/border, place a branch %cut between the current pixel and the edge/border. branch_cuts(row_current, col_current) = 1; %make this point a branchcut to the edge connected2edge(row_current, col_current) = 1; %flag that the residue is connected to the edge (helps eliminate extra cuts to edge) residues_balanced(row_current, col_current) = 1; %mark the current residue as balanced continue; %pass to next iteration of for loop; i.e., go to next unbalanced residue end %if active residue is on image edge/border %designate the current residue as balanced and active residues_active(row_current, col_current) = 1; %mark this residue as active residues_balanced(row_current, col_current) = 1; %mark this residue as balanced %add current pixel coordinates to active list by updating the entire list row_coord_act = row_current; %add the row coordinate of the current residue to the row-coordinate active list col_coord_act = col_current; %add the col coordinate of the current residue to the col-coordinate active list num_active = length(row_coord_act); %count the active residues %num_active = sum(sum(residues_active)); %count the active residues (should = 1 at this point) charge = residues(row_current, col_current); %store the initial residue charge for the search box %set the search box and search the box perimeter for radius = 1:max_searchbox_rad %loop through search box radii % disp(['radius = ', int2str(radius),';']) ka = 1; while num_active ~= 0 %for each residue in the active list row_active = row_coord_act(ka); %row coordinate of pixel corresponding to the active residue col_active = col_coord_act(ka); %column coordinate of the pixel corresponding to the active residue % disp(['The search box is centered at the active pixel, (', int2str(row_active), ',', int2str(col_active),').']) %Define the n x n (n = 2*radius + 1) search box such that it is %centered on the active pixel (m -> rows, n-> columns, as usual). m_min = max(row_active - radius, 1); %if the active row is the first row, (row_active - radius) will return an invalid (too small) row index m_max = min(row_active + radius, num_row); %if the active row is the last row, (row_active + radius) will return an invalid (too large) row index n_min = max(col_active - radius, 1); %if the active column is the first, (col_active - radius) will return an invalid (too small) column index n_max = min(col_active + radius, num_col); %if the active column is the last, (col_active + radius) will return an invalid (too large) column index %Note: it is only necessary to search the box perimeter, as the %center has already been evaluated. Therefore, search all elements %of the top and bottom rows, and only end elements of intermediate %rows. for m = m_min:m_max %loop through search box rows, m if (m == m_min || m == m_max) %if in top or bottom rows of search box n_step = 1; %search all elements/columns else %if in intermediate rows of search box n_step = n_max - n_min; %search only edge elements/columns end %if (m == m_min || m == m_max) for n = n_min:n_step:n_max %loop through search box columns, n % disp(['Searching box pixel (m,n) = (', int2str(m), ',', int2str(n), ').']) if (m < 1 || m > num_row || n < 1 || n > num_col) %if current box pixel is outside the image range, skip this box pixel continue; else if (m == 1 || m == num_row || n == 1 || n == num_col) %if current search box pixel is on image edge/border % disp(['(m,n) = (', int2str(m), ',', int2str(n), ') is on the image edge/border.']) if (connected2edge(row_active, col_active) == 0) %if the current center pixel has not been connected to the image edge/border [~, r_nearedge, c_nearedge] = Dist2Edge(row_active, col_active, num_row, num_col, border); branch_cuts = PlaceBranchCut(branch_cuts, row_active, col_active, r_nearedge, c_nearedge); %place a branch cut between the active pixel (box center) and the nearest image edge/border % disp(['A branch cut was placed between (', int2str(row