dynamicSPECT/slicerScripts/convertToNRRD.py

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

def main(configFile):
   print('Running with {}'.format(configFile))
   with open(configFile) as f:
      config=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


   net=labkeyInterface.labkeyInterface()
   fnet=os.path.join(os.path.expanduser('~'),'.labkey',config['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=[]  
   ds=db.selectRows(config['project'],config['schemaName'],config['queryName'],qFilter)


   try:
      rows=ds['rows']
   except KeyError:
      print('No rows returned')
      return
   for r in rows:
      print("Loading {}/{}".format(patientId(r,config),visitId(r,config)))
      nNames=nodeNames(r,config)
      rPath=fb.formatPathURL(config['project'],outputDir(r,config))
      rPath+='/'+nNames['CT']+'.nrrd'
      entryDone=fb.entryExists(rPath)
      if entryDone:
         try:
            if not config['recalculate']:
               print('Entry done')
               continue
         except KeyError:
            print('Entry done')
            continue

         print('Forced recalculation')

            
      print('{} available:{}'.format(nNames['CT'],fb.entryExists(rPath)))
      #loadPatient into slicer
      patient=loadPatient(ofb,r,config)
      #convert to nodes
      addCT(r,patient,config)
      addFrames(r,patient,config)
      addDummyInputFunction(r,patient,config)
      
      nodes=slicer.mrmlScene.GetNodesByClass('vtkMRMLScalarVolumeNode')
      print('Nodes')
      for node in nodes:
         print('\t{}'.format(node.GetName()))
         storeNode(fb,r,config,node)

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

      
def patientId(row,config):
   return row[config['ParticipantField']]

def visitId(row,config):
   return row['visitName']

def code(row,config):
   return '{}_{}'.format(patientId(row,config),visitId(row,config))

def pathList(row,config):
   return [config['baseDir'],patientId(row,config),visitId(row,config)]

def outputDir(row,config):
   return '/'.join(pathList(row,config))

def localDir(row,config):
   return os.path.join(config['tempDir'],code(row,config))

def nodeNames(row,config):
   names={}
   x=code(row,config)
   names['NM']=['{}_Volume{}'.format(x,i) for i in range(20)]
   names['CT']='{}_CT'.format(x)
   return names

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

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

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

   pd=parseDicom.parseDicom()

   masterPath=downloadAndUnzip(ofb,r,"nmMaster",config)
   pd.readMasterDirectory(masterPath,False)
   print('Time [{} .. {}]'.format(pd.frame_start,pd.frame_stop))
   nmPath=downloadAndUnzip(ofb,r,"nmCorrData",config)
   frame_data, frame_time, frame_duration, frame_origin, \
       frame_pixel_size, frame_orientation=\
       pd.readNMDirectory(nmPath,False)
   print('Frame time {}'.format(frame_time))
   ctPath=downloadAndUnzip(ofb,r,"ct",config)
   ct_data,ct_origin,ct_pixel_size, \
      ct_orientation=pd.readCTDirectory(ctPath,False)
   print('CT pixel {}'.format(ct_pixel_size))

   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 downloadAndUnzip(ofb,r,xcode,config):
   orthancId=r[xcode+'OrthancId']
   fileCode='{}_{}'.format(code(r,config),xcode)
   zipFile=os.path.join(config['tempDir'],fileCode+'.zip')
   ofb.getZip('series',orthancId,zipFile)
   zipDir=localDir(r,config)
   try:
      os.mkdir(zipDir)
   except FileExistsError:
      shutil.rmtree(zipDir)
   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,config):
   sys.path.append('../pythonScripts')
   import vtkInterface
   ct=patient['CT']
   vtkData=vtkInterface.numpyToVTK(ct['data'],ct['data'].shape)
   nodeName=code(r,config)+'_CT'
   addNode(nodeName,vtkData,ct['origin'],ct['pixel_size'],ct['orientation'],0)

def addFrames(r,patient,config):
   sys.path.append('../pythonScripts')
   import vtkInterface
   #convert data from numpy.array to vtkImageData
   #use time point it
   x=code(r,config)
   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=x+'_Volume'+str(it)
      vtkData=vtkInterface.numpyToVTK(frame_data,frame_data.shape)
      addNode(nodeName,vtkData,nm['origin'],nm['pixel_size'],nm['orientation'],1)



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,config):

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

   dnsNodeName=code(r,config)+'_Dummy'
   dns = slicer.mrmlScene.GetNodesByClassByName('vtkMRMLDoubleArrayNode',dnsNodeName)
   if dns.GetNumberOfItems() == 0:
      dn = slicer.mrmlScene.AddNode(slicer.vtkMRMLDoubleArrayNode())
      dn.SetName(dnsNodeName)
   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 storeNode(fb,row,config,node):

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

   fileName=node.GetName()+suffix

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

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


   






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()