A Stochastic estimation approach to emission tomography

Abstract

Emission Computed Tomography (ECT) is a medical imaging modality which is potentially useful in the study of human physiology and organ functions. In ECT, a physiologically active compound is radio-labelled and introduced into the body. These labelled compounds get distributed within the body in quantities proportional to the regional physiological activities. ECT imaging aims at obtaining a quantitative map of the spatial and temporal distributions of these compounds, from which information about the regional physiology can be obtained.

The physical processes involved in ECT imaging are intrinsically stochastic in nature and are described mathematically by appropriate stochastic models. With the emission process well modelled as a spatially independent Poisson process, the measurement process has also been modelled as an independent Poisson process. The observed data is considered a sample from this measurement stochastic process. In this case, the image reconstruction problem is essentially the estimation of the mean parameters of the emission process. This thesis deals with various stochastic estimation approaches for emission tomographic image reconstruction and consists of seven chapters. The beginning chapter is introductory in nature and discusses the potential use of ECT imaging in medical diagnostics and as a powerful medical research tool.

In Chapter 2, the physical and mathematical fundamentals, and instrumentation aspects of both Positron Emission Tomography (PET) and Single Photon Emission Computed Tomography (SPECT) are discussed. Assuming a Poisson model for the emission process, and taking into account the physical and geometrical aspects of the ECT imaging mode used, the probability function for observing a given measurement data (the likelihood function) has been derived. Maximisation of this function is the basis for the maximum likelihood estimation of emission densities.

In Chapter 3, a computer simulation of a PET system is described. In this system, an array of detectors in a ring geometry is used to collect the coincidence photon counts. With a square pixel decomposition of the object space, the pixel-detector probabilities are computed as the angle of intersection of rays from the center of a pixel to a pair of detectors. Using a mathematically generated phantom and the Hoffman brain phantom, which is a realistic simulation of the brain cross-section, measurement data are obtained on this simulated PET system. These sets of data have been used throughout this thesis for the evaluation of the reconstruction techniques developed.

Starting with the likelihood function for the measurement data, the conditions for the existence and uniqueness of a maximum likelihood solution are examined in Chapter 4. In this chapter, an iterative algorithm for computing the maximum likelihood estimate of emission densities is derived. This algorithm is based on the Expectation Maximisation (EM) formulation for ML parameter estimation from incomplete data. It is shown that the images reconstructed using this EM algorithm have much better quantitative accuracy and resolution than those produced by the Convolution Backprojection (CBP) method, which is a well-known deterministic reconstruction algorithm. However, the computation time required for each iteration of the EM algorithm is excessive and the iterations have very slow convergence. The EM algorithm has, therefore, been impractical for use in routine diagnostic and research applications. In this thesis, a major emphasis has been placed on developing faster algorithms for maximum likelihood estimation of emission densities.

Some of the currently available schemes for accelerating the convergence of the EM algorithm are implemented and evaluated using simulation studies. These include the over-relaxation method and the EM Search algorithm. As a major contribution in this thesis, a new class of fast maximum likelihood estimation algorithms for emission tomographic image reconstruction has been developed. These new algorithms are based on the acceleration of convergence by vector extrapolation. In these cyclic iterative algorithms, vector extrapolation techniques are integrated with the iterations in gradient-based MLE algorithms, with the objective of accelerating the convergence of the base iterations. The mathematical theory behind the minimal polynomial and reduced-rank vector extrapolation techniques, in the context of emission tomography, is presented. With the EM and EM Search algorithms in the base iterations, these extrapolation techniques are implemented and their convergence properties are studied. It is shown that with minimal additional computations, the proposed new approach substantially reduces the number of iterations required for obtaining the maximum likelihood estimate.

In Chapter 5, a three-dimensional maximum likelihood image reconstruction algorithm for a real SPECT system is developed. The SPECT system consists of a Dynascan 5C camera which counts photons on a 64 × 64 grid of square pixels. The measurements in SPECT are degraded by photon attenuation, scattering, and collimator blurring. These are modelled by two-dimensional, position-dependent response functions. By taking a series of point source acquisitions on the SPECT system with the Technetium-99m radionuclide, the distance-dependent primary response and the depth-dependent scatter response are modelled as two-dimensional Gaussian functions. These response functions are used in computing the pixel-detector probabilities. By using a simplifying assumption that the attenuation is uniform over a cross-section, a projector back-projector pair, which incorporates position-dependent blurring through pixel-detector probabilities, is developed. This projector back-projector pair is used in the implementation of the three-dimensional EM maximum likelihood reconstruction algorithm.

For testing the performance of the proposed three-dimensional SPECT MLE algorithm with real data, measurements have been taken using the SPECT system on a physical phantom and a human patient. The physical phantom used was the Jaszczak phantom, which consists of a ring of spheres and a number of rods in water medium. Technetium-99m activity was distributed in the water medium leaving the volume in the spheres and rods as cold spots. With the human patient, measurements are taken with Technetium-99m-labelled colloids introduced into the liver and spleen. With these real measurement data, the images reconstructed using the proposed three-dimensional EM ML algorithm have been of excellent quality. In our implementation with 64 × 64 × 64 pixel projection images reconstructed on a 64 × 64 × 64 image space, each EM iteration takes about 40 minutes on an IBM RISC 6000 computer. To obtain an acceptable image, about 40 such iterations are required. In order to accelerate the convergence of these iterations, the vector extrapolation techniques are implemented on this three-dimensional EM algorithm. This acceleration reduced the number of iterations required to about 15, which is a saving of about two-thirds of the computation effort of the non-accelerated EM ML algorithm.

The maximum likelihood estimate of emission densities is found to turn increasingly noisier as the ML iterations are continued and the estimates increase in likelihood. This has been known to be an example of the dimensional instability problem, which is common to maximum likelihood estimation of parameters of a continuous probability distribution function on the basis of a finite set of measurement data. An ad hoc approach for dealing with this problem is to stop the ML iterations prematurely before the images begin to turn noisy. In this approach, stopping the iterations based on a statistical criterion is examined in Chapter 6.

It can be assumed a priori that the radionuclide concentrations within subregions of common tissue types and common physiological activity are fairly homogeneous. The noisy appearance of the ML estimated images is an artifact and unacceptable. In cases where some such prior knowledge about the object is known, it can be incorporated in the reconstruction procedure to an advantage. The use of prior knowledge in the form of distribution functions, in a Bayesian framework, is the main topic of Chapter 6. The use of Gaussian, gamma, and Gibbs distributions as priors are examined. The Gaussian and gamma priors are pixel-based, in the sense that a mean estimate of the whole image is assumed to be known a priori and is used as the mean of the prior distribution. In practice, such a mean image is hardly available a priori. On the other hand, the Gibbsian prior achieves smoothness by incorporating spatial interactions between neighbouring pixels. In this case, abrupt changes in values for neighbouring pixels are penalised. A study on the use of various potential functions which model the spatial interactions is also done in this chapter.

The concluding chapter reviews the present status of the research in stochastic estimation methods for ECT image reconstruction. Future research in ECT imaging should be directed towards incorporating prior knowledge about the object in more fruitful ways. Incorporation of features from other modes of imaging, such as Magnetic Resonance Imaging (MRI) and X-ray CT, in emission tomographic image reconstruction will be of potential research interest in the future.