parseDicom.py 11 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341
  1. import os
  2. import sys
  3. import pydicom
  4. import numpy as np
  5. import re
  6. import slicer
  7. from slicer.ScriptedLoadableModule import *
  8. import pathlib
  9. #rom os import listdir
  10. #from os.path import isfile, join
  11. #onlyfiles = [f for f in listdir(mypath) if isfile(join(mypath, f))]
  12. #import Tkinter as tk
  13. #from Tkinter import filedialog
  14. #root = tk.Tk()
  15. #root.withdraw()
  16. #file_path = filedialog.askopenfilename()
  17. class parseDicom(ScriptedLoadableModule):
  18. def __init__(self, parent):
  19. ScriptedLoadableModule.__init__(self, parent)
  20. parent.title = "parseDicom"
  21. parent.categories = ["dynamicSPECT"]
  22. parent.dependencies = []
  23. parent.contributors = ["Andrej Studen (FMF/JSI)"] # replace with "Firstname Lastname (Org)"
  24. parent.helpText = """
  25. Parse dynamic SPECT DICOM files
  26. """
  27. parent.acknowledgementText = """
  28. This module was developed within the frame of the ARRS sponsored medical
  29. physics research programe to investigate quantitative measurements of cardiac
  30. function using sestamibi-like tracers
  31. """ # replace with organization, grant and thanks.
  32. self.parent = parent
  33. class parseDicomWidget(ScriptedLoadableModuleWidget):
  34. def setup(self):
  35. ScriptedLoadableModuleWidget.setup(self)
  36. self.logic=parseDicomLogic(self)
  37. class parseDicomLogic(ScriptedLoadableModuleLogic):
  38. def __init__(self,parent):
  39. ScriptedLoadableModuleLogic.__init__(self, parent)
  40. def setFileBrowser(self,fb):
  41. self.fb=fb
  42. def setTempBase(self,tb):
  43. self.tempBase=tb
  44. def filelist(self,mypath,remote=True):
  45. if remote:
  46. print("Using labkey")
  47. #not sure if labkey is available, so try it
  48. ok, files=self.fb.listRemoteDir(mypath)
  49. if not ok:
  50. print("Error accessing path")
  51. return []
  52. #files=[self.net.GetRelativePathFromLabkeyPath(f) for f in files]
  53. else:
  54. print("Using local files")
  55. #localPath=re.sub('file://','',mypath)
  56. localPath=mypath
  57. files = [os.path.join(localPath,f) for f in os.listdir(localPath)
  58. if os.path.isfile(os.path.join(localPath, f))]
  59. return files
  60. def getfile(self,f,remote=True):
  61. if remote:
  62. try:
  63. #not sure if labkey is available, but try it
  64. print("Using labkey")
  65. p=pathlib.Path(f)
  66. localPath=os.path.join(self.tempBase,p.name)
  67. self.fb.readFileToFile(f,localPath)
  68. return [open(localPath,'rb'),1]
  69. except:
  70. print('Could not access labkey. Exiting')
  71. return ['NULL',0]
  72. else:
  73. print("Using local directory")
  74. return [open(f,'rb'),1]
  75. return ['NULL',0]
  76. def readMasterFile(self,g):
  77. #this is the "master" file where data on other files can be had
  78. #here we found out the duration of the frame and their distribution through
  79. #phases and cycles
  80. try:
  81. plan = pydicom.dcmread(g)
  82. except:
  83. print ("{}: Not a dicom file".format(g))
  84. return False
  85. try:
  86. self.nframe=plan[0x0019,0x10a5].value;
  87. except:
  88. print ("{}: Not a master file".format(g))
  89. return False
  90. if not (type(self.nframe) is list) :
  91. print("nframe not a list")
  92. return False
  93. #nframe now holds for index i total number of frames collected up
  94. #to the end of each phase
  95. for i in range(1,len(self.nframe)):
  96. self.nframe[i]+=self.nframe[i-1]
  97. self.frame_start=plan[0x0019,0x10a7].value
  98. self.frame_stop=plan[0x0019,0x10a8].value
  99. self.frame_duration=plan[0x0019,0x10a9].value
  100. self.frame_time=np.zeros(self.nframe[-1]);
  101. self.frame_data=np.empty([1,1,1,self.nframe[-1]])
  102. self.center = [0,0,0]
  103. self.pixel_size =[0,0,0]
  104. self.frame_orientation=[0,0,0,0,0,0]
  105. return True
  106. def readNMFile(self,g):
  107. try:
  108. plan = pydicom.dcmread(g)
  109. except:
  110. print ("{}: Not a dicom file".format(g))
  111. return False
  112. try:
  113. pf=plan[0x0018,0x5020]
  114. except:
  115. print("Not a NM file. Exiting")
  116. return False
  117. try:
  118. phase=plan[0x0035,0x1005].value
  119. cycle=plan[0x0035,0x1004].value
  120. except:
  121. print("Missing phase/cycle values")
  122. return False
  123. #convert phase/cycle to frame index
  124. off=0
  125. if phase > 1:
  126. off=self.nframe[phase-2]
  127. ifi=off+cycle-1
  128. #from values in the master file determine frame time
  129. #(as the mid point between starting and ending the frame)
  130. self.frame_time[ifi]=0.5*(self.frame_start[ifi]+self.frame_stop[ifi]); #in ms
  131. print("({},{}) converted to {} at {} for {}".format(\
  132. phase,cycle,ifi,self.frame_time[ifi],self.frame_duration[ifi]))
  133. #play with pixel data
  134. if self.frame_data.shape[0] == 1:
  135. sh=np.transpose(plan.pixel_array,self.axisShift).shape;
  136. sh=list(sh)
  137. sh.append(self.nframe[-1])#add number of time slots
  138. self.frame_data=np.empty(sh)
  139. print(" Setting frame_data to {}".format(sh))
  140. #check & update pixel size
  141. pixel_size_read=[plan.PixelSpacing[0],plan.PixelSpacing[1],
  142. plan.SliceThickness]
  143. for i in range(0,3):
  144. if self.pixel_size[i] == 0:
  145. self.pixel_size[i] = float(pixel_size_read[i])
  146. if abs(self.pixel_size[i]-pixel_size_read[i]) > 1e-3:
  147. print('Pixel size mismatch {.2f}/{.2f}'.format(self.pixel_size[i],\
  148. pixel_size_read[i]))
  149. center_read=plan.DetectorInformationSequence[0].ImagePositionPatient
  150. print("Stored center at ({0},{1},{2})".format(self.center[0],self.center[1],self.center[2]))
  151. print("Read center at ({0},{1},{2})".format(center_read[0],center_read[1],center_read[2]))
  152. for i in range(0,3):
  153. if self.center[i] == 0:
  154. self.center[i] = float(center_read[i])
  155. if abs(self.center[i]-center_read[i]) > 1e-3:
  156. print('Image center mismatch {.2f}/{.2f}'.format(self.center[i],\
  157. center_read[i]))
  158. frame_orientation_read=plan.DetectorInformationSequence[0].ImageOrientationPatient
  159. for i in range(0,6):
  160. if self.frame_orientation[i] == 0:
  161. self.frame_orientation[i] = float(frame_orientation_read[i])
  162. if abs(self.frame_orientation[i]-frame_orientation_read[i]) > 1e-3:
  163. print('Image orientation mismatch {.2f}/{.2f}'.format(
  164. self.frame_rotation[i], frame_orientation_read[i]))
  165. self.frame_data[:,:,:,ifi]=np.transpose(plan.pixel_array,self.axisShift)
  166. return True
  167. def readCTFile(self,g):
  168. try:
  169. plan = pydicom.dcmread(g)
  170. except:
  171. print ("{}: Not a dicom file".format(g))
  172. return False
  173. if plan.Modality != 'CT':
  174. print ('{}: Not a CT file'.format(g))
  175. return False
  176. #this doesn't work in 2019 data version
  177. #if re.match("AC",plan.SeriesDescription) == None:
  178. # print (plan.SeriesDescription)
  179. # print ('Not a AC file')
  180. # continue
  181. try:
  182. iType=plan.ImageType
  183. except:
  184. print("Image type not found")
  185. return False
  186. if iType[3].find("SPI")<0:
  187. print("Not a spiral image")
  188. return False
  189. self.ct_data.append(\
  190. pydicom.pixel_data_handlers.util.apply_modality_lut(\
  191. plan.pixel_array,plan))
  192. self.ct_idx.append(plan.InstanceNumber)
  193. self.ct_z.append(plan.ImagePositionPatient[2])
  194. pixel_size_read=[plan.PixelSpacing[0],plan.PixelSpacing[1],
  195. plan.SliceThickness]
  196. for i in range(0,3):
  197. if self.ct_pixel_size[i] == 0:
  198. self.ct_pixel_size[i] = float(pixel_size_read[i])
  199. if abs(self.ct_pixel_size[i]-pixel_size_read[i]) > 1e-3:
  200. print('Pixel size mismatch {.2f}/{.2f}'.format(self.ct_pixel_size[i],
  201. pixel_size_read[i]))
  202. for i in range(0,2):
  203. if self.ct_center[i] == 0:
  204. self.ct_center[i] = float(plan.ImagePositionPatient[i])
  205. if abs(self.ct_center[i]-plan.ImagePositionPatient[i]) > 1e-3:
  206. print('Image center mismatch {.2f}/{.2f}'.format(self.ct_center[i],
  207. plan.ImagePositionPatient[i]))
  208. #not average, but minimum (!) why??
  209. if plan.ImagePositionPatient[2]<self.ct_center[2]:
  210. self.ct_center[2]=plan.ImagePositionPatient[2]
  211. for i in range(0,6):
  212. if self.ct_orientation[i] == 0:
  213. self.ct_orientation[i] = float(plan.ImageOrientationPatient[i])
  214. if abs(self.ct_orientation[i]-plan.ImageOrientationPatient[i]) > 1e-3:
  215. print('Image orientation mismatch {0:.2f}/{1:.2f}'.format(self.ct_orientation[i],\
  216. plan.ImageOrientationPatient[i]))
  217. return True
  218. def readMasterDirectory(self,mypath,remote=True):
  219. self.axisShift=(2,1,0)
  220. print("Reading master from {}".format(mypath))
  221. filelist=self.filelist(mypath,remote)
  222. for f in filelist:
  223. print('{}:'.format(f))
  224. g,ok=self.getfile(f,remote)
  225. if not(ok):
  226. return
  227. if self.readMasterFile(g):
  228. break
  229. def readNMDirectory(self,mypath,remote=True):
  230. files=self.filelist(mypath,remote)
  231. for f in files:
  232. g,ok=self.getfile(f,remote)
  233. if not(ok):
  234. continue
  235. self.readNMFile(g)
  236. return [self.frame_data,self.frame_time,self.frame_duration,self.center,
  237. self.pixel_size,self.frame_orientation]
  238. def readCTDirectory(self,mypath,remote=True):
  239. onlyfiles=self.filelist(mypath,remote)
  240. self.ct_data = []
  241. self.ct_idx = []
  242. self.ct_z = []
  243. self.ct_pixel_size = [0,0,0]
  244. self.ct_center = [0,0,0]
  245. self.ct_center[2]=1e30
  246. self.ct_orientation=[0,0,0,0,0,0]
  247. for f in onlyfiles:
  248. print('{}:'.format(f))
  249. g,ok=self.getfile(f,remote)
  250. if not(ok):
  251. return
  252. self.readCTFile(g)
  253. nz=len(self.ct_idx)
  254. #not average, again
  255. #ct_center[2]/=nz
  256. sh=self.ct_data[-1].shape
  257. sh_list=list(sh)
  258. sh_list.append(nz)
  259. data_array=np.zeros(sh_list)
  260. for k in range(0,nz):
  261. kp=int(np.round((self.ct_z[k]-self.ct_center[2])/self.ct_pixel_size[2]))
  262. data_array[:,:,kp]=np.transpose(self.ct_data[k])
  263. return data_array,self.ct_center,self.ct_pixel_size,self.ct_orientation