在Python的3D圖像進行中，通常用numpy array來進行非常友善的計算或者轉化，這裡記錄一下numpy資料的VTK可視化基本流程，包括面繪制（Surfase Rendering）和體繪制（Volume Rendering）。除去資料格式轉化，面繪制和體繪制在C++中也是類似的處理方法。
numpy資料轉成 vtkImageData
vtkImageData
首先得把numpy資料轉成vtk裡可以用的格式：numpy array -> vtkIImageData。這裡的numpy array是一個離散的三維空間資料場，0代表背景，非0代表前景點。
方法一： numpy_to_vtk
numpy_to_vtk
可以直接使用
vtk.util.numpy_support
裡的轉換方法，非常友善。這是它的說明文檔：
numpy_to_vtk(num_array, deep=0, array_type=None)
Converts a contiguous real numpy Array to a VTK array object.
This function only works for real arrays that are contiguous.
Complex arrays are NOT handled. It also works for multi-component
arrays. However, only 1, and 2 dimensional arrays are supported.
This function is very efficient, so large arrays should not be a
problem.
If the second argument is set to 1, the array is deep-copied from
from numpy. This is not as efficient as the default behavior
(shallow copy) and uses more memory but detaches the two arrays
such that the numpy array can be released.
WARNING: You must maintain a reference to the passed numpy array, if
the numpy data is gc'd and VTK will point to garbage which will in
the best case give you a segfault.
Parameters
----------
- num_array : a contiguous 1D or 2D, real numpy array.
對于3D array，可以用
flatten
或者
ravel
先轉成1D array就可以使用了。這裡的1D array得是C order（row-major order），最好使用deep copy以免出現一些記憶體管理的問題。
import numpy as np
import vtk
from vtk.util import numpy_support
# numpy_data is a 3D numpy array
shape = numpy_data.shape[::-1]
vtk_data = numpy_support.numpy_to_vtk(numpy_data.ravel(), 1, vtk.VTK_SHORT)
vtk_image_data = vtk.vtkImageData()
vtk_image_data.SetDimensions(shape)
vtk_image_data.SetSpacing(spacing)
vtk_image_data.SetOrigin(origin)
vtk_image_data.GetPointData().SetScalars(vtk_data)
# vtk_image_data: ready to use
方法二： vtkImageImport
vtkImageImport
使用Python也可以用
vtkImageImport
來直接導入C array，就是要先把numpy array變成string。
import numpy as np
import vtk
numpy_str = numpy_data.astype(np.int16).tostring()
x_extent = numpy_data.shape[2]
y_extent = numpy_data.shape[1]
z_extent = numpy_data.shape[0]
image_import = vtk.vtkImageImport()
image_import.SetImportVoidPointer(numpy_str, len(numpy_str))
# 也可以使用CopyImportVoidPointer() 會copy一份numpy_str
image_import.SetWholeExtent(0, x_extent-1, 0, y_extent-1, 0, z_extent-1)
image_import.SetDataExtent(0, x_extent-1, 0, y_extent-1, 0, z_extent-1)
image_import.SetDataScalarTypeToShort() # 根據需求指定資料類型
image_import.SetNumberOfScalarComponents(1)
# 如果是RGB資料的話，SetNumberOfScalarComponents(3)
image_import.Update()
vtk_image_data = vtk.vtkImageData()
vtk_image_data.SetSpacing(spacing)
vtk_image_data.SetOrigin(origin)
# vtk_image_data: ready to use
繪制
參考：https://www.cnblogs.com/XDU-Lakers/p/10822840.html
這個作者總結的很全面，還有代碼示例👍
面繪制
面繪制的意思是根據資料建立三角網格模型，再渲染網格。通俗來講就是根據算法進行輪廓識别和提取，生成了三維物體的表面，這個表面本質是由非常多的小三角面構成。
如何抽取物體表面（輪廓/等值面）？VTK提供了一些算法（filter）：
-
vtkMarchingCubes
- 必須指定一個或多個輪廓值
- 2D版本是
vtkMarchingSquares
- 多标簽圖像是
vtkDiscreteMarchingCubdes
vtkThreshold
-
vtkContourFilter
- 3D生成輪廓面，2D生成輪廓線，1D生成輪廓點
- 必須指定一個或多個輪廓值
- 不會自動生成normal，可以用
vtkPolyDataNormals
-
vtkFlyingEdges3D
- 提速版輪廓提取算法，适用于較大的資料
- 2D版本是
vtkFlyingEdges2D
什麼是輪廓值？其實就是根據這個值來提取輪廓。比如對于一個二值圖像，所有值為1的點是前景點（要顯示的物體），所有值為0的點是背景點。此時輪廓值是1，算法會找到值為1的點所構成物體的輪廓。如果是多标簽圖像，可以把多個标簽的值當作輪廓值。如果是其它醫學圖像，可以根據需求取不同的值，比如人體皮膚所對應的value值為500，人體骨骼所對應的value值為1150。設定多個輪廓值，這樣多個等值面都可以被提取出來。
經過這些filter，資料就被處理成了表現三維物體表面的
vtkPolyData
，再經過
vtkPolyDataMapper
，
vtkActor
、
vtkRenderer
、
vtkRenderWindow
、
vtkRenderWindowInteractor
顯示出來。
官方使用範例：(多标簽資料 面繪制 上色 / multilabel data)
# https://kitware.github.io/vtk-examples/site/Python/Modelling/SmoothDiscreteMarchingCubes/
import vtk
def main():
n = 20
radius = 8
blob = make_blob(n, radius)
discrete = vtk.vtkDiscreteMarchingCubes()
discrete.SetInputData(blob)
discrete.GenerateValues(n, 1, n)
smoothing_iterations = 15
pass_band = 0.001
feature_angle = 120.0
smoother = vtk.vtkWindowedSincPolyDataFilter()
smoother.SetInputConnection(discrete.GetOutputPort())
smoother.SetNumberOfIterations(smoothing_iterations)
smoother.BoundarySmoothingOff()
smoother.FeatureEdgeSmoothingOff()
smoother.SetFeatureAngle(feature_angle)
smoother.SetPassBand(pass_band)
smoother.NonManifoldSmoothingOn()
smoother.NormalizeCoordinatesOn()
smoother.Update()
lut = make_colors(n)
mapper = vtk.vtkPolyDataMapper()
mapper.SetInputConnection(smoother.GetOutputPort())
mapper.SetLookupTable(lut)
mapper.SetScalarRange(0, lut.GetNumberOfColors())
# Create the RenderWindow, Renderer and both Actors
#
ren = vtk.vtkRenderer()
ren_win = vtk.vtkRenderWindow()
ren_win.AddRenderer(ren)
ren_win.SetWindowName('SmoothDiscreteMarchingCubes')
iren = vtk.vtkRenderWindowInteractor()
iren.SetRenderWindow(ren_win)
actor = vtk.vtkActor()
actor.SetMapper(mapper)
ren.AddActor(actor)
colors = vtk.vtkNamedColors()
ren.SetBackground(colors.GetColor3d('Burlywood'))
ren_win.Render()
iren.Start()
def make_blob(n, radius):
blob_image = vtk.vtkImageData()
max_r = 50 - 2.0 * radius
random_sequence = vtk.vtkMinimalStandardRandomSequence()
random_sequence.SetSeed(5071)
for i in range(0, n):
sphere = vtk.vtkSphere()
sphere.SetRadius(radius)
x = random_sequence.GetRangeValue(-max_r, max_r)
random_sequence.Next()
y = random_sequence.GetRangeValue(-max_r, max_r)
random_sequence.Next()
z = random_sequence.GetRangeValue(-max_r, max_r)
random_sequence.Next()
sphere.SetCenter(int(x), int(y), int(z))
sampler = vtk.vtkSampleFunction()
sampler.SetImplicitFunction(sphere)
sampler.SetOutputScalarTypeToFloat()
sampler.SetSampleDimensions(100, 100, 100)
sampler.SetModelBounds(-50, 50, -50, 50, -50, 50)
thres = vtk.vtkImageThreshold()
thres.SetInputConnection(sampler.GetOutputPort())
thres.ThresholdByLower(radius * radius)
thres.ReplaceInOn()
thres.ReplaceOutOn()
thres.SetInValue(i + 1)
thres.SetOutValue(0)
thres.Update()
if i == 0:
blob_image.DeepCopy(thres.GetOutput())
max_value = vtk.vtkImageMathematics()
max_value.SetInputData(0, blob_image)
max_value.SetInputData(1, thres.GetOutput())
max_value.SetOperationToMax()
max_value.Modified()
max_value.Update()
blob_image.DeepCopy(max_value.GetOutput())
return blob_image
def make_colors(n):
"""
Generate some random colors
:param n: The number of colors.
:return: The lookup table.
"""
lut = vtk.vtkLookupTable()
lut.SetNumberOfColors(n)
lut.SetTableRange(0, n - 1)
lut.SetScaleToLinear()
lut.Build()
lut.SetTableValue(0, 0, 0, 0, 1)
random_sequence = vtk.vtkMinimalStandardRandomSequence()
random_sequence.SetSeed(5071)
for i in range(1, n):
r = random_sequence.GetRangeValue(0.4, 1)
random_sequence.Next()
g = random_sequence.GetRangeValue(0.4, 1)
random_sequence.Next()
b = random_sequence.GetRangeValue(0.4, 1)
random_sequence.Next()
lut.SetTableValue(i, r, g, b, 1.0)
return lut
if __name__ == '__main__':
main()
體繪制
體繪制是為每一個體素指定一個不透明度，并考慮每一個體素對光線的透射、發射和反射作用，實作三維重建。簡單來說就是能夠更完整的展示出整個物體，而不僅僅是表面。
VTK采用的是光線投射算法。光線投射算法（Ray-casting）原理：從圖像平面的每個像素都沿着視線方向發出一條射線，此射線穿過體資料集，按一定步長進行采樣，由内插計算每個采樣點的顔色值和不透明度，然後由前向後或由後向前逐點計算累計的顔色值和不透明度值，直至光線完全被吸收或穿過物體。該方法能很好地反映物質邊界的變化，使用Phong模型，引入鏡面反射、漫反射和環境反射能得到很好的光照效果，在醫學上可将各組織器官的性質屬性、形狀特征及互相之間的層次關系表現出來，進而豐富了圖像的資訊。（百度百科）
這個算法內建在
vtkVolumeMapper
，包括
vtkVolumeRayCastMapper
、
vtkFixedPointVolumeRayCastMapper
、
vtkGPUVolumeRayCastMapper
等。把
vtkImageData
經過volumeMapper，再放入
vtkVolume
，經過
vtkRenderer
、
vtkRenderWindow
、
vtkRenderWindowInteractor
顯示出體繪制圖像。
官方使用範例：
# https://kitware.github.io/vtk-examples/site/Python/VolumeRendering/SimpleRayCast/
import vtk
def main():
fileName = get_program_parameters()
colors = vtk.vtkNamedColors()
# This is a simple volume rendering example that
# uses a vtkFixedPointVolumeRayCastMapper
# Create the standard renderer, render window
# and interactor.
ren1 = vtk.vtkRenderer()
renWin = vtk.vtkRenderWindow()
renWin.AddRenderer(ren1)
iren = vtk.vtkRenderWindowInteractor()
iren.SetRenderWindow(renWin)
# Create the reader for the data.
reader = vtk.vtkStructuredPointsReader()
reader.SetFileName(fileName)
# Create transfer mapping scalar value to opacity.
opacityTransferFunction = vtk.vtkPiecewiseFunction()
opacityTransferFunction.AddPoint(20, 0.0)
opacityTransferFunction.AddPoint(255, 0.2)
# Create transfer mapping scalar value to color.
colorTransferFunction = vtk.vtkColorTransferFunction()
colorTransferFunction.AddRGBPoint(0.0, 0.0, 0.0, 0.0)
colorTransferFunction.AddRGBPoint(64.0, 1.0, 0.0, 0.0)
colorTransferFunction.AddRGBPoint(128.0, 0.0, 0.0, 1.0)
colorTransferFunction.AddRGBPoint(192.0, 0.0, 1.0, 0.0)
colorTransferFunction.AddRGBPoint(255.0, 0.0, 0.2, 0.0)
# The property describes how the data will look.
volumeProperty = vtk.vtkVolumeProperty()
volumeProperty.SetColor(colorTransferFunction)
volumeProperty.SetScalarOpacity(opacityTransferFunction)
volumeProperty.ShadeOn()
volumeProperty.SetInterpolationTypeToLinear()
# The mapper / ray cast function know how to render the data.
volumeMapper = vtk.vtkFixedPointVolumeRayCastMapper()
volumeMapper.SetInputConnection(reader.GetOutputPort())
# The volume holds the mapper and the property and
# can be used to position/orient the volume.
volume = vtk.vtkVolume()
volume.SetMapper(volumeMapper)
volume.SetProperty(volumeProperty)
ren1.AddVolume(volume)
ren1.SetBackground(colors.GetColor3d('Wheat'))
ren1.GetActiveCamera().Azimuth(45)
ren1.GetActiveCamera().Elevation(30)
ren1.ResetCameraClippingRange()
ren1.ResetCamera()
renWin.SetSize(600, 600)
renWin.SetWindowName('SimpleRayCast')
renWin.Render()
iren.Start()
def get_program_parameters():
import argparse
description = 'Volume rendering of a high potential iron protein.'
epilogue = '''
This is a simple volume rendering example that uses a vtkFixedPointVolumeRayCastMapper.
'''
parser = argparse.ArgumentParser(description=description, epilog=epilogue,
formatter_class=argparse.RawDescriptionHelpFormatter)
parser.add_argument('filename', help='ironProt.vtk.')
args = parser.parse_args()
return args.filename
if __name__ == '__main__':
main()
![TestLink導出用例轉換工具(XML2Excel)[圖]](data:image/gif;base64,R0lGODlhAQABAIAAAP///wAAACwAAAAAAQABAAACAkQBADs=)