pistep.m

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

function apha = pistep(iapha,eta,lambda,sdmu,disch,discv,mask,icmiter)
%   apha = pistep(eta,lambda,sdmu,disch,discv,mask,icmiter);
%   Absolute phase unwrapping using a Bayesiean approach
%   iapha    - estimated previous phase (first spep iapha = eta)
%   eta      - observed interferogram;
%   lambda   - matrix 2*coer*power|x1x2|/|Q|
%   sdmu     - StdDev of the GaussMarkov process
%   icmiter  - number of ICM iterations
%   discv    - vertical discontinuity matrix 
%              discv(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 
%              disch(i,j) = 1 means a horizontal discontinuity 
%                           between site (i,j) and site (i-1,j)
%                           discv \in [0 1]
%   mask     - float matrix with mask (0 - don't visit; 1-visit)
%         
%   NOTE:    - Discontinuity field must disconnect pixels no visited
%              form those visited
%   Author J.M. Bioucas Dias 2000
%   Topic - Interferometry

[M N] = size(eta);
mu = 1/(sdmu.^2);


apha = iapha;

phaeta = apha-2*pi*floor((apha+pi)/(2*pi)); % principal phase values 
phak   = apha - phaeta;     % 2*pi multiples


for itericm = 1:icmiter   % icm iter
   for nextsite=1:M*N  
      if mask(nextsite) ~= 0
         % find nextsite neighbors 
         c= floor((nextsite-1)/M)+1;  % column
         l= nextsite - (c-1)*M;       % line         
         % set 1  moreneighs(i) to 1 if there is no discontinuity
         % to the correspondent neighbor
         neighs = [0 0 0 0]; %[(l,c+1) (l-1,c)(l,c-1)(l+1,c)]
         discweight = [0 0 0 0]; 
         if ((c+1) <= N) 
            neighs(1) = nextsite+M;
            discweight(1) = 1 - discv(l,c+1);
         end
         if ((l-1) >= 1) 
            neighs(2) = nextsite-1;
            discweight(2) = 1 - disch(l,c);
         end
         if ((c-1) >= 1) 
            neighs(3) = nextsite-M;
            discweight(3) = 1 - discv(l,c);
         end
         if ((l+1) <= M) 
            neighs(4) = nextsite+1;         
            discweight(4) = 1 - disch(l+1,c);
         end
         discsum = sum(discweight);
         meanpha = 0;
         for j=1:4
            if neighs(j) ~= 0
               meanpha  = meanpha  + apha(neighs(j))*discweight(j);
            end
         end
         if discsum ~= 0.0
            meanpha  = meanpha /discsum;
            beta = 2*lambda(nextsite)/(mu*discsum);
            %icm steep
            apha(nextsite) = icmstep(eta(nextsite),meanpha-phak(nextsite),beta)+phak(nextsite);
         end
      end % end if mask
  end  % end nextsite
end % end icmiter