wavelet.m

对视频序列进行三维小波变换

function varargout = wavelet(WaveletName,Level,X,Ext,Dim)
%WAVELET  Discrete wavelet transform.
%   Y = WAVELET(W,L,X) computes the L-stage discrete wavelet transform
%   (DWT) of signal X using wavelet W.  The length of X must be
%   divisible by 2^L.  For the inverse transform, WAVELET(W,-L,X)
%   inverts L stages.  Choices for W are
%     'Haar'                                      Haar
%     'D1','D2','D3','D4','D5','D6'               Daubechies'
%     'Sym1','Sym2','Sym3','Sym4','Sym5','Sym6'   Symlets
%     'Coif1','Coif2'                             Coiflets
%     'BCoif1'                                    Coiflet-like [2]
%     'Spline Nr.Nd' (or 'bior Nr.Nd') for        Splines
%       Nr = 0,  Nd = 0,1,2,3,4,5,6,7, or 8
%       Nr = 1,  Nd = 0,1,3,5, or 7
%       Nr = 2,  Nd = 0,1,2,4,6, or 8
%       Nr = 3,  Nd = 0,1,3,5, or 7
%       Nr = 4,  Nd = 0,1,2,4,6, or 8
%       Nr = 5,  Nd = 0,1,3, or 5
%     'RSpline Nr.Nd' for the same Nr.Nd pairs    Reverse splines
%     'S+P (2,2)','S+P (4,2)','S+P (6,2)',        S+P wavelets [3]
%     'S+P (4,4)','S+P (2+2,2)'
%     'TT'                                        "Two-Ten" [5]
%     'LC 5/3','LC 2/6','LC 9/7-M','LC 2/10',     Low Complexity [1]
%     'LC 5/11-C','LC 5/11-A','LC 6/14',
%     'LC 13/7-T','LC 13/7-C'
%     'Le Gall 5/3','CDF 9/7'                     JPEG2000 [7]
%     'V9/3'                                      Visual [8]
%     'Lazy'                                      Lazy wavelet
%   Case and spaces are ignored in wavelet names, for example, 'Sym4'
%   may also be written as 'sym 4'.  Some wavelets have multiple names,
%   'D1', 'Sym1', and 'Spline 1.1' are aliases of the Haar wavelet.
%
%   WAVELET(W) displays information about wavelet W and plots the
%   primal and dual scaling and wavelet functions.
%
%   For 2D transforms, prefix W with '2D'.  For example, '2D S+P (2,2)'
%   specifies a 2D (tensor) transform with the S+P (2,2) wavelet.
%   2D transforms require that X is either MxN or MxNxP where M and N
%   are divisible by 2^L.
%
%   WAVELET(W,L,X,EXT) specifies boundary handling EXT.  Choices are
%     'sym'      Symmetric extension (same as 'wsws')
%     'asym'     Antisymmetric extension, whole-point antisymmetry
%     'zpd'      Zero-padding
%     'per'      Periodic extension
%     'sp0'      Constant extrapolation
%
%   Various symmetric extensions are supported:
%     'wsws'     Whole-point symmetry (WS) on both boundaries
%     'hshs'     Half-point symmetry (HS) on both boundaries
%     'wshs'     WS left boundary, HS right boundary
%     'hsws'     HS left boundary, WS right boundary
%
%   Antisymmetric boundary handling is used by default, EXT = 'asym'.
%
%   WAVELET(...,DIM) operates along dimension DIM.
%
%   [H1,G1,H2,G2] = WAVELET(W,'filters') returns the filters
%   associated with wavelet transform W.  Each filter is represented
%   by a cell array where the first cell contains an array of
%   coefficients and the second cell contains a scalar of the leading
%   Z-power.
%
%   [X,PHI1] = WAVELET(W,'phi1') returns an approximation of the
%   scaling function associated with wavelet transform W.
%   [X,PHI1] = WAVELET(W,'phi1',N) approximates the scaling function
%   with resolution 2^-N.  Similarly,
%   [X,PSI1] = WAVELET(W,'psi1',...),
%   [X,PHI2] = WAVELET(W,'phi2',...),
%   and [X,PSI2] = WAVELET(W,'psi2',...) return approximations of the
%   wavelet function, dual scaling function, and dual wavelet function.
%
%   Wavelet transforms are implemented using the lifting scheme [4].
%   For general background on wavelets, see for example [6].
%
%
%   Examples:
%   % Display information about the S+P (4,4) wavelet
%   wavelet('S+P (4,4)');
%
%   % Plot a wavelet decomposition
%   t = linspace(0,1,256);
%   X = exp(-t) + sqrt(t - 0.3).*(t > 0.3) - 0.2*(t > 0.6);
%   wavelet('RSpline 3.1',3,X);        % Plot the decomposition of X
%
%   % Sym4 with periodic boundaries
%   Y = wavelet('Sym4',5,X,'per');    % Forward transform with 5 stages
%   R = wavelet('Sym4',-5,Y,'per');   % Invert 5 stages
%
%  ################ % 2D transform on an image
%   t = linspace(-1,1,128); [x,y] = meshgrid(t,t);
%   X = ((x+1).*(x-1) - (y+1).*(y-1)) + real(sqrt(0.4 - x.^2 - y.^2));
%   Y = wavelet('2D CDF 9/7',2,X);    % 2D wavelet transform
%   R = wavelet('2D CDF 9/7',-2,Y);   % Recover X from Y
%   imagesc(abs(Y).^0.2); colormap(gray); axis image;
%
%   % Plot the Daubechies 2 scaling function
%   [x,phi] = wavelet('D2','phi');
%   plot(x,phi);
%
%   References:
%   [1] M. Adams and F. Kossentini.  "Reversible Integer-to-Integer
%       Wavelet Transforms for Image Compression."  IEEE Trans. on
%       Image Proc., vol. 9, no. 6, Jun. 2000.
%
%   [2] M. Antonini, M. Barlaud, P. Mathieu, and I. Daubechies.  "Image
%       Coding using Wavelet Transforms."  IEEE Trans. Image Processing,
%       vol. 1, pp. 205-220, 1992.
%
%   [3] R. Calderbank, I. Daubechies, W. Sweldens, and Boon-Lock Yeo.
%       "Lossless Image Compression using Integer to Integer Wavelet
%       Transforms."  ICIP IEEE Press, vol. 1, pp. 596-599.  1997.
%
%   [4] I. Daubechies and W. Sweldens.  "Factoring Wavelet Transforms
%       into Lifting Steps."  1996.
%
%   [5] D. Le Gall and A. Tabatabai.  "Subband Coding of Digital Images
%       Using Symmetric Short Kernel Filters and Arithmetic Coding
%       Techniques."  ICASSP'88, pp.761-765, 1988.
%
%   [6] S. Mallat.  "A Wavelet Tour of Signal Processing."  Academic
%       Press, 1999.
%
%   [7] M. Unser and T. Blu.  "Mathematical Properties of the JPEG2000
%       Wavelet Filters." IEEE Trans. on Image Proc., vol. 12, no. 9,
%       Sep. 2003.
%
%   [8] Qinghai Wang and Yulong Mo.  "Choice of Wavelet Base in
%       JPEG2000."  Computer Engineering, vol. 30, no. 23, Dec. 2004.

% Pascal Getreuer 2005-2006

if nargin < 1, error('Not enough input arguments.'); end
if ~ischar(WaveletName), error('Invalid wavelet name.'); end

% Get a lifting scheme sequence for the specified wavelet
Flag1D = isempty(findstr(lower(WaveletName),'2d'));
[Seq,ScaleS,ScaleD,Family] = getwavelet(WaveletName);

if isempty(Seq)
   error(['Unknown wavelet, ''',WaveletName,'''.']);
end

if nargin < 2, Level = ''; end
if ischar(Level)
   [h1,g1] = seq2hg(Seq,ScaleS,ScaleD,0);
   [h2,g2] = seq2hg(Seq,ScaleS,ScaleD,1);

   if strcmpi(Level,'filters')
      varargout = {h1,g1,h2,g2};
   else
      if nargin < 3, X = 6; end

      switch lower(Level)
      case {'phi1','phi'}
         [x1,phi] = cascade(h1,g1,pow2(-X));
         varargout = {x1,phi};
      case {'psi1','psi'}
         [x1,phi,x2,psi] = cascade(h1,g1,pow2(-X));
         varargout = {x2,psi};
      case 'phi2'
         [x1,phi] = cascade(h2,g2,pow2(-X));
         varargout = {x1,phi};
      case 'psi2'
         [x1,phi,x2,psi] = cascade(h2,g2,pow2(-X));
         varargout = {x2,psi};
      case ''
         fprintf('\n%s wavelet ''%s'' ',Family,WaveletName);

         if all(abs([norm(h1{1}),norm(h2{1})] - 1) < 1e-11)
            fprintf('(orthogonal)\n');
         else
            fprintf('(biorthogonal)\n');
         end
         
         fprintf('Vanishing moments: %d analysis, %d reconstruction\n',...
            numvanish(g1{1}),numvanish(g2{1}));
         fprintf('Filter lengths: %d/%d-tap\n',...
            length(h1{1}),length(g1{1}));         
         fprintf('Implementation lifting steps: %d\n\n',...
            size(Seq,1)-all([Seq{1,:}] == 0));
         
         fprintf('h1(z) = %s\n',filterstr(h1,ScaleS));
         fprintf('g1(z) = %s\n',filterstr(g1,ScaleD));
         fprintf('h2(z) = %s\n',filterstr(h2,1/ScaleS));
         fprintf('g2(z) = %s\n\n',filterstr(g2,1/ScaleD));

         [x1,phi,x2,psi] = cascade(h1,g1,pow2(-X));
         subplot(2,2,1);
         plot(x1,phi,'b-');
         if diff(x1([1,end])) > 0, xlim(x1([1,end])); end
         title('\phi_1');
         subplot(2,2,3);
         plot(x2,psi,'b-');
         if diff(x2([1,end])) > 0, xlim(x2([1,end])); end
         title('\psi_1');
         [x1,phi,x2,psi] = cascade(h2,g2,pow2(-X));
         subplot(2,2,2);
         plot(x1,phi,'b-');
         if diff(x1([1,end])) > 0, xlim(x1([1,end])); end
         title('\phi_2');
         subplot(2,2,4);
         plot(x2,psi,'b-');
         if diff(x2([1,end])) > 0, xlim(x2([1,end])); end
         title('\psi_2');
         set(gcf,'NextPlot','replacechildren');
      otherwise
         error(['Invalid parameter, ''',Level,'''.']);
      end
   end

   return;
elseif nargin < 5
   % Use antisymmetric extension by default
   if nargin < 4
      if nargin < 3, error('Not enough input arguments.'); end

      Ext = 'asym';
   end

   Dim = min(find(size(X) ~= 1));
   if isempty(Dim), Dim = 1; end
end

if any(size(Level) ~= 1), error('Invalid decomposition level.'); end

NumStages = size(Seq,1);
EvenStages = ~rem(NumStages,2);

if Flag1D   % 1D Transfrom
   %%% Convert N-D array to a 2-D array with dimension Dim along the columns %%%
   XSize = size(X);    % Save original dimensions
   N = XSize(Dim);
   M = prod(XSize)/N;
   Perm = [Dim:max(length(XSize),Dim),1:Dim-1];
   X = double(reshape(permute(X,Perm),N,M));

   if M == 1 & nargout == 0 & Level > 0
      % Create a figure of the wavelet decomposition
      set(gcf,'NextPlot','replace');
      subplot(Level+2,1,1);
      plot(X);
      title('Wavelet Decomposition');
      axis tight; axis off;

      X = feval(mfilename,WaveletName,Level,X,Ext,1);

      for i = 1:Level
         N2 = N;
         N = 0.5*N;
         subplot(Level+2,1,i+1);
         a = max(abs(X(N+1:N2)))*1.1;
         plot(N+1:N2,X(N+1:N2),'b-');
         ylabel(['d',sprintf('_%c',num2str(i))]);
         axis([N+1,N2,-a,a]);
      end

      subplot(Level+2,1,Level+2);
      plot(X(1:N),'-');
      xlabel('Coefficient Index');
      ylabel('s_1');
      axis tight;
      set(gcf,'NextPlot','replacechildren');
      varargout = {X};
      return;
   end

   if rem(N,pow2(abs(Level))), error('Signal length must be divisible by 2^L.'); end
   if N < pow2(abs(Level)), error('Signal length too small for transform level.'); end

   if Level >= 0           % Forward transform
      for i = 1:Level
         Xo = X(2:2:N,:);
         Xe = X(1:2:N,:) + xfir(Seq{1,1},Seq{1,2},Xo,Ext);

         for k = 3:2:NumStages
            Xo = Xo + xfir(Seq{k-1,1},Seq{k-1,2},Xe,Ext);
            Xe = Xe + xfir(Seq{k,1},Seq{k,2},Xo,Ext);
         end

         if EvenStages
            Xo = Xo + xfir(Seq{NumStages,1},Seq{NumStages,2},Xe,Ext);
         end

         X(1:N,:) = [Xe*ScaleS; Xo*ScaleD];
         N = 0.5*N;
      end
   else                     % Inverse transform
      N = N * pow2(Level);

      for i = 1:-Level
         N2 = 2*N;
         Xe = X(1:N,:)/ScaleS;
         Xo = X(N+1:N2,:)/ScaleD;

         if EvenStages
            Xo = Xo - xfir(Seq{NumStages,1},Seq{NumStages,2},Xe,Ext);
         end

         for k = NumStages - EvenStages:-2:3
            Xe = Xe - xfir(Seq{k,1},Seq{k,2},Xo,Ext);
            Xo = Xo - xfir(Seq{k-1,1},Seq{k-1,2},Xe,Ext);
         end

         X([1:2:N2,2:2:N2],:) = [Xe - xfir(Seq{1,1},Seq{1,2},Xo,Ext); Xo];
         N = N2;
      end
   end

   X = ipermute(reshape(X,XSize(Perm)),Perm);   % Restore original array dimensions
else        % 2D Transfrom
   N = size(X);

   if length(N) > 3 | any(rem(N([1,2]),pow2(abs(Level))))
      error('Input size must be either MxN or MxNxP where M and N are divisible by 2^L.');
   end

   if Level >= 0   % 2D Forward transform
      for i = 1:Level
         Xo = X(2:2:N(1),1:N(2),:);
         Xe = X(1:2:N(1),1:N(2),:) + xfir(Seq{1,1},Seq{1,2},Xo,Ext);

         for k = 3:2:NumStages
            Xo = Xo + xfir(Seq{k-1,1},Seq{k-1,2},Xe,Ext);
            Xe = Xe + xfir(Seq{k,1},Seq{k,2},Xo,Ext);
         end