A Bayesian View of Reconstruction

The EMML iterative algorithm maximizes the likelihood function for the case in which the entries of the data vector y = ( y ₁ , , y _I) ^T are assumed to be samples of independent Poisson random variables with mean values ( P x) _i; here, P is an I by J matrix with nonnegative entries and x = ( x ₁ , , x _J) ^T is the vector of nonnegative parameters to be estimated. Equivalently, it minimizes the Kullback-Leibler distance KL( y , P( x)). This situation arises in single photon emission tomography, where the y _i are the number of photons counted at each detector i, x is the vectorized image to be reconstructed and its entries x _j are (proportional to) the radionuclide intensity levels at each voxel j. When the signal-to-noise ratio is low, which is almost always the case in medical applications, maximizing likelihood can lead to unacceptably noisy reconstructions, particularly when J is larger than I. One way to remedy this problem is simply to halt the EMML algorithm after a few iterations, to avoid over-fitting the x to the noisy data. A more mathematically sophisticated remedy is to employ a Bayesian approach and seek a maximum a posteriori (MAP) estimate of x.

In the Bayesian approach we view x