Do, Minh N.

Vetterli, Martin

Infoscience - École polytechnique fédérale de Lausanne

FRAME RECONSTRUCTION OF THE LAPLACIAN PYRAMIDMinh N. Doy and Martin VetterliyxyLaboratory for Audio-Visual CommunicationsSwiss Federal Institute of Technology Lausanne (EPFL), 1015 Lausanne, SwitzerlandxDepartment of EECS, University of California at Berkeley, Berkeley CA 94720, USAEmail: fMinh.Do, Martin.Vetterlig@epfl.ch; Web: lcavwww.epfl.chProc. of IEEE Intl. Conf. on Acoustics, Speech, and Signal Processing, Salt Lake City, 2001ABSTRACTWe study the Laplacian pyramid (LP) as a frame operator, and thisreveals that the usual reconstruction is suboptimal. With orthog-onal filters, the LP is shown to be a tight frame, thus the optimallinear reconstruction using the dual frame operator has a simplestructure as symmetrical with the forward transform. For moregeneral cases, we propose an efficient filter bank for reconstructionin the LP that is shown to perform better than the usual method.Numerical results indicate that gains of more than 1 dB are actu-ally achieved.1. INTRODUCTIONMultiscale data representation is a powerful idea. In particular,it captures data in hierarchical structures where each level corre-sponds to a reduced-resolution approximation. One of the earlyexamples of such a scheme is the Laplacian pyramid (LP), pro-posed by Burt and Adelson [1] for image coding. The basic ideaof the LP is the following. First, derive a coarse approximation ofthe original signal by lowpass filtering and downsampling. Basedon this coarse version, predict the original (by upsampling and fil-tering) and then calculate the difference as the prediction error.The process can be iterated on the coarse version. A diagram forthe LP is shown in Fig. 1.+_M MH Gxx0y0y1Fig. 1. Laplacian pyramid scheme.A drawback of the LP is the implicit oversampling. Therefore,in compression applications it is normally replaced by subbandcoding or wavelet transform which is a critically sampled schemeand often an orthogonal decomposition. However, the LP has theadvantage over the critically sampled wavelet scheme in the sensethat each pyramid level generates only one bandpass image (evenfor multidimensional cases) which does not have “scrambled” fre-quencies. This frequency scrambling happens in the wavelet filterbank when a highpass channel after downsampling is folded backThis work was supported by a Department of Communication Systems- EPFL PhD Fellowship and the Swiss National Science Foundation grantnumber 21-52439.97.into the low frequency, and thus its spectrum is reflected. In theLP, this effect is avoided by only downsampling the lowpass chan-nel. Therefore the LP permits further subband decomposition to beapplied on its bandpass images. A possible scheme is a pyramidaldecomposition where bandpass images of the LP are fed into di-rectional filter banks. The end result is a set of directional subbandimages at multiple scales [2].For many applications like compression and denoising, the co-efficients in the transform domain are processed further and thiscan introduce errors due to quantization or thresholding. The pro-cessed coefficients are then used to reconstruct the original data.For the LP, the usual reconstruction – adding the prediction fromthe coarse version with the difference, produces perfect recon-structed image in the absence of noise but turns out to be sub-optimal otherwise.Our key observation is that the LP is a frame operator and thatone should use the dual frame operator for reconstruction. Whilethis seems a somewhat trivial observation, it has not been used inpractice, probably because the usual reconstruction, while subop-timal, is so simple. Yet, we will show that gains of more than 1 dBare actually possible.This paper is organized as follows. Section 2 introduces thenotations and sets up the LP in the frame analysis where the dualframe operator or the pseudo-inverse is defined. Section 3 consid-ers an important case where the LP is in fact a tight frame, andthus the dual frame operator has a simple structure as symmetricalwith the forward transform. In Section 4, inspired by the structureof the new reconstruction algorithm for the tight frame case, westudy a more general LP where it has better inverse operator withan efficient filter bank. Section 5 compares the two reconstructionalgorithms for multilevel LP’s. The performance improvement ofthe new reconstruction algorithm is illustrated with experimentalresults in Section 6.2. PRELIMINARIES2.1. Burt and Adelson’s Laplacian PyramidWe consider the general d-dimensional case where a discrete-timesignal is a sequence of real-valued numbers defined on the inte-ger lattice Zd, e.g. x[n]; n 2 Zd. The sampling operation isrepresented by a d  d nonsingular integer matrix M [3]. Signalswith finite energy belong to a Hilbert space l2(Zd) with the innerproduct is defined ashx; yi =Xn2Zdx[n]y[n];and thus the l2-norm is kxk =phx; xi =qPn2Zdx[n]2.With the notation as shown in Fig. 1, the filtering and down-sampling operation for the LP yieldsy0[n] =Xk2Zdx[k]h[Mn   k] = hx[];~h[  Mn]i; (1)where we denote ~h[n] = h[ n].Similarly, the upsampling and filtering operation results inx0[n] =Xk2Zdy0[k]g[n Mk]: (2)Writing signals as column vectors, for example x = (x[n] :n 2 Zd)T, we can express those operations as left matrix multi-plicationsy0= Hx and x0= Gy0;where H and G correspond to (# M)H and G(" M), respec-tively. In particular, H hasn~h[n Mk]on2Zdas its rows andG has fg[n Mk]gn2Zdas its columns. In the sequel, we denoteI as the identity matrices with appropriate sizes depending on thecontext.With this matrix notation, the output of the LP in the highpasschannels isy1= x  y0= x  GHx = (I  GH)x:By combining the previous relations, we can write the analysisoperator of the LP asy0y1| {z }y=HI  GH| {z }Tax: (3)The usual inverse transform of the LP compute: ^x = Gy0+y1, or^x = G I| {z }Tsy0y1: (4)It is easy to check that TsTa= I for any H and G. This agreeswith the well-known fact that the LP can be perfectly reconstructedwith any pair of filters H and G.2.2. Frame AnalysisThe frame analysis is a powerful concept in analyzing redundantrepresentations [4]. For a detailed introduction to frames, readersare referred to [5] (Chapter 2). It can be shown that a linear op-erator is a frame if and only if it is invertible on its range with abounded inverse. For the Laplacian pyramid, there always exists abounded reconstruction algorithm (i.e. the usual algorithm), hencethe LP is guaranteed to be a frame. In this case, the frame operatoris represented by the left matrix multiplication with Ta.Since the LP is a redundant transform, its frame operator ad-mits an infinite number of left inverses. Among those, the mostimportant one is called the dual frame operator, which is repre-sented by the pseudo inverse of Ta[6]Tya= (TTaTa) 1TTa: (5)When there is additive “noise” in the frame coefficients, thepseudo inverse eliminates the influence of errors that are orthog-onal to the range of the frame operator. Therefore, if instead ofhaving access to y = Tax we have ^y = y + e, then the pseudoinverse provides the solution ^x = Tya^y that minimizes the residualkTa^x ^yk. This is called the least-squares solution.Let T 1abe an arbitrary left inverse of Ta. With the noisemodel setup as above, we can write the reconstruction error whenusing T 1aas^x  x = T 1a^y   x = T 1a(y + e)  x= T 1a(Tax+ e)  x = T 1ae: (6)3. A TIGHT FRAME CASEFirst, let us consider the special case where filters in the LP areorthogonal filters with respect to the sampling matrix M. Moreprecisely, the orthogonal condition requires thathg[]; g[  Mn]i = [n]; and (7)h[n] = g[ n]; or H = GT: (8)Such filters can be designed using well-known methods, in-cluding separable and non-separable ones [3, 7]. The LP withorthogonal filters possess an important property as stated in thefollowing proposition.Proposition 1 The Laplacian pyramid with orthogonal filters, i.e.satisfying (7) and (8), is a tight frame.Proof: Under the orthogonality condition, (2) becomesx0[n] =Xk2Zdhx[]; g[  Mk]ig[n Mk]:Therefore if we denote V as the space spanned by the set oforthogonal vectors fg[   Mk]gk2Zdthen x0is the orthogonalprojection of x onto V . Together with the fact that y1= x  x0,using the Pythagorean theorem leads tokxk2= kx0k2+ ky1k2= ky0k2+ ky1k2: (9)where the equality kx0k = ky0k comes from the fact that y0are the coefficients in the orthonormal expansion (2) of x0in V . 2Since Tais a tight frame with TTaTa= I, its pseudo inverseis simply the transposed matrix. Thus we haveTya=GTI  GGTT= G I  GGT: (10)So the optimal (least-square) reconstruction using the pseudoinverse is^x = Tyay = Gy0+ (I  GGT)y1= G(y0  Hy1) + y1:The last expression above is derived in order to reduce thecomputational complexity of the pseudo inverse Tya. It leads to anefficient filter bank structure for reconstruction of the LP that isshown in Fig. 2.It is worth to stress two important facts here. First, the usualinverse is different from the pseudo inverse, and thus it is subopti-mal. Secondly, the pseudo inverse in this case has a symmetricalstructure with the forward transform hence it has the same orderof complexity.++_M MH Gx^y0y1Fig. 2. The proposed reconstruction scheme for the Laplacianpyramid in Fig. 1. It is the pseudo inverse when the filters areorthogonal.4. GENERAL CASESIn this section we consider the more general cases where H and Gare arbitrary filters. Even though any frame operator has a pseudoinverse, for complexity reason we would like to only consider thecases where the inverse operator has a fast structured transform.Inspired by the tight frame case, we restrict our attention to thereconstruction algorithm that has a structure as shown in Fig. 2.We then turn the problem around by asking for which filters suchan algorithm is indeed an inverse or pseudo inverse. This has thesame flavor as the filter design problem for perfect reconstructionfilter banks. The following proposition provides an answer to ourquestion.Proposition 2 ([8]) 1. The reconstruction shown in Figure 2is an inverse transform of the LP if and only if two filtersH and G are biorthogonal with respect to the sampling lat-tice M, which means the upper branch of the LP (GH) is aprojector, or HG = I.2. Furthermore, that reconstruction is the pseudo inverse ifand only if the upper branch of the LP (GH) is an orthogo-nal projector.Recall that given a Hilbert space H , a linear operator P map-ping H onto itself is called a projector if P2 = P . Furthermore,if P is self-adjoint or P = PT then P is called an orthogonalprojector.Remark: It is interesting to note that the two conditions forthe LP in the above proposition, i.e. projection and orthogonal pro-jection, are exactly the same with the conditions for the improvedLP’s that was studied in [9]. Those conditions lead to LP with in-terpolation and least squares LP, respectively. The motivation forthose modifications in [9] is to minimize the prediction error y1of the LP, whereas our motivation in this work is to have a betterreconstruction algorithm for the LP.So the minimum requirement for the filter bank shown inFig. 2 to be a bona fide inverse of the LP is the biorthogonalitycondition of the filters H and G. In this case, the LP has an inter-esting geometrical interpretation. Let us define two subspaces Vand ~V that are spanned by fg[ Mk]gk2Zdand f~h[ Mk]gk2Zd,respectively. These are also the column and row spaces of G andH. For all x in l2(Zd), the result of the approximation channel ofthe LP: x0= GHx, is a projection of x onto V . Furthermore,since HG = I the approximation error y1= x  x0has propertythatHy1= H(x  GHx) = Hx  Hx = 0:In other words, the approximation error is perpendicular to ~V .This fact is illustrated in Fig. 3. The upper channel in the LP (GH)can be called an oblique projector [10] and denoted by ~PV.V~VWxx0=~PVxy1=~PWxFig. 3. Graphical representation of a oblique projector. x0=~PVxis the projection of x onto V such that the approximation errorx ~PVx is perpendicular to ~V .Let us define W the orthogonal complementary subspace ofV , then it is easy to verify that y1= x ~PVx is a projection of xonto W such that the error is parallel to V (refer to Fig. 3). Denotethis projection as ~PW,~PW= I  GH.Now let us compare the reconstruction using the usual method,denoted REC-1:x1= Tsy = Gy0+ y1; (11)and the one using the filter bank in Fig. 2, denoted REC-2:x2= Try = Gy0+ (I  GH)y1: (12)Those two reconstruction algorithms are different in the wayof handling the prediction error y1. More specifically, the REC-1 method adds directly y1while the REC-2 method adds the ~PWprojection of y1to the reconstructed signal. So when there is noisein the LP coefficients, the REC-2 method eliminates the influenceof errors in ^y1that are parallel to V .In order to have more quantitative measurements in the per-formance of the two methods, suppose we wish to approximate xgiven ^y = Tax+ e. With no further information about the errorin the LP coefficients e, it makes sense to chose ^x to minimize theresidual kTa^x ^yk. As mentioned before, the optimal linear so-lution to this problem is the pseudo inverse of Ta. Using this as themeasurement for the performance in reconstruction, the followingresult states that REC-2 always performs better than REC-1.Proposition 3 ([8]) Assume that H and G are biorthogonal fil-ters. Let ^x1and ^x2be the results of reconstruction from noisy LPcoefficients ^y using REC-1 and REC-2, respectively. Then we havekTa^x1 ^yk  kTa^x2 ^yk ; (13)where equality holds if and only if ^y1= 0.5. MULTILEVEL LAPLACIAN PYRAMIDSWe now consider the situation where the LP scheme is iteratedon the coarse version. We concentrate, for the sake of simplicity,on the one-dimensional case where M = 2. By separating thesequences y1into two subsequences y10and y11of even and oddnumbered samples, respectively, we can re-formulate the LP as athree channels filter banks (see Fig. 4).Then, using polyphase-domain analysis, we can write the newanalysis filters as [8]K0(z) = 1 H(z)G0(z2);K1(z) = z  H(z)G1(z2); (14)2 22222+H GK0K1F0F1xx^y0y10y11Fig. 4. Laplacian pyramid as an oversampled filter bank.where Gi(z) are the polyphase components of G(z), orG(z) = G0(z2) + z 1G1(z2):Since H(z) and G(z) are both lowpass filters, K0(z) andK1(z) are high pass filter.On the synthesis side, for the REC-1 reconstruction we simplyhave_F0(z) = 1; and _F1(z) = z: (15)Whereas, for the REC-2 reconstruction, the synthesis filtersareF0(z) = K0(z); and F1(z) = K1(z): (16)When the filter bank in Fig. 4 is iterated on the lowpass chan-nel y0, then by using the multirate identity which says that G(z)followed by upsampling by 2 is equivalent to upsampling by 2followed by G(z2) [7], we have equivalent synthesis filters at then-level of a multilevel LP asF(n)i(z) = Fi(z2n 1)n 2Yk=0G(z2k); i = 0; 1: (17)preceded by upsampling by 2n. The impulse responses of F(n)i(z)are the synthesis functions for the LP.Now, consider what happen when (15) and (16) are substi-tuted into (17). In the REC-1 method, we see that the synthesisfunctions for the LP are all low frequency signals. Thus, the errorsfrom highpass subbands of a multilevel LP do not remain in thesesubbands but appear as broadband noise in the reconstructed sig-nal. In [11], this effect was noted as the most serious disadvantagesof the LP for coding applications.On the other hand, in the REC-2 method, the synthesis func-tions have the similar frequency characteristics as the analysisfunctions. Therefore, reconstruction using REC-2 method reme-dies the previous mentioned problem of the LP.6. EXPERIMENTAL RESULTSIn order to test the actual performance of the new reconstructionmethod with “real” signals we set up the experiments as follows.An input signal is transformed with a Laplacian pyramid. TheLP coefficients are uniformly quantized, and this introduces errorsdepending to the quantization step. The quantized LP coefficientsare then used to reconstruct the original signal using the usual andthe new inverses. This is the typical scenario in lossy compression.We use the signal-to-noise ratios (SNR) between the original andthe reconstructed signals as the measures of quality.Experiments with different test signals and filters (usingboth orthogonal and biorthogonal ones) indicate that the newreconstruction provides significant improvements over the usualmethod. Fig. 5 plots a typical result. In this case the input signal isthe 512 512 “Lena” image and the LP is computed with 5 levelsusing the Daubechies’ orthogonal maximally flat 8-taps filters [5].As can be seen, gains of more than 1 dB are actually obtained ateach quantization step.0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 0.05222426283032343638Quantization stepSNR (dB)With the usual inverseWith the new inverse  Fig. 5. Comparison of reconstruction from quantized Laplacianpyramid coefficients using two different inverses.Acknowledgments: The authors would like to thank Prof.Michael Unser, Dr. Vivek Goyal and Irena Maravic for their help-ful comments and suggestions.7. REFERENCES[1] P. J. Burt and E. H. Adelson, “The Laplacian pyramid as acompact image code,” IEEE Trans. Commun., vol. 31, no. 4,pp. 532–540, April 1983.[2] M. N. Do and M. Vetterli, “Pyramidal directional filter banksand curvelets,” in Proc. IEEE Int. Conf. on Image Proc.,2001, submitted.[3] P. P. Vaidyanathan, Multirate Systems and Filter Banks,Prentice-Hall, Englewood Cliffs, NJ, 1993.[4] V. K. Goyal and J. Kovacˇevic´, “Quantized frame expansionswith erasures,” Applied and Computational Harmonic Anal-ysis, 2001, to appear.[5] I. Daubechies, Ten Lectures on Wavelets, SIAM, Philadel-phia, PA, 1992.[6] R. A. Horn and C. R. Johnson, Matrix Analysis, CambridgeUniversity Press, 1985.[7] M. Vetterli and J. Kovacˇevic´, Wavelets and Subband Coding,Prentice-Hall, Englewood Cliffs, NJ, 1995.[8] M. N. Do and M. Vetterli, “Framing pyramids,” IEEE Trans.Signal Proc., 2001, submitted.[9] M. Unser, “An improved least squares Laplacian pyramidfor image compression,” Signal Proc., vol. 27, pp. 187–203,May 1992.[10] M. Unser and A. Aldroubi, “A general sampling theory fornon-ideal acquisition devices,” IEEE Trans. Signal Proc., pp.2915–2925, 1994.[11] E. P. Simoncelli and E. H. Adelson, “Subband transforms,”in Subband Image Coding, J. W. Woods, Ed., pp. 143–192.Kluwer Academic Publishers, Inc., Boston, MA, 1991.

Frame reconstruction of the Laplacian pyramid

Abstract

Similar works

Full text

Available Versions

Infoscience - École polytechnique fédérale de Lausanne