PET/CT Attenuation Correction

General
1

Hi everyone,
I am trying to write code that creates an attenuation map of the CT dicom and then performs attenuation correction on the PET dicom. I assumed the input CT and PET are already registered.
I also used the bilinear transformation that was mentioned in the this article. but I couldn’t get the desired result. I will be so glad If you give me suggestions to make my code work better.
I am also aware of two drawbacks in my code:

  1. my attenuation map is in 140 Kev, but I need it for 511 Kev. I didn’t make adjustment to accurately reflect the attenuation at 511 keV.
  2. I applied the correction directly on the PET pixel array. I don’t know wether I need a sinogram for my calculations or not.
    I don’t know how can I fix these issues.
    note that the energy of the CT images that I used for AC is 140 kev.

import pydicom
import cv2 as cv
import matplotlib.pyplot as plt
import numpy as np

class AttenuationCorretion:

def __init__(self, dicom_ct, dicom_pet):
    self.dicom_ct = pydicom.dcmread(dicom_ct)
    self.dicom_pet = pydicom.dcmread(dicom_pet)
    self.pixel_pet_uint16 = self.dicom_pet.pixel_array.astype(np.uint16)


def hounsfield_unit(self):
    self.pixel_ct_hu  = self.dicom_ct.pixel_array * self.dicom_ct.RescaleSlope + \
    self.dicom_ct.RescaleIntercept



def att_coe(self):
    self.AC_coef = np.zeros((512,512), dtype=float)
    for i, hu in enumerate(self.pixel_ct_hu.flatten()):
        if hu < 0:
            self.AC_coef.flat[i] = (hu * (0.151 - 0.000225)) / 1000
        elif hu > 0:
            self.AC_coef.flat[i] = 0.151 + (hu * 0.151 * (0.217 - 0.151)) / (1000 * (0.217 - 0.151))
    return self.AC_coef.reshape(pixel_ct_hu.shape)

def down_sample(self):
    self.pixel_ct_am_resized = cv.resize(self.AC_coef, (128,128))
    self.pixel_ct_am_resized_gfilter = cv.GaussianBlur(self.pixel_ct_am_resized, (3,3), 0)

def nor_corr(self):
    self.pixel_ac_pet = self.pixel_pet_uint16 * np.exp(self.pixel_ct_am_resized_gfilter)
    self.normalized_pixel_ac_pet = (self.pixel_ac_pet - self.pixel_ac_pet.min()) / \
        (self.pixel_ac_pet.max() - self.pixel_ac_pet.min())

def visualization(self):
    # create figure 
    fig = plt.figure(figsize=(10, 7)) 

    # setting values to rows and column variables 
    rows = 2
    columns = 2

    # Adds a subplot at the 1st position 
    fig.add_subplot(rows, columns, 1) 

    # showing image 

    plt.imshow(self.dicom_ct.pixel_array, cmap = 'gray') 
    plt.axis('off') 
    plt.title("CT") 

    # Adds a subplot at the 2nd position 
    fig.add_subplot(rows, columns, 2) 

    # showing image 


    plt.imshow(self.pixel_pet_uint16, cmap = 'gray') 
    plt.axis('off') 
    plt.title("nonAC PET") 

    # Adds a subplot at the 3rd position 
    fig.add_subplot(rows, columns, 3) 

    # showing image 

    
    plt.imshow(self.pixel_ct_am_resized_gfilter, cmap = 'gray') 
    plt.axis('off') 
    plt.title("Attenuation Map") 

    # Adds a subplot at the 4th position 
    fig.add_subplot(rows, columns, 4) 


    

    # showing image 
    plt.imshow(self.normalized_pixel_ac_pet, cmap = 'gray') 
    plt.axis('off') 
    plt.title("AC PET") 

    return fig

def process(self):
    self.hounsfield_unit()
    self.att_coe()
    self.down_sample()
    self.nor_corr()
    self.visualization()

attenuation_map = AttenuationCorretion(‘ctdcm/ct_file.dcm’,
‘petdcm/pet_file.dcm’)

attenuation_map.process()

AC
AC687×581 49.2 KB