


Iterative Methods in Image Reconstructionby Paola Favati, Grazia Lotti and Ornella Menchi A research group of scientists from the Departments of Informatics and Mathematics of Pisa University and from the Institute IITCNR, is working in image reconstruction problems in the areas of medicine and astronomy. It was found that practical iterative solutions are possible by using suitable preconditioners which increase the rate of convergence without mixing the signal (ie the original image) and the noise (ie small errors on the data). Image reconstruction is an inverse problem consisting of the identification of the input of a given instrument from the knowledge of its output. Inverse problems are generally illposed: they might not have a solution in the strict sense, solutions might not be unique and/or might not depend continuously on the data. The discrete version of the image reconstruction problem is a linear algebraic system, that is wellposed in the sense that it has a unique solution. However, this solution is completely corrupted by the noise, and is physically unacceptable because it generally does not reproduce the data within the experimental errors. This difficulty can be overcome by ‘regularization’, ie, by incorporating some expected physical properties of the object in terms, for example, of constraints on the size or on the smoothness of the solution. Starting from this idea, Tikhonov regularization methods solve leastsquared problems depending on a regularizing parameter. As an alternative some iterative methods, which enjoy an interesting regularization property known as semiconvergence, can be used. A semiconvergent method starts reconstructing the lowfrequency components of the solution; then, as the iteration progresses, the highfrequency components are reconstructed together with the noise components. Hence the method must be stopped before it starts to reconstruct the noise. The classical conjugate gradient (CG) method, which applies to symmetric positive definite systems, has the regularizing property. Due to illconditioning of the problem, the number of iterations required by CG to obtain a satisfactory result can be large and a preconditioning technique is required to increase the rate of convergence. This rate depends on the number of distinct eigenvalues of the matrix associated with the problem, and general purpose preconditioners are designed to cluster all eigenvalues around 1. In the present context, this type of preconditioner would be harmful, mixing the signal and the noise, whereas only the eigenvalues corresponding to the lowfrequency components of the signal should be clustered. In many applications, the function which describes how the imaging system affects the points of the original image (the socalled point spread function) is space invariant with respect to translations (in this case it is determined by the image of a single point source) and is bandlimited (ie it has a local action, a single source of light is blurred into a small spot, wherever the source is located).
These properties result in a strong structure of the matrix of the linear system, which turns out to be a twolevel band Toeplitz matrix. A band matrix has non zero elements only on a few diagonals around the principal one, and a Toeplitz matrix has equal elements on each diagonal. A twolevel band Toeplitz matrix is a block matrix which presents these two structures (band and Toeplitz form) both at the block level and inside each block. Circulant preconditioners, frequently applied to Toeplitz matrices, can be easily modified in order to cope with the noise. However, in the case of band Toeplitz matrices, a band preconditioner would be preferable, since it could be inverted with the same cost of a CG iteration. We have recently proposed a twolevel band preconditioner, which is effective for image reconstruction problems with the above properties and has a computational cost per iteration linear with respect to the number of pixels of the image. Figure 2 shows an example of a synthetic medical image (the 2D Hoffman phantom) blurred by the instrument used for acquisition and corrupted by noise. Figure 3 shows the image reconstructed by applying a few iterations of a preconditioned CG method. Please contact: 