TUTORIAL ON USING THE NIPET PACKAGE FOR PROCESSING THE LIST-MODE DATA

Download raw PET data

Example phantom data acquisition can be downloaded from this link: The list-mode data of phantom acquisition

Under the above link should be four files:

1. Two norm files (*.dcm, *.bf)
2. Two list-mode files (*.dcm, *.bf)

Download those files to folder used for storing and accessing the PET data.

Import required Python packages

import matplotlib.pyplot as plt
import numpy as np

also import the PET package, nipet:

import nipet

Initialisation of all scanner constants

Initialise all the mMR scanner constants, the transaxial and axial look-up tables:

Cnt, txLUT, axLUT = nipet.mmraux.mmrinit()

Provide the path to folder:

folderin = '/path/to/rawpet'

Then automatically read and recognise the input data by calling:

datain = nipet.mmraux.explore_input(folderin, False)

Processing the list-mode data

For histogramming the data list-mode data, you can choose your own time framing. The information is stored in Numpy array of 16-bit unsigned integers describing the duration of each consecutive frame. If only one element of '0' is present that means that the data will be histogrammed into single static sinogram, i.e.:

frames = np.array([0],  dtype=np.uint16)

The histogramming can be done with or without bootstrapping, depending on integer value of btp, i.e., btp = 0: normal histogramming; btp = 1: non-parametric bootstrapping; btp = 2: parametric bootstrapping (recommended over the non-parametric). Let's choose histogramming without bootstrapping:

btp = 0

The start and stop time of histogramming can also be specified apart from the above time frames. If Tstart and Tstop are qual, than the whole dataset is histogrammed.

Tstart = 0
Tstop  = 0

Then the histogramming is called by:

hst = nipet.lm.mmrhist.hist(datain, txLUT, axLUT, Cnt, frms=frames, nbtp=btp, tstart=Tstart, tstop=Tstop)

The above is the default histogramming and can be called without the optional variables:

hst = nipet.lm.mmrhist.hist(datain, txLUT, axLUT, Cnt)

It is important to check if the normalisation file with all the necessary normalisation components is up-to-date relative to the list-mode data (the normalisation file should be from the same day QC acquisition as the PET data). The time difference can be checked by this function:

nipet.mmraux.time_diff_norm_acq(datain)

The noise-reduced random events sinogram

The random event sinogram with reduced variance can be found by using the fan sums of delayed events outputted during the above histogramming, i.e.:

rsino, cmap = nipet.lm.mmrhist.rand(hst['fansums'], txLUT, axLUT, Cnt)

where rsino is the estimated random event sinogram and cmap is the crystal map of estimated single rates.

Full normalisation sinogram (including detector dead-time)

The normalisation sinogram (which includes the correction for detector dead-time) can be obtained using this function:

nsino = nipet.mmrnorm.get_sino(datain, hst, axLUT, txLUT, Cnt)

Plotting the processed list-mode data

The head curve (total counts per second for prompt and delayed event data) can be plotted this way:

plt.figure()
# plot prompt head curve
plt.plot(hst['phc'])
# plot delayeds head curve
plt.plot(hst['dhc'])
plt.show()

Any sinogram has the following dimensions: [sinoIdx, angleIdx, binIdx]. For span-11 they are [837, 252, 344]. For span-1 maximum sinoIdx = 4084 instead of 837. The prompt sinogram for sinogram #63 (this is a direct sinogram, roughly in the middle of the axial field of view) can be viewed without interpolation using this function:

plt.matshow(hst['psino'][63,:,:])
plt.show()

for delayed events:

plt.matshow(hst['dsino'][63,:,:])
plt.show()

which can be compared with the estimated random sinogram:

plt.matshow(rsino[63,:,:])
plt.show()

The prompt single slice rebinned sinograms can also be viewed:

plt.matshow(hst['pssr'][63,:,:])
plt.show()