wavecdf97.m

function c = wavecdf97(x, nlevel)
%WAVECDF97: Multi-level discrete 2-D wavelet transform 
%with the Cohen-Daubechies-Feauveau (CDF) 9/7 wavelet. 
%
% c = wavecdf97(x, nlevel) does the follows according to the value of 
%   nlevel:
%   nlevel > 0:   decomposes 2-dimension matrix x up to nlevel level;
%   nlevel < 0:   does the inverse transform to nlevel level;
%   nlevel = 0:   sets c equal to x;
%   omitted:      does the same as nlevel=5.  
%
% The boundary handling method is symmetric extension. 
%
% x may be of any size; it need not have size divisible by 2^L.
%   For example, if x has length 9, one stage of decomposition
%   produces a lowpass subband of length 5 and a highpass subband
%   of length 4.  Transforms of any length have perfect
%   reconstruction (exact inversion).
%   NOTE: the 5 lines above are quoted directly form [3].
%   
% If nlevel is so large that the approximation coefficients become 
%   a 1-D array, any further decomposition will be performed as for 1-D 
%   decomposition until the approximation coefficients be a scale number.  
%
% Lifting algorithm is not used here; we use subband filters directly.
%   Lifting algorithm and spline 5/3 wavelets and other jpeg2000 related 
%   codes will be available soon. 
%
% Example:
%   Y = wavecdf97(X, 5);    % Decompose image X up to 5 level
%   R = wavecdf97(Y, -5);   % Reconstruct from Y
%
% You can test wavecdf97.m with the following lines:     
%   % get a 2-D uint8 image 
%   x=imread('E:\study\jpeg2000\images\lena.tif');
%   % decompose
%   y=wavecdf97(x,2);
%   % show decomposed result 
%   figure;imshow(mat2gray(y));
%   % reconstruct without change of anything
%   ix=wavecdf97(y,-2);
%   % show and compare the original and reconstructed images
%   figure;subplot(1,2,1);imshow(x);subplot(1,2,2);imshow(uint8(ix));
%   % look at the MSE difference 
%   sum(sum((double(x)-ix).^2))/numel(x)
%
% Reference:
%   [1] D.S.Taubman et al., JPEC2000 Image Compression: F. S. & P.,
%       Chinese Edition, formula 10.6-10.9 in section 10.3.1 
%       and formula 10.13 in section 10.4.1.
%   [2] R.C.Gonzalez et al., Digital Image Processing Using MATLAB, 
%       Chinese Edition, function wavefast in section 7.2.2.
%   [3] Pascal Getreuer, waveletcdf97.m from Matlab file Exchange website
%   [4] Matlab files: biorwavf.m, wavdec2.m, wawrec2.m, etc.
%   
% Contact information: 
%   Email/MSN messenger:  wangthth@hotmail.com
%
% Tianhui Wang at Beijing, China,   July, 2006
%                  Last Revision:   Aug 5, 2006
%---------------------- input arguments checking  ----------------------%
error(nargchk(1,2,nargin));
if nargin == 1
    nlevel = 5; % default level
end
% check x
if ~isreal(x) || ~isnumeric(x) || (ndims(x) > 2)
    error('WAVELIFT:InArgErr', ['The first argument must' ...
        ' be a real, numeric 2-D or 1-D matrix.']);
end
if isinteger(x)
    x = double(x);
end
% check nlevel
if ~isreal(nlevel) || ~isnumeric(nlevel) || round(nlevel)~=nlevel
    error('WAVELIFT:InArgErr', ['The 2nd argument shall be ' ...
        'a real and numeric integer.']);
end
%---------------- forming low-pass and high-pass filters ---------------%
% CDF 9/7 filters: decomposition low-pass lp and high-pass hp
%                  reconstruction low-pass lpr and high-pass hpr
% The filter coefficients have several forms.
% What D.S.Taubman et al. suggest in [1] are used here:
lp = [.026748757411 -.016864118443 -.078223266529 .266864118443];
lp = [lp .602949018236 fliplr(lp)];
hp = [.045635881557 -.028771763114 -.295635881557];
hp = [hp .557543526229 fliplr(hp)];
lpr = hp .* [-1 1 -1 1 -1 1 -1] * 2;
hpr = lp .* [1 -1 1 -1 1 -1 1 -1 1] * 2;
% Matlab 'bior4.4' use the varied version (see Matlab's biorwavf.m):
%  lp=lp*sqrt(2);hp=hp*(-sqrt(2));lpr=lpr*(1/sqrt(2));hpr=hpr*(-1/sqrt(2));
% P.Getreuer's waveletcdf97.m [3] alters the Taubman's version as follows:
%  lp=lp*sqrt(2);hp=hp*sqrt(2);lpr=lpr*(1/sqrt(2));hpr=hpr*(1/sqrt(2));
% while R.C.Gonzalez et al in [2] alter the Taubman's version as follows:
%  lp=lp;hp=hp*(-2);lpr=lpr;hpr=hpr*(-1/2);
%----------------  remain unchanged when nlevel = 0  -------------------%
if nlevel == 0
    c = x;
%--------------------  decomposition,  if nlevel < 0  ------------------%
elseif nlevel > 0
    c = zeros(size(x));
    x = double(x);
    for i = 1:nlevel
        % [ll, hl; lh, hh]: 1-level FWT for x 
        temp = symconv2(x, hp, 'col');    % high filtering
        temp = temp(2:2:end, :);          % down sampling
        hh = symconv2(temp, hp, 'row');   % high filtering 
        hh = hh(:, 2:2:end);              % down sampling
        lh = symconv2(temp, lp, 'row');   % low filtering
        lh = lh(:, 1:2:end);              % down sampling
        
        temp = symconv2(x, lp, 'col');    % low filtering
        temp = temp(1:2:end, :);          % down sampling
        hl = symconv2(temp, hp, 'row');   % high filtering
        hl = hl(:, 2:2:end);              % down sampling
        ll = symconv2(temp, lp, 'row');   % low filtering
        ll = ll(:, 1:2:end);              % down sampling
        % update coefficient matrix
        c(1:size(x,1), 1:size(x,2)) = [ll, hl; lh, hh];
        % replace x with ll for next level FWT
        x = ll;
        % give a warning if nlevel is too large
        if size(x,1)<=1 && size(x,2)<=1 && i~=nlevel
            warning('WAVECDF97:DegradeInput', ['Only decompose to ' ...
                num2str(i) '-level instead of ' num2str(nlevel) ...
                ', \nas the approximation coefficients at ' num2str(i) ...
                '-level has row or/and column of length 1.']);
            break
        end
    end
%--------------------  reconstruction,  if nlevel < 0  -----------------%
else
    sx = size(x);
    % find reconstruction level
    nl = -nlevel;
    while sx(1)/2^nl<=1/2 && sx(2)/2^nl<=1/2,  nl = nl-1;  end
    if nl ~= -nlevel 
        warning('WAVECDF97:DegradeInput', ['Only reconstruct to ' ...
            num2str(nl) '-level instead of ' num2str(-nlevel) ...
            ', \nas the approximation coefficients at ' num2str(nl) ...
            '-level has row or/and column of length 1.']);
    end
    % nl-level reconstruction
    for i = 1:nl
        % find the target ll hl lh hh blocks
        sLL = ceil(sx/2^(nl-i+1));
        sConstructed = ceil(sx/2^(nl-i));
        sHH = sConstructed - sLL;
        lrow = sConstructed(1); lcol = sConstructed(2);
        ll = x(1:sLL(1), 1:sLL(2));
        hl = x(1:sLL(1), sLL(2)+1:sLL(2)+sHH(2));
        lh = x(sLL(1)+1:sLL(1)+sHH(1), 1:sLL(2));    
        hh = x(sLL(1)+1:sLL(1)+sHH(1), sLL(2)+1:sLL(2)+sHH(2));
        % upsample rows and low filter
        temp = zeros(sLL(1), lcol); temp(:, 1:2:end) = ll;
        ll = symconv2(temp, lpr, 'row'); 
        % upsample rows and high filter    
        temp = zeros(sLL(1), lcol); temp(:, 2:2:end) = hl;
        hl = symconv2(temp, hpr, 'row');
        % upsample columns and low filter      
        temp = zeros(lrow, lcol); temp(1:2:end, :) = ll + hl;
        l = symconv2(temp, lpr, 'col');  
        % upsample rows and high filter       
        temp = zeros(sHH(1), lcol); temp(:, 1:2:end) = lh;
        lh = symconv2(temp, lpr, 'row');
        % upsample rows and high filter       
        temp = zeros(sHH(1), lcol); temp(:, 2:2:end) = hh;
        hh = symconv2(temp, hpr, 'row');
        % upsample rows and high filter       
        temp = zeros(lrow, lcol); temp(2:2:end, :) = lh + hh;
        h = symconv2(temp, hpr, 'col');
        % update x with the new ll, ie. l+h
        x(1:lrow, 1:lcol) = l + h;
    end    
    % output
    c = x;
end
%------------------------- internal function  --------------------------%
%       2-dimension convolution with edges symmetrically extended       %
%-----------------------------------------------------------------------%
function y = symconv2(x, h, direction)
% symmetrically extended convolution(see section 6.5.2 in [1]):
%    x[n], E<=n<=F-1, is extended to x~[n] = x[n], 0<=n<=N-1;
%                                  x~[E-i] = x~[E+i], for all integer i
%                                x~[F-1-i] = x~[F-1+i], for all integer i
%    For odd-length h[n], to convolve x[n] and h[n], we just need extend x 
%    by (length(h)-1)/2  for both left and right edges. 
% The symmetric extension handled here is not the same as in Matlab 
%  wavelets toolbox nor in [2]. The last two use the following method:
%    x[n], E<=n<=F-1, is extended to x~[n] = x[n], 0<=n<=N-1;
%                                  x~[E-i] = x~[E+i-1], for all integer i
%                                x~[F-1-i] = x~[F+i], for all integer i 
l = length(h); s = size(x);
lext = (l-1)/2; % length of h is odd 
h = h(:)'; % make sure h is row vector 
y = x;
if strcmp(direction, 'row') % convolving along rows
    if ~isempty(x) && s(2) > 1 % unit length array skip convolution stage
        for i = 1: lext
            x = [x(:, 2*i), x, x(:, s(2)-1)]; % symmetric extension
        end
        x = conv2(x, h);
        y = x(:, l:s(2)+l-1); 
    end
elseif strcmp(direction, 'col') % convolving along columns
    if ~isempty(x) && s(1) > 1 % unit length array skip convolution stage
        for i = 1: lext 
            x = [x(2*i, :); x; x(s(1)-1, :)]; % symmetric extension
        end
        x = conv2(x, h');
        y = x(l:s(1)+l-1, :);
    end
end    
% EOF