Statistical algorithms for large systems of linear equations in medical imaging

In tomography, activity recorded on cameras provides acquired counts from which the medical image of interest is reconstructed. These acquired data are only indirectly related to the source distribution within the body. This latter distribution can only be inferred indirectly (by mathematical analysis) from the camera projection data. Both transmission and emission tomography are well modelled by a description(1) based on a spatially inhomogeneous Poisson point process. If f(s) is the (unknown) distribution, then the distribution of recorded activity (on a camera) is expressible as g(t) = integral(2 is an element of chi) a(t \ s)f(s)ds, for t is an element of y or, in suitably discretized form, g = Af, with A = (a(ts)). Reconstruction (i.e. estimating f from observations on g) constitutes an inverse problem. Conditional probabilities a(t \ s) are determined by characteristics defining the acquisition system and the physics of photon transport and may be regarded as known. Issues for effective reconstruction are the high dimensionality and statistical variability. Both row and column dimensions of A are very large, and usually unequal. Algorithms utilising the sparse structure of A can be effective. Statistical variability arises because g is only realized in an empirical distribution function (e.g. histogram counts y). We present algorithms suitable for solving linear systems y = Ax, for x greater than or equal to 0 with a general sparse, nonnegative matrix A. Ill-conditioned systems will be common, and statistical variability requires attention. The linear system will be studied under different conditions. In consistent linear systems a feasible solution x greater than or equal to 0 exists. If the linear system is inconsistent, estimation of parameters of a General Linear Model (GLM) specifying a Poisson sampling distribution is of interest. The model for Poisson distributed counts y (the bin counts of data from g) then involves a linear predictor eta and link function specifying the relationship between eta and mu = Ey. We formulate our goal as minimizing a suitable objective function. Cases of interest include: NNLS: non-negative least squares seeks x greater than or equal to 0 which minimizes //y - Ax//(2). A unique minimizer exists in a linear system of full rank. Emission: the minimizer x greater than or equal to 0 of deviance 2[Sigma y(i) log y(i)/mu(i) - Sigma (y(i) - mu(i))] (1) where mu = eta (identity link) and eta = Ax; Transmission: the minimizer x greater than or equal to 0 of the above deviance, but with fitted value mu = mu(0) exp(-eta) (log-link), with eta = Ax, and known offset eta(0) = log mu(0). Note that for consistent (i.e. noise-free) transmission data, log transformation (y) over tilde = - log y/mu(0) produces the linear system (y) over tilde = Ax. Iterative statistical algorithms exist for solving the constrained linear and non-linear optimizations above. Algorithms include generalized iterative scaling, Shepp and Vardi's(2) application of the EM algorithm in emission tomography, and accelerated forms like our ordered subsets EM (OSEM). Properties in consistent (no noise) and Poisson noise conditions will be considered. We also introduce applications in some newly emerging applications in medical emission tomography of the flexible algorithmic technology described.