preprocess.py 7.5 KB

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