Transform voxel coordinates to world coordinates of DICOM files

import numpy as np
import SimpleITK as sitk
import csv
import os
import linecache
import inspect

seriesuid_label = ‘seriesuid’
coordX_label = ‘coordX’
coordY_label = ‘coordY’
coordZ_label = ‘coordZ’
CADProbability_label = ‘probability’
SpacingDict = {}

import fnmatch

def getDicomDir(dir, uid, resultList, dataFormat=’mhd’):
for root, dirnames, filenames in os.walk(dir):
if len(filenames) != 0:
for filename in fnmatch.filter(filenames, ‘*.’ + dataFormat):
if dataFormat.lower() == ‘mhd’:
if (filename == uid + ‘.mhd’):
resultList.append(os.path.join(dir, filename))
break
elif dataFormat.lower() == ‘dcm’:
if (os.path.basename(dir) == uid):
resultList.append(dir)
break
for dirname in dirnames:
getDicomDir(os.path.join(dir, dirname), uid, resultList, dataFormat)

def readCSV(filename):
lines = []
with open(filename, “rb”) as f:
csvreader = csv.reader(f)
for line in csvreader:
lines.append(line)
return lines

def getSpacingInfo(filepath, dataFormat=’mhd’):
if dataFormat.lower() == ‘mhd’:
inputimage = sitk.ReadImage(filepath)
elif dataFormat.lower() == ‘dcm’:
reader = sitk.ImageSeriesReader()
filenames = reader.GetGDCMSeriesFileNames(filepath)
reader.SetFileNames(filenames)
inputimage = reader.Execute()
else:
print “Unsupport file type: %s. Code line %d in function getSpacingInfo().” % (dataFormat, inspect.stack()[0][2])
exit()
spacing = inputimage.GetSpacing()
origin = inputimage.GetOrigin()
direction = inputimage.GetDirection()
return spacing + origin + direction

def mhdReader(inputpath):
dicomPath = ‘/home/kong/4T/LIDC-IDRI/LIDC-IDRI-0038/1.3.6.1.4.1.14519.5.2.1.6279.6001.127335250316932247803669656678/1.3.6.1.4.1.14519.5.2.1.6279.6001.117383608379722740629083782428’
reader = sitk.ImageSeriesReader()
filenames = reader.GetGDCMSeriesFileNames(dicomPath)
reader.SetFileNames(filenames)
inputimage = reader.Execute()
spacing = inputimage.GetSpacing()
origin = inputimage.GetOrigin()
direction = inputimage.GetDirection()
voxelCoord = [1, 1, 1]
result = voxelCoordToWorld(voxelCoord, spacing, origin, direction)
return spacing, origin, direction

def voxelCoordToWorld(voxelCoord, outputSpacing, origin, direction):
direcArray = np.array(direction).reshape((3, 3))
spacingArray = np.diag(outputSpacing)
voxelArray = np.array(voxelCoord).reshape((3, 1))
orginArray = np.array(origin).reshape((3, 1))
result = np.matmul(np.matmul(direcArray, spacingArray), voxelArray) + orginArray
return result.reshape(3)

def worldToVoxelCoord(worldCoord, origin, outputSpacing):
stretchedVoxelCoord = np.absolute(worldCoord - origin)
voxelCoord = stretchedVoxelCoord / outputSpacing
return voxelCoord

def initialSpacingDict(spacing_filename):
spacingLines = readCSV(spacing_filename)[1:]
for spacingLine in spacingLines:
SpacingDict[spacingLine[0]] = np.array(spacingLine[1:], dtype=’float’)

def transformFile_PixelToWorld(csvfile, ImgDir, spacing_filename, outfile, dataFormat=’mhd’):
initialSpacingDict(spacing_filename)
results = readCSV(csvfile)
header = results[0]
i = 1
with open(outfile, “wb”) as f:
csvwriter = csv.writer(f, delimiter=’,’)
csvwriter.writerow(header)
for result in results[1:]:
nodule_seriesuid = result[header.index(seriesuid_label)]
if i % 10 == 0:
print “progress: %d/%d” % (i, len(results) - 1)
i += 1
if not SpacingDict.has_key(nodule_seriesuid):
resultList = []
getDicomDir(ImgDir, nodule_seriesuid, resultList, dataFormat=dataFormat)
if (len(resultList) > 1):
print “more than one dicom folds was found. “ \
“The first one will be used. “ \
“All the found folds are showed as follows:”
print resultList
if (len(resultList) == 0):
print ‘No dicom file of case %s is found’ % nodule_seriesuid
print ‘This line will be ignored. Line data: ‘ + str(result)
continue
resultList = resultList[0]
SpacingDict[nodule_seriesuid] = getSpacingInfo(resultList, dataFormat)

  # SpacingDict\[nodule_seriesuid\] = \[spacing,origin,direction\]
coordX = float(result\[header.index(coordX_label)\])
coordY = float(result\[header.index(coordY_label)\])
coordZ = float(result\[header.index(coordZ_label)\])
worldCoord = voxelCoordToWorld(\[coordX, coordY, coordZ\],  SpacingDict\[nodule_seriesuid\]\[
                               0:3\], SpacingDict\[nodule\_seriesuid\]\[3:6\],  SpacingDict\[nodule\_seriesuid\]\[6:\])
csvwriter.writerow(\[nodule_seriesuid, str(worldCoord\[0\]), str(
    worldCoord\[1\]), str(worldCoord\[2\]), result\[header.index(CADProbability_label)\]\])

def startTransform():
csvfile = ‘/home/kong/4T/evaluationScript/results_own.csv’
dicomDir = ‘/home/kong/4T/LIDC-IDRI’
outfile = ‘/home/kong/4T/evaluationScript/results_own_world.csv’
spacing_filename = ‘/home/kong/4T/evaluationScript/annotations/spacing.csv’
transformFile_PixelToWorld(csvfile, dicomDir, spacing_filename, outfile, ‘mhd’)

if __name__ == ‘__main__‘:
startTransform()