GE Discovery MI (Listmode Reconstruction; With Time of Flight)

This tutorial demonstrates how to reconstruct data listmode data in HDF5 data obtain from GE scanners. In this tutorial, we’ll reconstruct using time-of-flight listmode. The data can be obtained from https://zenodo.org/records/8404015

from pytomography.metadata import ObjectMeta
from pytomography.metadata.PET import PETLMProjMeta
from pytomography.projectors.PET import PETLMSystemMatrix
from pytomography.algorithms import OSEM, BSREM
from pytomography.priors import RelativeDifferencePrior
from pytomography.likelihoods import PoissonLogLikelihood
from pytomography.transforms.shared import GaussianFilter
import matplotlib.pyplot as plt
from pytomography.io.PET import clinical

We first need to get the required PET geometry information dictionary from the scanner name. Currently, we only support discovery_MI, but feel free to make a commit to support other scanners (add to the src/data/pet_scanner_info.txt file)

scanner_name = 'discovery_MI'
info = clinical.get_detector_info(scanner_name)

Now we obtain the required data from the listmode files.

# Get listmode events
detector_ids = clinical.get_detector_ids_hdf5('/disk1/ge_discovery_tof_mi_pet/dmi_nema_lm/LIST0000.BLF', scanner_name)
# Get scatter/random correction term at each event
additive_term = clinical.get_additive_term_hdf5('/disk1/ge_discovery_tof_mi_pet/dmi_nema_lm/corrections.h5')
# Get multiplicative weights (attenuation/normalization) at each event
weights = clinical.get_weights_hdf5('/disk1/ge_discovery_tof_mi_pet/dmi_nema_lm/corrections.h5')
# Get ALL valid detector-pairs and corresponding multiplicative weights
detector_ids_sensitivity, weights_sensitivity = clinical.get_sensitivity_ids_and_weights_hdf5('/disk1/ge_discovery_tof_mi_pet/dmi_nema_lm/corrections.h5', scanner_name)

First we’ll get the required object space and projection space metadata.

• The object_meta specifies the size of each voxel, and the number of voxels in each direction

• The proj_meta specifies the list of all measured listmode events (detector_ids), the time-of-flight metadata, and the detector IDs / weights used for computing the sensitivity image.

# Define object space reconstruction matrix
object_meta = ObjectMeta(
    dr=(2.78,2.78,2.78), #mm
    shape=(192,192,71) #voxels
)
# Get time-of-flight metadata and define projection space metadata
tof_meta = clinical.get_tof_meta(scanner_name)
proj_meta = PETLMProjMeta(
    detector_ids,
    info,
    tof_meta=tof_meta,
    detector_ids_sensitivity=detector_ids_sensitivity,
    weights_sensitivity=weights_sensitivity
    )

Now we’ll define the system matrix for the PET system. In this case we’ll use 4.5mm PSF modeling.

psf_transform = GaussianFilter(4.5) # 4.5mm gaussian psf
system_matrix = PETLMSystemMatrix(
    object_meta,
    proj_meta,
    obj2obj_transforms = [psf_transform],
    N_splits=8,
)

Now we’ll define the likelihood function (PoissonLogLikelihood for PET systems). The additive term contains scatters/randoms and needs to be normalized by the weights.

likelihood = PoissonLogLikelihood(
    system_matrix,
    additive_term=additive_term/weights,
)

We can now reconstruct using a variety of reconstruction algorithms:

recon_algorithm = OSEM(likelihood)
recon_OSEM = recon_algorithm(n_iters=2, n_subsets=34)

# This is probably very similar to "Q.Clear" is ;)
prior = RelativeDifferencePrior(beta=25, gamma=2)
recon_algorithm = BSREM(likelihood, prior=prior)
recon_BSREM = recon_algorithm(n_iters=20, n_subsets=34)

Plot axial slice of reconstruction:

fig, ax = plt.subplots(1,2,figsize=(5.5,3), gridspec_kw={'wspace': 0.0})
plt.subplot(121)
plt.imshow(recon_OSEM[:,:,20].T.cpu(), cmap='magma', interpolation='gaussian', vmax=0.5)
plt.axis('off')
plt.text(0.03, 0.97, 'OSEM', color='white', fontsize=15, fontweight='bold', transform=plt.gca().transAxes, ha='left', va='top')
plt.xlim(35,157)
plt.ylim(45,152)
plt.subplot(122)
plt.imshow(recon_BSREM[:,:,20].T.cpu(), cmap='magma', interpolation='gaussian', vmax=0.5)
plt.axis('off')
plt.text(0.03, 0.97, 'BSREM (RDP)', color='white', fontsize=15, fontweight='bold', transform=plt.gca().transAxes, ha='left', va='top')
plt.xlim(35,157)
plt.ylim(45,152)
fig.tight_layout()
plt.show()