Version that plots segment data

cardiacSPECT/cardiacSPECT.py

@@ -61,15 +61,15 @@ class cardiacSPECTWidget(ScriptedLoadableModuleWidget):
     # Add vertical spacer
     self.layout.addStretch(1)
 
-    addFrameButton=qt.QPushButton("Add Frame")
-    addFrameButton.toolTip="Add frame to VTK"
-    dataFormLayout.addWidget(addFrameButton)
-    addFrameButton.connect('clicked(bool)',self.onAddFrameButtonClicked)
+    #addFrameButton=qt.QPushButton("Add Frame")
+    #addFrameButton.toolTip="Add frame to VTK"
+    #dataFormLayout.addWidget(addFrameButton)
+    #addFrameButton.connect('clicked(bool)',self.onAddFrameButtonClicked)
 
-    addCTButton=qt.QPushButton("Add CT")
-    addCTButton.toolTip="Add CT to VTK"
-    dataFormLayout.addWidget(addCTButton)
-    addCTButton.connect('clicked(bool)',self.onAddCTButtonClicked)
+    #addCTButton=qt.QPushButton("Add CT")
+    #addCTButton.toolTip="Add CT to VTK"
+    #dataFormLayout.addWidget(addCTButton)
+    #addCTButton.connect('clicked(bool)',self.onAddCTButtonClicked)
 
     #
     # Parameters Area
@@ -85,10 +85,72 @@ class cardiacSPECTWidget(ScriptedLoadableModuleWidget):
     #
     # check box to trigger taking screen shots for later use in tutorials
     #
-    self.time_frame_select=qt.QLineEdit()
-    self.time_frame_select.setText("2")
+    hbox1=qt.QHBoxLayout()
+
+    frameLabel = qt.QLabel()
+    frameLabel.setText("Select frame")
+    hbox1.addWidget(frameLabel)
+
+    self.time_frame_select=qt.QSlider(qt.Qt.Horizontal)
+    self.time_frame_select.valueChanged.connect(self.onTimeFrameSelect)
+
+    #self.time_frame_select.connect('valueChanged()', self.onTimeFrameSelect)
+    self.time_frame_select.setMinimum(0)
+    self.time_frame_select.setMaximum(0)
+    self.time_frame_select.setValue(0)
+    self.time_frame_select.setTickPosition(qt.QSlider.TicksBelow)
+    self.time_frame_select.setTickInterval(5)
     self.time_frame_select.toolTip = "Select the time frame"
-    parametersFormLayout.addRow(self.time_frame_select)
+    hbox1.addWidget(self.time_frame_select)
+
+    parametersFormLayout.addRow(hbox1)
+
+    hbox2 = qt.QHBoxLayout()
+
+    meanROILabel = qt.QLabel()
+    meanROILabel.setText("MeanROI")
+    hbox2.addWidget(meanROILabel)
+
+    self.meanROIVolume = qt.QLineEdit()
+    self.meanROIVolume.setText("testVolume15")
+    hbox2.addWidget(self.meanROIVolume)
+
+    self.meanROISegment = qt.QLineEdit()
+    self.meanROISegment.setText("Segment_1")
+    hbox2.addWidget(self.meanROISegment)
+
+    computeMeanROI = qt.QPushButton("Compute mean ROI")
+    computeMeanROI.connect('clicked(bool)',self.onComputeMeanROIClicked)
+    hbox2.addWidget(computeMeanROI)
+
+    self.meanROIResult = qt.QLineEdit()
+    self.meanROIResult.setText("0")
+    hbox2.addWidget(self.meanROIResult)
+
+
+    parametersFormLayout.addRow(hbox2)
+
+    #row 3
+    hbox3 = qt.QHBoxLayout()
+
+    drawTimePlot=qt.QPushButton("Draw ROI time plot")
+    drawTimePlot.connect('clicked(bool)',self.onDrawTimePlotClicked)
+    hbox3.addWidget(drawTimePlot)
+
+    parametersFormLayout.addRow(hbox3)
+    #dataFormLayout.addWidget(hbox)
+
+    #row 4
+    hbox4 = qt.QHBoxLayout()
+
+    countSegments=qt.QPushButton("Count segmentation segments")
+    countSegments.connect('clicked(bool)',self.onCountSegmentsClicked)
+    hbox4.addWidget(countSegments)
+
+    self.countSegmentsDisplay=qt.QLineEdit()
+    self.countSegmentsDisplay.setText("0")
+    hbox4.addWidget(self.countSegmentsDisplay)
+    parametersFormLayout.addRow(hbox4)
 
 
     #
@@ -109,6 +171,8 @@ class cardiacSPECTWidget(ScriptedLoadableModuleWidget):
 
     self.logic=cardiacSPECTLogic()
 
+    self.resetPosition=1
+
   def cleanup(self):
     pass
 
@@ -118,14 +182,86 @@ class cardiacSPECTWidget(ScriptedLoadableModuleWidget):
     #imageThreshold = self.imageThresholdSliderWidget.value
 
   def onDataLoadButtonClicked(self):
-      self.logic.loadData()
+      self.logic.loadData(self)
 
-  def onAddFrameButtonClicked(self):
-      it=int(self.time_frame_select.text)
-      self.logic.addFrame(it)
-
-  def onAddCTButtonClicked(self):
-      self.logic.addCT()
+  def onTimeFrameSelect(self):
+       it=self.time_frame_select.value
+       selectionNode = slicer.app.applicationLogic().GetSelectionNode()
+       print("Propagating CT volume")
+       node=slicer.mrmlScene.GetFirstNodeByName("testCT")
+       selectionNode.SetReferenceActiveVolumeID(node.GetID())
+       if self.resetPosition==1:
+          self.resetPosition=0
+          slicer.app.applicationLogic().PropagateVolumeSelection(1)
+       else:
+          slicer.app.applicationLogic().PropagateVolumeSelection(0)
+       print("Propagating SPECT volume")
+       nodeName='testVolume'+str(it)
+       node=slicer.mrmlScene.GetFirstNodeByName(nodeName)
+       selectionNode.SetSecondaryVolumeID(node.GetID())
+       slicer.app.applicationLogic().PropagateForegroundVolumeSelection(0)
+       node.GetDisplayNode().SetAndObserveColorNodeID('vtkMRMLColorTableNodeRed')
+       lm = slicer.app.layoutManager()
+       sliceLogic = lm.sliceWidget('Red').sliceLogic()
+       compositeNode = sliceLogic.GetSliceCompositeNode()
+       compositeNode.SetForegroundOpacity(0.5)
+
+       #make sure the viewer is matched to the volume
+       print("Done")
+        #to access sliceLogic (slice control) use
+        #lcol=slicer.app.layoutManager().mrmlSliceLogics() (vtkCollection)
+        #vtkMRMLSliceLogic are named by colors (Red,Green,Blue)
+
+  def onComputeMeanROIClicked(self):
+        s=self.logic.meanROI(self.meanROIVolume.text,self.meanROISegment.text)
+        self.meanROIResult.setText(str(s))
+
+  def onDrawTimePlotClicked(self):
+        n=self.time_frame_select.maximum
+        ft=self.logic.frame_time
+
+        #find number of segments
+        ns = self.logic.countSegments()
+
+        #add the chart node
+        cn = slicer.mrmlScene.AddNode(slicer.vtkMRMLChartNode())
+
+        for j in range(0,ns):
+            segment="Segment_"+str(j+1)
+            #add node for data
+            dn = slicer.mrmlScene.AddNode(slicer.vtkMRMLDoubleArrayNode())
+            a = dn.GetArray()
+            a.SetNumberOfTuples(n)
+            for i in range(0,n):
+                vol="testVolume"+str(i)
+                fx=ft[i]
+                fy=self.logic.meanROI(vol,segment)
+                a.SetComponent(i, 0, fx)
+                a.SetComponent(i, 1, fy)
+                a.SetComponent(i, 2, 0)
+                print("({0:.2f}:{1:.2f})".format(fx,fy))
+
+            cn.AddArray(segment, dn.GetID())
+
+        cn.SetProperty('default', 'title', 'ROI time plot')
+        cn.SetProperty('default', 'xAxisLabel', 'time [ms]')
+        cn.SetProperty('default', 'yAxisLabel', 'Activity (arb)')
+
+        #update the chart node
+        cvns = slicer.mrmlScene.GetNodesByClass('vtkMRMLChartViewNode')
+        cvns.InitTraversal()
+        cvn = cvns.GetNextItemAsObject()
+        cvn.SetChartNodeID(cn.GetID())
+
+  def onCountSegmentsClicked(self):
+        self.countSegmentsDisplay.setText(self.logic.countSegments())
+
+  #def onAddFrameButtonClicked(self):
+#      it=int(self.time_frame_select.text)
+#      self.logic.addFrame(it)
+
+ # def onAddCTButtonClicked(self):
+#      self.logic.addCT()
 #
 #
 # cardiacSPECTLogic
@@ -141,7 +277,7 @@ class cardiacSPECTLogic(ScriptedLoadableModuleLogic):
   https://github.com/Slicer/Slicer/blob/master/Base/Python/slicer/ScriptedLoadableModule.py
   """
 
-  def loadData(self):
+  def loadData(self,widget):
     startDir=os.environ['HOME']
     inputDir=qt.QFileDialog.getExistingDirectory(None,
         'Select DICOM directory',startDir)
@@ -151,30 +287,45 @@ class cardiacSPECTLogic(ScriptedLoadableModuleLogic):
     sys.path.append(mypath)
     import parseDicom as pd
 
-    self.frame_data, self.frame_time, self.frame_center, \
+    self.frame_data, self.frame_time, self.frame_origin, \
         self.frame_pixel_size, self.frame_orientation=pd.read_dynamic_SPECT(inputDir)
 
-    self.ct_data,self.ct_center,self.ct_pixel_size, \
+    self.ct_data,self.ct_origin,self.ct_pixel_size, \
         self.ct_orientation=pd.read_CT(inputDir)
 
     self.ct_orientation=vi.completeOrientation(self.ct_orientation)
     self.frame_orientation=vi.completeOrientation(self.frame_orientation)
+
+    self.addCT()
+    self.addFrames()
+
+    widget.time_frame_select.setMaximum(self.frame_data.shape[3]-1)
+
     #additional message via qt
     qt.QMessageBox.information(
         slicer.util.mainWindow(),
         'Slicer Python','Data loaded')
 
 
-  def addNode(self,nodeName,v, orientation):
+  def addNode(self,nodeName,v, origin, pixel_size, orientation, dataType):
+
+       # if dataType=0 it is CT data, which gets propagated to background an is
+       #used to fit the view field dimensions
+
+       # if dataType=1, it is SPECT data, which gets propagated to foreground
+       #and is not fit; keeping window set from CT
+
        #nodeName='testVolume'+str(it)
        newNode=slicer.vtkMRMLScalarVolumeNode()
        newNode.SetName(nodeName)
 
        #pixel_size=[0,0,0]
-       pixel_size=v.GetSpacing()
-       print(pixel_size)
+       #pixel_size=v.GetSpacing()
+       #print(pixel_size)
        #origin=[0,0,0]
-       origin=v.GetOrigin()
+       #origin=v.GetOrigin()
+       v.SetOrigin([0,0,0])
+       v.SetSpacing([1,1,1])
        ijkToRAS = vtk.vtkMatrix4x4()
        #think how to do this with image orientation
 
@@ -187,27 +338,92 @@ class cardiacSPECTLogic(ScriptedLoadableModuleLogic):
        newNode.SetIJKToRASMatrix(ijkToRAS)
        newNode.SetAndObserveImageData(v)
        slicer.mrmlScene.AddNode(newNode)
-       selectionNode = slicer.app.applicationLogic().GetSelectionNode()
-       selectionNode.SetReferenceActiveVolumeID(newNode.GetID())
-       slicer.app.applicationLogic().PropagateVolumeSelection(0)
 
-  def addFrame(self,it):
+  def addFrames(self):
        #convert data from numpy.array to vtkImageData
        #use time point it
-       frame_data=self.frame_data[:,:,:,it];
-       nodeName='testVolume'+str(it)
-       self.addNode(nodeName,
-         vi.numpyToVTK3D(frame_data,self.frame_center,
-           self.frame_pixel_size),self.frame_orientation)
+       for it in range(0,self.frame_data.shape[3]):
+           frame_data=self.frame_data[:,:,:,it];
+           nodeName='testVolume'+str(it)
+           self.addNode(nodeName,
+              vi.numpyToVTK3D(frame_data,self.frame_origin,
+                self.frame_pixel_size),
+              self.frame_origin,
+              self.frame_pixel_size,
+              self.frame_orientation,1)
 
   def addCT(self):
        nodeName='testCT'
        self.addNode(nodeName,
             vi.numpyToVTK3D(self.ct_data,
-                self.ct_center,self.ct_pixel_size),
-            self.ct_orientation)
-
-
+                self.ct_origin,self.ct_pixel_size),
+            self.ct_origin,self.ct_pixel_size,
+            self.ct_orientation,0)
+
+  def meanROI(self, volName1, segment):
+    s=0
+
+    #get the segmentation mask
+    fNode=slicer.mrmlScene.GetNodesByClass("vtkMRMLSegmentationNode").GetItemAsObject(0)
+    segNode=slicer.vtkMRMLSegmentationNode.SafeDownCast(fNode)
+
+    #no python bindings for vtkSegmentation
+
+    #if segNode.GetSegmentation().GetNumberOfSegments()==0 :
+    #    print("No segments available")
+    #    return 0
+
+    mask = segNode.GetBinaryLabelmapRepresentation(segment)
+    if mask==None:
+        print("Segment {} not found".format(segment))
+        return s
+
+      #segNode.GetSegmentation().GetNthSegmentID(0))
+
+    #get mask at (x,y,z)
+    #mask.GetPointData().GetScalars().GetTuple1(mask.FindPoint([x,y,z]))
+
+    #get the image data
+    dataNode=slicer.mrmlScene.GetFirstNodeByName(volName1)
+    dataImage=dataNode.GetImageData()
+    # use IJK2RAS to get global coordinates
+    ijkToRas = vtk.vtkMatrix4x4()
+    dataNode.GetIJKToRASMatrix(ijkToRas)
+
+    #iterate over volume pixelData
+    n=dataImage.GetNumberOfPoints()
+
+    for i in range(0,n):
+
+      #get global coordinates of point i
+      [ix,iy,iz]=dataImage.GetPoint(i)
+      position_ijk=[ix, iy, iz, 1]
+      #ras are global coordinates (in mm)
+      position_ras=ijkToRas.MultiplyPoint(position_ijk)
+
+      #find point in mask with the same global coordinates
+      fr=[position_ras[0],position_ras[1],position_ras[2]]
+      maskValue=mask.GetPointData().GetScalars().GetTuple1(mask.FindPoint(fr))
+
+      if maskValue == 0:
+          continue
+
+      #use maskValue to project ROI data
+      s+=maskValue*dataImage.GetPointData().GetScalars().GetTuple1(i)
+
+    return s/n
+
+  def countSegments(self):
+    fNode=slicer.mrmlScene.GetNodesByClass("vtkMRMLSegmentationNode").GetItemAsObject(0)
+    segNode=slicer.vtkMRMLSegmentationNode.SafeDownCast(fNode)
+    i=1
+    while 1:
+        segName="Segment_"+str(i)
+        mask = segNode.GetBinaryLabelmapRepresentation(segName)
+        if mask==None:
+           break
+        i+=1
+    return i-1
 
 class cardiacSPECTTest(ScriptedLoadableModuleTest):
   """

cardiacSPECT/cardiacSPECT/cardiacSPECT.py

@@ -1,228 +0,0 @@
-import os
-import sys
-import unittest
-import vtk, qt, ctk, slicer
-from slicer.ScriptedLoadableModule import *
-import logging
-
-#
-# cardiacSPECT
-#
-
-class cardiacSPECT(ScriptedLoadableModule):
-  """Uses ScriptedLoadableModule base class, available at:
-  https://github.com/Slicer/Slicer/blob/master/Base/Python/slicer/ScriptedLoadableModule.py
-  """
-  def __init__(self, parent):
-    ScriptedLoadableModule.__init__(self, parent)
-    parent.title = "Cardiac SPECT"
-    parent.categories = ["Examples"]
-    parent.dependencies = []
-    parent.contributors = ["Andrej Studen (FMF/JSI)"] # replace with "Firstname Lastname (Org)"
-    parent.helpText = """
-    Load dynamic cardiac SPECT data to Slicer
-    """
-    parent.acknowledgementText = """
-    This module was developed within the frame of the ARRS sponsored medical
-    physics research programe to investigate quantitative measurements of cardiac
-    function using sestamibi-like tracers
-    """ # replace with organization, grant and thanks.
-    self.parent.helpText += self.getDefaultModuleDocumentationLink()
-    self.parent = parent
-
-#
-# cardiacSPECTWidget
-#
-
-class cardiacSPECTWidget(ScriptedLoadableModuleWidget):
-  """Uses ScriptedLoadableModuleWidget base class, available at:
-  https://github.com/Slicer/Slicer/blob/master/Base/Python/slicer/ScriptedLoadableModule.py
-  """
-  def __init__(self, parent):
-      ScriptedLoadableModuleWidget.__init__(self, parent)
-      self.logic=cardiacSPECTLogic()
-
-  def setup(self):
-    ScriptedLoadableModuleWidget.setup(self)
-
-    # Instantiate and connect widgets ...
-    dataButton = ctk.ctkCollapsibleButton()
-    dataButton.text = "Data"
-    self.layout.addWidget(dataButton)
-
-    # Layout within the sample collapsible button
-    dataFormLayout = qt.QFormLayout(dataButton)
-
-    dataLoadButton = qt.QPushButton("Load")
-    dataLoadButton.toolTip="Load data from DICOM"
-    dataFormLayout.addWidget(dataLoadButton)
-    dataLoadButton.connect('clicked(bool)',self.onDataLoadButtonClicked)
-
-    # Set local var as instance attribute
-    self.dataLoadButton = dataLoadButton
-
-    # Add vertical spacer
-    self.layout.addStretch(1)
-
-    dataImportButton=qt.QPushButton("Import")
-    dataImportButton.toolTip="Import time frame to VTK"
-    dataFormLayout.addWidget(dataImportButton)
-    dataImportButton.connect('clicked(bool)',self.onImportButtonClicked)
-
-    #
-    # Parameters Area
-    #
-
-    parametersCollapsibleButton = ctk.ctkCollapsibleButton()
-    parametersCollapsibleButton.text = "Parameters"
-    self.layout.addWidget(parametersCollapsibleButton)
-
-    # Layout within the dummy collapsible button
-    parametersFormLayout = qt.QFormLayout(parametersCollapsibleButton)
-
-    #
-    # check box to trigger taking screen shots for later use in tutorials
-    #
-    self.time_frame_select=qt.QLineEdit()
-    self.time_frame_select.setText("2")
-    self.time_frame_select.toolTip = "Select the time frame"
-    parametersFormLayout.addRow(self.time_frame_select)
-
-
-    #
-    # Apply Button
-    #
-    self.applyButton = qt.QPushButton("Apply")
-    self.applyButton.toolTip = "Run the algorithm."
-    self.applyButton.enabled = False
-    parametersFormLayout.addRow(self.applyButton)
-
-    # connections
-    self.applyButton.connect('clicked(bool)', self.onApplyButton)
-    #self.inputSelector.connect("currentNodeChanged(vtkMRMLNode*)", self.onSelect)
-    #self.outputSelector.connect("currentNodeChanged(vtkMRMLNode*)", self.onSelect)
-
-    # Add vertical spacer
-    self.layout.addStretch(1)
-
-
-  def cleanup(self):
-    pass
-
-  def onApplyButton(self):
-    #logic = cardiacSPECTLogic()
-    #imageThreshold = self.imageThresholdSliderWidget.value
-
-  def onDataLoadButtonClicked(self):
-      self.logic.load()
-
-  def onImportButtonClicked(self):
-      it=int(self.time_frame_select.text())
-      self.logic.setVolume(it)
-#
-# cardiacSPECTLogic
-#
-
-class cardiacSPECTLogic(ScriptedLoadableModuleLogic):
-  """This class should implement all the actual
-  computation done by your module.  The interface
-  should be such that other python code can import
-  this class and make use of the functionality without
-  requiring an instance of the Widget.
-  Uses ScriptedLoadableModuleLogic base class, available at:
-  https://github.com/Slicer/Slicer/blob/master/Base/Python/slicer/ScriptedLoadableModule.py
-  """
-
-  def load(self):
-    startDir=os.environ['HOME']
-    inputDir=qt.QFileDialog.getExistingDirectory(None,
-        'Select DICOM directory',startDir)
-
-    #use another script from the same file
-    mypath=os.environ['HOME']+'/software/build/dynamicSPECT'
-    sys.path.append(mypath)
-    import parseData as pd
-
-    self.frame_data,self.frame_time=pd.read_dynamic_SPECT(inputDir)
-
-    #additional message via qt
-    qt.QMessageBox.information(
-        slicer.util.mainWindow(),
-        'Slicer Python','Data loaded')
-
-  def setVolume(self,it):
-       #convert data from numpy.array to vtkImageData
-       #use time point it
-       frame_data=self.frame_data;
-
-       v=vtk.vtkImageData()
-       v.SetDimensions(
-            frame_data.shape[0],
-            frame_data.shape[1],
-            frame_data.shape[2])
-       v.SetOrigin(0.0*frame_data.shape[0],
-                0.0*frame_data.shape[1],
-                0.0*frame_data.shape[2])
-       sp = 1.0
-       v.SetSpacing(sp, sp, sp)
-       scalars = vtk.vtkShortArray()
-       for k in range(0,frame_data.shape[2]):
-           kOffset = k * frame_data.shape[1]*self.frame_data.shape[0]
-           for j in range(0,frame_data.shape[1]):
-               jOffset = j * frame_data.shape[0]
-               for i in range(0,frame_data.shape[0]):
-                   offset = i + jOffset + kOffset;
-                   scalars.InsertTuple1(offset,frame_data[i,j,k,it])
-       v.GetPointData().SetScalars(scalars);
-       self.v=v
-
-       nodeName='testVolume'+str(it)
-       newNode=slicer.vtkMRMLScalarVolumeNode()
-       newNode.SetName(nodeName)
-       ijkToRAS = vtk.vtkMatrix4x4()
-       ijkToRAS.Identity()
-       newNode.SetIJKToRASMatrix(ijkToRAS)
-       newNode.SetAndObserveImageData(v)
-       slicer.mrmlScene.AddNode(newNode)
-       selectionNode = slicer.app.applicationLogic().GetSelectionNode()
-       selectionNode.SetReferenceActiveVolumeID(newNode.GetID())
-       slicer.app.applicationLogic().PropagateVolumeSelection(0)
-
-
-
-class cardiacSPECTTest(ScriptedLoadableModuleTest):
-  """
-  This is the test case for your scripted module.
-  Uses ScriptedLoadableModuleTest base class, available at:
-  https://github.com/Slicer/Slicer/blob/master/Base/Python/slicer/ScriptedLoadableModule.py
-  """
-
-  def setUp(self):
-    """ Do whatever is needed to reset the state - typically a scene clear will be enough.
-    """
-    slicer.mrmlScene.Clear(0)
-
-  def runTest(self):
-    """Run as few or as many tests as needed here.
-    """
-    self.setUp()
-    self.test_cardiacSPECT1()
-
-  def test_cardiacSPECT1(self):
-    """ Ideally you should have several levels of tests.  At the lowest level
-    tests should exercise the functionality of the logic with different inputs
-    (both valid and invalid).  At higher levels your tests should emulate the
-    way the user would interact with your code and confirm that it still works
-    the way you intended.
-    One of the most important features of the tests is that it should alert other
-    developers when their changes will have an impact on the behavior of your
-    module.  For example, if a developer removes a feature that you depend on,
-    your test should break so they know that the feature is needed.
-    """
-
-    self.delayDisplay("Starting the test")
-    #
-    # first, get some data
-    #
-
-    self.delayDisplay('Test passed!')

parseDicom.py

@@ -181,6 +181,7 @@ def read_CT(mypath):
             print ('Not a AC file')
             continue
         #a slice of pure CT
+        print '.',
         ct_data.append(plan.pixel_array)
         ct_idx.append(plan.InstanceNumber)
         #ct_center.append(plan.ImagePositionPatient)
@@ -214,7 +215,7 @@ def read_CT(mypath):
                 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

vtkInterface.py

@@ -27,12 +27,14 @@ def numpyToVTK3D(numpy_array, origin, spacing):
     v.SetSpacing(spacing[0], spacing[1], spacing[2])
     scalars = vtk.vtkShortArray()
     for k in range(0,numpy_array.shape[2]):
+       #print '.',
        kOffset = k * numpy_array.shape[1]*numpy_array.shape[0]
        for j in range(0,numpy_array.shape[1]):
            jOffset = j * numpy_array.shape[0]
            for i in range(0,numpy_array.shape[0]):
                offset = i + jOffset + kOffset;
                scalars.InsertTuple1(offset,numpy_array[i,j,k])
+    #print
     v.GetPointData().SetScalars(scalars)
     return v