From 0efc775370d2ff91927d1b383a99eab78dc5538f Mon Sep 17 00:00:00 2001 From: leochanj Date: Wed, 21 Oct 2020 01:52:54 -0400 Subject: debug libextra and remove matlab FLUSH_CACHES --- .../src/matlabPyrTools/TUTORIALS/pyramids.m | 903 --------------------- 1 file changed, 903 deletions(-) delete mode 100755 SD-VBS/benchmarks/texture_synthesis/src/matlabPyrTools/TUTORIALS/pyramids.m (limited to 'SD-VBS/benchmarks/texture_synthesis/src/matlabPyrTools/TUTORIALS/pyramids.m') 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 " 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: -- cgit v1.2.2