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

TL;DR

This paper proposes a novel PET image reconstruction method using Gaussian mixture models and robust L1 estimation, improving resolution and reducing noise compared to traditional pixel-based approaches.

Contribution

It introduces a robust segmentation algorithm based on Gaussian mixture models and an iterative EM-like algorithm for PET image reconstruction.

Findings

Enhanced image resolution and noise robustness.

Reduced computational complexity and convergence time.

Demonstrated effectiveness on PET imaging data.

Abstract

An image or volume of interest in positron emission tomography (PET) is reconstructed from pairs of gamma rays emitted from a radioactive substance. Many image reconstruction methods are based on estimation of pixels or voxels on some predefined grid. Such an approach is usually associated with limited resolution of the reconstruction, high computational complexity due to slow convergence and noisy results. This paper explores reconstruction of PET images using the underlying assumption that the originals of interest can be modeled using Gaussian mixture models. A robust segmentation method based on statistical properties of the model is presented, with an iterative algorithm resembling the expectation-maximization algorithm. Use of parametric models for image description instead of pixels circumvent some of the mentioned limitations.

Figures

Tables

Table 1. TABLE I: Mean estimation error, K = 1 𝐾 1 K=1

$N=1000$	${𝒔}_1$	${𝒔}_2$	${𝒔}_3$
$L_1$ method	13.78%	11.88%	9.28%
$L_2$ method	8.27%	7.61%	7.6%

Equations

$$p(\mathbf{x}|{\tau }_k,{\mathbf{\mu }}_k,{\mathbf{\Sigma }}_k)=\sum _{k=1}^K{\tau }_{k}g(\mathbf{x}|{\mathbf{\mu }}_k,{\mathbf{\Sigma }}_k),$$

$$g(\bm{x}|\bm{\mu}_{k},\bm{\Sigma}_{k})=\frac{1}{\sqrt{(2\pi)^{d}|\bm{\Sigma}_{k}|}}\exp\Bigg{(}-\frac{1}{2}(\bm{x}-\bm{\mu}_{k})^{T}\bm{\Sigma}_{k}^{-1}(\bm{x}-\bm{\mu}_{k})\Bigg{)},$$

$${\mathbf{a}}_i^T\mathbf{x}+l_i=0,{\mathbf{a}}_i=[\tan {\psi }_i-1{]}^T,i=1,\dots ,N.$$

$$d^2({\mathbf{v}}_1,{\mathbf{v}}_2)=({\mathbf{v}}_1-{\mathbf{v}}_2{)}^T\mathbf{W}({\mathbf{v}}_1-{\mathbf{v}}_2),$$

$$\min d(\mathbf{x},\mathbf{\mu })\text{ s.t. }{\mathbf{a}}_i\mathbf{x}+l_i=0,$$

$$\left[\begin{array}{cc}2\mathbf{W}& {\mathbf{a}}_i\\ {\mathbf{a}}_i^T& 0\end{array}\right]\cdot \left[\begin{array}{c}{\mathbf{x}}_i^m\\ \lambda \end{array}\right]=\left[\begin{array}{c}2\mathbf{W}\mathbf{\mu }\\ -l_i\end{array}\right].$$

$$\min \sum _{i=1}^N({\mathbf{x}}_i^m-\mathbf{\mu }{)}^T\mathbf{W}({\mathbf{x}}_i^m-\mathbf{\mu }),$$

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

$$\tilde{\mathbf{W}}=\left[\begin{array}{cc}\mathbf{W}& 0\\ 0& 0\end{array}\right]\in M_3,{\mathbf{B}}_i={\left({\left[\begin{array}{cc}2\mathbf{W}& {\mathbf{a}}_i\\ {\mathbf{a}}_i^T& 0\end{array}\right]}^{-1}\right)}^T\cdot \tilde{\mathbf{W}}\in M_3,$$

$$\left[\begin{array}{c}{\mathbf{M}}_i\\ {\mathbf{m}}_i\end{array}\right]={\mathbf{B}}_i\left(\left[\begin{array}{cc}2\mathbf{W}& {\mathbf{a}}_i\\ {\mathbf{a}}_i^T& 0\end{array}\right]\left[\begin{array}{c}2\mathbf{W}\\ 0\end{array}\right]-\left[\begin{array}{cc}1& 0\\ 0& 1\\ 0& 0\end{array}\right]\right),$$

$$\left[\begin{array}{c}{\mathbf{N}}_i\\ 0\end{array}\right]={\mathbf{B}}_i^T\left[\begin{array}{c}2\mathbf{W}\\ 0\end{array}\right]-\tilde{\mathbf{W}}\left[\begin{array}{cc}1& 0\\ 0& 1\\ 0& 0\end{array}\right],{\mathbf{M}}_i,{\mathbf{N}}_i\in M_2,{\mathbf{m}}_i\in M_{12}.$$

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

$$(R\mathbf{\mu },R\mathbf{\Sigma }{R}^T),\text{ where }R=\left[\begin{array}{cc}\cos \phi & -\sin \phi \\ \sin \phi & \cos \phi \end{array}\right].$$

$$(R\mathbf{\mu }{)}_1=\cos \phi {\mu }_x-\sin \phi {\mu }_y,$$

$$(R\mathbf{\Sigma }{R}^T{)}_{11}={\cos }^2\phi {\Sigma }_{11}-2\cos \phi \sin \phi {\Sigma }_{12}+{\sin }^2\phi {\Sigma }_{22},$$

$$\mathbf{A}\mathbf{s}=\mathbf{b},\text{ where }\mathbf{s}=\left[\begin{array}{c}{\Sigma }_{11}\\ {\Sigma }_{12}\\ {\Sigma }_{22}\end{array}\right],$$

$$A=\left[\begin{array}{ccc}{\cos }^2{\phi }_1& -2\sin {\phi }_1\cos {\phi }_1& {\sin }^2{\phi }_1\\ {\cos }^2{\phi }_2& -2\sin {\phi }_2\cos {\phi }_2& {\sin }^2{\phi }_2\\ ⋮& ⋮& ⋮\\ {\cos }^2{\phi }_N& -2\sin {\phi }_1\cos {\phi }_N& {\sin }^2{\phi }_N\end{array}\right],\text{and }\mathbf{b}=\left[\begin{array}{c}(\cos {\phi }_1{\mu }_x-\sin {\phi }_1{\mu }_y+l_1\sin {\phi }_1{)}^2\\ ⋮\\ (\cos {\phi }_N{\mu }_x-\sin {\phi }_N{\mu }_y+l_N\sin {\phi }_N{)}^2\end{array}\right],$$

$$\underset{\mathbf{s}}{\min }\Vert \mathbf{A}\mathbf{s}-\mathbf{b}\Vert .$$

$${\mathbf{A}}^{\mathbf{T}}\mathbf{A}\mathbf{s}={\mathbf{A}}^{\mathbf{T}}\cdot \mathbf{b},\text{ i.e. }\mathbf{s}=({\mathbf{A}}^{\mathbf{T}}\mathbf{A}{)}^{-1}{\mathbf{A}}^{\mathbf{T}}\cdot \mathbf{b}.$$

$$\mathbf{A}\mathbf{s}=k\cdot \mathbf{b},$$

$$k=(\Phi (0.75){)}^2\approx 1.4826^2,$$

$$\mathbf{W}=(1-\alpha )\mathbf{I}+\alpha {\hat{\mathbf{\Sigma }}}^{-1},$$

$${\mathbf{\Sigma }}_1=0.05\mathbf{I},{\mathbf{\Sigma }}_2=\left[\begin{array}{cc}0.02& -0.01\\ -0.01& 0.05\end{array}\right],{\mathbf{\Sigma }}_3=\left[\begin{array}{cc}0.01& 0.02\\ 0.02& 0.05\end{array}\right].$$

$$\mathbf{\Sigma }=\left[\begin{array}{cc}{\sigma }_x^2& \rho {\sigma }_x{\sigma }_y\\ \rho {\sigma }_x{\sigma }_y& {\sigma }_y^2\end{array}\right],$$

$$\mathbf{s}=[{\Sigma }_{11}{\Sigma }_{12}{\Sigma }_{22}{]}^T=[{\sigma }_x^2\rho {\sigma }_x{\sigma }_y{\sigma }_y^2{]}^T$$

$${\mathbf{s}}_1=\left[\begin{array}{c}0.05\\ 0\\ 0.05\end{array}\right],{\mathbf{s}}_2=\left[\begin{array}{c}0.02\\ -0.01\\ 0.05\end{array}\right],{\mathbf{s}}_3=\left[\begin{array}{c}0.01\\ 0.02\\ 0.05\end{array}\right].$$

$$\Vert \mathbf{e}\Vert =\frac{\Vert \hat{\mathbf{s}}-\mathbf{s}\Vert }{\Vert \mathbf{s}\Vert },$$

$$\left[\begin{array}{c}{\mathbf{x}}_i^m\\ \lambda \end{array}\right]={\left[\begin{array}{cc}2\mathbf{W}& {\mathbf{a}}_i\\ {\mathbf{a}}_i^T& 0\end{array}\right]}^{-1}\cdot \left[\begin{array}{c}2\mathbf{W}\mathbf{\mu }\\ -l_i\end{array}\right].$$

$$d^2=\sum _{i=1}^N({\mathbf{x}}_i^m-\mathbf{\mu }{)}^T\mathbf{W}({\mathbf{x}}_i^m-\mathbf{\mu })$$

$$=\sum _{i=1}^N{\left(\left[\begin{array}{c}{\mathbf{x}}_i^m\\ \lambda \end{array}\right]-\left[\begin{array}{c}\mathbf{\mu }\\ \lambda \end{array}\right]\right)}^T\left[\begin{array}{cc}\mathbf{W}& 0\\ 0& 0\end{array}\right]\left(\left[\begin{array}{c}{\mathbf{x}}_i^m\\ \lambda \end{array}\right]-\left[\begin{array}{c}\mathbf{\mu }\\ \lambda \end{array}\right]\right),$$

$${\mathbf{x}}_{i\lambda }=\left[\begin{array}{c}{\mathbf{x}}_i^m\\ \lambda \end{array}\right],{\mathbf{\mu }}_{\lambda }=\left[\begin{array}{c}\mathbf{\mu }\\ \lambda \end{array}\right]\text{ and }\tilde{\mathbf{W}}=\left[\begin{array}{cc}\mathbf{W}& 0\\ 0& 0\end{array}\right].$$

Full text

I Introduction

Positron emission tomography (PET) scanners detect the annihilation photon pairs arising from the positron emissions. Image reconstruction implies generating a two- or three-dimensional image of a radiotracer’s concentration to estimate physiologic parameters for volumes of interest in vivo. In three-dimensional scanners, we can consider the tube (parallelepiped) joining any two detector elements as a volume of response (VOR). In the absence of physical effects and noise, the total number of coincidence events detected will be proportional to the total amount of tracer contained in the volume of response. In the two-dimensional case, we consider only 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). For a general overview of standard PET image reconstruction methodology, see e.g. [1, 2] or [3].

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 recently introduced. Since MLEM or OSEM estimation of pixels or voxels is usually noisy [1], 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 (see e.g. [4], [5]). The Gaussian mixture model (GMM) is a well-known and widely used model in a variety of segmentation and classification problems ([6],[7]).

In this paper, we propose a new robust algorithm for efficient reconstruction of an object from a PET image. Our method utilizes a novel $L_{1}$ minimizing algorithm to estimate the parameters of a Gaussian mixture model within framework similar to the standard expectation-maximization (EM) algorithm. Instead of voxels, we estimate parameters of the object’s model. We focus on the two-dimensional case for clarity, but an extension to three dimensions would follow in a similar fashion. This paper is organized into six sections. Section II gives an overview of traditional methodology in 2D PET imaging. In Section  III an introduction to Gaussian mixture models is presented. Section IV describes the estimation of GMM parameters from PET data. Section V shows the implementation of iterative estimates in the EM-like algorithm. Finally, experimental results with discussion conclude the paper.

II Two-Dimensional PET Imaging

In the two-dimensional case, the acquired data are collected along LORs through a two-dimensional object. Traditionally, data are organized into sets of projections, integrals along all LORs for a fixed direction $\phi$. The collection of all projections for $0\leq\phi<2\pi$ forms a two dimensional function of the distance of the LOR from the origin, denoted by $s$, and angle $\phi$. The line-integral transform of $f(x,y)\to p(s,\phi)$ is called the X-ray transform [8], which in 2D is the same as the Radon transform [9]. Since a fixed point traces a sinusoidal path in the projection space, the superposition of all sinusoids corresponding to each point of activity in a general object is called a sinogram. This is illustrated in Fig. 1.

A classical image reconstruction method is the filtered backprojection (FBP), which is based on the projection theorem of Fourier analysis. The algorithm reconstructs the image by calculating the inverse Fourier transform of the 2D Fourier transform of the backprojected image (see e.g. [2], [8], [10]). This technique does 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.

III Gaussian Mixture Models

A Gaussian mixture model [11] is a weighted sum of $K$ component Gaussian densities as given by the equation:

[TABLE]

where $\bm{x}$ is a $d$-dimensional observation, $\tau_{i}~{}(i=1,...,K)$ are the weights of each Gaussian component, and $g(\bm{x}|\bm{\mu}_{k},\bm{\Sigma}_{k})$ are the Gaussian component densities. We assume that a given observation $\bm{x}$ is a realization from exactly one of the $K$ Gaussian mixture components, and each component density is a $d$-variate Gaussian function of the form:

[TABLE]

with $|\bm{\Sigma}_{k}|$ denoting the determinant of $\bm{\Sigma}_{k}$. The set of probabilities $\{\tau_{k}\}$ such that $\sum_{k=1}^{K}\tau_{k}=1$ defines the probabilities that $\bm{x}$ belongs to the corresponding Gaussian component. The complete Gaussian mixture model is parameterized by the mean vectors $\{\bm{\mu}_{k}\}$, covariance matrices $\{\bm{\Sigma}_{k}\}$ and mixture weights $\{\tau_{k}\}$ for all Gaussian component densities.

Traditionally, in image segmentation GMMs are used to model values at points of observation, e.g. activity concentration in voxels [12] or pixel values [13]. Those models do not take into account the spatial correlation between observations, which can be corrected by introducing Markov random fields [12, 13].

We propose applying the GMM directly to the spatial data, focusing only on the locations and not the values at those locations.

In this scenario, the $K$ Gaussian components represent potential points of origin for activity. The points $\bm{x}$ that are realizations from these components are latent (unobserved). Our observations are lines (LORs) through these points at random angles $\psi$, however using convenient properties of Gaussian distributions we will still be able to accurately estimate the parameters.

IV Estimating Gaussian Parameters

For clarity, in this section we focus only on one Gaussian component. Assuming we know a set of $N$ LORs whose points of origin come from a Gaussian distribution with parameters $(\bm{\mu},\bm{\Sigma})$, we can estimate those parameters.

IV-A Estimating $\bm{\mu}$ using minimal distance

The mean vector $\bm{\mu}=[\mu_{x}~{}\mu_{y}]^{T}$ can be estimated in multiple ways. One method is to find the point in space whose total squared distance from all LORs is minimal. Clearly, the solution depends on our definition of distance.

Each LOR in two dimensions is uniquely given by its slope $k=\tan\psi$ and an intercept $l$ (or by any two equivalent parameters) and we can write the LOR equations as

[TABLE]

In general, we can define the squared distance between $d$-dimensional vectors $\bm{v}_{1}$ and $\bm{v}_{2}$ as

[TABLE]

for some $d\times d$ weight matrix $\bm{W}$. In particular, when $\bm{W}=\bm{I}$ the distance in (4) is the Euclidian distance, and for $\bm{W}=\bm{\Sigma}^{-1}$ we get the Mahalanobis distance. In our case $d=2$ and a planar line is given by one equation, but similar results would follow for higher dimensions. For a given $\bm{\mu}$, the point on the $i$-th line that is nearest to $\bm{\mu}$, i.e. the solution to

[TABLE]

denoted by $\bm{x}_{i}^{m}$, can be found by solving [14]

[TABLE]

Now the point $\bm{\mu}$ that is nearest to all lines is the solution of

[TABLE]

where $\bm{x}_{i}^{m}$ is given in (5). This can be shown to be (see Appendix A) the solution to:

[TABLE]

where

[TABLE]

and

[TABLE]

Note that the zeros and null-vectors added in expressions above appear for calculation purposes (see Appendix A) and do not change the original problem.

It can also be shown that, for $\bm{\Sigma}$ known, the estimate obtained using the Mahalanobis distance is also the maximum likelihood estimate for $\bm{\mu}$.

IV-B Estimating $\bm{\Sigma}$ using 1D projections

In the two-dimensional setting, since $\bm{\Sigma}$ is a symmetric matrix,

[TABLE]

we need to estimate three parameters: $\Sigma_{11}$, $\Sigma_{12}$ and $\Sigma_{22}$.

Observe a single line of response, and recall that one of the parameters determining the LOR is the slope $\tan\psi$, where $\psi$ is the angle between the line and the $x$-axis. That means that by rotating the $xy$-coordinate system by $\psi-\frac{\pi}{2}$, the LOR would be parallel to the new $y$-axis. This is equivalent to the rotation of the Gaussian distribution by $\varphi=\frac{\pi}{2}-\psi$. This is illustrated in Fig. 2. Since a rotated Gaussian distribution is again Gaussian, in this new coordinate system the LOR is parallel to the $y$-axis, and the Gaussian distribution has parameters

[TABLE]

Given the original LOR parameters, the coordinates of the point at which it intersects the new $x$-axis are $(-l\sin\varphi,0)$. Since marginal distributions of a Gaussian are again Gaussian, a 1D projection onto the new $x$-axis is a Gaussian random variable with expectation

[TABLE]

and variance

[TABLE]

i.e. we obtain the one-dimensional mean and variance simply by omitting rows and columns from their multidimensional counterparts.

Therefore, each LOR gives us a one-dimensional projection whose squared (Euclidian) distance from the mean, $(\cos\varphi\mu_{x}-\sin\varphi\mu_{y}+l\sin\varphi)^{2}$, can be used to estimate the variance in (8). This gives us a system of equations:

[TABLE]

[TABLE]

Since there are (dozens of) thousands of measurements, the problem in (9) is overdetermined. It does not have an exact solution, but we can find the best approximation, i.e. the solution to

[TABLE]

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

[TABLE]

As we will show in Sec. VI, OLS performs poorly when we introduce more than one component, so we will need to use alternative methods. It can be shown that in some cases $L_{1}$ minimization is preferred to the more traditional $L_{2}$ minimization because it is less sensitive, i.e. more resistant to gross and systematic errors [15]. The main argument against $L_{1}$ minimisation would be computational complexity, which can be alleviated by using iterative methods. In this paper, the solution is obtained by using the $L_{1}$ minimization algorithm proposed in [16]. We use it to solve a modified problem

[TABLE]

where $k$ is a constant corrective scale factor. In a sufficiently large sample, one would have enough data points in each direction $\phi$ 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 $\sigma$ [17], i.e.

[TABLE]

where $\Phi$ denotes the distribution function of the standard normal random variable.

V Iterative EM-like Algorithm

The expectation-maximization (EM) algorithm, first explained in [18], is an iterative method used to find maximum likelihood or maximum a posteriori (MAP) estimates of parameters in statistical models. The expectation (E) step creates a function for the expectation of the log-likelihood evaluated using the current estimate of parameters. The maximization (M) step then computes parameters maximizing the expected log-likelihood from the E step. Conversely, these estimates are then used in the next E step. In traditional applications of the EM algorithm to GMMs, the E step assigns each data point its membership probabilities $\{\tau_{k}\}$, i.e. the probabilities that the point belongs to each of the mixture components. In the M step, parameters of each component are estimated using the points ”belonging” to that component. As already mentioned, in the PET setting observations are lines, however we can replicate the iterative steps using estimates from Sec. IV. Alternatively, this could also be considered a Lloyd-like algorithm, where the difference from the conventional Lloyd’s algorithm [19] is that we allow different distance functions of the form (4). The algorithm initializes parameters arbitrarily, and then alternates between the following steps:

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

Estimate component parameters. Either from a hard or soft classification, where each LOR participates with its proportional share, parameters of each component are estimated using methodology from Sec. IV.

The $L_{1}$ minimization algorithm recursively reduces and increases dimensionality of the observed subspace and uses weighted median to efficiently find the global minimum, and has shown to overperform state-of-the-art competitive methods when there are relatively few parameters to be estimated from a very high number of equations. For details, see Appendix B and [16].

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.

[TABLE]

where $\alpha$ increases from [math] to $1$. Therefore, in later iterations we obtain an MLE-like estimate of $\bm{\mu}$.

VI Results and remarks

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 Sec. IV provides good estimates with both $L_{1}$ and $L_{2}$ minimization, for several covariance matrices with varying corresponding correlation coefficients:

[TABLE]

Since the 2D covariance matrix is of the form

[TABLE]

it is determined by three parameters – $\sigma_{x}$, $\sigma_{y}$ and $\rho(=\rho_{xy})$. This corresponds to the three-dimensional vector

[TABLE]

in (9). Instead of observing true and estimated $\bm{\Sigma}$ matrices we will calculate the error in estimation of $\bm{s}$ for each $\bm{s}$ that corresponds to matrices above:

[TABLE]

Note that the corresponding correlation coefficients can be calculated from these vectors, and they are $\rho_{1}=0$, $\rho_{2}\approx-0.3$, $\rho_{3}\approx 0.9$.

For each of these types we simulated a measurement from $N=1000$ and $N=10000$ points and calculated the average estimate for $L_{1}$ and $L_{2}$ estimates separately. We repeated the experiment $1000$. Accuracy of an estimate can be assessed in many ways, for illustrative purposes we chose relative error i.e.

[TABLE]

where $\|\cdot\|$ denotes the standard Euclidian norm in both expressions. Results are given in Table I.

Given that the variance of the traditional standard deviation estimator (from points) equals $\frac{\sigma^{4}}{n-1}$, estimations within $10\%$ from $N=1000$ LORs seem acceptable, which justifies the methodology described in Sec. IV. Significantly more accurate estimation from $N=10000$ LORs further confirms this. It is also notable that the accuracy of the estimate increases as $|\rho||$ increases from [math] to $1$.

For $K=2$ components we repeated the experiment as described in Sec. V for various combinations of types of vector $\bm{s}$. We used hard classification, where we assign each line to at most one component. We experimented with various constraints, from simply assigning LORs 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 $L_{1}$ minimization and the scaling factor $k$ correcting the bias. The $L_{2}$ minimization method proved inefficient, since wrongly assigned LORs would cause unstable estimations and ”breaks” in the algorithm. An illustration of the results for $K=2$, $N=4000$ is shown in Fig. 3, along with the corresponding classical FBP method.

At the end, we would like to draw the readers attention to the fact that the reconstructed image is given by its parametric model: mean vectors $\{\bm{\mu}_{k}\}$, covariance matrices $\{\bm{\Sigma}_{k}\}$ and mixture weights $\{\tau_{k}\}$. 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 the successful object representation. Hence, future research will be oriented to compressed sensing approach ([20, 7, 21]) for reducing the number of projections, in this case reduced radiotracer’s concentration. Due to robustness of the proposed reconstruction method, post-filtering step is not needed.

Appendix A

We will find the solution to (6) using classical matrix calculus techniques. First note that the solution to (5) is

[TABLE]

To accommodate this, we will denote the expression in (6) by $d^{2}$ and expand it:

[TABLE]

For simplicity, denote

[TABLE]

Now $d^{2}$ equals

[TABLE]

which we will differentiate piecewise to find the minimum. We have:

[TABLE]

where

[TABLE]

• 2)

[TABLE]

• 3)

[TABLE]

• 4)

[TABLE]

Note that in all calculations we use the fact that $\bm{W}$ and, by extension, $\widetilde{\bm{W}}$ are symmetric. By plugging these equations into $\frac{d}{d\bm{\mu}}$ and equating that with [math] to obtain the minimum, we get:

[TABLE]

Define

[TABLE]

From the previous equation we now have

[TABLE]

from which it follows that

[TABLE]

It remains to verify that this stationary point $\bm{\mu}$ is also a turning point. However, since a point whose sum of squared distances from all lines is minimal must exist from a geometrical perspective, the solution in (7) is indeed the (global) minimum.

Appendix B

The $L_{1}$ minimization algorithm recursively reduces and increases dimensionality of the observed subspace and uses weighted median to efficiently find the global minimum. Reduction of dimensionality is achieved by extracting of parameters and inserting them into remaining equations in (11). If $[A_{i1}~{}A_{i2}~{}A_{i3}]$ is the $i$-th row of matrix $\bm{A}$, and $b_{i}$ the $i$-th element of vector $k\cdot\bm{b}$, the $i$-th equation of the system is

[TABLE]

We choose equation $j_{1}$ from the set $i=1,...,N$ and extract one of its parameters, e.g. $\Sigma_{11}$:

[TABLE]

We insert it into all other equations and get a new system with only two unknown parameters:

[TABLE]

for $i\neq j_{1}.$ Now, we choose some other equation $j_{2},j_{2}\neq j_{1}$ and extract one of the remaining parameters, e.g. $\Sigma_{12}$:

[TABLE]

We insert $\Sigma_{12}$ into all other equations and get the system of $N-2$ equations with only one unknown parameter $\Sigma_{22}$:

[TABLE]

where $i\neq j_{1},j_{2}$,

[TABLE]

[TABLE]

Parameter $\Sigma_{22}$ is calculated by minimizing $L_{1}$ norm

[TABLE]

The value of parameter $\Sigma_{22}$ is given by the weighted median (MED):

[TABLE]

where $j_{3}$ is an ordinal number of the concomitant equation and $\Diamond$ is the replication operator. The weighted median can be obtained using the algorithm given in [16]. The value of parameter $\Sigma_{22}$ is an element of set $\{b^{(1)}_{i}/A^{(1)}_{i}\}$. Chosen equations $\{j_{1},j_{2},j_{3}\}$ define a local minimum in 1D, a vertex.

Return to a higher dimension is achieved by putting calculated parameter value $\Sigma_{22}$ in (13). To further descending in $L_{1}$ cost surface, we fix equations $j_{1}$ and $j_{3}$, and try to find new $j_{4}$ using (14)-(16). If $j_{4}\neq j_{2}$, new vertex is defined by $\{j_{1},j_{2},j_{3}\}=\{j_{1},j_{3},j_{4}\}$ and we repeat the same procedure. If $j_{4}=j_{2}$, we conclude that the vertex $\{j_{1},j_{2},j_{3}\}$ is a local minimum in observed 2D subspace, thus we return to 3D by fixing $j_{2}$ and $j_{3}$ and try to find new $j_{4}$ instead of $j_{1}$. We repeat the previous procedure until we cannot find new equation. Since the $L_{1}$ cost surface is convex, the global minimum is reached. Parameters $\Sigma_{11}$, $\Sigma_{12}$ and $\Sigma_{22}$ are given by (12), (13) and (16).