"""scan module.
This module defines pyCERR data object for images (CT, MR, PT, US).
Metadata can be imported from various file formats such as DICOM, NifTi.
It also provides methods to transform the Scan object to other formats such NifTi, SimpleITK
and for converting images to real world units and SUV calculation.
"""
from dataclasses import dataclass, field
import numpy as np
from pydicom import dcmread
from pydicom.dataset import Dataset, FileDataset, FileMetaDataset
import cerr.dataclasses.scan_info as scn_info
import SimpleITK as sitk
import json
import os
[docs]
def get_empty_list():
"""Return an empty list, used as a default factory for dataclass fields.
Returns:
list: An empty list.
"""
return []
[docs]
def get_empty_np_array():
"""Return an empty 3-D numpy array, used as a default factory for dataclass fields.
Returns:
np.ndarray: A zero-element array with shape (0, 0, 0).
"""
return np.empty((0,0,0))
[docs]
@dataclass
class Scan:
"""This class defines data object for volumetric images such as CT, MR, PET or derived image type.
Attributes:
scanArray (np.ndarray): numpy array for the image.
scanType (str): Type of scan. e.g. 'CT SCAN'
scanInfo (cerr.dataclasses.scan_info.ScanInfo): scan_info object containing metadata for each scan slice
scanUID (str): unique identifier for each scan.
assocDeformUID (str): optional, UID of associated deformation object that was used to generate this scan.
assocTextureUID (str): optional, UID of associated texture object that was used to generate this scan.
assocBaseScanUID (str): optional, UID of associated base scan in the deformation that was used to generate this scan.
assocMovingScanUID (str): optional, UID of associated moving scan in the deformation that was used to generate this scan.
Image2PhysicalTransM (np.ndarray): Transformation matrix to convert pyCERR row,col,slc to DICOM physical coordinates.
Image2VirtualPhysicalTransM (np.ndarray): Transformation matrix to convert pyCERR's scan row,col,slc to pyCERR virtual coordinates.
cerrToDcmTransM (np.ndarray): Transformation matrix to convert pyCERR virtual x,y,z coordinates to DICOM physical coordinates.
"""
scanArray: np.ndarray = field(default_factory=get_empty_np_array)
scanType: str = ''
scanInfo: scn_info.ScanInfo = field(default_factory=list)
uniformScanInfo: scn_info.UniformScanInfo = field(default_factory=get_empty_list)
scanArraySuperior: np.ndarray = field(default_factory=get_empty_np_array)
scanArrayInferior: np.ndarray = field(default_factory=get_empty_np_array)
thumbnails: np.ndarray = field(default_factory=get_empty_np_array)
transM: np.ndarray = field(default_factory=get_empty_np_array)
scanUID: str = ''
assocDeformUID: str = ''
assocTextureUID: str = ''
assocBaseScanUID: str = ''
assocMovingScanUID: str = ''
Image2PhysicalTransM: np.ndarray = field(default_factory=get_empty_np_array)
Image2VirtualPhysicalTransM: np.ndarray = field(default_factory=get_empty_np_array)
cerrToDcmTransM: np.ndarray = field(default_factory=get_empty_np_array)
def __getitem__(self, key):
"""Return the value of the named attribute.
Args:
key (str): Attribute name to retrieve.
Returns:
Any: Value of the requested attribute.
"""
return getattr(self, key)
def __setitem__(self, key, value):
"""Set the value of the named attribute.
Args:
key (str): Attribute name to set.
value (Any): Value to assign to the attribute.
"""
return setattr(self, key, value)
[docs]
class json_serialize(json.JSONEncoder):
[docs]
def default(self, obj):
"""Serialize a Scan instance to a JSON-compatible dictionary.
Args:
obj (Any): Object to serialize. When the object is an instance
of :class:`Scan` its ``scanUID`` is returned; all other types
fall back to an empty string.
Returns:
str | dict: ``{'scan': obj.scanUID}`` for Scan instances,
otherwise an empty string.
"""
if isinstance(obj, Scan):
return {'scan':obj.scanUID}
return "" #json.JSONEncoder.default(self, obj)
[docs]
def getScanArray(self):
""" Routine to obtain image in the units defined in planC.scan[scanNum].scanInfo[slcNum].imageUnits
Returns:
np.ndarray: CTOffset is added to to scanArray such that the resulting array is in
real world units such as HU, SUV
"""
scan3M = self.scanArray - self.scanInfo[0].CTOffset
return scan3M
[docs]
def getNiiAffine(self):
""" Routine for affine transformation of pyCERR scan object for storing in NifTi format
Returns:
np.ndarray: 3x3 affine matrix
"""
# https://neurostars.org/t/direction-orientation-matrix-dicom-vs-nifti/14382/2
affine3M = self.Image2PhysicalTransM.copy()
affine3M[0,:] = -affine3M[0,:] * 10 #nii row is reverse of dicom, cm to mm
affine3M[1,:] = -affine3M[1,:] * 10 #nii col is reverse of dicom, cm to mm
affine3M[2,:] = affine3M[2,:] * 10 # cm to mm
return affine3M
[docs]
def saveNii(self, niiFileName):
""" Routine to save pyCERR Scan object to NifTi file
Args:
niiFileName (str): File name including the full path to save the pyCERR scan object to NifTi file.
Returns:
int: 0 when NifTi file is written successfully.
"""
img = self.getSitkImage()
sitk.WriteImage(img, niiFileName)
return 0
# affine3M = self.getNiiAffine()
# scan3M = self.getScanArray()
# scan3M = np.moveaxis(scan3M,[0,1],[1,0])
# #scan3M = np.flip(scan3M,axis=[0,1]) # negated affineM to take care of reverse row/col compared to dicom
# # Determine whether CERR slice order matches DICOM
# # dcmImgOri = self.scanInfo[0].imageOrientationPatient
# # slice_normal = dcmImgOri[[1,2,0]] * dcmImgOri[[5,3,4]] \
# # - dcmImgOri[[2,0,1]] * dcmImgOri[[4,5,3]]
# # slice_normal = slice_normal.reshape((1,3))
# # zDiff = np.matmul(slice_normal, self.scanInfo[1].imagePositionPatient) - np.matmul(slice_normal, self.scanInfo[0].imagePositionPatient)
# # ippDiffV = self.scanInfo[1].imagePositionPatient - self.scanInfo[0].imagePositionPatient
# if flipSliceOrderFlag(self): #np.all(np.sign(zDiff) < 0):
# scan3M = np.flip(scan3M,axis=2) # CERR slice ordering is opposite of DICOM
# img = nib.Nifti1Image(scan3M, affine3M)
# success = nib.save(img, niiFileName)
# return success
[docs]
def getSitkImage(self):
""" Routine to convert pyCERR Scan object to SimpleITK Image object
Returns:
sitk.Image: SimpleITK Image
"""
#sitkArray = np.moveaxis(self.getScanArray(),[0,1,2],[1,2,0])
sitkArray = np.transpose(self.getScanArray(), (2, 0, 1)) # z,y,x order
# CERR slice ordering is opposite of DICOM
if flipSliceOrderFlag(self):
sitkArray = np.flip(sitkArray, axis = 0)
originXyz = list(np.matmul(self.Image2PhysicalTransM, np.asarray([0,0,0,1]).T)[:3] * 10)
xV, yV, zV = self.getScanXYZVals()
dx = np.abs(xV[1] - xV[0]) * 10
dy = np.abs(yV[1] - yV[0]) * 10
dz = np.abs(zV[1] - zV[0]) * 10
spacing = [dx, dy, dz]
img_ori = self.scanInfo[0].imageOrientationPatient
slice_normal = img_ori[[1,2,0]] * img_ori[[5,3,4]] \
- img_ori[[2,0,1]] * img_ori[[4,5,3]]
# Get row-major directions for ITK
dir_cosine_mat = np.hstack((img_ori.reshape(3,2,order="F"),slice_normal.reshape(3,1)))
direction = dir_cosine_mat.reshape(9,order='C')
img = sitk.GetImageFromArray(sitkArray)
img.SetOrigin(originXyz)
img.SetSpacing(spacing)
img.SetDirection(direction)
return img
[docs]
def getScanXYZVals(self):
""" Routine to obtain pyCERR scan object's x,y,z grid coordinates. The coordinates are in pyCERR's
virtual coordinate system.
Returns:
tuple: x, y, z coordinates corresponding to the columns, rows, slices of scan voxels
"""
scan_info = self.scanInfo[0]
sizeDim1 = scan_info.sizeOfDimension1-1
sizeDim2 = scan_info.sizeOfDimension2-1
# Add gridUnits/2 to the last value to account for numerical noise
# Calculate xVals
xvals = np.arange(scan_info.xOffset - (sizeDim2 * scan_info.grid2Units) / 2,
scan_info.xOffset + (sizeDim2 * scan_info.grid2Units) / 2 + scan_info.grid2Units/2,
scan_info.grid2Units)
# Calculate yVals (flipped left-right)
yvals = np.arange(scan_info.yOffset + (sizeDim1 * scan_info.grid1Units) / 2,
scan_info.yOffset - (sizeDim1 * scan_info.grid1Units) / 2 - scan_info.grid1Units/2,
-scan_info.grid1Units)
# Extract zValues from the scanStruct dictionary or object
zvals = np.asarray([si.zValue for si in self.scanInfo])
return (xvals,yvals,zvals)
[docs]
def getScanSize(self):
""" Routine to get scan dimensions.
Returns:
np.array: numRows, numCols, numSlcs of pyCERR scan object
"""
numRows, numCols, numSlcs = self.scanInfo[0].sizeOfDimension1, self.scanInfo[0].sizeOfDimension2, \
len(self.scanInfo)
return np.asarray([numRows, numCols, numSlcs])
[docs]
def getScanOrientation(self):
""" Routine to get orientation of sacn w.r.t. patient.
Returns:
str: 3-character String representing the orientation of Scans's row, column and slice.
"""
orientPos = ['L', 'P', 'S']
orientNeg = ['R', 'A', 'I']
flipDict = {}
for i in range(len(orientPos)):
flipDict[orientPos[i]] = orientNeg[i]
flipDict[orientNeg[i]] = orientPos[i]
img_ori = self.scanInfo[0].imageOrientationPatient
img_ori = img_ori.reshape(6,1)
slice_normal = img_ori[[1,2,0]] * img_ori[[5,3,4]] \
- img_ori[[2,0,1]] * img_ori[[4,5,3]]
slice_normal = slice_normal.reshape((1,3))
# img_ori = np.vstack((img_ori, slice_normal.reshape((3,1))))
# dir_cosine_mat = img_ori.reshape(3, 3,order="F")
# itk_direction = dir_cosine_mat.reshape(9, order="C")
itk_direction = getITKDirection(self)
itk_orient_str = sitk.DICOMOrientImageFilter_GetOrientationFromDirectionCosines(itk_direction)
zDiff = np.matmul(slice_normal, self.scanInfo[1].imagePositionPatient) - np.matmul(slice_normal, self.scanInfo[0].imagePositionPatient)
ippDiffV = self.scanInfo[1].imagePositionPatient - self.scanInfo[0].imagePositionPatient
if np.all(np.sign(zDiff) < 0):
# cerr slice direction is opposite to ITK/DICOM order. Hence, flip.
zOri = flipDict[itk_orient_str[-1]]
else:
# cerr slice direction is opposite to ITK/DICOM order
zOri = itk_orient_str[-1]
orientString = itk_orient_str[:2]
orientString = orientString + zOri
return orientString
[docs]
def getScanSpacing(self):
""" Routine to get voxel spacing in cm.
Returns:
np.array: 3-element array containing dx, dy, dz of scan
"""
x_vals_v, y_vals_v, z_vals_v = self.getScanXYZVals()
if y_vals_v[0] > y_vals_v[1]:
y_vals_v = np.flip(y_vals_v)
dx = abs(np.median(np.diff(x_vals_v)))
dy = abs(np.median(np.diff(y_vals_v)))
dz = abs(np.median(np.diff(z_vals_v)))
spacing_v = np.array([dx, dy, dz])
return spacing_v
[docs]
def convertDcmToCerrVirtualCoords(self):
"""Routine to get scan from DICOM to pyCERR virtual coordinates. More information
about virtual coordinates is on the Wiki https://github.com/cerr/pyCERR/wiki/Coordinate-system
"""
# Construct DICOM Affine transformation matrix
# To construct DICOM affine transformation matrix it is necessary to figure out
# whether CERR slice direction matches DICOM to get the position of the 1st slice
# according to DICOM convention. Since slices are sorted according to decreasing order of
# dot product between ImagePositionPatient and ImageOrientationPatient.
#
# Determining order of scanArray slices
# If (slice_normal . ipp_2nd_slice - slice_normal . ipp_1st_slice) > 0,
# then DICOM slice order is reverse of CERR scanArray and scanInfo.
# i.e. the 1st slice in DICOM will correspond to the last slice in
# scanArray and the last element in scanInfo.
# Compute slice normal
dcmImgOri = self.scanInfo[0].imageOrientationPatient
dcmImgOri = dcmImgOri.reshape(6,1)
# slice_normal = dcmImgOri[[1,2,0]] * dcmImgOri[[5,3,4]] \
# - dcmImgOri[[2,0,1]] * dcmImgOri[[4,5,3]]
# slice_normal = slice_normal.reshape((1,3))
# zDiff = np.matmul(slice_normal, self.scanInfo[1].imagePositionPatient) - np.matmul(slice_normal, self.scanInfo[0].imagePositionPatient)
# ippDiffV = self.scanInfo[1].imagePositionPatient - self.scanInfo[0].imagePositionPatient
if flipSliceOrderFlag(self): # np.all(np.sign(zDiff) < 0):
info1 = self.scanInfo[-1]
info2 = self.scanInfo[-2]
else:
info1 = self.scanInfo[0]
info2 = self.scanInfo[1]
pos1V = info1.imagePositionPatient / 10 # cm
pos2V = info2.imagePositionPatient / 10 # cm
deltaPosV = pos2V - pos1V
pixelSpacing = [info1.grid2Units, info1.grid1Units]
# Transformation for DICOM Image to DICOM physical coordinates
# Pt coordinate to DICOM image coordinate mapping
# Based on ref: https://nipy.org/nibabel/dicom/dicom_orientation.html
position_matrix = np.hstack((np.matmul(dcmImgOri.reshape(3, 2,order="F"),np.diag(pixelSpacing)),
np.array([[deltaPosV[0], pos1V[0]], [deltaPosV[1], pos1V[1]], [deltaPosV[2], pos1V[2]]])))
position_matrix = np.vstack((position_matrix, np.array([0, 0, 0, 1])))
#positionMatrixInv = np.linalg.inv(position_matrix)
self.Image2PhysicalTransM = position_matrix
# Get DICOM x,y,z coordinates of the center voxel.
# This serves as the reference point for the image volume.
sizV = self.scanArray.shape
xyzCtrV = position_matrix * np.array([(sizV[1] - 1) / 2, (sizV[0] - 1) / 2, 0, 1])
xOffset = np.sum(np.matmul(np.transpose(dcmImgOri[:3,:]), xyzCtrV[:3]))
yOffset = -np.sum(np.matmul(np.transpose(dcmImgOri[3:,:]), xyzCtrV[:3])) # (-)ve since CERR y-coordinate is opposite of column vector.
for i in range(len(self.scanInfo)):
self.scanInfo[i].xOffset = xOffset
self.scanInfo[i].yOffset = yOffset
xs, ys, zs = self.getScanXYZVals()
dx = xs[1] - xs[0]
dy = ys[1] - ys[0]
slice_distance = zs[1] - zs[0]
# Transformation for DICOM Image to CERR physical coordinates
# DICOM 1st slice is CERR's last slice (i.e. zs[-1]
if flipSliceOrderFlag(self): #np.all(np.sign(zDiff) < 0):
virPosMtx = np.array([[dx, 0, 0, xs[0]], [0, dy, 0, ys[0]], [0, 0, -slice_distance, zs[-1]], [0, 0, 0, 1]])
else:
virPosMtx = np.array([[dx, 0, 0, xs[0]], [0, dy, 0, ys[0]], [0, 0, slice_distance, zs[0]], [0, 0, 0, 1]])
self.Image2VirtualPhysicalTransM = virPosMtx
# Construct transformation matrix to convert cerr-xyz to dicom-xyz
self.cerrToDcmTransM = np.matmul(self.Image2PhysicalTransM, np.linalg.inv(self.Image2VirtualPhysicalTransM))
self.cerrToDcmTransM[:,:3] = self.cerrToDcmTransM[:,:3] * 10 # cm to mm
[docs]
def convertDcmToRealWorldUnits(self, opts={}):
""" Routine to convert pixel array from DICOM storage units to real world units.
Args:
opts (dict): Dictionary of options to convert to real world units. Currrently, only one option
if supported - importMRPreciseValueFlag (yes or no) to specify whether to convert MR image from
Philips scanner to precise values.
"""
importMRPreciseValueFlag = 'no'
if 'importMRPreciseValueFlag' in opts:
importMRPreciseValueFlag = opts['importMRPreciseValueFlag']
# Apply ReScale Intercept and Slope
scanArray3M = np.zeros(self.scanArray.shape, dtype=np.float32)
numSlcs = self.scanArray.shape[2]
rescaleSlopeV = np.ones(numSlcs)
realWorldImageFlag = False
for slcNum in range(numSlcs):
rescaleSlope = self.scanInfo[slcNum].rescaleSlope
rescaleIntrcpt = self.scanInfo[slcNum].rescaleIntercept
realWorldValueSlope = self.scanInfo[slcNum].realWorldValueSlope
realWorldValueIntercept = self.scanInfo[slcNum].realWorldValueIntercept
realWorldMeasurCodeMeaning = self.scanInfo[slcNum].realWorldMeasurCodeMeaning
philipsImageUnits = self.scanInfo[slcNum].philipsImageUnits
rescaleType = self.scanInfo[slcNum].rescaleType
manufacturer = self.scanInfo[slcNum].manufacturer
if manufacturer.lower() in ['philips'] and \
realWorldValueSlope is not None and \
not np.isnan(realWorldValueSlope) and \
self.scanInfo[slcNum].imageType.lower() in ['mr scan'] and \
realWorldMeasurCodeMeaning is not None:
realWorldImageFlag = True
scanArray3M[:, :, slcNum] = \
self.scanArray[:, :, slcNum] * realWorldValueSlope + realWorldValueIntercept
self.scanInfo[slcNum].imageUnits = realWorldMeasurCodeMeaning
else:
scanArray3M[:, :, slcNum] = \
self.scanArray[:, :, slcNum] * rescaleSlope + rescaleIntrcpt
if len(self.scanInfo[slcNum].imageUnits) == 0 and \
self.scanInfo[slcNum].imageType.lower() not in ['pt scan', 'nm scan']:
self.scanInfo[slcNum].imageUnits = rescaleType
rescaleSlopeV[slcNum] = rescaleSlope
minScanVal = np.min(scanArray3M)
ctOffset = max(0, -minScanVal)
scanArray3M += ctOffset
# Decommissioned conversion to unsigned int. Need to update logic to handle various data types - dicom, nii etc.
# minScanVal = np.min(scanArray3M)
# maxScanVal = np.max(scanArray3M)
# if not realWorldImageFlag and not np.any(np.abs(rescaleSlopeV - 1) > np.finfo(float).eps * 1e5):
# # Convert to uint if rescale slope is not 1
# if minScanVal >= -32768 and maxScanVal <= 32767:
# scanArray3M = scanArray3M.astype(np.uint16)
# else:
# scanArray3M = scanArray3M.astype(np.uint32)
for slcNum in range(numSlcs):
self.scanInfo[slcNum].CTOffset = ctOffset
self.scanArray = scanArray3M
# Convert Philips MR to precise values
if self.scanInfo[slcNum].imageType.lower() == 'mr scan' and \
importMRPreciseValueFlag.lower() == 'yes':
# Ref: Chenevert, Thomas L., et al. "Errors in quantitative image analysis due to platform-dependent image scaling."
manufacturer = self.scanInfo[0].manufacturer
if 'philips' in manufacturer.lower() and \
self.scanInfo[0].scaleSlope is not None and \
not realWorldImageFlag:
scaleSlope = self.scanInfo[0].scaleSlope
self.scanArray = self.scanArray.astype(np.float32) / (rescaleSlope * scaleSlope)
[docs]
def convertToSUV(self, suvType=None):
""" Routine to convert pixel array for PET scan from DICOM storage to SUV
Args:
suvType (str): optional, type of SUV. When not specified, the suvType is read from DICOM if available.
When not specified and not available in DIOCM, a default value of 'BW' is used. Currently supported
options are 'BW', 'BSA', 'LBM', 'LBMJANMA'
"""
scan3M = self.scanArray
headerS = self.scanInfo
scanSiz = scan3M.shape
suv3M = np.zeros(scanSiz)
numSlcs = scan3M.shape[2]
acqTimeV = np.empty(numSlcs,dtype=float)
# Allow None input for suvType
if suvType == "" or suvType == None:
suvType = headerS[0].suvType
if suvType == "":
suvType = 'BW'
for slcNum in range(numSlcs):
headerSlcS = headerS[slcNum]
if headerSlcS.acquisitionTime:
acqTimeV[slcNum] = dcm_hhmmss(headerSlcS.acquisitionTime)[0]
seriesTime = dcm_hhmmss(headerSlcS.seriesTime)[0]
seriesDate = np.nan
if headerSlcS.seriesDate:
seriesDate = dcm_to_np_date(headerSlcS.seriesDate)
injectionDate = np.nan
if headerSlcS.injectionDate:
injectionDate = dcm_to_np_date(headerSlcS.injectionDate)
acqStartTime = np.nan
if not np.any(np.isnan(acqTimeV)):
acqStartTime = np.min(acqTimeV)
for slcNum in range(scan3M.shape[2]):
headerSlcS = headerS[slcNum]
imgM = scan3M[:, :, slcNum] - headerSlcS.CTOffset
imgUnits = headerSlcS.imageUnits
imgMUnits = imgM.copy()
if imgUnits == 'CNTS':
activityScaleFactor = headerSlcS.philipsActivityConcentrationScaleFactor
suvScaleFactor = headerSlcS.philipsSUVScaleFactor
if activityScaleFactor != "":
imgMUnits = imgMUnits * activityScaleFactor
imgMUnits = imgMUnits * 1000 # Bq/L
imgUnits = 'BQL'
elif suvScaleFactor != "":
suvM = imgMUnits * suvScaleFactor # GML
imageUnits = 'GML'
elif imgUnits in ['BQML', 'BQCC']:
imgMUnits = imgMUnits * 1000 # Bq/L
imgUnits = 'BQL'
elif imgUnits in ['KBQCC', 'KBQML']:
imgMUnits = imgMUnits * 1e6 # Bq/L
imgUnits = 'BQL'
elif imgUnits in ['GML', 'CM2ML']:
suvM = imgMUnits
imageUnits = imgUnits
else:
#raise ValueError('SUV calculation is supported only for imageUnits BQML and CNTS')
import warnings
warnings.warn("'SUV calculation is supported only for imageUnits BQML and CNTS'")
return
if imgUnits == 'BQL':
decayCorrection = headerSlcS.petDecayCorrection
if len(headerSlcS.petDecayCorrectionDateTime) > 8:
scantime = dcm_hhmmss(headerSlcS.petDecayCorrectionDateTime[8:])[0]
elif len(headerSlcS.gePETDecayCorrectionDateTime) > 8:
scantime = dcm_hhmmss(headerSlcS.gePETDecayCorrectionDateTime[8:])[0]
elif len(headerSlcS.siemensPETDecayCorrectionDateTime) > 8:
scantime = dcm_hhmmss(headerSlcS.siemensPETDecayCorrectionDateTime[8:])[0]
elif decayCorrection == 'START':
scantime = seriesTime
if not np.isnan(acqStartTime) and acqStartTime < scantime:
scantime = acqStartTime
elif decayCorrection == 'ADMIN':
scantime = dcm_hhmmss(headerSlcS.injectionTime)[0]
elif decayCorrection == 'NONE':
scantime = np.nan
else:
scantime = np.nan
# Start Time for Radiopharmaceutical Injection
injection_time = dcm_hhmmss(headerSlcS.injectionTime)[0]
if not np.isnan(seriesDate) and not np.isnan(injectionDate):
date_diff = seriesDate - injectionDate
if date_diff < 5: # check whether it is a reasonable value
injection_time = injection_time - date_diff.item().total_seconds()
# Half Life for Radionuclide
half_life = headerSlcS.halfLife
# Total dose injected for Radionuclide
injected_dose = headerSlcS.injectedDose
# Modality
modality = headerSlcS.imageType
if modality.upper() == 'NM SCAN':
injected_dose = injected_dose * 1e6 # Convert MBq to Bq
# Fix issue where IOD is PT and injected_dose units are in MBq
if injected_dose < 1e5:
injected_dose = injected_dose * 1e6
# Calculate the decay
# The injected dose used to calculate suvM is corrected for the decay that
# occurs between the time of injection and the time of scan.
# decayFactor = e^(t1-t2/halflife)
if decayCorrection.upper() == 'NONE':
decay = 1
else:
decay = np.exp(-np.log(2) * (scantime - injection_time) / half_life)
# Calculate the dose decayed during procedure
injected_dose_decay = injected_dose * decay # in Bq
# Patient Weight
ptWeight = headerSlcS.patientWeight
# Calculate SUV based on type
# reference: http://dicom.nema.org/medical/Dicom/2017e/output/chtml/part16/sect_CID_85.html
# SUVbw and SUVbsa equations are taken from Kim et al. Journal of Nuclear Medicine. Volume 35, No. 1, January 1994. pp 164-167.
suvType = suvType.upper()
if suvType == 'BW': # Body Weight
suvM = imgMUnits * ptWeight / injected_dose_decay # pt weight in grams
imageUnits = 'GML'
elif suvType == 'BSA': # body surface area
# Patient height
# (BSA in m2) = [(weight in kg)^0.425 * (height in cm)^0.725 * 0.007184].
# SUV-bsa = (PET image Pixels) * (BSA in m2) * (10000 cm2/m2) / (injected dose).
ptHeight = headerSlcS.patientSize # units of meter
bsaMm = ptWeight**0.425 * (ptHeight * 100)**0.725 * 0.007184
suvM = imgMUnits * bsaMm / injected_dose_decay
imageUnits = 'CM2ML'
elif suvType == 'LBM': # lean body mass by James method
ptGender = headerSlcS.patientSex
ptHeight = headerSlcS.patientSize
if ptGender.upper() == 'M':
# LBM in kg = 1.10 * (weight in kg) - 120 * [(weight in kg) / (height in cm)]^2.
lbmKg = 1.10 * ptWeight - 120 * (ptWeight / (ptHeight * 100))**2
else:
# if gender == female
# LBM in kg = 1.07 * (weight in kg) - 148 * [(weight in kg) / (height in cm)]^2.
lbmKg = 1.07 * ptWeight - 148 * (ptWeight / (ptHeight * 100))**2
suvM = imgMUnits * lbmKg / injected_dose_decay
imageUnits = 'GML'
elif suvType == 'LBMJAMES128': # lean body mass by James method
imageUnits = 'GML'
elif suvType == 'LBMJANMA': # lean body mass by Janmahasatian method
ptHeight = headerSlcS['patientSize']
bmi = (ptWeight * 2.20462 / (ptHeight * 39.3701)**2) * 703
ptGender = headerSlcS['patientSex']
if ptGender.upper() == 'M':
lbmKg = (9270 * ptWeight) / (6680 + 216 * bmi) # male
else:
lbmKg = (9270 * ptWeight) / (8780 + 244 * bmi) # female
suvM = imgMUnits * lbmKg / injected_dose_decay
imageUnits = 'GML'
elif suvType == 'IBW': # ideal body weight
imageUnits = 'GML'
else:
return
suv3M[:, :, slcNum] = suvM
self.scanInfo[slcNum].imageUnits = imageUnits
self.scanInfo[slcNum].suvType = suvType
self.scanArray = suv3M
return
[docs]
def getScanDict(self):
""" Routine to get dictionary representation of scan metadata
Returns:
dict: fields of the dictionary are attributes of the Scan object.
"""
scanDict = self.__dict__.copy()
sInfoList = []
for sInfo in scanDict['scanInfo']:
sInfoDict = sInfo.__dict__.copy()
sInfoList.append(sInfoDict)
scanDict['scanInfo'] = sInfoList
return scanDict
[docs]
def getDcmScanInfo(self) -> list:
"""This routine return a list of scanInfo from scan object. The order corresponds to scanArray slice dimension
Returns:
list: list of scanInfo dictionaries
"""
scnDict = self.getScanDict()
scanInfoList = scnDict['scanInfo']
cerrSpecificKeyList = ['xOffset', 'yOffset', 'CTScale', 'distrustAbove',
'imageSource', 'transferProtocol', 'studyNumberOfOrigin',
'tapeOfOrigin', 'scanFileName', 'unitNumber', 'CTAir', 'CTWater',
'zValue','imageNumber','caseNumber','CTOffset','scanType',
'numberRepresentation', 'numberOfDimensions', 'voxelThickness',
'headInOut', 'positionInScan', 'patientAttitude', 'frameAcquisitionDuration',
'frameReferenceDateTime', 'scanNumber', 'scanID', 'patientBirthDate',
'scaleSlope', 'scaleIntercept', 'siteOfInterest','scanDate', 'patientPosition',
'philipsImageUnits', 'philipsRescaleSlope', 'philipsRescaleIntercept']
for scanInfoDict in scanInfoList:
for cerrKey in cerrSpecificKeyList:
del scanInfoDict[cerrKey]
# Create a new dictionary with DICOM names for keys
scanDir = os.path.dirname(__file__)
mappingFile = os.path.join(scanDir, 'dcm_cerr_name_map.json')
with open(mappingFile, 'r') as nameMapFile:
dcmCerrNameMap = json.load(nameMapFile)
dcmNameScanInfoList = []
for scanInfoDict in scanInfoList:
dcmNameScanInfoDict = dict('')
for key in scanInfoDict.keys():
if key in dcmCerrNameMap:
dcmNameScanInfoDict[dcmCerrNameMap[key]] = scanInfoDict[key]
if isinstance(dcmNameScanInfoDict[dcmCerrNameMap[key]], np.ndarray):
dcmNameScanInfoDict[dcmCerrNameMap[key]] = dcmNameScanInfoDict[dcmCerrNameMap[key]].tolist()
dcmNameScanInfoDict["PixelSpacing"] = [scanInfoDict['grid2Units'], scanInfoDict['grid1Units']]
dcmNameScanInfoList.append(dcmNameScanInfoDict)
return dcmNameScanInfoList #scanInfoList
[docs]
def flipSliceOrderFlag(scan):
""" Routine to determine slice order for determining the origin for conversion to NifTi and SimpleITK formats.
Args:
scan (cerr.dataclasses.scan.Scan): pyCERR scan object
Returns:
bool: True when dot product of slice normal and imagePositionPatient increases with slice order
"""
dcmImgOri = scan.scanInfo[0].imageOrientationPatient
slice_normal = dcmImgOri[[1,2,0]] * dcmImgOri[[5,3,4]] \
- dcmImgOri[[2,0,1]] * dcmImgOri[[4,5,3]]
slice_normal = slice_normal.reshape((1,3))
zDiff = np.matmul(slice_normal, scan.scanInfo[1].imagePositionPatient) - np.matmul(slice_normal, scan.scanInfo[0].imagePositionPatient)
ippDiffV = scan.scanInfo[1].imagePositionPatient - scan.scanInfo[0].imagePositionPatient
return np.all(np.sign(zDiff) < 0)
[docs]
def getITKDirection(scan):
"""
Args:
scan (cerr.dataclasses.scan.Scan): pyCERR scan object
Returns:
np.ndarray: 9-element array of direction cosines of row, column and slice w.r.t. patient.
"""
img_ori = scan.scanInfo[0].imageOrientationPatient
img_ori = img_ori.reshape(6,1)
slice_normal = img_ori[[1,2,0]] * img_ori[[5,3,4]] \
- img_ori[[2,0,1]] * img_ori[[4,5,3]]
slice_normal = slice_normal.reshape((1,3))
img_ori = np.vstack((img_ori, slice_normal.reshape((3,1))))
dir_cosine_mat = img_ori.reshape(3, 3,order="F")
itk_direction = dir_cosine_mat.reshape(9, order="C")
return itk_direction
[docs]
def dcm_hhmmss(time_str):
"""Parse a DICOM time string (HHMMSS) into its components and total seconds.
Args:
time_str (str): DICOM-format time string with at least 6 characters in
HHMMSS order (fractional seconds are accepted but currently ignored).
Returns:
tuple: A 5-element tuple ``(totSec, hh, mm, ss, fract)`` where
``totSec`` (int) is the total number of seconds since midnight,
``hh`` (int) is hours, ``mm`` (int) is minutes, ``ss`` (int) is
seconds, and ``fract`` is ``None`` (reserved for future use).
"""
hh = int(time_str[0:2])
mm = int(time_str[2:4])
ss = int(time_str[4:6])
fract = None # You can implement this part if needed
totSec = hh * 3600 + mm * 60 + ss
return totSec, hh, mm, ss, fract
[docs]
def dcm_to_np_date(dateStr):
"""Convert a DICOM date string (YYYYMMDD) to a numpy datetime64 object.
Args:
dateStr (str): DICOM-format date string, exactly 8 characters in
YYYYMMDD order.
Returns:
np.datetime64 | None: A ``numpy.datetime64`` with day precision when
``dateStr`` is 8 characters long, otherwise ``None``.
"""
dateObj = None
if len(dateStr) == 8:
dateObj = np.datetime64(dateStr[:4] + "-" + dateStr[4:6] + "-" + dateStr[6:], 'D')
return dateObj
[docs]
def get_slice_position(scan_info_item):
"""Extract the z-position from an enumerated scan-info pair for sort key use.
Args:
scan_info_item (tuple): A ``(index, ScanInfo)`` pair as produced by
``enumerate``. The second element must have a ``zValue`` attribute.
Returns:
float: The ``zValue`` of the scan-info slice, in centimetres.
"""
return scan_info_item[1].zValue
[docs]
def populateScanInfoFields(s_info, ds):
"""
Args:
s_info (cerr.dataclasses.scan_info.ScanInfo): pyCERR's scanInfo object for storing metadata per slice.
ds (pydicom.dataset.Dataset): pydicom dataset object
Returns:
cerr.dataclasses.scan_info.ScanInfo: scanInfo object with attributes populated from metadata from input ds.
"""
s_info.frameOfReferenceUID = ds.FrameOfReferenceUID
s_info.imageType = ds.Modality
if not "SCAN" in s_info.imageType.upper():
s_info.imageType = s_info.imageType + " SCAN"
if hasattr(ds,"SeriesDescription"): s_info.seriesDescription = ds.SeriesDescription
if hasattr(ds,"ManufacturerModelName"): s_info.scannerType = ds.ManufacturerModelName
if hasattr(ds,"Manufacturer"): s_info.manufacturer = ds.Manufacturer
s_info.scanFileName = ds.filename
s_info.sopInstanceUID = ds.SOPInstanceUID
s_info.sopClassUID = ds.SOPClassUID
s_info.seriesInstanceUID = ds.SeriesInstanceUID
s_info.studyInstanceUID = ds.StudyInstanceUID
s_info.bitsAllocated = ds.BitsAllocated
s_info.bitsStored = ds.BitsStored
s_info.pixelRepresentation = ds.PixelRepresentation
s_info.sizeOfDimension1 = ds.Rows
s_info.sizeOfDimension2 = ds.Columns
if hasattr(ds,"PatientName"): s_info.patientName = str(ds.PatientName)
if hasattr(ds,"PatientID"): s_info.patientID = ds.PatientID
if hasattr(ds,"AcquisitionDate"): s_info.acquisitionDate = ds.AcquisitionDate
if hasattr(ds,"AcquisitionTime"): s_info.acquisitionTime = ds.AcquisitionTime
if hasattr(ds,"SeriesDate"): s_info.seriesDate = ds.SeriesDate
if hasattr(ds,"SeriesTime"): s_info.seriesTime = ds.SeriesTime
if hasattr(ds,"StudyDate"): s_info.studyDate = ds.StudyDate
if hasattr(ds,"StudyTime"): s_info.studyTime = ds.StudyTime
if hasattr(ds,"StudyDescription"): s_info.studyDescription = ds.StudyDescription
if hasattr(ds,"PatientWeight"): s_info.patientWeight = ds.PatientWeight
if hasattr(ds,"PatientSize"): s_info.patientSize = ds.PatientSize
if ("0010","1022") in ds: s_info.patientBmi = ds["0010","1022"].value
if hasattr(ds,"PatientSex"): s_info.patientSex = ds.PatientSex
if hasattr(ds,"SeriesType"): s_info.petSeriesType = ds.SeriesType
if hasattr(ds,"Units"): s_info.imageUnits = ds.Units
if hasattr(ds,"CountsSource"): s_info.petCountSource = ds.CountsSource
if hasattr(ds,"NumberOfSlices"): s_info.petNumSlices = ds.NumberOfSlices
if hasattr(ds,"DecayCorrection"): s_info.petDecayCorrection = ds.DecayCorrection
if hasattr(ds,"CorrectedImage"): s_info.petCorrectedImage = ds.CorrectedImage
if ("0054","1006") in ds: s_info.suvType = ds["0054","1006"].value
if hasattr(ds,"WindowCenter"): s_info.windowCenter = ds.WindowCenter
if hasattr(ds,"WindowWidth"): s_info.windowWidth = ds.WindowWidth
if ("2005","140B") in ds: s_info.philipsImageUnits = ds["2005","140B"].value
if ("2005","140A") in ds: s_info.philipsRescaleSlope = ds["2005","140A"].value
if ("2005","1409") in ds: s_info.philipsRescaleIntercept = ds["2005","1409"].value
if hasattr(ds,"PatientIdentityRemoved"): s_info.patientIdentityRemoved = ds.PatientIdentityRemoved
if hasattr(ds,"DeidentificationMethod"): s_info.deIdentificationMethod = ds.DeidentificationMethod
if hasattr(ds,"DeidentificationMethodCodeSequence"):
for deIdMethod in ds.DeidentificationMethodCodeSequence:
methodStr = deIdMethod.CodeValue + ': ' + deIdMethod.CodeMeaning
s_info.deidentificationMethodDescription = np.append(s_info.deidentificationMethodDescription, methodStr)
if ("0018","0010") in ds: s_info.contrastBolusAgent = ds["0018","0010"].value
return s_info
[docs]
def populateRealWorldFields(s_info, perFrameSeq):
"""
Args:
s_info (cerr.dataclasses.scan_info.ScanInfo): pyCERR's scanInfo object for storing metadata per slice.
perFrameSeq (pydicom.dataset.Dataset): pydicom dataset object or ds.PerFrameFunctionalGroupsSequence
for multiFrameFlg images.
Returns:
cerr.dataclasses.scan_info.ScanInfo: scanInfo object with attributes populated from metadata from input ds.
"""
if 'RealWorldValueMappingSequence' in perFrameSeq:
RealWorldValueMappingSeq = perFrameSeq.RealWorldValueMappingSequence[0]
if hasattr(RealWorldValueMappingSeq,'RealWorldValueSlope'):
s_info.realWorldValueSlope = RealWorldValueMappingSeq.RealWorldValueSlope
if hasattr(RealWorldValueMappingSeq,'RealWorldValueIntercept'):
s_info.realWorldValueIntercept = RealWorldValueMappingSeq.RealWorldValueIntercept
if ("0040","08EA") in RealWorldValueMappingSeq:
if ("0008","0100") in RealWorldValueMappingSeq["0040","08EA"][0]:
s_info.realWorldMeasurCodeMeaning = RealWorldValueMappingSeq["0040","08EA"][0]["0008","0100"].value
elif ("0008","0119") in RealWorldValueMappingSeq["0040","08EA"][0]:
s_info.realWorldMeasurCodeMeaning = RealWorldValueMappingSeq["0040","08EA"][0]["0008","0119"].value
if ("0008","0104") in RealWorldValueMappingSeq["0040","08EA"][0]:
s_info.realWorldMeasurCodeMeaning = RealWorldValueMappingSeq["0040","08EA"][0]["0008","0104"].value
return s_info
[docs]
def populateRadiopharmaFields(s_info, seq):
"""
Args:
s_info (cerr.dataclasses.scan_info.ScanInfo): pyCERR's scanInfo object for storing metadata per slice.
seq (pydicom.dataset.Dataset): dataset containing radiopharma metadata for PET scan.
Returns:
cerr.dataclasses.scan_info.ScanInfo: scanInfo object with attributes populated from metadata from input ds.
"""
# populate radiopharma info
if ("0054","0016") in seq:
radiopharmaInfoSeq = seq["0054","0016"].value[0]
if hasattr(radiopharmaInfoSeq,"RadiopharmaceuticalStartDateTime"):
s_info.injectionDate = radiopharmaInfoSeq.RadiopharmaceuticalStartDateTime[:8]
s_info.injectionTime = radiopharmaInfoSeq.RadiopharmaceuticalStartDateTime[8:]
elif hasattr(radiopharmaInfoSeq,"RadiopharmaceuticalStartTime"):
s_info.injectionTime = radiopharmaInfoSeq.RadiopharmaceuticalStartTime
s_info.injectedDose = float(radiopharmaInfoSeq.RadionuclideTotalDose)
s_info.halfLife = float(radiopharmaInfoSeq.RadionuclideHalfLife)
if ("7053","1009") in seq: s_info.philipsActivityConcentrationScaleFactor = seq["7053","1009"].value
if ("0018", "9701") in seq: s_info.petDecayCorrectionDateTime = seq["0018", "9701"].value
if ("0071","1022") in seq: s_info.siemensPETDecayCorrectionDateTime = seq["0071","1022"].value # Siemens
if ("0009","100D") in seq: s_info.gePETDecayCorrectionDateTime = seq["0009","100D"].value # GE
if ("7053","1000") in seq: s_info.philipsSUVScaleFactor = seq["7053","1000"].value
return s_info
[docs]
def parseScanInfoFields(ds, multiFrameFlg=False) -> (scn_info.ScanInfo, Dataset.pixel_array, str):
"""
Args:
ds (pydicom.dataset.Dataset): Dataset object read from DICOM file
multiFrameFlg (bool): True when dataset is multiFrame image, otherwise False.
Returns:
cerr.dataclasses.scan_info.ScanInfo: scanInfo object with attributes populated from metadata from input ds.
"""
#numberOfFrames = ds.NumberOfFrames.real
# s_info.frameOfReferenceUID = ds.FrameOfReferenceUID
#s_info.seriesDescription = ds.SeriesDescription
if not multiFrameFlg: #numberOfFrames == 1:
scan_info = scn_info.ScanInfo()
scan_info = populateScanInfoFields(scan_info, ds)
if hasattr(ds,'RescaleSlope'): scan_info.rescaleSlope = ds.RescaleSlope
if hasattr(ds,'RescaleIntercept'): scan_info.rescaleIntercept = ds.RescaleIntercept
if hasattr(ds,'RescaleType'): scan_info.rescaleType = ds.RescaleType
if ("2005","100E") in ds: scan_info.scaleSlope = ds["2005","100E"].value
if ("2005","100D") in ds: scan_info.scaleIntercept = ds["2005","100D"].value
scan_info = populateRealWorldFields(scan_info, ds)
scan_info.grid1Units = ds.PixelSpacing[1] / 10
scan_info.grid2Units = ds.PixelSpacing[0] / 10
scan_info.sliceThickness = ds.SliceThickness / 10
scan_info.imageOrientationPatient = np.array(ds.ImageOrientationPatient)
scan_info.imagePositionPatient = np.array(ds.ImagePositionPatient)
slice_normal = scan_info.imageOrientationPatient[[1,2,0]] * scan_info.imageOrientationPatient[[5,3,4]] \
- scan_info.imageOrientationPatient[[2,0,1]] * scan_info.imageOrientationPatient[[4,5,3]]
scan_info.zValue = - np.sum(slice_normal * scan_info.imagePositionPatient) / 10
bVal1 = ("0043", "1039") # GE
bVal2 = ("0018", "9087") # Philips
bVal3 = ("0019", "100C") # Siemens
temporalPos = ("0020", "0100")
triggerTime = ("0018", "1060")
TR = ("0018", "0080")
FA = ("0018", "1314")
if bVal1 in ds: scan_info.bValue = ds["0043", "1039"].value
if bVal2 in ds: scan_info.bValue = ds["0018", "9087"].value
if bVal3 in ds: scan_info.bValue = ds["0019", "100C"].value
if temporalPos in ds: scan_info.temporalPositionIndex = ds["0020", "0100"].value
if triggerTime in ds: scan_info.triggerTime = ds["0018", "1060"].value
if TR in ds: scan_info.repetitionTime = float(ds["0018", "0080"].value)
if FA in ds: scan_info.flipAngle = float(ds["0018", "1314"].value)
scan_info = populateRadiopharmaFields(scan_info, ds)
else:
numberOfFrames = ds.NumberOfFrames.real
scan_info = np.empty(numberOfFrames, dtype=scn_info.ScanInfo)
for iFrame in range(numberOfFrames):
s_info = scn_info.ScanInfo()
s_info = populateScanInfoFields(s_info, ds)
if 'PerFrameFunctionalGroupsSequence' in ds:
perFrameSeq = ds.PerFrameFunctionalGroupsSequence[iFrame]
s_info = populateRealWorldFields(s_info, perFrameSeq)
if 'PixelValueTransformationSequence' in perFrameSeq:
PixelValueTransformSeq = perFrameSeq.PixelValueTransformationSequence[0]
if hasattr(PixelValueTransformSeq,'RescaleSlope'): s_info.rescaleSlope = PixelValueTransformSeq.RescaleSlope
if hasattr(PixelValueTransformSeq,'RescaleIntercept'): s_info.rescaleIntercept = PixelValueTransformSeq.RescaleIntercept
if ("2005","100E") in PixelValueTransformSeq: s_info.scaleSlope = PixelValueTransformSeq["2005","100E"].value
if ("2005","100D") in PixelValueTransformSeq: s_info.scaleIntercept = PixelValueTransformSeq["2005","100D"].value
s_info.imagePositionPatient = np.array(perFrameSeq.PlanePositionSequence[0].ImagePositionPatient)
s_info.imageOrientationPatient = np.array(perFrameSeq.PlaneOrientationSequence[0].ImageOrientationPatient)
slice_normal = s_info.imageOrientationPatient[[1,2,0]] * s_info.imageOrientationPatient[[5,3,4]] \
- s_info.imageOrientationPatient[[2,0,1]] * s_info.imageOrientationPatient[[4,5,3]]
s_info.zValue = - np.sum(slice_normal * s_info.imagePositionPatient) / 10
if 'FrameVOILUTSequence' in perFrameSeq:
s_info.windowWidth = float(perFrameSeq.FrameVOILUTSequence[0].WindowWidth)
s_info.windowCenter = float(perFrameSeq.FrameVOILUTSequence[0].WindowCenter)
else: # NM scans
s_info = populateRealWorldFields(s_info, ds)
if 'DetectorInformationSequence' in ds:
imagePositionPatientStart = np.array(ds.DetectorInformationSequence[0].ImagePositionPatient)
s_info.imageOrientationPatient = np.array(ds.DetectorInformationSequence[0].ImageOrientationPatient)
sliceSpacing = ds.SpacingBetweenSlices * iFrame
slice_normal = s_info.imageOrientationPatient[[1,2,0]] * s_info.imageOrientationPatient[[5,3,4]] \
- s_info.imageOrientationPatient[[2,0,1]] * s_info.imageOrientationPatient[[4,5,3]]
s_info.imagePositionPatient = imagePositionPatientStart + slice_normal * sliceSpacing
s_info.zValue = - np.sum(slice_normal * s_info.imagePositionPatient) / 10
if 'SharedFunctionalGroupsSequence' in ds:
PixelSpacing = ds.SharedFunctionalGroupsSequence[0].PixelMeasuresSequence[0].PixelSpacing
s_info.sliceThickness = ds.SharedFunctionalGroupsSequence[0].PixelMeasuresSequence[0].SliceThickness / 10
else:
PixelSpacing = ds.PixelSpacing
if isinstance(ds.SliceThickness, (float, int)):
s_info.sliceThickness = ds.SliceThickness / 10
elif 'SpacingBetweenSlices' in ds and isinstance(ds.SpacingBetweenSlices, (float, int)):
s_info.sliceThickness = ds.SpacingBetweenSlices / 10
s_info.grid1Units = PixelSpacing[1] / 10
s_info.grid2Units = PixelSpacing[0] / 10
# MR-specific tags
TR = ("0018", "0080")
FA = ("0018", "1314")
if TR in ds: scan_info.repetitionTime = float(ds["0018", "0080"].value)
if FA in ds: scan_info.flipAngle = float(ds["0018", "1314"].value)
s_info = populateRadiopharmaFields(s_info, ds)
scan_info[iFrame] = s_info
return (scan_info, ds.pixel_array, ds.SeriesInstanceUID)
[docs]
def loadSortedScanInfo(file_list):
"""
Args:
file_list (list): list of files to read into pyCERR's Scan object
Returns:
cerr.daatclasses.scan.Scan: pyCERR scan object containing metadata from the file_list.
"""
scan = Scan()
#scan_info = [] #scn_info.ScanInfo()
#scan_array = []
scan_array = [] #np.empty(len(file_list))
scan_info = np.empty(len(file_list),dtype=scn_info.ScanInfo)
count = 0
multiFrameFlag = False
for file in file_list:
ds = dcmread(file)
if np.any(ds.Modality == np.array(["CT","PT", "MR", "NM"])): #hasattr(ds, "pixel_array"):
if len(file_list) == 1 and 'NumberOfFrames' in ds:
multiFrameFlag = True
si_pixel_data = parseScanInfoFields(ds, multiFrameFlag)
scan_array = np.transpose(si_pixel_data[1], (1,2,0))
scan_info = si_pixel_data[0]
count = len(scan_info)
else:
si_pixel_data = parseScanInfoFields(ds)
#scan_info.append(si_pixel_data[0])
#scan_array.append(si_pixel_data[1])
scan_info[count] = si_pixel_data[0]
if not isinstance(scan_array, np.ndarray) and not scan_array:
imgSiz = list(si_pixel_data[1].shape)
imgSiz.append(len(file_list))
scan_array = np.empty(imgSiz, dtype=np.float32)
scan_array[:,:,count] = si_pixel_data[1]
count += 1
if count < scan_array.shape[2]:
scan_array = np.delete(scan_array,np.arange(count,scan_array.shape[2]),axis=2)
scan_info = np.delete(scan_info,np.arange(count,scan_array.shape[2]),axis=0)
# Filter out duplicate SOP Instances
if np.any(ds.Modality == np.array(["CT","PT", "MR", "NM"])) and not multiFrameFlag:
allSOPs = [s.sopInstanceUID for s in scan_info]
uniqSOPs, uniqInds = np.unique(allSOPs, return_index=True)
duplicateIDs = list(set(range(len(scan_info))) - set(uniqInds))
scan_array = np.delete(scan_array,duplicateIDs,axis=2)
scan_info = np.delete(scan_info,duplicateIDs,axis=0)
#sorted_indices = scan_info.sort(key=get_slice_position, reverse=False)
sort_index = [i for i,x in sorted(enumerate(scan_info),key=get_slice_position, reverse=False)]
#scan_array = np.array(scan_array)
#scan_array = np.moveaxis(scan_array,[0,1,2],[2,0,1])
#scan_info = np.array(scan_info)
scan_info = scan_info[sort_index]
scan_array = scan_array[:,:,sort_index]
scan_info = scn_info.deduce_voxel_thickness(scan_info)
scan.scanInfo = scan_info
scan.scanArray = scan_array
scan.scanUID = "CT." + si_pixel_data[2]
return scan
[docs]
def parseScanInfoFromDB(scanObj, scanInfoList):
"""Assign scanInfo from the list of dictionaries to scanObj by matching zValue per slice
Args:
scanObj (cerr.dataclasses.scan.Scan): pyCERR's Scan object whose scanInfo needs to be populated
scanInfoList (list): list of dictionaries with fields corresponding to scanInfo
Returns:
0 when field assignment is successful
"""
# Create a new dictionary with DICOM names for keys
scanDir = os.path.dirname(__file__)
mappingFile = os.path.join(scanDir, 'dcm_cerr_name_map.json')
with open(mappingFile, 'r') as nameMapFile:
dcmCerrNameMap = json.load(nameMapFile)
cerrDcmNameMap = {v: k for k, v in dcmCerrNameMap.items()}
zValsScanV = [s.zValue for s in scanObj.scanInfo]
# Get z-values for scanInfoList and assign scanInfo
numberOfFrames = len(scanInfoList) #ds.NumberOfFrames.real
for iFrame in range(numberOfFrames):
imageOri = np.array(scanInfoList[iFrame]['ImageOrientationPatient'])
imagePos = np.array(scanInfoList[iFrame]['ImagePositionPatient'])
slice_normal = imageOri[[1,2,0]] * imageOri[[5,3,4]] \
- imageOri[[2,0,1]] * imageOri[[4,5,3]]
zValue = - np.sum(slice_normal * imagePos) / 10
indSlc = np.argmin((zValsScanV - zValue)**2)
fieldNames = scanInfoList[iFrame].keys()
for fieldName in fieldNames:
if fieldName in cerrDcmNameMap:
cerrField = cerrDcmNameMap[fieldName]
if hasattr(scanObj.scanInfo[indSlc], cerrField):
setattr(scanObj.scanInfo[indSlc], cerrField, scanInfoList[iFrame][fieldName])
scanObj.scanInfo[indSlc].imageOrientationPatient = imageOri
scanObj.scanInfo[indSlc].imagePositionPatient = imagePos
scanObj.scanInfo[indSlc].grid1Units = scanInfoList[iFrame]['PixelSpacing'][1]
scanObj.scanInfo[indSlc].grid2Units = scanInfoList[iFrame]['PixelSpacing'][0]
return 0
[docs]
def getScanNumFromUID(assocScanUID,planC) -> int:
"""
Args:
assocScanUID (str): UID of scan.
planC (cerr.plan_container.planC): pyCERR's plan container object.
Returns:
int: index within planC.scan that matches input assocScanUID.
"""
uid_list = [s.scanUID for s in planC.scan]
if assocScanUID in uid_list:
return uid_list.index(assocScanUID)
else:
return None
[docs]
def getCERRScanArrayFromITK(itkImage, assocScanNum, planC):
""" This routine returns a numpy array in pyCERR coordinate system (orientation) from a SimpleITK Image.
Args:
itkImage (SimpleITK.Image): SimpleITK's Image object
assocScanNum (int): Scan index to associate orientation of itkImage in pyCERR.
planC (cerr.planC_container.planC): pyCERR's plan container object.
Returns:
np.ndarray: array in CERR virtual coordinates.
"""
if isinstance(itkImage, sitk.Image):
itkImage = sitk.GetArrayFromImage(itkImage)
cerrArray = np.transpose(itkImage, (1, 2, 0))
# flip slices in CERR z-slice order which increases from head to toe
if flipSliceOrderFlag(planC.scan[assocScanNum]):
cerrArray = np.flip(cerrArray, axis=2)
return cerrArray