IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 9, NO. 9, SEPTEMBER 2000 1573 A Multilevel Domain Decomposition Algorithm for Fast Reprojection of Tomographic Images Amir Boag, Senior Member, IEEE, Yoram Bresler, Fellow, IEEE, and Eric Michielssen, Senior Member, IEEE Abstract—A novel algorithm for fast computation of tomo- graphic image projections is presented. The method comprises a decomposition of an image into subimages followed by an aggre- gation of projections computed for the subimages. The multilevel domain decomposition algorithm is formulated as a recursive procedure. The computational cost of the proposed algorithm is comparable to that of FFT-based techniques while it appears to be more flexible than the latter. Numerical results demonstrate the effectiveness of the method. I. INTRODUCTION T HE reconstruction problem in two-dimensional (2-D) axial computed tomography (CT) is to recover an image from a set of its line-integral projections at different angles [1], [2]. The reverse operation, of computing a set of projections of a given 2-D image, is called reprojection, or computation of (angular samples of) the Radon transform. Reprojection is of interest in several applications. In X-ray CT, it is used [3] in iterative beam-hardening correction algorithms [4], in streak suppression algorithms [5], in algorithms for the removal of artifacts caused by the presence of metallic implants in the subject [6] or other high density objects [7], and in correcting for missing data [8], and partial volume effects [9]. In PET and SPECT imaging, reprojection has been used to compensate for attenuation [10]. Reprojection is also used in detection and compensation of various acquisition errors [11], [12] including those caused by patient motion [13], [14], and in accounting for Poisson noise statistics [15]. Reprojection is also a key component in iterative tomographic reconstruction algorithms, which offer various advantages over direct recon- struction methods in imaging modalities such as PET, SPECT, nondestructive testing, and more generally, in all situations of limited data (cf., [16], [17] and the references therein) and nonparallel-beam geometries [18]. Finally, reprojection is an important component in the efficient implementation of a wide Manuscript received March 24, 1999; revised April 6, 2000. This work was supported in part by grants from AFOSR via the MURI program (F49620-96-1- 0025) and by the NSF under Grant 95-02-138. The associate editor coordinating the review of this manuscript and approving it for publication was Prof. Jeffrey A. Fessler. A. Boag is with the Department of Electrical Engineering-Physical Elec- tronics, Tel-Aviv University, Ramat-Aviv 69978, Tel-Aviv, Israel (e-mail: boag@eng.tau.ac.il). Y. Bresler is with the Coordinated Science Laboratory and Department of Electrical and Computer Engineering, University of Illinois at Urbana-Cham- paign, Urbana, IL 61801 USA (e-mail: ybresler@uiuc.edu). E. Michielssen is with the Department of Electrical and Computer Engi- neering, University of Illinois at Urbana-Champaign, Urbana, IL 61801 USA (e-mail: emichiel@uiuc.edu). Publisher Item Identifier S 1057-7149(00)06977-3. variety of image processing algorithms [19], [20], including that for computing the well known Hough transform. Direct algorithms for reprojection [12], [18], [21], [22] re- main, in spite of various recent improvements, computationally very expensive, requiring operations for generating projections of an -pixel image. Likewise, the so-called “fast reprojection algorithms,” including those based on back- projection and implemented on fast hardware backprojectors [23], and those implemented on other parallel architectures [19], [20], still require operations. Their cost therefore scales as , and they also offer little advantage in implementations on general-purpose computers. These various algorithms, whether in software or hardware, are very expensive for typ- ical state-of-the-art image sizes. For example, using one of these methods to compute the reprojection of a -pixel image requires times the computation needed for a -pixel image. There are, however, two families of reprojection algorithms with cost. The first of these is the family of FFT- based algorithms, cf. [3], [24], and the references therein. These algorithms are based on the well-known Fourier Slice-Projec- tion Theorem [1], [2] (see Property 3 in the next section). The primary problem with this approach is the required step of inter- polation between the rectangular grid in Fourier space on which the transform of the image is computed using the FFT, and the polar grid on which the Fourier transform of the projections must be evaluated. The modern versions of these algorithms [24] overcome many of the associated difficulties, and offer good ac- curacy at a considerable speedup. The second family of reprojection algorithms [25]–[30] is based on a common principle. It uses a hierar- chical decomposition of the line integral in a given direction into shorter line integral in the same direction called “segments” or “links.” Adjacent segments are added up to create a seg- ment of double length. Furthermore, line integrals at adjacent directions have some segments in common, or the segments of one such line integral can be accurately interpolated from those of the other. This reuse of segments in successive dou- bling steps accounts for the computational efficiency of these algorithms. Starting with segments corresponding to line in- tegrals through two pixels, after such doubling steps, the complete line integrals are available. The earlier versions of these methods [25]–[28] used nearest-neighbor interpolation and computed the discrete Radon transform (partial sums of pixels whose center lies within a strip of predetermined width). However, versions that are more recent introduced linear [29] or 1057–7149/00$10.00 © 2000 IEEE