- Email: [email protected]
Computationally efficient filtered-backprojection algorithm for tomographic image reconstruction using Walsh transform
J. Vis. Commun. Image R. 17 (2006) 581–588 www.elsevier.com/locate/jvci
Computationally efficient filtered-backprojection algorithm for tomographic image reconstruction using Walsh transform Gylson Thomas *, V.K. Govindan Department of Computer Engineering, National Institute of Technology, Calicut 673 601, Kerala, India Received 31 March 2005; accepted 20 February 2006 Available online 11 April 2006
Abstract In this paper, we discuss the implementation of the filtered-backprojection (FBP) algorithm for tomographic image reconstruction using Walsh transform to exploit its fast computational ability. Walsh transform is the fastest unitary transform known so far. The major advantage of Walsh transform is that it involves only real additions and subtractions whereas Fourier transform involves complex multiplications and additions. Implementation of the proposed algorithm necessitates the design of an appropriate filter in Walsh domain. In this research, the known Fourier filter coefficients have been transformed into Walsh domain, thereby the 1 · N Fourier filter coefficients were converted into an N · N sparse matrix with nonzero elements in a special pattern. The proposed algorithm has been implemented by taken into account of the special nature of the Walsh domain filter coefficients and tested for its performance using the well-known ‘SheppLogan head phantom’ test image. The results demonstrate that the reconstruction strategy has comparable performance with a significant reduction of computing time. For example, with a 128 · 128-pixel image and 180 views, the speedup achieved is fourfold, with reconstructions qualitatively and visually the same as that of FBP algorithm in the Fourier domain. 2006 Elsevier Inc. All rights reserved. Keywords: Tomography; Walsh transform; Fast algorithm; Filtered-backprojection algorithm
1. Introduction The problem of reconstructing a two/three-dimensional function from measurements of line integrals has applications in such diverse areas as computerized tomography (CT), electron microscopy, radioastronomy, optical interferometry, earth resources exploration, nondestructive testing, radar, sonar, and ultrasonic imaging. The inverse Radon transform is an exact solution to this problem given by Radon in 1917 [1], for the case of continuous exact parallel line integrals taken over a continuous rotation about an origin. In CT-scanning the number of X-ray photons transmitted through an object of interest is measured by an array of detectors. This array is rotated by a fixed angular increment between each set of measurements. From these *
Corresponding author. Fax: +91 0495 2287250. E-mail addresses: [email protected], [email protected] (G. Thomas).
1047-3203/$ - see front matter 2006 Elsevier Inc. All rights reserved. doi:10.1016/j.jvcir.2006.02.003
582
G. Thomas, V.K. Govindan / J. Vis. Commun. Image R. 17 (2006) 581–588
measurements the total attenuation along each measurement paths is derived and an image of the linear attenuation coefficients may be reconstructed. A detailed description of this application and several approach to solutions are discussed in [2,3]. In many of the potential applications of tomographic medical imaging, it is of the interest to achieve ‘realtime’ imaging. As always, the meaning of ‘real-time’ is application dependent. In medical image reconstruction, where the data collection devices tend to be very expensive, it is important that the utilization of the devices is not impaired by the time required to compute the images. It is also often desirable to have these images available while the patient is still in the device, to avoid having to recall the patient due to the initial imaging procedure not meeting the diagnostic need. Even in applications where real-time imaging is not of the essence, fast computation of images is desired so that the imaging procedure can be evaluated and optimized without outrageous demands on computational resources [4]. The filtered-backprojection (FBP) algorithm [5] is the most common method to reconstruct a two-dimensional image from its parallel line projections as the detector rotates 180 around the object. In this algorithm, the projection data are first convolved with a ramp-filter kernel at each projection angle, then backprojected into the image space. This algorithm is efficient and provides a stable image with stationary noise propagation. In well-controlled environments, it is generally believed that the FBP algorithm provides a reconstructed image of high quality with normally available computer capacity and computation times [6,7]. Traditionally, for medical applications, backprojection has been accomplished by special hardware that exploits parallelism of the backprojection process to achieve near real-time reconstruction of a single slice [8]. However, there still exists a lag in reconstruction time that is becoming increasingly important with the emergence of new technologies able to acquire data at ever-faster rates such as in fifth generation CT-scanners. With the increasing data rates, the FBP algorithm becomes the bottleneck [9] in the reconstruction and display process, and presents a barrier to real-time imaging. Hence computational speed of the reconstruction algorithms is of great importance in the tomographic imaging field [10]. In this aspect, we have attempted to implement the FBP algorithm using Walsh transform in place of the widely accepted Fourier transform to exploit its fast computational ability. Walsh transform is the fastest unitary transform known so far. The transform involves only real additions and subtractions and can be computed using fast algorithm of complexity O (N log2N) [11,12]. A filter has been generated in the Walsh domain using the known Fourier filter coefficients [13]. Computer simulation of the FBP algorithm has been carried out both in the Fourier and Walsh domains and the results are compared. We limit ourselves to the case with parallel beams. Modern machines use the so called fan-beam projections [14], where many nonparallel beams emitted from one point in the form of a fan are considered. The fan-beam geometry has the advantage over parallel-beam geometry that it can collect multiple angles simultaneously in a short time interval. The fan-beam algorithm in the Walsh domain is left as a future work. The organization of the paper is as follows: In Section 2, we discuss the theoretical background of the tomographic image reconstruction problem along with the widely accepted FBP algorithm. Section 3 deals with the implementation of FBP algorithm in the Walsh domain. The comparative study of the algorithms both in the Fourier and Walsh domain is presented in Section 4 followed by a conclusion in Section 5. 2. Filtered-backprojection algorithm for parallel projection data Let a two-dimensional function f (x, y) represents a cross-section of the human body (Fig. 1). Any line running through the cross-section is called a ray and the integral of f (x, y) along a ray is called a ray integral and a set of ray integrals forms a projection [2]. The equation of line AB in Fig. 1 is given by x cos h þ y sin h ¼ t1 ;
ð1Þ
where t1 is the perpendicular distance of the line from the origin. The integral of the function f (x, y) along this line may be expressed as Z Z 1 Z 1 P h ðt1 Þ ¼ f ðx; yÞ ds ¼ f ðx; yÞdðx cos h þ y sin h t1 Þ dx dy; ð2Þ ray AB
1
1
G. Thomas, V.K. Govindan / J. Vis. Commun. Image R. 17 (2006) 581–588
583
A
Fig. 1. This figure illustrates the variables used in the discussion on the filtered-backprojection algorithm for parallel projection data. The function Ph(t) is a parallel projection of f (x, y).
where d(Æ) is the Dirac delta function. The function Ph(t) as a function of t (for a given value of h) defines the parallel projection of f (x, y) for angle h. The two-dimensional function Ph(t) is also called the Radon transform of f (x, y). Here, the problem is to recover the function f (x, y) from the known function Ph(t). Among the number of reconstruction algorithms existing to solve this problem, the FBP algorithm is the widely accepted one because of its high accuracy and speed. Fundamental to a number of reconstruction algorithms is the Fourier slice theorem which states that the one-dimensional Fourier transform of a projection is equal to the center cross-section of the two-dimensional Fourier transform of the image. Derived from this theorem [2] the FBP algorithm is expressed as Z p f ðx; yÞ ¼ Qh ðx cos h þ y sin hÞ dh; ð3Þ 0
where Qh(t) is called the filtered projection and are related to the projections Ph(t), by Z 1 Qh ðtÞ ¼ P h ðaÞhðt aÞ da;
ð4Þ
1
where the filter impulse response h(t) is the inverse Fourier transform of H(x) = |x|, as shown in Fig. 2, which is a ramp function in the frequency domain over the bandwidth (W) of the system given by H( ) 1/2
1/2
1/2
Frequency ( )
Fig. 2. The ideal filter response for the FBP algorithm, bandlimited to W = (1/2s).
584
G. Thomas, V.K. Govindan / J. Vis. Commun. Image R. 17 (2006) 581–588
hðtÞ ¼
Z
w
jxjej2pxt dx.
ð5Þ
w
These equations suggest the following steps for a digital implementation of the algorithm. 2.1. Step 1 (filtering) It is assumed that each projection is sampled with a sampling interval of ‘s’ cm. To avoid aliasing problems due to sampling, the value of W is taken as 1/2s. Substituting this value of W in Eq. (5), we get for the impulse response, 2 1 sinð2pt=2sÞ 1 sinðpt=2sÞ hðtÞ ¼ 2 2 . ð6Þ 2s 2pt=2s 4s pt=2s Now since the data are measured with a sampling interval of s, and is correspondingly band limited, for digital processing the impulse response need only be known with the same bandwidth and hence the same sampling interval. So from Eq. (6), we have 8 1 n ¼ 0; > < 4s2 ; n is even, ð7Þ hðnsÞ ¼ 0; > : 1 n2 p2 s2 ; n is odd, where n takes both negative and positive integer values. The frequency domain implementation may be expressed as Qh ðnsÞ ¼ s IFFTfFFTðP h ðnsÞ with ZPÞ FFTðhðnsÞ with ZPÞg;
ð8Þ
where FFT denotes fast Fourier transform, IFFT the inverse of FFT and ZP stands for zero padding. 2.2. Step 2 (backprojection) The second step deals with reconstruction of the image from the filtered projections using a digital approximation to the integral in Eq. (3). When the number of projections Mproj is large enough and uniformly distributed over 180, Eq. (3) may be approximated by f^ ðx; yÞ ¼
M proj p X Q ðx cos hi þ y sin hi Þ. M proj i¼1 hi
ð9Þ
The function f^ ðx; yÞ is the reconstructed approximation to the original function f (x, y). 3. Filtered-backprojection algorithm in the Walsh domain Each discrete transform has its own advantages and disadvantages over other transforms. Fourier domain filters are easier to design than Walsh domain filters. But Walsh transform can be computed much faster than Fourier transforms [13]. By making use of the superior computational properties of each transforms, image reconstruction algorithms with better performance can be obtained. Walsh functions are orthogonal functions defined over the interval (0, 1). These are rectangular pulses having only two discrete values +1 and 1. The Walsh transform is significantly more efficient to compute than the Fourier transform because each operation involves only real addition or subtraction, compared to the complex multiplication necessary for the Fourier transform. The discrete Walsh transform (DWT) [11] of the data sequence, X(n) = {X(0), X(1), . . . , X(N 1)} is defined as W x ðnÞ ¼
1 H w ðnÞX ðnÞ; N
ð10Þ
G. Thomas, V.K. Govindan / J. Vis. Commun. Image R. 17 (2006) 581–588
585
where Wx(k) is the kth DWT coefficient and Hw(n) is the N · N Walsh-ordered Hadamard matrix. Since Hw(n) is orthogonal and symmetric, the inverse transform IDWT is given by, X(n) = Hw(n)Wx(n). As per Eq. (8), the filtering in Fourier domain can be extended to the Walsh domain as Qh ðnsÞ ¼ s IFWTfFWTðP h ðnsÞ with ZPÞ Gw g;
ð11Þ
where FWT denotes fast Walsh transform, IFWT the inverse of FWT, Gw the Walsh domain filter coefficients and ZP stands for zero padding. The major concern here is the design of Walsh domain filter coefficients. We have adopted here a domain transformation method for generating Walsh filter coefficients. 3.1. Transformation of Fourier filter coefficients into Walsh domain In Figs. 3A and B, x and x 0 represent the input and output discrete signals, respectively, where x; x0 2 RN and N = 2q (or q = log2N), q is a positive integer. In Fig. 3A, F is the N · N discrete Fourier transform matrix and F1 is its inverse [15]. The Fourier gain matrix Gf is a complex diagonal N · N matrix of the form Gf ¼ diagfg0 ; g1 ; . . . ; gN 1 g;
ð12Þ
where gi ¼ gN i ; i = 1, 2, . . . , (N/2) 1 holds. The asterisk means complex conjugate. Also, g0 and gN/2 are real. In Fig. 3B W is the N · N discrete Walsh transform (DWT) matrix and W1 is its inverse. Gw, the Walsh gain matrix, is analogous to Gf. The problem here is to find Gw in terms of Gf. It is obvious that the two systems of Figs. 3A and B produce the same x 0 for the same x. This allows us to write F 1 Gf F ¼ W 1 Gw W
ð13Þ
and, upon solving for Gw, Gw ¼
1 WF 1 Gf FW N
ð14Þ
since W 1 ¼ N1 W : In the case of FBP algorithm, since Gf is a ramp function, Gw can be easily computed using Eq. (13). During the transformation process, the one-dimensional Fourier filter coefficients were replaced by a twodimensional filter coefficient matrix in Walsh domain which is the major hurdle in our algorithm (Tables 1a and 1b). It is observed that most of the off-diagonal elements are zeros and the nonzero elements are in a special pattern as shown in Table 1b. This factor is also considered during the development of the algorithm to minimize computations. The Walsh filter coefficients can be computed off-line and stored once and for all. After filtering, the backprojection operation can be performed in the same manner as in Fourier case.
A
x
F
Gf
F-1
FFT
FOURIER GAIN MATRIX
IFFT
W
Gw
W-1
FWT
WALSH GAIN MATRIX
IFWT
x′
B
x
Fig. 3. (A) Spectral filtering via FFT. (B) Spectral filtering via FWT.
x′
586
G. Thomas, V.K. Govindan / J. Vis. Commun. Image R. 17 (2006) 581–588
Table 1a Fourier domain ramp-filter coefficients (Gf) for N = 16 0
0.125
0.25
0.375
0.5
0.625
0.75
0.875
1
0.875
0.75
0.625
0 0 0.09 0 0 0 0.08 0 0 0 0.59 0 0 0 0.11 0
0 0 0 0.18 0 0 0 0 0 0 0 0.68 0 0 0 0
0.5
0.375
0.25
0.125
Table 1b Equivalent Walsh domain filter coefficients (Gw) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0.2 0 0 0 0.11 0 0 0 0.09 0 0 0 0.11 0 0
0 0 0.2 0 0 0 0.11 0 0 0 0.09 0 0 0 0.11 0
0 0 0 0.32 0 0 0 0 0 0 0 0.18 0 0 0 0
0 0 0 0 0.32 0 0 0 0 0 0 0 0.18 0 0 0
0 0.11 0 0 0 0.41 0 0 0 0.08 0 0 0 0.09 0 0
0 0 0.11 0 0 0 0.41 0 0 0 0.08 0 0 0 0.09 0
0 0 0 0 0 0 0 0.5 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0.5 0 0 0 0 0 0 0
0 0.09 0 0 0 0.08 0 0 0 0.59 0 0 0 0.11 0 0
0 0 0 0 0.18 0 0 0 0 0 0 0 0.68 0 0 0
0 0.11 0 0 0 0.09 0 0 0 0.11 0 0 0 0.8 0 0
0 0 0.11 0 0 0 0.09 0 0 0 0.11 0 0 0 0.8 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1
4. Simulations and performance evaluation The proposed algorithm is implemented in MATLAB7.0TM environment to test the performance of it using the ‘Shepp-Logan head phantom’ as the test image [5]. The test image resembles the human head cross-section. It is actually a composition of 10 ellipses of different parameters. The advantage of using such an image is that the projection data can be easily generated using analytical expressions. The performance of the proposed algorithm is compared with the widely accepted FBP algorithm in Fourier domain for various sizes of the head phantom image with 180 projections. Table 2 shows the computational time taken by both the methods for various sizes of the test image. From Table 2 it is obvious that the Walsh domain algorithm requires approximately one-fourth of computing time compared to the Fourier domain algorithm. Since the Walsh domain filter coefficients are two-dimensional in nature, the performance becomes poor as the image size increases. Table 3 demonstrates that the root mean square error (RMSE) and peak signal to noise ratio (PSNR) remain the same in both domains which ensure the quality of the image. Fig. 4 shows the original head phantom image and the reconstructed images using the FBP algorithm in both the Fourier and Walsh domains. It shows that the visual perception is the same in both algorithms. Hardware implementation may further improve the computational speed because Walsh transform involves only real additions and subtractions whereas Fourier transform involves the complex additions and multiplications.
Table 2 Comparison of computational time Image size
(A) FBPa algorithm Fourier domain (time, s)
(B) FBPa algorithm Walsh domain (time, s)
(A)/(B)
64 · 64 128 · 128 256 · 256 512 · 512
2.188 4.875 11.829 30.109
0.453 1.062 3.203 10.984
4.83 4.59 3.69 2.74
a
Filtered-backprojection.
G. Thomas, V.K. Govindan / J. Vis. Commun. Image R. 17 (2006) 581–588
587
Table 3 Comparison of reconstruction error Image size
128 · 128
a b
No. of projections (127 rays per projection)
180 90 45 18
Fourier domain
Walsh domain
RMSEa
PSNRb (dB)
RMSEa
PSNRb (dB)
0.1222 0.1329 0.1469 0.1659
+18.26 +17.21 +16.93 +15.61
0.1222 0.1329 0.1469 0.1659
+18.26 +17.21 +16.93 +15.61
Root mean square error. Peak signal to noise ratio.
Fig. 4. (A) The 128 · 128 pixel original head phantom image. Reconstructed image from 180 projections with 127 rays per projection using FBP algorithm (B) in Fourier domain. (C) in Walsh domain.
5. Conclusion In this paper, we have addressed the implementation of filtered-backprojection algorithm using Walsh transform in place of the widely accepted Fourier transform. Numerical results illustrate the superiority of the cost of computation in Walsh domain compared to the Fourier domain, with an approximate speedup of fourfold. It is demonstrated that the root mean square error and peak signal to noise ratio remain the same in both Walsh and Fourier domains, which ensure the quality of the reconstructed image. Also demonstrated that the visual perception of the reconstructed images in both domains the same. The major concern of this paper is the design of an appropriate filter in the Walsh domain. In this research, a domain transformation method has been adopted which results in a two-dimensional filter coefficients in the Walsh domain. This two-dimensional nature of the Walsh domain filter coefficients adversely affect the speed of the image
588
G. Thomas, V.K. Govindan / J. Vis. Commun. Image R. 17 (2006) 581–588
reconstruction process as the image size increases. A parallel computing procedure or a one-dimensional filter in Walsh domain could solve this problem. The above work can be carried as a future work as an extension of the present research work. References [1] J. Radon, On the determination of functions from their integral values along certain manifolds, IEEE Transactions on Medical Imaging 5 (1986) 170–176 (Translated by P.C. Parks from the original German text). [2] A.C. Kak, M. Slaney, Principles of Computerized Tomographic Imaging, IEEE Press, New York, 1988. [3] G.T. Herman, Image Reconstruction From Projections—The Fundamentals of Computerized Tomography, Academic, New York, 1980. [4] G.T. Herman, Image reconstruction from projections, Real-Time Imaging 1 (1995) 3–18. [5] L.A. Shepp, B.F. Logan, The Fourier reconstruction of a head section, IEEE Transactions on Nuclear Science 21 (1974) 21–42. [6] H.J. Scudder, Introduction to computer aided tomography, Proceedings of the IEEE 66 (6) (1978) 628–637. [7] M. Robert, R.M. Lewitt, Reconstruction algorithms: transform methods, Proceedings of IEEE 71 (1983) 390–408. [8] B.K. Gilbert, S.K. Kenue, R.A. Robb, A. Chu, A.H. Lent, E.S. Earl, Rapid execution of fan beam image reconstruction algorithms using efficient computational techniques and special-purpose processors, IEEE Transactions on Biomedical Engineering 28 (2) (1981) 98–116. [9] S. Matej, J.A. Fessler, G.K. Ivan, Iterative tomographic image reconstruction using Fourier-based forward and back-projectors, IEEE Transactions on Medical Imaging 23 (4) (2004) 401–412. [10] B. Samit, B. Yorams, O (N2 log2N) filtered-backprojection reconstruction algorithm for tomography, IEEE Transactions on Image Processing 9 (10) (2000) 1760–1773. [11] N. Ahmed, K.R. Rao, Orthogonal Transformations for Digital Signal Processing, Spring Verlag, New York, 1975. [12] W.K. Pratt, J. Kane, H.C. Andrews, Hadamard transform image coding, Proceedings of the IEEE 57 (1) (1969) 58–68. [13] E.L. Hall, Sequency domain design of frequency filters, IEEE Transactions on Computers 23 (5) (1974) 976–981. [14] B.K.P. Horn, Fan-beam reconstruction methods, Proceedings of the IEEE 67 (1979) 1616–1623. [15] J.Z. Christopher, Y. Maurice, Spectral filtering using the fast Walsh transform, IEEE Transactions on Acoustic Speech and Signal Processing 33 (4) (1985) 1246–1252.
Computationally efficient filtered-backprojection algorithm for tomographic image reconstruction using Walsh transform Gylson Thomas *, V.K. Govindan Department of Computer Engineering, National Institute of Technology, Calicut 673 601, Kerala, India Received 31 March 2005; accepted 20 February 2006 Available online 11 April 2006
Abstract In this paper, we discuss the implementation of the filtered-backprojection (FBP) algorithm for tomographic image reconstruction using Walsh transform to exploit its fast computational ability. Walsh transform is the fastest unitary transform known so far. The major advantage of Walsh transform is that it involves only real additions and subtractions whereas Fourier transform involves complex multiplications and additions. Implementation of the proposed algorithm necessitates the design of an appropriate filter in Walsh domain. In this research, the known Fourier filter coefficients have been transformed into Walsh domain, thereby the 1 · N Fourier filter coefficients were converted into an N · N sparse matrix with nonzero elements in a special pattern. The proposed algorithm has been implemented by taken into account of the special nature of the Walsh domain filter coefficients and tested for its performance using the well-known ‘SheppLogan head phantom’ test image. The results demonstrate that the reconstruction strategy has comparable performance with a significant reduction of computing time. For example, with a 128 · 128-pixel image and 180 views, the speedup achieved is fourfold, with reconstructions qualitatively and visually the same as that of FBP algorithm in the Fourier domain. 2006 Elsevier Inc. All rights reserved. Keywords: Tomography; Walsh transform; Fast algorithm; Filtered-backprojection algorithm
1. Introduction The problem of reconstructing a two/three-dimensional function from measurements of line integrals has applications in such diverse areas as computerized tomography (CT), electron microscopy, radioastronomy, optical interferometry, earth resources exploration, nondestructive testing, radar, sonar, and ultrasonic imaging. The inverse Radon transform is an exact solution to this problem given by Radon in 1917 [1], for the case of continuous exact parallel line integrals taken over a continuous rotation about an origin. In CT-scanning the number of X-ray photons transmitted through an object of interest is measured by an array of detectors. This array is rotated by a fixed angular increment between each set of measurements. From these *
Corresponding author. Fax: +91 0495 2287250. E-mail addresses: [email protected], [email protected] (G. Thomas).
1047-3203/$ - see front matter 2006 Elsevier Inc. All rights reserved. doi:10.1016/j.jvcir.2006.02.003
582
G. Thomas, V.K. Govindan / J. Vis. Commun. Image R. 17 (2006) 581–588
measurements the total attenuation along each measurement paths is derived and an image of the linear attenuation coefficients may be reconstructed. A detailed description of this application and several approach to solutions are discussed in [2,3]. In many of the potential applications of tomographic medical imaging, it is of the interest to achieve ‘realtime’ imaging. As always, the meaning of ‘real-time’ is application dependent. In medical image reconstruction, where the data collection devices tend to be very expensive, it is important that the utilization of the devices is not impaired by the time required to compute the images. It is also often desirable to have these images available while the patient is still in the device, to avoid having to recall the patient due to the initial imaging procedure not meeting the diagnostic need. Even in applications where real-time imaging is not of the essence, fast computation of images is desired so that the imaging procedure can be evaluated and optimized without outrageous demands on computational resources [4]. The filtered-backprojection (FBP) algorithm [5] is the most common method to reconstruct a two-dimensional image from its parallel line projections as the detector rotates 180 around the object. In this algorithm, the projection data are first convolved with a ramp-filter kernel at each projection angle, then backprojected into the image space. This algorithm is efficient and provides a stable image with stationary noise propagation. In well-controlled environments, it is generally believed that the FBP algorithm provides a reconstructed image of high quality with normally available computer capacity and computation times [6,7]. Traditionally, for medical applications, backprojection has been accomplished by special hardware that exploits parallelism of the backprojection process to achieve near real-time reconstruction of a single slice [8]. However, there still exists a lag in reconstruction time that is becoming increasingly important with the emergence of new technologies able to acquire data at ever-faster rates such as in fifth generation CT-scanners. With the increasing data rates, the FBP algorithm becomes the bottleneck [9] in the reconstruction and display process, and presents a barrier to real-time imaging. Hence computational speed of the reconstruction algorithms is of great importance in the tomographic imaging field [10]. In this aspect, we have attempted to implement the FBP algorithm using Walsh transform in place of the widely accepted Fourier transform to exploit its fast computational ability. Walsh transform is the fastest unitary transform known so far. The transform involves only real additions and subtractions and can be computed using fast algorithm of complexity O (N log2N) [11,12]. A filter has been generated in the Walsh domain using the known Fourier filter coefficients [13]. Computer simulation of the FBP algorithm has been carried out both in the Fourier and Walsh domains and the results are compared. We limit ourselves to the case with parallel beams. Modern machines use the so called fan-beam projections [14], where many nonparallel beams emitted from one point in the form of a fan are considered. The fan-beam geometry has the advantage over parallel-beam geometry that it can collect multiple angles simultaneously in a short time interval. The fan-beam algorithm in the Walsh domain is left as a future work. The organization of the paper is as follows: In Section 2, we discuss the theoretical background of the tomographic image reconstruction problem along with the widely accepted FBP algorithm. Section 3 deals with the implementation of FBP algorithm in the Walsh domain. The comparative study of the algorithms both in the Fourier and Walsh domain is presented in Section 4 followed by a conclusion in Section 5. 2. Filtered-backprojection algorithm for parallel projection data Let a two-dimensional function f (x, y) represents a cross-section of the human body (Fig. 1). Any line running through the cross-section is called a ray and the integral of f (x, y) along a ray is called a ray integral and a set of ray integrals forms a projection [2]. The equation of line AB in Fig. 1 is given by x cos h þ y sin h ¼ t1 ;
ð1Þ
where t1 is the perpendicular distance of the line from the origin. The integral of the function f (x, y) along this line may be expressed as Z Z 1 Z 1 P h ðt1 Þ ¼ f ðx; yÞ ds ¼ f ðx; yÞdðx cos h þ y sin h t1 Þ dx dy; ð2Þ ray AB
1
1
G. Thomas, V.K. Govindan / J. Vis. Commun. Image R. 17 (2006) 581–588
583
A
Fig. 1. This figure illustrates the variables used in the discussion on the filtered-backprojection algorithm for parallel projection data. The function Ph(t) is a parallel projection of f (x, y).
where d(Æ) is the Dirac delta function. The function Ph(t) as a function of t (for a given value of h) defines the parallel projection of f (x, y) for angle h. The two-dimensional function Ph(t) is also called the Radon transform of f (x, y). Here, the problem is to recover the function f (x, y) from the known function Ph(t). Among the number of reconstruction algorithms existing to solve this problem, the FBP algorithm is the widely accepted one because of its high accuracy and speed. Fundamental to a number of reconstruction algorithms is the Fourier slice theorem which states that the one-dimensional Fourier transform of a projection is equal to the center cross-section of the two-dimensional Fourier transform of the image. Derived from this theorem [2] the FBP algorithm is expressed as Z p f ðx; yÞ ¼ Qh ðx cos h þ y sin hÞ dh; ð3Þ 0
where Qh(t) is called the filtered projection and are related to the projections Ph(t), by Z 1 Qh ðtÞ ¼ P h ðaÞhðt aÞ da;
ð4Þ
1
where the filter impulse response h(t) is the inverse Fourier transform of H(x) = |x|, as shown in Fig. 2, which is a ramp function in the frequency domain over the bandwidth (W) of the system given by H( ) 1/2
1/2
1/2
Frequency ( )
Fig. 2. The ideal filter response for the FBP algorithm, bandlimited to W = (1/2s).
584
G. Thomas, V.K. Govindan / J. Vis. Commun. Image R. 17 (2006) 581–588
hðtÞ ¼
Z
w
jxjej2pxt dx.
ð5Þ
w
These equations suggest the following steps for a digital implementation of the algorithm. 2.1. Step 1 (filtering) It is assumed that each projection is sampled with a sampling interval of ‘s’ cm. To avoid aliasing problems due to sampling, the value of W is taken as 1/2s. Substituting this value of W in Eq. (5), we get for the impulse response, 2 1 sinð2pt=2sÞ 1 sinðpt=2sÞ hðtÞ ¼ 2 2 . ð6Þ 2s 2pt=2s 4s pt=2s Now since the data are measured with a sampling interval of s, and is correspondingly band limited, for digital processing the impulse response need only be known with the same bandwidth and hence the same sampling interval. So from Eq. (6), we have 8 1 n ¼ 0; > < 4s2 ; n is even, ð7Þ hðnsÞ ¼ 0; > : 1 n2 p2 s2 ; n is odd, where n takes both negative and positive integer values. The frequency domain implementation may be expressed as Qh ðnsÞ ¼ s IFFTfFFTðP h ðnsÞ with ZPÞ FFTðhðnsÞ with ZPÞg;
ð8Þ
where FFT denotes fast Fourier transform, IFFT the inverse of FFT and ZP stands for zero padding. 2.2. Step 2 (backprojection) The second step deals with reconstruction of the image from the filtered projections using a digital approximation to the integral in Eq. (3). When the number of projections Mproj is large enough and uniformly distributed over 180, Eq. (3) may be approximated by f^ ðx; yÞ ¼
M proj p X Q ðx cos hi þ y sin hi Þ. M proj i¼1 hi
ð9Þ
The function f^ ðx; yÞ is the reconstructed approximation to the original function f (x, y). 3. Filtered-backprojection algorithm in the Walsh domain Each discrete transform has its own advantages and disadvantages over other transforms. Fourier domain filters are easier to design than Walsh domain filters. But Walsh transform can be computed much faster than Fourier transforms [13]. By making use of the superior computational properties of each transforms, image reconstruction algorithms with better performance can be obtained. Walsh functions are orthogonal functions defined over the interval (0, 1). These are rectangular pulses having only two discrete values +1 and 1. The Walsh transform is significantly more efficient to compute than the Fourier transform because each operation involves only real addition or subtraction, compared to the complex multiplication necessary for the Fourier transform. The discrete Walsh transform (DWT) [11] of the data sequence, X(n) = {X(0), X(1), . . . , X(N 1)} is defined as W x ðnÞ ¼
1 H w ðnÞX ðnÞ; N
ð10Þ
G. Thomas, V.K. Govindan / J. Vis. Commun. Image R. 17 (2006) 581–588
585
where Wx(k) is the kth DWT coefficient and Hw(n) is the N · N Walsh-ordered Hadamard matrix. Since Hw(n) is orthogonal and symmetric, the inverse transform IDWT is given by, X(n) = Hw(n)Wx(n). As per Eq. (8), the filtering in Fourier domain can be extended to the Walsh domain as Qh ðnsÞ ¼ s IFWTfFWTðP h ðnsÞ with ZPÞ Gw g;
ð11Þ
where FWT denotes fast Walsh transform, IFWT the inverse of FWT, Gw the Walsh domain filter coefficients and ZP stands for zero padding. The major concern here is the design of Walsh domain filter coefficients. We have adopted here a domain transformation method for generating Walsh filter coefficients. 3.1. Transformation of Fourier filter coefficients into Walsh domain In Figs. 3A and B, x and x 0 represent the input and output discrete signals, respectively, where x; x0 2 RN and N = 2q (or q = log2N), q is a positive integer. In Fig. 3A, F is the N · N discrete Fourier transform matrix and F1 is its inverse [15]. The Fourier gain matrix Gf is a complex diagonal N · N matrix of the form Gf ¼ diagfg0 ; g1 ; . . . ; gN 1 g;
ð12Þ
where gi ¼ gN i ; i = 1, 2, . . . , (N/2) 1 holds. The asterisk means complex conjugate. Also, g0 and gN/2 are real. In Fig. 3B W is the N · N discrete Walsh transform (DWT) matrix and W1 is its inverse. Gw, the Walsh gain matrix, is analogous to Gf. The problem here is to find Gw in terms of Gf. It is obvious that the two systems of Figs. 3A and B produce the same x 0 for the same x. This allows us to write F 1 Gf F ¼ W 1 Gw W
ð13Þ
and, upon solving for Gw, Gw ¼
1 WF 1 Gf FW N
ð14Þ
since W 1 ¼ N1 W : In the case of FBP algorithm, since Gf is a ramp function, Gw can be easily computed using Eq. (13). During the transformation process, the one-dimensional Fourier filter coefficients were replaced by a twodimensional filter coefficient matrix in Walsh domain which is the major hurdle in our algorithm (Tables 1a and 1b). It is observed that most of the off-diagonal elements are zeros and the nonzero elements are in a special pattern as shown in Table 1b. This factor is also considered during the development of the algorithm to minimize computations. The Walsh filter coefficients can be computed off-line and stored once and for all. After filtering, the backprojection operation can be performed in the same manner as in Fourier case.
A
x
F
Gf
F-1
FFT
FOURIER GAIN MATRIX
IFFT
W
Gw
W-1
FWT
WALSH GAIN MATRIX
IFWT
x′
B
x
Fig. 3. (A) Spectral filtering via FFT. (B) Spectral filtering via FWT.
x′
586
G. Thomas, V.K. Govindan / J. Vis. Commun. Image R. 17 (2006) 581–588
Table 1a Fourier domain ramp-filter coefficients (Gf) for N = 16 0
0.125
0.25
0.375
0.5
0.625
0.75
0.875
1
0.875
0.75
0.625
0 0 0.09 0 0 0 0.08 0 0 0 0.59 0 0 0 0.11 0
0 0 0 0.18 0 0 0 0 0 0 0 0.68 0 0 0 0
0.5
0.375
0.25
0.125
Table 1b Equivalent Walsh domain filter coefficients (Gw) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0.2 0 0 0 0.11 0 0 0 0.09 0 0 0 0.11 0 0
0 0 0.2 0 0 0 0.11 0 0 0 0.09 0 0 0 0.11 0
0 0 0 0.32 0 0 0 0 0 0 0 0.18 0 0 0 0
0 0 0 0 0.32 0 0 0 0 0 0 0 0.18 0 0 0
0 0.11 0 0 0 0.41 0 0 0 0.08 0 0 0 0.09 0 0
0 0 0.11 0 0 0 0.41 0 0 0 0.08 0 0 0 0.09 0
0 0 0 0 0 0 0 0.5 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0.5 0 0 0 0 0 0 0
0 0.09 0 0 0 0.08 0 0 0 0.59 0 0 0 0.11 0 0
0 0 0 0 0.18 0 0 0 0 0 0 0 0.68 0 0 0
0 0.11 0 0 0 0.09 0 0 0 0.11 0 0 0 0.8 0 0
0 0 0.11 0 0 0 0.09 0 0 0 0.11 0 0 0 0.8 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1
4. Simulations and performance evaluation The proposed algorithm is implemented in MATLAB7.0TM environment to test the performance of it using the ‘Shepp-Logan head phantom’ as the test image [5]. The test image resembles the human head cross-section. It is actually a composition of 10 ellipses of different parameters. The advantage of using such an image is that the projection data can be easily generated using analytical expressions. The performance of the proposed algorithm is compared with the widely accepted FBP algorithm in Fourier domain for various sizes of the head phantom image with 180 projections. Table 2 shows the computational time taken by both the methods for various sizes of the test image. From Table 2 it is obvious that the Walsh domain algorithm requires approximately one-fourth of computing time compared to the Fourier domain algorithm. Since the Walsh domain filter coefficients are two-dimensional in nature, the performance becomes poor as the image size increases. Table 3 demonstrates that the root mean square error (RMSE) and peak signal to noise ratio (PSNR) remain the same in both domains which ensure the quality of the image. Fig. 4 shows the original head phantom image and the reconstructed images using the FBP algorithm in both the Fourier and Walsh domains. It shows that the visual perception is the same in both algorithms. Hardware implementation may further improve the computational speed because Walsh transform involves only real additions and subtractions whereas Fourier transform involves the complex additions and multiplications.
Table 2 Comparison of computational time Image size
(A) FBPa algorithm Fourier domain (time, s)
(B) FBPa algorithm Walsh domain (time, s)
(A)/(B)
64 · 64 128 · 128 256 · 256 512 · 512
2.188 4.875 11.829 30.109
0.453 1.062 3.203 10.984
4.83 4.59 3.69 2.74
a
Filtered-backprojection.
G. Thomas, V.K. Govindan / J. Vis. Commun. Image R. 17 (2006) 581–588
587
Table 3 Comparison of reconstruction error Image size
128 · 128
a b
No. of projections (127 rays per projection)
180 90 45 18
Fourier domain
Walsh domain
RMSEa
PSNRb (dB)
RMSEa
PSNRb (dB)
0.1222 0.1329 0.1469 0.1659
+18.26 +17.21 +16.93 +15.61
0.1222 0.1329 0.1469 0.1659
+18.26 +17.21 +16.93 +15.61
Root mean square error. Peak signal to noise ratio.
Fig. 4. (A) The 128 · 128 pixel original head phantom image. Reconstructed image from 180 projections with 127 rays per projection using FBP algorithm (B) in Fourier domain. (C) in Walsh domain.
5. Conclusion In this paper, we have addressed the implementation of filtered-backprojection algorithm using Walsh transform in place of the widely accepted Fourier transform. Numerical results illustrate the superiority of the cost of computation in Walsh domain compared to the Fourier domain, with an approximate speedup of fourfold. It is demonstrated that the root mean square error and peak signal to noise ratio remain the same in both Walsh and Fourier domains, which ensure the quality of the reconstructed image. Also demonstrated that the visual perception of the reconstructed images in both domains the same. The major concern of this paper is the design of an appropriate filter in the Walsh domain. In this research, a domain transformation method has been adopted which results in a two-dimensional filter coefficients in the Walsh domain. This two-dimensional nature of the Walsh domain filter coefficients adversely affect the speed of the image
588
G. Thomas, V.K. Govindan / J. Vis. Commun. Image R. 17 (2006) 581–588
reconstruction process as the image size increases. A parallel computing procedure or a one-dimensional filter in Walsh domain could solve this problem. The above work can be carried as a future work as an extension of the present research work. References [1] J. Radon, On the determination of functions from their integral values along certain manifolds, IEEE Transactions on Medical Imaging 5 (1986) 170–176 (Translated by P.C. Parks from the original German text). [2] A.C. Kak, M. Slaney, Principles of Computerized Tomographic Imaging, IEEE Press, New York, 1988. [3] G.T. Herman, Image Reconstruction From Projections—The Fundamentals of Computerized Tomography, Academic, New York, 1980. [4] G.T. Herman, Image reconstruction from projections, Real-Time Imaging 1 (1995) 3–18. [5] L.A. Shepp, B.F. Logan, The Fourier reconstruction of a head section, IEEE Transactions on Nuclear Science 21 (1974) 21–42. [6] H.J. Scudder, Introduction to computer aided tomography, Proceedings of the IEEE 66 (6) (1978) 628–637. [7] M. Robert, R.M. Lewitt, Reconstruction algorithms: transform methods, Proceedings of IEEE 71 (1983) 390–408. [8] B.K. Gilbert, S.K. Kenue, R.A. Robb, A. Chu, A.H. Lent, E.S. Earl, Rapid execution of fan beam image reconstruction algorithms using efficient computational techniques and special-purpose processors, IEEE Transactions on Biomedical Engineering 28 (2) (1981) 98–116. [9] S. Matej, J.A. Fessler, G.K. Ivan, Iterative tomographic image reconstruction using Fourier-based forward and back-projectors, IEEE Transactions on Medical Imaging 23 (4) (2004) 401–412. [10] B. Samit, B. Yorams, O (N2 log2N) filtered-backprojection reconstruction algorithm for tomography, IEEE Transactions on Image Processing 9 (10) (2000) 1760–1773. [11] N. Ahmed, K.R. Rao, Orthogonal Transformations for Digital Signal Processing, Spring Verlag, New York, 1975. [12] W.K. Pratt, J. Kane, H.C. Andrews, Hadamard transform image coding, Proceedings of the IEEE 57 (1) (1969) 58–68. [13] E.L. Hall, Sequency domain design of frequency filters, IEEE Transactions on Computers 23 (5) (1974) 976–981. [14] B.K.P. Horn, Fan-beam reconstruction methods, Proceedings of the IEEE 67 (1979) 1616–1623. [15] J.Z. Christopher, Y. Maurice, Spectral filtering using the fast Walsh transform, IEEE Transactions on Acoustic Speech and Signal Processing 33 (4) (1985) 1246–1252.