pyramids.m

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% IMAGE PYRAMID TUTORIAL 
%%%
%%% A brief introduction to multi-scale pyramid decompositions for image 
%%% processing.  You should go through this, reading the comments, and
%%% executing the corresponding MatLab instructions.  This file assumes 
%%% a basic familiarity with matrix algebra, with linear systems and Fourier
%%% theory, and with MatLab.  If you don't understand a particular
%%% function call, execute "help <functionName>" to see documentation.
%%%
%%% EPS, 6/96.  
%%% Based on the original OBVIUS tutorial.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Determine a subsampling factor for images, based on machine speed:
oim = pgmRead('einstein.pgm');
tic; corrDn(oim,[1 1; 1 1]/4,'reflect1',[2 2]); time = toc;
imSubSample = min(max(floor(log2(time)/2+3),0),2);
im = blurDn(oim, imSubSample,'qmf9');
clear oim;
clf; showIm(im, 'auto2', 'auto', 'im');

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

%%% LAPLACIAN PYRAMIDS: 

%% Images may be decomposed into information at different scales.  
%% Blurring eliminates the fine scales (detail):

binom5 = binomialFilter(5);
lo_filt = binom5*binom5';
blurred = rconv2(im,lo_filt);
subplot(1,2,1); showIm(im, 'auto2', 'auto', 'im');
subplot(1,2,2); showIm(blurred, 'auto2', 'auto', 'blurred');

%% Subtracting the blurred image from the original leaves ONLY the
%% fine scale detail:
fine0 = im - blurred;
subplot(1,2,1); showIm(fine0, 'auto2', 'auto', 'fine0');

%% The blurred and fine images contain all the information found in
%% the original image.  Trivially, adding the blurred image to the
%% fine scale detail will reconstruct the original.  We can compare
%% the original image to the sum of blurred and fine using the
%% "imStats" function, which reports on the statistics of the
%% difference between it's arguments:
imStats(im, blurred+fine0);

%% Since the filter is a lowpass filter, we might want to subsample
%% the blurred image.  This may cause some aliasing (depends on the
%% filter), but the decomposition structure given above will still be
%% possible.  The corrDn function correlates (same as convolution, but
%% flipped filter) and downsamples in a single operation (for
%% efficiency).  The string 'reflect1' tells the function to handle
%% boundaries by reflecting the image about the edge pixels.  Notice
%% that the blurred1 image is half the size (in each dimension) of the
%% original image.
lo_filt = 2*binom5*binom5';  %construct a separable 2D filter
blurred1 = corrDn(im,lo_filt,'reflect1',[2 2]);
subplot(1,2,2); showIm(blurred1,'auto2','auto','blurred1');

%% Now, to extract fine scale detail, we must interpolate the image
%% back up to full size before subtracting it from the original.  The
%% upConv function does upsampling (padding with zeros between
%% samples) followed by convolution.  This can be done using the
%% lowpass filter that was applied before subsampling or it can be
%% done with a different filter.
fine1 = im - upConv(blurred1,lo_filt,'reflect1',[2 2],[1 1],size(im));
subplot(1,2,1); showIm(fine1,'auto2','auto','fine1');

%% We now have a technique that takes an image, computes two new
%% images (blurred1 and fine1) containing the coarse scale information
%% and the fine scale information.  We can also (trivially)
%% reconstruct the original from these two (even if the subsampling of
%% the blurred1 image caused aliasing):

recon = fine1 + upConv(blurred1,lo_filt,'reflect1',[2 2],[1 1],size(im));
imStats(im, recon);

%% Thus, we have described an INVERTIBLE linear transform that maps an
%% input image to the two images blurred1 and fine1.  The inverse
%% transformation maps blurred1 and fine1 to the result.  This is
%% depicted graphically with a system diagram:
%%
%% IM --> blur/down2 ---------> BLURRED1 --> up2/blur --> add --> RECON
%%  |                   |                                  ^
%%  |	                |                                  |
%%  |	                V                                  |
%%  |	             up2/blur                              |
%%  |	                |                                  |
%%  |	                |                                  |
%%  |	                V                                  |
%%   --------------> subtract --> FINE1 -------------------
%% 
%% Note that the number of samples in the representation (i.e., total
%% samples in BLURRED1 and FINE1) is 1.5 times the number of samples
%% in the original IM.  Thus, this representation is OVERCOMPLETE.

%% Often, we will want further subdivisions of scale.  We can
%% decompose the (coarse-scale) BLURRED1 image into medium coarse and
%% very coarse images by applying the same splitting technique:
blurred2 = corrDn(blurred1,lo_filt,'reflect1',[2 2]);
showIm(blurred2)

fine2 = blurred1 - upConv(blurred2,lo_filt,'reflect1',[2 2],[1 1],size(blurred1));
showIm(fine2)

%% Since blurred2 and fine2 can be used to reconstruct blurred1, and
%% blurred1 and fine1 can be used to reconstruct the original image,
%% the set of THREE images (also known as "subbands") {blurred2,
%% fine2, fine1} constitute a complete representation of the original
%% image.  Note that the three subbands are displayed at the same size,
%% but they are actually three different sizes.

subplot(1,3,1); showIm(fine1,'auto2',2^(imSubSample-1),'fine1');
subplot(1,3,2); showIm(fine2,'auto2',2^(imSubSample),'fine2');
subplot(1,3,3); showIm(blurred2,'auto2',2^(imSubSample+1),'blurred2');

%% It is useful to consider exactly what information is stored in each
%% of the pyramid subbands.  The reconstruction process involves
%% recursively interpolating these images and then adding them to the
%% image at the next finer scale.  To see the contribution of ONE of
%% the representation images (say blurred2) to the reconstruction, we
%% imagine filling all the other subbands with zeros and then
%% following our reconstruction procedure.  For the blurred2 subband,
%% this is equivalent to simply calling upConv twice:
blurred2_full = upConv(upConv(blurred2,lo_filt,'reflect1',[2 2],[1 1],size(blurred1)),...
    lo_filt,'reflect1',[2 2],[1 1],size(im));
subplot(1,3,3); showIm(blurred2_full,'auto2',2^(imSubSample-1),'blurred2-full');

%% For the fine2 subband, this is equivalent to calling upConv once:
fine2_full = upConv(fine2,lo_filt,'reflect1',[2 2],[1 1],size(im));
subplot(1,3,2); showIm(fine2_full,'auto2',2^(imSubSample-1),'fine2-full');

%% If we did everything correctly, we should be able to add together
%% these three full-size images to reconstruct the original image:
recon = blurred2_full + fine2_full + fine1;
imStats(im, recon)

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

%%% FUNCTIONS for CONSTRUCTING/MANIPULATING LAPLACIAN PYRAMIDS

%% We can continue this process, recursively splitting off finer and
%% finer details from the blurred image (like peeling off the outer
%% layers of an onion).  The resulting data structure is known as a
%% "Laplacian Pyramid".  To make things easier, we have written a
%% MatLab function called buildLpyr to construct this object.  The
%% function returns two items: a long vector containing the subbands
%% of the pyramid, and an index matrix that is used to access these
%% subbands.  The display routine showLpyr shows all the subbands of the
%% pyramid, at the their correct relative sizes.  It should now be
%% clearer why these data structures are called "pyramids".
[pyr,pind] = buildLpyr(im,5-imSubSample); 
showLpyr(pyr,pind);

%% There are also "accessor" functions for pulling out a single subband:
showIm(pyrBand(pyr,pind,2));

%% The reconLpyr function allows you to reconstruct from a laplacian pyramid.
%% The third (optional) arg allows you to select any subset of pyramid bands
%% (default is to use ALL of them).
clf; showIm(reconLpyr(pyr,pind,[1 3]),'auto2','auto','bands 1 and 3 only');

fullres = reconLpyr(pyr,pind);
showIm(fullres,'auto2','auto','Full reconstruction');
imStats(im,fullres);

%% buildLpyr uses 5-tap filters by default for building Laplacian
%% pyramids.  You can specify other filters:
namedFilter('binom3')
[pyr3,pind3] = buildLpyr(im,5-imSubSample,'binom3');
showLpyr(pyr3,pind3);
fullres3 = reconLpyr(pyr3,pind3,'all','binom3');
imStats(im,fullres3);

%% Here we build a "Laplacian" pyramid using random filters.  filt1 is
%% used with the downsampling operations and filt2 is used with the
%% upsampling operations.  We normalize the filters for display
%% purposes.  Of course, these filters are (almost certainly) not very
%% "Gaussian", and the subbands of such a pyramid will be garbage!
%% Nevertheless, it is a simple property of the Laplacian pyramid that
%% we can use ANY filters and we will still be able to reconstruct
%% perfectly.

filt1 = rand(1,5); filt1 = sqrt(2)*filt1/sum(filt1)
filt2 = rand(1,3); filt2 = sqrt(2)*filt2/sum(filt2)
[pyrr,pindr] = buildLpyr(im,5-imSubSample,filt1,filt2);
showLpyr(pyrr,pindr);
fullresr = reconLpyr(pyrr,pindr,'all',filt2);
imStats(im,fullresr);

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

%%% ALIASING in the Gaussian and Laplacian pyramids:

%% Unless one is careful, the subsampling operations will introduce aliasing
%% artifacts in these pyramid transforms.  This is true even though the
%% Laplacian pyramid can be used to reconstruct the original image perfectly.
%% When reconstructing, the pyramid is designed in such a way that these
%% aliasing artifacts cancel out.  So it's not a problem if the only thing we
%% want to do is reconstruct.  However, it can be a serious problem if we
%% intend to process each of the subbands independently.

%% One way to see the consequences of the aliasing artifacts is by
%% examining variations that occur when the input is shifted.  We
%% choose an image and shift it by some number of pixels.  Then blur
%% (filter-downsample-upsample-filter) the original image and blur the
%% shifted image.  If there's no aliasing, then the blur and shift
%% operations should commute (i.e.,
%% shift-filter-downsample-upsample-filter is the same as
%% filter-downsample-upsample-filter-shift).  Try this for 2 different
%% filters (by replacing 'binom3' with 'binom5' or 'binom7' below),
%% and you'll see that the aliasing is much worse for the 3 tap
%% filter.

sig = 100*randn([1 16]);
sh = [0 7];  %shift amount
lev = 2; % level of pyramid to look at
flt = 'binom3';  %filter to use: 

shiftIm = shift(sig,sh);
[pyr,pind] = buildLpyr(shiftIm, lev, flt, flt, 'circular');
shiftBlur = reconLpyr(pyr, pind, lev, flt, 'circular');

[pyr,pind] = buildLpyr(sig, lev, flt, flt, 'circular');
res = reconLpyr(pyr, pind, lev, flt, 'circular');
blurShift = shift(res,sh);

subplot(2,1,1); r = showIm(shiftBlur,'auto2','auto','shiftBlur');
subplot(2,1,2); showIm(blurShift,r,'auto','blurShift');
imStats(blurShift,shiftBlur);

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

%%% PROJECTION and BASIS functions:

%% An invertible, linear transform can be characterized in terms
%% of a set of PROJECTION and BASIS functions.  In matlab matrix
%% notation:
%
%%     c = P' * x
%%     x = B * c
%
%% where x is an input, c are the transform coefficients, P and B
%% are matrices.  The columns of P are the projection functions (the
%% input is projected onto the the columns of P to get each successive
%% transform coefficient).  The columns of B are the basis
%% functions (x is a linear combination of the columns of B).

%% Since the Laplacian pyramid is a linear transform, we can ask: what
%% are its BASIS functions?  We consider these in one dimension for
%% simplicity.  The BASIS function corresponding to a given
%% coefficient tells us how much that coefficient contributes to each
%% pixel in the reconstructed image.  We can construct a single basis
%% function by setting one sample of one subband equal to 1.0 (and all
%% others to zero) and reconstructing. To build the entire matrix, we
%% have to do this for every sample of every subband:
sz = min(round(48/(sqrt(2)^imSubSample)),36);
sig = zeros(sz,1);
[pyr,pind] = buildLpyr(sig);
basis = zeros(sz,size(pyr,1));
for n=1:size(pyr,1)
  pyr = zeros(size(pyr));
  pyr(n) = 1;
  basis(:,n) = reconLpyr(pyr,pind);
end
clf; showIm(basis)

%% The columns of the basis matrix are the basis functions.  The
%% matrix is short and fat, corresponding to the fact that the
%% representation is OVERCOMPLETE.  Below, we plot the middle one from
%% each subband, starting with the finest scale.  Note that all of
%% these basis functions are lowpass (Gaussian-like) functions.
locations = round(sz * (2 - 3./2.^[1:max(4,size(pind,1))]))+1;
for lev=1:size(locations,2)
  subplot(2,2,lev);
  showIm(basis(:,locations(lev)));
  axis([0 sz 0 1.1]);
end

%% Now, we'd also like see the inverse (we'll them PROJECTION)
%% functions. We need to ask how much of each sample of the input
%% image contributes to a given pyramid coefficient.  Thus, the matrix
%% is constructed by building pyramids on the set of images with
%% impulses at each possible location.  The rows of this matrix are
%% the projection functions.
projection = zeros(size(pyr,1),sz);
for pos=1:sz
  [pyr,pind] = buildLpyr(mkImpulse([1 sz], [1 pos]));
  projection(:,pos) = pyr;
end
clf; showIm(projection);

%% Building a pyramid corresponds to multiplication by the projection
%% matrix.  Reconstructing from this pyramid corresponds to
%% multiplication by the basis matrix.  Thus, the product of the two
%% matrices (in this order) should be the identity matrix:
showIm(basis*projection);

%% We can plot a few example projection functions at different scales.
%% Note that all of the projection functions are bandpass functions,