zpimalgh.m

Matlab code for encoding an unwrapping phase InSAR image based on Markov Random field

function  apha = zpimAlgh(lambda,eta,iapha,sdmu,disch,discv,qual,mask,iter)
%         apha = zpim(lambda,eta,iapha,smdu,disch,discv,qual,mask,iter);
%         zpim algorithm
%         Absolute phase unwrapping using a Bayesiean approach
%   lambda   - matrix 2*coer*power|x1x2|/|Q|
%   eta      - wrapped image (ML solution)
%   iapha    - estimated previous phase (first spep iapha = eta)
%   sdmu     - StdDev of the GaussMarkov process
%   iter     - number of iterations
%   discv    - vertical discontinuity matrix 
%              disch(i,j) = 1 means a horizontal discontinuity 
%                           between site (i,j) and site (i,j-1);
%                           disch(i,j) \in [0,1] 
%   disch    - horizontal discontinuity matrix 
%              discv(i,j) = 1 means a horizontal discontinuity 
%                           between site (i,j) and site (i-1,j)
%                           discv \in [0 1]
%   qual     - quality map  (0 - interpolate; ~=0 quality map)
%   mask     - float matrix with mask (0 - don't visit; 1-visit)
%   NOTE:    - masked pixels is treated as a non-observed pixel.
%   Author J.M. Bioucas Dias 2000
%   Topic - Interferometry


mu = 1/(sdmu.^2);
[M N]=size(lambda);
riter = 20;				% interpolation iterations
							% (masked sites)
icmiter = 4;			% icm iterations

% isolate masked&qual pixels (those set zero)
[dh dv] = linefg(mask&qual);

%----------------------------------------
%  built discontinuity matrix for z-step z()
%  disch and discv are compacted in a single float.
%  latter dicsh and discv are to be recovered (scaled by 1000) 
%  Build discontinuity matrix.
%  dh and dv are compacted in a single float by
%
%  disc = 1000*(1-dh)+1000*1001*(1-dv);
%
%  To recover dh and dv (scaled by 1000) apply:
%  dh = rem(disc,1001); 
%  dv = floor(disc/1001);
disc=dbuild(0.99.*(disch|dh),0.99.*(discv|dv));
% 0.99 assures thar ther is no pathes completely disconnected

% ------------------------
% initialize
apha = iapha;

for i=1:iter
   fprintf('Iteration - %d\n',i);
   %phaeta = apha-2*pi*floor((apha+pi)/(2*pi)); % principal phase values 
   % unwrapping step == solve for 2*pi multiples
   fprintf('Z - step no. %d\n', i);
   apha = double(zs(single(apha),single(disc)));
   fprintf('pi - step no. %d\n', i');
   apha = pistep(apha, eta, lambda, sdmu, disch|dh, discv|dv, mask, icmiter);
end

% interpol masked+disconected sites
if (sum(dh(:)+dv(:)) > 0)
   apha = reconst(apha,qual,mask,disch,discv,riter);
end