Iterative Methods

1. Investigation of Image Reconstruction Methods
2. The Mathematics of Image Reconstruction
3. Our research at the MIRG
4. Demonstration of the Reconstruction Process

Investigation of Image Reconstruction Methods

The problem of reconstructing medical images from measurements of the radiation around the body of a patient belong to the class of inverse problems. Inverse problems are characterized by the fact that the information of interest (e.g. the distribution of radioactivity inside the patient, for SPECT) is not directly available. The imaging device (the camera) provides measurements of a transformation of this information. In practice, these measurements are both incomplete (sampling) and inaccurate (statistical noise). In practice, this means that one must give up recovering the exact image. Indeed, aiming for full recovery of the information usually results in unstable solutions. This means that the reconstructed image is very sensitive to inevitable measurement error. Otherwise expressed, slightly different data would have produced a significantly different image. One of the main problems in reconstructing a SPECT image is that the photons hitting the camera result from scatter and background radiation, as well as from the radioactivity distribution of interest. Hopefully, the data obtained by the SPECT camera are mostly the result of emission events which carry relevant information, but one must account for the scatter and background effects as well. That means the data received by the SPECT camera are a noisy version of the ideal data. In order to cope with this difficulty (and some others), we are led to define the reconstructed image as the solution of an optimization problem. This is a way to be faithful enough to the measured data without being too sensitive to their inherent noise. If no control over the reconstruction error is exercised, reconstruction artifacts may be misinterpreted, leading to erroneous diagnosis. Thus, image reconstruction is full-time work, and the computations we do here at MIRG are cutting-edge.

The Mathematics of Image Reconstruction (for the strong-hearted)

The relationship between the unknown distribution (or object) and the physical quantity which can be measured (the projections) is referred to as the forward relationship. For several imaging techniques, such as X-ray tomography or SPECT, the simplest model of forward relationship involves the Radon transform

or some of its generalizations (attenuated Radon transform, X-ray transform, etc.). If denotes the unknown distribution and the quantity measured by the imaging device, we have

The discrete version of the latter equation can be written as

(1)

The vector is a sampling of : for SPECT, each component of represents the number of counts obtained at a given bin while the camera was position at a given angle. The vector is a finite representation of the unknown object: in the case of SPECT, each component of is the level of activity in a given pixel (2D imaging) or voxel (3D imaging). The matrix is a discretized version of the Radon transform. In general, is not invertible. This suggests that to solve (1) in the least square sense, i.e. to solve the optimization problem

This amounts to solving the normal equation :

The latter equation often has infinitely many solutions. One may then consider applying a selection rule in order to define a unique solution. For example, one may choose, among all vectors minimizing the norm of , the vector having minimum norm (which is referred to as the minimum norm least square solution) or the vector having maximum entropy. This amounts to solving the optimization problem

(2)

Here,

is the norm of

or the negative of the entropy of

or any other criteria leading to a unique solution. Unfortunately, such solutions are usually unstable. In other words, a small perturbation of the data leads to non negligible perturbation of the solution. We must therefore resort to regularization theory, in order to define a unique and stable solution. Formally, a slight modification of (2) will provide simultaneously uniqueness and stability : let us define the reconstructed object (which we shall denote

in the sequel) as the solution of the following optimization problem :

The neg-entropy plays the role of a stabilizer. It is sometimes called regularizer. If is carefully selected, the corresponding solution will meet all physical and mathematical requirements. In fact, choosing is the main challenge that one faces when dealing with inverse problems of image reconstruction. Another challenge is to compute the solution (numerically), when the latter has been defined. In medical imaging, the dimension of (i.e. the number of pixels or voxels) ranges from thousand to hundreds of thousands. It is therefore crucial to be able to handle large scale optimization problems.

Our research at the MIRG

An important part of our work, in the Medical Imaging Research Group, has been to implement optimization methods which utilize the most modern techniques developed by applied mathematicians. Combining results from the theory of convex duality (one of the major conceptual tools of modern optimization, which originated in the theory of Lagrange multipliers) and the power of quasi-Newtonian optimization methods leads to significant improvement of reconstruction algorithms (in terms of speed). Another aspect of our work concerns the analysis of the stability of the solution, that is, of the sensitivity of
to perturbations in the data vector . The stability analysis can be performed using the notion of sensitivity matrix. The first order perturbation of the reconstructed object is then given by

(3)

where is the reconstruction error corresponding to an error in the data. The matrix is referred to as sensitivity matrix. From Equation (3), various types of error analysis can be performed. In particular, the (spectral) norm of will provide an upper bound on the size of the reconstruction error. As a matter of fact we have, for all

and equality holds when is a singular vector corresponding to the highest singular value. The sensitivity matrix also enables us to see how statistical features are carried from the data domain to the object domain. For example, if the covariance matrix of the noise is available, we can obtain the covariance matrix of the corresponding perturbation in the object domain as follows:

Here, is the transpose of . Each diagonal element of is the variance of the corresponding pixel value. The square roots of these elements are the components of the standard deviation map.

Demonstration of the Reconstruction Process:

Primal Reconstruction Algorithm