In positron emission tomography (PET), image reconstruction implies generating an image of a radiotracer’s concentration to estimate physiologic parameters for objects (volumes) of interest. Pairs of photons arise from emissions of annihilated positrons, and crystals placed in the scanner detect these pairs. When two are activated at the same time, the device records an event. Each possible pair of crystals is connected by a line or volume of response (LOR, VOR), depending on whether the scan is 2D or 3D. Assuming there are no contaminating physical effects or noise, the total number of coincidence events detected and the total amount of tracer contained in the tube are proportional. In the 2D case, we observe events along lines of response (LORs) joining two detector elements, lying within a specified imaging plane. The data are recorded as event histograms (sinograms or projected data) or as a list of recorded photon-pair events (list-mode data). Classic PET image reconstruction methodology is explained in more detail in e.g. Tong et al. (2010), Alessio and Kinahan (2006) or Reader and Zaidi (2007).

Modern image reconstruction methods are mostly based on maximum-likelihood expectation-maximization (MLEM) iterations. Maximum likelihood is used as the optimization criterion, combined with the expectation-maximization algorithm for finding its solution. To overcome the computational complexity and slow convergence of the MLEM, the ordered subsets expectation-maximization (OSEM) algorithm has been introduced (Hudson and Larkin, 1994). Since MLEM or OSEM estimation of pixels or voxels is usually noisy (Tong et al., 2010), use of post-filtering methods is necessary. On the one hand, image reconstruction from its projections is a mature research field with well known methods and proven results. On the other hand, known limitations made a challenge for a different approach presented in this paper.

In image segmentation, a number of algorithms based on model-based techniques utilizing prior knowledge have been proposed to model uncertainty, cf. Zhang et al. (2001, 1994). The Gaussian mixture model (GMM) is well-known and widely used in a variety of segmentation and classification problems (Friedman and Russell, 1997; Ralašić et al., 2018), as many observed quantities exhibit behaviour congruent with the model. A good overview of the application of GMMs and their generalizations to problems in image classification, image annotation and image retrieval can be found, e.g. in Tian (2018). However, in those problems the sample is from the GMM itself, whereas in this case the observed data is lower-dimensional (i.e. lines) and the points of origin are unknown.

In this paper, we propose a new robust method for reconstructing an object from a simulated PET image. The challenge of estimating Gaussian parameters from lower-dimensional data is solved by utilizing a novel L 1 minimizing algorithm. With it, our framework becomes similar to the standard Expectation-Maximization (EM) algorithm. Instead of values in individual pixels or voxels, we estimate parameters of the object’s model. This enables us to model one or more objects of interest, each as a mixture of one or more components, assuming that the entire object is a source of radiation with varying intensity. We focus on the two-dimensional case for clarity, but an extension to three dimensions would follow in a similar fashion. To the best of our knowledge, this is the first attempt to use parametric models in PET image reconstruction in this way. Parameters of Gaussian mixture models have not been estimated from lower-dimensional data due to computational complexity, which is circumvented here by the algorithm that efficiently finds global minimums. This paper is organized into six sections. Section 2 gives an overview of traditional methodology in 2D PET imaging. Section 3 describes the estimation of Gaussian mixture model parameters from PET data. Section 4 shows the implementation of iterative estimates in the EM-like algorithm. Finally, experimental results in Section 5 and discussion conclude the paper.

Data are collected along lines lying within a specific imaging plane. Traditionally, data are organized into sets of projections, integrals along all lines (LORs) for a fixed direction ϕ. The collection of all projections for 0 ⩽ ϕ < 2 π forms a two-dimensional function of the distance of the LOR from the origin, denoted by s, and angle ϕ. The line-integral transform of f ( x , y ) → p ( s , ϕ ) is called the X-ray transform (Natterer, 1986), which in 2D is the same as the Radon transform (Helgason, 1999). The path for any fixed point in the projection space resembles a sinusoidal curve. This is traditionally represented by a sinogram, a graph where these paths are drawn for all points of activity simultaneously. Figure 1 illustrates one synthetic measurement and the corresponding sinogram. Classical reconstruction methods, most notably the filtered backprojection (FBP), rely on the integral transform, cf. Natterer (1986), O’Sullivan and Pawitan (1993), Alessio and Kinahan (2006). They do not include or depend on any assumptions about the data. Assumptions such as form, volume and dependence could be utilized to obtain more precise estimates.

A Gaussian mixture model (Reynolds, 2015) is a weighted sum of K component Gaussian densities as given by the equation: (1) p ( x | τ k , μ k , Σ k ) = ∑ k = 1 K τ k f ( x | μ k , Σ k ) , where x is a d-dimensional observation, τ i ( i = 1 , … , K ) are the weights of each Gaussian component, and f ( x | μ k , Σ k ) are the Gaussian component densities. We assume that a given observation x is a realization from exactly one of the K Gaussian mixture components, and each component density is a d-variate Gaussian function of the form: (2) f ( x | μ k , Σ k ) = 1 ( 2 π ) d | Σ k | exp ( − 1 2 ( x − μ k ) T Σ k − 1 ( x − μ k ) ) . The set of probabilities { τ k } such that ∑ k = 1 K τ k = 1 defines the probabilities that x belongs to the corresponding Gaussian component. The complete Gaussian mixture model is parameterized by the mean vectors { μ k }, covariance matrices { Σ k } and mixture weights { τ k } for all Gaussian component densities.

Traditionally in this setting, the GMMs are used to model activity concentration in voxels (Layer et al., 2015) or pixel values (Nguyen and Wu, 2013). Those models do not take into account the spatial correlation between observations, which can be corrected by introducing Markov random fields (Layer et al., 2015; Nguyen and Wu, 2013). In other image segmentation problems, this can be addressed by modelling the mixture weights as probability vectors, thereby creating a spatially varying finite mixture model (SVFMM) (Xiong and Huang, 2015).

Here, we apply the GMM directly to the spatial data, focusing only on the locations and not the values at those locations. In other words, the location of the tracer, i.e. the particles that originate the gamma rays, is represented by the the K Gaussian components. The points x that are realizations from these components are latent (unobserved). Our observations, i.e. events, are lines through these points at random angles φ, however, using convenient properties of Gaussian distributions we will still be able to accurately estimate the parameters.

As mentioned earlier, Gaussian mixture models have naturally appeared in many signal processing (Yu and Sapiro, 2011; Léger et al., 2011) and biometrics problems (Reynolds, 2015; Reynolds et al., 2000). Due to the flexibility of the mixtures, GMMs are convenient and effective for modelling various types of signals, as well as image inverse and missing/noisy data problems.

Note that in our simulations we assume that K, the number of mixture components, is known. In practical applications it should also be estimated, which is a further (and not insignificant) problem. For mixture components that are sufficiently separated in space, such as in Fig. 1, the sinogram itself could give insight into the most adequate number of components, otherwise ideas suggested in, e.g. Leroux (1992), Chen (1995), Cheng (2005), should be explored.

For clarity, in this section we focus only on one Gaussian component. Assuming we know a set of N events whose points of origin come from a Gaussian distribution with parameters ( μ , Σ ), we can estimate those parameters.

In many statistical models, maximum likelihood parameter estimation can not be performed directly since most equations do not have an explicit closed solution. Several iterative methods have been developed to combat this, most notably the expectation-maximization (EM) algorithm. It had been proposed and used in many different circumstances before its proper introduction in Dempster et al. (1977). The name comes from its two-step setup, namely the expectation (E) and maximization (M) step which interchange until an acceptable solution is found. In the context of GMMs, and generally in emission tomography (Shepp and Vardi, 1982), the issue is that the true mixture membership is unknown (latent). The E step for given mixture parameters estimates probabilities { τ k } from (2) for each data point. The M step then (re)estimates the mixture parameters from the points assigned to that mixture with their corresponding probabilities. As already mentioned, our observations are lines and the true data points are latent, however, we can replicate the iterative steps using estimates from Section 3. Alternatively, this could also be considered a Lloyd-like algorithm, where the difference from the conventional Lloyd’s algorithm (Lloyd, 1982) is that we allow different distance functions of the form (4). The algorithm initializes parameters arbitrarily, and then alternates between the following steps: 1.

Compute class membership probabilities. For each event, compute the squared distance from each component mean. We distinguish between a hard classification where we assign the event to its nearest component, and a soft classification where membership probability is inversely proportional to the squared distance.

This approach allows for different variants, depending on the type of distance used (Tafro and Seršić, 2019).

In this paper, initial steps of the iterative algorithm use Euclidian distance in (4). Since later iterations improve the estimates, the distance gradually transforms into the Mahalanobis distance, i.e. (13) W = ( 1 − α ) I + α Σ ˆ − 1 , where α increases from 0 to 1. It is known that, for Σ given, the estimate obtained using the Mahalanobis distance is also the maximum likelihood estimate for μ. Therefore, in later iterations we obtain an MLE-like estimate of μ.

To evaluate the methodology, we experimented in the two-dimensional setting with K = 1 and K = 2 components.

First, for proof of concept we show that the method in Section 3 provides good estimates with both L 1 and L 2 minimization, for several covariance matrices with varying corresponding correlation coefficients: Σ 1 = 0.05 I , Σ 2 = 0.02 − 0.01 − 0.01 0.05 , Σ 3 = 0.01 0.02 0.02 0.05 . These matrices represent different shapes, i.e. different types of dependence, from independent ( Σ 1 ) to strongly dependent ( Σ 3 ) variables. Since the 2D covariance matrix is of the form Σ = σ x 2 ρ σ x σ y ρ σ x σ y σ y 2 , it is determined by three parameters – σ x , σ y and ρ ( = ρ x y ). This corresponds to the three-dimensional vector s = [ Σ 11 Σ 12 Σ 22 ] T = [ σ x 2 ρ σ x σ y σ y 2 ] T in (10). Instead of observing true and estimated Σ matrices we will calculate the error in estimation of s for each s that corresponds to matrices above: s 1 = 0.05 0 0.05 , s 2 = 0.02 − 0.01 0.05 , s 3 = 0.01 0.02 0.05 . Note that the corresponding correlation coefficients can be calculated from these vectors, and they are ρ 1 = 0, ρ 2 ≈ − 0.3, ρ 3 ≈ 0.9.

These results show proof of concept that it is possible to reconstruct data from latent mixture models using lower-dimensional observations and state-of-the-art computational techniques. Estimation from lower-dimensional measurements had not been thoroughly developed for Gaussian mixture models in general, and especially in this context. This paper shows that the main obstacle, slow computational speed, can be evaded by an alternative approach to minimization, which opens the door to various new procedures. Further work includes other metrics and membership calculation and utilizing traditional methods to obtain the number of mixing components and optimal initial values.

The fact that the reconstructed image is given by its parametric model (mean vectors { μ k }, covariance matrices { Σ k } and mixture weights { τ k }) has many benefits. This approach would allow for accurate reconstruction from significantly fewer measurements, thus decreasing patients’ exposure to radiation. It is virtually of infinite resolution, since the Gaussian components can be evaluated at each spatial point. The model is sparse: it consists from only a few parameters needed for successful object representation. Hence, future research will be oriented to compressed sensing approach (Rani et al., 2018; Ralašić et al., 2018, 2019) for reducing the number of projections, in this case reduced radiotracer’s concentration. Due to robustness of the proposed reconstruction method, a post-filtering step is not needed. This novel point of view aims to match the accuracy of iterative parametric models with the computational speed of the analytic approach by combining advantageous statistical models with state-of-the-art techniques from other image reconstruction problems.