1. Calculate count in given ROI. 2. Calculate parameter of certain ROI (chosen by threshold)

reconstructin_ROI.py

import dicom
from skimage.transform import iradon, iradon_sart
import numpy
import pylab

# readind data
ds = dicom.read_file("img.dcm") # read a dicom of pre-reconstructed SPECT
pixel_size= ds.PixelSpacing[1]/10 # read pixel size from header of dicom
print("Pixel size =", pixel_size, "cm")
pixel_volume= pixel_size**3
print("Pixel volume =", pixel_volume, "cm^3")

array = ds.pixel_array # read the array

# reconstruction
recon_3d = numpy.zeros((64,64,64)) # create a blank 3D imaging

for i in range(0,63):
    recon_3d[i,:,:] = iradon(numpy.transpose(array[:,i,:]), output_size=64) # create axial imaging

for i in range(28,39): #coronary view
    pylab.imshow(recon_3d[:,i,:])
    pylab.show()

ROI = numpy.zeros((64,64,64)) # create a blank ROI array

ROI[29:43, 28:39, 25:39] = 1 # assign certain area as selected
select = recon_3d * ROI

for i in range(28,39):
    pylab.imshow(select[:,i,:])
    pylab.show()

print("total count of selected =", numpy.sum(select[:,:,:]))
print("total volume of selected =", numpy.sum(ROI[:,:,:])*pixel_volume, "cm^3")
print("mean count per cm^3 =", numpy.sum(select[:,:,:])/(numpy.sum(ROI[:,:,:])*pixel_volume))


# threshold
maximal = numpy.max(select)

for i in range(0,63):
    for j in range(0,63):
        for k in range(0,63):
            if select[i,j,k] < 0.5*maximal:  # only select more than half of maximal
                ROI[i,j,k]= 0

select = recon_3d * ROI
                
for i in range(28,39):
    pylab.imshow(select[:,i,:])
    pylab.show()

print("Total volume of myocardium =", numpy.sum(ROI[:,:,:])*pixel_volume, "cm^3")
print("Total heart uptake =", numpy.sum(select[:,:,:]))