2D PET Image Reconstruction Using Robust L1 Estimation of the Gaussian Mixture Model

1Introduction

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 L1 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.

2Two-Dimensional PET Imaging

Fig. 1

(a) Simulated measurements for N=400 and K=2. (b) The corresponding sinogram.

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.

3Estimating the Gaussian Mixture Model

A Gaussian mixture model (Reynolds, 2015) is a weighted sum of K component Gaussian densities as given by the equation:

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.

3.2Estimation of Σ Using 1D Projections

In the two-dimensional setting, since Σ is a symmetric matrix,

Σ=Σ11Σ12Σ12Σ22,

we need to estimate three parameters: Σ11, Σ12 and Σ22.

Figure 2 depicts a single event and its corresponding LOR. Recall that the line is determined by two parameters, one of which is the slope tanφ, where φ is the angle between the line and the x-axis. If we rotate the coordinate system by φ−π2, in the new coordinate system the y-axis is parallel to the event. If we represent the Gaussian distribution as in the image, in this new setup it rotated by ψ=π2−φ. Since a rotated Gaussian distribution is again Gaussian, in this new coordinate system the event is parallel to the y-axis, and the Gaussian distribution has parameters

(8)

(Tμ,TΣTT),whereT=cosψ−sinψsinψcosψis the rotation matrix.

Fig. 2

Illustration of the rotation. The event (LOR) is parallel to the new y-axis.

Given the original event parameters, the coordinates of the intersection of the event and the new x-axis are (−lsinψ,0). The projection of the Gaussian distribution onto the new x-axis remains Gaussian, and the parameters of this one-dimensional distribution are components of the original two-dimensional Gaussian. From (8) we see that the expectation is given by

(Tμ)1=cosψμx−sinψμy,

and variance is

(9)

(TΣTT)11=cos2ψΣ11−2cosψsinψΣ12+sin2ψΣ22.

We can use these one-dimensional projections and their squared (Euclidian) distance from the mean, to estimate the variance in (9). This squared distance is equal to (cosψμx−sinψμy+lsinψ)2, which gives us a system of equations:

(10)

As=b,wheres=Σ11Σ12Σ22,A=cos2ψ1−2sinψ1cosψ1sin2ψ1cos2ψ2−2sinψ2cosψ2sin2ψ2⋮⋮⋮cos2ψN−2sinψ1cosψNsin2ψN,andb=(cosψ1μx−sinψ1μy+l1sinψ1)2⋮(cosψNμx−sinψNμy+lNsinψN)2.

This problem is clearly overdetermined, since there are (dozens of) thousands of measurements. There is no exact solution and the best approximation is found by solving

(11)

mins‖As−b‖.

Note that the solution will depend on the type of norm used in (11). Following classical methodology, the ordinary least squares (OLS) method gives the solution that minimizes the L2 norm. This solution can be found by solving the equivalent problem

ATAs=AT·b,i.e.s=(ATA)−1AT·b.

As we will show in Section 5, OLS performs poorly when we introduce more than one component, so we will need to use alternative methods. L1 minimization is a popular approach in engineering problems because it gives sparse solutions, but it had also been shown that it is less sensitive than the more traditional L2 minimization (Bektaş and Şişman, 2010). The main argument against L1 minimization would be computational complexity, which can be alleviated by using iterative methods. In this paper, the L1 minimization algorithm proposed in Sović Kržić and Seršić (2018) is used to solve a modified problem

(12)

As=k·b,

where k is a constant corrective scale factor. In a sufficiently large sample, one would have enough data points in each direction ϕ to obtain a one-dimensional variance estimate. In the absence of that, we are able to obtain a robust estimator for the parameters of the covariance matrix following the median absolute deviation (MAD) estimator of deviation σ (Rousseeuw and Croux, 1993), i.e.

k=(Φ(0.75))2≈1.48262,

where Φ denotes the distribution function of the standard normal random variable.

4EM-Like Algorithm

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:

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.

5Results

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 L1 and L2 minimization, for several covariance matrices with varying corresponding correlation coefficients:

Fig. 3

Left: Relative error in estimation of μ from a single measurement (N=1000 events). Right: Average relative error for increasing sample size, for three types of covariance.

5.2Two-Dimensional Estimation

For K=2 components we repeated the experiment as described in Section 4 for various combinations of types of vector s. We used hard classification, where we assign each line to at most one component. We experimented with various constraints, from simply assigning events to more likely components to assignations only when the probability of belonging is above a certain threshold.

For synthetic measurements from a variety of original (real) GMMs the algorithm proved robust regardless of the values of initial parameters, with estimation using L1 minimization and the scaling factor k correcting the bias. The L2 minimization method proved inefficient, since wrongly assigned events would cause unstable estimations and “breaks” in the algorithm.

Table 2 shows the percentage of correctly classified events for several combinations of covariance matrices, where a total of N=4000 events were simulated from two mixtures (N1=2500, N2=1500) for each combination. Events were classified with no probability thresholds, i.e. each event was assigned to the more likely component. We can conclude overall that the algorithm has a very good classification rate which, combined with the results for individual components, gives reason to expect accurate estimations for multiple components.

Table 2

Correct classification, K=2.

Cov. 1	Cov. 2	Comp. 1	Comp. 2	Total
Σ1	Σ2	87.63%	93.49%	89.83%
Σ2	Σ3	93.52%	92.13%	93.00%
Σ3	Σ1	91.26%	95.83%	92.98%

An illustration of the results for K=2, N=4000 is shown in Fig. 4, along with the corresponding classical FBP method.

Fig. 4

(a) Classical FBP reconstruction. (b) Proposed reconstruction using L1 minimization.

6Conclusion

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.