Hello everyone,
I was wondering if I could receive some guidance in terms of CT attenuation correction.
I went through the reconstruction tutorials without any issues but when I wanted to try it myself with some Lu-177 phantom data from a SPECT/CT system (Siemens, Symbia Intevo) I ran into issues.
The issue occurs when I try use SPECT attenuation transformation, and I get this error message
Warning: Could not find cortical bone peak: defaulting to 75% kVp value
warnings.warn(“Could not find cortical bone peak: defaulting to 75% kVp value”, category=Warning)
Meaning that I can’t obtain the effective energy from the CT images. I do see that my kVp is 130 in the DICOM header. Should I choose 2/3 of the kVp for my average energy as the effective energy?
Is there a way can I force PyTomography to use my 130 kVp and convert to effective energy, to convert to attenuation coefficient so I can have a mu map?
I tried to use the HU_to_mu found in pytomography.io.CT.attenuation_map , with the arguments files_CT, kvp=130 and energy of SPECT=208, but I still encounter the same error message.
It looks like the effective energy can be obtained by plotting a histogram over the HU values. Any ideas on how I’ll be able to get an effective energy from that? Won’t the histogram only show the number of occurrence on the vertical axis and HU value on the horizontal. Suggestions on how I can interpret the histogram into energy?
1 Like
I made an attempt to create my attenuation map.
Any comments on if this is the correct way to do it or not?
How can we confirm that the attenuation map is correct?
I was able to obtain new images after saving, however it looks like I only got the HU converted images.
Any suggestions on how I should correct the script?
Might have the wrong attenuation coefficient?
Attached are the histograms for the supposedly HU and attenuation coefficient converted histograms. The histogram looks like they match each other beside less frequent.
import os
import numpy as np
import pydicom
import matplotlib.pyplot as plt
#%%
ct_slices_dir = “Directory location for CT files”
#output_dir = “Directory location for the results, attenuation map/”
#%%
HU_air = -1024 # HU for air when converting the pixel values to HU in CT images
HU_water = 0 # HU for water when converting the pixel values to HU in CT images
attenuation_air = 0.001 # Attenuation coefficient for air at 130 keV?
attenuation_water = 0.2 # Attenuation coefficient for water at 130 keV?
slope_attenuation_coeff=(attenuation_water - attenuation_air) / (HU_water - HU_air)
#%%
attenuation_slices =
for filename in os.listdir(ct_slices_dir):
ds = pydicom.dcmread(os.path.join(ct_slices_dir, filename))
%Converting CT files to HU
HU_slices=((ds.pixel_array * ds.RescaleSlope) + ds.RescaleIntercept)
%Converting image pixels to attenuation coefficient
attenuation_slices=(attenuation_air + slope_attenuation_coeff * (HU_slices - HU_air))
Rewriting dicom file so it can be saved and visualized .
ds.PixelArray=attenuation_slices.astype(np.uint16).tobytes()
ds.save_as(‘Output name’)
Any ides/suggestions or comments are welcomed.
Hi @jakopbsf
I believe you have the right idea in creating the attenuation map, but the issue is that the relationship tends to change for high HU values. i.e. for -1024HU->0HU the relationship is characterized by one linear curve, while for 0HU->High HU, the relationship tends to be characterized by a different curve. All the stuff with effective kVp, etc, only applies to the HU above water region in PyTomography. Basically, there’s no reference point for a high HU region (unlike air/water), so we attempt to use cortical bone as a reference point within patients. If your phantom doesn’t have any high HU regions it shouldn’t really make any difference.
Also, for air/water, I think you want the attenuation coefficients at your isotope energy, not your CT energy.
Also, the issue is that the CT photon energy isnt constant (130kVp specifies the upper bound) and hence we have to take a sort-of average value.
1 Like