SlicerLabkeyExtension/DICOMtools/importDicom.py

import dicom
import vtkInterface as vi
import os
import vtk, qt, ctk, slicer
from slicer.ScriptedLoadableModule import *
import slicerNetwork
import json
import loadDicom
import DICOMLib

class importDicom(slicer.ScriptedLoadableModule.ScriptedLoadableModule):
    def __init__(self,parent):
        slicer.ScriptedLoadableModule.ScriptedLoadableModule.__init__(self, parent)
        self.className="importDICOM"
        self.parent.title="importDICOM"
        self.parent.categories = ["LabKey"]
        self.parent.dependencies = []
        self.parent.contributors = ["Andrej Studen (UL/FMF)"] # replace with "Firstname Lastname (Organization)"
        self.parent.helpText = """
                utilities for parsing dicom entries
                """
        self.parent.acknowledgementText = """
                Developed within the medical physics research programme of the Slovenian research agency.
                """ # replace with organization, grant and thanks.

class importDicomWidget(slicer.ScriptedLoadableModule.ScriptedLoadableModuleWidget):
  """Uses ScriptedLoadableModuleWidget base class, available at:
  https://github.com/Slicer/Slicer/blob/master/Base/Python/slicer/ScriptedLoadableModule.py
  """

  def setup(self):
      
    slicer.ScriptedLoadableModule.ScriptedLoadableModuleWidget.setup(self)
    self.logic=importDicomLogic(self)
    self.network=slicerNetwork.labkeyURIHandler()
    
    connectionCollapsibleButton = ctk.ctkCollapsibleButton()
    connectionCollapsibleButton.text = "Connection"
    self.layout.addWidget(connectionCollapsibleButton)

    connectionFormLayout = qt.QFormLayout(connectionCollapsibleButton)
    
    self.loadConfigButton=qt.QPushButton("Load configuration")
    self.loadConfigButton.toolTip="Load configuration"
    self.loadConfigButton.connect('clicked(bool)',self.onLoadConfigButtonClicked)
    connectionFormLayout.addRow("Connection:",self.loadConfigButton)

    self.DICOMDirectory=qt.QLineEdit("Test/Temp/%40files/TEST/MLEM")
    connectionFormLayout.addRow("LabKey directory:",self.DICOMDirectory)


    #loadDICOMButton=qt.QPushButton("Load")
    #loadDICOMButton.toolTip="Load DICOM"
    #loadDICOMButton.clicked.connect(self.onLoadDICOMButtonClicked)
    #connectionFormLayout.addRow("DICOM:",loadDICOMButton)

    self.DICOMFilter=qt.QLineEdit('{"seriesNumber":"SeriesLabel"}')
    connectionFormLayout.addRow("Filter(JSON):",self.DICOMFilter)

    loadDICOMFilterButton=qt.QPushButton("Load with filter")
    loadDICOMFilterButton.toolTip="Load DICOM with filter"
    loadDICOMFilterButton.clicked.connect(self.onLoadDICOMFilterButtonClicked)
    connectionFormLayout.addRow("DICOM:",loadDICOMFilterButton)

  def onLoadConfigButtonClicked(self):
       filename=qt.QFileDialog.getOpenFileName(None,'Open configuration file (JSON)',
            os.path.join(os.path.expanduser('~'),'.labkey'), '*.json')
       self.network.parseConfig(filename)
       self.network.initRemote()
       self.loadConfigButton.setText(os.path.basename(filename))
  
  def onLoadDICOMFilterButtonClicked(self):
      filter=json.loads(self.DICOMFilter.text)
      #print("Filter is {}".format(filter))
      self.logic.loadVolumes(self.network,self.DICOMDirectory.text,filter)
    
  def onLoadDICOMButtonClicked(self):
      self.logic.load(self.network,self.DICOMDirectory.text)
       


#equivalent of loadable + labkey interface
class dicomSeries():
    def __init__(self):
        self.data = []
        self.idx = []
        self.z = []
        self.pixel_size = [0,0,0]
        self.lpsOrigin = [0,0,0]
        self.lpsOrigin[2]=1e30
        self.lpsOrientation=[0,0,0,0,0,0]
        self.local=False


    def getfile(self,net,relativePath):
        if self.local:
                return open(relativePath,'rb')
        return net.readFileToBuffer(relativePath)

    def addFile(self,f):
        try:
            self.files.append(f)
        except:
            self.files=[f]

    def setLabel(self,label):
        self.label=label

    def getLabel(self):
        try:
            return self.label
        except:
            return None

    def setMetadata(self,key,value):
        try:
            self.metadata[key]=value
        except:
            self.metadata={key:value}

    def getMetadata(self):
        try:
            return self.metadata
        except:
            return {}

            

    def load(self,net):

        for f in self.files:
            print '{}:'.format(f)
            fileBuffer=self.getfile(net,f)
            self.loadFile(fileBuffer)

        nz=len(self.idx)
        sh=self.data[-1].shape
        sh_list=list(sh)
        sh_list.append(nz)
        data_array=np.zeros(sh_list)

        for k in range(0,nz):
            kp=int(np.round((self.z[k]-self.center[2])/self.pixel_size[2]))
            data_array[:,:,kp]=np.transpose(self.data[k])

        try:
            nodeName='Series'+self.label
        except:
            print('Could not set series label')
            nodeName='UnknownSeries'

        newNode=slicer.vtkMRMLScalarVolumeNode()
        newNode.SetName(nodeName)

        ijkToRAS = vtk.vtkMatrix4x4()
        #think how to do this with image orientation

        rasOrientation=[-self.lpsOrientation[i] if (i%3 < 2) else self.lpsOrientation[i]
          for i in range(0,len(self.lpsOrientation))]
        rasOrigin=[-self.lpsOrigin[i] if (i%3<2) else self.lpsOrigin[i] for i in range(0,len(self.lpsOrigin))]

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

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

        newNode.SetIJKToRASMatrix(ijkToRAS)

        v=vtk.vtkImageData()
        v.GetPointData().SetScalars(
            vtk.util.numpy_support.numpy_to_vtk(
                np.ravel(self.data,order='F'),deep=True, array_type=vtk.VTK_FLOAT))
        v.SetOrigin(0,0,0)
        v.SetSpacing(1,1,1)
        v.SetDimensions(self.data.shape)

        newNode.SetAndObserveImageData(v)
        slicer.mrmlScene.AddNode(newNode)
        volume={'node':newNode,'metadata':self.metadata}
        return volume

    def loadFile(self,fileBuffer):

        plan=dicom.read_file(fileBuffer)

        self.data.append(plan.pixel_array)
        self.idx.append(plan.InstanceNumber)
        self.z.append(plan.ImagePositionPatient[2])


        #pixelSize
        pixel_size=[plan.PixelSpacing[0],plan.PixelSpacing[1],
                    plan.SliceThickness]

        for i in range(0,3):
            if self.pixel_size[i] == 0:
                self.pixel_size[i] = float(pixel_size[i])
            if abs(self.pixel_size[i]-pixel_size[i]) > 1e-3:
                print 'Pixel size mismatch {.2f}/{.2f}'.format(self.pixel_size[i],
                        pixel_size[i])

        #origin
        for i in range(0,2):
            if self.lpsOrigin[i] == 0:
                self.lpsOrigin[i] = float(plan.ImagePositionPatient[i])
            if abs(self.lpsOrigin[i]-plan.ImagePositionPatient[i]) > 1e-3:
                print 'Image center mismatch {.2f}/{.2f}'.format(self.lpsOrigin[i],
                    plan.ImagePositionPatient[i])
        #not average, but minimum (!) why??

        if plan.ImagePositionPatient[2]<self.lpsOrigin[2]:
            self.lpsOrigin[2]=plan.ImagePositionPatient[2]


        #orientation
        for i in range(0,6):
            if self.lpsOrientation[i] == 0:
                self.lpsOrientation[i] = float(plan.ImageOrientationPatient[i])
            if abs(self.lpsOrientation[i]-plan.ImageOrientationPatient[i]) > 1e-3:
                print 'Image orientation mismatch {0:.2f}/{1:.2f}'.format(self.lpsOrientation[i],
                        plan.ImageOrientationPatient[i])

        return True

class importDicomLogic(slicer.ScriptedLoadableModule.ScriptedLoadableModuleLogic):
    def __init__(self,parent):
        slicer.ScriptedLoadableModule.ScriptedLoadableModuleLogic.__init__(self, parent)
        self.local=False
        self.tag={
        'studyInstanceUid':0x0020000d,
        'seriesInstanceUid':0x0020000e,
        'patientId':0x00100020,
        'patientName':0x00100010,
        'sequenceName':0x00180024,
        'seriesNumber':0x00200011,
        'percentPhaseFieldOfView':0x00180094,
        'modality': 0x00080060,
        'patientSex': 0x00100040,
        'patientBirthDate': 0x00100030,
        'patientComments': 0x00104000,
        'studyDescription': 0x00081030,
        'studyDate': 0x00080020,
        'studyId': 0x00200010,
        'studyTime': 0x00080030,
        'frameOfReferenceInstanceUid':0x00200052}
    


    def setLocal(self, basePath):
        self.local=True
        self.basePath=basePath

        
    def loadVolumes(self,net,directory,filter):
        #mimic examineForImport
        seriesList=self.examineForImport(net,directory,filter)
        print("Got {} series").format(len(seriesList))
        volumes=[]

        for s in seriesList:
            try:
                volumes.append(s.load(net))
                #often fails, e.g. JPEGLossles
            except:
                loadable=DICOMLib.DICOMLoadable()
                loadable.name='Series'+str(s.getLabel())
                print("Loading for {} number of files (pre-load) {}").format(loadable.name,len(s.files))
                loadable.files=[net.DownloadFileToCache(f) for f in s.files]
                print("Loading for {} number of files (pre-sort) {}").format(loadable.name,len(loadable.files))

                loadable.files,distances,loadable.warning=DICOMLib.DICOMUtils.getSortedImageFiles(loadable.files,1e-3)

                print("Loading for {} number of files {}").format(loadable.name,len(loadable.files))
                try:
                    volumeNode=self.volumePlugin.load(loadable)
                except:
                    self.volumePlugin=slicer.modules.dicomPlugins['DICOMScalarVolumePlugin']()
                    volumeNode=self.volumePlugin.load(loadable)
                volume={'node':volumeNode,'metadata':s.getMetadata()}
                volumes.append(volume)


        return volumes


    def listdir(self,net,relativeDirectory):
        if self.local:
            dirs=os.listdir(os.path.join(self.basePath,relativeDirectory))
            return [os.path.join(relativeDirectory,dir) for dir in dirs]
        return net.listRelativeDir(relativeDirectory)

    def getfile(self,net,relativePath):
        if self.local:
            return open(os.path.join(self.basePath,relativePath),'rb')
        return net.readFileToBuffer(relativePath)

    def examineForImport(self,net,directory,filter):
        #split by series
        seriesList=[]

        files=self.listdir(net,directory)
        if len(files)==0:
            print("No input found in {}".format(directory))
            return seriesList


        for f in files:
            fileBuffer=self.getfile(net,f)

            #validate
            try:
                plan = dicom.read_file(fileBuffer)
            except:
                print ("{}: Not a dicom file")
                continue

            #determine validity first
            fileValid=True
            for key in filter:
                if filter[key]==None:
                    continue
                if filter[key]=='SeriesLabel':
                    seriesTag=self.tag[key]
                    continue

                v=plan[self.tag[key]].value
                if not v==filter[key]:
                    print('Filter mismatch {}{:x}: {}/{}').format(key,self.tag[key],v,filter[key])
                    fileValid=False

            if not fileValid:
                continue

            #determine serieslabel second
            seriesLabel=plan[seriesTag].value
                
            try:
                if series.getLabel()==seriesLabel:
                    series.addFile(f)
                    continue
            except NameError:
                pass

            #add new series
            seriesList.append(dicomSeries())
            series=seriesList[-1]
            series.local=self.local

            #set series parameters
            series.addFile(f)
            series.setLabel(seriesLabel)
            for key in filter:
                if not filter[key]==None:
                    continue
                v=plan[self.tag[key]].value
                series.setMetadata(key,v)

        return seriesList