Basic PET image reconstruction

For basic Siemens Biograph mMR image reconstruction, NiftyPET requires:

1. PET list-mode data;

2. component-based normalisation file(s);

3. the μ-map image (the linear attenuation map as a 3D image).

An example of downloadable raw PET data of amyloid brain scan is provided in Raw brain PET data [8].

Initialisation

Prior to dealing with the raw input data, required packages need to be imported and the Siemens Biograph mMR scanner constants and parameters (transaxial and axial lookup tables, LUTs) have to be loaded:

import sys, os, logging
import numpy as np

# NiftyPET image reconstruction package (nipet)
from niftypet import nipet
# NiftyPET image manipulation and analysis (nimpa)
from niftypet import nimpa

# control the logging output
logging.basicConfig(level=logging.INFO)

# get all the scanner parameters
mMRpars = nipet.get_mmrparams()

Sorting and classification of input data

All the part of the input data is aimed to be automatically recognised and sorted in Python dictionary. This can be obtained by providing a path to the folder containing the unzipped file of the freely provided raw PET data in Raw brain PET data. The data is then automatically explored and sorted in the output dictionary datain:

# Enter the path to the input data folder
folderin = '/path/to/input/data/folder/'

# automatically categorise the input data
datain = nipet.classify_input(folderin, mMRpars)

The output of datain for the above PET data should be as follows:

In [5]: datain
Out[5]:
{'#mumapDCM': 192,
 'corepath': '/data/amyloid_brain',
 'lm_bf': '/data/amyloid_brain/LM/17598013_1946_20150604155500.000000.bf',
 'lm_dcm': '/data/amyloid_brain/LM/17598013_1946_20150604155500.000000.dcm',
 'mumapDCM': '/data/amyloid_brain/umap',
 'nrm_bf': '/data/amyloid_brain/norm/17598013_1946_20150604082431.000000.bf',
 'nrm_dcm': '/data/amyloid_brain/norm/17598013_1946_20150604082431.000000.dcm'}

The parts of the recognised and categorised input data include:

Type	Path
corepath	the core path of the input folder
#mumapDCM	the number of DICOM files for the μ-map, usually 192
mumapDCM	path to the MR-based DICOM μ-map
lm_bf	path to the list-mode binary file
lm_dcm	the path to the DICOM header of the list-mode binary file
nrm_bf	path to the binary file for component based normalisation
nrm_dcm	path to the DICOM header for the normalisation

Note

The raw mMR PET data can also be represented by a single *.ima DICOM file instead of the *.dcm and *.bf pairs for list-mode and normalisation data, which will be reflected in datain.

Specifying output folder

The path to the output folder where the products of NiftyPET go, as well as the verbose mode can be specified as follows:

# output path
opth = os.path.join( datain['corepath'], 'output')

# switch on verbose mode
logging.getLogger().setLevel(logging.DEBUG)

With the setting as above, the output folder output will be created within the input data folder.

Obtaining the hardware and object μ-maps

Since MR cannot image the scanner hardware, i.e., the patient table, head and neck coils, etc., the high resolution CT-based mu-maps are provided by the scanner manufacturer. These then have to be appropriately resampled to the table and coils position as used in any given imaging setting. The hardware and object μ-maps are obtained as follow:

# obtain the hardware mu-map (the bed and the head&neck coil)
muhdct = nipet.hdw_mumap(datain, [1,2,4], mMRpars, outpath=opth, use_stored=True)

# obtain the MR-based human mu-map
muodct = nipet.obj_mumap(datain, mMRpars, outpath=opth, store=True)

The argument [1,2,4] for Obtaining the hardware μ-map correspond to the hardware bits used in imaging, i.e.:

1. Head and neck lower coil

2. Head and neck upper coil

3. Spine coil

4. Table

Currently, the different parts have to be entered manually (they are not automatically recognised which are in use).

The option use_stored=True allows to reuse the already created hardware μ-map, without recalculating it (the resampling can take more than a minute).

Both output dictionaries muhdct and muodct will contain images among other parameters, such as the image affine matrix and image file paths.

In order to check if both μ-maps were properly loaded, the maps can be plotted together transaxially by choosing the axial index iz along the \(z\)-axis, as follows:

# axial index
iz = 60

# plot image with a colour bar
matshow(muhdct['im'][iz,:,:] + muodct['im'][iz,:,:], cmap='bone')
colorbar()

This will produce the following image:

Fig. 7 Composite of the hardware and object μ-maps. Observed can be the human head between the upper and lower head&neck coils, and the patient table below.

The sagittal image can be generated in a similar way, but choosing the slice along the \(x\)-axis, i.e.:

# axial index
ix = 170

# plot image with a colour bar
matshow(muhdct['im'][:,:,ix] + muodct['im'][:,:,ix], cmap='bone')
colorbar()

Fig. 8 Sagittal view of the composite of the hardware and object μ-maps. Observed can be the human head between the upper and lower head&neck coils, and the patient table on the right of the head.

Static image reconstruction

The code below provides full image reconstruction for the last 10 minutes of the acquisition to get an estimate of the amyloid load through the ratio image (SUVr).

recon = nipet.mmrchain(
    datain, mMRpars,
    frames = ['timings', [3000, 3600]],
    mu_h = muhdct,
    mu_o = muodct,
    itr=4,
    fwhm=0.0,
    outpath = opth,
    fcomment = 'niftypet-recon',
    store_img = True)

The input arguments are as follows:

argument	description
datain	input data (list-mode, normalisation and the μ-map)
mMRpars	scanner parameters (scanner constants and LUTs)
frames	definitions of time frame(s);
mu_h	hardware μ-map
mu_o	object μ-map
itr	number of iterations of OSEM (14 subsets).
fwhm	full width at half-maximum for the image post-smoothing
outpath	path to the output folder
fcomment	prefix for all the generated output files
store_img	store images (yes/no)

• the argument timings indicates that the start/stop times in the following sublist is user-specified and can be done for multiple time frames (see section Dynamic time frame definitions).

The reconstructed image can be viewed as follow:

matshow(recon['im'][60,:,:], cmap='magma')
colorbar()

Fig. 9 The transaxial slice of the amyloid PET reconstructed image. Voxel intensities are in Bq.