compute suv from pet

README.md

** PLEASE NOTE ***
there are many other solutions for computing SUV, see below link for more info.
https://qibawiki.rsna.org/index.php/Standardized_Uptake_Value_(SUV)

This gist contains code to perform the below with python (3.7).

+ download data from TCIA via tcia-rest api. (download.py, tciaclient.py)

+ compute SUV from PET. (utils.py)

See suv-demo.ipynb for example usage.

download.py

import sys,os
from tciaclient import TCIAClient

if __name__ == '__main__':
    
    series_instance_uid = sys.argv[1]
    file_path = sys.argv[2]

    folder = os.path.dirname(file_path)
    basename = os.path.basename(file_path)
    if not basename.endswith('.zip'):
        basename = basename+'.zip'

    tcia_client = TCIAClient(apiKey=None, baseUrl="https://services.cancerimagingarchive.net/services/v3",resource="TCIA")
    tcia_client.get_image(seriesInstanceUid=series_instance_uid,downloadPath=folder,zipFileName=basename)

'''
!rm -rf img
!mkdir -p img/pet;mkdir -p img/ct;
!python single_download.py 1.3.6.1.4.1.14519.5.2.1.5099.8010.427264300850965737262860055580 img/ct/ct.zip
!python single_download.py 1.3.6.1.4.1.14519.5.2.1.5099.8010.308184765901558710285007064772 img/pet/pet.zip
!unzip img/ct/ct.zip -d img/ct;unzip img/pet/pet.zip -d img/pet;

SeriesInstanceUID,StudyInstanceUID,Modality,ProtocolName,SeriesDate,SeriesDescription,BodyPartExamined,SeriesNumber,Collection,Manufacturer,ManufacturerModelName,SoftwareVersions,Visibility,ImageCount,PatientID,PatientName,PatientSex,StudyDate,StudyDescription,PatientAge,SeriesCount
1.3.6.1.4.1.14519.5.2.1.5099.8010.427264300850965737262860055580,1.3.6.1.4.1.14519.5.2.1.5099.8010.199920086920823171706454903251,CT,HeadNeckPETCT,1999-08-23,CT 5.0 H30s,HEADNECK,2.000000,Head-Neck Cetuximab,SIEMENS,Emotion,VA40C,1,113,0522c0001,0522c0001,F,1999-08-23,Neck^HeadNeckPETCT,000D,3
1.3.6.1.4.1.14519.5.2.1.5099.8010.308184765901558710285007064772,1.3.6.1.4.1.14519.5.2.1.5099.8010.199920086920823171706454903251,PT,,1999-08-23,PET WB,HEADNECK,606.000000,Head-Neck Cetuximab,CPS,1062,,1,112,0522c0001,0522c0001,F,1999-08-23,Neck^HeadNeckPETCT,000D,3
'''

suv-demo.ipynb

tciaclient.py

# https://raw.githubusercontent.com/TCIA-Community/TCIA-API-SDK/master/tcia-rest-client-python/src/tciaclient.py
import os
import urllib3
import urllib
import sys
import math
import traceback
#
# Refer https://wiki.cancerimagingarchive.net/display/Public/REST+API+Usage+Guide for complete list of API
#
class TCIAClient:
    GET_IMAGE = "getImage"
    GET_MANUFACTURER_VALUES = "getManufacturerValues"
    GET_MODALITY_VALUES = "getModalityValues"
    GET_COLLECTION_VALUES = "getCollectionValues"
    GET_BODY_PART_VALUES = "getBodyPartValues"
    GET_PATIENT_STUDY = "getPatientStudy"
    GET_SERIES = "getSeries"
    GET_PATIENT = "getPatient"
    GET_SERIES_SIZE = "getSeriesSize"
    CONTENTS_BY_NAME = "ContentsByName"

    def __init__(self, apiKey, baseUrl, resource, maxsize=5):
        self.apiKey = apiKey
        self.baseUrl = baseUrl + "/" + resource
        self.pool_manager = urllib3.PoolManager(maxsize=maxsize)
    def execute(self, url, queryParameters={},preload_content=True):
        queryParameters = dict((k, v) for k, v in queryParameters.items() if v)
        
        headers = {}
        if self.apiKey is not None:
            headers = {"api_key" : self.apiKey }

        queryString = "?%s" % urllib.parse.urlencode(queryParameters)
        requestUrl = url + queryString        
        request = self.pool_manager.request(method='GET', url=requestUrl , headers=headers, preload_content=preload_content)
        return request

    def get_modality_values(self,collection = None , bodyPartExamined = None , modality = None , outputFormat = "json" ):
        serviceUrl = self.baseUrl + "/query/" + self.GET_MODALITY_VALUES
        queryParameters = {"Collection" : collection , "BodyPartExamined" : bodyPartExamined , "Modality" : modality , "format" : outputFormat }
        resp = self.execute(serviceUrl , queryParameters)
        return resp

    def get_series_size(self, SeriesInstanceUID = None, outputFormat = "json"):
        serviceUrl = self.baseUrl + "/query/" + self.GET_SERIES_SIZE
        queryParameters = {"SeriesInstanceUID" : SeriesInstanceUID, "format" :
                           outputFormat}
        resp = self.execute(serviceUrl, queryParameters)
        return resp

    def contents_by_name(self, name = None):
        serviceUrl = self.baseUrl + "/query/" + self.CONTENTS_BY_NAME
        queryParameters = {"name" : name}
        resp = self.execute(serviceUrl,queryParameters)
        return resp

    def get_manufacturer_values(self,collection = None , bodyPartExamined = None, modality = None , outputFormat = "json"):
        serviceUrl = self.baseUrl + "/query/" + self.GET_MANUFACTURER_VALUES
        queryParameters = {"Collection" : collection , "BodyPartExamined" : bodyPartExamined , "Modality" : modality , "format" : outputFormat }
        resp = self.execute(serviceUrl , queryParameters)
        return resp

    def get_collection_values(self,outputFormat = "json" ):
        serviceUrl = self.baseUrl + "/query/" + self.GET_COLLECTION_VALUES
        queryParameters = { "format" : outputFormat }
        resp = self.execute(serviceUrl , queryParameters)
        return resp

    def get_body_part_values(self,collection = None , bodyPartExamined = None , modality = None , outputFormat = "json" ):
        serviceUrl = self.baseUrl + "/query/" + self.GET_BODY_PART_VALUES
        queryParameters = {"Collection" : collection , "BodyPartExamined" : bodyPartExamined , "Modality" : modality , "format" : outputFormat }
        resp = self.execute(serviceUrl , queryParameters)
        return resp

    def get_patient_study(self,collection = None , patientId = None , studyInstanceUid = None , outputFormat = "json" ):
        serviceUrl = self.baseUrl + "/query/" + self.GET_PATIENT_STUDY
        queryParameters = {"Collection" : collection , "PatientID" : patientId , "StudyInstanceUID" : studyInstanceUid , "format" : outputFormat }
        resp = self.execute(serviceUrl , queryParameters)
        return resp

    def get_series(self, collection = None , modality = None , studyInstanceUID = None, seriesInstanceUID = None , outputFormat = "json" ):
        serviceUrl = self.baseUrl + "/query/" + self.GET_SERIES
        queryParameters = {"Collection" : collection , "StudyInstanceUID": studyInstanceUID, "SeriesInstanceUID" : seriesInstanceUID  , "Modality" : modality , "format" : outputFormat }
        resp = self.execute(serviceUrl , queryParameters)
        return resp

    def get_patient(self,collection = None , outputFormat = "json" ):
        serviceUrl = self.baseUrl + "/query/" + self.GET_PATIENT
        queryParameters = {"Collection" : collection , "format" : outputFormat }
        resp = self.execute(serviceUrl , queryParameters)
        return resp

    def get_image(self , seriesInstanceUid , downloadPath, zipFileName):
        serviceUrl = self.baseUrl + "/query/" + self.GET_IMAGE
        queryParameters = { "SeriesInstanceUID" : seriesInstanceUid }
        #os.umask(0002)
        try:
            file = os.path.join(downloadPath, zipFileName)
            resp = self.execute( serviceUrl , queryParameters, preload_content=False)
            downloaded = 0
            CHUNK = 256 * 10240
            with open(file, 'wb') as fp:
                for chunk in resp.stream(CHUNK):
                    fp.write(chunk)
        except urllib3.exceptions.HTTPError as e:
            print("HTTP Error:",e.code , serviceUrl)
            return False
        except:
            traceback.print_exc()
            return False

        return True

utils.py

import os
import traceback
import datetime
import numpy as np
import pydicom
import SimpleITK as sitk

def sort_by_instance_number(image_file_list):
    data = []
    for row in image_file_list:        
        f=pydicom.dcmread(row)
        data.append({'f':row,'n':f.InstanceNumber})
    data=sorted(data,key=lambda x: x['n'])
    return [x['f'] for x in data]

def imread(fpath):
    if isinstance(fpath,list):
        image_file_list = fpath
        image_file_list = sort_by_instance_number(image_file_list)
        reader = sitk.ImageSeriesReader()
        reader.SetFileNames(image_file_list)
    elif fpath.endswith('.list'):
        with open(fpath,'r') as f:
            dicom_names = [x for x in f.read().split('\n') if len(x) > 0]
        if not os.path.exists(dicom_names[0]):
            image_file_list = [os.path.join(os.path.dirname(fpath),x) for x in dicom_names]
            image_file_list = sort_by_instance_number(image_file_list)
        reader = sitk.ImageSeriesReader()
        reader.SetFileNames(image_file_list)
    else:
        reader= sitk.ImageFileReader()
        reader.SetFileName(fpath)
    img = reader.Execute()
    arr = sitk.GetArrayFromImage(img)    
    spacing = img.GetSpacing()
    origin = img.GetOrigin()
    direction = img.GetDirection()    
    return arr,spacing,origin,direction

'''

One method to compute suv.
translated from https://github.com/mvallieres/radiomics by Martin Vallières
more refs.
https://wiki.cancerimagingarchive.net/display/Public/Head-Neck-PET-CT#0b323d6250cc42fa8fa09821fabe0bd7    
https://qibawiki.rsna.org/index.php/Standardized_Uptake_Value_(SUV)

'''
def compute_suv(image_file_list):
    
    estimated = False

    raw,spacing,origin,direction = imread(image_file_list)
    
    f=pydicom.dcmread(image_file_list[0])
    
    try:
        weight_grams = float(f.PatientWeight)*1000
    except:
        traceback.print_exc()
        weight_grams = 75000
        estimated = True
        
    try:
        # Get Scan time
        scantime = datetime.datetime.strptime(f.AcquisitionTime,'%H%M%S.%f')
        # Start Time for the Radiopharmaceutical Injection
        injection_time = datetime.datetime.strptime(f.RadiopharmaceuticalInformationSequence[0].RadiopharmaceuticalStartTime,'%H%M%S.%f')
        # Half Life for Radionuclide # seconds
        half_life = float(f.RadiopharmaceuticalInformationSequence[0].RadionuclideHalfLife) 
        # Total dose injected for Radionuclide
        injected_dose = float(f.RadiopharmaceuticalInformationSequence[0].RadionuclideTotalDose)

        # Calculate decay
        decay = np.exp(-np.log(2)*((scantime-injection_time).seconds)/half_life);
        # Calculate the dose decayed during procedure
        injected_dose_decay = injected_dose*decay; # in Bq        
    except:
        traceback.print_exc()
        decay = np.exp(-np.log(2)*(1.75*3600)/6588); # 90 min waiting time, 15 min preparation
        injected_dose_decay = 420000000 * decay; # 420 MBq
        estimated = True
    
    # Calculate SUV # g/ml
    suv = raw*weight_grams/injected_dose_decay
    
    return suv, estimated, raw,spacing,origin,direction

updated version of the similar code for ONLY image download is available below (commited as a code repo).
https://github.com/pangyuteng/tcia-image-download-python