Appendix One

Poisson noise¶

Assume that the expected amount of measured photons per unit time equals $r$. Let us now divide this unit of time in a $k$ time intervals. If $k$ is sufficiently large, the time intervals are small, so the probability of detecting a photon in such an interval is small, and the probability of detecting two or more is negligible. Consequently, we can assume that in every possible measurement, only zero or one photon is detected in every interval. The probability of detecting one photon in an interval is then {math}r / k. A measurement of $n$ photons must consist of $n$ intervals with one photon, and $k - n$ intervals with no photons. The probability of such a measurement is:

The first factor is the probability of detecting $n$ photons in $n$ intervals. The second factor is the probability of detecting no photons in $n - k$ intervals. The third factor is the amount of ways that these $n$ successful and $n-k$ unsuccessful intervals can be combined in different measurements.

As mentioned above, equation (1) becomes better when $k$ is larger. To make it exact, we simply have to let $k$ go to infinity.

It turns out that computing the limit of the logarithm is easier, so we have

So it follows that

Convolution and PSF¶

In section Mechanical collimation in the gamma camera, the collimator point spread function (PSF) was computed. The collimator PSF tells us which image is obtained for a point source at distance $H$ from the collimator. What happens if two point sources are positioned in front of the camera, both at the same distance $H$? Since the sources and the photons don’t interact with each other, all what was said above still applies, for each of the sources. The resulting image will consist of two PSFs, each centered at the detector point closest to the point source. Where the PSFs overlap, they must be added, since the detector area in the overlap region gets photons from both sources. The same is true for three, four, or one million point sources, all located at the same distance $H$ from the collimator. To predict the image for a set of one million point sources, simply calculate the corresponding PSFs centered at the corresponding positions on the detector, and sum everything.

The usual mathematical description of this can be considered as a two step approach:

1. Assume that the system is perfect: the image of a point source is a point, located on the perpendicular projection line through the point source. Mathematicians would call that “point” in the image a “Dirac impulse”. The image of two or more point sources contains simply two or more Dirac impulses, located at corresponding projection lines.

   Let $f(x,y)$ represent the image value at position $(x,y)$. This image can be regarded as the sum of an infinite number of Dirac impulses $\delta(x,y)$, one at every location $(x,y)$:

   $$\begin{align} f(x,y) &= (f \otimes \delta) (x,y)\\ &= \int_{-\infty}^{\infty} du \int_{-\infty}^{\infty} dv \, f(u,v) \, \delta(x-u, y-v) \end{align}$$

   (5)

2. Correction of the first step: the expected image of a point is not a Dirac impulse, but a PSF. Therefore, replace each of the Dirac impulses in the image by the corresponding PSF, centered at the position of the Dirac impulse.

   $$\begin{align} g(x,y) &= (f \otimes \mbox{PSF}) (x,y)\\ &= \int_{-\infty}^{\infty} du \int_{-\infty}^{\infty} dv \, f(u,v) \, \mbox{PSF}(x-u, y-v) \end{align}$$

   (6)

The second operation is called the convolution operation. Assume that a complex flat (e.g. a inhomogeneous radioactive sheet) tracer distribution would be put in front of the camera, parallel to the collimator (at distance $H$). What image will be obtained? Regard the distribution as consisting of an infinite number of point sources. This does not change the distribution, it only changes the way you look at it. Project all of these sources along an ideal perpendicular projection line into the image. You will now obtain an image consisting of an infinite number of Dirac impulses. Replace each of these impulses with the PSF and sum the overlapping values to obtain the expected image.

If the distance to the collimator is not the same for every point source, then things get more complicated. Indeed, the convolution operator assumes that the PSF is the same for all point sources. Therefore, to calculate the expected image, the previous procedure has to be applied individually for every distance, using the corresponding distance dependent PSF. Finally, all convolved images have to be summed to yield the expected image.

Combining resolution effects: convolution of two Gaussians¶

Very often, the PSF can be well approximated as a Gaussian. This fact comes in handy if we want to combine two PSFs. For example: what is the combined effect of the intrinsic resolution (PSF of scintillation detection) and the collimator resolution (collimator PSF)?

How can two PSFs be combined? The solution is given in appendix Convolution and PSF: one of the PSFs is regarded as a collection of Dirac impulses. The second PSF must be applied to each of these pulses. So we must compute the convolution of both PSFs. This appendix shows that if both are approximately Gaussian, the convolution is easy to compute.

Let us represent the first and second PSFs as follows:

Thus, $\sigma_1 = a / \sqrt{2}$ and $A = 1 / (\sqrt{2 \pi} \sigma_1)$, and similar for the other PSF. The convolution is then written as:

The exponentiation contains a term in $u^2$ and a term in $u$. We can get rid of the $u$ by requiring that both terms result from squaring something of the form $(u + C)^2$, and rewriting everything as $(u + C)^2 - C^2$. The $C^2$ is not a function of $u$, so it can be put in front of the integral sign.

The integrand is a Gaussian. The center is a function of $x$, but the standard deviation is not. The integral from $-\infty$ to ∞ of a Gaussian is a finite value, only dependent on its standard deviation. Consequently, the integral is not a function of $x$. Working out the factor in front of the integral sign and combining all constants in a new constant $D$, we obtain

So the convolution of two Gaussians is again a Gaussian. The variance of the resulting Gaussian is simply the sum of the input variances (by definition, the variance is the square of the standard deviation).

The FWHM of a Gaussian is proportional to the standard deviation, so we obtain a very simple expression to compute the FWHM resulting from the convolution of multiple PSFs:

Error propagation¶

The mean and the variance of a distribution $a$ are defined as:

where $E$ is the expectation, $a_i$ is a sample from the distribution and the summation is over the entire distribution. By definition, $\sigma_a$ is the standard deviation. When computations are performed on samples from a distribution (e.g. Poisson variables) the noise propagates into the result. Consequently, one can regard the result also as a sample from some distribution which can be characterized by its (usually unknown) mean value and its variance or standard deviation. We want to estimate the variance of that distribution to have an idea of the precision on the computed result. We will compute the variance of the sum or difference and of the product of two independent variables. We will also show how the variance on any function of independent samples can be easily estimated using a first order approximation.

Keep in mind that these derivations only hold for independent samples. If some of them are dependent, you should first eliminate them using the equations for the dependencies.

Sum or difference of two independent variables¶

We have two variables $a$ and $b$ with mean $\overline{a}$ and $\overline{b}$ and variance $\sigma_a^2$ and $\sigma_b^2$. We compute $a \pm b$ and estimate the corresponding variance $\sigma_{a \pm b}^2$.

because the expectation of $(a - \overline{a})$ is zero. The expectation $E\left[(a-\overline{a})(b-\overline{b})\right]$ is the covariance of $a$ and $b$. The expectation of the product is the product of the expectations if the variables are independent samples, and therefore, the covariance of independent variables is zero.

So in linear combinations the noise adds up, even if the variables are subtracted!

Product of two independent variables¶

For independent variables, the expectation of the product is the product of the expectations, so we have:

This expression is not very useful, it must be rewritten as a function of $\overline{a}$, $\overline{b}$, $\sigma_a$ and $\sigma_b$. To obtain that, we rewrite $a$ as $\overline{a} + (a - \overline{a})$:

Substituting this result for $E[a^2]$ and $E[b^2]$ in (16) we obtain:

The last line is a first order approximation which is acceptable if the relative errors are small. We conclude that when two variables are multiplied the relative variances must be added.

Any function of independent variables¶

If we can live with first order approximations, the variance of any function of one or more variables can be easily estimated. Consider a function $f(x_1, x_2, \ldots x_n)$ where the $x_i$ are independent samples from distributions characterized by $\overline{x_i}$ and $\sigma_i$. Applying first order equations means that $f$ is treated as a linear function:

The first step is a first order approximation: the expectation of a linear function is the function of the expectations. Similarly, the second line is a Taylor expansion, assuming all higher derivatives are zero. The third step is the application of (14).

With this approach you can easily verify that the variance on a product or division is obtained by summing the relative variances.

Central slice theorem¶

This appendix gives a proof of eq (62) for any angle θ. The projection $q(s,\theta)$ can be written as

The 1D Fourier transform of $q(s,\theta)$ along $s$ equals:

Now consider the 2D Fourier transform of $\lambda(x,y)$:

We now switch to a rotated coordinate system by setting:

It follows that $\nu_x x + \nu_y y = \nu_s s + \nu_r r$. The Jacobian of the transformation equals unity, so $dx dy = dr ds$. Thus, we obtain:

Comparison with (22) reveals that setting $\nu_r = 0$ in Λ produces $Q_1$:

A central profile through Λ along θ yields the Fourier transform of the projection $q(s,\theta)$.

The inverse Fourier transform of the ramp filter¶

To compute the inverse Fourier transform of the ramp filter, it is useful to consider it as the difference between a rectangular and a triangular filter, as illustrated in Figure 1. In practical implementations, the ramp filter must be broken off at some point, which we will call $W$ in the following. The highest value of $W$ in a discrete implementation is the Nyquist frequency, which is the highest frequency that can be represented with pixels. It equals 0.5, meaning that its period is 0.5 pixels long: a (co)sine at the Nyquist frequence has a value of 1 in one pixel and a value of -1 in the next.

Consider the rectangular filter ${\cal R}(\omega)$ of Figure 1: its value equals $W$ for frequency ω between $-W$ and $W$, and it is zero elsewhere. Its inverse Fourier transform equals:

The third equality follows from the fact that the integral of a sine over an interval symmetrical about zero is always zero. The inverse Fourier of a rectangular filter is a so-called sinc function.

The convolution of a rectangular function with itself is a triangular function. More precisely, consider the rectangular function ${\cal R}_2(\omega)$ and the triangular function ${\cal T}(\omega)$ defined as:

Then it is easy to verify that the convolution of ${\cal R}_2(\omega)$ with itself equals ${\cal T}(\omega)$:

Note that the rectangular filters ${\cal R}_2(\omega)$ and ${\cal R}(\omega)$ have a different width and amplitude.

Convolution in the Fourier transform corresponds to a product in the spatial domain, so the inverse Fourier transform of a triangular filter is the square of a sinc function. Therefore, the inverse Fourier transform of the difference between the rectangular and triangular filters, i.e. the ramp filter, equals:

Figure 2 shows a plot of (31). In a discrete implementation, where one only needs the function values at integer positions $x = -N, -N+1, \ldots,N$, the highest possible value for $W$ is 0.5, the Nyquist frequency. The red dots in the figure show these discrete function values for $W = 0.5$.

The Laplace transform¶

The Laplace transform is defined as:

The Laplace transform is very useful in computations involving differential equations, because integrals and derivatives with respect to $t$ are transformed to very simple functions of $s$. Some of its interesting features are listed below (most are easy to prove). The functions of the time are at the left, the corresponding Laplace transforms are at the right:

By combining (34) and (35) one obtains