pyramids.m

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

%% the original kernel is aligned with the left end of the shifted
%% one.  Thus, the inner product of these two will be the square of
%% the end tap, which will be non-zero.

%% However, one can easily create wavelet filters of length 2 that
%% will do the job.  This is the oldest known wavelet, known as the
%% "Haar".  The two kernels are [1,1]/sqrt(2) and [1,-1]/sqrt(2).
%% These are trivially seen to be orthogonal to each other, and shifts
%% by multiples of two are also trivially orthogonal.  The projection
%% functions of the Haar transform are in the rows of the following
%% matrix, constructed by applying the transform to impulse input
%% signals, for all possible impulse locations:

haarLo = namedFilter('haar')
haarHi = modulateFlip(haarLo)
subplot(2,1,1); lplot(haarLo); axis([0 3 -1 1]); title('lowpass');
subplot(2,1,2); lplot(haarHi); axis([0 3 -1 1]); title('highpass');

M = [corrDn(eye(32), haarLo, 'reflect1', [2 1], [2 1]); ...
    corrDn(eye(32), haarHi, 'reflect1', [2 1], [2 1])];
clf; showIm(M)
showIm(M*M') %identity!

%% As before, the filters are power-complementary (although the
%% frequency isolation is rather poor, and thus the subbands will be
%% heavily aliased):
plot(pi*[-32:31]/32,abs(fft(haarLo,64)).^2,'--',...
     pi*[-32:31]/32,abs(fft(haarHi,64)).^2,'-');

sig = mkFract([1,64],0.5);
[pyr,pind] = buildWpyr(sig,4,'haar','reflect1');
showWpyr(pyr,pind);

%% check perfect reconstruction:
res = reconWpyr(pyr,pind, 'haar', 'reflect1');
imStats(sig,res)

%% If you want perfect reconstruction, but don't like the Haar
%% transform, there's another option: drop the symmetry requirement.
%% Ingrid Daubechies developed one of the earliest sets of such
%% perfect-reconstruction wavelets.  The simplest of these is of
%% length 4:

daub_lo = namedFilter('daub2');
daub_hi = modulateFlip(daub_lo);

%% The daub_lo filter is constructed to be orthogonal to 2shifted
%% copy of itself.  For example:
[daub_lo;0;0]'*[0;0;daub_lo]

M = [corrDn(eye(32), daub_lo, 'circular', [2 1], [2 1]); ...
    corrDn(eye(32), daub_hi, 'circular', [2 1], [2 1])];
clf; showIm(M)
showIm(M*M') % identity!

%% Again, they're power complementary:
plot(pi*[-32:31]/32,abs(fft(daub_lo,64)).^2,'--',...
     pi*[-32:31]/32,abs(fft(daub_hi,64)).^2,'-');
 
%% The sum of the power spectra is again flat
plot(pi*[-32:31]/32,...
    fftshift(abs(fft(daub_lo,64)).^2)+fftshift(abs(fft(daub_hi,64)).^2));

%% Make a pyramid using the same code as before (except that we can't
%% use reflected boundaries with asymmetric filters):
[pyr,pind] = buildWpyr(sig, maxPyrHt(size(sig),size(daub_lo)), daub_lo, 'circular');
showWpyr(pyr,pind,'indep1');

res = reconWpyr(pyr,pind, daub_lo,'circular');
imStats(sig,res);

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

%%% ALIASING IN WAVELET TRANSFORMS

%% All of these orthonormal pyramid/wavelet transforms have a lot
%% of aliasing in the subbands.  You can see that in the frequency
%% response plots since the frequency response of each filter
%% covers well more than half the frequency domain.  The aliasing
%% can have serious consequences...

%% Get one of the basis functions of the 2D Daubechies wavelet transform:
[pyr,pind] = buildWpyr(zeros(1,64),4,daub_lo,'circular');
lev = 3;
pyr(1+sum(pind(1:lev-1,2))+pind(lev,2)/2,1) = 1;
sig = reconWpyr(pyr,pind, daub_lo,'circular');
clf; lplot(sig)

%% Since the basis functions are orthonormal, building a pyramid using this
%% input will yield a single non-zero coefficient.
[pyr,pind] = buildWpyr(sig, 4, daub_lo, 'circular');
figure(1);
nbands = size(pind,1)
for b=1:nbands
  subplot(nbands,1,b); lplot(pyrBand(pyr,pind,b));
  axis([1 size(pyrBand(pyr,pind,b),2) -0.3 1.3]);
end

%% Now shift the input by one sample and re-build the pyramid.
shifted_sig = [0,sig(1:size(sig,2)-1)];
[spyr,spind] = buildWpyr(shifted_sig, 4, daub_lo, 'circular');

%% Plot each band of the unshifted and shifted decomposition
nextFig(2);
nbands = size(spind,1)
for b=1:nbands
  subplot(nbands,1,b); lplot(pyrBand(spyr,spind,b));
  axis([1 size(pyrBand(spyr,spind,b),2) -0.3 1.3]);
end
nextFig(2,-1);

%% In the third band, we expected the coefficients to move around
%% because the signal was shifted.  But notice that in the original
%% signal decomposition, the other bands were filled with zeros.
%% After the shift, they have significant content.  Although these
%% subbands are supposed to represent information at different scales,
%% their content also depends on the relative POSITION of the input
%% signal.

%% This problem is not unique to the Daubechies transform.  The same
%% is true for the QMF transform.  Try it...  In fact, the same kind
%% of problem occurs for almost any orthogonal pyramid transform (the
%% only exception is the limiting case in which the filter is a sinc
%% function).

%% Orthogonal pyramid transforms are not shift-invariant.  Although
%% orthogonality may be an important property for some applications
%% (e.g., data compression), orthogonal pyramid transforms are
%% generally not so good for image analysis.

%% The overcompleteness of the Laplacian pyramid turns out to be a
%% good thing in the end.  By using an overcomplete representation
%% (and by choosing the filters properly to avoid aliasing as much as
%% possible), you end up with a representation that is useful for
%% image analysis.

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

%%% The "STEERABLE PYRAMID" 

%% The steerable pyramid is a multi-scale representation that is
%% translation-invariant, but that also includes representation of
%% orientation.  Furthermore, the representation of orientation is
%% designed to be rotation-invariant. The basis/projection functions
%% are oriented (steerable) filters, localized in space and frequency.
%% It is overcomplete to avoid aliasing.  And it is "self-inverting"
%% (like the QMF/Wavelet transform): the projection functions and 
%% basis functions are identical.  The mathematical phrase for a 
%% transform obeying this property is "tight frame".

%% The system diagram for the steerable pyramid (described in the
%% reference given below) is as follows:
%
% IM ---> fhi0 -----------------> H0 ---------------- fhi0 ---> RESULT
%     |                                                     |
%     |                                                     |
%     |-> flo0 ---> fl1/down2 --> L1 --> up2/fl1 ---> flo0 -|
%               |                                 |
%               |----> fb0 -----> B0 ----> fb0 ---|
%               |                                 |
%               |----> fb1 -----> B1 ----> fb1 ---|
%               .                                 .
%               .                                 .
%               |----> fbK -----> BK ----> fbK ---|
%
%% The filters {fhi0,flo0} are used to initially split the image into
%% a highpass residual band H0 and a lowpass subband.  This lowpass
%% band is then split into a low(er)pass band L1 and K+1 oriented
%% subbands {B0,B1,...,BK}.  The representatation is substantially
%% overcomplete.  The pyramid is built by recursively splitting the
%% lowpass band (L1) using the inner portion of the diagram (i.e.,
%% using the filters {fl1,fb0,fb1,...,fbK}).  The resulting transform is
%% overcomplete by a factor of 4k/3.

%% The scale tuning of the filters is constrained by the recursive
%% system diagram.  The orientation tuning is constrained by requiring
%% the property of steerability.  A set of filters form a steerable
%% basis if they 1) are rotated copies of each other, and 2) a copy of
%% the filter at any orientation may be computed as a linear
%% combination of the basis filters.  The simplest examples of
%% steerable filters is a set of N+1 Nth-order directional
%% derivatives.

%% Choose a filter set (options are 'sp0Filters', 'sp1Filters',
%% 'sp3Filters', 'sp5Filters'):
filts = 'sp3Filters';
[lo0filt,hi0filt,lofilt,bfilts,steermtx,harmonics] = eval(filts);
fsz = round(sqrt(size(bfilts,1))); fsz =  [fsz fsz];
nfilts = size(bfilts,2);
nrows = floor(sqrt(nfilts));

%% Look at the oriented bandpass filters:
for f = 1:nfilts
  subplot(nrows,ceil(nfilts/nrows),f);
  showIm(conv2(reshape(bfilts(:,f),fsz),lo0filt));
end

%% Try "steering" to a new orientation (new_ori in degrees):
new_ori = 360*rand(1)
clf; showIm(conv2(reshape(steer(bfilts, new_ori*pi/180 ), fsz), lo0filt));

%% Look at Fourier transform magnitudes:
lo0 = fftshift(abs(fft2(lo0filt,64,64)));
fsum = zeros(size(lo0));
for f = 1:size(bfilts,2)
  subplot(nrows,ceil(nfilts/nrows),f);
  flt = reshape(bfilts(:,f),fsz);
  freq = lo0 .* fftshift(abs(fft2(flt,64,64)));
  fsum = fsum + freq.^2;
  showIm(freq);
end

%% The filters sum to a smooth annular ring:
clf; showIm(fsum);

%% build a Steerable pyramid:
[pyr,pind] = buildSpyr(im, 4-imSubSample, filts);
 
%% Look at first (vertical) bands, different scales:
for s = 1:min(4,spyrHt(pind))
  band = spyrBand(pyr,pind,s,1);
  subplot(2,2,s); showIm(band);
end

%% look at all orientation bands at one level (scale):
for b = 1:spyrNumBands(pind)
  band = spyrBand(pyr,pind,1,b);
  subplot(nrows,ceil(nfilts/nrows),b);
  showIm(band);
end

%% To access the high-pass and low-pass bands:
low = pyrLow(pyr,pind);
showIm(low);
high = spyrHigh(pyr,pind);
showIm(high);

%% Display the whole pyramid (except for the highpass residual band),
%% with images shown at proper relative sizes:
showSpyr(pyr,pind);

%% Spin a level of the pyramid, interpolating (steering to)
%% intermediate orienations:

[lev,lind] = spyrLev(pyr,pind,2);
lev2 = reshape(lev,prod(lind(1,:)),size(bfilts,2));
figure(1); subplot(1,1,1); showIm(spyrBand(pyr,pind,2,1));
M = moviein(16);
for frame = 1:16
  steered_im = steer(lev2, 2*pi*(frame-1)/16, harmonics, steermtx);
  showIm(reshape(steered_im, lind(1,:)),'auto2');
  M(:,frame) = getframe;
end

%% Show the movie 3 times:
movie(M,3);

%% Reconstruct.  Note that the filters are not perfect, although they are good
%% enough for most applications.
res = reconSpyr(pyr, pind, filts); 
showIm(im + i * res);
imStats(im,res);

%% As with previous pyramids, you can select subsets of the levels
%% and orientation bands to be included in the reconstruction.  For example:

%% All levels (including highpass and lowpass residuals), one orientation:
clf; showIm(reconSpyr(pyr,pind,filts,'reflect1','all', [1]));

%% Without the highpass and lowpass:
clf; showIm(reconSpyr(pyr,pind,filts,'reflect1',[1:spyrHt(pind)], [1]));

%% We also provide an implementation of the Steerable pyramid in the
%% Frequency domain.  The advantages are perfect-reconstruction
%% (within floating-point error), and any number of orientation
%% bands.  The disadvantages are that it is typically slower, and the
%% boundary handling is always circular.

[pyr,pind] = buildSFpyr(im,4,4); % 4 levels, 5 orientation bands
showSpyr(pyr,pind);
res = reconSFpyr(pyr,pind);
imStats(im,res);  % nearly perfect

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The steerable pyramid transform given above is described in:
%
%   E P Simoncelli and W T Freeman. 
%   The Steerable Pyramid: A Flexible Architecture for Multi-Scale 
%   Derivative Computation.  IEEE Second Int'l Conf on Image Processing. 
%   Washington DC,  October 1995.
%
% Online access:
% Abstract:  http://www.cis.upenn.edu/~eero/ABSTRACTS/simoncelli95b-abstract.html
% Full (PostScript):  ftp://ftp.cis.upenn.edu/pub/eero/simoncelli95b.ps.Z
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Local Variables:
%% buffer-read-only: t 
%% End: