3D圖像可視化與numpy/vtk

Question

我試圖使用numpy/vtk顯示進一步的圖像（ct掃描）如本示例代碼（http://www.vtk.org/Wiki/VTK/Examples/Python/vtkWithNumpy）中所述，但我不明白爲什麼。 如果有人能幫助我，這將是善良的。

這裏是我的代碼：

import vtk 
import numpy as np 
import os 
import cv, cv2 
import matplotlib.pyplot as plt 
import PIL 
import Image 

DEBUG =True 
directory="splitted_mri/" 
w = 226 
h = 186 
d = 27 
stack = np.zeros((w,d,h)) 

k=-1 #add the next picture in a differente level of depth/z-positions 
for file in os.listdir(directory): 
    k+=1 
    img = directory + file 
    im = Image.open(img) 
    temp = np.asarray(im, dtype=int) 
    stack[:,k,:]= temp 
print stack.shape 

#~ plt.imshow(test) 
#~ plt.show() 

print type(stack[10,10,15]) 

res = np.amax(stack) 
res1 = np.amin(stack) 
print res, type(res) 
print res1, type(res1) 

#~ for (x,y,z), value in np.ndenumerate(stack): 
    #~ stack[x,y,z]=np.require(stack[x,y,z],dtype=np.int16) 
    #~ print type(stack[x,y,z]) 

stack = np.require(stack,dtype=np.uint16) 
print stack.dtype 


if DEBUG : print stack.shape 
dataImporter = vtk.vtkImageImport() 
data_string = stack.tostring() 

dataImporter.CopyImportVoidPointer(data_string, len(data_string)) 
dataImporter.SetDataScalarTypeToUnsignedChar() 
dataImporter.SetNumberOfScalarComponents(1) 
dataImporter.SetDataExtent(0, w-1, 0, 1, 0, h-1) 
dataImporter.SetWholeExtent(0, w-1, 0, 1, 0, h-1) 
essai = raw_input() 
alphaChannelFunc = vtk.vtkPiecewiseFunction() 
colorFunc = vtk.vtkColorTransferFunction() 
for i in range (0,255): 
    alphaChannelFunc.AddPoint(i, 0.9) 
    colorFunc.AddRGBPoint(i,i,i,i) 

volumeProperty = vtk.vtkVolumeProperty() 
volumeProperty.SetColor(colorFunc) 
#volumeProperty.ShadeOn() 
volumeProperty.SetScalarOpacity(alphaChannelFunc) 

# 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.SetVolumeRayCastFunction(compositeFunction) 
volumeMapper.SetInputConnection(dataImporter.GetOutputPort()) 

# The class vtkVolume is used to pair the preaviusly 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) 

#to quit, press q 
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() 

Answer 1

我終於發現什麼是錯的 這裏是我的新代碼

import vtk 
import numpy as np 
import os 
import matplotlib.pyplot as plt 
import PIL 
import Image 

DEBUG =False 
directory="splitted_mri/" 

l = [] 

k=0 #add the next picture in a differente level of depth/z-positions 
for file in os.listdir(directory): 
    img = directory + file 
    if DEBUG : print img 
    l.append(img) 
# the os.listdir function do not give the files in the right order 
#so we need to sort them 
l=sorted(l) 

temp = Image.open(l[0]) 
h, w = temp.size 
d = len(l)*5 #with our sample each images will be displayed 5times to get a better view 
if DEBUG : print 'width, height, depth : ',w,h,d 

stack = np.zeros((w,d,h),dtype=np.uint8) 

for i in l: 
    im = Image.open(i) 
    temp = np.asarray(im, dtype=int) 
    for i in range(5): 
     stack[:,k+i,:]= temp 
    k+=5 
    #~ stack[:,k,:]= temp 
    #~ k+=1 

if DEBUG : 
    res = np.amax(stack) 
    print 'max value',res 
    res1 = np.amin(stack) 
print 'min value',res1 

#convert the stack in the right dtype 
stack = np.require(stack,dtype=np.uint8) 

if DEBUG :#check if the image have not been modified 
test = stack [:,0,:] 
plt.imshow(test,cmap='gray') 
plt.show() 

if DEBUG : print 'stack shape & dtype' ,stack.shape,',',stack.dtype 

dataImporter = vtk.vtkImageImport() 
data_string = stack.tostring() 

dataImporter.CopyImportVoidPointer(data_string, len(data_string)) 
dataImporter.SetDataScalarTypeToUnsignedChar() 
dataImporter.SetNumberOfScalarComponents(1) 

#vtk uses an array in the order : height, depth, width which is 
#different of numpy (w,h,d) 
w, d, h = stack.shape 
dataImporter.SetDataExtent(0, h-1, 0, d-1, 0, w-1) 
dataImporter.SetWholeExtent(0, h-1, 0, d-1, 0, w-1) 

alphaChannelFunc = vtk.vtkPiecewiseFunction() 
colorFunc = vtk.vtkColorTransferFunction() 
for i in range(256): 
    alphaChannelFunc.AddPoint(i, 0.2) 
    colorFunc.AddRGBPoint(i,i/255.0,i/255.0,i/255.0) 
# for our test sample, we set the black opacity to 0 (transparent) so as 
#to see the sample 
alphaChannelFunc.AddPoint(0, 0.0) 
colorFunc.AddRGBPoint(0,0,0,0) 

volumeProperty = vtk.vtkVolumeProperty() 
volumeProperty.SetColor(colorFunc) 
#volumeProperty.ShadeOn() 
volumeProperty.SetScalarOpacity(alphaChannelFunc) 

# 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() 
# function to reduce the spacing between each image 
volumeMapper.SetMaximumImageSampleDistance(0.01) 

volumeMapper.SetVolumeRayCastFunction(compositeFunction) 
volumeMapper.SetInputConnection(dataImporter.GetOutputPort()) 

# The class vtkVolume is used to pair the preaviusly 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) 

#to auit, press q 
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()