Operational version

cardiacSPECT/cardiacSPECT.py

@@ -94,6 +94,25 @@ class cardiacSPECTWidget(ScriptedLoadableModuleWidget):
 
     self.dataLoadButton = dataLoadButton
 
+    self.patientId=qt.QLineEdit();
+    dataFormLayout.addRow('Patient ID', self.patientId)
+
+    patientLoadButton = qt.QPushButton("Load")
+    patientLoadButton.toolTip="Load data from DICOM"
+    dataFormLayout.addRow("Load Patient",patientLoadButton)
+    patientLoadButton.clicked.connect(self.onPatientLoadButtonClicked)
+
+    saveVolumeButton = qt.QPushButton("Save Volume to NRRD")
+    saveVolumeButton.toolTip="Save to NRRD"
+    dataFormLayout.addRow("Save to NRRD",saveVolumeButton)
+    saveVolumeButton.clicked.connect(self.onSaveVolumeButtonClicked)
+
+    transformNodeButton = qt.QPushButton("Transform Node")
+    transformNodeButton.toolTip="Transform node with patient based transform"
+    dataFormLayout.addRow("Transform Nodes",transformNodeButton)
+    transformNodeButton.clicked.connect(self.onTransformNodeButtonClicked)
+
+
     # Add vertical spacer
     self.layout.addStretch(1)
 
@@ -235,7 +254,8 @@ class cardiacSPECTWidget(ScriptedLoadableModuleWidget):
        it=self.time_frame_select.value
        selectionNode = slicer.app.applicationLogic().GetSelectionNode()
        print("Propagating CT volume")
-       node=slicer.mrmlScene.GetFirstNodeByName("testCT")
+       nodeName=self.patientId.text+'CT'
+       node=slicer.mrmlScene.GetFirstNodeByName(nodeName)
        selectionNode.SetReferenceActiveVolumeID(node.GetID())
        if self.resetPosition==1:
           self.resetPosition=0
@@ -243,7 +263,7 @@ class cardiacSPECTWidget(ScriptedLoadableModuleWidget):
        else:
           slicer.app.applicationLogic().PropagateVolumeSelection(0)
        print("Propagating SPECT volume")
-       nodeName='testVolume'+str(it)
+       nodeName=self.patientId.text+'Volume'+str(it)
        node=slicer.mrmlScene.GetFirstNodeByName(nodeName)
        selectionNode.SetSecondaryVolumeID(node.GetID())
        slicer.app.applicationLogic().PropagateForegroundVolumeSelection(0)
@@ -278,13 +298,13 @@ class cardiacSPECTWidget(ScriptedLoadableModuleWidget):
         for j in range(0,ns):
             #add node for data
             dn = slicer.mrmlScene.AddNode(slicer.vtkMRMLDoubleArrayNode())
-            dn.SetName(self.logic.getSegmentName(j))
+            dn.SetName(self.patientId.text+'_'+self.logic.getSegmentName(j))
             a = dn.GetArray()
             a.SetNumberOfTuples(n)
             dt=0;
             t0=0;
             for i in range(0,n):
-                vol="testVolume"+str(i)
+                vol=self.patientId.text+"Volume"+str(i)
                 fx=ft[i]
                 fy=self.logic.meanROI(vol,j)
                 dt=2*ft[i]-t0
@@ -314,7 +334,19 @@ class cardiacSPECTWidget(ScriptedLoadableModuleWidget):
   def onCountSegmentsClicked(self):
         self.countSegmentsDisplay.setText(self.logic.countSegments())
 
-  #def onAddFrameButtonClicked(self):
+  def onPatientLoadButtonClicked(self):
+      self.logic.loadPatient(self,self.patientId.text)
+
+  def onSaveVolumeButtonClicked(self):
+      self.logic.storeVolumeNodes(self.patientId.text,
+            self.time_frame_select.minimum,self.time_frame_select.maximum)
+
+  def onTransformNodeButtonClicked(self):
+      self.logic.applyTransform(self.patientId.text,
+          self.time_frame_select.minimum,self.time_frame_select.maximum)
+
+
+#def onAddFrameButtonClicked(self):
 #      it=int(self.time_frame_select.text)
 #      self.logic.addFrame(it)
 
@@ -345,17 +377,18 @@ class cardiacSPECTLogic(ScriptedLoadableModuleLogic):
   def loadData(self,widget):
 
     inputDir=str(widget.dataPath.text)
+    self.pd.readMasterDirectory(inputDir)
     self.frame_data, self.frame_time, self.frame_origin, \
-        self.frame_pixel_size, self.frame_orientation=self.pd.read_dynamic_SPECT(inputDir)
+        self.frame_pixel_size, self.frame_orientation=self.pd.readNMDirectory(inputDir)
 
     self.ct_data,self.ct_origin,self.ct_pixel_size, \
-        self.ct_orientation=self.pd.read_CT(inputDir)
+        self.ct_orientation=self.pd.readCTDirectory(inputDir)
 
     self.ct_orientation=vi.completeOrientation(self.ct_orientation)
     self.frame_orientation=vi.completeOrientation(self.frame_orientation)
 
-    self.addCT()
-    self.addFrames()
+    self.addCT('test')
+    self.addFrames('test')
 
     widget.time_frame_select.setMaximum(self.frame_data.shape[3]-1)
 
@@ -365,6 +398,47 @@ class cardiacSPECTLogic(ScriptedLoadableModuleLogic):
         'Slicer Python','Data loaded')
 
 
+  def loadPatient(self,widget,patientId):
+      print("Loading {}").format(patientId)
+      ds=self.net.loadDataset("dinamic_spect/Patients","Imaging")
+      for r in ds['rows']:
+        print r
+        if r['aliasID']==patientId:
+            visit=r
+
+      print visit
+      dicoms=self.net.loadDataset("Test/Transfer","Imaging")
+      for r in dicoms['rows']:
+          if not r['PatientId']==visit['aliasID']:
+              continue
+          if abs(r['SequenceNum']-float(visit['nmMaster']))<0.1:
+              path=r['_labkeyurl_Series']
+              path=path[:path.rfind('/')]
+              masterPath="labkey://"+path
+          if abs(r['SequenceNum']-float(visit['nmCorrData']))<0.1:
+              path=r['_labkeyurl_Series']
+              path=path[:path.rfind('/')]
+              nmPath="labkey://"+path
+          if abs(r['SequenceNum']-float(visit['ctData']))<0.1:
+              path=r['_labkeyurl_Series']
+              path=path[:path.rfind('/')]
+              ctPath="labkey://"+path
+
+
+      self.pd.readMasterDirectory(masterPath)
+      self.frame_data, self.frame_time, self.frame_origin, \
+          self.frame_pixel_size, self.frame_orientation=self.pd.readNMDirectory(nmPath)
+      self.ct_data,self.ct_origin,self.ct_pixel_size, \
+         self.ct_orientation=self.pd.readCTDirectory(ctPath)
+
+      self.ct_orientation=vi.completeOrientation(self.ct_orientation)
+      self.frame_orientation=vi.completeOrientation(self.frame_orientation)
+
+      self.addCT(patientId)
+      self.addFrames(patientId)
+
+      widget.time_frame_select.setMaximum(self.frame_data.shape[3]-1)
+
   def addNode(self,nodeName,v, lpsOrigin, pixel_size, lpsOrientation, dataType):
 
        # if dataType=0 it is CT data, which gets propagated to background an is
@@ -400,21 +474,21 @@ class cardiacSPECTLogic(ScriptedLoadableModuleLogic):
        newNode.SetAndObserveImageData(v)
        slicer.mrmlScene.AddNode(newNode)
 
-  def addFrames(self):
+  def addFrames(self,patientId):
        #convert data from numpy.array to vtkImageData
        #use time point it
        print "NFrames: {}".format(self.frame_data.shape[3])
        for it in range(0,self.frame_data.shape[3]):
            frame_data=self.frame_data[:,:,:,it];
-           nodeName='testVolume'+str(it)
+           nodeName=patientId+'Volume'+str(it)
            self.addNode(nodeName,
               vi.numpyToVTK(frame_data,frame_data.shape),
               self.frame_origin,
               self.frame_pixel_size,
               self.frame_orientation,1)
 
-  def addCT(self):
-       nodeName='testCT'
+  def addCT(self,patientId):
+       nodeName=patientId+'CT'
        self.addNode(nodeName,
             #vi.numpyToVTK3D(self.ct_data,
             #    self.ct_origin,self.ct_pixel_size),
@@ -422,6 +496,26 @@ class cardiacSPECTLogic(ScriptedLoadableModuleLogic):
             self.ct_origin,self.ct_pixel_size,
             self.ct_orientation,0)
 
+  def rFromI(i,volumeNode):
+      ijkToRas = vtk.vtkMatrix4x4()
+      volumeNode.GetIJKToRASMatrix(ijkToRas)
+      vImage=volumeNode.GetImageData()
+      i1=list(vImage.GetPoint(i))
+      i1=i1.append(1)
+      #ras are global coordinates (in mm)
+      position_ras=ijkToRas.MultiplyPoint(i1)
+      return position_ras[0:3]
+
+  def IfromR(pos,volumeNode):
+      fM=vtk.vtkMatrix4x4()
+      volumeNode.GetRASToIJKMatrix(fM)
+      fM.MultiplyPoint(pos)
+      vImage=volumeNode.GetImageData()
+      #nearest neighbor
+      return vImage.FindPoint(pos[0:3])
+
+
+
   def meanROI(self, volName1, i):
     s=0
 
@@ -460,7 +554,7 @@ class cardiacSPECTLogic(ScriptedLoadableModuleLogic):
     extent=mask.GetExtent()
     fM=vtk.vtkMatrix4x4()
     mask.GetWorldToImageMatrix(fM)
-
+    ns=0
     for i in range(0,n):
 
       #get global coordinates of point i
@@ -487,8 +581,9 @@ class cardiacSPECTLogic(ScriptedLoadableModuleLogic):
 
       #use maskValue to project ROI data
       s+=maskValue*dataImage.GetPointData().GetScalars().GetTuple1(i)
+      ns+=1
 
-    return s/n
+    return s/ns
 
   def countSegments(self):
     segNodeList=slicer.mrmlScene.GetNodesByClass("vtkMRMLSegmentationNode")
@@ -510,6 +605,90 @@ class cardiacSPECTLogic(ScriptedLoadableModuleLogic):
           return "NONE"
       return segNode.GetSegmentation().GetSegment(segNode.GetSegmentation().GetNthSegmentID(i)).GetName()
 
+  def storeNodeRemote(self,relativePath,nodeName):
+      labkeyPath=self.pd.net.GetLabkeyWebdavUrl()+'/'+relativePath
+      print ("Remote: {}").format(labkeyPath)
+      #checks if exists
+      self.pd.net.mkdir(labkeyPath)
+
+      localPath=self.pd.net.GetLocalCacheDirectory()+'/'+relativePath
+
+      node=slicer.mrmlScene.GetFirstNodeByName(nodeName)
+      if node==None:
+          print("Node {} not found").format(nodeName)
+          return
+
+      suffix=".nrrd"
+      if node.__class__.__name__=="vtkMRMLDoubleArrayNode":
+          suffix=".mcsv"
+      if node.__class__.__name__=="vtkMRMLTransformNode":
+          suffix=".h5"
+      fileName=nodeName+suffix
+      file=localPath+'/'+fileName
+      slicer.util.saveNode(node,file)
+      print("Stored to: {}").format(file)
+      f=open(file)
+      remoteFile=labkeyPath+'/'+fileName
+      self.pd.net.put(remoteFile,f.read())
+
+
+  def storeVolumeNodes(self,patientId,n1,n2):
+      #n1=self.time_frame.minimum;
+      #n2=self.time_frame.maximum
+      project="dinamic_spect/Patients"
+      relativePath=project+'/@files/'+patientId
+
+      print("Store CT")
+      nodeName=patientId+'CT'
+      self.storeNodeRemote(relativePath,nodeName)
+
+      print("Storing NM from {} to {}").format(n1,n2)
+      n=n2-n1+1
+      for i in range(n):
+          it=i+n1
+          nodeName=patientId+'Volume'+str(it)
+          self.storeNodeRemote(relativePath,nodeName)
+
+      segNodeName="Heart"
+      self.storeNodeRemote(relativePath,segNodeName)
+
+      doubleArrayNodeName=patientId+'_Ventricle'
+      self.storeNodeRemote(relativePath,doubleArrayNodeName)
+
+      transformNodeName=patientId+"_DF"
+      self.storeNodeRemote(relativePath,transformNodeName)
+
+  def applyTransform(self, patientId,n1,n2):
+      transformNodeName=patientId+"_DF"
+      transformNode=slicer.util.getFirstNodeByName(transformNodeName)
+      if transformNode==None:
+          print("Node {} not found").format(transformNodeName)
+          return
+      try:
+          if self.transformApplied:
+              print("Transform already applied")
+              return
+      except AttributeError:
+          print("Not defined yet")
+
+      self.transformApplied=True
+
+      n=n2-n1+1
+      for i in range(n):
+          it=i+n1
+          nodeName=patientId+'Volume'+str(it)
+          node=slicer.util.getFirstNodeByName(nodeName)
+          if node==None:
+              continue
+          node.ApplyTransform(transformNode.GetTransformToParent())
+
+      nodeName=patientId+'CT'
+      node=slicer.util.getFirstNodeByName(nodeName)
+      if not node==None:
+          node.ApplyTransform(transformNode.GetTransformToParent())
+
+
+
 
 class cardiacSPECTTest(ScriptedLoadableModuleTest):
   """

cardiacSPECT/parseDicom.py

@@ -83,155 +83,234 @@ class parseDicomLogic(ScriptedLoadableModuleLogic):
 
         return ['NULL',0]
 
-    def read_dynamic_SPECT(self,mypath):
-        axisShift=(2,1,0)
+    def readMasterFile(self,g):
+        #this is the "master" file where data on other files can be had
+        #here we found out the duration of the frame and their distribution through
+        #phases and cycles
+        try:
+            plan = dicom.read_file(g)
+        except:
+            print ("{}: Not a dicom file").format(g)
+            return False
+
+        try:
+            self.nframe=plan[0x0019,0x10a5].value;
+        except:
+            print ("{}: Not a master file").format(g)
+            return False
+        if not (type(self.nframe) is list) :
+            print("nframe not a list")
+            return False
+
+        #nframe now holds for index i total number of frames collected up
+        #to the end of each phase
+
+        for i in range(1,len(self.nframe)):
+            self.nframe[i]+=self.nframe[i-1]
+
+        self.frame_start=plan[0x0019,0x10a7].value
+        self.frame_stop=plan[0x0019,0x10a8].value
+        self.frame_duration=plan[0x0019,0x10a9].value
+
+        self.frame_time=np.zeros(self.nframe[-1]);
+        self.frame_data=np.empty([1,1,1,self.nframe[-1]])
+        self.center = [0,0,0]
+        self.pixel_size =[0,0,0]
+        self.frame_orientation=[0,0,0,0,0,0]
+        return True
+
+    def readNMFile(self,g):
+
+        try:
+            plan = dicom.read_file(g)
+        except:
+            print ("{}: Not a dicom file").format(g)
+            return False
+
+        try:
+            pf=plan[0x0018,0x5020]
+        except:
+            print("Not a NM file. Exiting")
+            return False
+
+        try:
+            phase=plan[0x0035,0x1005].value
+            cycle=plan[0x0035,0x1004].value
+        except:
+            print("Missing phase/cycle values")
+            return False
+
+        #convert phase/cycle to frame index
+        off=0
+        if phase > 1:
+            off=self.nframe[phase-2]
+        ifi=off+cycle-1
+
+        #from values in the master file determine frame time
+        #(as the mid point between starting and ending the frame)
+        self.frame_time[ifi]=0.5*(self.frame_start[ifi]+self.frame_stop[ifi]); #in ms
+
+        print "({},{}) converted to {} at {} for {}".format(
+            phase,cycle,ifi,self.frame_time[ifi],self.frame_duration[ifi])
+
+
+#play with pixel data
+        if self.frame_data.shape[0] == 1:
+            sh=np.transpose(plan.pixel_array,self.axisShift).shape;
+            sh=list(sh)
+            sh.append(self.nframe[-1])#add number of time slots
+            self.frame_data=np.empty(sh)
+            print "Setting frame_data to",sh
+
+        #check & update pixel size
+        pixel_size_read=[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_read[i])
+            if abs(self.pixel_size[i]-pixel_size_read[i]) > 1e-3:
+                print 'Pixel size mismatch {.2f}/{.2f}'.format(self.pixel_size[i],
+                pixel_size_read[i])
+
+        center_read=plan.DetectorInformationSequence[0].ImagePositionPatient
+        print "Stored center at ({0},{1},{2})".format(self.center[0],self.center[1],self.center[2])
+        print "Read   center at ({0},{1},{2})".format(center_read[0],center_read[1],center_read[2])
+        for i in range(0,3):
+            if self.center[i] == 0:
+                self.center[i] = float(center_read[i])
+            if abs(self.center[i]-center_read[i]) > 1e-3:
+                print 'Image center mismatch {.2f}/{.2f}'.format(self.center[i],
+                    center_read[i])
+
+        frame_orientation_read=plan.DetectorInformationSequence[0].ImageOrientationPatient
+        for i in range(0,6):
+            if self.frame_orientation[i] == 0:
+                self.frame_orientation[i] = float(frame_orientation_read[i])
+            if abs(self.frame_orientation[i]-frame_orientation_read[i]) > 1e-3:
+                print 'Image orientation mismatch {.2f}/{.2f}'.format(
+                    self.frame_rotation[i], frame_orientation_read[i])
+
+
+
+
+        self.frame_data[:,:,:,ifi]=np.transpose(plan.pixel_array,self.axisShift)
+
+        return True
+
+    def readCTFile(self,g):
+
+        try:
+            plan = dicom.read_file(g)
+        except:
+            print ("{}: Not a dicom file").format(g)
+            return False
+
+
+        if plan.Modality != 'CT':
+            print ('{}: Not a CT file').format(g)
+            return False
 
-        origin=re.sub('([^:/])://(.*)$',r'\1',mypath)
-        onlyfiles=self.filelist(mypath)
-        for f in onlyfiles:
-            print '{}:'.format(f)
+        #this doesn't work in 2019 data version
+        #if re.match("AC",plan.SeriesDescription) == None:
+        #    print (plan.SeriesDescription)
+        #    print ('Not a AC file')
+        #    continue
+        try:
+            iType=plan.ImageType
+        except:
+            print "Image type not found"
+            return False
 
-            g,ok=self.getfile(origin,f)
-            if not(ok):
-                return
+        if iType[3].find("SPI")<0:
+            print "Not a spiral image"
+            return False
 
-            try:
-                plan = dicom.read_file(g)
-            except:
-                print ("Not a dicom file")
-                continue
-            try:
-                nframe=plan[0x0019,0x10a5].value;
-            except:
-                print ("Tag not found;")
-                continue
-            if not (type(nframe) is list) :
-                print("nframe not a list")
-                continue
 
-    #this is the "master" file where data on other files can be had
-    #here we found out the duration of the frame and their distribution through
-    #phases and cycles
-            print('Found master file')
 
-            for i in range(1,len(nframe)):
-                nframe[i]+=nframe[i-1]
+        print '.',
+        self.ct_data.append(plan.pixel_array)
+        self.ct_idx.append(plan.InstanceNumber)
+        self.ct_z.append(plan.ImagePositionPatient[2])
+
+        pixel_size_read=[plan.PixelSpacing[0],plan.PixelSpacing[1],
+            plan.SliceThickness]
+
+
+        for i in range(0,3):
+            if self.ct_pixel_size[i] == 0:
+                self.ct_pixel_size[i] = float(pixel_size_read[i])
+            if abs(self.ct_pixel_size[i]-pixel_size_read[i]) > 1e-3:
+                print 'Pixel size mismatch {.2f}/{.2f}'.format(self.ct_pixel_size[i],
+                    pixel_size_read[i])
 
-            print(nframe)
+        for i in range(0,2):
+            if self.ct_center[i] == 0:
+                self.ct_center[i] = float(plan.ImagePositionPatient[i])
+            if abs(self.ct_center[i]-plan.ImagePositionPatient[i]) > 1e-3:
+                    print 'Image center mismatch {.2f}/{.2f}'.format(self.ct_center[i],
+                    plan.ImagePositionPatient[i])
+    #not average, but minimum (!) why??
 
-    #nframe now holds for index i total number of frames collected up
-    #to the end of each phase
+        if plan.ImagePositionPatient[2]<self.ct_center[2]:
+            self.ct_center[2]=plan.ImagePositionPatient[2]
 
-            frame_start=plan[0x0019,0x10a7].value
-            frame_stop=plan[0x0019,0x10a8].value
-            frame_duration=plan[0x0019,0x10a9].value
-            break
-    #print "rep [{}] start [{}] stop [{}] duration [{}]".format(
-    #len(rep),len(rep_start),len(rep_stop),len(rep_duration))
+        for i in range(0,6):
+            if self.ct_orientation[i] == 0:
+                self.ct_orientation[i] = float(plan.ImageOrientationPatient[i])
+            if abs(self.ct_orientation[i]-plan.ImageOrientationPatient[i]) > 1e-3:
+                print 'Image orientation mismatch {0:.2f}/{1:.2f}'.format(self.ct_orientation[i],
+                plan.ImageOrientationPatient[i])
 
-#select AC reconstructed data
-        frame_time=np.zeros(nframe[-1]);
-        frame_data=np.empty([1,1,1,nframe[-1]])
-        center = [0,0,0]
-        pixel_size =[0,0,0]
-        frame_orientation=[0,0,0,0,0,0]
+        return True
+
+    def readMasterDirectory(self,mypath):
+        self.axisShift=(2,1,0)
+
+        origin=re.sub('([^:/])://(.*)$',r'\1',mypath)
+
+        onlyfiles=self.filelist(mypath)
         for f in onlyfiles:
+            print '{}:'.format(f)
 
             g,ok=self.getfile(origin,f)
             if not(ok):
-                continue
+                return
 
-            try:
-                plan = dicom.read_file(g)
-            except:
-                print ("Not a dicom file")
-                continue
+            if self.readMasterFile(g):
+                break
 
-            try:
-                pf=plan[0x0018,0x5020]
-            except:
-                print ("ProcessingFunction not found")
-                continue
 
-            try:
-                phase=plan[0x0035,0x1005].value
-                cycle=plan[0x0035,0x1004].value
-            except:
-                print ("Phase/Cycle tag not found")
-                continue
-
-            #convert phase/cycle to frame index
-            off=0
-            if phase > 1:
-                off=nframe[phase-2]
-            ifi=off+cycle-1
-
-    #from values in the master file determine frame time
-    #(as the mid point between starting and ending the frame)
-            frame_time[ifi]=0.5*(frame_start[ifi]+frame_stop[ifi]); #in ms
-
-            print "({},{}) converted to {} at {} for {}".format(
-                phase,cycle,ifi,frame_time[ifi],frame_duration[ifi])
-
-    #play with pixel data
-            if frame_data.shape[0] == 1:
-                sh=np.transpose(plan.pixel_array,axisShift).shape;
-                sh=list(sh)
-                sh.append(nframe[-1])
-                frame_data=np.empty(sh)
-                print "Setting frame_data to",sh
-
-            #check & update pixel size
-            pixel_size_read=[plan.PixelSpacing[0],plan.PixelSpacing[1],
-                        plan.SliceThickness]
-
-            for i in range(0,3):
-                if pixel_size[i] == 0:
-                    pixel_size[i] = float(pixel_size_read[i])
-                if abs(pixel_size[i]-pixel_size_read[i]) > 1e-3:
-                    print 'Pixel size mismatch {.2f}/{.2f}'.format(pixel_size[i],
-                    pixel_size_read[i])
+    def readNMDirectory(self,mypath):
 
-            center_read=plan.DetectorInformationSequence[0].ImagePositionPatient
-            print "Stored center at ({0},{1},{2})".format(center[0],center[1],center[2])
-            print "Read   center at ({0},{1},{2})".format(center_read[0],center_read[1],center_read[2])
-            for i in range(0,3):
-                if center[i] == 0:
-                    center[i] = float(center_read[i])
-                if abs(center[i]-center_read[i]) > 1e-3:
-                    print 'Image center mismatch {.2f}/{.2f}'.format(center[i],
-                        center_read[i])
+        onlyfiles=self.filelist(mypath)
+        origin=re.sub('([^:/])://(.*)$',r'\1',mypath)
+        for f in onlyfiles:
 
-            frame_orientation_read=plan.DetectorInformationSequence[0].ImageOrientationPatient
-            for i in range(0,6):
-                if frame_orientation[i] == 0:
-                    frame_orientation[i] = float(frame_orientation_read[i])
-                if abs(frame_orientation[i]-frame_orientation_read[i]) > 1e-3:
-                    print 'Image orientation mismatch {.2f}/{.2f}'.format(
-                        frame_rotation[i], frame_orientation_read[i])
+            g,ok=self.getfile(origin,f)
+            if not(ok):
+                continue
 
+            self.readNMFile(g)
 
 
 
-            frame_data[:,:,:,ifi]=np.transpose(plan.pixel_array,axisShift)
 
-    #print('Orientation: ({0:.2f},{1:.2f},{2:.2f}),({3:.2f},{4:.2f},{5:.2f})').format( \
-    #    frame_orientation[0],frame_orientation[1],frame_orientation[2], \
-    #    frame_orientation[3],frame_orientation[4],frame_orientation[5])
 
-        return [frame_data,frame_time,center,pixel_size,frame_orientation]
+        return [self.frame_data,self.frame_time,self.center,
+            self.pixel_size,self.frame_orientation]
 
-    def read_CT(self,mypath):
+    def readCTDirectory(self,mypath):
         onlyfiles=self.filelist(mypath)
         origin=re.sub('([^:/])://(.*)$',r'\1',mypath)
 
-        ct_data = []
-        ct_idx = []
-        ct_z = []
-        ct_pixel_size = [0,0,0]
-        ct_center = [0,0,0]
-        ct_center[2]=1e30
-        ct_orientation=[0,0,0,0,0,0]
+        self.ct_data = []
+        self.ct_idx = []
+        self.ct_z = []
+        self.ct_pixel_size = [0,0,0]
+        self.ct_center = [0,0,0]
+        self.ct_center[2]=1e30
+        self.ct_orientation=[0,0,0,0,0,0]
         for f in onlyfiles:
             print '{}:'.format(f)
 
@@ -239,78 +318,18 @@ class parseDicomLogic(ScriptedLoadableModuleLogic):
             if not(ok):
                 return
 
-            try:
-                plan = dicom.read_file(g)
-            except:
-                print ("Not a dicom file")
-                continue
-
-            if plan.Modality != 'CT':
-                print ('Not a CT file')
-                continue
+            self.readCTFile(g)
 
-        #this doesn't work in 2019 data version
-        #if re.match("AC",plan.SeriesDescription) == None:
-        #    print (plan.SeriesDescription)
-        #    print ('Not a AC file')
-        #    continue
-            try:
-                iType=plan.ImageType
-            except:
-                print "Image type not found"
-                continue;
-
-            if iType[3].find("SPI")<0:
-                print "Not a spiral image"
-                continue;
-
-
-    #a slice of pure CT
-            print '.',
-            ct_data.append(plan.pixel_array)
-            ct_idx.append(plan.InstanceNumber)
-            ct_z.append(plan.ImagePositionPatient[2])
-
-            pixel_size_read=[plan.PixelSpacing[0],plan.PixelSpacing[1],
-                plan.SliceThickness]
-
-
-            for i in range(0,3):
-                if ct_pixel_size[i] == 0:
-                    ct_pixel_size[i] = float(pixel_size_read[i])
-                if abs(ct_pixel_size[i]-pixel_size_read[i]) > 1e-3:
-                    print 'Pixel size mismatch {.2f}/{.2f}'.format(ct_pixel_size[i],
-                        pixel_size_read[i])
-
-            for i in range(0,2):
-                if ct_center[i] == 0:
-                    ct_center[i] = float(plan.ImagePositionPatient[i])
-                if abs(ct_center[i]-plan.ImagePositionPatient[i]) > 1e-3:
-                        print 'Image center mismatch {.2f}/{.2f}'.format(ct_center[i],
-                        plan.ImagePositionPatient[i])
-        #not average, but minimum (!) why??
-
-            if plan.ImagePositionPatient[2]<ct_center[2]:
-                ct_center[2]=plan.ImagePositionPatient[2]
-
-            for i in range(0,6):
-                if ct_orientation[i] == 0:
-                    ct_orientation[i] = float(plan.ImageOrientationPatient[i])
-                if abs(ct_orientation[i]-plan.ImageOrientationPatient[i]) > 1e-3:
-                    print 'Image orientation mismatch {0:.2f}/{1:.2f}'.format(ct_orientation[i],
-                    plan.ImageOrientationPatient[i])
-
-        print
-        nz=len(ct_idx)
-    #not average, again
-    #ct_center[2]/=nz
-        sh=ct_data[-1].shape
+        nz=len(self.ct_idx)
+        #not average, again
+        #ct_center[2]/=nz
+        sh=self.ct_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((ct_z[k]-ct_center[2])/ct_pixel_size[2]))
-            data_array[:,:,kp]=np.transpose(ct_data[k])
+            kp=int(np.round((self.ct_z[k]-self.ct_center[2])/self.ct_pixel_size[2]))
+            data_array[:,:,kp]=np.transpose(self.ct_data[k])
 
-        return data_array,ct_center,ct_pixel_size,ct_orientation
+        return data_array,self.ct_center,self.ct_pixel_size,self.ct_orientation