convertToNRRD.py 8.4 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298
  1. import os
  2. import sys
  3. import json
  4. import pathlib
  5. import slicer
  6. import subprocess
  7. import shutil
  8. def main(configFile):
  9. print('Running with {}'.format(configFile))
  10. with open(configFile) as f:
  11. xconfig=json.load(f)
  12. fsetup=os.path.join(os.path.expanduser('~'),'.labkey','setup.json')
  13. with open(fsetup) as f:
  14. setup=json.load(f)
  15. sys.path.append(setup['paths']['nixWrapper'])
  16. import nixWrapper
  17. nixWrapper.loadLibrary('labkeyInterface')
  18. import labkeyInterface
  19. import labkeyDatabaseBrowser
  20. import labkeyFileBrowser
  21. nixWrapper.loadLibrary('orthancInterface')
  22. import orthancInterface
  23. import orthancDatabaseBrowser
  24. import orthancFileBrowser
  25. sys.path.append('../pythonScripts')
  26. import config
  27. net=labkeyInterface.labkeyInterface()
  28. fnet=os.path.join(os.path.expanduser('~'),'.labkey',xconfig['network'])
  29. net.init(fnet)
  30. net.getCSRF()
  31. fb=labkeyFileBrowser.labkeyFileBrowser(net)
  32. db=labkeyDatabaseBrowser.labkeyDB(net)
  33. onet=orthancInterface.orthancInterface()
  34. onet.init(fnet)
  35. ofb=orthancFileBrowser.orthancFileBrowser(onet)
  36. odb=orthancDatabaseBrowser.orthancDB(onet)
  37. qFilter=[]
  38. ds=db.selectRows(xconfig['project'],xconfig['schemaName'],xconfig['queryName'],qFilter)
  39. try:
  40. rows=ds['rows']
  41. except KeyError:
  42. print('No rows returned')
  43. return
  44. for r in rows:
  45. print("Loading {}/{}".format(config.getPatientId(r,xconfig),config.getVisitId(r,xconfig)))
  46. rPath=fb.formatPathURL(xconfig['project'],config.getOutputDir(r,xconfig))
  47. rPath+='/'+config.getNodeName(r,xconfig,'CT')+'.nrrd'
  48. entryDone=fb.entryExists(rPath)
  49. if entryDone:
  50. try:
  51. if not xconfig['recalculate']:
  52. print('Entry done')
  53. continue
  54. except KeyError:
  55. print('Entry done')
  56. continue
  57. print('Forced recalculation')
  58. print('{} available:{}'.format(config.getNodeName(r,xconfig,'CT'),fb.entryExists(rPath)))
  59. #loadPatient into slicer
  60. patient=loadPatient(ofb,r,xconfig)
  61. #convert to nodes
  62. addCT(r,patient,xconfig)
  63. addFrames(r,patient,xconfig)
  64. addDummyInputFunction(r,patient,xconfig)
  65. nodes=slicer.mrmlScene.GetNodesByClass('vtkMRMLScalarVolumeNode')
  66. print('Nodes')
  67. for node in nodes:
  68. print('\t{}'.format(node.GetName()))
  69. storeNode(fb,r,xconfig,node)
  70. nodes=slicer.mrmlScene.GetNodesByClass('vtkMRMLDoubleArrayNode')
  71. print('Nodes (double array)')
  72. for node in nodes:
  73. print('\t{}'.format(node.GetName()))
  74. storeNode(fb,r,xconfig,node)
  75. clearNodes(r,xconfig)
  76. def clearNodes(row,xconfig):
  77. nodes=slicer.mrmlScene.GetNodesByClass('vtkMRMLScalarVolumeNode')
  78. nodes1=slicer.mrmlScene.GetNodesByClass('vtkMRMLDoubleArrayNode')
  79. for n in nodes1:
  80. nodes.AddItem(n)
  81. res=[slicer.mrmlScene.RemoveNode(f) for f in nodes]
  82. def loadPatient(ofb,r,xconfig):
  83. sys.path.append('../pythonScripts')
  84. import parseDicom
  85. import vtkInterface
  86. pd=parseDicom.parseDicom()
  87. masterPath=downloadAndUnzip(ofb,r,"nmMaster",xconfig)
  88. pd.readMasterDirectory(masterPath,False)
  89. print('Time [{} .. {}]'.format(pd.frame_start,pd.frame_stop))
  90. clearUnzipDir(r,xconfig)
  91. nmPath=downloadAndUnzip(ofb,r,"nmCorrData",xconfig)
  92. frame_data, frame_time, frame_duration, frame_origin, \
  93. frame_pixel_size, frame_orientation=\
  94. pd.readNMDirectory(nmPath,False)
  95. print('Frame time {}'.format(frame_time))
  96. clearUnzipDir(r,xconfig)
  97. ctPath=downloadAndUnzip(ofb,r,"ct",xconfig)
  98. ct_data,ct_origin,ct_pixel_size, \
  99. ct_orientation=pd.readCTDirectory(ctPath,False)
  100. print('CT pixel {}'.format(ct_pixel_size))
  101. clearUnzipDir(r,xconfig)
  102. ct_orientation=vtkInterface.completeOrientation(ct_orientation)
  103. frame_orientation=vtkInterface.completeOrientation(frame_orientation)
  104. ct={'data':ct_data,'origin':ct_origin,'pixel_size':ct_pixel_size,
  105. 'orientation':ct_orientation}
  106. nm={'data':frame_data,'time':frame_time,'duration':frame_duration,
  107. 'origin':frame_origin,'pixel_size':frame_pixel_size,
  108. 'orientation':frame_orientation}
  109. return {'CT':ct,'NM':nm}
  110. print('Done')
  111. def clearUnzipDir(r,xconfig):
  112. import config
  113. zipDir=config.getLocalDir(r,xconfig)
  114. try:
  115. os.mkdir(zipDir)
  116. except FileExistsError:
  117. shutil.rmtree(zipDir)
  118. return zipDir
  119. def downloadAndUnzip(ofb,r,code,xconfig):
  120. import config
  121. orthancId=r[code+'OrthancId']
  122. fileCode='{}_{}'.format(config.getCode(r,xconfig),code)
  123. zipFile=os.path.join(xconfig['tempDir'],fileCode+'.zip')
  124. ofb.getZip('series',orthancId,zipFile)
  125. zipDir=clearUnzipDir(r,xconfig)
  126. try:
  127. outTxt=subprocess.check_output(["unzip","-d",zipDir,"-xj",zipFile])
  128. except subprocess.CalledProcessError:
  129. print("unzip failed for {}".format(zipFile))
  130. return ""
  131. return zipDir
  132. def addCT(r,patient,xconfig):
  133. sys.path.append('../pythonScripts')
  134. import config
  135. import vtkInterface
  136. ct=patient['CT']
  137. vtkData=vtkInterface.numpyToVTK(ct['data'],ct['data'].shape)
  138. nodeName=config.getNodeName(r,xconfig,'CT')
  139. addNode(nodeName,vtkData,ct['origin'],ct['pixel_size'],ct['orientation'],0)
  140. def addFrames(r,patient,xconfig):
  141. import vtkInterface
  142. import config
  143. #convert data from numpy.array to vtkImageData
  144. #use time point it
  145. nm=patient['NM']
  146. print("NFrames: {}".format(nm['data'].shape[3]))
  147. for it in range(0,nm['data'].shape[3]):
  148. frame_data=nm['data'][:,:,:,it];
  149. nodeName=config.getNodeName(r,xconfig,'NM',it)
  150. vtkData=vtkInterface.numpyToVTK(frame_data,frame_data.shape)
  151. addNode(nodeName,vtkData,nm['origin'],nm['pixel_size'],nm['orientation'],1)
  152. def addNode(nodeName,v,origin,pixel_size,orientation,dataType):
  153. #origin,orientation in lps
  154. #dataType=0 is CT (to background)
  155. #dataType=1 is SPECT, view not adjusted, foreground,
  156. newNode=slicer.vtkMRMLScalarVolumeNode()
  157. newNode.SetName(nodeName)
  158. v.SetOrigin([0,0,0])
  159. v.SetSpacing([1,1,1])
  160. ijkToRAS = vtk.vtkMatrix4x4()
  161. #think how to do this with image orientation
  162. #orientation from lps to ras
  163. rasOrientation=[-orientation[i] if (i%3 < 2) else orientation[i]
  164. for i in range(0,len(orientation))]
  165. #origin from lps to ras
  166. rasOrigin=[-origin[i] if (i%3<2) else origin[i] for i in range(0,len(origin))]
  167. for i in range(0,3):
  168. for j in range(0,3):
  169. ijkToRAS.SetElement(i,j,pixel_size[i]*rasOrientation[3*j+i])
  170. ijkToRAS.SetElement(i,3,rasOrigin[i])
  171. newNode.SetIJKToRASMatrix(ijkToRAS)
  172. newNode.SetAndObserveImageData(v)
  173. slicer.mrmlScene.AddNode(newNode)
  174. def addDummyInputFunction(r,patient,xconfig):
  175. import config
  176. nm=patient['NM']
  177. n=nm['data'].shape[3]
  178. dnsNodeName=config.getNodeName(r,xconfig,'Dummy')
  179. dns = slicer.mrmlScene.GetNodesByClassByName('vtkMRMLDoubleArrayNode',dnsNodeName)
  180. if dns.GetNumberOfItems() == 0:
  181. dn = slicer.mrmlScene.AddNode(slicer.vtkMRMLDoubleArrayNode())
  182. dn.SetName(dnsNodeName)
  183. else:
  184. dn = dns.GetItemAsObject(0)
  185. dn.SetSize(n)
  186. ft=nm['time']
  187. dt=nm['duration']
  188. for i in range(0,n):
  189. fx=ft[i]
  190. fy=dt[i]
  191. dn.SetValue(i, 0, fx)
  192. dn.SetValue(i, 1, fy)
  193. dn.SetValue(i, 2, 0)
  194. print('{} ({},{})'.format(i,fx,fy))
  195. def storeNode(fb,row,xconfig,node):
  196. import config
  197. suffix=".nrrd"
  198. if node.__class__.__name__=="vtkMRMLDoubleArrayNode":
  199. suffix=".mcsv"
  200. if (node.__class__.__name__=="vtkMRMLTransformNode" or \
  201. node.__class__.__name__=="vtkMRMLGridTransformNode"):
  202. suffix=".h5"
  203. fileName=node.GetName()+suffix
  204. localPath=os.path.join(config.getLocalDir(row,xconfig),fileName)
  205. slicer.util.saveNode(node,localPath)
  206. print("Stored to: {}".format(localPath))
  207. labkeyPath=fb.buildPathURL(xconfig['project'],config.getPathList(row,xconfig))
  208. print ("Remote: {}".format(labkeyPath))
  209. #checks if exists
  210. remoteFile=labkeyPath+'/'+fileName
  211. fb.writeFileToFile(localPath,remoteFile)
  212. def readPatient(fb,localDir,project,patientId):
  213. rDir=fb.formatPathURL(project,'/'.join([patientId]))
  214. lDir=os.path.join(localDir,patientId)
  215. if not os.path.isdir(lDir):
  216. os.makedirs(lDir)
  217. ok,files=fb.listRemoteDir(rDir)
  218. locFiles=[]
  219. for f in files:
  220. print(f)
  221. p=pathlib.Path(f)
  222. localFile=os.path.join(lDir,p.name)
  223. fb.readFileToFile(f,localFile)
  224. locFiles.append(localFile)
  225. print('Done')
  226. return locFiles
  227. if __name__=='__main__':
  228. main(sys.argv[1])
  229. quit()