preprocess.py 7.1 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252
  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. imageResampledField={"CT":"ctResampled","PET":"petResampled"}
  50. #projectNIfTIBase=os.path.join(labkeyBase,'files',project,'@files/nifti')
  51. #use webdav to transfer file (even though it is localhost)
  52. def getPatientLabel(row,participantField='PatientId'):
  53. return row[participantField].replace('/','_')
  54. def getVisitLabel(row):
  55. return 'VISIT_'+str(int(row['SequenceNum']))
  56. def getStudyLabel(row,participantField='PatientId'):
  57. return getPatientLabel(row,participantField)+'-'+getVisitLabel(row)
  58. def runPreprocess_DM(matlab,generalCodes,niftiTools,studyDir):
  59. #run after all directories have been assembled
  60. script="addpath('"+generalCodes+"');"
  61. script+="addpath('"+niftiTools+"');"
  62. script+="preprocess_DM('"+studyDir+"',0,0);"
  63. script+="test;"
  64. script+="quit();"
  65. #outText=subprocess.check_output(["/bin/echo",script])
  66. try:
  67. outText=subprocess.check_output([matlab,"-nojvm","-r",script])
  68. except subprocess.CalledProcessError as e:
  69. print("Failed with:\n{}".format(e.output.decode('utf-8')))
  70. return False
  71. print(outText.decode('utf-8'))
  72. return True
  73. def getDicom(ofb,row,zipDir,rawDir,im,imageSelector,\
  74. participantField='PatientId'):
  75. #Load the dicom zip file and unzips it. If zip file is already at the expected path, it skips the loading step
  76. #Return True for valid outcome and False for troubles in row formating or unzip failures
  77. seriesId=row[imageSelector[im]];
  78. if seriesId=="0":
  79. return False
  80. print("{}: {}".format(im,seriesId))
  81. fname=os.path.join(zipDir,\
  82. getStudyLabel(row,participantField)+'_'+im+".zip");
  83. #copy data from orthanc
  84. if os.path.isfile(fname):
  85. print("Data already loaded. Skipping")
  86. else:
  87. print("Loading data from orthanc")
  88. ofb.getZip('series',seriesId,fname)
  89. #unzip the zipped dicom series
  90. unzipDir=os.path.join(rawDir,im)
  91. if os.path.isdir(unzipDir):
  92. print("Data already unzipped")
  93. return True
  94. try:
  95. os.mkdir(unzipDir)
  96. except FileExistsError:
  97. shutil.rmtree(unzipDir)
  98. try:
  99. outTxt=subprocess.check_output(["unzip","-d",unzipDir,"-xj",fname])
  100. except subprocess.CalledProcessError:
  101. print("unzip failed for {}".format(fname))
  102. return False
  103. return True
  104. def updateRow(project,dataset,row,imageResampledField,gzFileNames,\
  105. participantField='PatientId'):
  106. row['patientCode']=getPatientLabel(row,participantField)
  107. row['visitCode']=getVisitLabel(row)
  108. for im in imageResampledField:
  109. row[imageResampledField[im]]=gzFileNames[im]
  110. db.modifyRows('update',project,'study',dataset,[row])
  111. i=0
  112. for row in ds["rows"]:
  113. #interesting files are processedDir/studyName_CT_notCropped_2mmVoxel.nii
  114. #asn processedDir/studyName_PET_notCropped_2mmVoxel.nii
  115. volumeFileNames={im:\
  116. getStudyLabel(row,participantField)+'_'+im+
  117. '_notCropped_2mmVoxel.nii'\
  118. for im in imageSelector}
  119. gzFileNames={im:f+".gz" \
  120. for (im,f) in volumeFileNames.items()}
  121. #build/check remote directory structure
  122. remoteDir=fb.buildPathURL(project,['preprocessedImages',\
  123. getPatientLabel(row,participantField),getVisitLabel(row)])
  124. gzRemoteFiles={im:remoteDir+'/'+f\
  125. for (im,f) in gzFileNames.items()}
  126. remoteFilePresent=[fb.entryExists(f)\
  127. for f in gzRemoteFiles.values()]
  128. for f in gzRemoteFiles.values():
  129. print("[{}]: [{}]".format(f,fb.entryExists(f)))
  130. if all(remoteFilePresent):
  131. print("Entry for row done.")
  132. updateRow(project,dataset,row,imageResampledField,\
  133. gzFileNames,participantField)
  134. continue
  135. #setup the directory structure for preprocess_DM
  136. studyDir=os.path.join(tempBase,getStudyLabel(row,participantField))
  137. if not os.path.isdir(studyDir):
  138. os.mkdir(studyDir)
  139. rawDir=os.path.join(studyDir,'Raw')
  140. if not os.path.isdir(rawDir):
  141. os.mkdir(rawDir)
  142. zipDir=os.path.join(studyDir,'Zip')
  143. if not os.path.isdir(zipDir):
  144. os.mkdir(zipDir)
  145. processedDir=os.path.join(studyDir,'Processed')
  146. if not os.path.isdir(processedDir):
  147. os.mkdir(processedDir)
  148. #specify local file names with path
  149. volumeFiles={im:os.path.join(processedDir,f)\
  150. for (im,f) in volumeFileNames.items()}
  151. gzFiles={im:f+".gz"\
  152. for (im,f) in volumeFiles.items()}
  153. filesPresent=[os.path.isfile(f) for f in gzFiles.values()]
  154. if not all(filesPresent):
  155. for im in imageSelector:
  156. #checks if raw files are already loaded
  157. getDicom(ofb,row,zipDir,rawDir,im,imageSelector,\
  158. participantField)
  159. #preprocess and zip
  160. ok=runPreprocess_DM(matlab,generalCodes,niftiTools,studyDir)
  161. if not ok:
  162. shutil.rmtree(studyDir)
  163. continue
  164. for f in volumeFiles.values():
  165. print("Running gzip {}".format(f))
  166. outText=subprocess.check_output(["/bin/gzip",f])
  167. print(outText.decode('utf-8'))
  168. #upload local files to remote
  169. for im in gzFiles:
  170. #for local,remote in zip(gzFiles,gzRemoteFiles):
  171. local=gzFiles[im]
  172. remote=gzRemoteFiles[im]
  173. print("Uploading {}".format(local))
  174. fb.writeFileToFile(local,remote)
  175. #update row and let it know where the processed files are
  176. updateRow(project,dataset,row,imageResampledField,gzFileNames,\
  177. participantField)
  178. #cleanup
  179. shutil.rmtree(studyDir)
  180. if i==-1:
  181. break
  182. i=i+1
  183. print("Done")