VTK：3D 矩阵的体积尺寸错误

Question

我想将 3-D 矩阵（X*Y 分辨率的 Z 切片）可视化为经典的 3-d 体素。 我在 MATLAB 中生成矩阵并在 Python 中导入它。 然后，按照代码here 和here，我想出了这个解决方案。 在这个演示中，我使用的矩阵应该生成一个包含 4 个 2*3 体素切片的 3D 图像。

在 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

在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]]]

接下来的代码会生成这个图像（为方便起见，我旋转了它）：

生成的形状不是 2*3*4 并且切片的大小不同，我做错了什么？我试图调整

dataImporter.SetDataExtent
dataImporter.SetWholeExtent
dataImporter.SetDataSpacing

如果我改变的话，还会以多种方式重塑矩阵 dataImporter.SetDataExtent(0,1,0,1,0,1) dataImporter.SetWholeExtent(0,1,0,1,0,1)

我按预期获得了一个 2x2x2 的立方体 但是如果我打电话

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

我得到一个 2x4x2 的实体（而不是 2x3x2）

如果我打电话

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

我得到一个 2x20x4 的实体

这似乎与 setDataExtent 和 SetWholeExtent 的文档相矛盾：

*数据的维度必须等于 (extent1-extent[0]+1) * (extent4-extent3+1) * (extent5-DataExtent@987654329 @+1)。例如，对于 2D 图像，使用 (0,width-1, 0,height-1, 0,0)。*

有什么想法吗？

完整代码如下 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()

我唯一的假设是这些体素不是立方体，但它们的一个维度是其他维度的两倍。但仍然无法解释为什么 4 个切片中只有 2 个受此影响。

[更新]： 似乎只有第一个和最后一个切片的大小是其他切片的一半。使用 20x30x40 矩阵可以看出，第一个和最后一个切片比其他切片更薄。

Answer 1

这是一个老问题，可能你已经找到答案了。

我的第一个猜测是，数据的存储方式与预期的读取方式之间存在某种不一致。可能，MATLAB 将 3D 矩阵存储为以列为主的数据结构，而 VTK 将此类数据恢复为以行为主的数据。 另一种可能性是在读取文件并获得 2x20x4 实体时交换尺寸。