VTK: объем 3D-матрицы имеет неправильные размеры

Question

Я хочу визуализировать трехмерную матрицу (срезы Z с разрешением X*Y) как классический трехмерный воксель. Я генерирую матрицу в MATLAB и импортирую ее в Python. Затем, следуя коду здесь и здесь, я придумываю это решение. В этой демонстрации я использую матрицу, которая должна генерировать трехмерное изображение, содержащее 4 среза 2*3 вокселей.

В 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

В питоне:

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:

* Размеры ваших данных должны быть равны (экстент 1- экстент [0]+1) * (экстент 4- экстент 3+1) * (экстент 5-DataExtent 4+1). Например, для 2D-изображения используйте (0, ширина-1, 0, высота-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()

Единственный гипотет, который у меня есть, состоит в том, что эти воксели не кубы, а одно из их размеров вдвое больше других. Но все равно не объяснил бы, почему это влияет только на 2 из 4 срезов.

[ОБНОВЛЕНИЕ]: Кажется, что только первый и последний кусочки в два раза меньше остальных. С матрицей 20x30x40 видно, что первый и последний срезы тоньше остальных.

Answer 1

Это старый вопрос, наверное, вы уже нашли ответ.

Мое первое предположение заключается в том, что существует какое-то несоответствие между способом хранения данных и тем, как ожидается их чтение. Вероятно, MATLAB хранит трехмерную матрицу как структуру данных основного столбца, а VTK восстанавливает такие данные как структуру данных основного столбца. Другая возможность заключается в том, что размеры меняются местами при чтении файла и получении твердого тела 2x20x4.