Reconstruction

Here we consider a 3D PET MLEM reconstruction algorithm, the derivation of which comes from these slides, based on the IAEA Nuclear Medicine Physics handbook.

We consider list-mode reconstruction, meaning that the input will be the measured lines of response (LORs), each defined by two interaction points, \((x_1,y_1,z_t,t_1)\) and \((x_2,y_2,z_2,t_2)\), and we seek to reconstruct a 3D image in a volume divided into \(N\) 3D voxels. To derive the MLEM procedure, we will consider several quantities. Note that in the following definitions, index \(i\) will run over the lines of response and index \(j\) will run over the voxels of the reconstructed image:

  • \(\lambda_j\): the value of the reconstructed image at voxel \(j\)

  • \(A_{ij}\): the probability that the emission of electron-positron annihilation radiation within voxel \(j\) will lead to detection in LOR \(i\) (referred to as the “system matrix”)

  • \(y_i\): the measured number of counts in LOR \(i\)

In practice, since LORs are not binned in any way (we consider a list of all the LORs observed), \(y_i = 1\) for each measured LOR and all other LORs (theoretically an infinite number) will have \(y_i = 0\).

We now consider a theoretical construction \(x_{ij}\), which contains the number of measured counts in LOR \(i\) that were emitted from voxel \(j\). We model these as Poisson distributed about a mean value \(\bar{x}_{ij} = A_{ij}\lambda_j\). The log-likelihood for observation of the image with voxels \(\lambda_j\) given the measured \(x_{ij}\) is

\[L(\lambda) = \sum_{ij}x_{ij}\ln(A_{ij}\lambda_j) - A_{ij}\lambda_j.\]

To maximize the likelihood, we look for \(\lambda_{j'}\) such that \(\partial L/\partial\lambda_{j'} = 0\).

\[\frac{\partial L}{\partial\lambda_{j'}} = \sum_{ij}\Bigl(\frac{x_{ij}}{\lambda_j}\delta_{jj'} - A_{ij}\delta_{jj'}\Bigr) = \frac{\sum_{i}x_{ij'}}{\lambda_j'} - \sum_{i}A_{ij'} = 0.\]

Therefore we require that for each \(\lambda_j\)

\[\lambda_j = \frac{\sum_{i}x_{ij}}{\sum_{i}A_{ij}}.\]

This will be solved iteratively, denoting the voxel values at step \(k\) as \(\lambda^{(k)}_j\) starting with an initial guess \(\lambda^{(0)}_j\) and solving for each successive \(\lambda^{(k+1)}_j\) given the measurements \(x^{(k)}_{ij}\) at step \(k\),

\[\lambda^{(k+1)}_j = \frac{\sum_{i}x^{(k)}_{ij}}{\sum_{i}A_{ij}}.\]

Since we don’t have all the information in \(x_{ij}\), we must now express it in terms of actual measurements that have been made. If we just use the expected value \(\langle x^{(k)}_{ij}\rangle = A_{ij}\lambda^{(k)}_j\), we will get \(\lambda^{(k+1)}_j = \lambda^{(k)}_j\). To actually incorporate the measured LORs, we write the measured number of counts in each LOR \(y_{i}\) as

\[y_{i} = \sum_{j}A_{ij}\lambda^{(k)}_{j} + b_{i},\]

where the term \(b_{i}\) contains additional contributions to the LOR, such as counts due to random scatters. Then we divide \(y_i\) by the entire right side of the equation above, writing

\[\frac{y_{i}}{\sum_{j}A_{ij}\lambda^{(k)}_{j} + b_{i}} = 1,\]

and we multiply \(\langle x^{(k)}_{ij}\rangle\) by 1,

\[\langle x^{(k)}_{ij}\rangle \rightarrow \langle x^{(k)}_{ij}\rangle\times 1 = A_{ij}\lambda^{(k)}_j \times \frac{y_{i}}{\sum_{j}A_{ij}\lambda^{(k)}_{j} + b_{i}}.\]

Now we can write the iterative equation for \(\lambda^{(k+1)}_j\) with the above expression for the expected value \(\langle x^{(k)}_{ij}\rangle\),

\[\lambda^{(k+1)}_j = \frac{\lambda^{(k)}_j}{\sum_{i}A_{ij}}\times \sum_{i}\Biggl(\frac{y_{i}A_{ij}}{\sum_{j}A_{ij}\lambda^{(k)}_{j} + b_{i}}\Biggr)\]

Since we will consider each LOR individually, we have for measured lines of response \(i_{m}\), \(y_{i_m} = 1\) and for all others \(y_{i} = 0\), so we can rewrite the sum over \(i\) as the sum over measured LORs \(i_m\),

\[\lambda^{(k+1)}_j = \frac{\lambda^{(k)}_j}{\sum_{i}A_{ij}}\times \sum_{i_{m}}\Biggl(\frac{A_{i_{m}j}}{\sum_{j}A_{i_{m}j}\lambda^{(k)}_{j} + b_{i_{m}}}\Biggr).\]

This is the final equation that must be implemented in the code. Some notes on the content of this equation:

  • \(A_{i_{m}j}\) will be calculated on-the-fly for each measured LOR \(i_{m}\). Note that this is done by considering the radiological path of the line of response through the active volume, that is, the set of active voxels encountered by the LOR from \((x_1,y_1,z_t)\) to \((x_2,y_2,z_2)\). These voxels can be determined using the Siddon algorithm [R. Siddon. Med. Phys. 12 (2), 252 (1985)]. Furthermore, the voxels can be weighted using the time of flight information, such that rather than considering the path as a uniform line in which the probability that the electron-positron emission came from any given voxel on the line is equal, the distribution can be a gaussian centered on the voxel corresponding to the recorded time-of-flight.

  • \(P_{i_{m}} \equiv \sum_{j}A_{i_{m}j}\lambda^{(k)}_{j} + b_{i_{m}}\) can then be calculated using the present image voxels \(\lambda^{(k)}_{j}\) and a given model for the additive contributions \(b_{i_{m}}\). Because \(P_{i_{m}}\) is a sum over contributions to LOR \(i_{m}\) from all voxels, it is called the projection (or forward projection) onto \(i_{m}\).

  • The sum \(BP_{j} = \sum_{i_{m}}(A_{i_{m}j}/P_{i_{m}})\) distributes the measured LORs back onto the image according to the probabilities \(A_{i_{m}j}\) that the observed LOR was due to emission from voxel \(j\). This is called the back projection.

  • The sum \(S_{j} \equiv \sum_{i}A_{ij}\) is the total probability of observing any LOR emitted from voxel \(j\). This is called the sensitivity matrix, as it describes how “sensitive” the detector is to emission from a specific location in the active region.

One iteration of the algorithm is then performed as follows.

  1. For each LOR \(i_{m}\):

    1. Compute the projection \(P_{i_{m}}\) of the current image.

    2. Compute the contribution to the backprojection \(A_{i_{m}j}/P_{i_{m}}\) and add it to the backprojection \(BP_{j}\) of voxel \(j\).

  2. Compute the \(\lambda^{(k+1)}_{j}\) by multiplying each value \(BP_{j}\) by the current voxel divided by the corresponding value in the sensitivity matrix, \(\lambda^{(k)}_{j}/S_{j}\).

A helper class MLEMReconstructor has been created to call the reconstruction algorithm, implemented in C, from Python.

class mlem.mlem_reconstruct.MLEMReconstructor(prefix='mlem', niterations=1, save_every=- 1, TOF=True, TOF_resolution=200.0, img_size_xy=180.0, img_size_z=180.0, img_nvoxels_xy=60, img_nvoxels_z=60, smatrix=None, libpath='lib/libMLEM.so')[source]

Python wrapper class to perform MLEM reconstruction. The class provides a method to call the C-based reconstruction using the configuration provided by the class variables.

The reconstructed image is stored in a 1D array of 32-bit (4-byte) floats with values corresponding to:

(x0,y0,z0), (x1,y0,z0) ... (xN,y0,z0)
(x0,y1,z0), (x1,y1,z0) ... (xN,y1,z0)
...
(x0,yN,z0), (x1,yN,z0) ... (xN,yN,z0)
(x0,y0,z1), (x1,y0,z1) ... (xN,y0,z1)
(x0,y1,z1), (x1,y1,z1) ... (xN,y1,z1)
...
(x0,yN,z1), (x1,yN,z1) ... (xN,yN,z1)

...

(x0,y0,zN), (x1,y0,zN) ... (xN,y0,zN)
(x0,y1,zN), (x1,y1,zN) ... (xN,y1,zN)
...
(x0,yN,zN), (x1,yN,zN) ... (xN,yN,zN)

To ensure the correct size of the sensitivity matrix, the class has to be initialized with the values of img_nvoxels_xy and img_nvoxels_z in the constructor method.

NOTE: the LOR points may need to be sorted for the reconstruction to be correct

Parameters

  • prefix (str) – the filename prefix for all saved files

  • niterations (int) – the number of iterations to perform on reconstruction

  • save_every (int) – save every specified number of iterations

  • TOF (bool) – boolean to enable TOF

  • TOF_resolution (float) – the TOF resolution in ps

  • img_size_xy (float) – the image size in the x and y dimensions (in mm)

  • img_size_z (float) – the image size in the z dimension (in mm)

  • img_nvoxels_xy (int) – the number of voxels in the x and y dimensions

  • img_nvoxels_z (int) – the number of voxels in the z dimension

  • smatrix (Optional[ndarray]) – the sensitivity matrix s[x,y,z]

  • libpath (str) – the path to the C++ reconstruction library

read_image(niter)[source]

Reads a reconstructed image stored as “(prefix)(niter).raw”.

Parameters

niter (int) – the iteration number to be read

Return type

ndarray

reconstruct(lor_x1, lor_y1, lor_z1, lor_t1, lor_x2, lor_y2, lor_z2, lor_t2)[source]

Performs reconstruction by calling a C-implemented function. The reconstruction is list-mode, so the inputs are the coordinates \((x,y,z,t)\) for each line of response (LOR). These coordinates are specified as separate lists.

Parameters

  • lor_x1 (list) – x-coordinates for the first point in the LOR

  • lor_y1 (list) – y-coordinates for the first point in the LOR

  • lor_z1 (list) – y-coordinates for the first point in the LOR

  • lor_t1 (list) – time coordinates for the first point in the LOR

  • lor_x2 (list) – x-coordinates for the second point in the LOR

  • lor_y2 (list) – y-coordinates for the second point in the LOR

  • lor_z2 (list) – y-coordinates for the second point in the LOR

  • lor_t2 (list) – time coordinates for the second point in the LOR

Return type

ndarray

Returns

array of shape [img_size_xy, img_size_xy,`img_size_z`] containing the reconstructed image