Accurate PET Reconstruction from Reduced Set of Measurements based on GMM

Tomislav Matulić*, Damir Seršić

* Corresponding author

Department of Electronic Systems and Information Processing, University of
Zagreb Faculty of Electrical Engineering and Computing, Unska 3, 10000 Zagreb,
Croatia

###### Abstract

In this paper, we provide a novel method for the estimation of unknown
parameters of the Gaussian Mixture Model (GMM) in Positron Emission Tomography
(PET). A vast majority of PET imaging methods are based on reconstruction
model that is defined by values on some pixel/voxel grid. Instead, we propose
a continuous parametric GMM model. Usually, Expectation-Maximization (EM)
iterations are used to obtain the GMM model parameters from some set of point-
wise measurements. The challenge of PET reconstruction is that the measurement
is represented by the so called lines of response (LoR), instead of points.
The goal is to estimate the unknown parameters of the Gaussian mixture
directly from a relatively small set of LoR-s. Estimation of unknown
parameters relies on two facts: the marginal distribution theorem of the
multivariate normal distribution; and the properties of the marginal
distribution of LoR-s. We propose an iterative algorithm that resembles the
maximum-likelihood method to determine the unknown parameters. Results show
that the estimated parameters follow the correct ones with a great accuracy.
The result is promising, since the high-quality parametric reconstruction
model can be obtained from lower dose measurements, and is directly suitable
for further processing.

Keywords: Positron emission tomography, Gaussian mixture model, Method of
moment, Iterative algorithm, Reduced measurements

##  1 Introduction

Positron emission tomography is a medical imaging modality that measures
metabolic activity of the observed tissue. It is based on electron-positron
annihilation. The annihilation happens due to the radioactive $\beta^{+}$
decay of the radioactive tracer that is injected into the measured tissue. The
result of the annihilation is two high-energy photons (511 keV each) that
travel along the same line but in opposite directions. If two photons hit
detectors in a short-coincidence time-window, then the data is being stored
and it is considered as a valid measurement event. Virtual path between the
two detectors is called line of response (LoR). The main challenge in medical
imaging is to reconstruct image from such measurement data. In PET imaging
systems, a list of all coincidence events corresponds to a measurement data
set.[31][3]

Analytical reconstruction algorithms [23][4][15][21] (filtered back-
projection, back-projection filtering) are rarely used nowadays. They are fast
and computationally inexpensive, but often do not include information about
the real PET imaging systems and, therefore, produce images of lesser quality.
Iterative methods are the golden standard in PET image reconstruction.[33]
[37] The core of iterative algorithms is the system matrix, which represents
the PET imaging system.[13][17] Maximum-likelihood expectation-maximization
(MLEM) algorithm [29] is the most widely used iterative method. A major
drawback of the MLEM algorithm is its slow convergence. To accelerate the
convergence, the ordered subset expectation-maximization (OSEM) was introduced
in [11]. It significantly reduces the time needed for image reconstruction
(more subsets = faster convergence), but it does not guarantee the maximum
likelihood solution. More recently, deep neural networks are exploited in
modeling of the PET imaging systems [9][2][32] and for image reconstruction
[7][35][8]. Pay attention that all mentioned algorithms result in a spatially-
discrete reconstruction model, namely, a grid of pixels or voxels.

The Gaussian mixture model [27] (GMM) is a well-investigated approach in a
variety of classification and segmentation problems.[6] [26] The application
of the GMM can be found in many problems in biometrics, signal processing, and
speech modeling. [27] [28] [20] [36] [18] [19] [30] [25] Two different ways
are often used for the estimation of parameters of GMM: the expectation-
maximization method [10][22] and the method of moments [12][14]. Method of
moments relies on tensor moments of higher order to estimate the unknown
parameters of the GMM. The theory of such a method is well known and
thoroughly investigated. The major drawbacks are:

  * •

high-order moments are computationally inefficient [24], since the $k$-th
moment of an $n$-dimensional random variable is a tensor of size $n^{k}$;

  * •

the method of moments leads to a multivariate polynomial system. Statistically
meaningful solutions of a such system often do not exist or are not
unique.[34]

A minimum order of the moment needed for the estimation of unknown parameters
increases with the number of components in the Gaussian mixture. Hence,
Expectation-Maximization (EM) algorithm is studied more intensely than the
method of moments. Although the EM algorithm converges slowly [1], its
drawbacks are far less concerning. Notice that the obtained GMM model is
spatially-continuous, and virtually of infinite resolution. It motivates our
research.

In this paper, we present a new method for estimation of the GMM parameters in
2D PET imaging. We can state this problem in another way - how to estimate the
unknown GMM parameters from the lines of response that originate from some
point that follows Gaussian distribution, but the line is fired under an
arbitrary angle? All the methods mentioned in the previous paragraph deal with
point-wise samples from the GMM distribution itself, while we deal with the
lines that originate from those point-wise sources. To the best of our
knowledge, the mentioned issue is unexplored or under-explored [30], [16], and
this paper gives a novel insight into it.

In this work, the problem and its solution is set in the projection domain.
Each point in the image gives a sine function in the projection domain. Thus,
the projection domain is often called sinogram: a collection of sines. Vice-
versa, a single point in the projection domain corresponds to a line of
response in the image. Notice that centers of each Gaussian component are
points in the image domain, too.

In our work, we exploit the following facts:

  1. 1.

the mean vector (i. e. the center $\bm{\mu}$) corresponds to the sine function
in the projection domain;

  2. 2.

projection under an arbitrary angle of a bivariate normal distribution is a
univariate normal distribution. This is a special case of the marginal
distribution theorem of the multivariate normal distribution;

  3. 3.

expressions for calculating higher moments of the Gaussians are well known,
and have shown to be useful for estimation of the GMM parameters.

In Section 3 we exploit the mentioned facts, and using some tricks arrange it
to get a solvable system of equations. Then, a two-step iterative process that
resembles the EM algorithm is exploited for determining all of the unknown
parameters of the GMM. The first step updates membership probabilities between
the lines of response and the GMM components. The second step reevaluates the
parameters of each Gaussian component according to the updated memberships.

Resulting GMM model has several advantages over the usual pixel or voxel based
approach. Essentially, it is an infinite resolution continuous model obtained
from a reduced set of measurements. Widespread digital processing methods are
based on difference equations that approximate the underlying differential
equations. In our case, we can apply them directly on the GMM model without
approximations.

This paper is divided into six sections. In Section 2 we explain a 2D PET
imaging system. We present an estimation method of the mean vector and the
covariance matrix for one component of the GMM in Section 3. In Section 4 we
show the proposed iterative, EM-like algorithm for the estimation of the
unknown parameters of the whole GMM model. The results are given in Section 5.
Finally, in Section 6 we conclude the paper.

##  2 Two-Dimensional PET Imaging

The mathematical background of PET imaging is based on the Radon transform.
The expression for the Radon transform is

$p(s,\theta)=\displaystyle\int\displaylimits_{-\infty}^{+\infty}\int\displaylimits_{-\infty}^{+\infty}f(x,y)\,\delta(x\cos(\theta)+y\sin(\theta)-s)\,dxdy,$
(1)

where function $f(x,y)$ represents metabolic activity of the object that is
being scanned by the PET system, and $p(s,\theta)$ models projection at the
angle $\theta$ \- an angle between the projection line and the x-axis. Value
$p(s,\theta)$ corresponds to the line integral that is perpendicular to the
projection line, and its distance from the origin is $s$. Described setup is
illustrated in Fig. 1.

Figure 1: Radon transform

In our experiment, function $f(x,y)$ is regarded as the probability density
function (PDF) of some GMM. Each sample of the GMM is a position of electron-
positron annihilation. The samples fire lines of response under an uniformly
distributed random angle. Actual positions of annihilation (samples) are
latent (hidden). In Fig. 2(a), we can see a simulation of the PET measurement.
We generated $N=2000$ lines firing in both directions from some hidden point-
wise samples that follow the GMM distribution. Contours are used to indicate
the two Gaussian components, and only a subset of LoR-s is displayed for
better visibility.

(a) 2D PET measurement

(b) LoRs represented as points $(s,\varphi)$. Solid curves depict means, while
dashed curves depict variances, both functions of angle $\varphi$.

Figure 2: 2D PET measurement

A line of response can be represented as a pair $(s,\varphi)$, where
$s\in\mathbb{R}$ is the oriented distance from the origin (see $s$-axis in
Fig. 1 or Fig. 3), and $\varphi\in\left[-\frac{\pi}{2},\frac{\pi}{2}\right]$
is the angle between the line and the x-axis. LoR-s correspond to points
$(s,\varphi)$ in the projection domain. The integral line and projection line
are perpendicular, i.e. the angles $\theta=\varphi\pm\frac{\pi}{2}$.
Measurement is illustrated as a set of points $(s,\varphi)$ in Fig. 2(b), or
as a set of lines of the response in Fig 2(a). The mean vector of the Gaussian
component traces a sine function in the projection domain, shown as a solid
curve. Dashed curves stand for $\pm 3\sigma$-neighborhood around the mean.

In the next section, we focus our attention on how to estimate the mean vector
$\bm{\mu}$ and covariance matrix $\bm{\Sigma}$ of a single component in the
mixture. Estimation of the mean vector is done by fitting _the best_ sinusoid
on pairs $(s,\varphi)$. Covariance matrix $\bm{\Sigma}$ is estimated from the
moments of the probability density function.

##  3 Estimation of Parameters of the Gaussian Mixture Model

The Gaussian mixture model with $K$ components is given by

$g(\bm{x};\tau_{k},\bm{\mu}_{k},\bm{\Sigma}_{k})=\displaystyle\sum\displaylimits_{k=1}^{K}\tau_{k}f_{G}(\bm{x};\bm{\mu}_{k},\bm{\Sigma}_{k}),$
(2)

where $\bm{\mu}_{k}$ is the mean vector, $\bm{\Sigma}_{k}$ is the covariance
matrix and $\tau_{k}$ is the weight, all associated to the $k$-th component.
In general, function $f_{G}(\bm{x};\bm{\mu}_{k},\bm{\Sigma}_{k})$ is a
$d$-variate normal (Gaussian) distribution:

$f_{G}(\bm{x};\bm{\mu}_{k},\bm{\Sigma}_{k})=\frac{1}{\sqrt{(2\pi)^{d}|\Sigma_{k}|}}\exp(-\frac{1}{2}\left(\bm{x}-\bm{\mu}_{k}\right)^{\intercal}\bm{\Sigma}_{k}^{-1}\left(\bm{x}-\bm{\mu}_{k}\right)).$
(3)

Weights $\tau_{k}>0$ should suffice the condition
$\sum\displaylimits_{k=1}^{K}\tau_{k}=1$, since $g$ defines the probability
density.

Now, we focus only on one component in the mixture. Since we deal with 2D PET
imaging, we have bivariate ($d=2$) normal distribution. Also, we assume that a
number of components $K$ in the mixture is known in advance. Estimation is
possible, but it is beyond the scope of this paper.

###  3.1 Estimation of Mean Vector

First, we take a look of the univariate normal distribution. Let
$(x_{i})_{i=1}^{N}$ be the samples of some univariate normal distribution.
Then, the least square problem for

$L_{uni}(c;x_{i})=\displaystyle\sum_{n=1}^{N}(c-x_{i})^{2}$ (4)

has the solution

$c_{LS}=\frac{1}{N}\displaystyle\sum_{n=1}^{N}x_{i}\quad=\hat{\mu}.$ (5)

Obviously, least square solution $c_{LS}$ corresponds to the known mean value
estimate of the univariate normal distribution $\hat{\mu}$.

(a) Estimation of $\mu_{y}$

(b) Estimation of $\mu_{x}$

Figure 3: Estimation of mean vector

The point $\bm{\mu}_{k}$ (i.e. the mean vector) corresponds to a sinusoidal
function in the projection domain:

$m(\varphi;A,B)=A\sin(\varphi)+B\cos(\varphi),$ (6)

where $A\in\mathbb{R}$ and $B\in\mathbb{R}$ are unknown parameters yet to be
determined. As depicted in Fig. 3(a), we have:

$\mu_{y}^{(k)}=m(\varphi=0;A,B)=B,$

and Fig. 3(b) gives:

$\mu_{x}^{(k)}=m(\varphi=-\frac{\pi}{2};A,B)=-A.$

Notice that both $\varphi=-\frac{\pi}{2}$ and $\varphi=\frac{\pi}{2}$
correspond to LoR-s parallel to $y$-axis, but differ in orientation of
$s$-axis. In Fig. 3(b), we denoted by $s$ direction when
$\varphi=-\frac{\pi}{2}$, and by $s^{\prime}$ when $\varphi=\frac{\pi}{2}$. We
can rewrite Eq. 6 as

$m(\varphi;\mu_{x}^{(k)},\mu_{y}^{(k)})=-\mu_{x}^{(k)}\sin(\varphi)+\mu_{y}^{(k)}\cos(\varphi),$
(7)

where $\bm{\mu}_{k}=\bigl{(}\begin{smallmatrix}\mu_{x}^{(k)}\\\
\mu_{y}^{(k)}\end{smallmatrix}\bigr{)}$ is the mean vector of $k$-th Gaussian
component.

We know that the projection under an arbitrary angle of a bivariate normal
distribution is a univariate normal distribution. Combining Eq. 7 and Eq. 4
results in a weighted least square problem:

$\begin{gathered}L(\mu_{x}^{(k)},\mu_{y}^{(k)};s_{i},\varphi_{i})=\displaystyle\sum_{i=1}^{N}p_{ij}(m(\varphi_{i};\mu_{x}^{(k)},\mu_{y}^{(k)})-s_{i})^{2}=\\\
\displaystyle\sum_{i=1}^{N}p_{ik}(-\mu_{x}^{(k)}\sin(\varphi_{i})+\mu_{y}^{(k)}\cos(\varphi_{i})-s_{i})^{2},\end{gathered}$
(8)

where $(s_{i},\varphi_{i})$ represent $i$-th LoR, $N$ is the total number of
LoR-s, and weights $p_{ik}$ are membership probabilities between $i$-th LoR
and $k$-th Gaussian component. The calculation of the membership probabilities
is given in Section 4.

To solve the least square problem, we set the gradient of the $L$ function to
zero. We get $2\times 2$ linear system $\bm{M}_{k}\bm{\mu}_{k}=\bm{b}_{k}$,
where:

$\begin{gathered}\bm{M}_{k}=\begin{pmatrix}-\displaystyle\sum_{i=1}^{N}p_{ik}\sin^{2}(\varphi_{i})&\displaystyle\sum_{i=1}^{N}p_{ik}\sin(\varphi_{i})\cos(\varphi_{i})\\\
-\displaystyle\sum_{i=1}^{N}p_{ik}\sin(\varphi_{i})\cos(\varphi_{i})&\displaystyle\sum_{i=1}^{N}p_{ik}\cos^{2}(\varphi_{i})\end{pmatrix}\\\
\bm{b}_{k}=\begin{pmatrix}\displaystyle\sum_{i=1}^{N}p_{ik}s_{i}\sin(\varphi_{i})\\\
\displaystyle\sum_{i=1}^{N}p_{ik}s_{i}\cos(\varphi_{i})\end{pmatrix}.\\\
\end{gathered}$ (9)

As in (5), the solution is estimated mean
$\bm{\hat{}}{\bm{\mu}_{k}}=\bigl{(}\begin{smallmatrix}\hat{\mu}_{x}^{(k)}\\\
\hat{\mu}_{y}^{(k)}\end{smallmatrix}\bigr{)}\mathrel{\mathop{\mathchar
58\relax}}$

$\bm{\hat{}}{\bm{\mu}}_{k}=\bm{M}^{-1}_{k}\bm{b}_{k}.$ (10)

###  3.2 Estimation of Covariance Matrix

Once we have the mean $\hat{\bm{\mu}}_{k}$, we need to estimate the covariance
matrix $\bm{\Sigma}_{k}$. Still, we focus only on one component in the
mixture, and our approach is based on the variance of the projected univariate
normal distribution.

We observe a centered ($\bm{\mu}=\bm{0}$) bivariate normal distribution whose
principal axes are the $x$ and $y$, without loss of generality. For any
bivariate normal distribution, we simply define a new coordinate system where
the new origin is the mean vector of the distribution (translation), and axes
that are the eigenvectors of the covariance matrix (rotation).

A centered bivariate normal distribution whose principal axes are the $x$ and
$y$ can be expressed as:

$f_{Gc}(x,y)=\frac{1}{2\pi\sigma_{1}\sigma_{2}}\exp(-(\frac{x^{2}}{2\sigma_{1}^{2}}+\frac{y^{2}}{2\sigma_{2}^{2}})).$

We want to find the distribution under some angle $\varphi$ as seen in Fig.
4(a). We must calculate integral:

$f_{p_{0_{\varphi}}}(s)=\int_{\gamma_{\varphi}}f_{Gc}(x,y)d\gamma_{\varphi},$
(11)

where $\gamma_{\varphi}(t)=(t,kt+l)$, $k=\tan(\varphi)$ and
$l=\frac{s}{\cos(\varphi)}$. Mapping $\gamma_{\varphi}$ parameterizes the set
of all lines under the angle $\varphi$. Fig. 4(b) explains the connection
between y-intercept ($l$) and $s$-axis. After integration we get

$f_{p_{0_{\varphi}}}(s)=\frac{1}{\sqrt{2\pi}\sigma_{p_{0}}(\varphi)}\exp(-\frac{s^{2}}{2\sigma^{2}_{p_{0}}(\varphi)}),$
(12)

where
$\sigma_{p_{0}}^{2}(\varphi)={\sigma_{1}^{2}\sin^{2}(\varphi)+\sigma_{2}^{2}\cos^{2}(\varphi)}$.
A detailed calculation can be found in A. Function $\sigma_{p_{0}}(\varphi)$
can be interpreted as the distance from the origin to one of the tangent lines
on the ellipse $\frac{x^{2}}{\sigma_{1}^{2}}+\frac{y^{2}}{\sigma_{2}^{2}}=1$,
under angle $\varphi$, as shown in Fig. 4(c). If we substitute $\varphi$ with
$\varphi-\varphi_{0}$ in $\sigma_{p_{0}}(\varphi)$, we get more general
expression for the variance of some rotated Gaussian (i.e. ellipse). It
corresponds to the rotation of the entire coordinate system by angle
$\varphi_{0}$.

(a) Line integral (projection) under angle $\varphi$

(b) Relation between $y$-intercept and $s$ axis

(c) Ellipse, $\sigma_{p_{0}}(\varphi)$

(d) Rotated ellipse, $\sigma_{p}(\varphi)$

Figure 4: Projection of bivariate normal distribution and elliptical
illustration of its variance

The expression is now:

$\sigma^{2}_{p}(\varphi)={\sigma_{1}^{2}\sin^{2}(\varphi-\varphi_{0})+\sigma_{2}^{2}\cos^{2}(\varphi-\varphi_{0})},$

as shown in Fig. 4(d). The probability density function of the projected and
centered ($\bm{\mu}=\bm{0}$) bivariate normal distribution under angle
$\varphi_{0}$ is:

$f_{p_{\varphi}}(s;\sigma_{1},\sigma_{2},\varphi_{0})=\frac{1}{\sqrt{2\pi}\sigma_{p}(\varphi)}\exp(-\frac{s^{2}}{2\sigma^{2}_{p}(\varphi)}).$
(13)

For $k$-th component, additionally we need to translate the coordinate system
by mean vector $\bm{\mu}_{k}$. In the projection space it is:

$(s_{c}^{(k)},\varphi)=(s-m(\varphi;\mu_{x}^{(k)},\mu_{y}^{(k)}),\varphi).$
(14)

A point $(s_{c}^{(k)},\varphi)$ is a LoR in the centered coordinate system. In
Fig. 5, we can see points $(s,\varphi)$ and the mean before (Fig. 5(a)) and
after translation (Fig. 5(b)).

Consider $N$ LoR-s $(s_{c_{i}}^{(k)},\varphi_{i})$. Samples $\varphi_{i}$
follow a uniform distribution on the interval
$[-\frac{\pi}{2},\frac{\pi}{2}]$, due to the nature of electron-positron
annihilation. Interestingly, $s_{c_{i}}^{(k)}$ follow the distribution whose
probability density function is:

$f_{s_{c}^{(k)}}(s_{c};\sigma_{1}^{(k)},\sigma_{2}^{(k)},\varphi_{0}^{(k)})=\displaystyle\frac{1}{\pi}\int\displaylimits_{-\frac{\pi}{2}}^{\frac{\pi}{2}}\frac{1}{\sqrt{2\pi}\sigma_{p}(\varphi)}e^{-\frac{s_{c}^{2}}{2\sigma^{2}_{p}(\varphi)}}d\varphi.$
(15)

This expression comes from the fact that for fixed angle $\varphi$ we have a
probability density function
$f_{p_{\varphi}}(s_{c};\sigma_{1}^{(k)},\sigma_{2}^{(k)},\varphi_{0}^{(k)})$.
By averaging all contributions for all angles $\varphi$, we get the expression
(15), which is the marginal distribution of $s_{c}^{(k)}$.

(a) Estimation of mean value in the projection domain

(b) Pairs $(s_{c},\varphi)$ in a new coordinate system

Figure 5: Estimation of mean value and LoR-s in the original (left) and
centered (right) coordinate systems

Similarly to the problem from Eq. 4, let $(x_{i})_{i=1}^{N}$ be the samples
that follow distribution $f_{p_{\varphi}}$ from Eq. 13 for some fixed angle
$\varphi$. The least square problem:

$L_{sq}(d;x_{i})=\displaystyle\sum_{n=1}^{N}(d-x_{i}^{2})^{2}$ (16)

has the solution

$d_{LS}=\frac{1}{N}\displaystyle\sum_{n=1}^{N}x_{i}^{2}\approx
E[X^{2}]=\sigma^{2}_{p}(\varphi),$ (17)

where $\sigma^{2}_{p}(\varphi)$ is the variance of the distribution with the
probability density function $f_{p_{\varphi}}$. Therefore, we can state the
weighted least square problem:

$\begin{gathered}L_{sq}(\varphi_{0}^{(k)},\sigma_{1}^{(k)},\sigma_{2}^{(k)})=\\\
\displaystyle\sum\displaylimits_{i=1}^{N}p_{ik}\cdot[(\sigma_{1}^{(k)})^{2}\sin^{2}(\varphi_{i}-\varphi_{0}^{(k)})+(\sigma_{2}^{(k)})^{2}\cos^{2}(\varphi_{i}-\varphi_{0}^{(k)})-(s_{c_{i}}^{(k)})^{2}]^{2}.\end{gathered}$
(18)

Calculating $\nabla
L_{sq}(\varphi_{0}^{(k)},\sigma_{1}^{(k)},\sigma_{2}^{(k)})=\bm{0}$ leads to a
nonlinear system of equations with no closed-form solution. To solve this
system, we need to find the roots of a polynomial with an order greater than
four, which can only be solved numerically.

Therefore, in this paper, we propose a different approach. If we assume that
$\sigma_{1}^{(k)}$ and $\sigma_{2}^{(k)}$ are known, and focus only on
minimizing (18) with respect to rotation $\varphi_{0}^{(k)}$, we get a simple
polynomial system with a closed-form solution. Similarly, once we know
$\varphi_{0}^{(k)}$, and want to find $\sigma_{1}^{(k)}$ and
$\sigma_{2}^{(k)}$ that solve the problem stated by the Eq. 18, we get a
linear $2\times 2$ system of equations.

Let us put the moments of the probability density function $f_{s_{c}}$ in
relation with the variances $\sigma_{1}^{(k)}$ and $\sigma_{2}^{(k)}$. Notice
that odd moments are zero since $f_{s_{c}}$ is an even function. We calculate
the second and the fourth moment. The second moment is

$\begin{gathered}E[S_{c}^{2}]=\frac{1}{\pi}\displaystyle\int\displaylimits_{-\infty}^{+\infty}\displaystyle\int\displaylimits_{-\frac{\pi}{2}}^{\frac{\pi}{2}}s_{c}^{2}f_{s_{c}}(s_{c};\sigma_{1}^{(k)},\sigma_{2}^{(k)},\varphi_{0}^{(k)})d\varphi
ds_{c}=\frac{1}{\pi}\displaystyle\int\displaylimits_{-\frac{\pi}{2}}^{\frac{\pi}{2}}\sigma_{p}^{2}(\varphi)=\\\
\frac{1}{\pi}\displaystyle\int\displaylimits_{-\frac{\pi}{2}}^{\frac{\pi}{2}}[(\sigma_{1}^{(k)})^{2}\cos^{2}(\varphi-\varphi_{0}^{(k)})+(\sigma_{2}^{(k)})^{2}\sin^{2}(\varphi-\varphi_{0}^{(k)})]d\varphi=\frac{(\sigma_{1}^{(k)})^{2}}{2}+\frac{(\sigma_{2}^{(k)})^{2}}{2},\end{gathered}$
(19)

while the fourth moment is

$\begin{gathered}E[S_{c}^{4}]=\frac{1}{\pi}\displaystyle\int\displaylimits_{-\infty}^{+\infty}\displaystyle\int\displaylimits_{-\frac{\pi}{2}}^{\frac{\pi}{2}}s_{c}^{4}f_{s_{c}}(s_{c};\sigma_{1}^{(k)},\sigma_{2}^{(k)},\varphi_{0}^{(k)})d\varphi
ds_{c}=\frac{1}{\pi}\displaystyle\int\displaylimits_{-\frac{\pi}{2}}^{\frac{\pi}{2}}3\sigma_{p}^{4}(\varphi)=\\\
\frac{3}{\pi}\displaystyle\int\displaylimits_{-\frac{\pi}{2}}^{\frac{\pi}{2}}[(\sigma_{1}^{(k)})^{2}\cos^{2}(\varphi-\varphi_{0}^{(k)})+(\sigma_{2}^{(k)})^{2}\sin^{2}(\varphi-\varphi_{0}^{(k)})]^{2}d\varphi=\\\
\frac{9(\sigma_{1}^{(k)})^{4}}{8}+\frac{3(\sigma_{1}^{(k)})^{2}(\sigma_{2}^{(k)})^{2}}{4}+\frac{9(\sigma_{2}^{(k)})^{4}}{8}.\end{gathered}$
(20)

But, second and fourth moments can be estimated from:

$\begin{gathered}M_{2w}^{(k)}=\displaystyle\frac{\sum_{i=1}^{N}p_{ik}(s_{c_{i}}^{(k)})^{2}}{\sum_{n=1}^{N}p_{ik}}\\\
M_{4w}^{(k)}=\displaystyle\frac{\sum_{i=1}^{N}p_{ik}(s_{c_{i}}^{(k)})^{4}}{\sum_{n=1}^{N}p_{ik}}\end{gathered},$
(21)

where $s_{c_{i}}^{(k)}$ are samples as described in Eq. 14 and $p_{ik}$ is the
probability that $i$-th LoR belongs to the $k$-th component of the mixture.

From Eq. 19, 20, and 21 we can estimate $(\sigma_{1}^{(k)})^{2}$ and
$(\sigma_{2}^{(k)})^{2}$. By taking the square root, we get $\sigma_{1}^{(k)}$
and $\sigma_{2}^{(k)}$. The solution of system:

$\begin{gathered}\frac{(\hat{\sigma}_{1}^{(k)})^{2}}{2}+\frac{(\hat{\sigma}_{2}^{(k)})^{2}}{2}=M_{2w}^{(k)}\\\
\frac{9(\hat{\sigma}_{1}^{(k)})^{4}}{8}+\frac{3(\hat{\sigma}_{1}^{(k)})^{2}(\hat{\sigma}_{2}^{(k)})^{2}}{4}+\frac{9(\hat{\sigma}_{2}^{(k)})^{4}}{8}=M_{4w}^{(k)}\end{gathered}$
(22)

gives estimates $(\hat{\sigma}_{1}^{(k)})^{2}$ and
$(\hat{\sigma}_{2}^{(k)})^{2}$. By expressing $(\hat{\sigma}_{2}^{(k)})^{2}$
from the first equation and inserting it into the second one, we get a
quadratic equation in variable $(\hat{\sigma}_{1}^{(k)})^{2}$:

$\begin{gathered}3(\hat{\sigma}_{1}^{(k)})^{4}-6M_{2w}(\hat{\sigma}_{1}^{(k)})^{2}+9M_{2w}^{(k)}-2M_{4w}^{(k)}=0\end{gathered}.$
(23)

Notice that the system in Eq. 22 is symmetric, so we can choose
$(\hat{\sigma}_{1}^{(k)})^{2}\geq(\hat{\sigma}_{2}^{(k)})^{2}$. The solution
of the system (22) is given by Eq. 23, and the final solution is:

$\begin{gathered}(\hat{\sigma}_{1}^{(k)})^{2}=M_{2w}^{(k)}+\sqrt{2}\cdot\sqrt{\frac{M_{4w}^{(k)}}{3}-(M_{2w}^{(k)})^{2}},\\\
(\hat{\sigma}_{2}^{(k)})^{2}=M_{2w}^{(k)}-\sqrt{2}\cdot\sqrt{\frac{M_{4w}^{(k)}}{3}-(M_{2w}^{(k)})^{2}}.\end{gathered}$
(24)

Now, we can calculate $\varphi_{0}^{(k)}$ by solving the weighed least square
problem in (18). By using trigonometric power-reduction formulae we can
rewrite Eq. 18 as:

$\begin{gathered}L_{sq}(\varphi_{0}^{(k)},\sigma_{1}^{(k)},\sigma_{2}^{(k)}))=\\\
\displaystyle\sum\displaylimits_{i=1}^{N}p_{ik}\cdot\bigg{(}\frac{(\sigma_{2}^{(k)})^{2}-(\sigma_{2}^{(k)})^{2}}{2}\cos(2\varphi_{0}^{(k)}-2\varphi_{i})+\frac{(\sigma_{1}^{(k)})^{2}+(\sigma_{2}^{(k)})^{2}}{2}-(s_{c_{i}}^{(k)})^{2}\bigg{)}^{2}.\end{gathered}$
(25)

We set the partial derivative with respect to $\varphi_{0}^{(k)}$ to zero,
i.e. $\frac{\partial L_{sq}}{\partial\varphi_{0}^{(k)}}(\varphi_{0}^{(k)})=0$:

$\begin{gathered}\displaystyle\sum\displaylimits_{i=1}^{N}(M_{ik}\cos(\alpha_{0}^{(k)}-\alpha_{i})+N_{ik})\sin(\alpha_{0}^{(k)}-\alpha_{i})=0,\end{gathered}$
(26)

where

$\begin{gathered}M_{ik}=p_{ik}((\sigma_{2}^{(k)})^{2}-(\sigma_{2}^{(k)})^{2}),\\\
N_{ik}=p_{ik}((\sigma_{1}^{(k)})^{2}+(\sigma_{2}^{(k)})^{2}-2(s_{c_{i}}^{(k)})^{2}),\\\
\alpha_{0}^{(k)}=2\varphi_{0}^{(k)},\text{and}\\\
\alpha_{i}=2\varphi_{i}.\end{gathered}$ (27)

The equation can be written in a more useful form by using angle difference
identities:

$\begin{gathered}A_{s2}\sin^{2}(\alpha_{0}^{(k)})+A_{c2}\cos^{2}(\alpha_{0}^{(k)})+A_{sc}\sin(\alpha_{0}^{(k)})\cos(\alpha_{0}^{(k)})+\\\
A_{s}\sin(\alpha_{0}^{(k)})+A_{c}\cos(\alpha_{0}^{(k)})=0,\end{gathered}$ (28)

where

$\begin{gathered}A_{s2}=\displaystyle\sum\displaylimits_{i=1}^{N}M_{ik}\sin(\alpha_{i})\cos(\alpha_{i}),\\\
A_{c2}=-\displaystyle\sum\displaylimits_{i=1}^{N}M_{ik}\sin(\alpha_{i})\cos(\alpha_{i})=-A_{s2},\\\
A_{sc}=\displaystyle\sum\displaylimits_{i=1}^{N}M_{ik}\cos^{2}(\alpha_{i})-\displaystyle\sum\displaylimits_{i=1}^{N}M_{ik}\sin^{2}(\alpha_{i}),\\\
A_{s}=\displaystyle\sum\displaylimits_{i=1}^{N}N_{ik}\cos(\alpha_{i}),\\\
A_{c}=-\displaystyle\sum\displaylimits_{i=1}^{N}N_{ik}\sin(\alpha_{i}).\\\
\end{gathered}$ (29)

Introducing substitution $x=\cos(\alpha_{0}^{(k)})$ and
$y=\sin(\alpha_{0}^{(k)})$, we get a system of two polynomials equations:

$\begin{cases}\begin{gathered}x^{2}+y^{2}=1,\\\
A_{s2}y^{2}-A_{s2}x^{2}+A_{sc}xy+A_{s}y+A_{c}x=0.\end{gathered}\end{cases}$
(30)

Now, we substitute $y=\pm\sqrt{1-x^{2}}$ in the second equation and calculate:

$\begin{gathered}\pm(A_{sc}x+A_{s})\sqrt{1-x^{2}}=2A_{s2}x^{2}-A_{c}x+A_{s2}\Big{/}^{2}\Rightarrow\\\
(4A_{s2}^{2}+A_{sc}^{2})x^{4}+(2A_{sc}A_{s}-4A_{s2}A_{c})x^{3}+(A_{c}^{2}+A_{s}^{2}-A_{sc}^{2}-4A_{s2}^{2})x^{2}+\\\
(2A_{s2}A_{c}-2A_{sc}A_{s})x+A_{s2}^{2}-A_{s}^{2}=0\end{gathered}$ (31)

Therefore, we get a quartic equation that has an algebraic (close-form)
solution as stated by Abel-Ruffini theorem. More details on how to find roots
of a polynomial of the forth order analytically can be found in [5]. There are
eight pairs of solutions. Four of them do not satisfy the second equation in
(30), leaving us with only four solutions. Two of the remaining candidates are
complex conjugates. Thus, only two real solutions remain - global minimum and
global maximum. We pick the solution
$\hat{\varphi}_{0}^{(k)}=\frac{\hat{\alpha}_{0}^{(k)}}{2}=\frac{1}{2}\arctan(\frac{y}{x})$
with the minimum value of
$L_{sq}(\hat{\varphi}_{0}^{(k)},\hat{\sigma}_{1}^{(k)},\hat{\sigma}_{2}^{(k)})$,
which is the solution of the least square problem stated in Eq. 25.

Pay attention that now we have all data needed for construction of the
covariance matrix, as shown in Eq. 33.

The last step is to reevaluate $\sigma_{1}^{(k)}$ and $\sigma_{2}^{(k)}$ from
the Eq. 18 when angle $\varphi_{0}^{(k)}$ is known. This step is optional, but
useful, especially when the components significantly overlap. To solve the
least square problem from Eq. 18, we set partial derivatives with respect to
$(\sigma_{1}^{(k)})^{2}$ and $(\sigma_{2}^{(k)})^{2}$ to zero (
$\frac{\partial L_{sq}}{\partial(\sigma_{1}^{(k)})^{2}}=\frac{\partial
L_{sq}}{\partial(\sigma_{2}^{(k)})^{2}}=0$ ). We get a linear $2\times 2$
system: $\underbrace{\begin{pmatrix}M_{11}&M_{12}\\\
M_{21}&M_{22}\end{pmatrix}}_{\bm{M}_{k}}\bm{\sigma}^{2}_{k}=\bm{b}_{k}$, where

$\begin{gathered}M_{11}=\displaystyle\sum_{i=1}^{N}p_{ik}\sin^{4}(\varphi_{0}^{(k)}-\varphi_{i}),\\\
M_{12}=M_{21}=\displaystyle\sum_{i=1}^{N}p_{ik}\sin^{2}(\varphi_{0}^{(k)}-\varphi_{i})\cos^{2}(\varphi_{0}^{(k)}-\varphi_{i}),\\\
M_{22}=\displaystyle\sum_{i=1}^{N}p_{ik}\cos^{4}(\varphi_{0}^{(k)}-\varphi_{i}),\\\
\bm{b}_{k}=\begin{pmatrix}\displaystyle\sum_{i=1}^{N}p_{ik}(s_{c_{i}}^{(k)})^{2}\sin^{2}(\varphi_{0}^{(k)}-\varphi_{i})\\\
\displaystyle\sum_{i=1}^{N}p_{ik}(s_{c_{i}}^{(k)})^{2}\cos^{2}(\varphi_{0}^{(k)}-\varphi_{i})^{2}\end{pmatrix},\\\
\bm{\sigma}^{2}_{k}=\begin{pmatrix}(\sigma_{1}^{(k)})^{2}\\\
(\sigma_{2}^{(k)})^{2}\end{pmatrix}.\\\ \end{gathered}$ (32)

The solution of the above system is
$\hat{\bm{\sigma}}^{2}_{k}=\bm{M}_{k}^{-1}\bm{b}_{k}$.

Finally, from known estimates $\hat{\sigma_{1}}^{(k)}$,
$\hat{\sigma}_{2}^{(k)}$ and $\hat{\varphi}_{0}^{(k)}$ we can construct
covariance matrix $\bm{\Sigma}_{k}$. The covariance matrix can be decomposed
as:

$\bm{\Sigma}_{k}=\bm{U}_{k}\bm{D}_{k}\bm{U}_{k}^{T}$ (33)

where
$\bm{D}_{k}=diag((\hat{\sigma}_{1}^{(k)})^{2},(\hat{\sigma}_{2}^{(k)})^{2})$
is a diagonal matrix and $\bm{U}_{k}=[\bm{u}_{1}\;\;\bm{u}_{2}]$ is a unitary
(rotation) matrix. Corresponding eigenvectors $\bm{u}_{1}$ and $\bm{u}_{2}$
are unit vectors in direction of the principal axes. Eigenvectors can be
expressed in terms of angle $\varphi_{0}^{(k)}$ as

$\begin{gathered}\bm{u}_{1}=\begin{pmatrix}\cos(\hat{\varphi}_{0}^{(k)})\\\
\sin(\hat{\varphi}_{0}^{(k)})\end{pmatrix},\\\
\bm{u}_{2}=\begin{pmatrix}\cos(\hat{\varphi}_{0}^{(k)}+\frac{\pi}{2})\\\
\sin(\hat{\varphi}_{0}^{(k)}+\frac{\pi}{2})\end{pmatrix}.\end{gathered}$ (34)

In Fig. 6, we can see eigenvectors and the corresponding square root of
eigenvalues $\hat{\sigma_{1}}^{(k)}$, $\hat{\sigma}_{2}^{(k)}$.

Figure 6: Interpretation of eigenvectors and eigenvalues of covariance matrix

##  4 Proposed Algorithm

At the beginning, we arbitrarily initialize mean vectors, and randomly assign
LoR-s to a subset of Gaussian components. We remind the reader that the number
of components $K$ is known in advance. The proposed algorithm consists of two
steps. Both steps are iterative, and the first step provides for a rough
estimate of the unknown parameters. In each iteration of the second step we
recalculate the mean vectors, the covariance matrices, as well as the
membership probabilities.

###  4.1 First step - initialization

Initially, we assign lines of response randomly to Gaussian components. Each
component should have a roughly equal number of associated LoR-s
$\approx\frac{N}{K}$. In this step, we do the hard (or modal) classification
of LoR-s. That means that each LoR belongs to only one of the components, i.e.
the membership probabilities are binary, and the number of LoR-s of $k$-th
component is

$\hat{L}_{k}=\sum_{i=1}^{N}p_{ik}.$ (35)

We calculate the mean vectors for each Gaussian component as described by Eq.
8 and 10. Then, we reassign the LoR-s by pairing them to the least distanced
mean vectors. Notice that the membership probabilities are still binary. We
repeat the iterative procedure (mean, reassignment), until the mean vectors
changes are sufficiently small.

From known probabilities $p_{ik}$ (zeros and ones) and mean vectors for each
component in the mixture, we calculate the initial estimate of the covariance
matrices. Thanks to the analytical results presented in the previous Section,
the calculation is straightforward. We estimate $\sigma_{1}^{(k)}$ and
$\sigma_{2}^{(k)}$ from moments for each component separately as described in
Eq. 21, 22 and 24. Afterward, we calculate $\hat{\varphi}_{0}^{(k)}$ as
discussed in the previous Section. Then we reevaluate $\hat{\sigma}_{1}^{(k)}$
and $\hat{\sigma}_{2}^{(k)}$ as stated in Eq. 32 and update
$\hat{\varphi}_{0}^{(k)}$ with new $\hat{\sigma}_{1}^{(k)}$ and
$\hat{\sigma}_{2}^{(k)}$. Finally, we calculate each covariance matrix via
eigenvalue decomposition (Eq. 33).

###  4.2 Second step

The second step is similar to the first step. The main difference is the way
we calculate probabilities $p_{ik}$. Instead of hard classification, in this
step we utilize soft (proportional) classification. Let us focus on only one
LoR with a total of $K$ components in the mixture. We want to determine
$p_{ik}$ \- the probability that $i$-th line belongs to $k$-th component. For
each component, we calculate the line integral:

$\begin{gathered}\tilde{p}_{ik}=\int_{\gamma_{i}}\tau_{k}f_{G}(\bm{x};\bm{\mu}_{k},\bm{\Sigma}_{k})d\gamma_{i},\end{gathered}$
(36)

where $f_{G}$ denoted a bivariate normal distribution and
$\tau_{k}=\frac{L_{k}}{N}$ is an estimated weight of the $k$-th component. The
mean vector $\bm{\mu}_{k}$ and covariance matrix $\bm{\Sigma}_{k}$ are known
from the previous step. Mapping $\gamma_{i}$ is a parametrization of the line
of response we are currently observing. The solution to this integral is a
point from the univariate normal distribution (see A and Eq. 13). Keep in mind
that, since $i$-th response have occurred, the overall probability must be 1.
Therefore, we define

$\begin{gathered}p_{ik}=\frac{\tilde{p}_{ik}}{\sum_{k=1}^{K}\tilde{p}_{ik}}.\end{gathered}$
(37)

We keep the same ratio as in $\tilde{p}_{ik}$, but force the probabilities to
add up to one for each LoR. Additionally, the estimation of the number of
LoR-s from Eq. 35 still holds.

After we updated all membership probabilities $p_{ik}$, we calculate a new
mean vector $\bm{\mu}_{k}$ as already described (Eq. 10). With newly estimated
mean vectors, we determine new covariance matrices $\bm{\Sigma}_{k}$ in the
same way as presented in the last subsection. Then we repeat the procedure,
i.e. we calculate $p_{ik}$ from $\bm{\mu}_{k}$ and $\bm{\Sigma}_{k}$. This can
be described in 8 steps:

  1. 1.

Calculate $p_{ik}$ (Eq. 36 and Eq. 37) for all lines $i$ and all components
$k$ that best fit the current model.

  2. 2.

Estimate the mean vector for each component according to newly obtained
$p_{ik}$.

  3. 3.

Calculate $\hat{\sigma}_{1}^{(k)}$ and $\hat{\sigma}_{2}^{(k)}$ from Eq. 24.

  4. 4.

Determine $\hat{\varphi}_{0}^{(k)}$ from Eq. 30 and 31.

  5. 5.

Reevaluate $\hat{\sigma}_{1}^{(k)}$ and $\hat{\sigma}_{2}^{(k)}$ from Eq. 32.

  6. 6.

Update $\hat{\varphi}_{0}^{(k)}$ from Eq. 30 and 31 with new
$\hat{\sigma}_{1}^{(k)}$ and $\hat{\sigma}_{2}^{(k)}$.

  7. 7.

Use eigenvalue decomposition to calculate the covariance matrix (Eq. 33).

  8. 8.

Do 3)-7) for every component in the mixture ($k=1,2,3,..,K$).

  9. 9.

Repeat 1)-8) until the change in estimated weight of all components is
sufficiently small.

We obtained all mixture parameters: $\bm{\mu}_{k}$ and $\bm{\Sigma}_{k}$ and
$\tau_{k}$. Such a set of parameters represents continuous model of the
reconstructed image. Pay attention that the proposed algorithm is numerically
efficient, in spite of the number of equations in the previous
section.111Implementation of described algorithm is available on
https://github.com/tm2005/Accurate-PET-Reconstruction-from-Reduced-Set-of-
Measurements-based-on-GMM.

##  5 Results

Our test mixture consists of three Gaussian components, in which two
components highly overlap.

Figure 7: The probability density function of the GMM (Ground truth)

The first component has mean vector $\bm{\mu}_{1}=\begin{pmatrix}0\\\
0\end{pmatrix}$ and covariance matrix
$\bm{\Sigma}_{1}=\begin{pmatrix}0.0625&0\\\ 0&0.0625\end{pmatrix}$. The second
component highly overlaps with the first component and its parameters are
$\bm{\mu}_{2}=\begin{pmatrix}-0.4\\\ -0.4\end{pmatrix}$ and
$\bm{\Sigma}_{2}=\begin{pmatrix}0.04&0.03\\\ 0.03&0.09\end{pmatrix}$. The
third component has the mean $\bm{\mu}_{3}=\begin{pmatrix}1.25\\\
-1\end{pmatrix}$ and covariance $\bm{\Sigma}_{3}=\begin{pmatrix}0.04&0.006\\\
0.006&0.01\end{pmatrix}$. During the simulation, we created $N_{1}=3500$ lines
that originate from the first component, $N_{2}=2500$ lines that originate
from the second component, and $N_{3}=1000$ from the last component. Hence,
the corresponding weights in the Gaussian mixture model are $\tau_{1}=0.5$,
$\tau_{2}=\frac{5}{14}\approx 0.36$, and $\tau_{3}=\frac{1}{7}\approx 0.14$.
The PDF of our GMM is (as depicted in Fig 7):

$\begin{gathered}g(x,y)=\tau_{1}\cdot
f_{G}(x,y;\bm{\mu}_{1},\bm{\Sigma}_{1})+\tau_{2}\cdot
f_{G}(x,y;\bm{\mu}_{2},\bm{\Sigma}_{2})\\\ +\tau_{3}\cdot
f_{G}(x,y;\bm{\mu}_{3},\bm{\Sigma}_{3}),\end{gathered}$ (38)

where $f_{G}$ is bivariate normal distribution as in Eq 3.

We have repeated the simulation $n=100$ times with the same parameters and
applied the proposed algorithm to determine the unknown parameters of the
Gaussian mixture model. To check our algorithm, we measure four errors. The
first one measures an error between the ideal and estimated mean vector

$\bm{\mu}_{er_{i}}=||\bm{\hat{}}{\bm{\mu}}_{i}-\bm{\mu}_{i}||_{2}.$ (39)

The second one evaluates the error between the ideal and estimated covariance
matrix

$\bm{\Sigma}_{er_{i}}=||\bm{\hat{}}{\bm{\Sigma}}_{i}-\bm{\Sigma}_{i}||_{F}.$
(40)

Also, we measure the error between the ideal and estimated weights in the GMM

$\tau_{er_{i}}=|\hat{\tau}_{i}-\tau_{i}|.$ (41)

(a) Mean vector error

(b) Covariance matrix error

(c) Mixture weight error

Figure 8: Errors of individual parameters of GMM in each simulation

Unlike the mentioned measures that deal with each parameter separately, the
last one is the Kullback-Leibler (KL) divergence. It measures an overall
error, since it compares the two PDF-s:

$D_{KL}(\hat{g}|g)=\int\displaylimits_{-\infty}^{+\infty}\int\displaylimits_{-\infty}^{+\infty}\hat{g}(x,y)\log\left(\frac{\hat{g}(x,y)}{g(x,y)}\right)dxdy,$
(42)

where $g(x,y)$ is the ideal PDF (Eq. 38) and $\hat{g}(x,y)$ is the estimate of
the PDF obtained via the proposed algorithm.

In Figure 8, we can see mean, covariance, and components weight error in each
simulation. Pay attention that the non-overlapping component has smaller
errors. Moreover, all errors are very small, and the averages are presented in
Table 1.

Component | Avg. Mean | Avg. Cov. | Avg. Weight
---|---|---|---
First | 0.035 | 0.014 | 0.019
Second | 0.029 | 0.021 | 0.018
Third | 0.011 | 0.004 | 0.002
Table 1: Average error of parameters of GMM in 100 simulations

The KL divergence of 100 simulations is depicted in Fig. 9. The calculation of
the KL divergence is done numerically. The KL divergence is always smaller
than $0.023$, while the mean value of all simulations is $0.013$.

Figure 9: KL divergence in each simulation

From all that has been said, we conclude that the estimated PDF of GMM highly
resembles the true PDF. A typical reconstruction of the observed setup can be
seen in Fig. 10. As expected, the difference between the true PDF (Fig. 7) and
the estimated one is not noticeable.

Figure 10: Estimated PDF of the GMM

##  6 Conclusion

We presented a novel method for estimation of unknown parameters of the
Gaussian mixture model in Positron Emission Tomography. In contrast to
competitive reconstruction methods based on pixel or voxel grid, the obtained
GMM is continuous, and virtually of infinite resolution. It can be directly
used for further analytical processing. There are several well-known GMM
estimation methods based on samples. In PET imaging, the samples are not
known. The challenge was to estimate the GMM parameters from projections,
namely from the lines of response that fires from unknown radioactive decay
points under some random angles. We solved the problem analytically, and
proposed a systematic approach for obtaining the unknown parameters, under the
assumption that the number of Gaussian components is known.

First, we estimate the mean vectors for each component in the mixture. For
each component, we get a $2\times 2$ linear system of equations. Then, the
covariance matrices of each component are estimated in four steps. In the
first step, we propose to use higher moments of each Gaussian component to
estimate variances in the principal axes. Then, we obtain the direction of the
larger principal axis. In the third step, we additionally tune the variance
estimates. Finally, we use the eigenvalue decomposition to get the desired
covariance matrix of each component. Weights of Gaussian components are given
by the membership probabilities. Described steps are integrated into an
iterative ML-like algorithm.

The results presented in this paper show that recovery of the unknown
parameters is possible even when two components significantly overlap. A
relatively small number of LoR-s is needed for an accurate reconstruction,
thus eventually leading to lower radiation doses in the PET imaging.

## Funding

This research was supported by the Croatian Science Foundation
[IP-2019-04-6703].

##  Appendix A Line integration of centered and not rotated bivariate normal
distribution

In this appendix we present a detailed procedure on how to calculate an
integral stated in Eq. 11:

$f_{p_{0_{\varphi}}}(s)=\int_{\gamma_{\varphi}}f_{Gc}(x,y)d\gamma_{\varphi}.$

We remind the reader that $|\varphi|\leq\frac{\pi}{2}$. Since we are dealing
with a line integral, we calculate $||\gamma_{\varphi}^{\prime}(t)||_{2}$:

$||\gamma_{\varphi}^{\prime}(t)||_{2}=||(1,k)||_{2}=\sqrt{1+k^{2}}=\frac{1}{\cos(\varphi)}.$

We calculate
$f_{p_{0_{\varphi}}}(s)=\int_{\gamma_{\varphi}}f_{Gc}(x,y)d\gamma_{\varphi}$:

$\begin{gathered}\int\displaylimits_{\gamma_{\varphi}}f_{Gc}(x,y)d\gamma=\int\displaylimits_{-\infty}^{+\infty}\frac{1}{2\pi\sigma_{1}\sigma_{2}}\exp(-(\frac{t^{2}}{\sigma_{1}^{2}}+\frac{(kt+l)^{2}}{\sigma_{2}^{2}}))\cdot||\gamma^{\prime}(t)||dt=\\\
\frac{\exp(-\frac{t^{2}}{2\sigma_{2}^{2}})}{2\pi\sigma_{1}\sigma_{2}\cos(\varphi)}\int\displaylimits_{-\infty}^{+\infty}\exp(-(t^{2}(\frac{1}{2\sigma_{1}^{2}}+\frac{k^{2}}{2\sigma_{2}^{2}})+\frac{kl}{\sigma_{2}^{2}}t))dt=\\\
\frac{\exp(-\frac{t^{2}}{2\sigma_{2}^{2}})\exp(\frac{k^{2}t^{2}}{4\sigma_{2}^{4}}\frac{1}{\frac{1}{2\sigma_{1}^{2}}+\frac{k^{2}}{2\sigma_{2}^{2}}})}{2\pi\sigma_{1}\sigma_{2}\cos(\varphi)}\int\displaylimits_{-\infty}^{+\infty}\exp(-(t\sqrt{\frac{1}{2\sigma_{1}^{2}}+\frac{k^{2}}{2\sigma_{2}^{2}})}+\frac{1}{\sqrt{\frac{1}{2\sigma_{1}^{2}}+\frac{k^{2}}{2\sigma_{2}^{2}})}}\frac{kl}{2\sigma_{2}^{2}})^{2})dt=\\\
\frac{1}{2\pi\sigma_{1}\sigma_{2}\cos(\varphi)}\frac{\sqrt{\pi}}{\sqrt{\frac{1}{2\sigma_{1}^{2}}+\frac{k^{2}}{2\sigma_{2}^{2}}}}\exp(\frac{-l^{2}}{2\sigma_{2}^{2}}(1-\frac{2\sigma_{1}^{2}\sigma_{2}^{2}}{\sigma_{2}^{2}+k^{2}\sigma_{1}^{2}}\frac{1}{2\sigma
2^{2}}))=\\\
\frac{1}{\sqrt{2\pi}\sigma_{1}\sigma_{2}\cos(\varphi)}\frac{\sigma_{1}\sigma_{2}}{\sqrt{\sigma_{2}^{2}+k^{2}\sigma_{1}^{2}}}e^{\frac{-l^{2}}{2(\sigma_{2}^{2}+k^{2}\sigma_{1}^{2})}}=\\\
\frac{1}{\sqrt{2\pi}}\frac{1}{\cos(\varphi)\sqrt{\sigma_{2}^{2}+\tan^{2}(\varphi)\sigma_{1}^{2}}}\exp(\frac{-s^{2}}{2\cos^{2}(\varphi)(\sigma_{2}^{2}+\tan^{2}(\varphi)\sigma_{1}^{2})})=\\\
\frac{1}{\sqrt{2\pi}\sqrt{\sigma_{2}^{2}\cos^{2}(\varphi)+\sigma_{1}^{2}\sin^{2}(\varphi)}}\exp(\frac{-s^{2}}{2(\sigma_{2}^{2}\cos^{2}(\varphi)+\sigma_{1}^{2}\sin^{2}(\varphi))})=\\\
\frac{1}{\sqrt{2\pi}\sigma_{p_{0}}(\varphi)}\exp(-\frac{s^{2}}{2\sigma_{p_{0}}(\varphi)^{2}}).\end{gathered}$

Here, we used known result:

$\int\displaylimits_{-\infty}^{+\infty}\exp(-(at+b)^{2})dt=\frac{\sqrt{\pi}}{|a|},\,\,\,a\neq
0.$

## References

  * [1] Animashree Anandkumar, Daniel Hsu, and Sham M. Kakade.  A method of moments for mixture models and hidden markov models.  In Proceedings of the 25th Annual Conference on Learning Theory, volume 23 of Proceedings of Machine Learning Research, pages 33.1–33.34. PMLR, 25–27 Jun 2012.
  * [2] Yannick Berker, Joscha Maier, and Marc Kachelries.  Deep scatter estimation in PET: Fast scatter correction using a convolutional neural network.  In 2018 IEEE Nuclear Science Symposium and Medical Imaging Conference Proceedings (NSS/MIC). IEEE, November 2018.
  * [3] Simon R. Cherry and Magnus Dahlbom.  PET: Physics, instrumentation, and scanners.  In PET, pages 1–117. Springer New York.
  * [4] Michel Defrise and Paul Kinahan.  Data acquisition and image reconstruction for 3d PET.  In The Theory and Practice of 3D PET, pages 11–53. Springer Netherlands, 1998.
  * [5] Amir Fathi and Nastaran Sharifan.  A classic new method to solve quartic equations.  Applied and Computational Mathematics, 2(2):24, 2013.
  * [6] Nir Friedman and Stuart Russell.  Image segmentation in video sequences: A probabilistic approach.  In UAI’97: Proceedings of the Thirteenth conference on Uncertainty in artificial intelligence, pages 175–181. Morgan Kaufmann Publishers Inc., August 1997.
  * [7] Kuang Gong, Jiahui Guan, Kyungsang Kim, Xuezhu Zhang, Jaewon Yang, Youngho Seo, Georges El Fakhri, Jinyi Qi, and Quanzheng Li.  Iterative PET image reconstruction using convolutional neural network representation.  IEEE Transactions on Medical Imaging, 38(3):675–685, March 2019\.
  * [8] Ida Häggström, C. Ross Schmidtlein, Gabriele Campanella, and Thomas J. Fuchs.  DeepPET: A deep encoder–decoder network for directly solving the PET image reconstruction inverse problem.  Medical Image Analysis, 54:253–262, May 2019.
  * [9] Xiao Han.  MR-based synthetic CT generation using a deep convolutional neural network method.  Medical Physics, 44(4):1408–1419, March 2017.
  * [10] Trevor Hastie, Robert Tibshirani, and Jerome Friedman.  The Elements of Statistical Learning.  Springer New York, 2009.
  * [11] H.M. Hudson and R.S. Larkin.  Accelerated image reconstruction using ordered subsets of projection data.  IEEE Transactions on Medical Imaging, 13(4):601–609, 1994.
  * [12] Daniel M. Kane.  Robust learning of mixtures of gaussians.  In Proceedings of the Thirty-Second Annual ACM-SIAM Symposium on Discrete Algorithms, pages 1246–1258. Society for Industrial and Applied Mathematics, 2021.
  * [13] Parisa Khateri, Jannis Fischer, Werner Lustermann, Charalampos Tsoumpas, and Günther Dissertori.  Implementation of cylindrical PET scanners with block detector geometry in STIR.  EJNMMI Physics, 6(1), July 2019.
  * [14] Rima Khouja, Pierre-Alexandre Mattei, and Bernard Mourrain.  Tensor decomposition for learning gaussian mixtures from moments.  Journal of Symbolic Computation, 113:193–210, November 2022.
  * [15] P.E. Kinahan and J.G. Rogers.  Analytic 3d image reconstruction using all detected events.  IEEE Transactions on Nuclear Science, 36(1):964–968, February 1989.
  * [16] Anita Gribl Koscevic and Davor Petrinovic.  Extra-low-dose 2d PET imaging.  In 2021 12th International Symposium on Image and Signal Processing and Analysis (ISPA). IEEE, September 2021.
  * [17] Richard M. Leahy and Jinyi Qi.  Statistics and Computing, 10(2):147–165, 2000.
  * [18] Flavien Leger, Guoshen Yu, and Guillermo Sapiro.  Efficient matrix completion with gaussian models.  In 2011 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, May 2011.
  * [19] Jichuan Li and Arye Nehorai.  Gaussian mixture learning via adaptive hierarchical clustering.  Signal Processing, 150:116–121, September 2018.
  * [20] Saikat Majumder.  A gaussian mixture model method for eigenvalue-based spectrum sensing with uncalibrated multiple antennas.  Signal Processing, 192:108404, March 2022.
  * [21] T. Matulic, R. Bagaric, and D. Sersic.  Enhanced reconstruction for PET scanner with a narrow field of view by using backprojection method.  In 2021 44th International Convention on Information, Communication and Electronic Technology (MIPRO). IEEE, September 2021.
  * [22] Kevin P. Murphy.  Machine Learning: A Probabilistic Perspective.  The MIT Press, 2012.
  * [23] D. Panetta and N. Camarlinghi.  3D Image Reconstruction for CT and PET: A Practical Guide with Python.  CRC Press.
  * [24] Donya Rahmani, Mahesan Niranjan, Damien Fay, Akiko Takeda, and Jacek Brodzki.  Estimation of gaussian mixture models via tensor moments with application to online learning.  Pattern Recognition Letters, 131:285–292, March 2020.
  * [25] M. Raitoharju, Á.F. García-Fernández, R. Hostettler, R. Piché, and S. Särkkä.  Gaussian mixture models for signal mapping and positioning.  Signal Processing, 168:107330, March 2020.
  * [26] Ivan Ralasic, Azra Tafro, and Damir Sersic.  Statistical compressive sensing for efficient signal reconstruction and classification.  In 2018 4th International Conference on Frontiers of Signal Processing (ICFSP). IEEE, September 2018.
  * [27] Douglas Reynolds.  Gaussian mixture models.  In Encyclopedia of Biometrics, pages 827–832. Springer US, 2015\.
  * [28] Douglas A. Reynolds, Thomas F. Quatieri, and Robert B. Dunn.  Speaker verification using adapted gaussian mixture models.  Digital Signal Processing, 10(1-3):19–41, January 2000.
  * [29] L. A. Shepp and Y. Vardi.  Maximum likelihood reconstruction for emission tomography.  IEEE Transactions on Medical Imaging, 1(2):113–122, October 1982\.
  * [30] Azra Tafro, Damir Seršić, and Ana Sović Kržić.  2d PET image reconstruction using robust l1 estimation of the gaussian mixture model.  Informatica, pages 653–669, 2022.
  * [31] Juan José Vaquero and Paul Kinahan.  Positron emission tomography: Current challenges and opportunities for technological advances in clinical and preclinical imaging systems.  Annual Review of Biomedical Engineering, 17(1):385–414, December 2015.
  * [32] Yan Wang, Biting Yu, Lei Wang, Chen Zu, David S. Lalush, Weili Lin, Xi Wu, Jiliu Zhou, Dinggang Shen, and Luping Zhou.  3d conditional generative adversarial networks for high-quality PET image estimation at low dose.  NeuroImage, 174:550–562, July 2018.
  * [33] Ying Wang, Xuezhu Zhang, Mengxi Zhang, Dong Liang, Xin Liu, Hairong Zheng, Yongfeng Yang, and Zhanli Hu.  An improved PET image reconstruction method based on super-resolution.  Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, 946:162677, December 2019.
  * [34] Yihong Wu and Pengkun Yang.  Optimal estimation of gaussian mixtures via denoised method of moments.  The Annals of Statistics, 48(4), August 2020.
  * [35] Zhaoheng Xie, Reheman Baikejiang, Tiantian Li, Xuezhu Zhang, Kuang Gong, Mengxi Zhang, Wenyuan Qi, Evren Asma, and Jinyi Qi.  Generative adversarial network based regularized image reconstruction for PET.  Physics in Medicine & Biology, 65(12):125016, June 2020.
  * [36] Guoshen Yu and Guillermo Sapiro.  Statistical compressed sensing of gaussian mixture models.  IEEE Transactions on Signal Processing, 59(12):5842–5858, December 2011.
  * [37] Song Zhang and Yuan Mei Wang.  An approach to positron emission tomography based on penalized cross-entropy minimization.  Signal Processing, 81(5):1069–1074, May 2001.