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

Question

我想把一个三维矩阵(X*Y分辨率的Z片)可视化为经典的三维voxel.我在MATLAB中生成矩阵，并在Python中导入它，然后，按照代码 此处 和 此处，我想出了这个解决方案.通过这个演示，我使用一个矩阵，应该生成一个包含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的文档相矛盾。

*你的数据的尺寸必须等于(extreme)。1-extent[0]+1) * (extent)4-程度3(+1) * (范围5-DataExtent4+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()

我唯一的假设是，这些voxels不是立方体，而是它们的一个尺寸是其他尺寸的两倍。但仍然无法解释为什么4个片子中只有2个受到影响。

更新]：似乎只有第一个和最后一个切片的尺寸是其他切片的一半。在20x30x40的矩阵中，可以看到第一个和最后一个片子比其他片子薄。

Answer 1

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

我的第一个猜测是，数据的存储方式和预期的读取方式之间存在某种不一致。可能是MATLAB将三维矩阵存储为列-主数据结构，而VTK将这类数据恢复为行-主。另一种可能是文件读取时尺寸被调换了，你得到的是2x20x4的实体。