pyramids.m

matlab的steel金字塔小波分解源代码

%% except for the coarsest subband which is lowpass.
for lev=1:size(locations,2)
  subplot(2,2,lev);
  showIm(projection(locations(lev),:));
  axis([0 sz -0.3 0.8]);
end
  
%% Now consider the frequency response of these functions, plotted over the
%% range [-pi,pi]:
for lev=1:size(locations,2)
  subplot(2,2,lev);
  proj = projection(locations(lev),:);
  plot(pi*[-32:31]/32,fftshift(abs(fft(proj',64))));
  axis([-pi pi -0.1 3]);
end

%% The first projection function is highpass, and the second is bandpass.  Both
%% of these look something like the Laplacian (2nd derivative) of a Gaussian.
%% The last is lowpass, as are the basis functions.  Thus, the basic operation
%% used to create each level of the pyramid involves a simple highpass/lowpass
%% split.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% QMF/WAVELET PYRAMIDS.

%% Two things about Laplacian pyramids are a bit unsatisfactory.
%% First, there are more pixels (coefficients) in the representation
%% than in the original image. Specifically, the 1-dimensional
%% transform is overcomplete by a factor of 4/3, and the 2-dimensional
%% transform is overcomplete by a factor of 2.  Secondly, the
%% "bandpass" images (fineN) do not segregate information according to
%% orientation.

%% There are other varieties of pyramid.  One type that arose in the
%% speech coding community is based on a particular pairs of filters
%% known as a "Quadrature Mirror Filters" or QMFs.  These are closely
%% related to Wavelets (essentially, they are approximate wavelet
%% filters).

%% Recall that the Laplacian pyramid is formed by simple hi/low
%% splitting at each level.  The lowpass band is subsampled by a
%% factor of 2, but the highpass band is NOT subsampled.  In the QMF
%% pyramid, we apply two filters (hi- and lo- pass) and subsample BOTH
%% by a factor of 2, thus eliminating the excess coefficients of the
%% Laplacian pyramid.

%% The two filters must have a specific relationship to each
%% other. In particular, let n be an index for the filter samples.
%% The highpass filter may be constructed from the lowpass filter by
%% (1) modulating (multiplying) by (-1)^n (equivalent to shifting by
%% pi in the Fourier domain), (2) flipping (i.e., reversing the order
%% of the taps), (3) spatially shifting by one sample.  Try to
%% convince yourself that the resulting filters will always be
%% orthogonal to each other (i.e., their inner products will be zero)
%% when shifted by any multiple of two.

%% The function modulateFlip performs the first two of these operations.  The
%% third (spatial shifting) step is built into the convolution code.
flo = namedFilter('qmf9')';
fhi = modulateFlip(flo)';
subplot(2,1,1); lplot(flo); axis([0 10 -0.5 1.0]); title('lowpass');
subplot(2,1,2); lplot(fhi); axis([0 10 -0.5 1.0]); title('highpass');

%% In the Fourier domain, these filters are (approximately)
%% "power-complementary": the sum of their squared power spectra is
%% (approximately) a constant.  But note that neither is a perfect
%% bandlimiter (i.e., a sinc function), and thus subsampling by a
%% factor of 2 will cause aliasing in each of the subbands.  See below
%% for a discussion of the effect of this aliasing.

%% Plot the two frequency responses:
freq = pi*[-32:31]/32;
subplot(2,1,1);
plot(freq,fftshift(abs(fft(flo,64))),'--',freq,fftshift(abs(fft(fhi,64))),'-');
axis([-pi pi 0 1.5]); title('FFT magnitudes');
subplot(2,1,2);
plot(freq,fftshift(abs(fft(flo,64)).^2)+fftshift(abs(fft(fhi,64)).^2));
axis([-pi pi 0 2.2]); title('Sum of squared magnitudes');

%% We can split an input signal into two bands as follows:
sig = mkFract([1,64],1.6);
subplot(2,1,1); showIm(sig,'auto1','auto','sig');
lo1 = corrDn(sig,flo,'reflect1',[1 2],[1 1]);
hi1 = corrDn(sig,fhi,'reflect1',[1 2],[1 2]);
subplot(2,1,2); 
showIm(lo1,'auto1','auto','low and high bands'); hold on; plot(hi1,'--r'); hold off; 

%% Notice that the two subbands are half the size of the original
%% image, due to the subsampling by a factor of 2.  One subtle point:
%% the highpass and lowpass bands are subsampled on different
%% lattices: the lowpass band retains the odd-numbered samples and the
%% highpass band retains the even-numbered samples.  This was the
%% 1-sample shift relating the high and lowpass kernels (mentioned
%% above).  We've used the 'reflect1' to handle boundaries, which
%% works properly for symmetric odd-length QMFs.

%% We can reconstruct the original image by interpolating these two subbands
%% USING THE SAME FILTERS:
reconlo = upConv(lo1,flo,'reflect1',[1 2]);
reconhi = upConv(hi1,fhi,'reflect1',[1 2],[1 2]);
subplot(2,1,2); showIm(reconlo+reconhi,'auto1','auto','reconstructed');
imStats(sig,reconlo+reconhi);

%% We have described an INVERTIBLE linear transform that maps an input
%% image to the two images lo1 and hi1.  The inverse transformation
%% maps these two images to the result.  This is depicted graphically
%% with a system diagram:
%%
%% IM ---> flo/down2 --> LO1 --> up2/flo --> add --> RECON
%%     |                                      ^
%%     |	                              |
%%     |	                              |
%%      -> fhi/down2 --> HI1 --> up2/fhi ----- 
%% 
%% Note that the number of samples in the representation (i.e., total
%% samples in LO1 and HI1) is equal to the number of samples in the
%% original IM.  Thus, this representation is exactly COMPLETE, or
%% "critically sampled".

%% So we've fixed one of the problems that we had with Laplacian
%% pyramid.  But the system diagram above places strong constraints on
%% the filters.  In particular, for these filters the reconstruction
%% is no longer perfect.  Turns out there are NO
%% perfect-reconstruction symmetric filters that are
%% power-complementary, except for the trivial case [1] and the
%% nearly-trivial case [1 1]/sqrt(2).

%% Let's consider the projection functions of this 2-band splitting
%% operation.  We can construct these by applying the transform to
%% impulse input signals, for all possible impulse locations.  The
%% rows of the following matrix are the projection functions for each
%% coefficient in the transform.
M = [corrDn(eye(32),flo','circular',[1 2]), ...
     corrDn(eye(32),fhi','circular',[1 2],[1 2])]';
clf; showIm(M,'auto1','auto','M');

%% The transform matrix is composed of two sub-matrices.  The top half
%% contains the lowpass kernel, shifted by increments of 2 samples.
%% The bottom half contains the highpass.  Now we compute the inverse
%% of this matrix: 
M_inv = inv(M);
showIm(M_inv,'auto1','auto','M_inv');

%% The inverse is (very close to) the transpose of the original
%% matrix!  In other words, the transform is orthonormal.
imStats(M_inv',M);

%% This also points out a nice relationship between the corrDn and
%% upConv functions, and the matrix representation.  corrDn is
%% equivalent to multiplication by a matrix with copies of the filter
%% on the ROWS, translated in multiples of the downsampling factor.
%% upConv is equivalent to multiplication by a matrix with copies of
%% the filter on the COLUMNS, translated by the upsampling factor.

%% As in the Laplacian pyramid, we can recursively apply this QMF 
%% band-splitting operation to the lowpass band:
lo2 = corrDn(lo1,flo,'reflect1',[1 2]);
hi2 = corrDn(lo1,fhi,'reflect1',[1 2],[1 2]);

%% The representation of the original signal is now comprised of the
%% three subbands {hi1, hi2, lo2} (we don't hold onto lo1, because it
%% can be reconstructed from lo2 and hi2).  Note that hi1 is at 1/2
%% resolution, and hi2 and lo2 are at 1/4 resolution: The total number
%% of samples in these three subbands is thus equal to the number of
%% samples in the original signal.
imnames=['hi1'; 'hi2'; 'lo2'];
for bnum=1:3
  band = eval(imnames(bnum,:));
  subplot(3,1,bnum); showIm(band); ylabel(imnames(bnum,:));
  axis([1 size(band,2) 1.1*min(lo2) 1.1*max(lo2)]);
end

%% Reconstruction proceeds as with the Laplacian pyramid: combine lo2 and hi2
%% to reconstruct lo1, which is then combined with hi1 to reconstruct the
%% original signal:
recon_lo1 = upConv(hi2,fhi,'reflect1',[1 2],[1 2]) + ...
            upConv(lo2,flo,'reflect1',[1 2],[1 1]);
reconstructed = upConv(hi1,fhi,'reflect1',[1 2],[1 2]) + ...
                upConv(recon_lo1,flo,'reflect1',[1 2],[1 1]);
imStats(sig,reconstructed);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% FUNCTIONS for CONSTRUCTING/MANIPULATING QMF/Wavelet PYRAMIDS

%% To make things easier, we have bundled these qmf operations and
%% data structures into an object in MATLAB.

sig = mkFract([1 64], 1.5);
[pyr,pind] = buildWpyr(sig);
showWpyr(pyr,pind);

nbands = size(pind,1);
for b = 1:nbands
  subplot(nbands,1,b); lplot(pyrBand(pyr,pind,b));
end
	
res = reconWpyr(pyr,pind);
imStats(sig,res);

%% Now for 2D, we use separable filters.  There are 4 ways to apply the two 
%% filters to the input image (followed by the relavent subsampling operation):
%%   (1) lowpass in both x and y
%%   (2) lowpass in x and highpass in y 
%%   (3) lowpass in y and highpass in x
%%   (4) highpass in both x and y.  
%% The pyramid is built by recursively subdividing the first of these bands
%% into four new subbands.

%% First, we'll take a look at some of the basis functions.
sz = 40;
zim = zeros(sz);
flo = 'qmf9'; edges = 'reflect1';
[pyr,pind] = buildWpyr(zim);

% Put an  impulse into the middle of each band:
for lev=1:size(pind,1)
  mid = sum(prod(pind(1:lev-1,:)'));
  mid = mid + floor(pind(lev,2)/2)*pind(lev,1) + floor(pind(lev,1)/2) + 1;
  pyr(mid,1) = 1;
end

% And take a look at the reconstruction of each band:
for lnum=1:wpyrHt(pind)+1
  for bnum=1:3
    subplot(wpyrHt(pind)+1,3,(wpyrHt(pind)+1-lnum)*3+bnum);
    showIm(reconWpyr(pyr, pind, flo, edges, lnum, bnum),'auto1',2,0);
  end
end

%% Note that the first column contains horizontally oriented basis functions at
%% different scales.  The second contains vertically oriented basis functions.
%% The third contains both diagonals (a checkerboard pattern).  The bottom row
%% shows (3 identical images of) a lowpass basis function.

%% Now look at the corresponding Fourier transform magnitudes (these
%% are plotted over the frequency range [-pi, pi] ):
nextFig(2,1);
freq = 2 * pi * [-sz/2:(sz/2-1)]/sz;
for lnum=1:wpyrHt(pind)+1
  for bnum=1:3
    subplot(wpyrHt(pind)+1,3,(wpyrHt(pind)+1-lnum)*3+bnum);
    basisFn = reconWpyr(pyr, pind, flo, edges, lnum, bnum);
    basisFmag = fftshift(abs(fft2(basisFn,sz,sz)));
    imagesc(freq,freq,basisFmag);
    axis('square'); axis('xy'); colormap('gray');
  end
end
nextFig(2,-1);

%% The filters at a given scale sum to a squarish annular region:
sumSpectra = zeros(sz);
lnum = 2;
for bnum=1:3
  basisFn = reconWpyr(pyr, pind, flo, edges, lnum, bnum);
  basisFmag = fftshift(abs(fft2(basisFn,sz,sz)));
  sumSpectra = basisFmag.^2 + sumSpectra;
end
clf; imagesc(freq,freq,sumSpectra); axis('square'); axis('xy'); title('one scale');

%% Now decompose an image:
[pyr,pind] = buildWpyr(im);

%% View all of the subbands (except lowpass), scaled to be the same size
%% (requires a big figure window):
nlevs = wpyrHt(pind);
for lnum=1:nlevs
  for bnum=1:3
    subplot(nlevs,3,(lnum-1)*3+bnum); 
    showIm(wpyrBand(pyr,pind,lnum,bnum), 'auto2', 2^(lnum+imSubSample-2));
  end
end

%% In addition to the bands shown above, there's a lowpass residual:
nextFig(2,1);
clf; showIm(pyrLow(pyr,pind));
nextFig(2,-1);

% Alternatively, display the pyramid with the subbands shown at their
% correct relative sizes:
clf; showWpyr(pyr, pind);

%% The reconWpyr function can be used to reconstruct the entire pyramid:
reconstructed = reconWpyr(pyr,pind);
imStats(im,reconstructed);

%% As with Laplacian pyramids, you can specify sub-levels and subbands
%% to be included in the reconstruction.  For example:
clf
showIm(reconWpyr(pyr,pind,'qmf9','reflect1',[1:wpyrHt(pind)],[1]));  %Horizontal only
showIm(reconWpyr(pyr,pind,'qmf9','reflect1',[2,3])); %two middle scales

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% PERFECT RECONSTRUCTION: HAAR AND DEBAUCHIES WAVELETS

%% The symmetric QMF filters used above are not perfectly orthogonal.
%% In fact, it's impossible to construct a symmetric filter of size
%% greater than 2 that is perfectly orthogonal to shifted copies
%% (shifted by multiples of 2) of itself.  For example, consider a
%% symmetric kernel of length 3.  Shift by two and the right end of