A multiscale spline derivative-based transform for image fusion and enhancement


Material Information

A multiscale spline derivative-based transform for image fusion and enhancement
Physical Description:
Koren, Iztok
Publication Date:

Record Information

Rights Management:
All applicable rights reserved by the source institution and holding location.
Resource Identifier:
aleph - 26349757
oclc - 36680276
System ID:

Full Text







Copyright 1996


Iztok Koren

:Iojimn starsem


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


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.



LIST OF TABLES ...............

LIST OF FIGI'RES ..............

ABSTRAC(T ...................


1 INTRODUCTION ...........

1.1 Shortcomings of Traditional \l,.w
1.2 Thesis Overview .........
1.3 Notation .............

1,,, of Wavelet Analysis


2.1 Central B-Splines: Definition and Properties . .
2.2 B-Spline Signal Interpolations . . . . . .
2.3 B-Spline Signal Approximations . . . . . .


3.1 1-D Discrete Dyadic Wavelet Transform Revisited .
3.2 Rem arks . . . . . . . . . . . .


4.1 2-D Discrete Dyadic Wavelet Transform Revisited .
4.2 Steerable Functions . . . . . . . . .
4.3 Steerable Dyadic Wavelet Transform . . . . .
4.4 Multiscale Spline Derivative-Based Transform . .


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


2.1 Transfer functions of direct B-spline 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


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 { 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 1-D 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 2-D discrete dyadic wavelet tranllsform. 35

4.2 Filter bank implementation of a spline derivative-based 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



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 B-spline. 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 one-dimensional discrete dyadic wavelet transform will serve

as a basis for deriving a steerable two-dimensional 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 spline-based approximation. The resulting algo-

rithm for computing a multiscale spline derivative-based 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 wavelet-based implementations and provides more flexibility

for different fusion and enhancement strategies by enabling orientation analysis.


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 rotation-invariant. By translation-invariant 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 translation-invariant by translations determined by the coarsest scale (e.g.. sixteen samples for
tlie analysis from Figure 1.1) [10].

transform that enables rotation-invariant processing. As an example of such a trans-

form, let us consider filtering with the first derivative of a two-dimensional Gaussian

probability density function in two directions, specifically, along x and along y-axis.

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 time-frequency analysis, (2) higher order derivatives of a Gaus-

sian can be. similar to the first derivative, used for rotation-invariant 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 scale-space "tracking" of emergent features.

1.2 Thesis Overview

We wish to construct an invertible, translation and rotation-invariant 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 translation-invariant 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 B-spline functions. In Chapter 2. we review

the basics of B-spline 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 scale-space filtering and reconstruction from edges.

In case of the first derivative of a Gaussian, a rotation-invariant transform

can be obtained by a tensor product extension of the one-dimensional 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 derivative-based transform, which in addition

enables rotation-invariant 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 mirror-extending 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,


L2(R) and L2(R2) denote the Hilbert spaces of measurable, square-integrable

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
(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 square-summable discrete signals f(n)
and f(t, ny). respectively.
The z-transform 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

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.


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 on the interval [.xO. .+i] is called a spline function of order p with the knot (i.e..
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 (B-splines). or more pre-
cisely,. cent ral B-splines having knolts at i E Z for p odd and at i + 1 for p even [34].
central l B-splines 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].
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 = P-1 (,"' + -) _P-1 ( X- 1) .

Another interesting property of B-splines is the fact that they converge to a
Gaussian as their order tends to infinity. Unser t al. [44] derived the Gaussian
-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 Battle-Lemari. wavelet bases
stem, can be built [10]. Here, we will not pursue constructions of orthogonal, semi-
orthogonal. or biorthogonal spline-wavelets any further-reader looking for a detailed
treatment of this subject may find a good starting point in books [8. 9].

'2.'- [B-Splin_,, Sit ,i;dl Thil,.r)l,,l iti'ii

Unser 0 al. developed a fast digital filtering scheme for B-spline signal pro-
cessing [43]. They defined a discrete B-spline 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


Interpolation of a discrete signal s(n) G 12(Z) by sp(x') G SP' using central


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 B-spline coefficients c(n) can

be obtained from samples s(n) as a "direct B-spline transform." Equation (2.2)

therefore represents "indirect B-spline transform" of a sequence {c(n)}.

After taking the z-transform of (2.2). the direct B-spline filters are found to be

p1;1(n) = 7-'{[Bp(z)]1}. Since B-spline functions ip(.r) are compactly supported,
indirect B-spline filters bp(n) are finite impulse response (FIR) filters, while direct

B-spline filters b-1(n) are infinite impulse response (IIR) filters. Aldroubi et al. [1]

showed that IIR filters b-1(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 B-splines 3p(x)

are symmetric.

Table 2.1 shows the z-transforms of direct B-spline filters for the first ten

orders. We postpone the discussion on implementation details of B-spline filters until

Chapter 5.

Table 2.1: Transfer functions of direct B-spline filters for orders from 0 to 9.
Z B# (z)

0 1

1 1

*) 8
_+6+ -1

3 6

4 384
2+76 -+230 -. -'t -

5 120
z2+26z+66+26z-1 +-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 B-spline 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 B-splines:

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.


68 -6

Figure 2.2: Fourier


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,


-8-. -6- -4,

4r 6,

-., -6. -

2.3 B i liin, q-iuji-i, \d | |l i i ,:..i i ,;fit I I -

Central B-splines are also simple to use when the goal is function approxima-

tion. Least-squares B-spline 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)

d(O) = (s(x),3,(x i)), (2.6)


W3(x) bp+= 1 + 3pi(x-

is the dual spline of order p [45]. Spline functions ,3(x.) and 3p(x) satisfy the biorthog-

onalitv condition

3p(xr M), 3p(x n) = 6(m n), ?, In e Z.

(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 L2-norm (i.e., least-squares) 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, shift-invariant 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
s,(x)= ,(i) 3(x- i) (2.7)


a(i) = 1 (s(,), 3,(x- i)).

where q-1 is the convolution inverse of the cross-correlation sequence

q12(i) = (3,(r -i),0-3(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 spaces--the 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 least-squares solution but

simpler to implement than the orthogonal projection (e.g., [47]).





$" 62 -6 6

(a) (14

Figure 2.3: Spline functions (a1) 3p(2) and (b) their respective Fourier transforms

,Jp( ) for {0.1,.2,3,4}.


In this chapter, we will augment the one-dimensional discrete dyadic wavelet

transform [28] to obtain a strong foundation for derivations of two-dimensional 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 translation-invariant, 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


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
^(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 2-1.

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


I ( ')2 t!(2' ) (2 ). (3.2)
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-


( + + {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

'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)

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(t-i) 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)

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.

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 one-dimensional 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


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)
while the one using the algorithm from [28] equals

-,-(n) = G- "(2) Im H-s(2'i). (3.14)
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).


1 5,



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
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(,) = 3-s (-ir+1 2 L+J Ccos (p)) (3.16)
m =1 /

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_____


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.01t2-1
2 O i ii.2,-',


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_____
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.i-2t. 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.I-i.1231 "- 1 7-"
5 -0.001953125 -0.0002441 t-i.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
(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).



the zero-crossings of I,n.(.r) for d = 2, a scale-space image of the representation sim-

ilar to [48] can be obtained. The only differences are that curves in the scale-space

plane are computed for dyadic scales only and that 0,,, is a close approximation of a

Gaussian instead of a true Gaussian.


In this chapter, we will extend the one-dimensional transform from Chapter

3 to two dimensions1 and derive several two-dimensional transforms with translation

and at least approximately rotation-invariant decomposition.

_1 D D1-,1.l.- D'.,-h- W ,.,L, I Tiaw-f-.ii R," i-it,.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) = 2-2,2 t,(2- mx. 2-y)

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
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, yi-i
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, 2-y).

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 two-dimensional analog of the one-dimensional

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 spline-based 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 one-dimensional 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 G-s(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 y) = I|(w .. )|2 o(22,). 2, )|2. (4.9)

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 one-dimensional case, a two-dimensional discrete dyadic wavelet

transform can be implemented as a fast hierarchical filtering scheme. To derive

such an implementation, we, similar to the one-dimensional case from Section 3.1,

use the definition of the two-dimensional 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


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
*p'- t.r) >3 g-s(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

\V|q---s(C2. .^y ) C"^q S(COr, ~y )B, -'(^.)B',7-1(C~y)L~p+r+1 (ol*) p++, +( 'y) (--s (Cr)

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" }14-s(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)

F{f 11 ,;'.(.r,.9) :, ,y= } (o;, oy) B, (B7 )(o 1)B ,7;(y B-)B,+,.+I(c,, )"
*Bp+,.+,(W-y) G_,(2-mwy) I H_,(2,w)H_,(2'1y). (4.13)

Equations (1.12) and (4.13) describe thlie decomposition part of the filter bank

implementation of a two-dimensional 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(/)

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


For d(= I il then follows that

lI',,, .x, y) = 2"'V(,s )(x.y). (1.14)

Thus for d(=2 we can write
Sll, .s(x y) = 22" A(, i), )(r. y). (4.15)

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 zero-crossings of (4.15) correspoiinds to the filtering carried out with a Marr-

Hildreth edge detector (Laplacian of Gaussian) [29]. (Note that both edge detectors




Hs(2N) J(2o )


Ks(o,) T,(o~y)

Br(Ox) Br,((Oy)

- H_(2(o ) H_(2(y)
Figure 4.1: Filter bank implementation of a two-dimensional discrete dyadic wavelet
transform (a) decomposition and (b) reconstruction for two levels of analysis. H',(.)
denotes the complex conjugate of H-L(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 zero-crossings of Z, 12 I. ,s(.r. y) is therefore an approximation to the (C'annv or

the NI rr-Hildreth 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 two-dimniensional 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
f(r.0 0o) = ci(0).f(r,, 0). (4.16)
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 B-splines can be used, as we saw in Section 4-1. to approximate
Laplacian of a Gaussian for approximately rotation-invariant, 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:
f(r. 0)= a a (r)d*" (.1.17)

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,(1.) -,, oo = i(0o)a,,(,/') (-7 i' ( .18)

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 one-dimensional transform

from Chapter 3. Steerable transforms are nothing new-quite 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 B-splines which, as the order of B-splines 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 one-dimensional 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


.(,.)= E EyIF 0(y). (-1.23)
m =-x. i=1
where \ (xr.y) denotes \m(0'x. y) rotated by Oi and \,,, (x. y,) = 2-2" \(2-m., 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 one-dimensional wavelets

from Section 3.1:`

/sin( P+-"+'
( ,... ) = 0(ju. cos(.'))d (sin( ) (1.25)

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).


o(2w,.) = Ht,(w,),)(.,), (4.26)

i"(<.,.o) = s(L,,, ;o ) o(,,.), (4.27)

,(2)) = t,,( r) (a,,.). (4.28)


'("-,..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
C,,(^,.. 9,s,- 0)A, ,, ,,- 0O) + Ht(.,)Lt(.,) = 1. (4.30)
Setting o(C',.) = I,(-,.), and employing (3.5) and (4.25) with (4.26) and (4.27).

we find that

t ( ,.) = /- (-,,-) (1.31)


,( ,., 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


tL r) =t(,.) (4.33)


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 -l-B ,- i T ,li i i..1 11i

Let us pause here for a brief assessment of tie two-dimensional steerable trans-

forml derived so far. \Ve have chosen steerable wavelets (4.25) which are equal to d-th

order derivatives of circularly symmetric spline functions in the direction of .r-axis

(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 x-y 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 B-splines.

Let us approximate wavelets from (4.25) with
d/ +^W^L) / /.
u(x -y) = '+-) 3p+d(Y)- (4.35)

Based oni the fact that B-splines 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., d-th order derivatives of the normalized Gaussian in the direction of .r-axis) 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 B-spline 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) fron-m (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.)

where i = (cos 0o. sin 0o) = (1s., IV).
Note that in case of Gaussian, which is both ,r-y 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, x-y 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+d-i( '.'/)" /' = 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 ( )3-p-2(,y )

KS (.r) Ks(-y),P-i.'-}p-(\

T-2(,,,)I\_2 (^}3-2(^)p^

3 (",'.&-y)

\ (t... y






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(^~),






\3(w) = (1 H(,2),


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


-(;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 )





wxv ere

14( .)= i- |I(.)I.


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 onie-dimiensional transform from Chapter 3 as much as possible. All

decoimposit ion filters (;d(,:) were first paired with corresponding reconstruction filters

41 4-32,,x -2,
\4(ua,',:. ,.

( '-.-^..v,) = 7: (+,,.)1,;( .y)) _=,(*J,.')*( ).

Kd(,). and then all other potential digital filters were specified as polynomials in

H-s(-:). 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 higher-degree 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+d-i( ty any) ddt dty,

where i = 0, 1 ..... d as in (4.37).

Again, we approximate since functions with r-order cardinal splines, then use

(2.4), and get

D FT r-f U) s'=y-(
17~.S,, Y } I -Y IY (Ux' C,,y) B -I (,) B -l( ,,y) Bl,+,+i+l(,')

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 derivative-based

transform for d = {1. 2.3,4}. For d = 1, we recognize (except for the prefiltering

and postfiltering) the filter bank implementation of a two-dimensional 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 y-axis 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 ones-a possible

advantage for certain applications. They have, on the other hand, little control over

the function from which derivatives are computed (to obtain a d-th 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 rx-y separable.

Br'(Cx)B;Coy) Bp+r+i(Owx) Bp+,,1(Co)

BR|1. )J Bp+r+i((Oy)


Hs(2m),) Hs(2m(Oy)


G.,(2mx)J GCs(2 (y)


Hs(2eo) H(2my)

Br((o),) Br(o)y)

Ks(2mmo) T(2m(oy)

Tl(2m%,,) K~s(mwy)

Ls(2mox) Ls(2moy)

K(2mo),) T(2mOy)

Kl(2 mox) Ks(2moy)

T2(2m(o) K(2 %))

Ls(2mx) Ls(2mwy)

Figure 4.2: Filter bank implementation of a multiscale spline derivative-based 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.

B-1 )

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.


Except for the steerable dyadic wavelet transform, all transforms presented in

Chapters 3 and 4 can be implemented as filter banks consisting of one-dimensional 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 two-dimensional 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 one-dimensional 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 even-length 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

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

Now. we are ready to examine the filter bank implementation of a one-dimensional
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
F!,,,u,(n) = f(0)uj,(n)+-y .f(')[Itjj(ii-2'"'i)+u11(i+2"'i)]. nC e [0. A--1]. (5.3)

(-- ) 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
}- -
(-,.OU(n) = g(i)[uII(n -i 1) u11(n (+i)]. + E [1, A- 1].

for d odd. (5.3) for d(I even.

H_,,oii(n) = h(i)[uII(n I 1) + uii(rn +i)]. EC [O. A].

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
/-- -
G,, u(n) ="( (i)[u(n- 2m(i +.s))- u((n, + 2"(i + s))], n E [1. N'- I],

for d odd and p even.
(G_.,,,u(n) = g(i)[un( 2"..(I + s)) uii(n + 2"(i + s))]. n C [0. A' 1].

for d and p odd.
F-.1u(n) = f(O)I((n) + > f(i)[uI(u 2"t1) + it(n + 2l)]. ( + [0. A], (5.5)

with /(n) = g(nt) for d and p even,
--s.mu(n) = h(i/)[u/i(n 2"1(i + .s)) + utI(n + 2m(i + .))], e [0. ]. (5.6)

for p even. and (5.3) for p odd, where

1 1(-n) ift E [ -. -1]
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
Kx.ou(n) = + k(i)[1n)( -i+)- utiI(n+ i)]. + C [0, N 1].

for d odd, (5.3) for d even,
L,,oa(,) = l(i)[ui(n i + 1) + u(n +i)]. n E [0, A 1].

for p even. and (5.3) for p odd, where

S-u(-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
A'2.,,1u() = k(i + 1)[aI(r? 2'(i + s)) -ilii(n + 2' (1 .s))]. a E [0, N'].

for d odd and )p even. (5.5) with f(ii) = k'(n) for d and ) even.

', ..,,,(,) = k A.(i + l)[u,(, 2"(i + .)) n -(- +- 2'(i 2+ .))]. c [0. A I].

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).

i-t(--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 one-dimensional case. whose filter

bank implementation is depicted in Figure 4.2. Two-dimensional transform filter

bank implementations (Figures 4.1 and 4.2) are not only comprised of .rx-y 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 two-dimensional 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 zero-phase 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


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 b-1(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 B-1(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- az-1 +l-acz

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 a-2
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.


c(n) = (c(;? + 1)- c+(n)) i = X-2, N-3.....0.


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(AV-1).

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.X--i 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)

N-1 \ N+ i
00+0 +I +
c(A -1) 1 (C+ -1) + 0 + (i))
-- (+(-+ I -i,(I)).
1= --- ^


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- n-1) 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

C,(0) = 2[ -1 111p()],, -

-^EE^ ^]t + [C2I^1}(
1 =(0) + E J uQ)
11=0 i=O


a |11(1

c(X-1) 0Nx1
C(A- 1) -- (C+ (N -) E ]+o u( (5.16)
o n=O i=0

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] 2-V-__,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 b-l(n) from Figures 1.2(h) and

4.2(j) with tin E [1, 31 1] and p even are
C+ (0) =^ [-Up?]
[o=--x,/ -

y a ,,, (0)+ o ,'--- "(k')+
1n=O'O +

+ 2 ] + [ u(l) (5.17)

C(N) (2c+(N) -u(N)).
1 a2

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.


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 wavelet-based 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 derivative-based trans-


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 derivative-based 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. Low-level 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 rotation-noninvariant

wavelet techniques than for the multiscale spline derivative-based 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


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 derivative-based transform.


modify wavelet transform coefficients, and Figure 6.2(c) shows the result of enhance-

inent using the multiscale spline derivative-based 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 rotation-invariant 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 rotation-invariant 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 derivative-based transform.



Figure 6.3: The magnitude of filter coefficients for
(three levels of analysis shown).

(;0( w.^/) at oi

{O. 7. 3}


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


We extended the one-dimnensional discrete dyadic wavelet transform to higher-

order derivatives and even-order 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 two-dimensional transforms. All of them were derived

with the goal of eliminating artifacts due to lack of translation and rotation in-

variance. \We presented a two-dimensional discrete dyadic wavelet transform with a

first-derivative wavelet as an extension of the one originally proposed in [28], a two-

dimensional discrete dyadic wavelet transform with a second-derivative wavelet that

can approximate Laplacian of a Gaussian, and a steerable dyadic wavelet transform

implemented in a near-perfect reconstruction filter bank which may be preferred when

there is a need for orientation processing along equally spaced angles. WNe derived

a multiscale splinie derivative-based transform which uses x-y separable filters in a

perfect reconstruction filter bank and enables fast translation and rotation-invariant

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.


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.


[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 scale-space filtering," IEEE Trans. Pathtrn A1nal. IMach.
Inh1Il.. vol. 8. no. 1. pp. 26-33, 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. 173-182.
[6] J. ('annv. "A computational approach to edge detection," IEEE Trans. Pattern
Anal. Mach. Intell.. vol. 8. no. 6, pp. 679-698, 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.
[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. 163--189.
[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. Hel-Or 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. Kronland-Martinet. J. Morlet, and Ph. Tchanmitchian, "A
real-time algorithm for signal analysis with the help of the wavelet transform." in
lai'arlts. 'ime-Frequei i' / Mthods and Phas( Spac(. J. M. Comibes, A. Gross-
mann. and Ph. Tchamitchian. Eds., Springer-Verlag. Berlin. GermanY. 1989. pp.
[16] A. J.. Jerri. "The Shannon sampling theorem-its 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 Timnr-Frcqucncy 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 3-d 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:32-2: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. 736-749.
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~ i-sensor image fusion using the
wavelet transform," in Proc. IEEE Int. ('onf. Image Process., Austin. TX, 1994,
vol. 1, pp. 51-55.
[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. 710-732. 1992.

[29] D. Marr and E. Hildreth, "Theory of edge detection." Proc R. .Soc. London
Ser. B. vol. 207. no. 1167. pp. 187-217, 1980.
ill] A. V. Oppenheim and R. W. Schafer, Discrete-Tin. 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. 569-586. 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. 112-141.
[35] I. J. Schoenberg. "C(ardinal interpolation and spline functions." .J. Appro.r.
Thory, vol. 2. no. 2, pp. 167-206, 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. 87-93, 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 multi-scale derivative computation." in Proc. IEEE Int. ('onf. Imnag
Process.. Washington, D.C., 1995, vol. 3. pp. 444-447.
[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. 587-607,
[40] M. .1. T. Smith and T. P. Barnwell, III, "Exact reconstruction techniques for
tree-structured subband coders." IEEE Trans. Acoust. cSpech S.5ignal Process..
vol. 34. no. 3 pp. 434-441, 1986.
[41] A. Toet. "Image fusion by a ratio of low-pass pyramid," Pattern Rcognlit. Lett.,
vol. 9. no. 4. pp. 245-25: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 B-spline transforms for continuous
image representation and interpolation," IEEE Trans. Pattern Anal. Mach.
Intll., vol. 13. no. 3. pp. 277-2-'., 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. -%, 1-872. 1'l1i1'.


[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. 95-103. 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. 3519-3523. 199-1.
[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. 1459-1462.
[48] A. Witkin. "'Scale space filtering," in Proc. Int. Joint Conf. Artif. Inf l.,
Karlsruhe. Germany. 1983, pp. 1019-1022.


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


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

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


3 1262 08555 0696

Full Text
xml version 1.0 standalone yes
PreviousPageID P110