Appendix Two

Expectation of Poisson data contributing to a measurement¶

Assume the configuration of Figure 12: two radioactive sources contribute to a single measurement $N$. We know the a-priori expected values $\bar{a}$ and $\bar{b}$ for each of the sources, and we know the measured count $N = a + b$. The question is to compute the expected values of $a$ and $b$ given $N$.

By definition, the expected value equals:

The denominator should be equal to 1, but if we keep it, we can apply the equation also if $p$ is known except for a constant factor.

The prior expectations for $a$ and $b$ are:

After the measurement, we know that the first detector can only have produced $a$ counts if the other one happened to contribute $N-a$ counts. Dropping some constants this yields:

Applying (1) yields:

Let us first look at the denominator:

A bit more work needs to be done for the numerator:

Combining numerator and denominator results in:

The convergence of the EM algorithm¶

This section explains why maximizing the likelihood of the complete variables $L_x$ is equivalent to maximizing the likelihood $L$ of the observed variables $Q$.

First some notations and definitions must be introduced. As in the text, Λ denotes the reconstruction, $Q$ the measured sinogram and $X$ the complete variables. We want to maximize the logarithm of the likelihood given the reconstruction:

The conditional likelihood of the complete variables $f(X | \Lambda)$ can be written as:

where $k(X | Q, \Lambda)$ is the conditional likelihood of the complete variables, given the reconstruction and the measurement. These definitions immediately imply that

The objective function we construct during the E-step is defined as

which means that we write the log-likelihood of the complete variables as a function of $\Lambda'$, and that we eliminate the unknown variables $X$ by computing the expectation based on the current reconstruction Λ. Combining (10) and (11) results in

Finally, we define a generalized EM (GEM) algorithm. This is a procedure which computes a new reconstruction from the current one. The procedure will be denoted as $M$. $M$ is a GEM-algorithm if

This means that we want $M$ to increase $h$. We can be more demanding and require that $M$ maximizes $h$; then we have a regular EM algorithm, such as the MLEM algorithm of section Iterative Reconstruction.

Now, from equation (12) we can compute what happens with $L$ if we apply a GEM-step to increase the value of $h$:

Because $M$ is a GEM-algorithm, we already know that $h(M(\Lambda) | \Lambda) - h(\Lambda, \Lambda)$ is positive. If we can also show that

then we have proven that every GEM-step increases the likelihood $L$.\

By definition, we have

Therefore, the left hand side of equation (15) can be rewritten as

which can be written as

with the additional requirement that $\int k(X | Q, \Lambda) dX = \int k(X |Q, M(\Lambda)) dX = 1$. It turns out that (18) is always positive, due to the convexity of $t \ln t$. We will have a look at that now.\

Consider the function $\psi(t) = t \ln t$. It is only defined for $t > 0$. Here are its derivatives:

The second derivative is continuous and strictly positive ($t \ln t$ is a convex function). Consequently, we can always find a value $u$ such that

Because $\psi(1) = 0$ and $\psi'(1) = 1$, this becomes:

Consider now an integral of the same form as in (18):

We can rework it such that we can exploit the convexity of $t \ln t$:

Now we have proved that the GEM-step increases $L$. It still remains to be shown that the algorithm indeed converges towards the maximum likelihood solution. If you want to know really everything about it, you should read the paper by Salvo and Defrise, “A Convergence Proof of MLEM and MLEM-3 With Fixed Background”, IEEE TMI, 38 (2019).

Backprojection of a projection of a point source¶

2D parallel projection¶

Consider a point source located in the center of the field of view. The projection is a sinogram $q(s, \theta)$, given by

where $\delta(x)$ is the delta function: $\delta(x) = 1$ if $x = 0$, and $\delta(x) = 0$ if $x \neq 0$. The backprojection $b(x,y)$ is then given by

To proceed, the argument of the delta function must be changed. This can be done with the following expression:

where $f'$ is the derivative of $f$, and $x_n, n=1 \ldots N$ are the zeros of $f(x)$, i.e. $f(x_n) = 0$. A simple example is $\delta(ax) = \delta(x) / |a|$. Applying this to (26) yields:

where $\theta_0$ is such that $x\cos\theta_0 + y\sin\theta_0 = 0$. This is satisfied if

which results in

2D parallel projection with TOF¶

Consider again the projection of a point source in the center of the field of view. With time-of-flight, the parallel projection along angle θ can be considered as a Gaussian blurring along lines with angle θ. That gives for the TOF-sinogram $q(x,y,\theta)$:

where σ represents the uncertainty of the TOF-measurement. The corresponding backprojection operator is obtained by applying the same blurring to the projections, followed by integration over θ. The convolution of a Gaussian with itself is equivalent to multiplying its standard deviation with $\sqrt{2}$ (see appendix Combining resolution effects: convolution of two Gaussians). One finds

which is equivalent to equation (41).