Display Issue with Volume Node Created Programmatically

Support
2

Creating a node is an expensive operation. You create 1183 volume nodes (and associated display and storage nodes) just to copy some voxels into them and delete them. This increases execution time from about 1 sec to 30sec.

Probably the issue was how you copied/moved value ranges between numpy arrays. As memory of image arrays managed by VTK, if you apply complex range operations, the safest is to make a copy of the numpy array and update the volume with the processed array. Only use getNode and array in testing code, not in actual production code, as looking up nodes by name is not reliable (multiple nodes may have the same name). After cleaning up and simplifying the code a bit, you can get this version that works correctly:

from __future__ import division
import numpy as np
import math
import SampleData

def createEmptyVolume(imageSize, imageSpacing, nodeName):
  voxelType=vtk.VTK_SHORT
  # Create an empty image volume
  imageData=vtk.vtkImageData()
  imageData.SetDimensions(imageSize)
  imageData.AllocateScalars(voxelType, 1)
  thresholder=vtk.vtkImageThreshold()
  thresholder.SetInputData(imageData)
  thresholder.SetInValue(0)
  thresholder.SetOutValue(0)
  thresholder.Update()
  # Create volume node
  volumeNode=slicer.mrmlScene.AddNewNodeByClass("vtkMRMLScalarVolumeNode", nodeName)
  volumeNode.SetSpacing(imageSpacing)
  volumeNode.SetAndObserveImageData(thresholder.GetOutput())
  volumeNode.CreateDefaultDisplayNodes()
  volumeNode.CreateDefaultStorageNode()
  return volumeNode

def performVesselness(inputNode, cutSpacing, createTempVolumes):
  inputImageArray = slicer.util.arrayFromVolume(inputNode)
  dim = inputNode.GetImageData().GetDimensions()
  spacing = inputNode.GetSpacing()
  numOfSlices = [0,0,0]
  for i in range(3):
    numOfSlices[i] = int(math.ceil(dim[i] / cutSpacing[i]))
  print(numOfSlices)
  outputNode = createEmptyVolume(dim, spacing, 'VesselnessFiltered')
  outputImageArray = slicer.util.arrayFromVolume(outputNode).copy()
  for ii in range(numOfSlices[0]):
    xMin = ii * cutSpacing[0]
    xMax = min((ii + 1) * cutSpacing[0], dim[0])
    for jj in range(numOfSlices[1]):
      yMin = jj * cutSpacing[1]
      yMax = min((jj + 1) * cutSpacing[1], dim[1])
      for kk in range(numOfSlices[2]):
        zMin = kk * cutSpacing[2]
        zMax = min((kk + 1) * cutSpacing[2], dim[2])
        tileDim = [xMax-xMin, yMax-yMin, zMax-zMin]
        if createTempVolumes:
          tempVolume = createEmptyVolume(tileDim, spacing, 'tempV')
          tempVolumeArray = slicer.util.arrayFromVolume(tempVolume)
          tempVolumeArray[:] = inputImageArray[zMin:zMax, yMin:yMax, xMin:xMax]
          outputImageArray[zMin:zMax, yMin:yMax, xMin:xMax] = tempVolumeArray
          slicer.mrmlScene.RemoveNode(tempVolume)
        else:
          outputImageArray[zMin:zMax, yMin:yMax, xMin:xMax] = inputImageArray[zMin:zMax, yMin:yMax, xMin:xMax]
  slicer.util.updateVolumeFromArray(outputNode, outputImageArray)
  # Copy origin, spacing, axis directions
  ijkToRAS = vtk.vtkMatrix4x4()
  inputNode.GetIJKToRASMatrix(ijkToRAS)
  outputNode.SetIJKToRASMatrix(ijkToRAS)
  return outputNode

inputVolume = SampleData.SampleDataLogic().downloadMRHead()
outputVolume = performVesselness(inputVolume, [20,20,20], True)  # slow
outputVolume = performVesselness(inputVolume, [20,20,20], False)  # fast
print np.array_equal(array('MR*'), array('VesselnessFiltered'))
3 Likes