To use the double-density discrete wavelet transform for 2-D signal processing, we must implement a two-dimensional analysis and synthesis filter bank structure. This can simply be done by alternatively applying the transform first to the rows, then to the columns of an image. This gives rise to nine 2-D subbands, one of which is the 2-D lowpass scaling filter, and the other eight of which make up the eight 2-D wavelet filters, as shown in Figure 7 below.

Figure 7. An Oversampled Filter Bank for 2-D Images.

We implement the 2-D analysis filter bank with the afb3_2D.m MATLAB function. This function returns the subband images as two variables: (1) lo, which is the lowpass subband image, and (2) hi, which is a cell array containing the eight other highpass subband images.

function [lo, hi] = afb3_2D(x, af1, af2)

% 2-D Analysis Filter Bank
%
% USAGE:
% [lo, hi] = afb3_2D(x, af1, af2);
% INPUT:
%    x - N x M matrix
%    af1 - analysis filter for the columns
%    af2 - analysis filter for the rows
%    af(:, 1) - lowpass filter
%    af(:, 2) - first highpass filter
%    af(:, 3) - second highpass filter
% OUTPUT:
%     lo - lowpass subband
%     hi - cell array containing the eight following highpass subbands:
%         1) hi{1} = 'lo hi1' subband
%         2) hi{2} = 'lo hi2' subband
%         3) hi{3} = 'hi1 lo' subband
%         4) hi{4} = 'hi1 hi1' subband
%         5) hi{5} = 'hi1 hi2' subband
%         6) hi{6} = 'hi2 lo' subband
%         7) hi{7} = 'hi2 hi1' subband
%         8) hi{8} = 'hi2 hi2' subband

if nargin < 3
   af2 = af1;
end

% filter along columns
[L, H1, H2] = afb3_2D_A(x, af1, 1);

% filter along rows
[lo,    hi{1}, hi{2}] = afb3_2D_A(L, af2, 2);
[hi{3}, hi{4}, hi{5}] = afb3_2D_A(H1, af2, 2);
[hi{6}, hi{7}, hi{8}] = afb3_2D_A(H2, af2, 2);

Similarly, we implement the 2-D synthesis filter bank with the sfb3_2D.m MATLAB function. This function reconstructs the output image from the nine subband images.

As with the 1-D transform, the 2-D double-density DWT of an input image x is implemented by recursively applying the 2-D analysis filter bank on the lowpass subband image. The MATLAB function double_f2D.m computes the double-density DWT, w, of an input image x, for any J number of scales. MATLAB's cell array data structure is used to store each of the nine subband signals such that for j = 1 - J, w{j}{d} is the high frequency subband signal produced at stage j, where d = 1 - 8. The low frequency subband signal produced at stage j is then stored in w{J+1}.

function w = double_f2D(x, J, af1, af2)

% Forward Double-Density Discrete 2-D Wavelet Transform
%
% USAGE:
%   w = double_f2D(x, J, af)
% INPUT:
%   x - N by M matrix
%      1) M, N are both even
%      2) min(M,N) >= 2^(J-1)*length(af)
%   J - number of stages
%  af - analysis filters
% OUPUT:
%   w - cell array of wavelet coefficients

if nargin < 4
    af2 = af1;
end

for k = 1:J
    [x w{k}] = afb3_2D(x, af1, af2);
end

w{J+1} = x;

The inverse double-density DWT is computed with the double_i2D.m MATLAB function. This function computes the inverse double-density DWT at the Jth scale by repeatedly calling the sfb3_2D.m function.

To verify the perfect reconstruction condition, we write another small piece of MATLAB code (as seen below), similar to the one in the previous section. First, we generate a random 2-D input image x of size 128x64. The analysis and synthesis filter bank is then computed first from the filters1.m function, and later from the filters2.m filter bank. The 3-scale forward double-density 2-D DWT is computed with the double_f2D.m function to generate a function w. We then apply the inverse double-density 2-D DWT to the wavelet coefficients in w using the double_i2D.m function to show that it reconstructs the original image x. To ensure that the perfect reconstruction condition has been met, we then find the difference between the input image x and the output image y. As you can see, the differences found here are again significantly close to zero, and will be considered as zero.

>> [af, sf] = filters1;		% analysis and synthesis filters
>> x = rand(128,64);		% create generic 2-D signal
>> w = double_f2D(x, 3, af);	% analysis filter bank (3 stages)
>> y = double_i2D(w, 3, sf);	% synthesis filter bank (3 stages)
>> err = x - y;			% compute error signal
>> max(max(abs(err)))		% verify that error is zero

ans =

  2.9310e-014

>> [af, sf] = filters2;		% analysis and synthesis filters
>> x = rand(128,64);		% create generic 2-D signal
>> w = double_f2D(x, 3, af);	% analysis filter bank (3 stages)
>> y = double_i2D(w, 3, sf);	% synthesis filter bank (3 stages)
>> err = x - y;			% compute error signal
>> max(max(abs(err)))		% verify that error is zero

ans =

  7.7827e-013

As previously mentioned, there are eight wavelets associated with the 2-D double-density wavelet transform. We can plot the wavelets with the MATLAB code fragment below (doubledwt2D_plots.m).

>> [af, sf] = filters1;		% analysis and synthesis filters
>> J = 5;			% five stages
>> L = 8*2^(J+1);		
>> N = L/2^J;
>> x = zeros(L,8*L);		% create zero signal
>> w = double_f2D(x,J,af);	% analysis filter bank (5 stages)
>> w{J}{1}(N/2,N/2+0*N) = 1;	% set single coefficient to 1 for each of 
>> w{J}{2}(N/2,N/2+1*N) = 1;	% the eight 2-D wavelets
>> w{J}{3}(N/2,N/2+2*N) = 1;
>> w{J}{4}(N/2,N/2+3*N) = 1;
>> w{J}{5}(N/2,N/2+4*N) = 1;
>> w{J}{6}(N/2,N/2+5*N) = 1;
>> w{J}{7}(N/2,N/2+6*N) = 1;
>> w{J}{8}(N/2,N/2+7*N) = 1;
>> y = double_i2D(w,J,sf);	% synthesis filter bank (5 stages)
>> figure(1)		
>> clf
>> imagesc(y);			% construct an image of all eight wavelets
>> axis equal	
>> axis off
>> colormap(gray(128))		% convert to a grayscale image

As in the 1-D case, we set all of the wavelet coefficients to zero, with the exception of one wavelet coefficient in each of the eight subbands. We then take the inverse wavelet transform and plot the images. The following figure illustrates the eight wavelets as gray scale images.

Figure 8. 2-D Double-Density Wavelets.

Note that the first two wavelets are oriented in the vertical direction, while the third and sixth are oriented in the horizontal direction. The remaining wavelets do not have any specific orientations, but rather combine the two diagonal orientations, which gives rise to the checkerboard affect. Because half of the wavelets have no dominating orientation, the 2-D double-density DWT is poor at isolating the two diagonal orientations.

Summary of Programs:

1. filters1.m - First set of asymmetric perfect reconstruction filters
2. filters2.m - Second set of symmetric perfect reconstruction filters
3. afb3_2D.m - 2-D analysis filter bank
4. sfb3_2D.m - 2-D synthesis filter bank
5. double_f1D.m - Forward 2-D double-density wavelet transform
6. double_i1D.m - Inverse 2-D double-density wavelet transform
7. cshift2D.m - 2-D circular shift
8. doubledwt2D_plots.m - Plot 2-D double-density wavelets