Uncertainty Estimation

Special thanks to Dr. Peyman Sheikhzadeh at the Tehran University Of Medical Science for providing this anonymized data

import os
import numpy as np
from pytomography.io.SPECT import dicom
from pytomography.transforms.SPECT import SPECTAttenuationTransform, SPECTPSFTransform
from pytomography.algorithms import OSEM, PGAAMultiBedSPECT
from pytomography.projectors.SPECT import SPECTSystemMatrix
from pytomography.likelihoods import PoissonLogLikelihood
from pytomography.utils import print_collimator_parameters
import matplotlib.pyplot as plt
import torch
from rt_utils import RTStructBuilder
from pytomography.callbacks import DataStorageCallback

save_path = '/mnt/mydisk2/pytomo_tutorial_data/SPECT'

The following code cell loads projection data and a corresponding CT file and creates a reconstruction algorithm; if you are unfamiliar with the code below or require more explanation, see that “DICOM Data Introduction” tutorial.

• The only difference is that we now get a variance estimate for the scatter in line 7 of the cell below, and then input that to the likelihood in line 22

• The collimator corresponding to this data is General Electric s medium energy general purpose collimator.

file_NM = os.path.join(save_path, 'Lu177-PSMA-GEDisc', 'bed2_projections.dcm')
path_CT = os.path.join(save_path, 'Lu177-PSMA-GEDisc', 'CT')
files_CT = [os.path.join(path_CT, file) for file in os.listdir(path_CT)]
object_meta, proj_meta = dicom.get_metadata(file_NM, index_peak=1)
projections = dicom.get_projections(file_NM)
photopeak = projections[1]
scatter, scatter_variance_estimate = dicom.get_energy_window_scatter_estimate(file_NM, index_peak=1, index_lower=3, index_upper=2, return_scatter_variance_estimate=True)
# Build system matrix
attenuation_map = dicom.get_attenuation_map_from_CT_slices(files_CT, file_NM, index_peak=1)
psf_meta = dicom.get_psfmeta_from_scanner_params('GI-MEGP', energy_keV=208)
att_transform = SPECTAttenuationTransform(attenuation_map)
psf_transform = SPECTPSFTransform(psf_meta)
system_matrix = SPECTSystemMatrix(
    obj2obj_transforms = [att_transform,psf_transform],
    proj2proj_transforms = [],
    object_meta = object_meta,
    proj_meta = proj_meta)
# Likelihood
likelihood = PoissonLogLikelihood(system_matrix, photopeak, additive_term=scatter, additive_term_variance_estimate=scatter_variance_estimate)
# Reconstruction algorithm
recon_algorithm = OSEM(likelihood)

Given photopeak energy 208.0 keV and CT energy 120 keV from the CT DICOM header, the HU->mu conversion from the following configuration is used: 208.0 keV SPECT energy, 120 keV CT energy, and scanner model optimact540

Up to this point, we’ve created the reconstruction algorithm but we haven’t performed reconstruction yet. If we plan on computing uncertainty in various VOIs, we are required to use the DataStorageCallback callback in the image reconstruction algorithm. This callback stores a copy of the reconstructed image at each iteration and subiteration: all this data is required when computing uncertainties later. It can be created as follows:

data_storage_callback = DataStorageCallback(likelihood, torch.clone(recon_algorithm.object_prediction))

Then we reconstruct using this callback:

recon_OSEM = recon_algorithm(n_iters = 4, n_subsets = 8, callback=data_storage_callback)

Here’s a coronal maximum intensity projection of the reconstruction

maximum_intensity_projection = recon_OSEM.max(axis=1)[0].cpu()

plt.figure(figsize=(4,4))
plt.imshow(maximum_intensity_projection.T, cmap='Greys', vmax=4, interpolation='gaussian', origin='lower')
plt.axis('off')
plt.show()

Suppose we want to estimate the uncertainty in activity in one of the kidneys. Let’s obtain it from an RTStruct file containing segmentations of all the organs. (The cell below is not required for this, but it will allow us to determine the mask names contained in the RTStruct file):

file_RT = os.path.join(save_path, 'Lu177-PSMA-GEDisc', 'segmentations.dcm')
rtstruct = RTStructBuilder.create_from(
        dicom_series_path=path_CT,
        rt_struct_path=file_RT
    )
print(rtstruct.get_roi_names())

['spleen', 'kidney_right', 'kidney_left', 'gallbladder', 'liver', 'stomach', 'pancreas', 'adrenal_gland_right', 'adrenal_gland_left', 'lung_upper_lobe_left', 'lung_lower_lobe_left', 'lung_upper_lobe_right', 'lung_middle_lobe_right', 'lung_lower_lobe_right', 'esophagus', 'trachea', 'thyroid_gland', 'small_bowel', 'duodenum', 'colon', 'urinary_bladder', 'sacrum', 'vertebrae_S1', 'vertebrae_L5', 'vertebrae_L4', 'vertebrae_L3', 'vertebrae_L2', 'vertebrae_L1', 'vertebrae_T12', 'vertebrae_T11', 'vertebrae_T10', 'vertebrae_T9', 'vertebrae_T8', 'vertebrae_T7', 'vertebrae_T6', 'vertebrae_T5', 'vertebrae_T4', 'vertebrae_T3', 'vertebrae_T2', 'vertebrae_T1', 'vertebrae_C7', 'vertebrae_C6', 'vertebrae_C5', 'vertebrae_C4', 'vertebrae_C3', 'vertebrae_C2', 'vertebrae_C1', 'heart', 'aorta', 'pulmonary_vein', 'brachiocephalic_trunk', 'subclavian_artery_right', 'subclavian_artery_left', 'common_carotid_artery_right', 'common_carotid_artery_left', 'brachiocephalic_vein_left', 'brachiocephalic_vein_right', 'atrial_appendage_left', 'superior_vena_cava', 'inferior_vena_cava', 'portal_vein_and_splenic_vein', 'iliac_artery_left', 'iliac_artery_right', 'iliac_vena_left', 'iliac_vena_right', 'humerus_left', 'humerus_right', 'scapula_left', 'scapula_right', 'clavicula_left', 'clavicula_right', 'femur_right', 'hip_left', 'hip_right', 'spinal_cord', 'gluteus_maximus_left', 'gluteus_maximus_right', 'gluteus_medius_left', 'gluteus_medius_right', 'gluteus_minimus_left', 'gluteus_minimus_right', 'autochthon_left', 'autochthon_right', 'iliopsoas_left', 'iliopsoas_right', 'brain', 'skull', 'rib_right_4', 'rib_right_3', 'rib_left_1', 'rib_left_2', 'rib_left_3', 'rib_left_4', 'rib_left_5', 'rib_left_6', 'rib_left_7', 'rib_left_8', 'rib_left_9', 'rib_left_10', 'rib_left_11', 'rib_left_12', 'rib_right_1', 'rib_right_2', 'rib_right_5', 'rib_right_6', 'rib_right_7', 'rib_right_8', 'rib_right_9', 'rib_right_10', 'rib_right_11', 'rib_right_12', 'sternum', 'costal_cartilages']

We can see the 'kidney_right' is listed as one of the available masks. We can obtain this mask as follows:

• file_RT is the path to the RTStruct file

• file_NM is the path the the SPECT projection data; this is required for aligning the mask

• dicom_series_path is the reference DICOM path of the RTStruct

• rt_struct_name is the name of the RTStruct we wish to open

mask_name = 'kidney_left'
kidney_mask = dicom.get_aligned_rtstruct(
    file_RT = file_RT,
    file_NM = file_NM,
    dicom_series_path = path_CT,
    rt_struct_name = mask_name
)

We can compute the relative and absolute uncertainty corresponding to counts in this region using the compute_uncertainty function of the reconstruction algorithm. We also need to pass in the DataStorageCallback which contains every iterative update to the objects.

uncertainty_abs, uncertainty_pct = recon_algorithm.compute_uncertainty(
    mask = kidney_mask,
    data_storage_callback = data_storage_callback,
    return_pct = True,
    include_additive_term=True # must have additive_term_variance_estimate in likelihood
)
print(f'Estimated uncertainty in {mask_name}: {uncertainty_pct:.2f}%')

Estimated uncertainty in kidney_left: 0.94%

This can be repeated for any of the masks above! In addition, you can play around the with iterations, subsets (or even try different reconstruction algorithms) and see how it effects the estimated uncertainty.

Multiple Bed Positions

PyTomography is also capable of uncertainty estimation in masks that are contained in two seperate bed positions. This uncertainty estimation takes into account the stitching weights between two different bed positions (if the mask is near the edge) when computing the organ variance in each seperate region.

• The code below is a summary of the multibed position tutorial, except that we now include a DataStorageCallback for each bed position

files_NM = [
    os.path.join(save_path, 'Lu177-PSMA-GEDisc', 'bed1_projections.dcm'),
    os.path.join(save_path, 'Lu177-PSMA-GEDisc', 'bed2_projections.dcm'),
]
path_CT = os.path.join(save_path, 'Lu177-PSMA-GEDisc', 'CT')
files_CT = [os.path.join(path_CT, file) for file in os.listdir(path_CT)]

projectionss = dicom.load_multibed_projections(files_NM)


def initialize_reconstruction_algorithm_singlebed(i):
    # Change these depending on your file:
    index_peak = 1
    index_lower = 3
    index_upper = 2
    projections = projectionss[i]
    file_NM = files_NM[i]
    object_meta, proj_meta = dicom.get_metadata(file_NM, index_peak=index_peak)
    photopeak = projections[index_peak]
    scatter = dicom.get_energy_window_scatter_estimate_projections(file_NM, projections, index_peak, index_lower, index_upper)
    # Build system matrix
    attenuation_map = dicom.get_attenuation_map_from_CT_slices(files_CT, file_NM, index_peak=index_peak)
    psf_meta = dicom.get_psfmeta_from_scanner_params('GI-MEGP', energy_keV=208)
    att_transform = SPECTAttenuationTransform(attenuation_map)
    psf_transform = SPECTPSFTransform(psf_meta)
    # Create system matrix
    system_matrix = SPECTSystemMatrix(
        obj2obj_transforms = [att_transform,psf_transform],
        proj2proj_transforms= [],
        object_meta = object_meta,
        proj_meta = proj_meta)
    likelihood = PoissonLogLikelihood(system_matrix, photopeak, additive_term=scatter)
    reconstruction_algorithm = OSEM(likelihood)
    # Return only the reconstruction algorithm initialization (not calling it for any iterations/subsets yet)
    return reconstruction_algorithm

In order to compute uncertainty, we need to use the PGAAMultiBedWrapper reconstruction algorithm, which takes in multiple initialized reconstruction algorithms for each bed position (this class can also be used to reconstruct multiple bed positions without uncertainty as well)

recon_algo_upper = initialize_reconstruction_algorithm_singlebed(0)
recon_algo_lower  = initialize_reconstruction_algorithm_singlebed(1)
recon_algo = PGAAMultiBedSPECT(files_NM, [recon_algo_upper, recon_algo_lower])
# Initialize callback for each bed position
callbacks = [DataStorageCallback(r.likelihood, r.object_prediction) for r in recon_algo.reconstruction_algorithms]
reconstructed_image_multibed = recon_algo(4, 8, callback=callbacks)

Given photopeak energy 208.0 keV and CT energy 120 keV from the CT DICOM header, the HU->mu conversion from the following configuration is used: 208.0 keV SPECT energy, 120 keV CT energy, and scanner model optimact540
Given photopeak energy 208.0 keV and CT energy 120 keV from the CT DICOM header, the HU->mu conversion from the following configuration is used: 208.0 keV SPECT energy, 120 keV CT energy, and scanner model optimact540

Lets plot the reconstruction (note how the liver overlaps two different bed positions, since it spans the center of the image)

maximum_intensity_projection = reconstructed_image_multibed.max(axis=1)[0].cpu()
plt.figure(figsize=(3,5))
plt.imshow(maximum_intensity_projection.T, cmap='Greys', vmax=5, origin='lower', interpolation='gaussian')
plt.axis('off')
plt.show()

When we obtain the organ mask, we need to

1. Specify the shape as the shape of the reconstructed image (not the batch dimension, so we use [1:])

2. Specify the file_NM in the most superior bed position

mask_name = 'liver'
liver_mask = dicom.get_aligned_rtstruct(
    file_RT = file_RT,
    file_NM = files_NM[0],
    dicom_series_path = path_CT,
    rt_struct_name = mask_name,
    shape=reconstructed_image_multibed.shape
)

Now we can compute the uncertainty in the liver, which spans both bed positions

uncertainty_abs, uncertainty_pct = recon_algo.compute_uncertainty(liver_mask, callbacks, return_pct=True)

print(f'Estimated uncertainty in {mask_name}: {uncertainty_pct:.2f}%')

Estimated uncertainty in liver: 1.56%