OSEM

OSEM (Ordered Subsets Expectation Maximization) is the main reconstruction algorithm.

Basic Usage

This is an example usage performing the OSEM reconstruction of a list-mode dataset.

import pyyrtpet as yrt

scanner = yrt.Scanner("<MyScannerFile.json>")

lm = yrt.ListModeLUTOwned(scanner, "<MyListMode.lmDat>")

# Create OSEM with scanner
osem = yrt.OSEM(scanner)

# Set reconstruction parameters
osem.num_MLEM_iterations = 10
osem.num_OSEM_subsets = 8

# Set image parameters
img_params = yrt.ImageParams(nx=128, ny=128, nz=64,
    length_x=128.0, length_y=128.0, length_z=64.0)
osem.setImageParams(img_params)

# Set data input
osem.setDataInput(lm)

# Set projector
osem.setProjector("Siddon")

# Generate sensitivity images (Will also save the sensitivity image into a NIfTI file)
[sens_image] = osem.generateSensitivityImages("sens_image.nii.gz")

# Set sensitivity image(s)
osem.setSensitivityImage(sens_image)

# Run reconstruction (Will also save the reconstructed image into a NIfTI file)
result = osem.reconstruct("recon_image.nii.gz")

Configuration

Reconstruction Parameters

#<>
osem.num_MLEM_iterations = 20  # Number of iterations
osem.num_OSEM_subsets = 8      # Number of subsets

Data Input

#<>
osem.setDataInput(projection_data)

Projector

#<>
osem.setProjector("DD")      # For the distance-driven projector
osem.setProjector("SIDDON")  # For the Siddon projector

osem.setNumRays(3)                 # Multi-ray (Siddon only)
osem.addProjPSF("<proj_psf.csv>")  # Adding projection-space PSF (DD only)

# Time-of-flight (TOF width in picoseconds, number of standard deviations to fill
#    from the TOF kernel)
osem.addTOF(500.0, 3)

Corrections

Sensitivity correction (sometimes also called normalization correction) can be done by providing a histogram (which can be of any format that inherits from Histogram, not necessarily Histogram3D) of the sensitivity coefficient for every detector pair.

Attenuation correction can be done by providing an attenuation image or by providing a histogram of Attenuation Correction Factors (ACFs).

Randoms correction can be performed by either providing a histogram of the randoms estimate for every detector pair, or by encoding the randoms estimate for every event (overloaded in the getRandomsEstimate member function)

Scatter correction (Not to be confused with scatter estimation) is perfomed by providing a histogram of the scatter estimate for every detector pair.

#<>
osem.setSensitivityHistogram(sensitivity_his)
osem.setAttenuationImage(att_image)   # Attenuation correction (From mu-map)
osem.setACFHistogram(acf_his)         # Attenuation correction (From histogram)
osem.setRandomsHistogram(randoms_his) # Randoms correction (From histogram)
osem.setScatterHistogram(scatter_his) # Scatter correction

Point Spread Function (Image-space)

#<>
osem.addImagePSF("<psf.csv>")                           # Uniform PSF
osem.addImagePSF("<psf_variant.csv>", yrt.ImagePSFMode.VARIANT)  # Spatially variant PSF

PSF Modes

yrt.ImagePSFMode.UNIFORM - Same PSF for all voxels
yrt.ImagePSFMode.VARIANT - Spatially variant PSF

Output Options

To save intermediate iteration, you must provide a RangeList object

# Save intermediate iterations
range_list = yrt.RangeList()
range_list.insertSorted(0,1) # Save iterations 0 and 1
range_list.insertSorted(5,8) # Save iterations 5,6,7,8
range_list.insertSorted(10,10) # Save iteration 10

to the OSEM object

#<>
# The filename given will be added a prefix to specify which iteration is was saved
osem.setSaveIterRanges(range_list, "./iterations/recon_image_intermediate.nii.gz")

# Custom initial estimate (instead of a uniform image containing the value 0.1)
osem.setInitialEstimate(initial_image)

# Reconstruction mask (To disable some voxels from the reconstruction)
osem.setMaskImage(mask_image)

Output

The reconstruct() method returns an ImageOwned object, but can also save the reconstructed image to disk if provided a filename.

#<>
result = osem.reconstruct() # Will not save to disk
result = osem.reconstruct("recon_image.nii.gz") # Will save to disk

# Access as numpy array
result_np = np.array(result, copy=False)