EEE TRANSACTIONS ON IMAGE PROCESSING. VOL I. NO. 2. APRIL 1992 Image Coding using Wavelet Transform Marc Antonini, Michel Barlaud, Member, IEEE, Pierre Mathieu, and Ingrid Daubechies, Member, IEEE together with its implementation as described by Mallat [27], satisfies each of these conditions. wvelet s; this new method involves two steps. First, we use aa wavelet transform and a vector quantization coding der to obtain a set of biorthogonal sub- scheme. The wavelet coefficients are coded considering a asses of images; the original image is decomposed at different noise shaping bit allocation procedure. This technique ex maintains constant the number of pixels required to describe in the image data, enabling bit rate reduction. %N s an Section lI describes the wavelet transforms used in this ory, the wavelet coefficients are vector quantized using a multi- paper. After a quick review of wavelets in general, we resolution codebook. Furthermore, to encode the wavelet coef- explain in more detail the properties and construction of hich assumes that details at high resolution are less visible to re regular biorthogonal wavelet bases. We then extend tI m ize a picture as quickiy as possible at minimum cost we scheme with separable filters. The new coding scher wavelet t ransforesi is tranimission scheme It is shown that the next presented in Section ll. We focus particularly in this transmission ients, on the asymptotic coding gain that can be achieved orthogonal wavelet, multiscale py. using vector quantization in the subimages, and on the algorithm, vector quantization, noise shaping, pro- optimal allocation across the subimages. Experimental re- sults are given in Section IV for images taken within and outside of the training Se ACTIc different fields, digitized images are replacing A. A Short Review of Wavelet Analysis onal analog images as photograph or x-rays Wavelets are functions generated from one single func costly. The information contained in the images must 1/ therefore, be compressed by extracting only the visible elements, which are then encoded, The quantity of data (For this introduction we assume t is a one-dimen- nvolved is thus reduced substantially sional variable). The mother wavelet v has to satisfy A fundamental goal of data compression is to reduce dx v(x)=0, which implies at least the bit rate for transmission or storage while maintaining (Technically speaking, the condition on y should an acceptable fidelity or image quality. Compression can dw |Y(w)l2-< oo, where y is the Fourier trans be achieved by transforming the data, projecting it on a form of y; if v(n) decays faster than t for t-o, then basis of functions, and then encoding this transform. Be- this condition is equivalent to the one above). The defi cause of the nature of the image signal and the mecha- nition of wavelets as dilates of one function means that nisms of human vision, the transform used must accept high frequency wavelets correspond to a I or narrow nonstationarity and be well localized in both the space and width, while low frequency wavelets have a>I or wider frequency domains. To avoid redundancy, which hinders width compression, the transform must be at least biorthogonal The basic idea of the wavelet transform is to represent and lastly, in order to save CPU time, the corresponding any arbitrary function f as a superposition of wavelets algorithm must be fast. The two-dimensional wavelet Any such superposition decomposes f into different scale transform defined by Meyer and Lemarie [31], [24], [25], levels, where each level is then further decomposed with a resolution adapted to the level. One way to achieve such ipt received February 7, 1990: revised March 26 a decomposition writes f as an integral over a and b of chies is with at&T Bell Laboratories Murray Hill. N) 07974. tice, one prefers to write f as a discrete superpose ition(sum IEEE Log Number 9106073 rather than integral). Therefore, one introduces a discre 1057714992s3.00@1992IEEE
IEEE TRANSACTIONS ON IMAGE PROCESSING, vOL L. NO 2. APRIL 199 tization,a=0,b= nboao, with m, nEZ, and ao >1, for decomposition. Such filters are well known since the bo >0 fixed. The wavelet decomposition is then work of Smith and Barnwell [35] and of Vetterli [37] .The ∫=∑cmn(f)ψnn (1) xtra ingredient in the orthonormal tion is that it writes the signal to be decomposed as a st compositions of this type were studied in [14], [15]. For blocks. The filters must satisfy the additional condition a0=2, bo= 1 there exist very special choices of y such that the ym. constitute an orthonormal basis, so that II H( Cm,n(f)=(Vm. n f )= dx ym. a(x)f(x) decay faster than C(1+5I in this case. Different bases of this nature were con- structed by Stromberg [36], Meyer [31], Lemarie [24 H(E)=2-1/2∑hne Battle [7, and Daubechies [16]. All these examples cor- respond to a multiresolution analysis, a mathematical tool This extra regularity requirement is usually not satisfied nvented by Mallat [27], which is particularly well adapted by the exact reconstruction filters in the ASSP literature to the use of wavelet bases in image analysis, and which gives rise to a fast computation algorithm B. Applications of Wavelet Bases to Image Analysis n a multiresolution analysis, one really has two func 1) Biorthogonal Wavelet Bases: Since images are tions: the mother wavelet y and a scaling function g. One mostly smooth(except for occasional edges)it seems ap- lso introduces dilated and translated versions of the scal ing function, m, n (x)=2-m/2p(2-x-n). For fixed m propriate that an exact reconstruction subband coding scheme for image analysis should correspond to an or the m n are orthonormal. We denote by m the space thonormal basis with a reasonably smooth mother wave spanned by the m, n; these spaces Vm describe successive let. In order to have fast computation, the filters should each with resolutic For each m,theψ be short(short filters lead to less smoothness, however a space wm which is exactly the orthogonal complement so they cannot be too short), On the other hand it is de- in Vm-1 of Vm; the coefficients (vm nf>, therefore. de- sirable that the FIR filters used be linear phase, since such scribe the information lost when going from an ap- without the need for phase compensation. Unfortunately ation of f with resolution 2"I to the coarser ap- there are no nontrivial orthonormal linear phase FIR fil nation with resolution 2m. All this is translated ters with the exact ruction property [351, regard into the following algorithm for the computation of the less of any regularity considerations. The only symmetric cmn(f)=〈ψm,n,f〉( for more details,se27]): xact reconstruction filters are those corresponding to the cm.(f)=282n-K@m-L&() Haar basis with all other h,, gn=0 ann(f)=∑h2Aam-1.(f) One can preserve linear phase(corresponding to sym- (2)metry for the wavelet)by relaxing the orthonormality re- quirement, and using biorthogonal bases. It is then still here g,=(-1)h-I+I and h=2/2dxo(x-n)o(2x). possible to construct examples where the mother wavelets In fact the am, n(f)are coefficients characterizing the pro- have arbitrarily high regularity ction of f onto Vm. If the functionf is given in sampled In such a scheme, we still decompose as in(2),but form, then one can take these samples for the highest or. reconstruction Ines der resolution approximation coefficients ao. n, and (2)de scribes a subband coding algorithm on these sampled val an-1.()=∑[2n-1an(f)+gm=1Cm.(f)(4) ues, with low-pass filter h and high-pass filter g. Because where the filters h, g may be different from h, g. In order of their association with orthonormal wave\nf oses, these to have exact reconstruction, we impose filters give exact reconstruction, i.e gn=(-1)"h an-1.(f)=2[hn-1am(f)+82n-1Cmn( 8n=(-1)hn+ Most of the orthonormal wavelet bases have infinitely g differently from supported y, corresponding to filters h and g with infi- the usual exact reconstruction ding scheme nitely many taps. The construction in [16] gives v with with synthesis filters different fr ecomposition fil finite support, and therefore, corresponds to FIR filters. ters. If the filters satisfy the additional condition that It follows that the orthonormal bases in [16] correspond to a subband coding scheme with exact reconstruction property, using the same FiR filters for reconstruction as II A(2-5)and II H(2-5)
ANTONINI et al.: IMAGE CODING USING WAVELET TRANSFORM decay faster than C(1+|)-sas|→∞, for some ∈>0, where H(o) B(5)=212∑hne)H()=2-12∑hne ) then we can give the following interpretation to(2)and (4). Define functions c and o by d(x)=∑hnd(2x-n)ando(x)=∑hnd(2x-n) Their Fourier transforms are exactly the infinite products Feauveau explores the construction from the point of view Fable functions, compactly supported if the filters h and ically one has two hierarchies of spaces in the bior- thogonal case, each corresponding to one pair of filters (x)=∑gn(2x-n)andv(x)=∑gn(2x-n) It is shown in [12] that arbitrarily high regularity can be achieved by both y and y. provided one chooses suf- Then, the a,mn(/) and cm. n() in(2)can be rewritten ficiently long filters. In particular, if the functions and Y are, respectively. (k-1)and(k- 1)times continu ously differentiable, then the trigonometric polynomials am. ()=(om. m/)=2 m/ dx om. n(a)f() (1 +e"f, respectively. so that the length of the corre m(f)=(vm, n,f>=2-m/ dx vom n(x)f(r) By(5). divisibility of A(E)by(1+e-5)means that y will have k consecutive moments zero and reconstruction is simply f=∑(ψmn,f)n,n xrv(x)=0, for I=0,I For more details concerning this discussion, see [12 The filter bank structure with the associating wavelets It is well known (and it can easily be checked by using nd scaling functions is depicted on the following sub band coding scheme(Fig. 1) Taylor expansions)that if y has k moments zero, then the If the infinite products in(6a)decay even faster than Coefficients (Vm. n, f)will represent functions f. which imposed above, then o and o and consequently y and y re k times differentiable, with a high compression poten will be reasonably smooth. Note that (7)is very similar tial(many coefficients will be negligibly small) o the orthonormal decomposition described in Section Many examples of biorthogonal wavelet bases with rea II-A; the only difference is that the expansion of f with sonably regular v and y can be constructed; for our ap- respect to the basis yo. n uses coefficients computed via ym n. which is linked to the number of zero moments of the dual basis tmn with y different from y. This interpre, 4, is more important than the regularity of the yo,. or the oding schemes; in particular, convergence of the infinite number of zero moments of y. Within the limits imposed tation is not possible for all exact reconstruction subband by the support widths, we will, therefore, try to choose as larg ∑hn=21/2and∑hn=21/2 、气 n terms of trigonometric polynomials H()and() he exact reconstruction requirement condition on h and Moreover, (7) can only hold h given in(5)reduces to(for symmetric filters) ∑(-1)yhn=0and∑(-1)yhn=0 H()B()+HE+)HE+丌)=1 Most exact reconstruction subband coding schemes do Together with divisibility of H and H, respectively. by ot satisfy these conditions (1+e-i 'and(1+e/5). this leads to(see [12]) Biorthogonal bases of wavelets hat Herley and Vetterli [38]. Reference [12] contains a de sin(/2)}+sin(/2)yR)(9) onditions stated above. the wavelets do indeed constitute merically stable bases(Riesz bases)and a discussion of where R(E) is an odd polynomial in cos(E), and where necessary and sufficient conditions for regularity. In [18] 2/=k+k(symmetry of h and h forces k+ k to be even)
EEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 1. No. 2. APRIL 1992 TABLE I FILTER COEFFICIENTS FOR THE SPLINE FILTERS W n 19/64 1/8 3/643/128 1/2 1/4 Many examples are possible. We have studied in par cular the following three examples, which belong to three different families 2)Spline Filters: One can choose, e. g,R=0,with A(E)=cos(E/2 e kE/2 where x=0 if k is even,K=I if k is odd. This corresponds to the filters called"spline filters"in[12](because the corresponding function is a B-spline function)or"binomial filters"'in[38](because the h are simply fficients). It then follows hat H(5)= ∑ We have looked at one example from this family; it responds to /=3, k=2. The coefficients h,and are listed in Table I; the corresponding scaling functions -1s and wavelets are plotted in Fig. 2 It is clear that the two filters in the first example have very uneven length. This is typical for all the examples in this family of"spline filters 3)A Spline variant with Less Dissimilar Lengths:This family still uses R =0 in(9), but factorizes the right hand side of (9), breaking up the polynomial of degree sin(E/2)with real coefficients, one to be allocated to H the other to H, so as to make the lengths of h and h as close as possible The example presented here is the"smallest',one in 4 this family(shortest h and h ); it corresponds to I-4 and k= 4. The filter coefficients are listed in Table II: the corresponding scaling functions and wavelets are plotted Note that, unlike examples I and 3 where the 2-1/h 2-1/h, are rational the entries in Table lI are truncated decimal expansions of irrational numbers. The functions in examples I and 2 look very similar(compare Figs 2(a)and 3(a)); a more detailed analysis shows that the one in example 2 is more regular, however. Both correspond to 4 vanishing moments for y 4)Filters Close to Orthonormal Filters: Finally, there Xist many e for which R*0. In particular ther very close to each other, and both very close to an or. Fiter 3 w ith ling function. 2 and wavelet un cti for exam ple in pline conormal wavelet filter Surprisingly, for the first example of this series, one of the two filters is a Laplacian pyramid filter pro- wavelets are plotted in Fig. 4. It is clear that the scaling rery similar R(E)=48 cos(5)/175. The filter coefficients are listed similar y and 4. Note that in this case, the filter coef in Table Ill; the corresponding scaling functions and cients are again rationa
ANTONINI et aL.: IMAGE CODING USING WAVELET TRANSFORM TABLE I FILTER COEFFICIENTS FOR THE SPLINE VARIANT WITH LESS DISSIMILAR ± 2-1/h,0.6029490.266864-0.078223-0.0168640.02674 2-120.5575430.295636-0.028 0.045636 d etsψ. ngths: /=4-kk=4).(a) Scaling function o.(b)Scaling function o(c)Wavelet 4. (d) Wavelet 2 various extensions of the one-dimensional wavelet FILTER COEFFICIENTS FOR EXAMPLE 3. THE ENTRIES ARE RATIONAL, AND transform to higher dimensions. We follow Mallat [27 LAPLACIAN PYRAMID FILTER PROPOSED IN [9]. IN THIS CASE and use a two-dimensional wavelet transform in which =2=k,E= horizontal and vertical orientations are considered pre 0 ±4 esential In two-dimensional wavelet analysis one introduces like o in the one-dimensional case, a scaling function o(x. y) 17/28 such that o(x, v)= o(x)o( v) (11) The two biorthogonal filters in this example are both close to an orthonormal wavelet filter of length 6 con- where o(r) is a one-dimensional scaling function structed in [17], where it was called a:coiflet. Being Let v(x)be the one-dimensional wavelet associated with an orthonormal wavelet filter, the coiflet is nonsymme- the scaling function (x). Then, the three two-dimer tric. The filters in this example are shorter than in exam- sional wavelets are defined as ples I and 2, but k is also smaller. The next example in this family corresponds to k= 4(and /=4); the filters h ψ"(x,y)=o(x)ψ(y) and h then have length 9 and 15: they are both close to a coiflet of length 12 (x,y)=ψ(x)q(y) 5) Extension to the Two-Dimensional Case: There ex ψ(x,y)=ψ(x)(y) (12)