VTK: 3D matrix's volume has wrong dimensions

Question

I want to visualize a 3-D matrix (Z slices of X*Y resolution) as the classic 3-d voxel. I generate the matrix in MATLAB and import it in Python. Then, following the code here and here, I come up with this solution. With this demonstration I am using a matrix that should generate a 3D image containing 4 slices of 2*3 voxels.

In MATLAB

C(:,:,1) =

   5   5   5
   5   5   5


C(:,:,2) =

   15   15   15
   15   15   15


C(:,:,3) =

   25   25   25
   25   25   25


C(:,:,4) =

   35   35   35
   35   35   35

In python:

Cmat = spio.loadmat('CMAT.mat')['C']
>>>print Cmat.shape
(2, 3, 4)

Cmat = np.ascontiguousarray(Cmat.T)
>>>print Cmat
[[[ 5  5]
  [ 5  5]
  [ 5  5]]

 [[15 15]
  [15 15]
  [15 15]]


 [[25 25]
  [25 25]
  [25 25]]

 [[35 35]
  [35 35]
  [35 35]]]

And the code that will follow produces this image (which I rotated for convenience):

The resulting shape is not a 2*3*4 and the slices have not the same size, what am I doing wrong? I tried to tweak

dataImporter.SetDataExtent
dataImporter.SetWholeExtent
dataImporter.SetDataSpacing

And reshaping the matrix in many ways, if I change dataImporter.SetDataExtent(0,1,0,1,0,1) dataImporter.SetWholeExtent(0,1,0,1,0,1)

I obtain a 2x2x2 cube as expected but if I call

dataImporter.SetDataExtent(0, 1, 0, 2, 0, 1)
dataImporter.SetWholeExtent(0, 1, 0, 2, 0, 1)

I obtain a 2x4x2 solid (instead of 2x3x2)

if I call

dataImporter.SetDataExtent(0, 1, 0, 10, 0, 2)
dataImporter.SetWholeExtent(0, 1, 0, 10, 0, 2)

I obtain a 2x20x4 solid

This appears to contradict the documentation for setDataExtent and SetWholeExtent:

*The dimensions of your data must be equal to (extent1-extent[0]+1) * (extent4-extent3+1) * (extent5-DataExtent4+1). For example, for a 2D image use (0,width-1, 0,height-1, 0,0).*

Any Idea?

Full Code below MATLAB:

C = zeros(2,3,4)
C(:,:,1) = 5;
C(:,:,2) = 15;
C(:,:,3) = 25;
C(:,:,4) = 35;

save Cmat C

Python:

import vtk
from numpy import *
import numpy as np
import scipy.io as spio

data_matrix = zeros([2, 3, 4], dtype=uint8)

Cmat = spio.loadmat('CMAT.mat')['C']
Cmat = np.ascontiguousarray(Cmat.T)
print Cmat
data_matrix = Cmat
# For VTK to be able to use the data, it must be stored as a VTK-image. This can be done by the vtkImageImport-class which
# imports raw data and stores it.
dataImporter = vtk.vtkImageImport()
# The previously created array is converted to a string of chars and imported.
data_string = data_matrix.tostring()
dataImporter.CopyImportVoidPointer(data_string, len(data_string))
# The type of the newly imported data is set to unsigned char (uint8)
dataImporter.SetDataScalarTypeToUnsignedChar()
# Because the data that is imported only contains an intensity value (it isn't RGB-coded or something similar), the importer
# must be told this is the case.
dataImporter.SetNumberOfScalarComponents(1)
# The following two functions describe how the data is stored and the dimensions of the array it is stored in.
dataImporter.SetDataExtent(0, 1, 0, 2, 0, 3)
dataImporter.SetWholeExtent(0, 1, 0, 2, 0, 3)

# This class stores color data and can create color tables from a few color points. For this demo, we want the three cubes
# to be of the colors red green and blue.
colorFunc = vtk.vtkColorTransferFunction()
colorFunc.AddRGBPoint(5, 1, 0.0, 0.0)  # Red
colorFunc.AddRGBPoint(15, 0.0, 1, 0.0) # Green
colorFunc.AddRGBPoint(25, 0.0, 0.0, 1) # Blue
colorFunc.AddRGBPoint(35, 0.0, 0, 0.0) # Black

# The previous two classes stored properties. Because we want to apply these properties to the volume we want to render,
# we have to store them in a class that stores volume properties.
volumeProperty = vtk.vtkVolumeProperty()
volumeProperty.SetColor(colorFunc)

volumeMapper = vtk.vtkOpenGLGPUVolumeRayCastMapper()
volumeMapper.SetInputConnection(dataImporter.GetOutputPort())

# The class vtkVolume is used to pair the previously declared volume as well as the properties to be used when rendering that volume.
volume = vtk.vtkVolume()
volume.SetMapper(volumeMapper)
volume.SetProperty(volumeProperty)

# With almost everything else ready, its time to initialize the renderer and window, as well as creating a method for exiting the application
renderer = vtk.vtkRenderer()
renderWin = vtk.vtkRenderWindow()
renderWin.AddRenderer(renderer)
renderInteractor = vtk.vtkRenderWindowInteractor()
renderInteractor.SetRenderWindow(renderWin)

# We add the volume to the renderer ...
renderer.AddVolume(volume)
# ... set background color to white ...
renderer.SetBackground(1, 1, 1)
# ... and set window size.
renderWin.SetSize(400, 400)


# A simple function to be called when the user decides to quit the application.
def exitCheck(obj, event):
    if obj.GetEventPending() != 0:
        obj.SetAbortRender(1)


# Tell the application to use the function as an exit check.
renderWin.AddObserver("AbortCheckEvent", exitCheck)

renderInteractor.Initialize()
# Because nothing will be rendered without any input, we order the first render manually before control is handed over to the main-loop.
renderWin.Render()
renderInteractor.Start()

The only hypotesis I have is that these voxels are not cubes but one of their dimensions is twice the others. But still would not explain why only 2 of the 4 slices are impacted by this.

[UPDATE]: It seems that only the first and last slices are half the size of the other ones. With a 20x30x40 matrix it can be seen that the first and last slices are thinner than the rest.

Answer 1

This is an old question, probably you have found out the answer.

My first guess is that there is some kind of inconsistency between the way data is stored and how is expected to be read. Probably, MATLAB is storing the 3D matrix as column-major data structure, and VTK is recovering such data as row-major. Another possibility would be that dimensions are swapped when file is read and you obtained the 2x20x4 solid.