3d与numpy / vtk的图象形象化(3d image visualisation with numpy/vtk)

3d与numpy / vtk的图象形象化(3d image visualisation with numpy/vtk)

我正在尝试使用numpy / vtk显示更多图像（ct-scan），如本示例代码中所述（ 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()

I'm trying to display further images (ct-scan) using numpy/vtk as describe in this sample code (http://www.vtk.org/Wiki/VTK/Examples/Python/vtkWithNumpy) but I don't get it and don't know why. If someone could help me it would be kind.

Here's my code :

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()

最满意答案

我终于发现这里的错误是我的新代码

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()

I finally find out what was wrong here's my new code

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()