mlem

Maximum Likelihood Expectation Maximization (MLEM) is an iterative algorithm used to estimate image intensities from projection data under a statistical model, typically Poisson noise. It is a specialization of the Expectation Maximization (EM) framework for emission tomography and was introduced by Shepp and Vardi in 1982 for PET imaging. The method seeks the image that maximizes the likelihood of obtaining the measured projections, given a forward model of the imaging system.

In the common formulation, the measured projection data y are modeled as Poisson with mean A x, where x is the image to reconstruct and A is the system matrix that describes the probability of detecting a photon from each voxel along each projection bin. At iteration k, the algorithm computes the forward projection p^{(k)} = A x^{(k)} and updates the image by x^{(k+1)}_j = x^{(k)}_j times the sum over all projections of a_{ij} y_i divided by p^{(k)}_i, all normalized by the corresponding column sums of A. In compact form, x^{(k+1)} = x^{(k)} ∘ [ A^T ( y ∘ (1 / (A x^{(k)})) ) ] / A^T 1, where ∘ denotes element-wise multiplication. The update guarantees non-decreasing likelihood under Poisson assumptions.

MLEM is commonly accelerated by methods such as Ordered Subsets EM (OS-EM), which uses subsets of data to speed convergence, though with potential oscillations. Extensions include MAP-EM, which adds prior information to regularize the solution and reduce noise. MLEM remains a foundational reconstruction technique in emission tomography, particularly PET and SPECT, where accurate system models and Poisson statistics are central.