I am doing vesselness analysis on some dataset that requires performing all the tasks from the python interpreter. The steps so far are:
- Splitting the 3D data volume into smaller 3D blocks (In the attached python code I used MRHead sample data as an example)
- 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'))