Wavelet based multiresolution expectation maximization image reconstruction algorithm for positron emission tomography

Computerized Medical Imaging and Graphics PERGAMON

Computerized Medical Imaging and Graphics 24 (2000) 359–376 www.elsevier.com/locate/compmedimag

Wavelet based multiresolution expectation maximization image reconstruction algorithm for positron emission tomography A. Raheja a, A.P. Dhawan b,* a

Department of Computer Science, Philadelphia University, Philadelphia, PA 19144, USA Department of Electrical & Computer Engineering, New Jersey Institute of Technology, Newark, NJ 07102, USA

b

Received 16 December 1999; accepted 23 May 2000

Abstract Maximum Likelihood (ML) estimation based Expectation Maximization (EM) [IEEE Trans Med Imag, MI-1 (2) (1982) 113] reconstruction algorithm has shown to provide good quality reconstruction for positron emission tomography (PET). Our previous work [IEEE Trans Med Imag, 7(4) (1988) 273; Proc IEEE EMBS Conf, 20(2/6) (1998) 759] introduced the multigrid (MG) and multiresolution (MR) concept for PET image reconstruction using EM. This work transforms the MGEM and MREM algorithm to a Wavelet based Multiresolution EM (WMREM) algorithm by extending the concept of switching resolutions in both image and data spaces. The MR data space is generated by performing a 2D-wavelet transform on the acquired tube data that is used to reconstruct images at different spatial resolutions. Wavelet transform is used for MR reconstruction as well as adapted in the criterion for switching resolution levels. The advantage of the wavelet transform is that it provides very good frequency and spatial (time) localization and allows the use of these coarse resolution data spaces in the EM estimation process. The MR algorithm recovers low-frequency components of the reconstructed image at coarser resolutions in fewer iterations, reducing the number of iterations required at finer resolution to recover high-frequency components. This paper also presents the design of customized biorthogonal wavelet filters using the lifting method that are used for data decomposition and image reconstruction and compares them to other commonly known wavelets. 䉷 2000 Elsevier Science Ltd. All rights reserved. Keywords: Multiresolution reconstruction; Wavelets; Expectation maximization; Positron emission tomography; Lifting scheme

1. Introduction Expectation Maximization (EM) algorithm has been used for Positron Emission Tomography (PET) image reconstruction with several variants [4–6]. Recently, wavelet processing has been used [7–12] to deduce a multiresolution (MR) representation to solve inverse problems similar to this problem of image reconstruction. Use of wavelets for MR tomographic reconstruction using the Filtered Backprojection operator has been investigated [9–12]. In image analysis and reconstruction approaches, MR methods have produced better results than conventional approaches [13]. A similar approach for representation of the EM operator in the wavelet domain is not so trivial because of the complex nature of the operator and the sparse nature of the probability matrix (weight matrix which entirely depends on the geometry of the imaging system). We have used standard EM estimation method [1] for MR reconstruction using wavelet analysis to improve efficiency * Corresponding author. Tel.: ⫹1-419-530-8267; fax: ⫹1-419-530-7392. E-mail address: [email protected] (A.P. Dhawan).

and reconstruction quality. The wavelets are used to construct a MR data space, which are then used in the estimation process. The beauty of the wavelet transform to provide localized frequency-space representation of the data allows us to perform the estimation using these decomposed components. The advantage of this method lies with the fact that the noise in the acquired data becomes localized in the high–high or diagonal frequency bands and thus not using these bands for estimation at coarser resolution helps speed up the recovery of various frequency components with reduced noise estimation. Hence, it is very important to have a wavelet basis that allows for good spatial frequency localization into the various subbands at different levels. To ensure the above, biorthogonal wavelet filters were custom-designed using the lifting scheme [14–16] that provides a robust method of creating new biorthogonal wavelets from existing biorthogonal wavelets. The new wavelet filters have a frequency response that ideally localizes the various frequency components of the tube data for reconstruction purposes. A new and robust wavelet based stopping criterion and a wavelet spline based interpolation [3] method used for grid

0895-6111/00/$ - see front matter 䉷 2000 Elsevier Science Ltd. All rights reserved. PII: S0895-611 1(00)00035-5

360

A. Raheja, A.P. Dhawan / Computerized Medical Imaging and Graphics 24 (2000) 359–376

mapping is presented for this iterative algorithm. This approach can be generalized to other image reconstruction methods that use the concept of EM. The main goal of this paper is to introduce the MR concept in the EM framework and to demonstrate the improved reconstruction by generating a wavelet basis appropriate for frequency localization of tomographic PET data. The rest of the paper is organized as follows. Section 2 reviews the wavelet theory in general and the design of biorthogonal wavelets using the method of lifting. The wavelet based MR approach to EM is described in Section 3. Section 4 gives a brief overview of the design of the custom wavelet filters for the problem of tomographic image reconstruction and compares the new wavelet filters with commonly known wavelet filters. Section 5 discusses the results of the reconstruction using phantom PET data and compares the effect of various wavelet bases in reconstruction. Finally, Section 6 presents conclusions and discusses the future direction of this work. 2. Wavelets and MR analysis

where Wj is a subspace spanned by c (2 jt ⫺ k) for all integers j, k 僆 Z. This basic requirement of a MR analysis is satisfied by nesting of the spanned subspaces similar to scaling functions i.e. c (t) can be expressed as a weighted sum of the shifted c (2t) as p X c…t† ˆ 2 g n f…2t ⫺ n† …6† n

p where gn is a set of filter coefficients and 2 maintains the norm of the wavelet function with the scale of two. The wavelet-spanned subspace is such that it satisfies the relation Vm⫹1 ˆ Vm 丣 Wm where L2 ˆ … 丣 W⫺2 丣 W⫺1 丣 W0 丣 W1 丣 W2 丣 …

(7)

Hence, the wavelets are functions whose dilations and translations form an orthonormal basis of L 2(R). Since, the wavelets span the “difference” or orthogonal complement spaces …W⫺1 ⬜ W⫺2 ⬜ W⫺3 ⬜ V⫺1 †; the orthogonality requires the scaling and wavelet filter coefficients to be related through the following:

2.1. Orthogonal and biorthogonal wavelets

gn ˆ …⫺1†n h1⫺n

Assuming most readers to be familiar with wavelets, only a brief introduction to the concept of discrete wavelet transform is presented here. Details on wavelets can be looked up in Refs. [17–19]. A scaling function f (t) in time t, also known as the “father wavelet,” can be defined as

In the case of orthogonal wavelets as described above, hn and gn are called a pair of quadrature-mirror low-pass and high-pass filters. The biorthogonal wavelet system generalizes the classical orthogonal case. For biorthogonal wavelets, a dual scaling function is similar to the MR approach of the scaling function is defined as p X f~ …t† ˆ 2 h~ n f~ …2t ⫺ n† …9†

fj;k …t† ˆ 2j=2 f…2j t ⫺ k†

…1†

where j is a scaling parameter and k a translation parameter and j; k 僆 Z; set of all integers. The scaling and translation of which generates a family of functions that span by using the following “dilation” equations: p X f…t† ˆ 2 hn f…2t ⫺ n† …2† n

p where hn is a set of filter coefficients and 2 maintains the norm of the wavelet function with the scale of two. To induce a MR analysis of L 2(R), where R is the space of all real numbers, it is required to have a nested chain of closed subspaces defined as … 傺 V⫺1 傺 V0 傺 V1 傺 V2 傺 … 傺 L2

…3†

Specifically, there exists a function c (t) in time t, the “mother wavelet,” such that

c j;k …t† ˆ 2j=2 c…2 j t ⫺ k†

…4†

form an orthonormal basis of L 2(R). The wavelet basis induces an orthogonal decomposition of L 2(R), i.e. one can write … 傺 W⫺1 傺 W0 傺 W1 傺 W2 傺 … 傺 L2

…5†

…8†

n

Similarly, one can define the dual wavelet as p X c~ …t† ˆ 2 g~ n f~ …2t ⫺ n†

…10†

n

The following relations can now relate the dual scaling and wavelet filters: gn ˆ …⫺1†n h~ 1⫺n

…11†

g~ n ˆ …⫺1†n h1⫺n

…12†

The subspaces spanned by the scaling functions, dual scaling functions and wavelet and dual wavelet functions follow the relationship Vj ⬜ W~ j

…13†

V~ j ⬜ Wj

…14†

2.2. Wavelet expansion and reconstruction of 1D and 2D signals Mallat [20] has showed using a pyramid algorithm how orthogonal wavelets could be used to create an MR

A. Raheja, A.P. Dhawan / Computerized Medical Imaging and Graphics 24 (2000) 359–376

(a)

2

h

h

2

aj-1

g

2

bj-1

aj

aj+1 g

361

2

bj

(b)

2

aj-1

~ h

aj 2

bj-1

2

~ h

2

g~

g~

aj+1 bj

Fig. 1. (a) Two-stage two-band analysis tree. h and g are the low-pass and high-pass decomposition wavelet filters as mentioned in the previous section. (b) Two-stage two-band synthesis tree. h˜ and g˜ are the low-pass and high-pass reconstruction wavelet filters as mentioned in the previous section.

g

1 2

3

Dj

2 1

g

j;k

h

1 2

Dj2

g

1 2

Dj

Aj+1

1

2 1

h

where the brackets denote the inner product. In the MR analysis, the function is successively approximated in the subspace Vj. From the representation in Fig. 1(a), one can write aj⫹1 …t† ˆ

h

1 2

Aj

(a) 2j0

representation of a signal and Unser et al. [21] extended this to the case of non-orthogonal wavelets using splines. Using the definitions in the section above, any signal aj⫹1 …t† 僆 L2 …R† can be written as E XD aj⫹1 …t†; c j;k …t† c j;k …t† …15† aj⫹1 …t† ˆ

J-1

X

aj0 …k†fj0 ;k …t† ⫹

⫺1 X JX

k

Fig. 2. (a) Single level decomposition of an image Aj⫹1. g and h are the highand low-pass filters, 2 # 1: keep one column out of two, and 1 # 2: keep one row out of two. (b) Schematic representation of the 2D-wavelet transform coefficients of an image. Aj ; refers to the approximation at resolution j while D1j ; D2j ; and D3j ; refer to the details at resolution j.

bj …k†c j;k …t†

…16†

jˆj0

where j0 ˆ 2 in Fig. 2(b). Similar to the analysis, the reverse process of synthesis to get a perfect reconstructed signal is seen in Fig. 1(b). The synthesis of any wavelet-decomposed signal should give a perfectly reconstructed signal if orthogonal or biorthogonal wavelet filters are used. In case of a 2D signal, the same concept discussed above can be extended to L 2(R 2). The wavelet MR analysis can be constructed using a separable scaling function that can be written as

f…x; y† ˆ f…x†f…y†

(b)

k

…17†

This implies that the wavelet transform can be performed as described above for the 1D case to the 2D signal. This is achieved by first applying the transform to the rows of the image and then to the columns of the image. 2.3. Multiresolution analysis Mallat [20] showed how separable functions can be used

362

A. Raheja, A.P. Dhawan / Computerized Medical Imaging and Graphics 24 (2000) 359–376

to create a MR representation of a 2D signal. In brief, a onestage 2D wavelet transform can be performed as seen in Fig. 2(a). Performing the wavelet transform on an image result in four sub-images, each one of them is half the size of the original because of the downsampling of two for both rows and columns. These subimages Aj, D1j , D2j and D3j are referred to as the low–low (LL), low–high (LH), high– low (HL) and high–high (HH) subbands, respectively. Unlike a complete analysis, where each subband is further decomposed to the desired level, a wavelet tree analysis decomposes only the LL band further to the desired level. Synthesis for an analyzed image is done similar to that of a signal in Fig. 1(b). 2.4. Design of biorthogonal wavelet filters using lifting The idea of lifting as a method of biorthogonal wavelet construction was introduced by Wim Sweldens [15]. Lifting is also a simple method that is used to increase the number of vanishing moments or dual vanishing moments of a wavelet or dual wavelet. This section defines and discusses briefly the design considerations while using the method of lifting. The reader is pointed towards Sweldens’s work [14– 16,22] for an in depth understanding of the lifting scheme from which this entire section has been derived. The lifting scheme was inspired by two theorems, one from the work of Cohen et al. [23] and the other from Chui [24]. Sweldens’s corollary [16] uses these two theorems to define the following. Take an initial set of finite biorthogonal filters ~ Then a new set of finite biorthogonal filters {h; h~0 ; g0 ; g}: ~ g; g} ~ can be found as {h; h; ~ v† ˆ h~ 0 …v† ⫹ g… ~ v†s…2v† h…

…18†

g…v† ˆ g0 …v† ⫺ h…v†s…2v†

k

X

h~ 0k f…2x ⫺ k† ⫹

k

c~ …x† ˆ 2

X

c~ …x† ˆ f~ …2x ⫺ 1† ⫺ af~ …x† ⫺ bf~ …x ⫺ 1†

X

s⫺k c~ …x ⫺ k†

…21†

…23†

The coefficients or lifting parameters a and b are chosen such that the resulting dual wavelet is symmetric and has a vanishing moment. sym_ab implies the new wavelet to be symmetric i.e. sym_ab ˆ a ˆ b

…24†

The pth vanishing moment, denoted by momp is defined as Z …25† momp ˆ x p c~ …x† ˆ 0 mom0 for Eq. (23) lead to the symmetry condition. Using condition mom1, the result is given as Z1 Z1 Z1 c~ …x† dx ˆ w~ …2x ⫺ 1† dx ⫺ a w~ …x† dx ⫺ 1=2

⫺ 1=2

⫺b

…19†

where s(2v ) is a trigonometric polynomial. As seen from the above definition, the scaling filter and the dual wavelet filter remain unchanged. The lifting scheme theorem is defined in terms of biorthogonal functions as follows. Take an initial set of finite biorthogonal scaling functions and wavelets {f; f~ 0 ; c 0 ; c~ 0 }: Then a new set {f; f~ ; c; c~ }; which is formally biorthogonal can be found as: X c…x† ˆ c 0 …x† ⫺ sk f…x ⫺ k† …20†

f~ …x† ˆ 2

The sequence s provides the ability to manipulate the wavelets and dual functions that can be constructed from a simple scaling function. This sequence can also be used to provide a wavelet with increased number of vanishing moments or have a desired shape. Sweldens’s work [14] provides examples using the finite impulse response (FIR) filter coefficients based upon the work of Deslauriers and Dubuc’s [25] as the initial set of biorthogonal filters that are derived from interpolating scaling functions. Sweldens’s definition [16] of an interpolating scaling function implied that every dual interpolating scaling function has as its dual, a Dirac function as its scaling function. This concept is combined with the lifting method to generate new wavelets with higher vanishing moments. In Ref. [14], the linear spline hat functions are used with two lifting parameters to create a new dual wavelet and a scaling function. The equation for the dual wavelet can be written as

⫺ 1=2

Z1 ⫺ 1=2

w~ …x ⫺ 1† dx

ˆ0

(26)

Solving this yields a ˆ 1=2 ⫺ b

…27†

Using Eqs. (24) and (27), the lifting coefficients are a ˆ b ˆ 1=4: These coefficients lead to the following representation of the combination of old wavelet with two scaling functions (reduced by a factor of 1/4 each) at the same level to form the new dual wavelet as seen in Fig. 3 [16]. This technique of finding the lifting parameters and then combining old wavelets with scaling functions from the same level results in design of new wavelets.

k

g~ k f~ …2x ⫺ k†

k

where the coefficients sk can be freely chosen.

…22†

3. Method: wavelet based MREM The MREM reconstruction method is an extension of the multigrid EM (MGEM) reconstruction concept introduced

A. Raheja, A.P. Dhawan / Computerized Medical Imaging and Graphics 24 (2000) 359–376

(a)

363

(b)

(c)

Combination

(d)

(e)

Fig. 3. (a) Piecewise linear spline hat dual scaling functions at level j. (b) Old dual wavelet function at level j. (c) Dual scaling functions at neighboring level j. (d) Combination of old dual wavelets and one-fourth times dual scaling function at the same level. (e) The new dual wavelet created as a result of lifting.

by Dhawan et al. [2] but merged with the MR approach of wavelets. A schematic approach to the wavelet based MREM method is presented and compared with the MG approach in Fig. 4. In the MGEM algorithm [2], a set of grids G k, k ˆ 1; …; K represented the same reconstruction space at different resolutions. G 0 was the coarsest grid while G K was the finest resolution grid. The reconstruction pyramid consisted of square grids of sizes S k, k ˆ 1; …; K such that the ratio, S k⫹1 : Sk ˆ 1 : 2: The maximum frequency, f k, of the image on grid G k is given by 1/2S k, which is twice the maximum frequency that can be represented by G k⫺1. Using the hierarchical structure of the reconstruction pyramid, the low-frequency components of the finest resolution were recovered quickly at the coarser grids. In the MREM algorithm [3], the authors introduced the concept of varying resolution in the detector space. This was done by combining groups of adjacent detectors i.e. rebinning the tube data to create new data spaces at coarser

resolutions. This process of rebinning increases the signal to noise ratio of the data at coarser resolutions, but using these new data spaces in the MR framework provided faster and better quality image reconstruction. In this work, wavelets are used to create MR data spaces in order to use the lower resolution data spaces without using much of the high-frequency components of the original tube data. The new Wavelet based MREM (WMREM) reconstruction algorithm aims at using the relatively noise free tube data at coarser resolutions to recover the low-frequency components of the image and recover most of the high-frequency components at the original resolution of the data. Thus, this algorithm merges the EM algorithm with the wavelet based MR approach, and also utilizes spatial-frequency localization analysis obtained from the wavelet processing in the image reconstruction. To optimize spatial-frequency localization analysis for better reconstruction, the wavelet is custom designed using the method of lifting [14–16].

364

A. Raheja, A.P. Dhawan / Computerized Medical Imaging and Graphics 24 (2000) 359–376 Tube data Wavelet Transform and Analysis p(b,d) weight matrix

Tube data at p(b,d) coarser resolutions weight matrix

Tube data

Wavelet Transform: Transition Criterion and Interpolation

Multigrid Estimation Process (MGEM)

Multiresolution Estimation Process (MREM)

Final Reconstructed Image = Final Estimations of the MGEM algorithm

Reconstructed Image = Inverse Wavelet Transform of the final estimations at various resolutions

(a)

(b)

Fig. 4. The difference of the concept between the (a) multigrid EM (MGEM) and the (b) wavelet based multiresolution EM (WMREM).

Level 0

Level 1

Level 2

2

2

2

h

n0*(d) Tube data 128x128

h

n1*(d)

32x32

2

g

n2*(d)

32x32

2

h

n3*(d)

32x32

2

g

n4*(d)

32x32

h

2

2

2

h

g

2

g

n5*(d)

64x64

2

h

n6*(d)

64x64

2

g

n7*(d)

64x64

g

(a)

n*0 Original Tube data

Multiresolution Representation

(128x128)

n*1 32x32

n*2 32x32

n*3 32x32

n*4 32x32

n*5 64x64

n*6

n*7

64x64

64x64

(b) Fig. 5. (a) Two-channel two-stage analysis of the tube data n0ⴱ(d) at the finest resolution using the wavelet low- and high-pass decomposition filters h and g, respectively. (b) MR representation of the data spaces.

A. Raheja, A.P. Dhawan / Computerized Medical Imaging and Graphics 24 (2000) 359–376

n1*(d)

32x32

EM

2

~ h

n2*(d)

32x32

EM

2

g~

n3*(d)

32x32

EM

2

~ h

365

Wavelet reconstruction

Initial estimation

n5*(d)

64x64

EM

2

g~

n6*(d)

64x64

EM

2

~ h

Reconstructed Image

Fig. 6. MR reconstruction scheme for WMREM. EM is the Expectation Maximization process. h˜ and g˜ are the low-pass and high-pass reconstruction filters.

Create data spaces DWT2( n *0 (d ) ) results in

n 1* , n *2 , n *3 , n *4 , n *5 , n *6 , n *7 Initialize Grid level, k = 0 (coarsest) Data space, j = 2 (coarsest) Iterations, i = 0

λ

0 0

Initial Solution * * * * * for n 1 , n 2 , n 3 , n 5 , n 6

λ j , k = EM( λ j , k ,n 1) i i+1 λ j , k = EM( λ j , k ,n 2) i i+1 λ j , k = EM( λ j , k ,n 3) i+1

i

1

1

2

2

3

3

i = i +1

Is intra-grid level performance measure satisfied ?

NO

Is current grid resolution > detector resolution?

YES

YES

Wavelet Decomposition

NO

Is grid optimization measure satisfied ?

Final Reconstructed Image

Final Reconstructed Image

YES

Is J=0 K =2 NO 0

YES

λ j ,k+1 = DWT -1(λj ,k, λ ji ,k , λ ji ,k)

λ j ,k+1 = DWT -1(λ j ,k , λ j ,k , λ j ,k )

k=k+1 j = j -1

k = k+1

i

4

1

2

Wavelet Synthesis

3

0

8

i

i

i

4

5

6

Wavelet Synthesis

i=0 Fig. 7. Flowchart for the WMREM algorithm. DWTn refers to a forward discrete wavelet transform and n refers to the level. When n ⬎ 0; it is a forward discrete wavelet transform and inverse wavelet transform when n ⬍ 0:

366

A. Raheja, A.P. Dhawan / Computerized Medical Imaging and Graphics 24 (2000) 359–376

Fig. 8. Shepp and Logan Phantom.

3.1. MR data spaces The MR data spaces can be constructed using the original tube data and performing a two-stage or two-level two-band analysis as seen in the Fig. 5(a). The reason for choosing two levels of decomposition is because it provides optimal data spaces in terms of information contained in them. The selection of two levels is discussed in detail under the results in a latter section. Using the notation introduced earlier in Section 2 for the separable wavelet transform of an image, we can write the following: nⴱ1 …d† ˆ A2 nⴱ0 …d†

…28†

nⴱ2 …d† ˆ D12 nⴱ0 …d†; nⴱ3 …d† ˆ D22 nⴱ0 …d†; nⴱ4 …d† ˆ D32 nⴱ0 …d† …29† nⴱ5 …d† ˆ D11 nⴱ0 …d†; nⴱ6 …d† ˆ D21 nⴱ0 …d†; nⴱ7 …d† ˆ D31 nⴱ0 …d† …30† Eq. (28) is the low-pass component of the tube data after two levels of decomposition also called the LLLL component. The three detail components at level 2 in Eq. (29) are usually referred to as the LLLH, LLHL and LLHH components, respectively. The three detail components at level 1 in Eq. (30) are referred as the LH, HL and HH components, respectively. The image representation in Fig. 5(b) displays the creation of 7 new data spaces from the original data using a two-level wavelet-decomposition. The data spaces nⴱ4 …d† and nⴱ7 …d† shall not be used in the reconstruction algorithm since they contain the possible high-frequency noise of the data collection process. The coarser resolution data created in this manner is relatively free of the high frequency of components of the original data. 3.2. WMREM algorithm In the WMREM method, the MREM is performed within

the framework of EM using the newly constructed data spaces and switching grid and data space resolution simultaneously to perform a faster and better quality reconstruction. The interesting property of spatial localization of the data after performing wavelet decomposition allows for the use of these data spaces. The LL data space nⴱ1 …d† is simply a low-pass filtered image and can be used without any modification. The HL and LH components i.e. nⴱ2 …d†; nⴱ3 …d†; nⴱ5 …d† and nⴱ6 …d† are subjected to thresholding or DC shifting to be used as data spaces for the estimation process. Both these methods, i.e. thresholding and DC shifting are evaluated in the framework of this algorithm. The WMREM reconstruction process using wavelet basis is shown in Fig. 6. It is interesting to note that the interpolation step (for mapping from a coarser grid level to a finer one) performed in the MREM [3] using a wavelet spline based method is now naturally performed by the wavelet reconstruction scheme until the finest data resolution is reached. The WMREM algorithm, however, does retain the concept of performing EM at one resolution finer than the desired grid level i.e. it performs EM at the desired resolution and one grid finer using the finest data space which is the original data itself. The mapping from the desired grid level to the next finer level is done using the wavelet spline method [3]. Also, the stopping criterion is similar to the one introduced for MREM. A flowchart of this algorithm is shown in Fig. 7. The algorithm can be described briefly as follows: Step 1. Create new data spaces using a separable 2D wavelet transform on the data. Step 2. Threshold/DC shift the data spaces that involve the LH or HL filtering process. Step 3. Initialize the starting solution to a constant value for the LL, LH and HL spaces at the coarsest level, i.e. say level ˆ J and the HL and LH components at other levels i.e. J ⫺ 1; …; 1: In this algorithm, J ˆ 2: Coarsest grid level starts from K ˆ 0 which corresponds to coarsest data level at J ˆ 2: Step 4. Perform EM using these 3 data spaces at level J till the stopping criterion is satisfied. Grid used has the same dimension as the data spaces. Step 5. Perform inverse wavelet transform using the three images as a result of step 4. The HH image is taken as zeros for this wavelet synthesis step. Step 6. The reconstructed image is used as the starting solution for the LL component at the next finer resolution, i.e. J ⫺ 1: The LH and HL components use a constant value image as the starting solution similar to Step 3. Grid resolution is the same as that of the data space, i.e. K ˆ 1 at data level J ˆ 1: Step 7. Steps 4, 5 and 6 are performed on these data spaces until the coarsest resolution of the data is reached i.e. J ˆ 0; which is the original data. Step 8. The reconstructed image is used as starting solution for the original data. EM is performed till the transition criterion is satisfied. The resulting solution is grid

A. Raheja, A.P. Dhawan / Computerized Medical Imaging and Graphics 24 (2000) 359–376

(a)

(b)

(c)

(d)

(e)

367

(f)

Fig. 9. Images reconstructed using different wavelet basis in the WMREM algorithm that uses threshold on the high-pass data spaces: (a) biorthogonal spline N~ ˆ 1; N ˆ 5; (b) Daubechies orthogonal N ˆ 3; (c) biorthogonal spline N~ ˆ 6; N ˆ 8; (d) least asymmetric N~ ˆ 4; (e) biorthogonal spline N~ ˆ 3; N ˆ 1; and (f) Daubechies orthogonal N ˆ 20:

mapped to the next finer resolution (i.e. K ˆ 3) using a wavelet spline based interpolation method. Step 9. The mapped image is used as the starting solution for the finest grid and EM is performed until the transition criterion is satisfied. The final reconstructed image is taken as the LL component of the wavelet transformed image reconstructed at level K ˆ 3 and J ˆ 0: 3.3. The transition and stopping criterion Decomposing the image created using the LL component

at level j after each iteration, the frequency content of the image reconstructed at that resolution can be monitored. The information contained in each decomposed frequency band is represented by the energy. The energy is defined as: Energy ˆ

I X

…band_pass_coefficient…i††2

…31†

iˆ0

where I is the total number of pixels in each frequency band. The energy in each frequency band is monitored with

368

A. Raheja, A.P. Dhawan / Computerized Medical Imaging and Graphics 24 (2000) 359–376

(a)

(b)

(c)

(d)

Fig. 10. Images reconstructed using different wavelet basis in the WMREM algorithm that used DC shifting on the high-pass data spaces: (a) biorthogonal spline N~ ˆ 1; N ˆ 5, (b) biorthogonal spline N~ ˆ 6; N ˆ 8; (c) least asymmetric N ˆ 4; and (d) biorthogonal spline N~ ˆ 3; N ˆ 1:

respect to the normalized RMS error defined as v u I uX u …l…i† ⫺ l^ …i††2 u u iˆ1 Error ˆ u u I X u t l…i†2

…32†

the image reconstructed at the jth iteration at a certain grid level. It is observed that the RMS error decreases and reaches a minimum and then starts increasing. The transition point corresponds to an iteration before the RMS error reaches its minimum and the change in the RMS error becomes very small.

iˆ1

where l (i) is the known emission density in box i, l^ …i† is the computed density, and I is the total number of pixels in the image. The energy in the high–high frequency band represents the high-frequency content of the image at that grid level. The grid level is switched when the energy in this band ceases to increase and starts to decrease. This is done in order to retrieve the high-frequency components of the desired grid resolution, which are the low-frequency components of a grid level finer than the desired grid resolution. Mathematically, the transition criterion [3] is written as E…u j⫹1 † ⫺ E…u j † ⬍ 0:0 E…u j † ⫺ E…u j⫺1 †

and

E…u j † ⫺ E…u j⫺1 † ⬎ 0:0

…33† where E is the energy operator defined in Eq. (31), and u j is

4. Custom design of biorthogonal wavelet filters As mentioned in Section 2, one can use a set of biortho~ that can be used to create gonal wavelet filters {h; h~ 0 ; g0 ; g} ~ g; g} ~ with only the dual wavelet and the a new set {h; h; scaling filter remaining the same. Cohen et al. [23] created a set of compactly supported biorthogonal wavelets that resemble spline hat functions. This Daubechies wavelet [26] which has N~ ˆ 2 and N ˆ 2 is a special case of the second-generation wavelet created by using the lifting coefficients a ˆ 1=4 and b ˆ 1=4 and Deslauriers and Dubuc’s [25] FIR filter coefficients for the first order …N~ ˆ 1† interpolating scaling function. In this paper, we use the second order interpolating spline function …N~ ˆ 3† as the dual scaling function. The reason for choosing this function is discussed later in Section 5. The family of biorthogonal wavelets constructed by Daubechies

A. Raheja, A.P. Dhawan / Computerized Medical Imaging and Graphics 24 (2000) 359–376

369

function and the scaling function can be written as

c~ …x† ˆ f~ …2x† ⫺ f~ …2x ⫺ 1† ⫺ af~ …x ⫹ 3† ⫺bf~ …x ⫹ 2† ⫺ cf~ …x ⫹ 1† ⫺ d f~ …x† ⫺ef~ …x ⫺ 1† ⫺ f f~ …x ⫺ 2† ⫺ gf~ …x ⫺ 3†

(36)

f…x† ˆ f…2x† ⫹ f…2x ⫺ 1† ⫹ ac…x ⫹ 3† ⫹bc…x ⫹ 2† ⫹ cc…x ⫹ 1† ⫹ d c…x† ⫹ec…x ⫺ 1† ⫹ f c…x ⫺ 2† ⫹ gc…x ⫺ 3†

(37)

Table 1 lists the lifting coefficients found using different symmetry and momentum conditions. Each set gives rise to a different set of wavelet filters. Three and five lifting parameters were also used but since the wavelet for seven lifting coefficients gave good results, only a few sets of lifting coefficients have been listed. The wavelet names as mentioned in the table above have been assigned arbitrarily. The bior3lift7.1 performed best for the WMREM algorithm. The frequency response of the low-pass and high-pass analysis filters for the new wavelet and for some other commonly used wavelets are shown in Fig. 8. 5. Results

Fig. 11. Frequency response for: (a) low-pass; and (b) high-pass decomposition filter for the bior3lift7.1 compared with other commonly used wavelets Daubechies biorthogonal bior4.6, bior5.5, bior1.3 and Daubechies orthogonal d20 and d6.

[26] are special cases of second generation wavelets that can be constructed using different orders of interpolating spline functions as dual scaling functions and lifting coefficients derived using symmetry and moment conditions. The set of FIR filter coefficients for the second order spline function to be used for the lifting method are given as p h~ ˆ {1=4; 3=4; 3=4; 1=4†= 2 and …34† p g ˆ { ⫺1=4; 3=4; ⫺3=4; 1=4}= 2 h ˆ old samples

and

g~ ˆ even samples

…35†

Using the lifting method, the equations for the dual wavelet

This new algorithm was evaluated quantitatively using the normalized RMS error and qualitatively by comparison of reconstructed images and profiles through the reconstructed and original phantom. The algorithm is contrasted with the traditional EM method and the Filtered Backprojection. Commonly known wavelets were also used in the WMREM algorithm and the results are compared with the best custom designed wavelet suitable for tomographic image reconstruction. The simulated emission phantom data was created using the method described by Shepp and Logan [27]. The results presented here were generated using the phantom of 128 × 128 pixels resolution as seen in Fig. 8. It was assumed to be 128 discrete detectors equally spaced around the circle p of radius 2 circumscribing the image boxes. The total number of emissions was set to 10 million. To demonstrate the transition criterion, the energy in the HH band of the wavelet decomposed image after each iteration is plotted for each level of reconstruction grid. Daubechies “least asymmetric” [26] compactly supported wavelet with maximum number of vanishing moments for N ˆ 4 was used for analysis of the image after each iteration to develop the transition criterion [3]. This wavelet has a spectral response that allows for better quantization of the high-frequency components of an image compared with some other wavelets investigated.

370

A. Raheja, A.P. Dhawan / Computerized Medical Imaging and Graphics 24 (2000) 359–376

Fig. 12. Decomposition: (a) scaling function; (b) wavelet function and reconstruction; (c) scaling function; and (d) wavelet function for the custom designed bior3lift7.1 wavelet.

The choice of the most suitable wavelet for use in this algorithm was not trivial. Different wavelets were used and evaluated for the analysis and synthesis in the WMREM algorithm. We initially tested the algorithm using the popular Daubechies orthogonal wavelets [26] of different orders. As seen in Fig. 9b and f, the system of orthogonal wavelets does not perform well for this algorithm. The Daubechies “least asymmetric” wavelets with different number of

vanishing moments were also evaluated in this framework and these provided better results as seen in Fig. 9d and d. The Cohen et al. [23] biorthogonal spline wavelet family was used in this framework and these provided the best results (reconstructed images) as seen in Fig. 9c and e among all the commonly known wavelets evaluated. Hence, this work also brings forth the effect of using different wavelet bases in the specific problem at hand.

Table 1 Lifting coefficients using seven lifting parameters (lifting coefficient derived using various conditions of symmetry and vanishing moments) Wavelet name

Conditions used for lifting coefficients

bior3lift7.1

sym_ag, sym_bf sym_ce, mom0, mom1, mom3, mom5 sym_ag, sym_bf sym_ce, mom0, mom1, mom3, mom5 sym_ag, sym_bf sym_ce, mom0, mom1, mom3, mom5 sym_ag, sym_bf sym_ce, mom0, mom1, mom3, mom5 sym_ag, sym_bf sym_ce, mom0, mom1, mom3, mom5

bior3lift7.2

bior3lift7.3

bior3lift7.4

bior3lift7.5

a

b

c

d

e

f

g

15 2048

⫺31 512

459 2048

0

⫺459 2048

31 512

⫺15 2048

225 31744

⫺30347 507904

3501 15872

225 253952

⫺3501 15872

29897 507904

⫺225 31744

12639 1715200

⫺104309 1715200

9639 42880

⫺153 1715200

⫺76653 343040

103391 1715200

⫺6243 857600

⫺1395 6381952

⫺1559533 102111232

⫺3989475 25527808

⫺1395 51055616

3989475 25527808

1562123 102111232

1395 6381952

⫺14103 8865248

⫺3525331 141843968

57015 316616

⫺14103 70921984

⫺57015 316616

⫺3553537 141843968

14103 8865248

A. Raheja, A.P. Dhawan / Computerized Medical Imaging and Graphics 24 (2000) 359–376

371

Fig. 13. Transition curves for different grid levels for the WMREM algorithm: (a) Grid level: 32, transition iteration: 1; (b) Grid level: 64, transition iteration: 9; (c) Grid level: 128, transition iteration: 18; and (d) Grid level: 256, transition iteration: 6.

5.1. Evaluation of various wavelet bases for use in WMREM algorithm The following two cases discuss the effect of subjecting high-pass filtered data spaces to thresholding and DC shifting to provide a positive distribution to the EM estimation process. Reconstruction results using both these methods are discussed. Case I Thresholding the coefficients of data spaces n2ⴱ(d), ⴱ n3 (d) n5ⴱ(d), and n6ⴱ(d). ᭙ n j …d† ⬍ 0:0 make nj …d† ˆ 0:0;

j ˆ 2; 3; 5 and 6

Performing this operation on the data spaces, some of the coefficients are set to zero, i.e. some information in these data spaces is lost but it is required to perform this in order to provide a distribution that can be used for the EM estimation process. Some of the reconstructed images in Fig. 9 are difficult to distinguish qualitatively, hence Table 2 provides

a comparison in terms of the normalized RMS error. As seen from the errors, the orthogonal wavelets like the Daubechies do not perform well in this scheme of tomographic reconstruction. Also, the least asymmetric wavelet for order N ˆ 4 performs well but with some image degradation along the edges. The biorthogonal spline wavelet family [23,26] for N~ ˆ 1; 2; 3; 4; 5 and 6 were also tested for use in the WMREM algorithm and the family of N~ ˆ 3 provided the Table 2 RMS errors using different wavelet bases for the algorithm. The high-pass filtered data spaces use the method of thresholding Wavelet name

RMS error

Wavelet name

RMS error

Bior N~ ˆ 1; N ˆ 5 Bior N~ ˆ 3; N ˆ 1 Bior N~ ˆ 6; N ˆ 8

0.151491 0.117206 0.118325

Daubechies N ˆ 3 Least asymmetric N ˆ 4 Bior3lift7.1 N~ ˆ 3; 7 lifting parameters

0.783103 0.116828 0.114550

372

A. Raheja, A.P. Dhawan / Computerized Medical Imaging and Graphics 24 (2000) 359–376

Fig. 14. NRMS error curves for different grid levels for the WMREM algorithm: (a) 32 × 32; (b) 64 × 64; (c) 128 × 128; and (d) 256 × 256.

best results. This initiated the use of second order interpolating spline function for the custom design using the method of lifting. Case II DC shifting the coefficients of data spaces n2ⴱ(d), ⴱ n3 (d) n5ⴱ(d), and n6ⴱ(d). nj …d† ˆ abs…min† ⫹ nj …d†;

…39†

where min ˆ minimum …nj …d††: This method finds the minimum value for each data space, and DC shifts the complete data space by adding the absolute of the minimum value. The idea behind this method is to preserve the statistics of the coefficients even though the values are being changed. The reconstruction was performed on phantom data using this method and the results are displayed in Fig. 10. The quantitative measure using the RMS errors is compared in Table 3. As seen from these results, the method of thresholding performs better than DC shifting. This is primarily due to the reason that DC shifting actually changes the coefficient

values in the data spaces, which results in distortion of the information contained in the data space. Hence, the reconstructed images are of poorer quality. In case of thresholding, less information is used from the data space as a result of equating negative coefficients to zero. The method of thresholding does not change the values of the coefficients used for the estimation process i.e. even though less information is used, but without distortion of the data. Hence, the thresholding method results in better Table 3 RMS errors using different wavelet bases for the algorithm. The high-pass filtered data spaces use the method of DC shifting Wavelet name

RMS error Wavelet name

Bior N~ ˆ 1; N ˆ 5 0.322867 Bior N~ ˆ 3; N ˆ 1 0.153096 Bior N~ ˆ 6; N ˆ 8 0.136050

RMS error

Daubechies N ˆ 3 0.897563 Least Asymmetric N ˆ 4 0.151780 Bior3lift7.1 N~ ˆ 3; 7 0.118339 lifting parameters

A. Raheja, A.P. Dhawan / Computerized Medical Imaging and Graphics 24 (2000) 359–376

(a)

(b)

(c)

(d)

373

Fig. 15. (a) Original phantom, Phantom data reconstructed using (b) WMREM (using the best custom designed wavelet bior3lift7.1), (c) ML-EM and (d) Filtered Backprojection algorithm.

quality reconstructed images as compared with DC shifting method. 5.2. Custom designed wavelet The importance of the custom design of the wavelet is better understood by evaluating the reconstructed images in Figs. 9 and 10. As seen from these figures, these images differ because the localization of spatial frequency in the various data spaces plays an important role in the MR concept. Hence, it becomes important to custom design wavelets for data of a specific nature. The wavelets are being used to localize a certain range of spatial frequency in the LL, LH and HL bands that are used for reconstruction. The reconstruction can be improved if appropriate frequencies can be localized in these data spaces, i.e. include the desired high frequencies excluding the noise. Many custom wavelets were designed using the lifting method as seen in Table 1. The custom design of wavelets is steered by the

frequency response of the wavelet filters. The frequency response of the decomposition low-pass and high-pass wavelet filters of a few commonly used wavelets and the best custom designed wavelet are compared in Fig. 11. The most suitable custom designed wavelet used for analysis and synthesis in the WMREM algorithm is shown in Fig. 12. This is the bior3lift7.1 wavelet of Table 1. Fig. 11 demonstrates that the cut-off frequency of the low-pass decomposition filter plays a crucial role in choosing the wavelet to generate MR data spaces. Comparing the frequency response of the low-pass decomposition wavelet filters in this figure, it is seen that the cut-off frequency for bior3lift7.1 is shifted to the right, hence including higher frequencies as compared to other wavelet filters. Also, the magnitude response for bior3lift7.1 indicates a lower magnification of the low frequencies, which is a highly desirable property for these data spaces because the low frequencies are already provided by the LL data space. The best performing WMREM algorithm uses the custom

374

A. Raheja, A.P. Dhawan / Computerized Medical Imaging and Graphics 24 (2000) 359–376

Fig. 16. Intensity profiles for a vertical line 75 pixels from the left side through the phantom fort the WMREM and the single grid EM (ML-EM).

designed bior3lift7.1 wavelet filters. The reconstructed image using WMREM algorithm that used custom designed filter is seen in Fig. 14 and is discussed in the next two subsections.

5.3. Transition criterion The coarsest resolution grid level is 32 × 32 which uses data spaces n1ⴱ(d), n2ⴱ(d) and n3ⴱ(d), also at the same resolution. Assigning every pixel, a constant value initializes the coarsest resolution grid. The grid level 64 × 64 used by data space n5ⴱ(d) and n6ⴱ(d) is also initialized in a similar fashion. The variation of energy in the HH band and RMS error with iterations is plotted for each grid level in Figs. 13 and 14, respectively. As seen from these plots for the energy and errors, the transition point seen in the energy curves lies just before the RMS error reaches a minimum and starts to increase. The region before the errors reach a minimum is ideal for the transition to the next finer resolution since the region beyond this iteration starts contributing to the noise to the reconstruction. Fig. 13 shows that the number of iterations used to reconstruct the phantom was 1 at grid level 32 × 32 using data spaces n1ⴱ(d), n2ⴱ(d) and n3ⴱ(d), 9 at grid level 64 × 64 using data spaces n4ⴱ(d), n5ⴱ(d) and n6ⴱ(d), 18 at grid level 128 × 128 using n0ⴱ(d) i.e. original data and 6 at grid level 256 × 256 using original data. This transition criterion is more stringent and robust as compared to the maximum likelihood.

5.4. Qualitative and quantitative comparison The reconstructed phantom image is compared with the Maximum Likelihood SGEM (ML-EM) algorithm and the Filtered Backprojection algorithm. The images in Fig. 15 show that the WMREM using custom designed wavelet filters provides a better reconstruction. The ML-EM algorithm was run for the same amount of CPU time as the WMREM for a fair comparison. The RMS errors for the various algorithms are tabulated in Table 4. The lowest RMS error for the WMREM is a quantitative measure that corroborates the performance of the algorithm. A profile through the reconstructed phantom and the original phantom for qualitative comparisons is seen in Fig. 16. 5.5. Real PET data reconstruction The WMREM reconstruction algorithm was tested on actual PET camera data. The data was is a brain slice collected from an ECAT scanner. The data was also reconstructed using the ML-EM algorithm and the Filtered Backprojection algorithms. Fig. 17 demonstrates that the WMREM performs well and gives better quality reconstructed images. Table 4 RMS errors for various reconstruction methods Phantom

WMREM

Shepp

0.114

SGEM (ML-EM) 0.143

Filtered backprojection 0.583

MGEM log-likelihood 0.165

A. Raheja, A.P. Dhawan / Computerized Medical Imaging and Graphics 24 (2000) 359–376

(a)

(b)

(c) Fig. 17. Brain slice data reconstructed using the (a) WMREM (using custom designed wavelet), (b) SGEM and (c) Filtered Backprojection.

6. Conclusions We have presented a new reconstruction algorithm that utilizes the wavelet representation and spatio-frequency localization in EM based tomographic image reconstruction. This algorithm can be extended to solving inverse problems that use the EM estimation method. Wavelets have been used primarily to form data sub-spaces by performing two-channel orthonormal filter analysis on the original tube data. The HH data sub-spaces are not used since they contain most of the possible data collection noise. Wavelet synthesis with the remaining data sub-spaces is used to

375

reconstruct the image with EM estimation performed using the data sub-spaces and grid levels at the same resolution. The use of wavelet MR approach in this manner eliminates the interpolation step that was used earlier in the MREM algorithm [3], and provides better reconstruction. Wavelet-decomposition of estimated image at each iteration is also used to develop a transition and stopping criterion. This method uses Daubechies least asymmetric wavelet …N ˆ 4†; since its frequency response provides a good localization of the high-frequency contents of an image. The data spaces should have appropriate frequency components of the data for better reconstruction and these data spaces should be more de-correlated than the data itself. This implies that the wavelet filters of the selected wavelet basis should have a frequency response ideal to localize the data into MR data spaces used for reconstruction. Hence, biorthogonal wavelet filters have been custom designed for use in the wavelet-analysis and synthesis in this algorithm. To have an idea regarding the design of wavelet filters better suited for the nature of tomographic data, commonly known wavelets like the Daubechies orthogonal, biorthogonal spline, least asymmetric etc. were used in the framework of this algorithm. The orthogonal wavelets did not perform well but the Daubechies biorthogonal spline family [26] performed quite well. This led to the use of second order interpolating spline function as the dual scaling function in the lifting method to create new biorthogonal wavelet filters. The new custom designed wavelet filters used for the WMREM algorithm have a frequency response ideal for tomographic data. These wavelet filters include the appropriate higher frequencies in the data spaces required for improved image reconstruction. Also, the filters have a lower gain for the low frequencies, which is helpful in localizing more high-frequency data in the high-pass filtered data space. The new algorithm presented in this paper provides better quality images for PET image reconstruction compared to earlier MG and single grid EM methods. The algorithm can be extended to incorporate other derivatives of EM method for PET image reconstruction such as the Penalized EM [4]. The algorithm, in general, can also be used for other applications involving solutions for inverse problems that use EM as the estimation method. The algorithm can be parallelized for more efficient implementation. The estimation process using the three data spaces at each grid level can be run in parallel to save time and make the algorithm faster. The MR part of the algorithm before the original data is utilized, helps provide a good initial estimate to the remaining MG part of algorithm.

Acknowledgements The authors wish to thank Dr Kevin Wheeler of the NASA Ames Research Labs for discussion regarding custom wavelet design. We also thank Sloan Kettering Medical Center, Dayton for providing PET data.

376

A. Raheja, A.P. Dhawan / Computerized Medical Imaging and Graphics 24 (2000) 359–376

References [1] Shepp LA, Vardi Y. Maximum Likelihood reconstruction for emission tomography. IEEE Trans Med Imag 1982;MI-1(2):113–22. [2] Ranganath MV, Dhawan AP, Mullani N. A multigrid expectation maximization reconstruction algorithm for positron emission tomography. IEEE Trans Med Imag 1988;7(4):273–8. [3] Raheja A, Dhawan AP. Multiresolution Expectation Maximization reconstruction algorithm for PET using wavelet processing. Proc IEEE EMBS Conf 1998;20(2/6):759–62. [4] Alvaro R De Pierro. A modified EM algorithm for penalized likelihood estimation in emission tomography. IEEE Trans Med Imag 1995;14(1):273–8. [5] Fessler JA, Hero AO. Penalized maximum likelihood image reconstruction using space-alternating generalized EM algorithms. IEEE Trans Imag Proces 1995;4:1417–29. [6] Kay Jim. The EM algorithm in medical imaging. Stat Meth Med Res 1997;6:55–75. [7] Wang G, Zhang J, Pan G. Solution of inverse problems in image processing by wavelet expansion. IEEE Trans Imag Proces 1995;4(5): 579–93. [8] Zhu W, Wang Y, Deng Y, Yao Y, Barbour RL. A wavelet-based multiresolution regularized least squares reconstruction approach for optical tomography. IEEE Trans Med Imag 1997;16(2):210–7. [9] Kolaczyk ED. A wavelet shrinkage approach to tomographic image reconstruction. J Am Stat Assoc 1996;91(435):1079–90. [10] Peyrin F, Zaim M, Goutte R. Construction of wavelet decompositions for tomographic images. J Math Imag Vision 1993;3:105–21. [11] Delaney AH, Bresler Y. Multiresolution tomographic reconstruction using wavelets. IEEE Trans Imag Proces 1995;4(6):799–813. [12] Bhatia M, Karl WC, Willsky AS. A wavelet-based method for multiscale tomographic reconstruction. IEEE Trans Med Imag 1996;15(1):92–101. [13] Hackbusch W. Multi-grid methods and applications. New York: Springer, 1985. [14] Sweldens W. The lifting scheme: a construction of second generation wavelets. Technical Report 1995:6, Department of Mathematics, University of South Carolina, 1995 (ftp://ftp.math.sc.edu/pub/ imi_95/imi95_6.ps). [15] Sweldens W. The lifting scheme: a custom-design construction biorthogonal wavelets. Appl Comput Harmon Anal 1996;3(2):186–200. [16] Sweldens W, Schroder P. Building your own wavelets at home. Technical Report 1995:5, Department of Mathematics, University of South Carolina, 1995 (ftp://ftp.math.sc.edu/pub/imi_95/imi95_6.ps). [17] Burrus CS, Gopinath RA, Guo Haitao. Introduction to wavelets and wavelet transforms: a primer. Englewood Cliffs, NJ: Prentice-Hall, 1998. [18] Vetterli M, Kovacevic J. Wavelets and subband coding. Englewood Cliffs, NJ: Prentice-Hall, 1995. [19] Chui CK. An introduction to wavelets. New York: Academic Press, 1992. [20] Mallat S. A theory for multiresolution signal decomposition: the wavelet representation. IEEE Trans Pattern Anal Machine Intell 1989;11:674–93. [21] Unser M, Aldroubi A, Eden M. A family of polyspline wavelet transforms. Signal Proces 1993;30(2):141–62. [22] Wheeler K. Smoothing non-uniform data samples using wavelets. PhD thesis, Department of Electrical Engineering, University of Cincinnati, 1996. [23] Cohen A, Daubechies I, Feauveau JC. Biorthogonal bases of compactly supported wavelets. Commun Pure Appl Math 1992;45(5):485–560. [24] Chui CK. An introduction to wavelets, wavelet analysis and its applications, vol. 1. New York: Academic Press, 1992. [25] Deslauriers G, Dubuc S. Symmetric iterative interpolation processes. Constructive Approx 1989;5(1):49–68.

[26] Daubechies I. Ten lectures on wavelets. Philadelphia, PA: SIAM, 1992. [27] Shepp LA, Logan BF. The Fourier reconstruction of a head section. Technical Report 92, Department of Computer Science, SUNY, Buffalo, NY, 1975.

Amar Raheja obtained BS in Physics in 1992 and MS in Physics in 1994, both from Indian Institute of Technology, Kharagpur, India. He was with the Department of Radiology, University of Texas Southwestern Medical Center, Dallas from 1994–1997 as a PhD candidate in Radiological Sciences. He obtained his PhD in Bioengineering at University of Toledo with a specialization in application of wavelets in medical imaging in 1999. He is currently an Assistant Professor of Computer Science at Philadelphia University. His current research interests are image and signal processing using wavelets, medical imaging and soft computing (neural networks and genetic algorithms). He is an associate member of IEEE, the IEEE engineering in medicine and biology society and the IEEE computer society.

Atam P. Dhawan obtained his BEng and MEng degrees in Electrical Engineering from the University of Roorkee, Roorkee, India. He was a Canadian Commonwealth Fellow at the University of Manitoba where he completed his PhD in Electrical Engineering with specialization in medical imaging and image analysis in 1985. In 1984, he won the first prize and the Martin Epstein Award in the Symposium of Computer Application in Medical Care Paper Competition at the Eighth SCAMC Annual Congress in Washington, DC, for his work on developing a 3D imaging technique to detect early skin-cancer called melanoma. From 1985 to 1988, he was an Assistant Professor in the Department of Electrical Engineering at the University of Houston. Later, in 1988, he joined the University of Cincinnati as an Assistant Professor where he became Professor of Electrical and Computer Engineering and Computer Science, and Radiology (joint appointment). From 1990–1996, he was the Director of Center for Intelligent Vision and Information System. From 1996–1998, he was Professor of Electrical Engineering at the University of Texas at Arlington, and Adjunct Professor of Radiology at the University of Texas Southwestern Medical Center at Dallas. He is currently Professor in the Department of Electrical & Computer Engineering at the New Jersey Institute of Technology and Director of Medical Imaging and Informatics Laboratory. Dr Dhawan has published more than 50 research articles in refereed journals, and edited books, and 85 research papers in refereed conference proceedings. Dr Dhawan is a recipient of Martin Epstein Award (1984), National Institutes of Health FIRST Award (1988), Sigma-Xi Young Investigator Award (1992), University of Cincinnati Faculty Achievement Award (1994) and the prestigious IEEE Engineering in Medicine and Biology Early Career Achievement Award (1995). He is an Associate Editor of IEEE Transactions on Biomedical Engineering, Associate Editor of IEEE Transactions on Rehabilitation Engineering, and Editor of International Journal of Computing Information and Technology. He has served on many IEEE EMBS professional committees and has delivered Workshops on Intelligent Biomedical Image Analysis in IEEE EMBS International Conferences (1996 and 1997). He is the Chair of the “New Frontiers in Biomedical Engineering” Symposium at the World Congress 2000 on Medical Physics and Biomedical Engineering. His current research interests are medical imaging, multimodality brain mapping, intelligent image analysis, multigrid image reconstruction, wavelets, genetic algorithms, neural networks, adaptive learning and pattern recognition.