3-D Discrete Wavelet Transform



1. 3D Filter Banks


To use the wavelet transform for volume and video processing we must implement a 3D version of the analysis and synthesis filter banks. In the 3D case, the 1D analysis filter bank is applied in turn to each of the three dimensions. If the data is of size N1 by N2 by N3, then after applying the 1D analysis filter bank to the first dimension we have two subband data sets, each of size N1/2 by N2 by N3. After applying the 1D analysis filter bank to the second dimension we have four subband data sets, each of size N1/2 by N2/2 by N3. Applying the 1D analysis filter bank to the third dimension gives eight subband data sets, each of size N1/2 by N2/2 by N3/2. This is illustrated in the diagram below.



Fig. 3.1 The resulution of a 3-D signal is reduced in each dimension


The block diagram of the 3D analysis filter bank is shown in the following diagram.

Fig. 3.2 Block diagram of a single level decomposition for the 3D DWT [4]


The 3D analysis filter bank is implemented with the Matlab function afb3D.m. This function calls a sub-function, afb3D_A.m, which applies the 1D analysis filter bank along one dimension only. The function afb3D.m returns two variables: lo is the lowpass subband data set; hi is a cell array containing the seven other subband data sets. Although not shown here, the function afb3D_A is contained in the same file as afb3D.

Table 3.1 Matlab function afb3D.m

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

% 3D Analysis Filter Bank
%
% [lo, hi] = afb3D(x, af1, af2, af3);
% INPUT:
%    x - N1*N2*N3 matrix, where min(Ni) > 2*length(filter)
%           (Ni are even)
%    afi - analysis filter for dimension i
%    afi(:, 1) - lowpass filter
%    afi(:, 2) - highpass filter
% OUTPUT:
%     lo, hi - lowpass, highpass subbands

if nargin < 3
   af2 = af1;
   af3 = af1;
end

% filter along dimension 1
[L, H] = afb3D_A(x, af1, 1);

% filter along dimension 2
[LL LH] = afb3D_A(L, af2, 2);
[HL HH] = afb3D_A(H, af2, 2);

% filter along dimension 3
[LLL LLH] = afb3D_A(LL, af3, 3);
[LHL LHH] = afb3D_A(LH, af3, 3);
[HLL HLH] = afb3D_A(HL, af3, 3);
[HHL HHH] = afb3D_A(HH, af3, 3);

lo    = LLL;
hi{1} = LLH;
hi{2} = LHL;
hi{3} = LHH;
hi{4} = HLL;
hi{5} = HLH;
hi{6} = HHL;
hi{7} = HHH;

The 3D synthesis filter bank, implemented with the Matlab function sfb3D.m and sfb3D_A.m, combines the eight subband data sets to obtain the original image of size N1 by N2 by N3. The following code fragment verifies the perfect reconstruction property of the 3D analysis/synthesis filter bank.


VERIFY PERFECT RECONSTRUNCTION

>> x = rand(32,64,16);
>> [af, sf] = farras;
>> [lo, hi] = afb3D(x, af, af, af);
>> y = sfb3D(lo, hi, sf, sf, sf);
>> err = x - y;
>> max(max(max(abs(err))))

ans =

   5.7510e-014
   

2. 3D Discrete Wavelet Transform


As in the 1D case, the separable 3D discrete wavelet transform of a signal x is implemented by iterating the 3D analysis filter bank on the lowpass subband data. In this case, at each scale there are seven highpass subbands instead of one. The following function, dwt3D.m, computes the J-scale 3D DWT of an data set x by repeatedly calling afb3D.m.

Table 3.2 Matlab function dwt3D.m

function w = dwt3D(x, J, af)

% discrete 3-D wavelet transform
%
% w = dwt3D(x, stages, af)
% INPUT:
%    x - N1*N2*N3 matrix (Ni even)
%    J - number of stages
%    af  - analysis filters
% OUPUT:
%    w - cell array of wavelet coefficients
%
% % Example:
% [af, sf] = farras;
% x = rand(128,64,64);
% w = dwt3D(x,3,af);
% y = idwt3D(w,3,sf);
% err = x-y; 
% max(max(max(abs(err))))

for k = 1:J
    [x w{k}] = afb3D(x, af, af, af);
end
w{J+1} = x;

Again, w is a Matlab cell array; for j = 1..J, d = 1..7, w{j}{d} is one of the seven subband data sets produced at stage j. w{J+1} is the lowpass subband data set produced at the last stage. The function idwt3D.m computes the inverse 3D DWT.


The following code fragment verifies the perfect reconstruction of the 3D DWT. We create a random input signal x of size 128 by 64 by 64, apply the 3-scale 3D DWT and its inverse, and show it reconstructs x from the wavelet coefficients in w.


VERIFY PERFECT RECONSTRUNCTION

 >> [af, sf] = farras;
 >> x = rand(128,64,64);
 >> w = dwt3D(x,3,af);
 >> y = idwt3D(w,3,sf);
 >> err = x-y;
 >> max(max(max(abs(err))))

 ans =

   1.2845e-013

There are seven wavelets associated with the separable 3D wavelet transform. It is a bit difficult to display the wavelets because they are functions of three variables. However, the isosurface provides a way to represent a function of three variables, similar to a contour plot of function of two variables. Shown below is an isosurface of one of the seven wavelets. To interpret the isosurface note that, like a contour plot, the points on the surfaces are points where the function is equal-valued.



The blue and red surfaces represent the points where the wavelet is equal-valued (blue: w(x, y, t) = 0.3; red: w(x, y, t) = -0.3).


Note that this wavelet has the same checkerboard artifact present in the separable 2-D case - it does not have a dominant orientation. For this reason, the separable 3-D DWT does a somewhat poor job at isolating in its subbands distinct orientations of 3D data sets. The isosurface figure was produced using the Matlab program dwt3D_plots.m.

Table 3.3 Matlab function dwt3D_plots.m

% dwt3D_plots
% DISPLAY 3D WAVELETS OF DWT3D.M

[af, sf] = farras;
J = 4;
L = 3*2^(J+1);
N = L/2^J;
x = zeros(L,L,L);
w = dwt3D(x,J,af);
w{J}{7}(N/2,N/2,N/2) = 1;
y = idwt3D(w,J,sf);
figure(1)
clf
v = 1:L;
S = 0.017;
p1 = patch(isosurface(v,v,v,y,S));
isonormals(v,v,v,y,p1);
set(p1,'FaceColor','red','EdgeColor','none'); 
hold on
p2 = patch(isosurface(v,v,v,y,-S));
isonormals(v,v,v,y,p2);
set(p2,'FaceColor','blue','EdgeColor','none'); 
hold off
daspect([1 1 1]);
view(-30,30); 
camlight;
lighting phong
grid
axis([32 48 32 48 32 48])
set(gca,'fontsize',7)
title('3-D WAVELET ISOSURFACE (SEPARABLE TRANSFORM)')
set(gcf,'paperposition',[0.5 0.5 0 0]+[0 0 3 3])
print -djpeg95 dwt3D_plots

As in the 1D and 2D case, we first set allcoefficients to zero for the exception of one wavelet coefficient, and then we then take the inverse wavelet transform. The isosurface command in Matlab is used to produce this figure. Different sets of filters will produce different wavelets, however, the checkerboard artifact will still be present for the separable 3D DWT.


Summary of programs:


  1. afb3D.m - Analysis filter bank (single stage)
  2. sfb3D.m - Synthesis filter bank (single stage)
  3. dwt3D.m - Forward wavelet transform (multi-stage)
  4. idwt3D.m - Inverse wavelet transform (multi-stage)
  5. cshift3D.m - Circular shift
  6. dwt3D_plots.m - Plot 3D wavelet isosurface


Separable DWT
  1-D DWT
  2-D DWT
  3-D DWT
Dual-tree DWT
  1-D DWT
  2-D DWT
  3-D DWT
Denoising
  Thresholding
  BiShrink