Gaussian Mixture Estimation from Lower-Dimensional Data with Application to PET Imaging

Abstract

In positron emission tomography (PET), the original points of emission are unknown, and the scanners record pairs of photons emitting from those origins and creating lines of response (LORs) in random directions. This presents a latent variable problem, since at least one dimension of relevant information is lost. This can be solved by a statistical approach to image reconstruction—modeling the image as a Gaussian mixture model (GMM). This allows us to obtain a high-quality continuous model that is not computationally demanding and does not require postprocessing. In this paper, we propose a novel method of GMM estimation in the PET setting, directly from lines of response. This approach utilizes some well-known and convenient properties of the Gaussian distribution and the fact that the random slopes of the lines are independent from the points of origin. The expectation–maximization (EM) algorithm that is most commonly used to estimate GMMs in the traditional setting here is adapted to lower-dimensional data. The proposed estimation method is unbiased, and simulations and experiments show that accurate reconstruction on synthetic data is possible from relatively small samples.

1. Introduction

Positron emission tomography (PET) is an imaging technique that uses radioisotopes (i.e., radiotracers) and rings of crystal detectors. The injected tracer, designed to attach itself to the object or process of interest, decomposes and emits pairs of photons in opposite directions, creating a line that the scanner detects and records as an event when two crystals are activated simultaneously [1]. In a fully three-dimensional scanner, all possible pairs of crystals create tubes, i.e., volumes of response (VORs), which capture the events. Alternatively, in the two-dimensional case, pairs of crystals form lines of response (LORs), and the scanner detects events along a single plane, which are later used to reconstruct the 3D image.

Traditional reconstruction methods in tomography generally fall under two categories—analytic or iterative algorithms [2,3]. Analytic methods are based on the Fourier transform and its inverse, using the projection-slice theorem (sometimes called the central section theorem) [4]. The most used is the filtered back projection, also known as the Radon inverse transform [5,6]. This method filters the projections to reduce noise and artifacts and then backprojects them along different angles, and is still widely used in CT scan reconstructions [7]; however, the resulting image is of lower resolution compared with other methods. Similar traditional methods, now outperformed by others, are back projection filtering, where the filtering step occurs after the back projection, and the Fourier transform method that transforms the data into the frequency domain, which manipulates it there and then transforms it back to the spatial domain [8]. These methods are computationally simple and fast; however, they assume that the data are noise-free and they do not account for randomness or stochastic variability. This means that they require a lot of pre- and postprocessing in order to obtain good image quality. By definition, the inverse problem is ill-posed and is therefore very sensitive to small perturbations in the original data. Before reconstruction, raw data are corrected for attenuation and noise, and filters can be applied in both the spatial and frequency domains. Filtering is also commonly applied at the postprocessing stage, as well as various regularization techniques that enforce smoothness and penalize undesired properties [9].

The main problem in PET image reconstruction is attenuation, i.e., the loss of detection of true coincidence events due to photons being absorbed in the body, scattering out of the field of view or changing course after collision. There are many attenuation correction techniques using either additional CT scans or properties of the PET scanner itself. In clinical settings, the choice of method depends on resources like computational power and availability of combined PET/CT scanners [10,11]. In recorded data, one of the consequences of attenuation is wrongly recorded or false events. When one or both of the photons in a pair change course, it activates the wrong crystal, and the resulting measurement is a line that does not contain the real point of origin. This new event differs from the real unrecorded one, and depending on the size of the error, i.e., the angle by which the course had changed, it is considered as photon noncollinearity (small) or scattered coincidence (large difference). Finally, given that the scanner records events at discrete times, it is possible for two or more events to occur “at the same time”. If more than one pair of crystals is activated at the same time, the scanner can pair them incorrectly, resulting in random coincidences [12].

With increased computational powers, iterative models became widespread at the end of the 20th century [13,14,15,16,17,18] and are still most widely used, often in combination with analytical techniques [3]. They are usually based on the expectation–maximization (EM) algorithm [19], traditionally using maximum likelihood estimation (MLEM) [20]. In general, these techniques discretize the image domain into distinct pixels or voxels and iteratively improve their estimation of pixel/voxel values. In the maximum likelihood approach, the Poisson distribution is most commonly used to model the number of times each detector pair would record an event. The MLEM was later updated to accelerate convergence with methods that use only parts of the data at each step, e.g., the ordered subsets EM algorithm (OSEM) [21]. These methods are more resistant to noise, but they still require postfiltering for denoising, as well as some regularization methods to counteract issues arising from ill-conditioned problems [22].

It should be mentioned that recent progress in medical imaging has been achieved by utilizing deep learning and AI-based solutions [23,24,25,26]. This new approach is proving to be most useful in solving previously unfeasible problems that are ill-posed or lacking the mathematical framework, but in standard image reconstruction problems, traditional methods are still competitive [27]. Deep learning methods can be used to improve existing methodology, e.g., by denoising the data [28], but a full implementation of AI-based models is limited by evaluation and interpretability issues [29,30] and computational memory [31].

Standard modern algorithms in image reconstruction include some sort of stochastic model that utilizes prior or additional knowledge. The Gaussian mixture model (GMM) has been proven suitable for a variety of segmentation problems [32,33,34]. In PET reconstruction problems, the original data (i.e., points of origin) are unobserved, so GMMs are used to model activity at points of observation, with spatial dependence modeled by Markov random fields [35,36,37]. An alternative approach is to model the unknown points of origin using GMMs directly by utilizing nice statistical properties of the Gaussian distribution to overcome missing information [38,39]. This novel strategy yields a continuous model of infinite resolution that can be reconstructed directly and evaluated at any point in the imaging plane. Furthermore, it requires a much smaller sample than traditional methods, which means that scanning requires lower doses of radiation and less time in the scanner.

This paper presents a new method for estimating Gaussian parameters from lower-dimensional data. Using synthetic measurements, we demonstrate how this approach can be used to estimate a full Gaussian mixture model on PET data using an iterative EM-like algorithm. In Section 2, we show how to estimate parameters of a single distribution, and in Section 3, we expand this idea to simultaneously estimate parameters of a mixture of Gaussian distributions. Section 4 demonstrates the approach on synthetic data, and a discussion concludes the paper.

2. Gaussian Estimation from PET Data

Original PET data are a set of events recorded when the scanner registers a pair of photons activate two crystals simultaneously. These events can be represented as lines connecting the pairs of crystals, while the points of origin remain unknown (Figure 1).

Traditionally, these data would be stored as a histogram of events across all angles, called a sinogram. The original location of the tracer that emitted the photons would then be reconstructed using analytic or iterative methods, and corrected for measurement errors. In this paper, we assume that the points of origin follow a known distribution and estimate the parameters of that distribution directly from the lines.

To start, assume that the location of any point $\mathit{x}$, where a positron is annihilated and emits the pair of photons, follows a single Gaussian distribution with the parameters $\mathit{\mu }$ and $\mathbf{\Sigma }$. Note that the following calculations are performed in the 2D setting for simplicity, but analogous claims are true in the three-dimensional case. An event, i.e., line, in the 2D setting can be represented by two parameters, and possible representations are of course interchangeable. In this work, we will observe the line slope $k=tan\phi$, where $\phi \in \left[-\frac{\pi }{2},\frac{\pi }{2}\right]$, and an intercept l. Therefore, N events give us a series of N equations:

$$y=tan{\phi }_i·x+l_i,i=1,ߪ,N.$$

(1)

2.1. Estimation of $\mu$

It has been shown that the mean $\mathit{\mu }$ can be estimated in many ways, e.g., from intersections of lines [40] or directly from the sinogram [39]. For this part, we follow the methodology from [38], which shows that the mean is accurately estimated as the point that minimizes the total squared Euclidean distance from all events. If we denote ${\mathit{a}}_i={[tan{\phi }_i-1]}^T$, the location of this point can be calculated as follows:

$${\mathit{\mu }}^T=\left(\sum _{i=1}^N{l}_i{\mathit{m}}_i\right)·{\left(\sum _{i=1}^N\left(2{\mathit{M}}_i-{\mathit{N}}_i\right)\right)}^{-1},$$

(2)

where

$$\begin{array}{c}\left[\begin{array}{c}{\mathit{M}}_i\\ {\mathit{m}}_i\end{array}\right]={\mathit{B}}_i\left(\left[\begin{array}{cc}2\mathit{I}& {\mathit{a}}_i\\ {\mathit{a}}_i^T& 0\end{array}\right]\left[\begin{array}{c}2\mathit{I}\\ \mathbf{0}\end{array}\right]-\left[\begin{array}{cc}1& 0\\ 0& 1\\ 0& 0\end{array}\right]\right),\left[\begin{array}{c}{\mathit{N}}_i\\ \mathbf{0}\end{array}\right]={\mathit{B}}_i^T\left[\begin{array}{c}2\mathit{I}\\ \mathbf{0}\end{array}\right]-\tilde{\mathit{I}}\left[\begin{array}{cc}1& 0\\ 0& 1\\ 0& 0\end{array}\right],\\ {\mathit{B}}_i={\left({\left[\begin{array}{cc}2\mathit{I}& {\mathit{a}}_i\\ {\mathit{a}}_i^T& 0\end{array}\right]}^{-1}\right)}^T·\tilde{\mathit{I}}\mathrm{and}\tilde{\mathit{I}}=\left[\begin{array}{cc}\mathit{I}& \mathbf{0}\\ \mathbf{0}& 0\end{array}\right].\end{array}$$

Once the center has been estimated, we can also find the points on each line that are nearest to it. Consider the i-th event and the original unknown point of origin ${\mathit{x}}_i={[x_i,y_i]}^T$. Equation (1) can also be written as follows:

$$cos{\phi }_i(y-y_i)=sin{\phi }_i(x-x_i).$$

(3)

Any line perpendicular to the i-th line has the slope $-\frac{1}{tan{\phi }_i}=-\frac{cos{\phi }_i}{sin{\phi }_i}$, and especially one that goes through $\mathit{\mu }$ satisfies the following:

$$-sin{\phi }_i(y-{\mu }_y)=cos{\phi }_i(x-{\mu }_x).$$

(4)

For simplicity, denote $s_i=sin{\phi }_i$ and $c_i=cos{\phi }_i$. If we omit the index i in the coordinates for legibility, denote the point on line i nearest to the center $\mathit{\mu }$ as ${\mathit{y}}_i={[x_s,y_s]}^T$. This point satisfies (3) and (4), which can be written in matrix form:

$$\left[\begin{array}{cc}{s}_i& -c_i\\ c_i& s_i\end{array}\right]·\left[\begin{array}{c}{x}_s\\ y_s\end{array}\right]=\left[\begin{array}{cccc}{s}_i& -c_i& 0& 0\\ 0& 0& c_i& s_i\end{array}\right]·\left[\begin{array}{c}{x}_i\\ y_i\\ {\mu }_x\\ {\mu }_y\end{array}\right],$$

and gives the following solution:

$${\mathit{y}}_i=\left[\begin{array}{c}{x}_s\\ y_s\end{array}\right]=\left[\begin{array}{c}{s}_i^2{x}_i---s_i{c}_i{y}_i+c_i^2{\mu }_x+s_i{c}_i{\mu }_y\\ -s_i{c}_i{x}_i+c_i^2{y}_i+s_i{c}_i{\mu }_x+s_i^2{\mu }_y\end{array}\right].$$

(5)

It is straightforward to show from (5) that $E\left({\mathit{y}}_i\right)=\mathit{\mu }$ as well.

2.2. Estimation of $\mathsf{\Sigma }$

Statistical properties of the set of “nearest points” found in the previous section by (5) can be connected to the statistical properties of the original set of (unknown) points, which will enable us to estimate the original covariance matrix:

$$\mathbf{\Sigma }=\left[\begin{array}{cc}{\mathsf{\Sigma }}_{11}& {\mathsf{\Sigma }}_{12}\\ {\mathsf{\Sigma }}_{12}& {\mathsf{\Sigma }}_{22}\end{array}\right].$$

Consider the ith point of origin ${\mathit{x}}_i={[x_i,y_i]}^T$ from a Gaussian distribution with the parameters $\left(\mathit{\mu }\mathbf{\Sigma }\right)$ and its corresponding point nearest to the center, $\mathit{y}={[x_s,y_s]}^T$. We will now show how their covariance matrix $\mathit{C}=E((\mathit{y}-\mathit{\mu })·{(\mathit{y}-\mathit{\mu })}^T)$ relates to the original covariance matrix $\mathbf{\Sigma }$. First, assume for simplicity that $\mathit{\mu }={[0,0]}^T$, i.e., that the measurements have been centered. Equation (5) can now be written as follows:

$$\left[\begin{array}{c}{x}_s\\ y_s\end{array}\right]=\left[\begin{array}{c}{s}_i^2{x}_i-s_i{c}_i{y}_i\\ -s_i{c}_i{x}_i+c_i^2{y}_i\end{array}\right]=\left[\begin{array}{cc}{s}_i^2& -s_i{c}_i\\ -s_i{c}_i& c_i^2\end{array}\right]·\left[\begin{array}{c}{x}_i\\ y_i\end{array}\right].$$

Then we have the following:

$$\mathit{C}=E({[x_s,y_s]}^T·[x_s,y_s])=E\left(\left[\begin{array}{cc}{s}_{11}& s_{12}\\ s_{21}& s_{22}\end{array}\right]\right),$$

where

$$\begin{array}{cc}\hfill s_{11}& =s_i^4{x}_i^2-2s_i^3{c}_i{x}_i{y}_i+s_i^2{c}_i^2{y}_i^2\hfill \\ \hfill s_{12}& =-s_i^3{c}_i{x}_i^2+2s_i^2{c}_i^2{x}_i{y}_i-c_i^3{s}_i{y}_i^2=s_{21}\hfill \\ \hfill s_{22}& =s_i^2{c}_i^2{x}_i^2-2s_i{c}_i^3{x}_i{y}_i+c_i^4{y}_i^2.\hfill \end{array}$$

This can be expressed in matrix form as follows:

$$\left[\begin{array}{c}{s}_{11}\\ s_{12}\\ s_{22}\end{array}\right]=\left[\begin{array}{ccc}{s}_i^4& -2s_i^3{c}_i& s_i^2{c}_i^2\\ -s_i^3{c}_i& 2s_i^2{c}_i^2& -c_i^3{s}_i\\ s_i^2{c}_i^2& -2s_i{c}_i^3& c_i^4\end{array}\right]·\left[\begin{array}{c}{x}_i^2\\ x_i{y}_i\\ y_i^2\end{array}\right].$$

(6)

The key observation here is that the first matrix on the right-hand side is a function only of the angles at which the lines appear. These angles are random, uniformly distributed, and independent from the points of origin. Therefore, when taking the expectation of (6), it can be separated into a product of two expectations:

$$E\left(\left[\begin{array}{c}{s}_{11}\\ s_{12}\\ s_{22}\end{array}\right]\right)=E\left(\left[\begin{array}{ccc}{s}_i^4& -2s_i^3{c}_i& s_i^2{c}_i^2\\ -s_i^3{c}_i& 2s_i^2{c}_i^2& -c_i^3{s}_i\\ s_i^2{c}_i^2& -2s_i{c}_i^3& c_i^4\end{array}\right]\right)·E\left(\left[\begin{array}{c}{x}_i^2\\ x_i{y}_i\\ y_i^2\end{array}\right]\right).$$

(7)

Since the angle $\phi$ is uniformly distributed over $\left[-\frac{\pi }{2},\frac{\pi }{2}\right]$, we can calculate, e.g., as follows:

$$E\left(s_i^4\right)={\int }_{-\frac{\pi }{2}}^{\frac{\pi }{2}}{sin}^4\phi d\phi =\frac{3}{8},E(-2s_i^3{c}_i)=-2{\int }_{-\frac{\pi }{2}}^{\frac{\pi }{2}}{sin}^3\phi cos\phi d\phi =0,$$

and similar calculations for other expressions in (7) yield the following:

$$E\left(\left[\begin{array}{ccc}{s}_i^4& -2s_i^3{c}_i& s_i^2{c}_i^2\\ -s_i^3{c}_i& 2s_i^2{c}_i^2& -c_i^3{s}_i\\ s_i^2{c}_i^2& -2s_i{c}_i^3& c_i^4\end{array}\right]\right)=\left[\begin{array}{ccc}\frac{3}{8}& 0& \frac{1}{8}\\ 0& \frac{1}{4}& 0\\ \frac{1}{8}& 0& \frac{3}{8}\end{array}\right]=:\mathit{M}.$$

(8)

To estimate the original covariance matrix, we begin by estimating the covariance matrix of the (centered) nearest points:

$$\widehat{\mathit{C}}=\sum _{i=1}^N{\mathit{y}}_i{\mathit{y}}_i^T=\left[\begin{array}{cc}{c}_{11}& c_{12}\\ c_{21}& c_{22}\end{array}\right].$$

From (7) and (8), the following is derived:

$$\begin{array}{cc}\hfill \left[\begin{array}{c}{\sigma }_x^2\\ {\sigma }_{xy}\\ {\sigma }_y^2\end{array}\right]& =E\left(\left[\begin{array}{c}{x}_i^2\\ x_i{y}_i\\ y_i^2\end{array}\right]\right)={\mathit{M}}^{-1}·E\left(\left[\begin{array}{c}{s}_{11}\\ s_{12}\\ s_{22}\end{array}\right]\right)\hfill \\ \hfill & =\left[\begin{array}{ccc}3& 0& -1\\ 0& 4& 0\\ -1& 0& 3\end{array}\right]·E\left(\left[\begin{array}{c}{s}_{11}\\ s_{12}\\ s_{22}\end{array}\right]\right),\hfill \end{array}$$

(9)

which means the following:

$$\left[\begin{array}{c}{\widehat{\sigma }}_x^2\\ {\widehat{\sigma }}_{xy}\\ {\widehat{\sigma }}_y^2\end{array}\right]=\left[\begin{array}{ccc}3& 0& -1\\ 0& 4& 0\\ -1& 0& 3\end{array}\right]·\left[\begin{array}{c}{c}_{11}\\ c_{12}\\ c_{22}\end{array}\right].$$

3. Gaussian Mixtures in PET Imaging

3.1. Gaussian Mixture Model

The Gaussian mixture model [41] assumes that a d-dimensional observation $\mathit{x}$ is generated from a mixture of K Gaussian densities, as follows:

$$p\left(\mathit{x}\right|{\pi }_k,{\mathit{\mu }}_k,{\mathbf{\Sigma }}_k)=\sum _{k=1}^K{\pi }_k{f}_G\left(\mathit{x}\right|{\mathit{\mu }}_k,{\mathbf{\Sigma }}_k),$$

(10)

where ${\pi }_i(i=1,...,K)$ are the weights of each Gaussian component, and $f_G\left(\mathit{x}\right|{\mathit{\mu }}_k,{\mathbf{\Sigma }}_k)$ are the component densities. Each density is a d-variate Gaussian function of the following form:

$$f_G\left(\mathit{x}\right|{\mathit{\mu }}_k,{\mathbf{\Sigma }}_k)=\frac{1}{\sqrt{{\left(2\pi \right)}^d\left|{\mathbf{\Sigma }}_k\right|}}exp\left(-\frac{1}{2}{(\mathit{x}-{\mathit{\mu }}_k)}^T{\mathbf{\Sigma }}_k^{-1}(\mathit{x}-{\mathit{\mu }}_k)\right),$$

and we assume that $\mathit{x}$ originates from exactly one of the components (clusters). This is a model with a latent (hidden) variable that records component membership. If we denote this variable by Z, we can write ${\pi }_k=P(Z=k)$ as the probability that $\mathit{x}$ belongs to the k-th Gaussian component, with ${\sum }_{k=1}^K{\pi }_k=1$. The model parameters that need to be estimated are the mean vectors $\left\{{\mathit{\mu }}_k\right\}$, covariance matrices $\left\{{\mathbf{\Sigma }}_k\right\}$, and mixture weights $\left\{{\pi }_k\right\}$ for each Gaussian component.

Even in the traditional setting, this is a model with latent variables. Since the component membership of each observation is unknown, the maximum likelihood approach will not provide a closed-form solution. Several estimation methods have been proposed, and the most common one is the iterative expectation–maximization (EM) algorithm [19]. In the following section, we will adapt this approach to our estimation method from Section 2.

3.2. EM-like Algorithm

Broadly, the original EM algorithm iterates between two steps:

E:

Given a set of parameters $\Theta =\{{\pi }_k,{\mathit{\mu }}_k,{\mathbf{\Sigma }}_k\}$, calculate the probability that each observation ${\mathit{x}}^{\left(i\right)}$ originated from the k-th component:

$$h_k^{\left(i\right)}=P(Z_i=k|{\mathit{x}}^{\left(i\right)})=\frac{p\left({\mathit{x}}^{\left(i\right)}|{\mathit{\mu }}_k,{\mathbf{\Sigma }}_k\right){\pi }_k}{{\sum }_{j=1}^{K}p\left({\mathit{x}}^{\left(i\right)}|{\mathit{\mu }}_j,{\mathbf{\Sigma }}_j\right){\pi }_j}.$$

These probabilities are sometimes also called responsibilities, or posterior probabilities.

M:

Given the set of probabilities $\left\{h_k^{\left(i\right)}\right\}$, estimate the parameters of each component. This step is named after the maximum likelihood method that is used to estimate the parameters.

The algorithm can be initialized either from a random set of parameters $\theta$ or by randomly assigning the probabilities h. Note that, in our approach, the estimation is not derived directly from max likelihood given that we do not have the original data; nevertheless, we can adapt this algorithm to achieve satisfactory results.

3.2.1. Initialization

To ensure a “smart” initialization that will minimize the number of iterations later, we adapt the approach from [39]. This part combines k-means(-like) clustering with parameter estimation from the previous section (Matlab (R2022b) code for all experiments can be made available on request from the corresponding author).

• Each line is initially randomly assigned to exactly one component, while ensuring that each component has approximately the same number of lines.
• Mean vectors of each component are estimated as described in Section 2.
• Lines are reassigned to the component whose estimated mean vectors are nearest to them in terms of Euclidean distance. In this step, we still adhere to the so-called hard classification; i.e., each line is assigned to only one component.
  Steps 2 and 3 are repeated until the changes in mean vectors are sufficiently small.

Once the hard classification and reassignment are completed, the initial covariances are estimated once from all lines assigned to them following the methodology from Section 2. Initial weights $\left\{{\pi }_k\right\}$ are calculated as the relative size of the components, i.e., the proportion of all lines that are assigned to a specific component.

3.2.2. E Step

As in the traditional algorithm, given a set of parameters $\theta =\{{\pi }_k,{\mathit{\mu }}_k,{\mathbf{\Sigma }}_k\}$, this step calculates the probabilities $\left\{h_k^{\left(i\right)}\right\}$ that the i-th observation belongs to the k-th component. If we denote by ${\gamma }_i$ our i-th observation (i.e., some parametrization of the line), we need to calculate the following:

$$h_k^{\left(i\right)}=P(Z_i=k|{\gamma }_i)=\frac{p\left({\gamma }_i|{\mathit{\mu }}_k,{\mathbf{\Sigma }}_k\right){\pi }_k}{{\sum }_{j=1}^{K}p\left({\gamma }_i|{\mathit{\mu }}_j,{\mathbf{\Sigma }}_j\right){\pi }_j}.$$

(11)

The term $p\left({\gamma }_i|{\mathit{\mu }}_k,{\mathbf{\Sigma }}_k\right)$ is the integral of the density along the following line:

$$p\left({\gamma }_i|{\mathit{\mu }}_k,{\mathbf{\Sigma }}_k\right)={\int }_{{\gamma }_i}{f}_G\left(\mathit{x}\right|{\mathit{\mu }}_k,{\mathbf{\Sigma }}_k)\mathrm{d}{\gamma }_i,$$

where $f_G$ again denotes the (bivariate) Gaussian density. For a quicker computation, we can utilize the following nice properties of the Gaussian distribution:

• A Gaussian distribution retains properties when rotated.
• Marginal distributions of a Gaussian are again Gaussian.

If an event is at an angle $\phi$, rotating the coordinate system by $\psi =\frac{\pi }{2}-\phi$ makes it parallel to the y-axis and defined by its abscissa coordinate $x=-lsin\psi$. The coordinates of the center $\mathit{\mu }$ and covariance $\mathbf{\Sigma }$ are now $R\mathit{\mu }$ and $R\mathbf{\Sigma }{R}^T$, where R denotes the rotation matrix. Given that the line is now parallel to the y-axis, the integral along the line simply gives us the marginal distribution in the first coordinate (in our new coordinate system):

$$p\left({\gamma }_i|{\mathit{\mu }}_k,{\mathbf{\Sigma }}_k\right)=f_G(-lsin\psi |{\left(R_i{\mathit{\mu }}_k\right)}_1,{\left(R_i{\mathbf{\Sigma }}_k{R}_i^T\right)}_{11}),$$

(12)

where $R_i$ denotes the operator of rotation by ${\psi }_i$ and $f_G$ is now the univariate Gaussian density. In the two-dimensional setting, this can be calculated as follows:

$$p\left({\gamma }_i|{\mathit{\mu }}_k,{\mathbf{\Sigma }}_k\right)=f_G(-lsin\psi |cos{\psi }_i{\mu }_x-sin{\psi }_i{\mu }_y,{cos}^2{\psi }_i{\mathsf{\Sigma }}_{11}-2cos{\psi }_{i}sin{\psi }_i{\mathsf{\Sigma }}_{12}+{sin}^2{\psi }_i{\mathsf{\Sigma }}_{22}).$$

For more details, see [38], and Figure 2 illustrates the idea.

3.2.3. M-like Step

Given the estimated posterior probabilities $\left\{h_k^{\left(i\right)}\right\}$ and previous component centers $\left\{{\mathit{\mu }}_k\right\}$, we re-estimate the parameters. For each component, we have the following:

• Find the nearest points ${\mathit{y}}_i$ to the (previous iterations’) center for each line.
• Calculate the weighted mean and covariance of ${\mathit{y}}_i^{\prime }s$ with probabilities h as weights:

  $${\mathit{\mu }}_k=\frac{{\sum }_i{h}_k^{\left(i\right)}{\mathit{y}}_i}{{\sum }_i{h}_k^{\left(i\right)}},{\mathit{C}}_k=\frac{{\sum }_i{h}_k^{\left(i\right)}({\mathit{y}}_i-{\mathit{\mu }}_{\mathit{k}}){({\mathit{y}}_i-{\mathit{\mu }}_{\mathit{k}})}^T}{{\sum }_i{h}_k^{\left(i\right)}}.$$
• Calculate ${\mathbf{\Sigma }}_k$ from ${\mathit{C}}_k$ as in Section 2.
• Mixture weights ${\pi }_k$ are calculated, as in the original EM algorithm, as the proportion of all lines (events) assigned to each component: ${\displaystyle {\pi }_k=\frac{1}{N}\sum _{i=1}^N{h}_k^{\left(i\right)}}$.

The E and M steps are repeated until there are only small changes between iterations, specifically until there is virtually no change in estimated weight of each component. The following section evaluates our approach, both for a single component and in the mixture setting.

4. Experiments and Results

4.1. Single-Component Estimation

To demonstrate that the approach described in Section 2 correctly estimates the covariance matrix, we experimented with five covariances representing different shapes and types:

$$\begin{array}{c}{\mathbf{\Sigma }}_1=\left[\begin{array}{cc}0.0625& 0\\ 0& 0.0625\end{array}\right],{\mathbf{\Sigma }}_2=\left[\begin{array}{cc}0.04& 0.03\\ 0.03& 0.09\end{array}\right],{\mathbf{\Sigma }}_3=\left[\begin{array}{cc}0.04& 0.006\\ 0.006& 0.01\end{array}\right]\\ {\mathbf{\Sigma }}_4=\left[\begin{array}{cc}0.09& -0.0135\\ -0.0135& 0.09\end{array}\right],{\mathbf{\Sigma }}_5=\left[\begin{array}{cc}0.04& 0.004\\ 0.004& 0.04\end{array}\right]\end{array}$$

Note that the corresponding correlation coefficients are as follows:

$${\rho }_1=0,{\rho }_2=0.5,{\rho }_3=0.3,{\rho }_4=-0.15,{\rho }_5=0.1.$$

Estimation was performed $n=100$ times for each covariance type and for sample sizes N varying from $N=1000$ to $N=100,000$. For each iteration, the relative error was calculated as follows:

$$\mathrm{Rel}.\mathrm{error}=\frac{\left|\right|\mathbf{\Sigma }-{\mathit{S}\left|\right|}_F}{{\left|\right|\mathbf{\Sigma }\left|\right|}_F},$$

(13)

where $\left|\right|·{\left|\right|}_F$ denotes the Frobenius norm of the matrix. Figure 3 shows the average relative error for all covariance types and all values of N. We see that accurate reconstruction is achieved even from small (in the context of PET imaging) samples, starting with approximately $N=2000$, regardless of the covariance type.

4.2. Gaussian Mixture Estimation

To illustrate the effectiveness of the proposed algorithm, several experiments were conducted for two and three components. First, for $K=2$, the following parameters were set:

$$N=N_1+N_2,\frac{N_2}{N_1}=\frac{5}{7},{\mathit{\mu }}_1=\left[\begin{array}{c}0\\ 1\end{array}\right],{\mathit{\mu }}_2=\left[\begin{array}{c}1\\ 0\end{array}\right],{\mathbf{\Sigma }}_1=\left[\begin{array}{cc}0.0625& 0\\ 0& 0.0625\end{array}\right],{\mathbf{\Sigma }}_2=\left[\begin{array}{cc}0.04& 0.03\\ 0.03& 0.09\end{array}\right].$$

The total sample size N varied from 3000 ($N_1=1750,N_2=1250$) to 60,000 ($N_1$ = 35,000, $N_2$ = 25,000). Each experiment was repeated 100 times, and the results are given in Figure 4.

The estimated component weight was calculated as follows:

$${\widehat{N}}_k=\sum _{i=1}^N{h}_k^{\left(i\right)},$$

where $h_k^{\left(i\right)}$ was calculated as in (11). The criterion for stoppage was that each component weight changed by less than 10 between iterations. Note that the algorithm performed only a few (approx. 4–5) iterations each time before there was virtually no change in component weights. Figure 5 shows the ratio of relative and actual component sizes $\widehat{N_k}/N_k$, calculated as an average over all simulations for each sample size. Clearly, component weights are well estimated regardless of the total sample size.

Similar experiments were performed for $K=3$ components with varying sizes and shapes. The following parameters were used:

$$\begin{array}{c}N=N_1+N_2+N_3,\frac{N_2}{N_1}=\frac{5}{7},\frac{N_3}{N_1}=\frac{2}{7},{\mathit{\mu }}_1=\left[\begin{array}{c}0\\ 1\end{array}\right],{\mathit{\mu }}_2=\left[\begin{array}{c}1\\ 0\end{array}\right],{\mathit{\mu }}_3=\left[\begin{array}{c}1.25\\ -1\end{array}\right]\\ {\mathbf{\Sigma }}_1=\left[\begin{array}{cc}0.0625& 0\\ 0& 0.0625\end{array}\right],{\mathbf{\Sigma }}_2=\left[\begin{array}{cc}0.04& 0.03\\ 0.03& 0.09\end{array}\right],{\mathbf{\Sigma }}_3=\left[\begin{array}{cc}0.04& 0.006\\ 0.006& 0.01\end{array}\right].\end{array}$$

The total sample size N varied from 3500 ($N_1=1750,N_2=1250,N_3=500$) to 105,000 ($N_1$ = 52,500, $N_2$ = 37,500, $N_3$ = 15,000). Each experiment was repeated 100 times, and Figure 6 depicts the mean relative errors calculated as in (13) using the Euclidean norm for centers. The centers are well estimated even for very small samples, as well as the covariance of the smallest (third) component. A slightly larger sample is required for accurate covariance estimation for components with larger proportions of weight. Still, a total sample size of ≈100,000 is still significantly smaller than the traditional PET measurement size.

As in the two-component case, Figure 7 shows that component sizes (weights) were correctly estimated even from small sample sizes. It can also be noted that the number of iterations increased approximately from 10 to 20 with growing sample sizes, with no simulation exceeding 22 iterations per component. The results presented here are representative of many simulations for different types and combinations of covariance matrices, which were omitted here for brevity. Some additional similar results can be found in Appendix A for comparison.

Finally, Figure 8 illustrates the similarity between the original Gaussian model, i.e., the unknown components, and the proposed estimated reconstruction.

4.3. Noise Resistance

To demonstrate how this approach would behave in a more realistic setting, two types of noise have been simulated. First, a noise effect was added to some points of origin. Points were generated from three mixture components as described in Section 4.2, and a certain number of them were offset by error:

$$\tilde{\mathit{x}}=\mathit{x}+\mathit{\epsilon },\mathit{\epsilon }\sim \mathcal{N}\left(\mathbf{0},\sigma \mathit{I}\right).$$

Recall that when a photon slightly changes course, the resulting measurement is an event that is “close” to the original point, but does not pass through it. Since line directions are random and independent from the points themselves, this type of noise mimics the occurrence of photon noncollinearity—an event is recorded that does not contain the real point of origin, but is close to it. The parameters $\sigma =0.005$ and $\sigma =0.01$ were selected to represent approximately 10% and 20% variability of the original data, respectively, and a randomly selected subset of 10% of data was “noisified”. The results are shown in Figure 9. Compared with Figure 6, we see that the error decays more slowly, but satisfactory precision is still achieved.

To assess how the estimation behaves with different quantities of photon noncollinearity, an error with $\sigma =0.005$ was added to an increasing percentage of data (from 1% to 20%). This experiment was performed on a sample size N = 105,000, corresponding to component sizes of $N_1$ = 52,500, $N_2$ = 37,500, $N_3$ = 15,000, and the results are shown in Figure 10. The precision deteriorates with added noise, but even with 20% of noisy data, the error terms are still below 5%.

To represent random coincidences, additional false measurements were added to the original measurements. The points of origin were generated randomly within the scanner field of view and added to the real measurements. The algorithm was slightly modified to accommodate potential outliers. One of the oldest and simplest ideas is the so-called “three-sigma rule”, where measurements are classified as outliers if they are more than three standard deviations away from the mean [42]. Here, estimation was first performed on the full dataset. Then, each line was projected to one dimension as described in Section 3.2.2, along with the component means and variances (see Figure 2). If the line was more than three standard deviations away from all three components, it was classified as an outlier. Finally, the parameters were estimated again with the outliers omitted. The experiment was performed with the parameters as in Section 4.2 and the total sample size ranging from $N=3500$ to N = 105,000, with 2% random events added. The results are shown in Figure 11. More advanced and precise outlier detection methods will be incorporated when the algorithm is applied to real data, but this demonstrates that the method on its own is robust for small amounts of random coincidences.

This section shows that the algorithm is resistant to a relatively small amount of noise and false measurements. Therefore, less pre- and postprocessing would be required to achieve accurate reconstruction, which is important given that excessive regularization risks losing important features of the data.

5. Conclusions

In this paper, we have presented a new way of estimating Gaussian parameters from lower-dimensional observations in the PET context. The EM algorithm is adapted for GMM estimation from projections, with experiments showing that accurate reconstruction is possible from relatively small samples. The method relies on the geometric properties of the Gaussian distribution, and simulations demonstrate that the obtained unbiased estimators quickly converge to real parameter values. Compared with traditional methods that estimate individual pixels (or voxels) and give quantized values, this approach results in a continuous model of the PET image.

The limitation of this work is that it was performed on synthetic measurements, and application to real data is the next step in ongoing research. This includes incorporating scanner properties [43], as well as some issues that traditionally arise when dealing with mixture models, such as determining the proper number of components [44]. As it stands, this manuscript serves as proof of concept and of the feasibility of the approach.

To the best of our knowledge, Gaussian mixture models have not yet been used in modeling PET images in this way, and estimating model parameters from data of this type is a new and unresearched topic. Given the significance and usefulness of Gaussian models in image reconstruction problems, these results are encouraging and imply efficient reconstruction from fewer measurements that are currently used. This would, of course, mean that accurate imaging could be achieved while exposing patients to significantly lower doses of radiation. Additionally, the method could be applied to an array of other mixture models where the distributions have similar properties to the Gaussian, thus expanding the types of images that can be reconstructed in this manner.