dynamicSPECT/slicerScripts/convertToNRRD.py

import os
import sys
import json
import pathlib
import slicer
import subprocess
import shutil

#for all rows in xconfig[queryName] import dicom and convert to nrrd

def main(configFile):
   print('Running with {}'.format(configFile))
   with open(configFile) as f:
      xconfig=json.load(f)

   fsetup=os.path.join(os.path.expanduser('~'),'.labkey','setup.json')
   with open(fsetup) as f:
      setup=json.load(f)

   sys.path.append(setup['paths']['nixWrapper'])
   import nixWrapper

   nixWrapper.loadLibrary('labkeyInterface')
   import labkeyInterface
   import labkeyDatabaseBrowser
   import labkeyFileBrowser
   

   nixWrapper.loadLibrary('orthancInterface')
   import orthancInterface
   import orthancDatabaseBrowser
   import orthancFileBrowser
   
#one should be smart to figure this out if pwd is not the same as the directory of the script

   pwd=os.path.dirname(os.path.abspath(__file__))
   pwdUp=os.path.dirname(pwd)
   pythonScripts=os.path.join(pwdUp,'pythonScripts')
   sys.path.append(pythonScripts)
   import config
   import getData
   print('Config loaded')

   net=labkeyInterface.labkeyInterface()
   fnet=os.path.join(os.path.expanduser('~'),'.labkey',xconfig['network'])
   net.init(fnet)
   net.getCSRF()
   fb=labkeyFileBrowser.labkeyFileBrowser(net)
   db=labkeyDatabaseBrowser.labkeyDB(net)


   onet=orthancInterface.orthancInterface()
   onet.init(fnet)
   ofb=orthancFileBrowser.orthancFileBrowser(onet)
   odb=orthancDatabaseBrowser.orthancDB(onet)

   qFilter=config.getFilter(xconfig)
   ds=db.selectRows(xconfig['project'],xconfig['schemaName'],xconfig['queryName'],qFilter)


   try:
      rows=ds['rows']
   except KeyError:
      print('No rows returned')
      return
   for r in rows:
      #check row sanity
      

      print("Loading {}/{}".format(config.getPatientId(r,xconfig),config.getVisitId(r,xconfig)))
      rPath=fb.formatPathURL(xconfig['project'],config.getOutputDir(r,xconfig))
      rPath+='/'+config.getNodeName(r,xconfig,'CT')+'.nrrd'
      entryDone=fb.entryExists(rPath)
      if entryDone:
         try:
            if not xconfig['recalculate']:
               print('Entry done')
               getData.updateStatus(db,r,xconfig,'convertToNRRD')
               continue
         except KeyError:
            print('Entry done')
            getData.updateStatus(db,r,xconfig,'convertToNRRD')
            continue

         print('Forced recalculation')

      code=config.getCode(r,xconfig)
      print('{} available:{}'.format(config.getNodeName(r,xconfig,'CT'),fb.entryExists(rPath)))
      #loadPatient into slicer
      patient=loadPatient(ofb,r,xconfig)
      if not patient:
         print(f'{code} failed to load from Orthanc')
         getData.updateStatus(db,r,xconfig,'convertToNRRD',status=False)
         continue
      #convert to nodes
      addCT(r,patient,xconfig)
      addFrames(r,patient,xconfig)
      addDummyInputFunction(r,patient,xconfig)
      
      nodes=slicer.mrmlScene.GetNodesByClass('vtkMRMLScalarVolumeNode')
      print('Nodes')
      for node in nodes:
         print('\t{}'.format(node.GetName()))
         storeNode(fb,r,xconfig,node)

      nodes=slicer.mrmlScene.GetNodesByClass('vtkMRMLDoubleArrayNode')
      print('Nodes (double array)')
      for node in nodes:
         print('\t{}'.format(node.GetName()))
         storeNode(fb,r,xconfig,node)

      nodes=slicer.mrmlScene.GetNodesByClass('vtkMRMLTableNode')
      print('Nodes (table)')
      for node in nodes:
         print('\t{}'.format(node.GetName()))
         storeNode(fb,r,xconfig,node)



      clearNodes(r,xconfig)
      #addCT and addFrames fill r['ct'] and r['spect']
      db.modifyRows('update',xconfig['project'],xconfig['schemaName'],xconfig['queryName'],[r])
      getData.updateStatus(db,r,xconfig,'convertToNRRD')

      
def clearNodes(row,xconfig):
   nodes=slicer.mrmlScene.GetNodesByClass('vtkMRMLScalarVolumeNode')
   nodes1=slicer.mrmlScene.GetNodesByClass('vtkMRMLDoubleArrayNode')
   nodes2=slicer.mrmlScene.GetNodesByClass('vtkMRMLTableNode')
   for n in nodes1:
      nodes.AddItem(n)
   for n in nodes2:
      nodes.AddItem(n)

   res=[slicer.mrmlScene.RemoveNode(f) for f in nodes]

def loadPatient(ofb,r,xconfig):
   sys.path.append('../pythonScripts')
   import parseDicom
   import vtkInterface

   pd=parseDicom.parseDicom()

   masterPath=downloadAndUnzip(ofb,r,"nmMaster",xconfig)
   if not masterPath:
      return None
   pd.readMasterDirectory(masterPath,False)
   print('Time [{} .. {}]'.format(pd.frame_start,pd.frame_stop))
   clearUnzipDir(r,xconfig)

   nmPath=downloadAndUnzip(ofb,r,"nmCorrData",xconfig)
   if not nmPath:
      return None
   frame_data, frame_time, frame_duration, frame_origin, \
       frame_pixel_size, frame_orientation=\
       pd.readNMDirectory(nmPath,False)
   print('Frame time {}'.format(frame_time))
   clearUnzipDir(r,xconfig)

   ctPath=downloadAndUnzip(ofb,r,"ct",xconfig)
   if not ctPath:
      return None
   ct_data,ct_origin,ct_pixel_size, \
      ct_orientation=pd.readCTDirectory(ctPath,False)
   print('CT pixel {}'.format(ct_pixel_size))
   clearUnzipDir(r,xconfig)

   ct_orientation=vtkInterface.completeOrientation(ct_orientation)
   frame_orientation=vtkInterface.completeOrientation(frame_orientation)

   ct={'data':ct_data,'origin':ct_origin,'pixel_size':ct_pixel_size,
      'orientation':ct_orientation}
   nm={'data':frame_data,'time':frame_time,'duration':frame_duration,
      'origin':frame_origin,'pixel_size':frame_pixel_size,
      'orientation':frame_orientation}

   return {'CT':ct,'NM':nm}
   
   print('Done')

def clearUnzipDir(r,xconfig):
   import config

   zipDir=config.getLocalDir(r,xconfig)
   try:
      os.mkdir(zipDir)
   except FileExistsError:
      shutil.rmtree(zipDir)
   return zipDir

def downloadAndUnzip(ofb,r,code,xconfig):
   import config 
   
   pathList=xconfig['tempDir'].split('/')
   pathList.insert(0,os.path.expanduser('~'))
   tempDir=os.path.join(*pathList)
   if not os.path.isdir(tempDir):
      print('Creating {}'.format(tempDir))
      os.makedirs(tempDir)

   orthancId=r[code+'OrthancId']
   if not orthancId:
      return None

   fileCode='{}_{}'.format(config.getCode(r,xconfig),code)
   zipFile=os.path.join(tempDir,fileCode+'.zip')
   ofb.getZip('series',orthancId,zipFile)
   zipDir=clearUnzipDir(r,xconfig)

   try:
      outTxt=subprocess.check_output(["unzip","-d",zipDir,"-xj",zipFile])
   except subprocess.CalledProcessError:
      print("unzip failed for {}".format(zipFile))
      return ""

   return zipDir

def addCT(r,patient,xconfig):
   import config
   import vtkInterface
   ct=patient['CT']
   vtkData=vtkInterface.numpyToVTK(ct['data'],ct['data'].shape)
   nodeName=config.getNodeName(r,xconfig,'CT')
   addNode(nodeName,vtkData,ct['origin'],ct['pixel_size'],ct['orientation'],0)
   r['ct']=f'{nodeName}.nrrd'


def addFrames(r,patient,xconfig):
   import vtkInterface
   import config
   #convert data from numpy.array to vtkImageData
   #use time point it
   nm=patient['NM']
   print("NFrames: {}".format(nm['data'].shape[3]))
   for it in range(0,nm['data'].shape[3]):
      frame_data=nm['data'][:,:,:,it];
      nodeName=config.getNodeName(r,xconfig,'NM',it)
      vtkData=vtkInterface.numpyToVTK(frame_data,frame_data.shape)
      addNode(nodeName,vtkData,nm['origin'],nm['pixel_size'],nm['orientation'],1)
      #last one will be kept
      r['spect']=f'{nodeName}.nrrd'


def addNode(nodeName,v,origin,pixel_size,orientation,dataType):
   #origin,orientation in lps
   #dataType=0 is CT (to background)
   #dataType=1 is SPECT, view not adjusted, foreground, 
   newNode=slicer.vtkMRMLScalarVolumeNode()
   newNode.SetName(nodeName)
   v.SetOrigin([0,0,0])
   v.SetSpacing([1,1,1])
   ijkToRAS = vtk.vtkMatrix4x4()
   #think how to do this with image orientation
   #orientation from lps to ras
   rasOrientation=[-orientation[i] if (i%3 < 2) else orientation[i]
        for i in range(0,len(orientation))]
   #origin from lps to ras
   rasOrigin=[-origin[i] if (i%3<2) else origin[i] for i in range(0,len(origin))]

   for i in range(0,3):
      for j in range(0,3):
         ijkToRAS.SetElement(i,j,pixel_size[i]*rasOrientation[3*j+i])

      ijkToRAS.SetElement(i,3,rasOrigin[i])

   newNode.SetIJKToRASMatrix(ijkToRAS)
   newNode.SetAndObserveImageData(v)
   slicer.mrmlScene.AddNode(newNode)

def addDummyInputFunction(r,patient,xconfig):

   import config
   
   nm=patient['NM']
   n=nm['data'].shape[3]

   dnsNodeName=config.getNodeName(r,xconfig,'Dummy')
   dns = slicer.mrmlScene.GetNodesByClassByName('vtkMRMLDoubleArrayNode',dnsNodeName)
   if dns.GetNumberOfItems() == 0:
      try:
         dn = slicer.mrmlScene.AddNode(slicer.vtkMRMLDoubleArrayNode())
      except AttributeError:
         addDummyTable(dnsNodeName,n,nm)
         return
   else:
      dn = dns.GetItemAsObject(0)

   dn.SetSize(n)

   ft=nm['time']
   dt=nm['duration']
   for i in range(0,n):
      fx=ft[i]
      fy=dt[i]
      dn.SetValue(i, 0, fx)
      dn.SetValue(i, 1, fy)
      dn.SetValue(i, 2, 0)
      print('{} ({},{})'.format(i,fx,fy))



def addDummyTable(dnsNodeName,n,nm): 
   
   #add vtkMRMLTableNodes
   dns = slicer.mrmlScene.GetNodesByClassByName('vtkMRMLTableNode',dnsNodeName)
   if dns.GetNumberOfItems() == 0:
      dn = slicer.mrmlScene.AddNode(slicer.vtkMRMLTableNode())
      dn.SetName(dnsNodeName)
   else:
      dn = dns.GetItemAsObject(0)
      dn.RemoveAllColumns()
      

   ft=nm['time']
   dt=nm['duration']
   tCol=vtk.vtkDoubleArray()
   dCol=vtk.vtkDoubleArray()
   for i in range(n):
      tCol.InsertNextValue(ft[i])
      dCol.InsertNextValue(dt[i])

   tcol=dn.AddColumn(tCol)
   tcol.SetName('time')
   dcol=dn.AddColumn(dCol)
   dcol.SetName('duration')

def storeNode(fb,row,xconfig,node):
   import config

   suffix=".nrrd"
   isTable=False
   if node.__class__.__name__=="vtkMRMLDoubleArrayNode":
      suffix=".mcsv"
   if node.__class__.__name__=="vtkMRMLTableNode":
      suffix=".mcsv"
      isTable=True      
   if (node.__class__.__name__=="vtkMRMLTransformNode" or \
         node.__class__.__name__=="vtkMRMLGridTransformNode"):
      suffix=".h5"

   fileName=node.GetName()+suffix

   localPath=os.path.join(config.getLocalDir(row,xconfig),fileName)

   if isTable:
      storeTable(node,localPath)
   else:
      slicer.util.saveNode(node,localPath)
   print("Stored to: {}".format(localPath))
   labkeyPath=fb.buildPathURL(xconfig['project'],config.getPathList(row,xconfig))
   print ("Remote: {}".format(labkeyPath))
   #checks if exists
   remoteFile=labkeyPath+'/'+fileName
   fb.writeFileToFile(localPath,remoteFile)


def storeTable(node,filename):
   #mimic old vtkMRMLDoubleArray format
   table=node.GetTable()
   ft=table.GetColumnByName('time')
   fd=table.GetColumnByName('duration')
   n=ft.GetNumberOfValues()
   print(f'Storing {n} values')
   with open(filename,'w') as f:
      f.write(f'# measurement file {filename}\n')
      f.write('# no labels\n')
      for i in range(n):
         _t=ft.GetTuple1(i)
         _d=fd.GetTuple1(i)
         print(f'{_t},{_d},0')
         f.write(f'{_t},{_d},0\n')

      






def readPatient(fb,localDir,project,patientId):

    rDir=fb.formatPathURL(project,'/'.join([patientId]))
    lDir=os.path.join(localDir,patientId)
    if not os.path.isdir(lDir):
        os.makedirs(lDir)

    ok,files=fb.listRemoteDir(rDir)
    locFiles=[]
    for f in files:
        print(f)
        p=pathlib.Path(f)
        localFile=os.path.join(lDir,p.name)
        fb.readFileToFile(f,localFile)
        locFiles.append(localFile)
    print('Done')
    return locFiles

if __name__=='__main__':
    main(sys.argv[1])
    quit()