UFDC Home  Search all Groups  UF Institutional Repository  UF Institutional Repository  UF Theses & Dissertations   Help 
Material Information
Record Information

Full Text 
A MULTISCALE SPLINE DERIVATIVEBASED TRANSFORM FOR IMAGE FUSION AND ENHANCE NIlINT By IZTOK KOREN A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFII.I.l NT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 1996 Copyright 1996 by Iztok Koren :Iojimn starsem AC(IK NOWLED(GEM 1ENTS I would like to thank my mentors Dr. Fred J. Taylor and Dr. Andrew F. Laine for their hell) and support. They gave me the freedom to explore research directions that interested me most an(t provided mine with the excellent work environment. I am grateful to Dr. Jose ('. Principe. Dr. John M. .I. A\nderson. and Dr. Kermit Sigmon for serving on my committee and for their comments on my disser tat ion. I would also like to thank all the members of the High Speed I)igit al Architec ture Laboratory and the Image Processing Research Group for their friendship and Shelp. Last but not least. I want to thank my family for t heir support and under standing. and for bearing with me when I was unbearable. TABLE OF CONTENTS ACKNOWLEDGEMENTS .......... LIST OF TABLES ............... LIST OF FIGI'RES .............. ABSTRAC(T ................... CHAPTERS 1 INTRODUCTION ........... 1.1 Shortcomings of Traditional \l,.w 1.2 Thesis Overview ......... 1.3 Notation ............. 1,,, of Wavelet Analysis 2 SIGNAL PROCESSING USING CENTRAL BSPLINES 2.1 Central BSplines: Definition and Properties . . 2.2 BSpline Signal Interpolations . . . . . . 2.3 BSpline Signal Approximations . . . . . . 3 TRANSFORM IN ONE DIMENSION . . . . . . 3.1 1D Discrete Dyadic Wavelet Transform Revisited . 3.2 Rem arks . . . . . . . . . . . . 4 TRANSFORMS IN TWO DIMENSIONS . . . . . 4.1 2D Discrete Dyadic Wavelet Transform Revisited . 4.2 Steerable Functions . . . . . . . . . 4.3 Steerable Dyadic Wavelet Transform . . . . . 4.4 Multiscale Spline DerivativeBased Transform . . 5 IMPLEMENTATION ISSUES . . . . . . . . 5.1 Finite Impulse Response Filters . . . . . . 5.2 Infinite Impulse Response Filters . . . . . . 6 APPLICATIONS . . . . . . . . . . . 6.1 Im age Fusion . . . . . . . . . . . 6.2 Image Enhancement . . . . . . . . . . . ... .. 61 7 CONCLUSION . . . . . . . . . . . . . ... .. 66 REFERENCES . . . . . . . . . . . . . . . . ... .. 68 BIOGRAPHICAL SKETCH . . . . . . . . . . . . ... .. 72 LIST OF TABLES 2.1 Transfer functions of direct Bspline filters for orders from 0 to 9. . 12 3.1 Impulse responses h (n), g(n), l(n). and k'(n) for p=0 and d E {1.2.3}. 2.5 3.2 Impulse responses h(n). y(it), l(,n), and k(n) for p= 1 and d E {1.2. 3}. 25 3.3 Impulse responses h( (n). g(n). 1(n). and Ak(n) for p=2 and d E { 12.3}. 2. 6 4.1 Imnpulse responses t (n) for e {0, 1,2} . . . . . . .... :32 LIST O FFIGITRES 1.1 Discrete wavelet transform of two signals translated to each other.. 2 1.2 Discrete wavelet transform and "algorithLine a trous." . . . . 3 2.1 Spline functions .,(x.r) and /,(;x) for p C {0. 1.2.3. } . . . . 10 2.2 Fourier transforms p(,) an(d [,(&') for /) C {0.1.2.3. 4} . . . 13 o 0 2.3 Splines 3,,(,x) and J,(,Jc) for p {0.1,2. 2,3,4) . . . . . . ... .. 16 3.1 Filter bank implementation of a 1D discrete dyadic wavelet transform. 21 3:1.2 Comparison of two discrete implementations using sinc(.x) as an input. 23 3.3 Wavelets .(.r)= ,) for p e {0.1,2} and d Ce {1 2.3}. . . . 27 4.1 Filter bank implementation of a 2D discrete dyadic wavelet tranllsform. 35 4.2 Filter bank implementation of a spline derivativebased transform... 48 6.1 (Cinomparison of image fusion results. . . . . . . . . . 62 6.2 (C'omparison of image enhancement results . . . . . . ..... .. 64 6.3 The magnitude of filter coefficients for '( (..). . . . ... .. 65 Abstract of Dissertation Presented to the Graduate School of the V'niversitv of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy A MI'LTIS(CALE SPLINE DERIVATIVEBASED TRANSFORM FOR I\1 GE FUSION AND ENIHANCE\II'NT By lztok Koren December 11'v, Chairman: Dr. Fred J. Taylor Cochairman: I)r. Andrew F. Laine l.,j,,r Department: Electrical and (Comniputer Engineering Analyzing images along distinct scales is advantageous in many image pro cessing and computational vision tasks. Traditional wavelet transforms have become a popular tool for multiscale image analysis and image compression suffer from arti facts, lack shift and rotation invariance, and may introduce aliasing and anisotropies. VWe propose a highly redundant representation tliat eliminates these artifacts. A new transform is constructed as an extension of the discrete (ldyadic wavelet transform with a wavelet being the first derivative of a central Bspline. We begin the construction in one dimension by providing an efficient initialization procedure for the comnputation of the discrete dyadic wavelet transform and by extending the transform to encompass higher order derivatives. This modified onedimensional discrete dyadic wavelet transform will serve as a basis for deriving a steerable twodimensional transform. Such a transform is advantageous in that it does not introduce any of the artifacts described above. How ever, exact reconstruction is problematic via implementation in the spatial domain. \We solve this problem by devising a splinebased approximation. The resulting algo rithm for computing a multiscale spline derivativebased transform is implemented efficiently as a filter bank consisting of separable filters alone. For both finite and infinite impulse response filters present in the filter bank. we specify implementa tion details of a realization that alleviates the boundary effects caused by the use of circular convolution. Finally, we use the derived transformin for multiscale image fusion and enhance ment. We show empirically that the transform does not introduce artifacts commonly reported for traditional waveletbased implementations and provides more flexibility for different fusion and enhancement strategies by enabling orientation analysis. CHAPTER 1 INTRODI( TION Analyzing images across multiple scales and resolution has become a power ful tool for solving compelling problems in computational vision, image processing, and( pattern recognition. Wavelet theory encompasses multiscale and multiresolution representations, such as subband filtering [40], image pyramids [4]. and scale space filtering [48]. into a unified mathematical framework. In the area of image processing, there remain few research areas to which wavelet analysis has not been applied. For example, problems in image compression, denoising. restoration, enhancement, reg ist rat ion. fusion, segmentation, and analysis, have all been approached wit li d(listinct kinds of wavelet processing. Though ubiquitous, wavelel analysis is not without problems of its own. Lack of translation invariance, one of the major problems of the wavelet transform [10]. is in multiple dimensions accompanied with lack of rotation invariance. 1.1 Sh ic I. uiii .,I I If 1 1 i J i ,i.il NM .f1 1k l Id n' \\.I% e A wavelet transform in its most commonly used orthogonal or biorthogonal forms is not translation and rotationinvariant. By translationinvariant transform. we mean a transform that commutes with a translation operator. Since we will deal primarily with discrete transforms in this work, we constrain the translation parameter to integer multiples of a sampling period. Lack of translation invariance of the discrete wavelet transform is illustrated in Figure 1.1. Here. we can clearly see how a translation of the input signal by one (a) (b) Figure 1.1: (a) Original signal and (b) signal translated one sample to the left with its discrete wavelet transform coefficients shown across dyadic scales 2". m EC {1.2,.3.4}. saminple results in a completely different set of trainsform coefficientls orthogonall wavelet DA1 B41 [10] was used in this experiment). Noninvariance under translations of an orthogonal and biorthogonal wavelet transform is (due to lower sampling density at coarser scales.2 A straightforward way of dealing with this problem is to consl ruct a redundant transform by using the same sampling frequency for the input signal and all scales of the transform. A filter bank implementation of such a Itransform, called "algorithme a trous" [15]. is based upon the fact that downsampling followed by filtering is equivalent to filtering with the upsamnpled filter before the downsampling, as shown in Figure 1.2. G(z) 2 G(z) G(z) 2 G(Z) H(z) 2t G(z) 2v H(z) G(z4) H(z) 2, H(z) H(z) 2Y H(z4) (a) (1)) Figure 1.2: Filter bank implementation for (a) a discrete wavelet transform and (b) "'algorithme a trous" decompositions for three levels of analysis. Lack of rotation invariance is another short coming of traditional (i.e.. orthog onal and biorthogonal) wavelet techniques. In defining rotation invariance, we are a bit less strict than with translation invariance. We do not require that the transform commutes with a rotation operator here. Even in the case of a simple filtering, this would limit us to circularly symmetric filters only. Our requirement for analysis is a 'The number in DAUB4 refers to twice the order of the wavelet (i.e., two in this case). 21n practice, since analysis is p)erformned over a finite range of scales, a discrete wavelet transform is translationinvariant by translations determined by the coarsest scale (e.g.. sixteen samples for tlie analysis from Figure 1.1) [10]. transform that enables rotationinvariant processing. As an example of such a trans form, let us consider filtering with the first derivative of a twodimensional Gaussian probability density function in two directions, specifically, along x and along yaxis. By linearly combining the results of filtering in these two directions, filtering with the first derivative of a Gaussian in any direction can be computed. This fact was used by Canny [6] for edge detection. A determined edge direction rotates as an input image is rotated. After choosing the fundamental properties of the transform, one must decide upon the basis functions to be applied. For our studies, we selected basis functions that well approximated derivatives of a Gaussian. because (1) the Gaussian proba bility density function is optimally concentrated in both time and frequency domain, and thus suitable for timefrequency analysis, (2) higher order derivatives of a Gaus sian can be. similar to the first derivative, used for rotationinvariant processing [12], and (3) the Gaussian function generates a causal (in a sense that a coarse scale de pends exclusively on the previous finer scale) scale space [2]. The last property makes possible scalespace "tracking" of emergent features. 1.2 Thesis Overview We wish to construct an invertible, translation and rotationinvariant two dimensional transform that decomposes an image using a set of basis functions which are approximations to derivatives of the Gaussian probability density function. To obtain a translationinvariant representation, we will use "algorithme a trous," or more precisely, we will extend the discrete dyadic wavelet transform [28], which uses "algorithme a trous." Derivatives of a Gaussian will be approximated by wavelets which are derivatives of central Bspline functions. In Chapter 2. we review the basics of Bspline processing necessary to derive transforms in later chapters. Next. Chapter 3 examines a dyadic wavelet transform in one dimension. We present extensions of the discrete dyadic wavelet transform to higher order derivatives and to even order spline functions. To compute a discrete transform from a contin uous one, the discrete computation must be properly initialized. We devise a new initialization procedure that is shown to be more accurate than the one i,_._,I ed by Mallat and Zhong [28]. After the derivation of the transform, we point out relevant connections to scalespace filtering and reconstruction from edges. In case of the first derivative of a Gaussian, a rotationinvariant transform can be obtained by a tensor product extension of the onedimensional transform to two dimensions. For second order derivatives, a tensor product approach can approximate Laplacian of a Gaussian, but cannot perform orientation analysis (due to its isotropic nature). In Chapter 4, we first examine the above two cases, and then develop a new multiscale spline derivativebased transform, which in addition enables rotationinvariant processing for derivatives of orders higher than two. Given that the developed transform is highly redundant, an efficient imple mentation is extremely advantageous. This issue is addressed in Chapter 5. When filtering of a signal is performed by circular convolution, boundary effects may result. A standard way of dealing with this problem is by mirrorextending the input signal to the filter. We combine such an extension with symmetry/antisymmetry of the filters, to achieve savings in both computation time and memory. Finally, in Chapter 6. we present two sample applications: image fusion and enhancement. WVe demonstrate that the new transform does not cause artifacts com monly associated with traditional wavelet techniques and is an attractive analytic tool for designing novel fusion and enhancement strategies. 1.3 Notation We use symbols N, Z, and R for the sets of naturals, integers, and reals, respectively. L2(R) and L2(R2) denote the Hilbert spaces of measurable, squareintegrable functions f(x) and f(x y/). respectively. The inner product of two functions f(x) E L2(R) and g(x) E L2(R) is given by (f(,r),g(.r)) = f(xr) g(x) d.r. The norm of a function f(x) E L2(R) is defined as 1f1 = f(rx) I dr. The convolution of functions f(xr) E L2(R) and g(x) E L2(R) is computed as f *g(xr) = f(t)g(r t)dt. and the convolution of two functions f(x, y) G L2(R2) and g(x. y) C L2(R2) equals .f g(X. ) f f(t.,. t,,)g(x t,. y ty) (dtx dt. The Fourier transform of a function f(x) E L2(R) is defined as /()= f(x)j' dx. and the Fourier transform of a function f(x,y) C L2(R2) is equal to f (xr .y)J() d+Y dy. 12(Z) and 12(Z2) stand for the spaces of squaresummable discrete signals f(n) and f(t, ny). respectively. The ztransform of a discrete signal f(n) E 12(Z) is defined as F(z)= > f(n)z" The convolution of discrete signals f(n) C 12(Z) and g(n) E 12(Z) is equal to f *g(n) = (i f(m)g(n m), and the convolution of discrete signals f(ii., ?Y) E 12(Z2) and g(nx. ny) E 12(Z2) is given by f* y(,r. 1y) P > .f(rx, my)g(nx II,, .y y). 'm, 7n = 'X' The Fourier transform of a discrete signal f(n) G I2(Z) is equal to the z transform evaluated on the unit circle F(^)'= .f',)J . and the Fourier transform of a= discrete signal f (, ) 1(Z2) is defined as aunl the Fourier transform of a discretee signal f( us., us) G /2( Z2) is (deflned as f(n~. n~) ( flr X For later use, we define the following functions: 1. the unit impulse function (.r+) := { ,(.r) := {1 2. the unit step function for x = 0 otherwise. for x for x' 3. the rectangular function rect(x) := for .1 for xI 4. the since function sin( rx) sinc(,r) := and 7T r 5. the unit impulse sequence I1 for n = 0 ( [) := 0 otherwise. where .r E R and n E Z. (CHAPTER 2 SIGNAL PROCESSING USING CENTRAL BSPLINES In this chapter we briefly review fundamentals of spline processing needed for derivations in (Chapters 3 and 4. 2.1 (C'ert ,,I 1 t l ,i,.. fl fi ,i_ ii '11 rid P i, 11i1  (liven real numbers c grid point) sequence Xr.. 2.....X,. if it is (1) a polynomial of degree p or less in each interval [xr.,rx+l]. i = 0.1,... m, and (2) continuous in its derivatives up to the order p I on the interval x[.o. +] (i.e., CP [xo, r'm+l]). Here. we will concentrate primarily on basis splines (Bsplines). or more pre cisely,. cent ral Bsplines having knolts at i E Z for p odd and at i + 1 for p even [34]. central l Bsplines of order ) (with p+ I knots) are defined as I,.r :+= P+1 + ) p+1 1+ P+ P.i=0 I Figure 2.1 shows 3(.r)and their Fourier transforms jj(,) for p E {0. 1.2.3. 41}. A family of functions {,p(r m)} ,,,Z forms a basis of S'". a space of order p spline functions with knots at i for p odd and at i + 1 for p even ( E Z) [31. 35]. 2 Except for p = 0. the basis functions {/.J3p(xr inm)} are not orthogonal. Let us list some properties of functions 3p(x) [3. 34]: 1. 3p(.r) are nonnegative functions with a support of length p+ 1. p+1 times 2. 3 3(x) = ,30(x), where "*" denotes the convolution operator, or, equivalently, in the Fourier domain: 3 (~Sin( )P2 \ 2 3. .3p(x) +((p+l ,) 3 (x + + X) pl(X 1 2 JpX) Xp i(x ( + ) + (P+ ( p 4. J = P1 (,"' + ) _P1 ( X 1) . Another interesting property of Bsplines is the fact that they converge to a Gaussian as their order tends to infinity. Unser t al. [44] derived the Gaussian approximation virp. 22 71 .p where ,p = Ratio .. for example, is already well below i 2'' Denoting by SP a spline function space spanned by { 3p(x inm) for p odd and by {3p(.r j i)} for p even subscriptt in SP refers to the fact thliat spline functions have knots at integers), spline spaces form a nested sequence ... C S), C ... C SP C S' C S'_i C ... C S'_, C .... By orthogonalizing this basis functions a multiresolution analysis of L2(R), from which the BattleLemari. wavelet bases stem, can be built [10]. Here, we will not pursue constructions of orthogonal, semi orthogonal. or biorthogonal splinewavelets any furtherreader looking for a detailed treatment of this subject may find a good starting point in books [8. 9]. '2.' [BSplin_,, Sit ,i;dl Thil,.r)l,,l iti'ii Unser 0 al. developed a fast digital filtering scheme for Bspline signal pro cessing [43]. They defined a discrete Bspline of order p) and expansion factor (spacing between knots) mi as P n \ bplm(n) := 1, 1. n ? G Z. Figure 2.1: Spline functions (a) 3p(x) and (b) q/(,r(x) for p E {0.1.2, 3.4}. 3 2 2 3 '* Henceforth, if the distance between knots equals one. we will write bp(n) instead of bpi(n). Interpolation of a discrete signal s(n) G 12(Z) by sp(x') G SP' using central Bsplines s(n)= .s(x) = c(i)3(i) (2.1) n i=^x=n can now be written as a convolution sum s(n)= c* bp(n). (2.2) If s(n) are samples of a function s(x) bandlimited to [7. 7r] (i.e., the support of its Fourier transform S'(u) is in [7r, w]). it can be shown that .p(.r) + s(x) as p oc [1, :',,]. In [43]. they refer to a linear operator by which Bspline coefficients c(n) can be obtained from samples s(n) as a "direct Bspline transform." Equation (2.2) therefore represents "indirect Bspline transform" of a sequence {c(n)}. After taking the ztransform of (2.2). the direct Bspline filters are found to be p1;1(n) = 7'{[Bp(z)]1}. Since Bspline functions ip(.r) are compactly supported, indirect Bspline filters bp(n) are finite impulse response (FIR) filters, while direct Bspline filters b1(n) are infinite impulse response (IIR) filters. Aldroubi et al. [1] showed that IIR filters b1(n) are stable (i.e., the region of convergence of B;1(z) includes the unit circle [30]) for any order p. Note that both indirect and direct B spline filters are symmetric, which follows from the fact that central Bsplines 3p(x) are symmetric. Table 2.1 shows the ztransforms of direct Bspline filters for the first ten orders. We postpone the discussion on implementation details of Bspline filters until Chapter 5. Table 2.1: Transfer functions of direct Bspline filters for orders from 0 to 9. Z B# (z) 0 1 1 1 *) 8 _+6+ 1 3 6 :+4+:I' 4 384 2+76 +230 . 't  5 120 z2+26z+66+26z1 +2 6 46080 + .''"; +" 1D l I:'.4* '+ *** t  7 5040 Z3+ .,, l i "I Il,.+ I I ' + I . T  8 +.: 10321920 : I I .' "'' + ''! " )III"."  I "'l I' ."' It :  l ." *.  q +362880 _________________________ 36280__+L_______"____' _ ,,*.,,* t 1 l,,, t  l',.1 lll s .'1+ I I,.11 + It ., ,.' +. Instead of using Bspline interpolation as given by (2.1) it is sometimes con venient to express the interpolating function Sp(x) in terms of discrete sampl)les s(n) )= E .()b(* ). (2.3) where ip(.r) is the cardinal spline of order p. In the frequency domain, cardinal splines converge to an ideal lowpass filter with cutoff frequency 7r (i.e.. qp(x) * sinc(x)) as p tends to infinity [1. 45]. which establishes the asymptotic equivalence with Shannon's sampling theorem [37]. Using (2.1) and (2.2) with (2.3) cardinal splines can be related to Bsplines: tip(X)= b '(i) (, )(2.4) Cardinal splines i;p(,r) and 4p(w) for p E {0, 1, 2, 3, 4} are shown in Figure 2.2. 2 68 6 Figure 2.2: Fourier {0,1,2,3.,4}. transforms of spline functions (a) p({.,) and (b) rjp(,) for p E 4 6 8 11 6 6 89 8 6" 2 4 6, S 6 B, 6 4 I 6, p=4 8. 6 4, 4r 6, ., 6.  2.3 B i liin, qiujii, \d  l i i ,:..i i ,;fit I I  Central Bsplines are also simple to use when the goal is function approxima tion. Leastsquares Bspline approximation of s(x) E L2(R) is achieved by computing the orthogonal projection of this function onto SP. We have (x)= = () = d(i)0.i(x i) (2.5) with 0 d(O) = (s(x),3,(x i)), (2.6) where W3(x) bp+= 1 + 3pi(x 0 is the dual spline of order p [45]. Spline functions ,3(x.) and 3p(x) satisfy the biorthog onalitv condition 0 3p(xr M), 3p(x n) = 6(m n), ?, In e Z. 0 (Note that. since both 3p(x) and 3p(x) form a basis of SP. they can be interchanged in (2.5) and (2.6).) 0 0 Figure 2.3 shows functions 3,(x) and their Fourier transforms 3p(w) for p C {0.1. 2. 3.4}. An interesting alternative to the minimum L2norm (i.e., leastsquares) ap proximation of a signal is obtained by computing oblique instead of orthogonal pro jection of the signal onto the spline function space. Unser and Aldroubi proposed an independent specification of the sampling and approximation spaces [42]: a linear op erator maps coefficients of the input signal expansion over sampling space basis into the coefficients of the approximation space basis expansion from which the signal's projection onto the approximation space is recovered. Constraining the entire system to be linear, shiftinvariant for integer translations, and consistent (i.e., the system acts as an identity operator for functions that belong to the approximation space), the obtained solution for the signal approximation is the projection of a signal onto the approximation space perpendicular to the sampling space [42]. Analogously to (2.5) and (2.6) this projection can be expressed as OCO s,(x)= ,(i) 3(x i) (2.7) with a(i) = 1 (s(,), 3,(x i)). where q1 is the convolution inverse of the crosscorrelation sequence q12(i) = (3,(r i),03(x)) = b,+,+, (,). When the sampling space S' and the approximation space S' are identical (i.e.. s = r), an orthogonal projection given by (2.5) and (2.6) results. (Note that the described oblique projection is not restricted to spline function spacesthe only requirement is that both sampling space basis and approximation space basis are Riesz bases of the corresponding function spaces [42].) Signal approximation (2.7) is particularly attractive in situations where the sampling space is given a prior (e.g., by the impulse response of the acquisition device [42]) or when such an approximation is close to the optimal leastsquares solution but simpler to implement than the orthogonal projection (e.g., [47]). 66 3> 22 5 I $" 62 6 6 (a) (14 aa 00 Figure 2.3: Spline functions (a1) 3p(2) and (b) their respective Fourier transforms ,Jp( ) for {0.1,.2,3,4}. C(IAPTER 3 TRANSFORM IN ONE DIMENSION In this chapter, we will augment the onedimensional discrete dyadic wavelet transform [28] to obtain a strong foundation for derivations of twodimensional trans forms in Chapter 1. 3. 1 I) D i i,.t. D 1'.,,, Id W x ,.l.i Ti, i',,r i ii R.'. I I o A discrete wavelet transform is obtained from a continuous representation by discretizing dilation and translation parameters such that the resulting set of wavelets constitutes a frame. The dilation parameter is typically discretized bv an exponential sampling with a fixed dilation step and the translation parameter )v integer multiples of a fixed step [10]. Unfortunately, the resulting transform is variant under translations, a property which makes it less attractive for image analysis (cf. Section 1.1). As we have already mentioned in (C'hapter 1. sampling the translation param eter with the same sampling period as the input function to the transform results in a translationinvariant, but redundant representation. The dyadic wavelet transform proposed by Mallat and Zhong [28] is one such representation. Let us begin with a brief review of properties of the dyadic wavelet transform as described in [28]. but included here for completeness. The dyadic wavelet transform of a function s(.x) E L2(R) is defined as a sequence of functions L{ 'z() (3.1) where 1Vms(,r) = s t'v,(.r). and (x() = 2"('(2 ') is a wavelet '(.r) expanded by a dilation parameter (or scale) 2"'. Note the use of convolution instead of an inner product. To ensure coverage of the frequency axis the requirement on the Fourier trans form of t',,(Jx) is the existence of A1 > 0 and B1 < oc such that A, < j ,(2'1c)2 K B_ 17 = N: is satisfied almost everywhere. The constraint on the Fourier transform of the (nonunique) reconstructing function \ (x) is N: ^(2 %) <(2"+ ) =1. m7 =  ,x. A function s(x) can then be completely reconstructed from its dyadic wavelet trans form using the identity .s(xr) = l 'm.' \m(.). m __  ix. where \,,(x) = 2" (2"'x). In numerical applications, processing is performed on discrete rather than continuous functions. When the function to be transformed is in the discrete form, the scale 2'" can no longer vary over all m G Z. Finite sampling rate prohibits the scale from being arbitrarily small, while computational resources restrict the use of an arbitrarily large scale. Let the finest scale be normalized to 1 and the coarsest scale set to 21. The smoothing of a function E(x) L2(R) is defined as ='(.) Z S o (x). where O,,(xr) = 2'o(2 ,x) with 0m E Z. and 0(x) is a smoothing function (i.e., its integral is equal to 1 an(d 0(x)  0 as xI + oc). In [28]. a real smoothing function )(x) was selected, whose Fourier transform satisfied I ( ')2 t!(2' ) (2 ). (3.2) m=l In addition, it was shown that any discrete function of finite energy s(n) E 12(Z) can be written as the uniform sampling of some function smoothed at scale 1. i.e., s(n) = So0f(n). where f(.r) C L2(R) is not unique. Thus. the discrete dyadic wavelet transform of s(n) for any coarse scale 2" was defined as a sequence of discrete func tions ( + + {IW(u s)},n1.n}z. where .s is a u(.r) dependent sampling shift. The above initialization s(n) = Sof(n) is rather standard in the discrete wavelet transform computation [10], although it yields correct results (i.e., the dis crete wavelet transform is equal to the samples of its continuous counterpart) only when s(n) = Sos(n). Here, we will concentrate on wavelets which are derivatives of spline functions and this will lead us to a simple initialization procedure [46] that alleviates the above problem. For a certain choice of wavelets the discrete dyadic wavelet transform can be implemented within a fast hierarchical digital filtering scheme. Next. we shall summarize the relations between filters, wavelets, and smoothing functions. First, let us introduce a real smoothing function ,(x) such that (3.2) can be rewritten as1 V7(^ f (2 ... v(2`,), (3.3) rn =0 and let us set o(.r) = 1,(x) (i.e., we restrict ourselves to wavelets which are spline functions). 'Note that the sum index determines the range of scales of the discrete transform: using (3.2) we have t'(2w) and \(2.') at the finest scale of the transform, while for (3.3) we get () and \(w). Computing (3.3) for the finest two scales shows that ,'(.) \(w.)= (. ) (.) .,,(2w) ?(2 ). (3.4) /,(2,,;) can be related to () by expressing ip(2,) as (cf. Proposition 1 of [46]) 1 sinl(2 2sin and using E .(1+= Ci ( 2) +6) = Si,,( ) ): J) (2L,) = cos ()) 3() (3.5) (Note that a relation similar to (3.5) can be derived for integer scales provided that the dilation parameter and the order p) are not both even '111.) Let F(,.') be a digital filter frequency response and let F(') = /*'*F(). If we choose ( = ( ('_ ) (w), (3.6) ( L) = (,) ( ). (3.7) and If(,), ; (cos ())+ (3.9) where ., E {0. } is a filter dependent sampling shift needed for g(n). l(,). k(n). and h(n) to be FIR filters, and insert Equations (3.5) (3.9) into (3.41). we observe the relation between frequency responses of the filters (;()K(,) + H(w)L() = 1. (3.10) Similar to orthogonal and biorthogonal discrete wavelet transforms, the dlis crete dyadic wavelet transform can be implemented within a hierarchical filtering scheme. To derive such a digital filtering scheme, let us assume that .(*.') from (3.1) is banidliinited to [Tr. ]. I'sing Shannon's sampling theorem [37] and (3.6) in the definition of the dyadic wavelet transform (3.1) with m =0 shows .O(.r) .(O)sinc(ti) Y g,(1()3p(x  ) (11. By making use of the fact that the cardinal spline functions tend to the since function as their order r approaches infinity, and employing (2.1) we can write 11js^) ^ W3 ~ S()B\^,{) 3p(&)G, ^). or. by using (3.5) and (3.9). 2F{ r),I } Br B.+'(I )^ G (2`,,;) ILJ I (2'&). (3.11) 1=0 Equation (3.11) entirely specifies the discrete dyadic wavelet transform de composition. while the reconstruction follows from (3. 1)(3.9). Three levels of a filter bank implementation are shown in Figure 3.1. (Note that the initialization is the same as the one proposed in [46] and that except for the prefiltering and postfilter ing. this scheme is implementing "algorithme a trous.") Noninteger shifts at scale I for filters with ,= are rounded to the nearest integer. 2 3.:'.,i Bp^,q(() G^(2(, K,(2..i) Bp ]{(O)) B,((0) H (o) G(,(4(i) K,(4,) L,((,) H (m2) L,(2o) HF(4() L,(4o) Figure 3.1: Filter bank implementation of a onedimensional discrete dyadic wavelet transform decomposition (left) and reconstruction (right) for three levels of analysis. \\e will now l)erform a siml)le experiment which will illustrate the difference between the implementation of the discrete dyadic wavelet transform as originally proposed in [28] (i.e., without prefiltering and postfiltering) and the one from Figure 3.1. Let .s(x) = sinc(.r). p = 2. and g(n) = 2o( + l)2e(i) (this particular choice for p) and 9(n) results in the same wavelet as was used by l.,ll.i and collaborators [27. 28]). The dyadic wavelet transform of .s(.r) at a scale 2'" (3.1) in tl he frequency domain is then l ( ;) = (_,(2"'% ) 2(2',) rect( ) (3.12) The Fourier Itransform of the discrete dyadic wavelet, transform of .s(?;) = i(n) at a scale 2'" using spline based initialization follows from (3.11) B{ I(,,s(u)} = 3 ( )B,+ ( ) ( 2... (2 .. ?) 1 H ,(2'.;). (3.13) i=0 while the one using the algorithm from [28] equals ,(n) = G "(2) Im Hs(2'i). (3.14) i=0 In Figure 3.2 a comparison of magnitudes of (3.13) and (3.14) versus (3.12) is shown: in Figure 3.2(a) magnitudes of (3.12) (solid) and (3.14) (dashed) are plotted for in E {0.1, 2.3}. while the dashed curves in 3.2(b) represent magnitudes of (3.13) with r=5. By choosing the appropriate order r. (3.13) can approximate (3.12) in tlie interval [7. 7] arbitrarily good. while originally proposed (3.14) has troubles at finer scales. Mallat and Zhong [28] noticed that there was a problem with their discrete transform computation, and introduced a set of constants associated with the discrete transform coefficients at dyadic scales. They chose the values of constants such that the transform coefficient modulus maxima remained constant over all (Idyadic scales for a step edge input signal. In relation to Figure 3.2(a) this is equivalent to multiplying F{ lIT',,;(n)} by a distinct constant for each m. Clearly, this can improve over the situation depicted by Figure 3.2(a). but can still not compete with the described spline based initialization. 0 5 Figure 3.2: (a) Fourier transform magnitudes of the dvadic wavelet transform of s(,r) = sinc(,r) (solid) and the originally proposed discrete dyadic wavelet transformn [28] of .s(n) = 6(n) (dashed). (b) Fourier transform magnitudes of the dyadic wavelet transform of ,.(r) (solid) and the discrete dyadic wavelet transform using quintic splines for interpolation of ,( n ) (dashed). n0 1 5, 0 25 Next. we will choose filters in the filter bank implementation of the discrete dyadic wavelet transform. As already mentioned, we are interested in wavelets which are derivatives of spline functions. G(uw) in (3.6) is therefore the Fourier transform of the difference operator centered around s (cf. Property 4 in Section 2): (;(,) = (2j sin ())d. (3.15) where d is the order of the derivative and the sampling shift for this filter is = dmod2 2 Since H(u,) was already given by (3.9). the remaining two filters to be deter mined are L(&c) and K(,;). Both of them are (as is true for )(x) and \(x)) nonunique. If we choose L(,.') such that we can express K(.') in terms of a finite geometric series having the smallest number of elements for an arbitrary p. we get Ld. J ) 1 m 1 1d+l  / 2\ (p+l)(2"n 1) L(,) = 3s (ir+1 2 L+J Ccos (p)) (3.16) m =1 / and 1 7 (:.N\drod 2 (CP / /P,\\ 2 1\2 ()) in (S co \)2 ) ( (:3.17) where [x] denotes the largest integer smaller than x, the sampling shift for L(L}) is the same as the one for H(.') (i.e.. =s mo 2 ). and the sampling shift for K(&,) is the same as the one for G(,c). Note that Equations (3.9) and (3.15)(3.17) work fine from the mathematical point of view, but in practice the reconstruction may become cumbersome when both p) and d are large (the lengths of impulse responses h(n), g(n), l(n). and k(n) are p+2, d+1., (p+l)(d(d+l) mod2)+l, and p(/+(p+l)(dmod2)+1. respectively, while for the frequency responses of the decomposition filters we observe that limp_,, Hs({) = 6,,(, + 2?o7) and limd j_(2) G,() = 6,(' + (2n+l)7r) with n C Z). It is also worth noting that K(w.') is a lowpass filter when p is even (i.e., the reconstruction function \(x) is a wavelet only for p odd). Tables 3.1. 3.2, and 3.3 list impulse responses of the four filters for p {0. 1. 2} and d E {1. 2. 3}. while Figure 3.3 shows wavelets .'(r) = +() for the same values of p and d. Wavelets from this family have a support of length d+p+ 1. regularity order p, and are either symmetric or antisvmmetric. ible 3.1: Impulse responses h(n), g(,'). 1(n), and k(n) for p=0 and d E {1.2,: Sd= 1 d = 2 d = 3 n h(n) g(n) l(n) k(n) g(n) 1(n) k(n) g(n) 1(n) k(n) 2 1 1 0.5 1 1 3 0.125 0 0.5 1 0.5 0.25 2 0.5 0.25 3 0.625 0.0625 1 0.5 0.25 1 0.5 1 0.625 0.0625 2 ___ _____________ 0.125_____ 3}. Table 3.2: Impulse responses h(in). g(in), l(n). and k(n) for p I and d E {1.2.,3}. As already discussed, wavelets with p= 2 and d= 1 from a family of wavelets with p even and d= 1 were used in [27, 28], whereas filters with p= 1 and d= 2 from a family of filters with p odd and d= 2 were employed by Laine and collaborators [11, 22, 23]. Here described transform puts no restrictions on the choice of p or d whatsoever, and uses a better initialization procedure than the one originally proposed in [28]. h(n) d = 1 d = 2 , l(n) g(n) 141) g(n) k(n) 1 0.25 1 0.0625 1 0.0625 0 0.5 1 0.3125 2 0.375 1 0.25 0.3125 1 0.01t21 2 O i ii.2,', TL nT h(ne) g(n) !(7?) k(nn foi p(.n E 1. i 2 0.125 0.015625 0.01.".,I 1 0.375 1 0.125 0.10'1i:7 1 0.125 0.125 0 0.375 1 0.375 0.34375 2 0.375 0.46875 1 0.125 0.375 0.34375 1 0.375 0.125 2 0.125 0.109375 0.125 0.015625 3 ___________0.015625_____ d=3~ n l1(n) g(n) 1(n) k(n) 4 0.001953125 0.000244140625 3 0.017578125 0.003662109375 2 0.125 1 0.0703125 (i.i2t. l.7,t875 1 0.375 3 0 N.I..37 0.I'I 2'1.13 25' 0 0.375 3 0.7. .i' 1,''5 0.13113710'.(1 7 1 0.125 1 0.'. .', 11 " 0.13037109375 2 0.0 ','i371 I r I'll 1,_'3125 3 0.0703125 0.I i2'i. I,71875 4 0.017578125 (I.Ii.1231 " 1 7" 5 0.001953125 0.0002441 ti.5 Table 3.3: Impulse responses b(??), g(1), 40!), and k(n) for p=2 and d C {1,2,3}. 2 3 ~j .... ... 1 3 2 1 0 1 2 3 (c) (x) d2 Op+2(X) Figure 3.3: (a) Wavelets (a) ;'(x) (= ' ). (b) t'(x) = d2+ and (c) ,(xr) '3+(J) for p=O (dashdotted), p=l (solid), and p=2 (dashed). Page Missing or Unavailable 29 the zerocrossings of I,n.(.r) for d = 2, a scalespace image of the representation sim ilar to [48] can be obtained. The only differences are that curves in the scalespace plane are computed for dyadic scales only and that 0,,, is a close approximation of a Gaussian instead of a true Gaussian. CHAPTER 4 TRANSFOR\IS IN TWO DIMENSIONS In this chapter, we will extend the onedimensional transform from Chapter 3 to two dimensions1 and derive several twodimensional transforms with translation and at least approximately rotationinvariant decomposition. _1 D D1,1.l. D'.,h W ,.,L, I Tiawf.ii R," iit,.1 The dyadic wavelet transform of a function s(x, .) E L2(R2) is defined as a set of functions [28] {t ( 2y), ( )} (4.1) where li(. y) = x '(.r, y,) for i = {1,2} and ,(x. y) = 22,2 t,(2 mx. 2y) are wavelets t y(x.y) expanded by a dilation parameter 2". To ensure coverage of the frequency space there must exist an A2 > 0 and B2 < x such that 2 A,2 < Z(2"';,2" )12 < B2 m=o; i=l is satisfied almost everywhere. If (nonunique) functions XI(x. y). \2(x. y) are chosen such that their Fourier transforms satisfy x, 2 V ^ \ ^ I'lVO"1' *~)m. ^ \< 9 m"*, ,' 2^ 2 f ' y\ 1, yii mn=,x i=1 the function s({ry) may be reconstructed from its dyadic wavelet transform by x 2 m= xX i=1 where X, (x, y) = 2 X (2.r, 2y). 1For extensions to higher dimensions, please refer to [18, 19]. However, when processing discrete functions the scale 2'm may no longer vary over all m E Z. Let thle finest scale be normalized to 1 and the coarsest scale set to be 2M. Let us. similar to [28], introduce a real smoothing function O(x. y.) such that its Fourier transform satisfies x< 2  _/ , \' :?' *" ., *" ^i(*gm .)m 11\) 0 Yr i) 12 =_ ', I: z^ M2}\ 'r^l (41.2) m=O i=1 Here, as in one dimension, a finite energy discrete function (s(n,,. n V) E 12(Z2)) can be written as the uniform sampling of some function smoothed at scale 1: .s(n, n ) f= n.(11. ?j.). where f(xr,y) E L2 (R2) is not unique, and 5', f(.r. y) f O, (.r. y). This led M.:ill,, and Zhong [28] to a twodimensional analog of the onedimensional definition of the discrete dyadic wavelet transform:2 { _./'!(,,l+. .. I + ). { Il p,1 ,, + ny +s), 11 f(n,+2 n +s + ),,[o. ][,,. ,Y, . We will use, as in Section 3.1. a splinebased initialization procedure. To implement a multidimensional discrete dyadic wavelet transform within a fast hierarchical digital filtering scheme, the wavelets were chosen to be separable products of onedimensional functions: ,(.r, ) = l .) 0 (y), (4.3) '2(X, y) = ( (y) 0(V), (4.4) where 6(x) and (.r) were chosen as described in Section 3.1 (i.e.. o(x) = 3p(x) and /'(w) specified by (3.6)). From (4.3), (4.4), and (3.6). we may write s i3 e io Gs(eutth f) ines(sleo he tran)sfor (4.5), "'2 in GS,ecio ).1 utef a (4o6) 2 'As ill Section 3. 1, we put the finest scale of the t~ransformi at in = 0. where G(,() is given by (3.15) for (I E {1,2}. Choosing = (4.7) where K(c) and 7Ti(c) are digital filter frequency responses, we may compute (4.2) for the finest two scales by 2 2' ,. 2 ,) <(2 ,, 2 y) = I(w .. )2 o(22,). 2, )2. (4.9) i=1 Inserting the terms defined by (4.5). (4.6), (4.7), (4.8), (3.5). and (3.9) with o(rC', Y) = 3p(wr) 3,9(wu) into (4.9) results in K(Y.)G(t:)T7(w,) + K(^,)G(^y)Ti(,)) + jH(.\)2H(',)j2 = 1. (4.10) Equation (4.10) represents a relation between the frequency responses of the digital filters used to implement a multidimensional discrete dyadic wavelet transform and is a multidimensional analog to (3.10). Solving (4.10) for Tl(,c) by substituting K(,)G(,,) from (3.10) yields the closed formula [28] T ) = (1 + IH(')12) (1.11) In Table 4.1 we provide the filter coefficients for T1(w') from (4.11) for p E {0.1.2). All other filters from (4.10) were already specified in Section 3.1. Table 4.1: Impulse responses t1(o) for p E {0, 1.24. n p=0 p=l p=2 3 nI 1i'7 [25 2 0.03125 0.046875 1 0.125 0.125 0.1171875 0 0.75 0.6875 0 1,,1,,_ 1 0.125 0.125 0.1171875 2 0.03125 0.046875 3 0.0078125 As in the onedimensional case, a twodimensional discrete dyadic wavelet transform can be implemented as a fast hierarchical filtering scheme. To derive such an implementation, we, similar to the onedimensional case from Section 3.1, use the definition of the twodimensional dyadic wavelet transform (4.1) and require (cov) = 0 for 1olj > rt or Iy  > 7t. Using Shannon's sampling theorem in two dimensions [16]. (4.5), (4.6). and m =0, we get 0/g(s rx. '(x. y)'' W ls(x~y) = / Z^ ,s(i .,iy)sinc(tr. i )sinc(ty y). x .. iz ,=,>5 i=,x s(,,,),^(.r ^ n'py ty) dt dt and W0's(.r.y)= J > s(i.xi)sinc(tx i)sinc(t iy). X .. J X *x  '  "3 ( ' 1' ) E 9 ("y( v (Xd'' t *X *p' t.r) >3 gs(I')Jp(y tfy m)d1(11/. As in one dimension, we make use of the fact that the cardinal spline functions converge to the since function as their order r tends to infinity, and write \Vqs(C2. .^y ) C"^q S(COr, ~y )B, '(^.)B',71(C~y)L~p+r+1 (ol*) p++, +( 'y) (s (Cr) and s (C. O ) ,4 ,5'(CO.r, 'y )r1 (C) 71 (C) p+rl (C4x)!)p+r+l (Cy) G(sy) or for in E [0, 1) and discrete signal processing .{Wf (.r" }14s(ox., W) B? 1 (,x) Br 1 (cy)Bl+r+l ()" r 1 .BP+r+Iy) ( Gs(2mwx) 11H i,(2.,cx)I_ (2',). (4.12) i=0 and F{f 11 ,;'.(.r,.9) :, ,y= } (o;, oy) B, (B7 )(o 1)B ,7;(y B)B,+,.+I(c,, )" nm1 *Bp+,.+,(Wy) G_,(2mwy) I H_,(2,w)H_,(2'1y). (4.13) i=0 Equations (1.12) and (4.13) describe thlie decomposition part of the filter bank implementation of a twodimensional discrete dyvadic wavelet transform. The recon struction part can be obtained from (1.7)(4.9) with (3.5) and (3.9). I he entire filter bank implementation of the transform is shown in Figure 4.1. Except for the pre filtering and post filtering, we readily recognize the implemental ion proposed in [28]. Using the fact thai a wavelet i'(xr) is equal to a first (d= 1) or a second (d/=2) derivative of a spline function 1,,+d(,r), (4.3) and (4.4) may be rewritten as , 0. 0) d',(,1 21) ,,'(..O ., i.d e {1,2}. w herie dI(.. Y/) = 3p+l(.) 3p(/) and l)2(r,.9) = 4,(.*) 3,,+,ff,') Let us devote I1,s(,.y) (1Vb.(..). llj(.y)). V ( ). A V2 ( + ). and assume tha t )'(,r.y.) can be approximated by d)(.r.y) for both i E {1.2}. For d(= I il then follows that lI',,, .x, y) = 2"'V(,s )(x.y). (1.14) Thus for d(=2 we can write 2 Sll, .s(x y) = 22" A(, i), )(r. y). (4.15) t=1 With )(.r. y) being a Gaussian, finding local extrema of (4.14) in the direction of gradient V corresponds to the filtering stage of a Canny edge detector [6], and finding zerocrossings of (4.15) correspoiinds to the filtering carried out with a Marr Hildreth edge detector (Laplacian of Gaussian) [29]. (Note that both edge detectors GS(WO Qs(2(o) Gs,(20y) Hs(2N) J(2o ) (a) Ks(o,) T,(o~y) Br(Ox) Br,((Oy)  H_(2(o ) H_(2(y) (b) Figure 4.1: Filter bank implementation of a twodimensional discrete dyadic wavelet transform (a) decomposition and (b) reconstruction for two levels of analysis. H',(.) denotes the complex conjugate of HL(s.c). 0 HBVO) ff)O Bp+r+l((O,) Bp+r+l((Oy) Hs(X) H,((Oy) E+r+#)') RP'+,+]((OY) Tj(ox) K(oy) K,(2O)x) T,(2(oy) T,(2wx) Ks(2oy) R*,(O)X) ft*s((Oy) involve postprocessing). Edge detection based on finding local extrema of 1',,s(r. q) or zerocrossings of Z, 12 I. ,s(.r. y) is therefore an approximation to the (C'annv or the NI rrHildreth edge detector over a range of dyadic scales. The differences stem from the fact that d(.r. y) is neither a Gaussian nor is U'(.. 1) equal 1to (.x, y). 4.2 Steerable Functions When extending the transform from (Chapter 3 to two dimensions, one of the first questions that come to mind is how to deal with the fact that derivatives can be defined in any direction of the plane.3 In case of a first derivative of a Gaussian, one can simply compute first derivatives of a Gaussian in x and y directions and then determine the derivative in any direction from these two direct ional derivatives [6]. For higher order derivatives of a Gaussian. Freeman and Adelson [12] showed that order+1 directional operators are needed for synthesizing the operator at any orientation. In fact. functions with the property of expressing their arbitrary rotations as linear combinations of some functions are not limited to derivatives of a Gaussian. Below, we briefly restate some of the results from [12]. A twodimniensional function is called steerablele" if its rotations generate a finite dimensional space. Rotations of a steerable function f(r, 0) can therefore be expressed as 1 f(r.0 0o) = ci(0).f(r,, 0). (4.16) i=1 where 00 denotes an arbitrary angle. {ci(Oo)} stands for a set of interpolating func tions. {fi(r. 0)} is a set of basis functions, and "r = 2 + y2 and 0 = arg(.r. y) are polar radius and angle, respectively. If f(r.0) represents a filter kernel, the result of filtering with a rotated fil ter /(r.0 00o) can be computed simply by {ci(Oo)} weighted linear combination of outputs from basis filters {/,(r, 0)}. Only the outputs from basis filters need to be 3Second derivatives of central Bsplines can be used, as we saw in Section 41. to approximate Laplacian of a Gaussian for approximately rotationinvariant, although not directional, processing. precomputed and then the output from a filter rotated by any angle 0()u can be found by interpolating between them. When a large number of rotations of a template filter is required, the above scheme can lead to substantial savings in both computational cost (time) and memory requirements (space). The necessary condition for a function ff(r. 0) to be steerable is that /(1r. 0) is bandlimited in its polar angle: N f(r. 0)= a a (r)d*" (.1.17) n=.N This can be verified by inserting (1.17) into (4.16) and by assuming, for con venience. that f.(i,, 0) = f( r. 0 0,): I (i,(1.) ,, oo = i(0o)a,,(,/') (7 i' ( .18) i=1 where {Oi} is a set of user defined angles and n C [N'. 0]. 4 Since only nonzero coefficients a,,(r) are of interest, both sides of (4.18) can be divided by av(r). This yields a matrix equation from which interpolating functions ci(Oo) can be determined: 1 1 1 1 l c(00) C0A1 ( 10.. 1 f'9O C2(00) S. : : (4.19) J O jNOI C. N 0^ ^ 2 ... f^ '_ .cI(^o) For coefficients a,, = 0 the rows corresponding to each n were removed from the matrix formulation shown in (4.19). For this system to have a solution, the angles {0,} must be chosen such that the columns of the matrix are linearly independent. In [12] they proved that the minimum munmber of basis functions fi( r. 0) needed to steer f(r. 0) according to (4.16) is equal to the number of nonzero coefficients a,/(r) in the Fourier series expansion (4.17). To conclude this brief description of steerability, let us only remark that func tions which are not steerable (i.e.. do not have a finite number of terms in (1.17)) 4Note that the constraints are the same for n E [.V,1] and n E [1, ]. so that a subset of all possible values for n E [.V, V] can be taken. can be approximated with steerable functions (a singular value decomposition was employed for approximating such functions efficiently [31]). and that the technique of expressing transformed versions of a function as linear combinations of a fixed set of basis functions is not limited to rotations (extensions to translations [39]. scalings [31. 39]. and general transformations [14] have been reported). 4.3 %l.',i. 1,l. P '..,1i, W ,, ,.1, I Ti..ii,fi II Returning to the question from the beginning of Section 4.2. the answer seems obvious: one needs to construct a steerable analog to the onedimensional transform from Chapter 3. Steerable transforms are nothing newquite a few [13. 17. 1,. 39] have been devised, some of them [17, .'j] exactly for the computation of directional derivatives. Here. we are not interested in any directional derivatives: we want to use derivatives of central Bsplines which, as the order of Bsplines increases, tend to derivatives of a Gaussian. W\e define a steerable dyadic wavelet transform of a function S(xa. !y) G L2(R2) at a scale 2'. in E Z. as [21] \ (x, 0) = (x, y). (4.20) where ,,'(.r. y) denotes ',, (.r, ) rotated by 0, ,,,(x.y.) = 2m 2'.(2"`x. 2'mj), i,'(a. y) is a steerable wavelet that can be steered with I basis functions, and 0i =  with i e {1.2 .... I}. Analogously to the onedimensional case (cf. Section 3.1) we require the two dimensional Fourier plane to be covered by the dyadic dilations of '(2'. 2",,): there must exist A3 > 0 and B3 < cc such that x, I A3 < 2 1)1(,0 2 < B3 (4.21) m=,x, =1 is satisfied almost everywhere. If (nonunique) reconstructing functions \' (x, y) are chosen such that their Fourier transforms satisfy xc I Z Z '(2n)1,2 j)\{(2nU."'..'y) = 1. L(4.22) 1=00 i=1 the function s(. y) may be reconstructed from its steerable dyadic wavelet transform by .(,.)= E EyIF 0(y). (1.23) m =x. i=1 where \ (xr.y) denotes \m(0'x. y) rotated by Oi and \,,, (x. y,) = 22" \(2m., 2"' y). To derive an algorithm for the fast computation of the transform, we, similar to (3.3). introduce two smoothing functions such that o(#,..q %v)(&,...q) 5 t i2(2m 2.2wr ) v)(2'rb.. 2>c). (4.24) m=0 i=1 We choose wavelets that are steerable analogs to the onedimensional wavelets from Section 3.1:` /sin( P+"+' ( ,... ) = 0(ju. cos(.'))d (sin( ) (1.25) S2 where L,,. = / + and co = arg(x,awy). These wavelets can be steered with d+1 + dYl basis functions (i.e.. I in (4.16) is equal to d+1l). Choosing o(2w,.) = Ht,(w,),)(.,), (4.26) i"(<.,.o) = s(L,,, ;o ) o(,,.), (4.27) ,(2)) = t,,( r) (a,,.). (4.28) and '(",..O) = KI( 0i)r(). (4.29) 'This choice can be viewed as an extension of the transform presented in [20. 21. 24]. and using (1.26) through (4.29) with (1.21) computed for the finest two scales, we obt a in I C,,(^,.. 9,s, 0)A, ,, ,, 0O) + Ht(.,)Lt(.,) = 1. (4.30) i=1 Setting o(C',.) = I,(,.), and employing (3.5) and (4.25) with (4.26) and (4.27). we find that t ( ,.) = / (,,) (1.31) and ,( ,., 0) = (cos(.0))< ( ,), (4.32) where IIH(') and G(.'() are specified by (3.9) and (3.15). respectively. By inserting (4.31) and (4.32) into (4.30). the missing two filters can be chosen as tL r) =t(,.) (4.33) and t. ,i( ,.. e) = (cos(,( ,.). (4.34) where L(.') and K(.') are given by (3.16) and (3.17). respectively, and (, = Z (cos(' 0, ))24. .I M u ll .... ,1 Dliii I I :, .i t 1 lB , i T ,li i i..1 11i Let us pause here for a brief assessment of tie twodimensional steerable trans forml derived so far. \Ve have chosen steerable wavelets (4.25) which are equal to dth order derivatives of circularly symmetric spline functions in the direction of .raxis (note that knots for thliese splines are circles) and we have laid a foundation for fil ter bank implementations in (4.30). An obvious shortcoming of this scheme is the fact that none of the filter kernels obtained from (4.31) through (1.31) is compactly supported on the rectangular grid. For implementations using digital filters, one is therefore forced to approximate these frequency responses, and by doing so. (4.30) may not hold anymore. Filters in filter bank implementations of steerable pyramids described in [17. s. 39], for example, were designed by using various techniques for approximating desired frequency responses. None of the reported filter banks achieved perfect reconstruction and all filter kernels were nonseparable. Here, we will take a different approach. We will approximate the wavelets in (4.25) in a way that will lead to xy separable filters in a. perfect reconstruction filter bank imple mentation of the transform such that the quality of approximation will improve by increasing the order of Bsplines. Let us approximate wavelets from (4.25) with d/ +^W^L) / /. u(x y) = '+) 3p+d(Y) (4.35) Based oni the fact that Bsplines tend to a Gaussian as their order increases, it is easy to see that both wavelets (4.25) and (4.35) converge to the same functions (i.e., dth order derivatives of the normalized Gaussian in the direction of .raxis) as p + 0. In order to steer wavelets ,'(x,,y) given by (4.35) (note that steering will be only approximate, since these wavelets are not steerable), we need to find basis functions that will approximately steer ,(x,y). Computing rotations, as we did in (4.20). is not practical here, because arbitrary rotations of (4.35) cannot be expressed exactly in terms of central Bspline functions from Chapter 2. Instead, we take advantage of the property of circularly symmetric functions that rotations of their directional derivatives are equal to directional derivatives in rotated directions: P {ad,(X. Y) *p P^(^7)1 d^ ^ where 71Z stands for rotation by 0o, 7Jx,) = iiVo(x, y), o,.(.r. y) is a circularly symmetric function, and 600 denotes vector i7 = (cos 0. sin 0) rotated bv 00. Let us choose P(X y) = p+d(X)!3p+d(Y), which is approximately circularly symmetric function for higher order splines. A rotation of t'(x. y) fronm (41.35) by 00 can therefore be approximated by roaino ,q) (") ,,;^ ^^ U0,( 1) 01t.r, y) 4 ( d d "i 3,+4(,r) d'3,+,1(y) (.6 *''~ 0 r? dfi Yo I rY dxrli dyi(1.) i=' where i = (cos 0o. sin 0o) = (1s., IV). Note that in case of Gaussian, which is both ,ry separable and circularly symmetric. (4.36) becomes exact (e.g., for o(x,y) = c(X2+y2). 00 = 0. and d = {2,4}, we obtain, up to a scaling factor, xy separable basis set for the second and fourth derivative of Gaussian from Tables III and VII of [12]). Having found a set of basis functions (4.36) that approximately steer wavelets (4.35). we want to construct a transform such that Equations (4.20) through (4.24) will be valid superscriptt i must be viewed now as an index, rather than rotation by O). In frequency domain, we can express basis functions from (4.36) as <'i+ (' y) = 's t "_7( )G'_s( v),^p+i( ')i^p+di( '.'/)" /' = 0. 1 ..... d; (1.37) where Gd(s) is given by (3.15) and G'0() = 1. Choosing appropriate '('U .y) to obtain a relation needed for the filter bank implementation of the transform is more complicated than in one dimension. Since we could not find a general solution for arbitrary d, we solve each case separately. Below, we present solutions for the first four orders. When d < 2, we impose 5( c,.,) = o(w'.,.w'y ) = .p( ,).p(,y), a constraint analogous to the one from Section 3.1. For d = 1. we write similar to [28] (.,.,uy) = A'( )T ( ) ,( ),pi(v ). (4.38) 2.,) = r (,)A.( ) ,)p() ( ). (4.39) where A'K(,) and T1(&) are given by (3.17) and (4.11), respectively. Computing (4.24) for the finest two scales and inserting (3.5). (4.37). (4.38), and (4.39) yields a relation between frequency responses ( J1(,)/'(&,)Tl(.9) + T1(.>)G1(,,)KI(.y) + jH(,)H(,)\2 = 1. For d = 2. we choose K 2 1,T2 lO ))p ( )3p2(,y ) KS (.r) Ks(y),Pi.'}p(\ T2(,,,)I\_2 (^}32(^)p^ 3 (",'.&y) \ (t... y (4.40) (4.41) (4.42) (4.43) where T72(0) = IH(,)2. Using (3.5). (4.37). and (4.40) through (4.42) with (4.24) results in G2(x)K2(,x)T2() _+ l(A()'1(,)G(l,)'i (+,) + T2(,.r)G2("y)K2(, )+ +IH( .)H( ,)12= 1. For orders d > 2. we require o(W.,u.'y) = i3p(.r)3p(,y) and 5(u,.,w ) = ('.,.) ,) where ^(w.) is specified by (3.7) and (3.16). For d = 3. we choose reconstructing functions \ ( ,'r .i ) __ A^ 2(,a 1 ..)I',(i^ ) 1_3 ( W^ )13(), r _1( '., 2 . A )(I)Aj (,.;y) 13(w)1 (w ) ) _2(.,) (W,). X W). J)1A(&q)^_3(W.)Y(^~), (4.44) (4.45) (4.46) (4.47) where 1 \3(w) = (1 H(,2), v2 (4.48) and Y_(w') E LP(R) denotes a function such that =(w) = ,}i(')w3_('). i E N. Employing (4.37). (3.5). (3.7). and (1.14) through (4.47) wit h (4.21) yields a relation (;a(&.)K 1 )A'( )1 ,) 2(,)(;2( y ) 32(w., )1 ( )c ) 17 + G(3("q) 3(& )+ G+ I( )L(,)II( )L( ) 1= G1. where L(,;) is specified by (3.16). For d = 41. our choices are (1. I1 ) (4.50) (4.51) (4.52) (4.53) wxv ere 14( .)= i I(.)I. (4.51) I.sing the above functions (4.49) through (4.53). (1.37). (3.5). and (3.7) in (4.24) computed for the finest two scales gives G( (, ) A4( '. ) X(4T, y) + Gl,) .,. )G3 ( ) )K' (3) G(2( )\2(,,)(:)G 2(2')I 1;2(,')V, +('y)+ (,1(,)A'K (,.) W(,y) ()+ +72(,,,)G<1(.^)A'( 1.) + t(&,)L(.,)H(+)L(y) = 1. Here. we have even more freedom for choosing the reconstructing functions than in one dimension. The above functions for d = {2,3.41} were found by trying to imitate lie oniedimiensional transform from Chapter 3 as much as possible. All decoimposit ion filters (;d(,:) were first paired with corresponding reconstruction filters 41 432,,x 2, \4(ua,',:. ,. ( '.^..v,) = 7: (+,,.)1,;( .y)) _=,(*J,.')*( ). Kd(,). and then all other potential digital filters were specified as polynomials in Hs(:). We inserted thus specified filters into a relation between their frequency responses and solved for the unknown polynomial coefficients. Since we allowed more filters with higherdegree polynomials than expected in the solution, the resulting system of linear equations was underdetermined. This allowed enough freedom for removal of undesired digital filters and for balance between degrees of polynomials. The described procedure for determination of reconstructing functions and filters involves quite a lot of heuristics to obtain the appropriate solution from the underdetermined linear system. Unfortunately, we are not aware of any systematic way (aside from numerical optimization, which may be pretty cumbersome) to obtain solutions comparable to the ones above. Next, we will derive a filter bank implementation of the transform. As in Section 4.1. we assume a bandlimited input signal: (,,,'y) = 0 for yj > > or ',l > 7r. Using Shannon's sampling theorem in two dimensions [16] with (4.20) and basis functions from (1.37), we can write I .s(.r,! f) = s(i ,iy)sinc(tf i)sinc(ty iy) G^J^,=00 i,=,X, . Z Sgdsi( '.)3p+i( t ) 9 g(0my)3V+di( ty any) ddt dty, where i = 0, 1 ..... d as in (4.37). Again, we approximate since functions with rorder cardinal splines, then use (2.4), and get D FT rf U) s'=y( 17~.S,, Y } I Y IY (Ux' C,,y) B I (,) B l( ,,y) Bl,+,+i+l(,') =0 Using (4.55) with an approximation Bp+,.+i+l (c)  Bp+,(,')Bi(,(), we can ob tain a filter bank implementation of the transform decomposition. The reconstruction part follows from (4.24), (4.37), and reconstructing functions for distinct values of d. Figure 1.2 shows filter bank implementations of a multiscale spline derivativebased transform for d = {1. 2.3,4}. For d = 1, we recognize (except for the prefiltering and postfiltering) the filter bank implementation of a twodimensional discrete dyadic wavelet transform from [28]. For d = 2, however, our transform differs from the filter bank presented in [22] (i.e., the corresponding transform described in Section 4.1): second derivative is computed only in the directions of ,r and yaxis in [18, 22]. which is not enough for steering. Although not particularly appropriate for orientation anal ysis. such a scheme may still, as we have seen in Section 4.1, efficiently approximate Laplacian of Gaussian across dyadic scales. A transform similar to the one described in this section, was presented in [17. 38.39]. Their filter bank implementation is less redundant (downsampling is used on the lowpass branch, while simultaneously keeping aliasing negligible by using a filter with insignificant amount of energy for l\,j > ) and uses reconstruction filters with same iii.i.,iitude frequency responses as the decoinposition onesa possible advantage for certain applications. They have, on the other hand, little control over the function from which derivatives are computed (to obtain a dth derivative, they multiply a circularly .. i liit'tric filter by (jcos 0)' with all filters being obtained by recursive minimization of a weighted combination of constraints), filter bank does not have perfect reconstruction, and none of the filters is rxy separable. Br'(Cx)B;Coy) Bp+r+i(Owx) Bp+,,1(Co) BR1. )J Bp+r+i((Oy) (b) Gs(2mox) Hs(2m),) Hs(2m(Oy) (c) G(2mox) G.,(2mx)J GCs(2 (y) G2(2m(Oy) Hs(2eo) H(2my) (e) Br((o),) Br(o)y) Ks(2mmo) T(2m(oy) Tl(2m%,,) K~s(mwy) Ls(2mox) Ls(2moy) (d) K(2mo),) T(2mOy) Kl(2 mox) Ks(2moy) T2(2m(o) K(2 %)) Ls(2mx) Ls(2mwy) (f) Figure 4.2: Filter bank implementation of a multiscale spline derivativebased trans form for inm E [0. 31 1]: (a) Prefiltering, (b) postfiltering, (c) decomposition and (d) reconstruction modules for d = 1, and (e) decomposition and (f) reconstruction modules for d = 2. B1 ) G.(2m(o) B2(2moy) K(2mco,) Bf(2oy) G2(2mO,) Gs(2(o,) K2(2mCo,) V3(2m(,x) K(2m(O,) V3(2m(o) G_,(2m0)) G2(2mcoy) Ks(2mCo,) V3(2mo,) K2(2moy) V3(2m(oy) B2(2m(o,) Ga(2mwy) 2'{2m(,,,) K3(2m(O) H,(2mo,) H,(2m(o,) Ls(2m(o,) Ls(2m0y) (g) (h) G4(2m(),) B3(2m, K4(2mo),x) B3(2mO)y) T2(2m(oy) Gs(2mj Q Gs(2m(ov) B(2mw))^) Ks(mw,) Ej2m(2 )K(2m(0y) G2(2m(,) G2(2mw,) K2(2m),) V4(2mco) K2(2mw ) V4(2'ko,,) G.s(2mo),) B2(2m(o,) G,(2moy) B(2mo,) Ks(2m(o) K3(2my) B3(2mo),) G4(2mwy) B93'(2mCo,) T2(2%m(o,) K4(2mwo) H(2mo,) H(2mo,) Ls(2mo,) Ls(2m)y) (i) (0) Figure 1.2: (Continued: (g) Decomposition and (h) reconstruction modules for d = 3, and (i) decomposition and (j) reconstruction modules for d = 1. Decomposition modules are recursively applied at I lthe locations of the filled circles. CHAPTER 5 IMPLEMENTATION ISSUES Except for the steerable dyadic wavelet transform, all transforms presented in Chapters 3 and 4 can be implemented as filter banks consisting of onedimensional fil ters only. In this chapter, we show how to take advantage of symmetry/antisymmetry of filters when combined with a mirror extended input signal. 5.1 Finite Impulse Response Filters Since all twodimensional filters used in the filter bank implementations of the transforms from Chapter 4 are xr/y separable, we will begin this section with a detailed description of FIR filter implementations for the onedimensional discrete dyadic wavelet transform implementation from Figure 3.1. The extension to two dimensions will then be straightforward. Let us refer to filters applied at scale 2' as filters at level nm+1, and let filters at level 1 (Equations (3.9). (3.15) through (3.17). (4.11). (4.43). (4.48). and (4.54)) be called "'original filters." to distinguish them from their upsampled versions. As an input to the filter bank from Figure 3.1, we consider a real signal s(n) E 12(Z). n C [0, A 1]. Depending on the length of each filter impulse response, filtering an input signal may be computed either by multiplying the discrete Fourier transforms of the two sequences or by circularly convolving s(n) with a filter's impulse response.1 Of course, such a periodically extended signal ii,,. change abruptly at the boundaries causing artifacts. A comminon remedy for such a problem is realized by constructing 1As is customary in image processing, we use circular, rather than linear, convolution. a mirror extended signal [18] s( 1? 1) if n e [ 1] . .;() { L( ) if n E [0, 1]. (5.1) where we chose the signal .s (n) to be supported in [NA, 1]. It will become evident shortly, that mirror extension is particularly elegant in conjunction with symmetric/antisymmetric filters. Let us first classify synmmetric/antisymmetric real evenlength signals into four types [30]: Type I f(n) = f(n), Type II f(n) =f(n 1). Type III f(n) =f(n), Type IV f'(n) = f(n 1). where i E [A.V 1]. Note that for Type I signals the values at .f(0) and f(N) are unique, and that for Type III signals the values at f(0) and f(NV) are equal to zero. Using properties of the Fourier transform it is easy to show that the convolu tion of symmetric/ai ,it. iii ,.tric real signals results in a symmetric/antisynmmnetric real signal. If a symmetric/antisymmetric real signal has an even length, then there always exists an integer shift such that the shifted signal belongs to one of the above types. Now. we are ready to examine the filter bank implementation of a onedimensional discrete dyadic wavelet transform from Figure 3.1 with filters given by (3.9) and (3.15) through (3.17) driven by a mirrored signal Ism,,(n) at the input. Let the number of levels 3M be restricted by 1 < 1 + log2, (5.2)1 Lmar I where Lr, is the length of the longest original FIR filter impulse response. Eachli FIR filter block in the filter bank consists of a filter and a circular shift operator. Equation (5.2) guarantees that the length of the filter impulse response does not exceed the length of the signal at any block. Since our input signal ...(n) is of Type II and noninteger shifts at level 1 are rounded to the nearest integer, it follows that a processed signal at any point in the filter bank belongs to one of the types defined al)ove. This means that filtering a signal of length 2A can be reduced to filtering a signal of approximately one half of its length. (For Types I and III. N + 1 samples are needed. However, for Type III one needs to store only N 1 values because zero values are always present at the zeroth and ( \)th sample position). Implementation is particularly simple for FIR filters designed with d even and p odd. Filters are of Type I in this case. so the signal at any point of tlie filter bank will be of Type II. An FIR filter block from the filter bank shown in Figure 3.1 can therefore be implemented by L1 F!,,,u,(n) = f(0)uj,(n)+y .f(')[Itjj(ii2'"'i)+u11(i+2"'i)]. nC e [0. A1]. (5.3) i=t whIere ( ) ifn C[ .I] /(,,) = u (,) if n C [0. X ] (5.1) S(2N I) if NC [ 3)']. a(n) is an input signal to a block. f(I?) is an impulse response of some original filter. L is the length of the filter, and AN is the length of an input signal s(ii) to thlie filter bank. Implemnientation of filters b1 (u) used for prefiltering and postfiltering represents a special case of (5.3) with min =0. A filter bank with the above implementation of blocks and signal s(n ) at the input yields equivalent results as circular convolution for s,,,(n) as defined by (5.1). In addition to requiring one half the amount of memory, the computational savings over a circular convolution implementation of blocks are. depending on the original filter length, three to four times fewer multiplications and one half as many additions. A similar approach can be used for other filters. However, things get slightly more complicated in this case. because the filters are not of the same type and the signal components within the filter bank are of distinct types. As a consequence, an implementation of blocks that use distinct original filters may not be thile same, and the implementation of blocks at level 1 may differ from the implementation of blocks at other levels of analysis. The decomposition blocks at level 1 can be implemented by }  2 (,.OU(n) = g(i)[uII(n i 1) u11(n (+i)]. + E [1, A 1]. i=o for d odd. (5.3) for d(I even. L1 2 H_,,oii(n) = h(i)[uII(n I 1) + uii(rn +i)]. EC [O. A]. i=0 for p even. and (5.3) for p odd, where ujj(1) is defined by (5.4). g(n) and h(n) are impulse responses of the filters computed from (3.15) and (3.9), respectively, and L is the length of the corresponding impulse response. The output from a block G_,(w) at level 1 is of Type III for d odd and of Type II for d even. while the output from H,(,) at the same level is of Type I for p even and of Type II for p odd. The decomposition blocks at subsequent levels nm E [1, .1 1] can be imple miented by /  2 G,, u(n) ="( (i)[u(n 2m(i +.s)) u((n, + 2"(i + s))], n E [1. N' I], i=o for d odd and p even. Li 21 (G_.,,,u(n) = g(i)[un( 2"..(I + s)) uii(n + 2"(i + s))]. n C [0. A' 1]. i=0 for d and p odd. L F.1u(n) = f(O)I((n) + > f(i)[uI(u 2"t1) + it(n + 2l)]. ( + [0. A], (5.5) i=1 with /(n) = g(nt) for d and p even, .L1 s.mu(n) = h(i/)[u/i(n 2"1(i + .s)) + utI(n + 2m(i + .))], e [0. ]. (5.6) i=0 for p even. and (5.3) for p odd, where 1 1(n) ift E [ . 1] 2" uI(o) = j (n) if a C [0, N] (5.7) a(2N a) if E [N + 1, ]. 2 " The outputs from blocks Gs(2'mw) are of Type III for d odd and p even, of Type IV for d and p odd. and of Type I for d and p even, whereas the outputs from H_,(2'2";) are of Type I for p even and of Type II for p odd. N .:i. the reconstruction blocks at level 1 can be implemented by L 2 Kx.ou(n) = + k(i)[1n)( i+) utiI(n+ i)]. + C [0, N 1]. i=1 for d odd, (5.3) for d even, L 2 L,,oa(,) = l(i)[ui(n i + 1) + u(n +i)]. n E [0, A 1]. i=1 for p even. and (5.3) for p odd, where Su(n) if I? [. 1] 0 if 7, = 0 tI(n(n) u(= () if I [1. N 1] (5.8) 0 if t = N a(2N n) if E [N+1,j. uI(n) is as defined by (5.7) and k(n) is an impulse response of the filter from (3.17). Note that both outputs from blocks K,(() and Ls(') are of Type II. The reconstruction blocks at subsequent levels can be implemented by L1 2 A'2.,,1u() = k(i + 1)[aI(r? 2'(i + s)) ilii(n + 2' (1 .s))]. a E [0, N']. i=o0 for d odd and )p even. (5.5) with f(ii) = k'(n) for d and ) even. L 21 ', ..,,,(,) = k A.(i + l)[u,(, 2"(i + .)) n ( + 2'(i 2+ .))]. c [0. A I]. i=0 for d and J) odd. L.m,,u(lu) = H_ ,,"(.). for p even, and (5.3) for /) odd. where iui2(l) is given by (5.8). it(1) ift E[E .1] iy(n,) = u(nI) if E [0. 1] u(2N it 1) ifn [X, 3]. and Its.,, u(n) is specified by (5.6). We observe that thei outputs from blocks K.5(2"',) and L,(2`',). m G [1. 1 1]. are of Type I for p even, and of Type II for p odd. When we compare the above implementation of blocks to circular convolution driven by a mirrored signal ,, (u) at the input, we observe that approximately twofold less memory space, three to four times fewer multiplications and one half as many additions are required. (F1or Type I signals an additional sample has to be saved because two values are without a pair). Until now. we have talked only about the onedimensional case. whose filter bank implementation is depicted in Figure 4.2. Twodimensional transform filter bank implementations (Figures 4.1 and 4.2) are not only comprised of .rxy separable filters solely, but also use all the filters from Section 3.1. Everything presented in this section so far is therefore directly applicable to the twodimensional case. Filters which have not been treated vet (i.e.. ti(n), t2(1) t'3(0), 1'4(0). and filters b,(n) from the decomposition modules of Figure 4.2) can all be realized by (5.3) for p) odd or m = 0 and (5.5) otherwise (,f(n) denotes an impulse response of any of the above mentioned zerophase filters in this case).2 2'In case of filters u'3(n) and v.((n). a slightly more efficient implementation can be obtained by precomputing new filters k t3(n) and k v40(), and then implementing then by K, ,u(n), in E [0. M 1]. The implementation presented( in this section performs all operations in the spatial domain, however, one could also implement the structures shown in Figures 3.1. 1.1. and 4.2 with anll input signal .s,, (Mn) (Equation (5.1)) in the frequency domain. For short filter impulse responses, such as those given in Tables 3.1. 3.2 and 3.3. the spatial implemnentat ion described in this section is certainly more efficient. For long filter impulse responses, however, filtering is faster if imlplenmented in the frequency domain. Additional (let ails on alternative implementation strategies c('an be found in [33]. 5.2 Infinite Impulse Response Filters Implelnentation of IMI filters b. (n1) which were introduced in Section 2.2 is a bit more involved than the one encountered in the previous section. Fortunately, the number of different cases is much smaller here: possible input to b1(n,) in filter banks from Figures 3.1, 4.1. and 4.2 is either of Type II or of Type I.3 We will use ideas and a few results from [43]. Let us first take a closer look at the system function B1(z). This function can be written as a cascade of terms 1 n )+ (1 +z (") which can be expressed in a parallel form as E()+ ( +1 1 i). (5.10) E )=1  \1 az1 +lacz where o and are poles of tlie causal and the anticausal filter, respectively. The impulse response of this term can be writ tenll as c(n) o)i . 1 a2 3Note that symmetry types in this section slightly differ from those defined in Section 5.1: here. mirror extended signals are periodically repeated, so that they stretch from x to xc. We choose to implement E(z) in a cascade form and therefore extract the difference equations from (5.9): c+(n) = u(n) +oc c+( 1) n = 1.2.....  1. (5.11) c(n) = (c(;? + 1) c+(n)) i = X2, N3.....0. (5.12) where u(n) denotes the input to tlie block, c+(n) is the output from the causal part, and c'(i) stands for the output from the block. To solve (5.11) and (5.12) we need boundary conditions c+(0) and c(AV1). Let us begin with filters I (n) in filter bank implementations from Figures 3.1. 4.1, 1.2(a). 4.2(b). and Figures 4.2(h) and 4.2(j) with n =0. We derive 0 .\'1 01+1 + o 2.Xi 10 C+(0)= 1 ' )l( = 1(O)+ 12 + (i) ( 0)+ +l(0. (5.13) =x i=0 i=O and. using parallel form (5.10) N1 \ N+ i 00+0 +I + c(A 1) 1 (C+ 1) + 0 + (i)) i=o  (+(+ I i,(I)). 1=  ^ I=Nli0 (5.11) w liere f uJJ( mod (2A)) if > 0 ul1p"() ~ { "/((n + 1) mod (2A')) if n < 0. 1f (n) if I G [O.x 1] uj"( u(2 n1) ifn C [N. 2X 1]. A is I lie length of an input signal to the filter bank, and i0 < A I is selected such that o' falls below a predefined precision threshold. For IIR filters from Figures 4.2(h) and 4.2(j) with ni C [1. M 1] and 1 odd. we get 0 C,(0) = 2[ 1 111p()],,  ^EE^ ^]t + [C2I^1}( 1 =(0) + E J uQ) 11=0 i=O (5.15) and a 11(1 c(X1) 0Nx1 C(A 1)  (C+ (N ) E ]+o u( (5.16) o n=O i=0 where 0 {" if E Z S 0 ot lierwise. If N is a power of two and 2" < N ((5.2) guarantees I hal the latter condition is always true) (5.15) becomes r ,+I] 2V__,l .\1 [o2,, + [o ,, ,() S+(o) = ,(0) + V L c, I" i=0 I 02 while (5.1(i) can be written as N r .V, 1] f .0j. L cA* a1)^ + 1 2N j c(.\l) 1 ^ ( E   k () =0 0 1 o5[ linallv. the boundary conditions for filters bl(n) from Figures 1.2(h) and 4.2(j) with tin E [1, 31 1] and p even are o C+ (0) =^ [Up?] [o=x,/  y a ,,, (0)+ o ,' "(k')+ 1n=O'O + 71=0 + 2 ] + [ u(l) (5.17) and O C(N) (2c+(N) u(N)). 1 a2 whIere Uy(1 ) = ul(n nrod (2N)). f u(n) if n e [0, N] it, ) u(2N n) if n E [. + 1.2N 1]. and c(:N) = c+(.') with c(n) denoting response of the anticausal filter from (5.10) was used. Again, if V is a power of two and 2' < N. (5.17) can be simplified (())+ ^AQ+E^1^ 4, + 2N, t(> u(0) + 2It'U(N) + i=1 0{ 21 [1 2,\'1 C,+(0) =  i  1 (n' Series in expressions for c+(0), c(A 1). and c(AV) for filters from Figures 4.2(h) and 1.2(j) with in E [1. A1 1] can be, similar to (5.13) and (5.14). truncated according to the desired precision. For orders p) greater than three, we implement B'(:) as a cascade of terms E(:) with different n's. Note that the output from block E( ) is always of the same type as the input to it. CHAPTER 6 APPLICATIONS In this chapter, we present a comparison between multiscale spline derivative based transform and traditional wavelet techniques on two sample applications. 6.1 Image Fusion Image fusion combines particular aspects of information from the same imag ing mnodalitv or from distinct imaging modalities and can be used to improve the reliability of a particular computational vision task or to provide a human observer with a deeper insight about the nature of observed data. Whether it is combining different sensors or extending the dynamic range of a single sensor, the goal is to achieve more accurate inferences that can be achieved by a single sensor or a single sensor setting. The simplest method of fusing images is accomplished by computing their average. Such a technique does combine features from input images in the fused image, however, the contrast of the original features can be significantly reduced. Among more sophisticated methods, multiscale and multiresolution analyses have become particularly popular. Different pyramids [5. 41] and waveletbased techniques [7, 20. 26, 32] have been applied to this problem. Fusion is performed in the transform space by computing local statistics across the decomposition scales. Typical size of neighborhood is between a single pixel and 5 x 5 area with some loss of contrast reported for the latter [5]. Some authors argue about of the particular criteria for fusion [5, 26]. while others propose a set of different combination rules to perform image fusion in a variety of ways [7]. Here. we are not trying to find the best criteria for fusion. We simply use a "maximum magnitude" rule (i.e., at each position and scale of the transform the corresponding transform coefficients are compared and those with greater magnitude included for reconstruction) and compare the results obtained by traditional wavelet techniques with the ones produced by the multiscale spline derivativebased trans form. Figure 6.1(a).(b) shows a pair of images with distinct (lower in Figure 6.1(a) and upper in Figure 6.1(b)) areas blurred. Figure 6.1(c) demonstrates the ideal result of image fusion obtained by manually cutting and pasting the upper part of Figure 6.1(a) and the lower part of Figure 6.1(b). Image fusion resulting fromin the use of discrete wavelet transforms is shown in Figure 6.1(d) for the orthogonal wavelet DAUB12 [7] and in Figure 6.1(e) for the linear phase biorthogonal wavelet Bior6.8. Figure 6.1(f) illustrates the result of image fusion with a multiscale spline derivative based transform being employed. Note that the image fused using the multiscale spline derivativebased transform is virtually indistinguishable from the one obtained manually. Both orthogonal and biorthogonal wavelet transform exhibit artifacts. On the average, biorthogonal transform performs better than orthogonal due to the linear phase of wavelets and filters, but still causes serious artifacts. Basing the decision rule on an area rather than on a single pixel certainly reduces artifacts of traditional wavelet methods, however, it reduces sharpness as well. As noted by Burt and Kolczynski [5], "1tli,'i' is a tradeoff between sharpness and shift invariance." Here, we argue the point that the proper choice of transform can eliminate artifacts due to translation and rotation invariance altogether. Another possible source of artifacts is misregistration. Lowlevel fusion al gorithms typically require input images that have already been properly aligned. Misregistration causes artifacts [7]. and it can be recognized from our discussion in Section 1.1 that they are more obvious for translation and rotationnoninvariant wavelet techniques than for the multiscale spline derivativebased transform. 6.2 lli,,_' F lil,,i li, ii.li Similar to the previous section, we will point out problems associated with traditional wavelet analysis on an example. Below, we briefly present the problem of contrast enhancement in digital maninography. In mammography, early detection of breast cancer relies upon the ability to distinguish between malignant and benign mammographic features. Contrast en hancement can make more obvious unseen or barely seen features of a mammograms without requiring additional radiation. The complexity of mammographic images and the subtlety of malignancies can present a challenge even to expert radiologists. In addition to dealing with barely visible mammographic features, human observers sometimes simply overlook abnormalities. Such mistakes affect the number of false negative cases considerable. Computer processing of digital mammograms can assist radiologists to reach a correct diagnosis more consistently. A variety of computer based techniques have been reported in almost three decades of research. The advent of direct digital mammography devices has made digital image processing techniques more attractive for screening. Contrast enhancement in a wavelet transform framework can be achieved by applying a (typically nonlinear) function to wavelet coefficients and then reconstruct ing an enhanced image with modified coefficients. We are not trying to find the best method for a specific type of images or noise here, rather we want to illustrate pos sible sources of artifacts when orthogonal and biorthogonal wavelet transforms are used. Figure 6.2(a) shows an original mammographic image. Figure 6.2(b) depicts the result of enhancement using piecewise linear enhancement function [11, 23] to (c) (d) (e) (f) Figure 6.1: (a) An image with lower part blurred. (b) An image with upper part blurred. (c) Fused image obtained by combining images from (a) and (b)) mian ually. Results of fusion performed by (d) discrete orthogonal wavelet transform (DA!'BI2). (e) discrete biortliogonal wavelet transform (Bior6.8). and (f) multiscale spline derivativebased transform. Ow modify wavelet transform coefficients, and Figure 6.2(c) shows the result of enhance inent using the multiscale spline derivativebased transformin and the same enhance ment function as Figure 6.2(b). In this figure, we. more than trying to point out better enhancement results in Figure 6.2(c). observe that Figure 6.2(b) sliowxs fea tures (artifacts) thai are not present in thlie original image (they are particularly obvious in the area right from the mass in Figure 6.2(b)). Here, the goal of enhance ment is making visible features that are present in the original image but obstructed by different structures within the same image. Generation of artificial features, while not being rare when nonlinearly processing coefficienlts of thlie orthogonal or biort hog onal wavelet transform, is certainly not desirable for this particular application. For more details on different enhancement strategies in minammnography. please refer to [11, 22.2:, 2.5]. Another problem related to enhancement and other kinds of processing of dig ital minamniogramis (and not unique to mammography) is that the same mammiogram could have been easily translated and/or rotated. From obvious reasons, it is not desirable that the result of processing can be significantly affected by positioning only one more reason for a translation and rotationinvariant transform. Figure 6.3 shows outputs from steerable basis filters applied to the image from Figure 6.2(a) at scales {1.2. 4}. By linearly combining these outputs analysis along arbitrary direction can be performed (cf. Chapter 4). This enables, as illustrated in Chapter 1. a rotationinvariant processing [21]. (b) (c) Figure 6.2: (a) An original mammographic image containing a spicular mass. An enhanced image using (b) discrete biorthogonal wavelet transform (Bior6.8) and (c) multiscale spline derivativebased transform. I:'  Figure 6.3: The magnitude of filter coefficients for (three levels of analysis shown). (;0( w.^/) at oi {O. 7. 3} CHAPTER 7 CONCLUSION In this thesis we constructed a new transform that does not introduce artifacts due to translation and rotation invariance, which are inherent to traditional wavelet analysis. We extended the onedimnensional discrete dyadic wavelet transform to higher order derivatives and evenorder spline functions and developed an improved initial ization procedure. Comparison to the originally employed initialization [28] showed significantly better performance of our procedure for finer scales of analysis. We developed several twodimensional transforms. All of them were derived with the goal of eliminating artifacts due to lack of translation and rotation in variance. \We presented a twodimensional discrete dyadic wavelet transform with a firstderivative wavelet as an extension of the one originally proposed in [28], a two dimensional discrete dyadic wavelet transform with a secondderivative wavelet that can approximate Laplacian of a Gaussian, and a steerable dyadic wavelet transform implemented in a nearperfect reconstruction filter bank which may be preferred when there is a need for orientation processing along equally spaced angles. WNe derived a multiscale splinie derivativebased transform which uses xy separable filters in a perfect reconstruction filter bank and enables fast translation and rotationinvariant directional analysis of images. We empirically showed the presence of artifacts in image fusion and enhance ment applications when traditional wavelet methods were used. Here developed trans forms did not exhibit similar artifacts, and, in case of fusion, yielded nice results without the use of postprocessing. 67 Future work will consent rate on a more thorough exploitation of the transformi space for the existing applications, and potential new applications (e.g.. denoising). In particular, we will try to use the causality property of a Gaussian kernel mentioned in Chapter 1 to trace features across scales and process them according to their mnult iscale behavior. REFERENCES [1] A. Aldroubi. M. Unser. and M. Eden, "Cardinal spline filters: Stability and convergence to the ideal since interpolator," Signal Process.. vol. 28, no. 2. pp. 127 138. 1992. [2] J. Babaud. A. P. Witkin, MN. Baudin, and R. 0. Duda. "Uniqueness of the Gaussian kernel for scalespace filtering," IEEE Trans. Pathtrn A1nal. IMach. Inh1Il.. vol. 8. no. 1. pp. 2633, 1986. [3] C. de Boor. A Practical G(uid to 5plincs, Springer Verlag. New York. NY, 1978. [4] P. J. Burt and E. I1. Adelson. "The Laplacian pyramid as a compact image code." IEEE Trans. Comm an., vol. 31, no. 4. pp. 532 540. 1983. [5] P. J. Burt and R. J. Kolczynski, "Enhanced image capture through fusion." in Proc. Int. ('onf. Comprnut. Vision, Berlin. Germany, 1993. pp. 173182. [6] J. ('annv. "A computational approach to edge detection," IEEE Trans. Pattern Anal. Mach. Intell.. vol. 8. no. 6, pp. 679698, 1986. [7] L. J. Chipman, T. M. Orr, and L. N. Graham. "'WVavelets and image fusion." in Proc. IEEE Int. ('of. II.,,ii, Process., Washington. D.C.. 1995. vol. 3, pp. 248251. [8] C. K. Chui, An Introduction to Wavclets, Academic Press, San Diego. CA. 1992. [9] C. K. Chui, Ed.. IVarvelts: A Tutorial in Theory and Applications. Academic Press. San Diego. C(A. 1992. [10] I. Daubechies, Ten Lectures on Wavelcts, SIAM. Philadelphia. PA. 1992. / [1 J. Fan and A. Laine. "Mlultiscale contrast enhancement and denoisig in digital v radiographs." in Wiavelts in Medicine and Bi,. A. Aldroubi and M. Unser, Eds.. CRC Press. Boca Raton, FL, 1996. pp. 163189. [12] W. T. Freeman and E. H. Adelson, "The design and use of steerable filters," IEEE Trans. Pattern Anal. Mach. Intell., vol. 13. no. 9. pp. 891 'iit,. 1991. [13] H. Greenspan, S. Belongie. R. Goodman. P. Perona., S. Rakshit. and C. H. Anderson. "Overcomplete steerable pyramid filters and rotation invariance." in Proc. IEEE ('omputt. Soc. C'onf. C('omnput. Vision Pattern Recognit.. Seattle, WA. 1994, pp. 222 228. [14] Y. HelOr and P. ('. Teo. "Canonical decomposition of steerable functions," in Proc. IEEE ('omptut. Soc. C'onf. Comnput. IVision Pattern I? .,i,.t., San Fran cisco, CA, 1996. pp. MIV'I s16. [15] M. Holschneider. R. KronlandMartinet. J. Morlet, and Ph. Tchanmitchian, "A realtime algorithm for signal analysis with the help of the wavelet transform." in lai'arlts. 'imeFrequei i' / Mthods and Phas( Spac(. J. M. Comibes, A. Gross mann. and Ph. Tchamitchian. Eds., SpringerVerlag. Berlin. GermanY. 1989. pp. 286297. [16] A. J.. Jerri. "The Shannon sampling theoremits various extensions and appli cations: A tutorial review," Proc. IEEE, vol. 65, no. 11. pp. 1565 1596. 1977. [17] A. Karasaridis and E. Simoncelli, "A filter design technique for steerable pyramid image transforms." in Proc. IEEE Int. Conf. Acoust. Speech Signlal Process.., Atlanta. GA. 1996. vol. 4. pp. *23.' 2:392. [18] I. Koren and A. Laine, "A discrete dyadic wavelet transform for multidimensional feature analysis." in TimnrFrcqucncy and lavelet Transformns in Bionmedical E,,,", .:,.' 1.q. M. Akay. Ed., IEEE Press. New York. NY. 1997. [19] I. Koren. A. F. Laine. J. Fan, and F. J. Taylor. "Edge detection in echocardio graphic image sequences by 3d multiscale analysis." in Proc. IEEE I0t. Conf. Imay( Process.. Austin. TX, 1994, vol. 1, pp. 288 292. [20] I. Koren. A. Laine. and F. Taylor, "Image fusion using steerable dyadic wavelet transform." in Proc. IEEE Int. Conf Imnag Process.. WVashington. D.C., 1995, vol. 3, pp. 2:322:35. [21L I. Koren. A. Laine. F. Taylor. and M. Lewis, "Interactive wavelet processing and techniques applied to digital mammography," in Proc. IEEE Int. Conf. Acoust. /Spe5ch Signal Procss.. Atlanta, GA, i',. vol. 3, pp. 1415 1418. ""] A. Laine, J. Fan. and S. Schuler, "A framework for contrast enhancement by v', dyadic wavelet analysis." in Digital Mammography. A. G. Gale. S. MI. Astley, D. R. Dance, and A. Y. Cairns, Eds.. Elsevier. Amsterdam, The Netherlands, /1991, pp. 91 100. "] A. Laine. J. Fan. and W. Yang, "Wavelets for contrast enhancement of digital mammography." IEEE Eui, Med. Biol. Mag., vol. 14. no. 5. pp. 536 550, '',. [24] A. Laine. I. Koren. WV. Yang, and F. Taylor. "A steerable dyadic wavelet transform and interval wavelets for enhancement of digital mammography," in Wa(Telet Applicationis II. Proc. SPIE, H. H. Szu. Ed.. Orlando. FL, I''., vol. 2491, pp. 736749. i_24 A. F. Laine. S. Schuler. J. Fan, and W. Huda, "Mi.iiiographic feature en hancement by multiscale analysis," IEEE Trans. Med. Imaging, vol. 13. no. 4, pp. 725 740. 1994. [26] H. Li, B. S. M.Liiijunath, and S. K. Mitra, '.ti~ isensor image fusion using the wavelet transform," in Proc. IEEE Int. ('onf. Image Process., Austin. TX, 1994, vol. 1, pp. 5155. [27] S. Mallat and W. L. Hwang, "Singularity detection and processing with wavelets." IEEE Trans. Inf. Theory, vol. 38, no. 2. pp. 617 643, 1992. [28] S. Mallat and S. Zhong, "'Characterization of signals from multiscale edges." IEEE Trans. Pattern Anal. Mach. Intll.. vol. 14, no. 7. pp. 710732. 1992. [29] D. Marr and E. Hildreth, "Theory of edge detection." Proc R. .Soc. London Ser. B. vol. 207. no. 1167. pp. 187217, 1980. ill] A. V. Oppenheim and R. W. Schafer, DiscreteTin. Signal Procs.sing, Prentice Hall, Englewood Cliffs, NJ, 1989. [31] P. Perona. "Deformable kernels for early vision." IEEE Trans. Pattern Anal. Mach. Intell., vol. 17, no. 5. pp. 4s. I'n', 1995. [32] T. Ranchin. L. Wald. and M. Mangolini, "Efficient data fusion using wavelet transform: the case of SPOT satellite images," in Mathematical Imaging: 'atvela t Applications in Signal and Image Pi .. f ,',. Proc. SPIE. A. F. Laine. Ed., San Diego. CA. '1,. vol. 2034, pp. 171 178. [33] 0. Rioul and P. Duliamel, "Fast algorithms for discrete and continuous wavelet transforms." IEEE Trans.. Inf T1, "ry. vol. 38. no. 2. pp. 569586. 1992. [34] I. J. Schoenberg. "'Contributions to the problem of approximation of equidistant data by analytic functions," Q. Appl. 11,/0, vol. 4, no. 1. pp. 45 99. 112141. 1946. [35] I. J. Schoenberg. "C(ardinal interpolation and spline functions." .J. Appro.r. Thory, vol. 2. no. 2, pp. 167206, 1969. 1i"" I. J. Schoenberg. "Notes on spline functions III: On the convergence of the interpolating cardinal splines as their degree tends to infinity," Isr. J. MAlath., vol. 16. no. 1. pp. 8793, 1973. [37] C. E. Shannon. "Communication in the presence of noise," Proc. IRE. vol. 37. no. 1, pp. 10 21. 'i. [38] E. P. Simoncelli and W. T. Freeman, "The steerable pyramid: A flexible archi tecture for multiscale derivative computation." in Proc. IEEE Int. ('onf. Imnag Process.. Washington, D.C., 1995, vol. 3. pp. 444447. [39] E. P. Simoncelli. W. T. Freeman, E. H. Adelson. and D. J. Heeger. "Shiftable multiscale transforms." IEEE Trans. Iinf. I/i.',,, vol. 38. no. 2. pp. 587607, 1992. [40] M. .1. T. Smith and T. P. Barnwell, III, "Exact reconstruction techniques for treestructured subband coders." IEEE Trans. Acoust. cSpech S.5ignal Process.. vol. 34. no. 3 pp. 434441, 1986. [41] A. Toet. "Image fusion by a ratio of lowpass pyramid," Pattern Rcognlit. Lett., vol. 9. no. 4. pp. 24525:3, 1989. [42] M. User and A. Aldroubi. "A general sampling theory for nonideal acquisition devices." IEEE Trans. Signal Process., vol. 42. no. 11, pp. 2915 2925. 1',i, . [43] M. Unser. A. Aldroubi. and NI Eden, "Fast Bspline transforms for continuous image representation and interpolation," IEEE Trans. Pattern Anal. Mach. Intll., vol. 13. no. 3. pp. 2772'., 1991. [44] M. User, A. Aldroubi, and M. Eden, "On the asymptotic convergence of 13B spline wavelets to Gabor functions," IEEE Trans. Inf Theory. vol. 38, no. 2, pp. %, 1872. 1'l1i1'. 71 [45] M. Unser, A. Aldroubi. and M. Eden, "Polynomial spline signal approximations: Filter design and asymptotic equivalence with Shannon's sampling theorem," IEEE Trans. IJf. Theory, vol. 38, no. 1. pp. 95103. 1992. [46] M. Unser. A. Aldroubi, and S. J. Schiff, "Fast implementation of the continuous wavelet transform with integer scales," IEEE Traits. Signal Process.. vol. 42, no. 12. pp. 35193523. 1991. [47] M. Vrhel. C. Lee. and M. Unser, "Fast computation of the continuous wavelet transform through oblique projections," in Proc. IEEE Int. ('of. Acoust. ,5ptch Signal Process., Atlanta. GA. "'l'i. vol. 3, pp. 14591462. [48] A. Witkin. "'Scale space filtering," in Proc. Int. Joint Conf. Artif. Inf l., Karlsruhe. Germany. 1983, pp. 10191022. BIOGRAPHICAL SKETCH Iztok Koren earned the B.S. and M.S. degrees in electrical engineering from the University of Ljubljana, Slovenia, in 1987 and 1991. respectively. He will receive the Ph.D. degree in electrical and computer engineering from the University of Florida. Gainesville. in 1' .. His research interests include image processing, digital signal processing, and wavelets. He is a member of Eta Kappa Nu, Tau Beta Pi, and the IEEE. I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Docor of P ilosophy. Fred J, T.h1 Chairman Professor of electrical and Computer Engineering I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy. "\ T c, .' _ Andrew F. Laine Cochairmani Associate Professor of Computer and Information Science and Engineering I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation ,, ll adequate, in scope and quality, as a dissertation for the degree of Dotorof/PF 'ophb S.l,,,, (,'. P n l C' ~of Eletrical and Computer I Er'iiieering I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doc of Philosophy. yJohn M. M. Anderson yAssistant Professor of Electrical and S Computer Engineering I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality. as a dissertation for the degree of otor ofPhisophy. Anmit Sigmon 1 Associate Professor of Mathematics This dissertation was submitted to the Graduate Faculty of the College of Engineering and to the Graduate School and was accepted as partial fulfillment of the requirements for the degree of Doctor of Philosophy. December 1996 Winfred M. Phillips Dean, College of Engineering Karen A. Holbrook Dean, Graduate School LO 1780 UNIVERSITY OF FLORIDA 3 1262 08555 0696 
Full Text 
xml version 1.0 standalone yes
Volume_Errors Unscanned PreviousPageID P110 