⽤VTK-Python实现nii格式腹部CT三维重建与⼈机交互
python-vtk做医学nii格式的三维重建
课题的开始
导师定了课题《医学影像CT的三维重建》，因为实验室的主流技术还是深度学习做医学影像和⾃然图像的处理（2D），从来没有做过三维的东西。抱着试⼀试的⼼态，开始查阅相关⽂献。⾸先说明，主流的三维重建问题分享中，多视⾓三维重建和多层⾯三维重建。前者input⼏张某物体个⾓度的jpg，output该物体的三维图像，⽬前计算机视觉顶会论⽂多有涉及，不管是单⽬还是多⽬的三维重建，⽤深度学习⽣成三维模型即为输出；后者input⼀个⼆维序列图像，output该物体的三维图像，接下来主要说说多层⾯三维重建
关于多层⾯三维重建（医学影像+⼯业CT）
了解到的相关论⽂分两种
第⼀种 打着三维重建的题⽬，去做插值，分割某⼀器官，或者融合⼆维图像（其实就是为三维重建⽽准备的，对⼆维序列数据的预处理。插值即是使数据饱满，分割即是使数据分类，融合则适⽤于MRI多模态融合在三维重建），最后调⽤VTK的Ray Casting 去完成⼀个三维重建的效果。当然，仅做到这⾥你
会得到⼀个根本没法⽤的三维可视化效果，需要调整颜⾊，透明度，窗宽窗位，然后你就会得到⼀个看起来像那么⼀回事的三维模型
第⼆种 1980s~1990s的经典体绘制⽅法论⽂（个⼈感觉：光线投影算法作为体绘制算法的⼀种，⽬前来说是效果最好的直接体绘制算法之⼀，只是20世纪末受限于GPU运算速度，在现在来说已经不太是⼀个问题）
所以，我只得先去使⽤vtk-python实现对⾃⼰⼿上腹部nii数据的三维重建，继续往下
使⽤vtk-Python完成腹部nii数据的三维重建与⼈机交互
// 3d reconstruction ud by vtk
from vtk.util.vtkImageImportFromArray import*
import vtk
import SimpleITK as sitk
import numpy as np
import time
# path ='../vtk/nii_data_low/1_1.nii' #gmentation volume
path ='../vtk/volume_2.nii' #gmentation volume糖醋
ds = sitk.ReadImage(path) #读取nii数据的第⼀个函数sitk.ReadImage
print('ds: ',ds)
data = sitk.GetArrayFromImage(ds) #把itk.image转为array
print('data: ',data)
print('shape_of_data',data.shape)
# 去掉Hu值⼩于x的点
# time_start=time.time()
# sum =0
# for i,iindex in enumerate(data):
# for j,jindex in enumerate(data[i]):
# for k,kindex in enumerate(data[j]):
# sum = sum +1
# if data[i][j][k]<1129:
# data[i][j][k]=-1024
# time_end=time.time()
# print('time cost',time_end-time_start,'s')
# sum
spacing = ds.GetSpacing() #三维数据的间隔
print('spacing_of_data',spacing)
# data = data[50:]
# data = data[:,:,300:]
srange =[np.min(data),np.max(data)]
print('shape_of_data_chenged',data.shape)
img_arr =vtkImageImportFromArray() #创建⼀个空的vtk类-----vtkImageImportFromArray
img_arr =vtkImageImportFromArray() #创建⼀个空的vtk类-----vtkImageImportFromArray
print('img_arr: ',img_arr)
img_arr.SetArray(data) #把array_data塞到vtkImageImportFromArray（array_data）img_arr.SetDataSpacing(spacing) #设置spacing
origin =(0,0,0)
img_arr.SetDataOrigin(origin) #设置vtk数据的坐标系原点
img_arr.Update()
# srange = img_arr.GetOutput().GetScalarRange()
print('spacing: ',spacing)
print('srange: ',srange)
# 键盘控制交互式操作
class KeyPressInteractorStyle(vtk.vtkInteractorStyleTrackballCamera):
def __init__(lf,parent=None):韩语常用语
lf.parent = vtk.vtkRenderWindowInteractor()
if(parent is not None):
lf.parent = parent
lf.AddObrver("KeyPressEvent",lf.keyPress)
如何计算圆的面积def keyPress(lf,obj,event):
key = lf.parent.GetKeySym()
if key =='Up':
gradtfun.AddPoint(-100,1.0)
gradtfun.AddPoint(10,1.0)
gradtfun.AddPoint(20,1.0)
三爱三节volumeProperty.SetGradientOpacity(gradtfun)
#下⾯这⼀⾏是关键，实现了actor的更新
renWin.Render()
if key =='Down':
tfun.AddPoint(1129,0)
tfun.AddPoint(1300.0,0.1)
tfun.AddPoint(1600.0,0.2)
tfun.AddPoint(2000.0,0.1)
tfun.AddPoint(2200.0,0.1)
化妆课tfun.AddPoint(2500.0,0.1)
tfun.AddPoint(2800.0,0.1)
tfun.AddPoint(3000.0,0.1)
#下⾯这⼀⾏是关键，实现了actor的更新
renWin.Render()
def StartInteraction():
renWin.SetDesiredUpdateRate(10)
def EndInteraction():
renWin.SetDesiredUpdateRate(0.001)
def ClipVolumeRender(obj):
obj.GetPlanes(planes)
volumeMapper.SetClippingPlanes(planes)
ren = vtk.vtkRenderer()
renWin= vtk.vtkRenderWindow()
renWin.AddRenderer(ren) #把⼀个空的渲染器添加到⼀个空的窗⼝上
renWin.AddRenderer(ren)
iren = vtk.vtkRenderWindowInteractor()
iren.SetRenderWindow(renWin) #把上⾯那个窗⼝加⼊交互操作
有关青春的素材iren.SetRenderWindow(renWin) #把上⾯那个窗⼝加⼊交互操作
iren.SetInteractorStyle(KeyPressInteractorStyle(parent = iren)) #在交互操作⾥⾯添加这个⾃定义的操作例如up,down min = srange[0]
max = srange[1]
# diff = max - min #体数据极差
# slope =4000/ diff
# inter =-slope * min
# shift = inter / slope
# print(min, max, slope, inter, shift) #这⼏个数据后⾯有⽤
diff = max - min #体数据极差
inter =4200/ diff团员学习心得
shift =-min
print(min, max, inter, shift) #这⼏个数据后⾯有⽤
# diffusion = vtk.vtkImageAnisotropicDiffusion3D()
# diffusion.SetInputData(img_arr.GetOutput())
# diffusion.SetNumberOfIterations(10)
# diffusion.SetDiffusionThreshold(5)
# diffusion.Update()
shifter = vtk.vtkImageShiftScale() # 对偏移和⽐例参数来对图像数据进⾏操作数据转换，之后直接调⽤shifter shifter.SetShift(shift)
shifter.SetScale(inter)
shifter.SetOutputScalarTypeToUnsignedShort()
shifter.SetInputData(img_arr.GetOutput())
shifter.ReleaDataFlagOff()
shifter.Update()
tfun = vtk.vtkPiecewiFunction() # 不透明度传输函数---放在tfun
tfun.AddPoint(1129,0)
tfun.AddPoint(1300.0,0.1)
tfun.AddPoint(1600.0,0.12)
tfun.AddPoint(2000.0,0.13)
tfun.AddPoint(2200.0,0.14)
tfun.AddPoint(2500.0,0.16)
tfun.AddPoint(2800.0,0.17)
tfun.AddPoint(3000.0,0.18)
gradtfun = vtk.vtkPiecewiFunction() # 梯度不透明度函数---放在gradtfun
gradtfun.AddPoint(-1000,9)
gradtfun.AddPoint(0.5,9.9)
gradtfun.AddPoint(1,10)
ctfun = vtk.vtkColorTransferFunction() # 颜⾊传输函数---放在ctfun
ctfun.AddRGBPoint(0.0,0.5,0.0,0.0)
ctfun.AddRGBPoint(600.0,1.0,0.5,0.5)
ctfun.AddRGBPoint(1280.0,0.9,0.2,0.3)
ctfun.AddRGBPoint(1960.0,0.81,0.27,0.1)
ctfun.AddRGBPoint(2200.0,0.9,0.2,0.3)
ctfun.AddRGBPoint(2500.0,1,0.5,0.5)
ctfun.AddRGBPoint(3024.0,0.5,0.5,0.5)
# ctfun.AddRGBPoint(0.0,0.5,0.0,0.0)
# ctfun.AddRGBPoint(600.0,1.0,255,0.5)
# ctfun.AddRGBPoint(1280.0,0.9,0.2,255)回忆过去美好的句子
# ctfun.AddRGBPoint(1960.0,255,0.27,0.1)
# ctfun.AddRGBPoint(2200.0,0.9,0.2,0.3)
# ctfun.AddRGBPoint(2500.0,1,0.5,0.5)
# ctfun.AddRGBPoint(3024.0,0.5,0.5,0.5)
volumeMapper = vtk.vtkGPUVolumeRayCastMapper() #映射器volumnMapper使⽤vtk的管线投影算法volumeMapper.SetInputData(shifter.GetOutput()) #向映射器中输⼊数据：shifter(预处理之后的数据) volumeProperty = vtk.vtkVolumeProperty() #创建vtk属性存放器，向属性存放器中存放颜⾊和透明度volumeProperty.SetColor(ctfun)
volumeProperty.SetScalarOpacity(tfun)
# volumeProperty.SetGradientOpacity(gradtfun)
# volumeProperty.SetGradientOpacity(gradtfun) volumeProperty.SetInterpolationTypeToLinear() # volumeProperty.ShadeOn()
newvol = vtk.vtkVolume() #演员
newvol.SetMapper(volumeMapper)
newvol.SetProperty(volumeProperty)
outline = vtk.vtkOutlineFilter()
outline.SetInputConnection(shifter.GetOutputPort())
outlineMapper = vtk.vtkPolyDataMapper()
outlineMapper.SetInputConnection(outline.GetOutputPort())
outlineActor = vtk.vtkActor()
outlineActor.SetMapper(outlineMapper)
ren.AddActor(outlineActor)
ren.AddVolume(newvol)
ren.SetBackground(0,0,0)
renWin.SetSize(600,600)
planes = vtk.vtkPlanes()
boxWidget = vtk.vtkBoxWidget()
boxWidget.SetInteractor(iren)
boxWidget.SetPlaceFactor(1.0)
boxWidget.PlaceWidget(0,0,0,0,0,0)
boxWidget.InsideOutOn()
boxWidget.AddObrver("StartInteractionEvent", StartInteraction) boxWidget.AddObrver("InteractionEvent", ClipVolumeRender) boxWidget.AddObrver("EndInteractionEvent", EndInteraction)
outlineProperty = boxWidget.GetOutlineProperty() outlineProperty.SetReprentationToWireframe() outlineProperty.SetAmbient(1.0)
outlineProperty.SetAmbientColor(1,1,1)
outlineProperty.SetLineWidth(9)
lectedOutlineProperty = boxWidget.GetSelectedOutlineProperty() lectedOutlineProperty.SetReprentationToWireframe() lectedOutlineProperty.SetAmbient(1.0) lectedOutlineProperty.SetAmbientColor(1,0,0) lectedOutlineProperty.SetLineWidth(3)
ren.RetCamera()
iren.Initialize()
renWin.Render()
iren.Start()
记录⾃⼰的学习过程， 如有问题，请多指正
