A Prospective Study on Algorithms Adapted to the Spatial Frequency in Tomography

The use of iterative algorithms in tomographic reconstruction always leads to a frequency adapted rate of convergence in that low frequencies are accurately reconstructed after a few iterations, while high frequencies sometimes require many more computations. In this paper, we propose to build frequency adapted (FA) algorithms based on a condition of incomplete backprojection and propose an FA simultaneous algebraic reconstruction technique (FA-SART) algorithm as an example. The results obtained with the FA-SART algorithm demonstrate a very fast convergence on a highly detailed phantom when compared to the original SART algorithm. Though the use of such an FA algorithm may seem difficult, we specify in which case it is relevant and propose several ways to improve the reconstruction process with FA algorithms.


INTRODUCTION
Iterative methods are increasingly used in computed tomography (image reconstruction from projections) where in most cases they advantageously replace analytical methods. The latter methods, whose success is based upon a great robustness and speed, unfortunately cannot take into account a certain number of physical and geometrical parameters which should appear in the point spread function (PSF). These parameters (intrinsic resolution of the detector, spatial variation of sensitivity, and others, related to specific applications) can be easily incorporated in the projection operator that models the process of imaging in iterative methods. This however does not change the nature of the problem to be solved, formulated as a large system of linear equations: where A is the projection operator, x is typically a set of (unknown) values taken by the voxels in the 3D space (attenuation coefficients, activity . . . ), and b is the set of measurements, often organized as a set of projection images. In addition, besides the exponential growth in computer processing power that obviously plays an important role in the feasibility of large computations, solving methods for (1) profit from a broad effort realized in general inverse prob-lems methods and from a continuous exchange between theory and applications.
One important ongoing study on the speed of convergence is contingent on the consideration of different algorithms (algebraic reconstruction technique (ART) [1], simultaneous algebraic reconstruction technique (SART) [2], expectation maximization (EM) [3], or conjugated gradient for the most famous), different types of implementation (sequential (Seq), simultaneous (Sim), or block-iterative (BI)), regularizations, or relaxation parameters (see [4] for details). A solid mathematical framework was formed, which in most cases tendered proof of convergence, including at times for the inconsistent case-generally that of computed tomography.
The idea of the work presented here has risen from the fact that the high frequencies (the details) are always reconstructed after the low frequencies in iterative methods for computed tomography. In the case of highly detailed objects (in a sense that we will later specify), an important number of iterations might be required to obtain the precision expected on the reconstructed image. This very general fact has been experienced by all who worked with iterative algorithms, and keys for the understanding of this phenomena can be found in [5], but in a framework of integral geometry resulting in a Dirac projection PSF and thus quite different from our context.

International Journal of Biomedical Imaging
In this article, we propose a way to accelerate the convergence of iterative methods in the specific case of highly detailed objects, based on the use of an incomplete backprojection operator. This can be roughly formulated as a weighting of each correction d k i, j of voxel j by equation i during the kth iteration with weights w k i, j such as i w k i, j = 1 for all j and k and the incomplete backprojection condition (IBC): To our knowledge, this condition has not been thus far reported in the use of reconstruction algorithms for computed tomography. Besides the formulation of a frequency adapted (FA) algorithm that realizes the reconstruction with IBC starting from the well-known SART algorithm, a major concern of this paper is to explain how the IBC permits to accelerate the convergence by updating the value of the voxels from the most relevant measurements and to give experimental results to support both the convergence of the algorithm and the gain obtained when compared to SART.
In Section 2, we describe the type of problems for which we believe the adaptation involved by the IBC can be useful and formulate the IBC in terms of projections onto convex sets. We explain the IBC for it to realize the expected acceleration and give an adaptation of the simultaneous algebraic reconstruction technique (SART) algorithm of Andersen and Kak [2] modified by the IBC. In Section 3, we show results on simulations, phantom, and small animal with a comparison between SART and FA-SART. Section 4 is a discussion on the proposed method which also suggests a number of different use for FA algorithms and Section 5 concludes the paper.

BUILDING AN ALGORITHM ADAPTED TO THE SPATIAL FREQUENCY
In computed tomography like in other applications leading to inverse problems, a set of discrete measurements is physically obtained, presumably by the application of some linear operator A to a studied object, discretized into a set of unknowns. For each unknown x j (a voxel in tomography), the point spread function (PSF) (we should actually talk about the voxel spread function) appears in the column A •, j of the matrix A whose coefficients are assumed to be nonnegative here, to simplify the notations. For some imaging techniques, the paradigmatic case of which is single photon emission computed tomography (SPECT), this PSF is always different from a Dirac function since it notably includes the detector components (collimator, scintillator, photo-multiplicators, electronic boards, etc.). Indeed, the photons emitted from a voxel j reach far more than one single pixel on the detector and one has to take this into account in the matrix A when modeling this problem. PSF in SPECT is generally modelled by discretized Gaussian functions but some authors proposed more elaborate models, especially for pinhole SPECT (see Metzler et al. [6]).
But whatever the chosen model is, the actual situation is that two close voxels have many measurements in com- Object space Figure 1: Overlapping between two close voxels' PSF. The darkest zone will be used to reconstruct voxel x j as well as x j+4 (and all those in between), thus leading to a smoothed reconstruction after only a few iterations. mon in their PSF (see Figure 1). This may result in a slow rate of convergence: the backprojection of these shared measurements leads to a smoothing for the voxels' values, that only tends to vanish as the number of iterations increases. This is not necessarily inadequate since when a large and homogeneous structure is to be imaged, this smoothing will mainly reduce the noise if one stops the iterations before the noise is reconstructed (Figure 2(a)). But when considering an image of projection with high frequencies, using these shared measurements in the backprojection, or at least the furthest away from the center of the PSF, will hold the voxels' values distant from the solution x * for a certain number of iterations ( Figure 2(b)). Throughout this paper, we will characterize higly detailed objects in an informal manner by reference to this figure, that is, when significant details of the image are of the same size or smaller than the PSF's width.
Our idea is to only use the most central measurements of the PSF to accelerate the convergence in the case of imaging with high spatial frequencies. To do so, we considered the SART algorithm whose iteration is where j = 1, . . . , N indexes the voxels and i = 1, . . . , M indexes the pixels so that x k ∈ R N for all k, b ∈ R M , A is an MxN matrix modelling the problem and λ is the relaxation parameter. This can also be written in a matrix form discussed by Jiang and Wang [7] (T denotes the transpose of a matrix or a vector): with  Figure 2: Comparison of the two main cases, when a profile corresponding to a large and homogeneous structure is measured with noise (a), and when the measured profile corresponds to a higly detailed object (b). On (a), using all the discrete measurements of voxel's PSF will smooth the reconstruction which will reduce the noise. On (b), "important" high frequencies are smoothed and more iterations will be needed to achieve an accurate reconstruction.

an NxN diagonal matrix and
Although this might sound antinomic with the simultaneous of SART, we can write a component Seq form of SART that can be derived from the analysis presented in [8] by Censor and Elfving. Then, applying this algorithm in the matrix form of (4), it is well known that the new estimate x k+1 is the projection of x k onto the hyperplane This is our starting point for a frequency adapted algorithm based on SART. For the purpose of the incomplete backprojection condition, given 0 ≤ ρ ≤ 1, let A ρ denote the matrix so that A ρ only keeps in each column the coefficients of A higher than ρ times the top of the PSF. We also define the ma- . Then, following (4) we define the next iterate corresponding to the treatement of a single projection i as It can been showen in the exact same manner as for (4) (in the sequential case) that P ρ i (x) projects x onto the hyperplane H i . But (7) is also designed to have P ρ i (x) in the affine subspace of R N : thus projecting the current estimate onto C k i := H i ∩ Ω k i , which is a closed convex set of R N . Now, we want to consider the scheme (7) in its natural case, that is, simultaneous. Thus we change the iteration (3) for the following. x k+1 The BI version requires some sophistication in the notations and for this purpose, we adopt pretty much the same ones as Censor and Elfving in [8]. Let then T be the number of blocks and, for t = 1, 2, . . . , T, let the blocks of indices B t ⊆ {1, 2, . . . , M} be an ordered subset of the form is the number of elements in B t . In our case, we define T = M and for each t, B t = {t}. For t = 1, 2, . . . , T, let A t denote the matrix formed by taking all the rows {a i } of A whose indices belong to the block of indices B t , that is, A t := a t . The same applies to A ρ with A ρ t the matrix formed by taking all the rows { a i } of A ρ whose indices belong to the block of indices B t and a The vector b is partitioned similarly with b t denoting the elements of b whose indices belong to the block of indices B t , that is, b t := b t . Let us also denote a j,t(k) c the jth column of A t , then the BI FA-SART writes as follows.

Algorithm 2.2 (BI FA-SART).
x k+1 These algorithms realize a backprojection characterized by a contracted cone compared to the original SART (see Figure 3). These results show improvements with FA-SART when the PSF is large (2.5 mm and 1.5 mm) and better images with SART when the data are noisy (Gaussian noise, 20% of the maximum of the data), especially for big structures. No significant difference is visible when the PSF is small, because the stage of deconvolution is of less importance during the reconstruction process.

Material and method
The applications shown here are simulations and pinhole SPECT acquisitions made on a small animal single head dedicated SPECT gamma camera (Gaede Medizinsysteme GMBH, Freiburg, Germany) with a 6.5 mm NaI(Tl) crystal 25 photomultipliers and a small field of view of 17 cm × 17 cm. The camera was equipped with a tungsten pinhole collimator of 120 mm in focal length and 1.5 mm in diameter for the phantom, 2 mm in diameter for the mouse study. A more complete description of the device, as well as results on small animal imaging can be found in [9]. Prior to reconstruction, we used a correction of the center of rotation applied on the projections after measurement of the defect on a single-line-source phantom as described in [10]. Both the FA-SART and SART algorithms were block iterative with each block corresponding to a different image of projection.

Simulation
The simulation performed here on a computed phantom illustrates the typical behavior of the FA-SART algorithm, which accelerates the convergence when the PSF is large (see Figures 4(a), 4(b) and 4(e), 4(f)). For a thin PSF, no visible difference appears since FA-SART acts on the deconvolution. Bigger structures should not be reconstructed with ρ close to 1 from noisy data to avoid having the high frequencies of the noise emphasized. Three iterations were computed for all reconstructions.

Phantom
A phantom, seen in Figure 5, with three types of cylindrical cavities (diameters 1 mm, 1.5 mm and 2 mm separated with the same distance), was filled with 99m Tc, 20 MBq activity.
The parameters for the acquisition were 30 mm radius of rotation, 60 projections of 64 × 64 pixels on 180 • and 1 minute per projection. For the comparison in Figures 5(b)-5(d), we used the FA-SART algorithm with ρ = 1 and the SART algorithm (which is FA-SART with ρ = 0). For this highly detailed object, the FA-SART reconstruction (5(b)) does better than 3 iterations of SART (5(c)) or even 10 iterations (5(d)), due to an important smoothing effect during the first iterations. Figure 5(e) provides a promising way to use FA-SART in a preconditioning of system (1). It shows a reconstruction obtained after 3 iterations of FA-SART (ρ = 1) followed by 2 iterations of SART.

Cardiac mouse imaging
For mouse heart perfusion, a normal adult female CD1 mouse (Mouse Clinical Institute, Ilkirch, France) weighting 30 g was injected with 400 MBq of 99mTc-Tetrofosmin (Amersham, General Electric Healthcare, USA). The radius of rotation of the camera was 2.5 cm corresponding to a zoom factor of 5 and 48 projections of 64 × 64 pixels were acquired over 180 • (see [9] for details). Figure 6 demonstrates a satisfying result with FA-SART and ρ = 0.8 that makes both the left and right ventricles visible.

DISCUSSION
From the obtained results, it seems that the FA-SART algorithm is able to reconstruct high details in less iterations than SART. While the reconstructions performed from physical structures (phantom and animal) show a faster convergence to what seems to be the true distribution, those performed from simulated data demonstrate that this acceleration increase with the PSF's width of the system, that is, when the intrinsic resolution decreases. These simulations also prove that the improvement obtained with the FA-SART algorithm Vincent Israel-Jost et al. is not due to an inadequacy of the matrix A that models the system. Thus, the method only applies to a class of problems with large PSF, that is, involving an important step of deconvolution like SPECT or optical tomography. It would make no sense, for those applications where the PSF is close to a Dirac function like in computed tomography scanning, to contract the backprojection cone since this cone is already almost a line. It should also be recalled that the gain obtained on highly detailed structures is balanced by a loss in large and homogeneous structures, for which the smoothing properties of the SART algorithm reduce the noise. On the contrary, using the FA-SART algorithm with a high threshold ρ would increase the noise by emphasizing the high frequencies in the reconstruction. This leads to the major difficulty of such an algorithm: how it should be used. We have shown results with a same value for ρ for all the voxels and all the iterations. This is a first way to accomodate to the FA-SART: the context of SPECT is that of functional imaging, in which in general only one structure is considered interesting. Thus, even if several structures appear in the same image, the parameter ρ can be set to a value corresponding to a given organ or a type of examination (myocardial perfusion produces a sharper image of the heart than blood-pool imaging so these two explorations of the same organ would require different values of ρ). This has been for the most part our way to proceed thus far, although other uses can be considered.
A use in preconditioning of the system (1) with two or three iterations of FA-SART followed by a solving with the SART algorithm seems to combine the best of these two algorithms by providing both a fast reconstruction of the high frequencies and an adequate rendering of the bigger structures. Also a value of ρ adapted to the local frequencies of the projection images might successfully change the local properties of the reconstruction, with a high value of ρ only when high-frequencies are detected. We started to investigate this way by creating local frequencies maps of the projection images composed of wavelets coefficients and this will constitute the base of our further works.
A noticeable positive point of the FA-SART algorithm is that it requires very few adaptations when SART has already been implemented. On a programming point of view, certain loops are just partially executed, which also means that one iteration of FA-SART is slightly faster than one of SART.
Since the backprojection operation significantly differs from other studied algorithms, it is by no means obvious to deduce that the FA-SART algorithm converges from an analysis of SART or of the general Landweber scheme. Results concerning the convergence have nevertheless been greatly improved by the suggestion of an anonymous referee and a paper dedicated to a mathematical analysis will be presented elsewhere for publication. This theoretical work will support the experimental evidences shown here, with a study of the convergence of the algorithm but also of its behavior in function of the main frequency.

CONCLUSION
The FA-SART algorithm has been designed to be able to change the usual rate of convergence of iterative algorithms. While low frequencies are generally the first reconstructed, we permitted to invert this phenomenon in order to have the details appear in the image after very few iterations. We evidenced this behavior of the FA-SART algorithm by showing applications and proposed several possible ways to use it. Theoretical questions will be discussed in a forthcoming paper. His work concerns the production or improvement of iterative algorithms in cone-beam SPECT. His research interests include iterative algorithms and image processing, as well as philosophy of science and image depiction. He is now pursuing a second Master's degree in the philosophy of science from the University of Paris I -La Sorbonne.
Philippe Choquet, veterinary since 1988, got a predoctoral degree in medical imaging and signal processing in 1991, and a Ph.D. degree dissertation in 1996. Since 2001, he has been an Associate Professor in the Biophysics and Nuclear Medicine Department, CHU Hautepierre, Strasbourg. He is the author and coauthor of 29 papers in peerreviewed journals as well as 76 oral and poster communications in national and international meetings in the field of low field MRI, NMR, and MRI of hyperpolarized gases and small animal scintigraphy.