Display Issue with Volume Node Created Programmatically

Support
1

I am doing vesselness analysis on some dataset that requires performing all the tasks from the python interpreter. The steps so far are:

  1. Splitting the 3D data volume into smaller 3D blocks (In the attached python code I used MRHead sample data as an example)
  2. For each of the smaller block I perform some analysis (for the sake of simplicity, I just copied the data in the attached python code)

So the resulting volume node (called ‘VesselnessFiltered’) should have exactly the same data array as the sample MRHead node. Indeed I checked numpy.array_equal(array(‘MR*’), array(‘VesselnessFiltered’))
and the return value is True. However, when I load the resulting ‘VesseslnessFiltered’ node into the volume rendering module, I see nothing. So I went to check the volume in the Volume module and I see that the scalar range is 0 to 0 and the histogram also gives nothing. So I assume that it might be because I didn’t properly create a volume node (but I followed this when creating an empty volume so I am really confused now). Also a side question is whether there is an efficient way to split a larger 3D array into smaller ones programmatically and in a custom way in 3D slicer. The attached py file merely copied the array data from a volume node to another volume node but it took around 40 seconds to run. Any help is much appreciated!

%%%%%%% The Python Code %%%%%%%%
from __future__ import division
import numpy as np
import math
import SampleData

def createEmptyVolume(imageSize, imageSpacing, nodeName):
	# imageSize=[512, 512, 512]
	# imageSpacing=[1.0, 1.0, 1.0]
	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)
	# Create volume node
	volumeNode=slicer.vtkMRMLScalarVolumeNode()
	volumeNode.SetSpacing(imageSpacing)
	volumeNode.SetImageDataConnection(thresholder.GetOutputPort())
	volumeNode.SetName(nodeName)
	# Add volume to scene
	slicer.mrmlScene.AddNode(volumeNode)
	displayNode=slicer.vtkMRMLScalarVolumeDisplayNode()
	slicer.mrmlScene.AddNode(displayNode)
	colorNode = slicer.util.getNode('Grey')
	displayNode.SetAndObserveColorNodeID(colorNode.GetID())
	volumeNode.SetAndObserveDisplayNodeID(displayNode.GetID())
	volumeNode.CreateDefaultStorageNode()

def performVesselness(inputNode, cutSpacing, outputNode=None):
	inputImageData = array(inputNode.GetName())
	xDim, yDim, zDim = inputNode.GetImageData().GetDimensions()
	xCutSpacing, yCutSpacing, zCutSpacing = cutSpacing
	xSpacing, ySpacing, zSpacing = inputNode.GetSpacing()
	createEmptyVolume([xDim, yDim, zDim], [xSpacing, ySpacing, zSpacing], 'VesselnessFiltered')
	xNumOfSlices = int(math.ceil(xDim / xCutSpacing))
	yNumOfSlices = int(math.ceil(yDim / yCutSpacing))
	zNumOfSlices = int(math.ceil(zDim / zCutSpacing))
	for ii in range(xNumOfSlices):
		if ii != xNumOfSlices - 1:
			xLocs = range(ii * xCutSpacing , (ii + 1) * xCutSpacing)
		else:
			xLocs = range(ii * xCutSpacing , xDim)

		for jj in range(yNumOfSlices):
			if jj != yNumOfSlices - 1:
				yLocs = range(jj * yCutSpacing , (jj + 1) * yCutSpacing)
			else:
				yLocs = range(jj * yCutSpacing , yDim)

			for kk in range(zNumOfSlices):
				if kk != zNumOfSlices - 1:
					zLocs = range(kk * zCutSpacing , (kk + 1) * zCutSpacing)
				else:
					zLocs = range(kk * zCutSpacing , zDim)

				createEmptyVolume([len(xLocs), len(yLocs), len(zLocs)], [xSpacing, ySpacing, zSpacing], 'tempV')
				array('tempV')[:] = inputImageData[np.ix_(zLocs, yLocs, xLocs)]
				array('VesselnessFiltered')[np.ix_(zLocs, yLocs, xLocs)] = array('tempV')[:]
				slicer.mrmlScene.RemoveNode(getNode('tempV'))

SampleData.SampleDataLogic().downloadMRHead()
performVesselness(getNode('MR*'), [20,20,20])
print np.array_equal(array('MR*'), array('VesselnessFiltered'))
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
3

Thank you Iassoan! Your code solves my problem perfectly. The reason why I created those dummy nodes is that in my real code I would provide these nodes to an internal function (‘computeVesselnessVolume’) of the Slicer-VMTK plugin, which takes volume node instead of data array as input argument. So in this case I guess I have to use expensive node creation procedure (please let me know if there is another way around). Still, I am a little bit curious why via my version of the code the new volume cannot be displayed as the underlying data array is exactly the same (as indicated by the output of numpy.array_equal). Is there something else in the volume node whose value I didn’t set properly that caused the display to fail? Thanks!

4

If you modify voxels using numpy, then when you finished with the changes you need no notify VTK about the changes by calling:

v=getNode('VesselnessFiltered')
v.GetImageData().GetPointData().GetScalars().Modified()
v.Modified()

To see the same volume, you need to set the same window/level in the display node and copy image origin, spacing, and directions from the original volume.

Speed: you definitely should not create 1183 temporary volume nodes and run VMTK filtering 1183 times. Why would you split the image into so many small pieces and run VMTK on these small pieces instead of the whole volume?

1 Like
5

Thank you for your reply! I only copied the data array previously and perhaps that is why the slicer could not display it. As for the speed, performing vmtk filtering on the whole volume would have a tradeoff between losing fine details of the vessel (using stricter filtering parameters) and including non-vessel voxels (using less strict parameters). So I was trying to split the whole volume into smaller ones, apply filtering with strict parameters on each of them, exclude those subvolumes that merely have any vessels and refilter the remaining subvolumes that do have vessels with less strict parameters so that the fine details could be revealed. And yes you are right, my true dataset has a dimension of 640x880x880 and I definitely have to find a better way to do this instead of creating an enormous amount of subvolumes…

6

Volume subdivision may indeed decrease the computation time but you have to keep the number of subvolumes small (up to a few ten or so).

You may be able to reduce number of voxels by a factor of 10-100x by cropping to a single rectangular region of interest, using Crop volume module.

7

Thank you for your advice. This might be a little bit deviated from this post, but am I right in that cropping the volume using crop volume module is not as efficient as directly taking blocks from the numpy array, at least on the programmatic level? Because it is hard to use the mouse to drag the same ROI every time, or even if there is a way to programmatically drag the same ROI every time, then isn’t it the same as directly taking blocks from the numpy array?

8

If you just crop without resampling then Crop volume module will not be much slower. What I propose is just to use a single, clinician-defined region of interest. In that case a slight speed difference is negligible. Depending on the shape and size of the vessel of interest, voxel count may be reduced by 10-100x, so you may not need more sophisticated subdivision.

9

Thank you for all of the help and time! I have selected your previous post as the solution so that others can see it easily.

1 Like
10

Hi Andras,
I need the following results after creating new scalar volume. The input array to update the volume is from fuzzy classification.
when I write using nrrd like below I get the required results.
nrrd.write(dir_path+“/brain_ventricle_membership.nrrd”, brain_ventricle_membership, header)

image
image493×520 44.1 KB

but within a newly created volume I do not get desired results as above.

instead I am getting

image
The section of code I am using is below:
brain_ventricle_membership = np.zeros(brainImg.shape)
brain_ventricle_membership[brainMask > 0] = u[0]

dim = inputVolume.GetImageData().GetDimensions()
spacing = inputVolume.GetSpacing()
outputNode = self.createEmptyVolume(dim, spacing, ‘newVolume’)
slicer.util.updateVolumeFromArray(outputNode, brain_ventricle_membership)
ijkToRAS = vtk.vtkMatrix4x4()
inputVolume.GetIJKToRASMatrix(ijkToRAS)
outputNode.SetIJKToRASMatrix(ijkToRAS)

11

What is this line supposed to do?

brain_ventricle_membership[brainMask > 0] = u[0]

As far as I can tell, you set all the voxels in the mask to a constant value, which is exactly what you show in the result image. Maybe you meant this:

brain_ventricle_membership[brainMask <= 0] = 0

You may also need to set appropriate window/level in the display node of the volume so that you can see the gray levels.

12

the line

brain_ventricle_membership[brainMask > 0] = u[0]

is assigning some values calculated by fuzzy classification to region where the brain ventricles exists. brainMask is nrrd image with binary values, white where the whole brain exists and black with no brain region its like a mask. How can I assign scalar values to the same region. u[0] contains the membership values from fuzzy classification.

how to set appropriate window levels?

13

If you set brain_ventricle_membership the same way for both the first and second example then it should be OK.

See example here: Documentation/Nightly/ScriptRepository - Slicer Wiki