效果图
废话少说先上效果图,右边红色的plane是想要获取3D模型对应切面的平面,左边是这个切面的切片的图像。
实现目标
定义一个任意角度的切面,都能把体绘制模型的这个切面的图像获取,并且能够把这个切面图像转为numpy格式供其他逻辑继续处理。MPR三维重建只能获取xyz轴的切面无法满足任意角度切片的需求。
安装依赖
pip install vtk numpy opencv-python scipy
Code
体渲染的知识可以参考我上一篇文章:Python使用VTK对容积超声图像进行体绘制(三维重建)
import vtk
import numpy as np
import cv2
from vtkmodules.util.numpy_support import numpy_to_vtk, vtk_to_numpy
from scipy.ndimage import rotate # 3D旋转
# 定义原点vtk对象原点位置
origin = [0, 0, 0]
# 定义xyz旋转角度
x_angle = 0
y_angle = 0
z_angle = 0
# 导入3D模型, 根据自己场景调整, 最终类型是vtkImageData就可以
label_arr = np.load("./your_model.npy")
label_arr = rotate(label_arr, 180, axes=(0, 1), reshape=True, order=1, mode='constant', cval=0, prefilter=True)
label_arr = rotate(label_arr, -90, axes=(0, 2), reshape=True, order=1, mode='constant', cval=0, prefilter=True)
# numpy转为vtkTypeUInt8Array
vtk_data = numpy_to_vtk(label_arr.ravel(), array_type=vtk.VTK_UNSIGNED_CHAR)
# 获取模型尺寸
model_dims = label_arr.shape
# 定义一个vtkImageData数据对象, 作为vtkImageReslice的输入
image_data = vtk.vtkImageData()
# 设置vtkImageData的三维尺寸
image_data.SetDimensions(model_dims[2], model_dims[1], model_dims[0])
# 设置像素间隔
image_data.SetSpacing(1, 1, 1)
# 设置原点
image_data.SetOrigin(origin)
# 将数据导入vtkImageData数据对象
image_data.GetPointData().SetScalars(vtk_data)
# 定义vtkRenderer
renderer = vtk.vtkRenderer()
# 设置renderer相机焦点(x, y, z)
renderer.GetActiveCamera().SetFocalPoint((64, 90, 64))
# 设置renderer相机与焦点的距离(x, y, z + distance)
renderer.GetActiveCamera().SetPosition((64, 90, 630))
# 在水平面上旋转镜头10度
renderer.GetActiveCamera().Azimuth(-5)
# 创建一个vtkRenderWindow对象,用于显示vtkRenderer对象中的内容
render_window = vtk.vtkRenderWindow()
render_window.AddRenderer(renderer)
# 定义显示窗口尺寸
render_window.SetSize(640, 480)
# 创建一个vtkRenderWindowInteractor对象,用于处理交互事件
interactor = vtk.vtkRenderWindowInteractor()
interactor.SetRenderWindow(render_window)
# 定义xyz轴显示
axes = vtk.vtkAxesActor()
# 定义xyz轴长度
axes.SetTotalLength(250, 250, 250)
renderer.AddActor(axes)
# 定义一个平面与vtkImageReslice处于相同位姿,方便观察
plane_source = vtk.vtkPlaneSource()
# 设置原点
plane_source.SetOrigin(origin)
# 设置另外两个点, 三个点确定一个平面。同时定义方向与尺寸
plane_source.SetPoint1(168, 0, 0)
plane_source.SetPoint2(0, 168, 0)
# 使用vtkImageReslice获取切片
reslice = vtk.vtkImageReslice()
# 设置vtkImageReslice的输入
reslice.SetInputData(image_data)
# 设置插值方法
reslice.SetInterpolationModeToCubic()
# 设置间隔
reslice.SetOutputSpacing(1, 1, 1)
# 设置切片范围
reslice.SetOutputExtent(0, 128, 0, 128, 0, 0)
# 设置切片原点
reslice.SetResliceAxesOrigin(origin)
# 设置切片旋转角度
transform = vtk.vtkTransform()
transform.RotateX(x_angle)
transform.RotateY(y_angle)
transform.RotateZ(z_angle)
matrix = transform.GetMatrix()
reslice.SetResliceAxes(matrix)
reslice_axes = reslice.GetResliceAxes()
# 将平面与切片旋转角度保持一致
transform = vtk.vtkTransform()
transform.SetMatrix(reslice_axes)
transform_filter = vtk.vtkTransformPolyDataFilter()
transform_filter.SetInputConnection(plane_source.GetOutputPort())
transform_filter.SetTransform(transform)
mapper2 = vtk.vtkPolyDataMapper()
mapper2.SetInputConnection(transform_filter.GetOutputPort())
actor = vtk.vtkActor()
actor.SetMapper(mapper2)
actor.GetProperty().SetColor(1, 0, 0)
renderer.AddActor(actor)
# 设置切片输出维度
reslice.SetOutputDimensionality(2)
# 必须执行Update(),否则无法获取切片
reslice.Update()
# 应用体绘制算法,生成三维模型
volume_mapper = vtk.vtkGPUVolumeRayCastMapper()
volume_mapper.SetInputData(image_data)
volume_property = vtk.vtkVolumeProperty()
opacity_transfer_function = vtk.vtkPiecewiseFunction()
opacity_transfer_function.AddPoint(15, 0.0)
opacity_transfer_function.AddPoint(255, 1.0)
volume_property.SetScalarOpacity(opacity_transfer_function)
color_transfer_function = vtk.vtkColorTransferFunction()
color_transfer_function.AddRGBPoint(0, 0.0, 0.0, 0.0)
color_transfer_function.AddRGBPoint(255, 1.0, 1.0, 1.0)
volume_property.SetColor(color_transfer_function)
volume_property.SetScalarOpacityUnitDistance(1)
volume = vtk.vtkVolume()
volume.SetMapper(volume_mapper)
volume.SetProperty(volume_property)
renderer.AddVolume(volume)
renderer.SetBackground(1.0, 1.0, 0.8)
# 将VTK图像数据转换为NumPy数组
vtk_image = reslice.GetOutput()
height, width, _ = vtk_image.GetDimensions()
sc = vtk_image.GetPointData().GetScalars()
numpy_array = np.array(vtk_to_numpy(sc))
img = numpy_array.reshape(height, width)
img = rotate(img, 180, axes=(0, 1), reshape=True, order=1, mode='constant', cval=0, prefilter=True)
img = cv2.resize(img, (480, 480), interpolation=cv2.INTER_CUBIC)
cv2.imshow('rotate', img)
interactor.Initialize()
render_window.Render()
interactor.Start()
cv2.waitKey(0)
其他需求
- Q: vtkImageReslice的切片图像能否显示在plane(vtkPlaneSource)上?
A: 将vtkImageReslice转为纹理vtkTexture,然后设置为plane对应VtkActor的纹理。
# 假设 reslice 是你的 vtkImageReslice 对象
reslice.Update()
# 创建纹理
texture = vtk.vtkTexture()
texture.SetInputData(reslice.GetOutput())
# 创建平面源
plane = vtk.vtkPlaneSource()
# 创建映射器并设置纹理坐标
mapper = vtk.vtkPolyDataMapper()
mapper.SetInputConnection(plane.GetOutputPort())
mapper.SetScalarModeToUsePointData()
# 创建演员并设置纹理
actor = vtk.vtkActor()
actor.SetMapper(mapper)
actor.SetTexture(texture)