morphbwt2.m

matlab程序

function X = morphbwt2(X, Level)
%MORPHBWT2  2D morphological binary wavelet transform.
%   Y = MORPHBWT2(X,L) computes the L level decomposition of binary
%   image X using the 2D median morphological wavelet [1].  The transform
%   coefficients are also binary.  The image dimensions must be divisible
%   by 2^L.
%
%   MORPHBWT2(X,-L) is the inverse transform, reversing L levels.
%
%   Note: Using MORPHBWT2 on multi-color images works approximately but
%   does not have perfect reconstruction.  Try MEDLIFT2 for a grayscale
%   median-based morphological wavelet transform.
%
%   Example:
%   Y = morphbwt2(X,3);   % Perform 3 levels of decomposition on X
%   R = morphbwt2(Y,-3);  % Recover X from Y
%
%   Reference:
%   [1] Heijmans and Goutsias.  ``Nonlinear Multiresolution Signal
%       Decomposition Schemes--Part II: Morphological Wavelets.''  IEEE
%       Transactions on Image Processing, Vol. 9, No. 11, Nov. 2000.
%
%   See also LINBWT2, MEDLIFT2.

% Pascal Getreuer 2005

if nargin < 2, error('Not enough input arguments.'); end

N = size(X);

if rem(N,pow2(abs(Level))), error('Invalid input size.'); end

if Level > 0
   for k = 1:Level
      N = size(X)*pow2(1-k);
      A = X(1:2:N(1),1:2:N(2));
      A(:,:,2) = X(1:2:N(1),1:2:N(2));
      A(:,:,3) = X(1:2:N(1),2:2:N(2));
      A(:,:,4) = X(2:2:N(1),1:2:N(2));
      A(:,:,5) = X(2:2:N(1),2:2:N(2));
      A = sort(A,3);
      X(1:N(1),1:N(2)) = [A(:,:,3), ...
         bitxor(X(1:2:N(1),1:2:N(2)),X(2:2:N(1),1:2:N(2)));
         bitxor(X(1:2:N(1),1:2:N(2)),X(1:2:N(1),2:2:N(2))), ...
         bitxor(X(1:2:N(1),1:2:N(2)),X(2:2:N(1),2:2:N(2)))];
   end
elseif Level < 0
   for k = Level:-1
      N = size(X)*pow2(k+1);
      T = bitxor(X(1:N(1)/2,1:N(2)/2),min(min(X(1:N(1)/2,N(2)/2+1:N(2)), ...
               X(N(1)/2+1:N(1),1:N(2)/2)),X(N(1)/2+1:N(1),N(2)/2+1:N(2))));
      Y([1:2:N(1),2:2:N(1)],1:2:N(2)) = [T;bitxor(T,X(1:N(1)/2,N(2)/2+1:N(2)))];
      Y([1:2:N(1),2:2:N(1)],2:2:N(2)) = [bitxor(T,X(N(1)/2+1:N(1),1:N(2)/2));
         bitxor(T,X(N(1)/2+1:N(1),N(2)/2+1:N(2)))];
      X(1:N(1),1:N(2)) = Y;
   end
end