preprocess.py 7.2 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254
  1. import os
  2. import json
  3. import re
  4. import subprocess
  5. import nibabel
  6. import shutil
  7. import sys
  8. shome=os.path.expanduser('~nixUser')
  9. fhome=os.path.expanduser('~')
  10. with open(os.path.join(fhome,".labkey","setup.json")) as f:
  11. setup=json.load(f)
  12. sys.path.insert(0,setup["paths"]["labkeyInterface"])
  13. import labkeyInterface
  14. import labkeyDatabaseBrowser
  15. import labkeyFileBrowser
  16. sys.path.insert(0,setup["paths"]["orthancInterface"])
  17. import orthancInterface
  18. import orthancFileBrowser
  19. #sys.path.insert(1,shome+'/software/src/IPNUMM/dicomUtils')
  20. #import loadDicom
  21. fconfig=os.path.join(fhome,'.labkey','network.json')
  22. #matlab=os.path.join("/","data","software","install","matlab","bin","matlab")
  23. matlab=setup["paths"]["matlab"]
  24. #os.path.join(fhome,"software","install","matlab","bin","matlab")
  25. generalCodes=setup["paths"]["generalCodes"]
  26. #ios.path.join(fhome,"software","src","generalCodes")
  27. niftiTools=setup["paths"]["niftiTools"]
  28. #niftiTools=os.path.join(fhome,"software","src","NifTiScripts")
  29. net=labkeyInterface.labkeyInterface()
  30. net.init(fconfig)
  31. db=labkeyDatabaseBrowser.labkeyDB(net)
  32. fb=labkeyFileBrowser.labkeyFileBrowser(net)
  33. onet=orthancInterface.orthancInterface()
  34. onet.init(fconfig)
  35. ofb=orthancFileBrowser.orthancFileBrowser(onet)
  36. hi=0
  37. project='iPNUMMretro/Study'
  38. #for prospective, set
  39. project='IPNUMMprospektiva/Study'
  40. dataset='Imaging1'
  41. tempBase=os.path.join(fhome,'temp')
  42. participantField='PatientId'
  43. #for prospective set
  44. participantField='ParticipantId'
  45. #all images from database
  46. ds=db.selectRows(project,'study',dataset,[])
  47. #imageSelector={"CT":"CT","PET":"PETWB"};
  48. imageSelector={"CT":"CT_orthancId","PET":"PETWB_orthancId"}
  49. #output
  50. imageResampledField={"CT":"ctResampled","PET":"petResampled","patientmask":"ROImask"}
  51. #projectNIfTIBase=os.path.join(labkeyBase,'files',project,'@files/nifti')
  52. #use webdav to transfer file (even though it is localhost)
  53. def getPatientLabel(row,participantField='PatientId'):
  54. return row[participantField].replace('/','_')
  55. def getVisitLabel(row):
  56. return 'VISIT_'+str(int(row['SequenceNum']))
  57. def getStudyLabel(row,participantField='PatientId'):
  58. return getPatientLabel(row,participantField)+'-'+getVisitLabel(row)
  59. def runPreprocess_DM(matlab,generalCodes,niftiTools,studyDir):
  60. #run after all directories have been assembled
  61. script="addpath('"+generalCodes+"');"
  62. script+="addpath('"+niftiTools+"');"
  63. script+="preprocess_DM('"+studyDir+"',0,0);"
  64. script+="test;"
  65. script+="quit();"
  66. #outText=subprocess.check_output(["/bin/echo",script])
  67. try:
  68. outText=subprocess.check_output([matlab,"-nojvm","-r",script])
  69. except subprocess.CalledProcessError as e:
  70. print("Failed with:\n{}".format(e.output.decode('utf-8')))
  71. return False
  72. print(outText.decode('utf-8'))
  73. return True
  74. def getDicom(ofb,row,zipDir,rawDir,im,imageSelector,\
  75. participantField='PatientId'):
  76. #Load the dicom zip file and unzips it. If zip file is already at the expected path, it skips the loading step
  77. #Return True for valid outcome and False for troubles in row formating or unzip failures
  78. seriesId=row[imageSelector[im]];
  79. if seriesId=="0":
  80. return False
  81. print("{}: {}".format(im,seriesId))
  82. fname=os.path.join(zipDir,\
  83. getStudyLabel(row,participantField)+'_'+im+".zip");
  84. #copy data from orthanc
  85. if os.path.isfile(fname):
  86. print("Data already loaded. Skipping")
  87. else:
  88. print("Loading data from orthanc")
  89. ofb.getZip('series',seriesId,fname)
  90. #unzip the zipped dicom series
  91. unzipDir=os.path.join(rawDir,im)
  92. if os.path.isdir(unzipDir):
  93. print("Data already unzipped")
  94. return True
  95. try:
  96. os.mkdir(unzipDir)
  97. except FileExistsError:
  98. shutil.rmtree(unzipDir)
  99. try:
  100. outTxt=subprocess.check_output(["unzip","-d",unzipDir,"-xj",fname])
  101. except subprocess.CalledProcessError:
  102. print("unzip failed for {}".format(fname))
  103. return False
  104. return True
  105. def updateRow(project,dataset,row,imageResampledField,gzFileNames,\
  106. participantField='PatientId'):
  107. row['patientCode']=getPatientLabel(row,participantField)
  108. row['visitCode']=getVisitLabel(row)
  109. for im in imageResampledField:
  110. row[imageResampledField[im]]=gzFileNames[im]
  111. db.modifyRows('update',project,'study',dataset,[row])
  112. i=0
  113. for row in ds["rows"]:
  114. #interesting files are processedDir/studyName_CT_notCropped_2mmVoxel.nii
  115. #asn processedDir/studyName_PET_notCropped_2mmVoxel.nii
  116. volumeFileNames={im:\
  117. getStudyLabel(row,participantField)+'_'+im+
  118. '_notCropped_2mmVoxel.nii'\
  119. for im in imageResampledField}
  120. gzFileNames={im:f+".gz" \
  121. for (im,f) in volumeFileNames.items()}
  122. #build/check remote directory structure
  123. remoteDir=fb.buildPathURL(project,['preprocessedImages',\
  124. getPatientLabel(row,participantField),getVisitLabel(row)])
  125. gzRemoteFiles={im:remoteDir+'/'+f\
  126. for (im,f) in gzFileNames.items()}
  127. remoteFilePresent=[fb.entryExists(f)\
  128. for f in gzRemoteFiles.values()]
  129. for f in gzRemoteFiles.values():
  130. print("[{}]: [{}]".format(f,fb.entryExists(f)))
  131. if all(remoteFilePresent):
  132. print("Entry for row done.")
  133. updateRow(project,dataset,row,imageResampledField,\
  134. gzFileNames,participantField)
  135. continue
  136. #setup the directory structure for preprocess_DM
  137. studyDir=os.path.join(tempBase,getStudyLabel(row,participantField))
  138. if not os.path.isdir(studyDir):
  139. os.mkdir(studyDir)
  140. rawDir=os.path.join(studyDir,'Raw')
  141. if not os.path.isdir(rawDir):
  142. os.mkdir(rawDir)
  143. zipDir=os.path.join(studyDir,'Zip')
  144. if not os.path.isdir(zipDir):
  145. os.mkdir(zipDir)
  146. processedDir=os.path.join(studyDir,'Processed')
  147. if not os.path.isdir(processedDir):
  148. os.mkdir(processedDir)
  149. #specify local file names with path
  150. volumeFiles={im:os.path.join(processedDir,f)\
  151. for (im,f) in volumeFileNames.items()}
  152. gzFiles={im:f+".gz"\
  153. for (im,f) in volumeFiles.items()}
  154. filesPresent=[os.path.isfile(f) for f in gzFiles.values()]
  155. if not all(filesPresent):
  156. #use imageSelector -> inputs
  157. for im in imageSelector:
  158. #checks if raw files are already loaded
  159. getDicom(ofb,row,zipDir,rawDir,im,imageSelector,\
  160. participantField)
  161. #preprocess and zip
  162. ok=runPreprocess_DM(matlab,generalCodes,niftiTools,studyDir)
  163. if not ok:
  164. shutil.rmtree(studyDir)
  165. continue
  166. for f in volumeFiles.values():
  167. print("Running gzip {}".format(f))
  168. outText=subprocess.check_output(["/bin/gzip",f])
  169. print(outText.decode('utf-8'))
  170. #upload local files to remote
  171. for im in gzFiles:
  172. #for local,remote in zip(gzFiles,gzRemoteFiles):
  173. local=gzFiles[im]
  174. remote=gzRemoteFiles[im]
  175. print("Uploading {}".format(local))
  176. fb.writeFileToFile(local,remote)
  177. #update row and let it know where the processed files are
  178. updateRow(project,dataset,row,imageResampledField,gzFileNames,\
  179. participantField)
  180. #cleanup
  181. shutil.rmtree(studyDir)
  182. if i==-1:
  183. break
  184. i=i+1
  185. print("Done")