2017-02-20 75 views
0

我正尝试使用numpy/vtk显示使用CT扫描获取的图像。为此,我遵循sample codethis question的回答,但我没有得到好的结果,我不知道原因。用于3D图像渲染和可视化的VTK/numpy

我已经检查出来,我正确加载数据,所以它似乎我做错了渲染时。任何类型的帮助,高度赞赏。这是我的结果到现在为止:提前

This is the 3D image that I obtain

感谢。

这是我的代码:

import os 
import sys 
import pylab 
import glob 

import vtk 
import numpy as np 

#We order all the directories by name 
path="data/Images/" 
tulip_files = [t for t in os.listdir(path)] 
tulip_files.sort() #the os.listdir function do not give the files in the right order so we need to sort them 

#Function that open all the images of a folder and save them in a images list 
def imageread(filePath):  
    filenames = [img for img in glob.glob(filePath)] 
    filenames.sort() 

    temp = pylab.imread(filenames[0]) 
    d, w = temp.shape 
    h = len(filenames) 
    print 'width, depth, height : ',w,d,h 

    volume = np.zeros((w, d, h), dtype=np.uint16) 
    k=0 
    for img in filenames: #assuming tif  
     im=pylab.imread(img) 
     assert im.shape == (500,500), 'Image with an unexpected size' 
     volume[:,:,k] = im 
     k+=1 
    return volume 

#We create the data we want to render. We create a 3D-image by a X-ray CT-scan made to an object. We store the values of each 
#slice and we complete the volume with them in the z axis 
matrix_full = imageread(path+'Image15/raw/reconstruction/*.tif') 

# 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 = matrix_full.tostring() 
dataImporter.CopyImportVoidPointer(data_string, len(data_string)) 
# The type of the newly imported data is set to unsigned short (uint16) 
dataImporter.SetDataScalarTypeToUnsignedShort() 
# Because the data that is imported only contains an intensity value (it isnt RGB-coded or someting 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. 
w, h, d = tulip_matrix_full.shape 
dataImporter.SetDataExtent(0, h-1, 0, d-1, 0, w-1) 
dataImporter.SetWholeExtent(0, h-1, 0, d-1, 0, w-1) 

# This class stores color data and can create color tables from a few color points. 
colorFunc = vtk.vtkPiecewiseFunction() 
colorFunc.AddPoint(0, 0.0); 
colorFunc.AddPoint(65536, 1); 

# The following class is used to store transparency-values for later retrieval. 

alphaChannelFunc = vtk.vtkPiecewiseFunction() 
#Create transfer mapping scalar value to opacity 
alphaChannelFunc.AddPoint(0, 0.0); 
alphaChannelFunc.AddPoint(65536, 1); 

# 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) 
volumeProperty.SetScalarOpacity(alphaChannelFunc) 
#volumeProperty.ShadeOn(); 

# This class describes how the volume is rendered (through ray tracing). 
compositeFunction = vtk.vtkVolumeRayCastCompositeFunction() 
# We can finally create our volume. We also have to specify the data for it, as well as how the data will be rendered. 
volumeMapper = vtk.vtkVolumeRayCastMapper() 
volumeMapper.SetMaximumImageSampleDistance(0.01) # function to reduce the spacing between each image 
volumeMapper.SetVolumeRayCastFunction(compositeFunction) 
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(550, 550) 
renderWin.SetMultiSamples(4) 

# 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() 
+0

你能比“渲染时出现错误”更具体吗?也许你有轴交换? – Benjamin

+0

@Benjamin感谢您的回放。根据上面的图像,我得到了一些与我试图想象的东西非常相似的东西,但只在一个平面上(那种圆圈)。我在“深度”轴上看不到任何东西,这就是为什么我倾向于认为在渲染过程中出现问题的原因。无论如何,我将再次检查轴! –

回答

0

最后,我找到了解决办法。我做了两个重要的变化:

  1. 变化的不透明度值。我有很多接近黑色的体素,所以我修改了不透明度,将它们视为黑色(0.0)。

    alphaChannelFunc.AddPoint(15000,0.0);
    alphaChannelFunc.AddPoint(65536,1);

  2. 更改阵列顺序。看来,在VTK阵列顺序是Fortran order,所以我改变下一个函数来定义正确的轴:

    dataImporter.SetDataExtent(0,H-1,0,d-1,0,W-1)
    dataImporter.SetWholeExtent(0,H-1,0,d-1,0,W-1)

而现在它的作品!