1 files changed, 0 insertions, 903 deletions

diff --git a/SD-VBS/benchmarks/texture_synthesis/src/matlabPyrTools/TUTORIALS/pyramids.m b/SD-VBS/benchmarks/texture_synthesis/src/matlabPyrTools/TUTORIALS/pyramids.m
deleted file mode 100755
index 2123c69..0000000
--- a/SD-VBS/benchmarks/texture_synthesis/src/matlabPyrTools/TUTORIALS/pyramids.m
+++ /dev/null
@@ -1,903 +0,0 @@
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%%% 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,
-%% 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
-%% 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: