zpim.m

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

function  apha = zpim(eta,varargin)
% ZPIM algorithm for  Matlab 5.x
% Version 3.0, June 20 2003
% Copyright (c) Jose Bioucas Dias.
%
% Absolute phase estimation in interderometic SAR (InSAR)
% pim algoritm w as proposed in IEEEIP, vol 11, no. 4, pp. 408-422, 2002.
% See also http://www.lx.it.pt/~bioucas/publications
%
%==========================================================================
%--Description of the parameters
%
%   lambda   - matrix 2*coer*power|x1x2|/|Q|
%   eta      - wrapped image (ML solution)
%   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 are treated as a non-observed pixel.
%
%============================================================================
%--Example1
%
% apha = zpim(eta); eta contains the interferogram (wrapped image). 
% 
% Defaults:
%	lambda           = 2*ones(M,N);
%  iapha				  = eta				% phase obtained in previous iterations
%	sdmu				  = 0.8;				% "a priori" phase standard deviation
%	disch            = zeros(M,N);	% horizontal discontinuity
%	discv				  = zeros(M,N);	% vertical discontinuity
%	qual				  = ones(M,N);    % quality map
%	mask				  = ones(M,N);		% defines where to estimate phase 
%
%============================================================================
%--Example2
%
% apha = zpim(eta, 'Lambda',lambda); 
% eta - contains the interferogram (wrapped image). 
% lambda   - matrix 2*coer*power|x1x2|/|Q|
% 
% Remaining pars: as in Example 1.
%           
%============================================================================
%--Example3
%
% apha = zpim(eta, 'Lambda',lambda,'Sdmu',sdmu,'Disch',disch,'Discv',discv,...
%             'Qual',qual,'Iter',iter); 
%           
%
%

[M,N] = size(eta);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Read the optional parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Default varargin parameters
lambda           = 2*ones(M,N);
sdmu				  = 0.8;				% "a priori" phase standard deviation
iapha            = eta;				% estimate obtained in previous iterations
disch            = zeros(M,N);	% horizontal discontinuity
discv				  = zeros(M,N);	% vertical discontinuity
qual				  = ones(M,N);    % quality map
mask				  = ones(M,N);		% defines where to estimate phase 
iter				  = 4;				% zpim iterations


if (rem(length(varargin),2)==1)
  error('Optional parameters should always go by pairs');
else
  for i=1:2:(length(varargin)-1)
    % change the value of parameter
    switch varargin{i}
     case 'Lambda'
        lambda = varargin{i+1};	
     case 'Iapha'
      iapha = varargin{i+1};
     case 'Sdmu'
      sdmu = varargin{i+1};
     case 'Disch'
      disch = varargin{i+1};
     case 'Discv'
      discv = varargin{i+1};
     case 'Qual'
      qual = varargin{i+1};
     case 'Mask'
      mask = varargin{i+1};
     case 'Iter'
      iter = varargin{i+1};      
     otherwise
      % Hmmm, something wrong with the parameter string
      error(['Unrecognized parameter: ''' varargin{i} '''']);
    end;	%end swich
  end; %end for
end %end if

% call zpimAlgh
apha = zpimAlgh(lambda,eta,iapha,sdmu,disch,discv,qual,mask,iter);