VTK图像处理之vtkImageReslice(Python 版本)

三维图像切面提取 切片（Slice）或切面是三维图像比较常用的概念，尤其在医学图像中。通过提取切面可以方便地浏览和分析图像内部组织结构。VTK中vtkImageReSlice类可以实现图像切面的提取。在实际开发中，四视图中冠状视面、矢状面和横断面（显示过图像内部一点且平行于XY、YZ、XZ平面的平面），需要用到此类。

在显示图像时，很多时候我们需要的不是矢状面，冠状面或者横断面这种与坐标轴平行的正交面，而是一些斜切面。这时候就可以利用vtkImageReslice 这个类来实现。

这个类的输入最主要的是SetResliceAxes()这个函数，这个函数的输入是vtkMatrix4x4 其代表的是新的坐标系在世界坐标系下的坐标值。可以把vtkMatrix4x4看成是4列，其中第一，二，三列分别为新的坐标系在世界坐标系下的向量坐标，第四列代表的是新的坐标原点在世界坐标系下的坐标值。如下面的四列值就表示新的坐标系原点为（100,200,10），x方向向量坐标为（1,0,0），y方向向量坐标为（0,0,-1），z方向向量坐标为（0,1,0）。

[1,0,0,100;0,0,1,200;0,-1,0,10;0,0,0,1]

话不多说先看结果！

import vtkmodules.all as vtk
import numpy as np
import cv2
file_path_t2 = "E:/ShoulderData/cai/reconstruction/cai_shoulder_T2_5.nii.gz"
file_path_t1 = "E:/ShoulderData/cai/reconstruction/cai_shoulder_T1_4.nii.gz"#数据路径

readerT2 = vtk.vtkNIFTIImageReader()
readerT2.SetFileName(file_path_t2)
readerT2.Update()

readerT1 = vtk.vtkNIFTIImageReader()
readerT1.SetFileName(file_path_t1)
readerT1.Update()

extentT2 = readerT2.GetOutput().GetExtent()
spacingT2 = readerT2.GetOutput().GetSpacing()
originT2 = readerT2.GetOutput().GetOrigin()

extentT1 = readerT1.GetOutput().GetExtent()
spacingT1 = readerT1.GetOutput().GetSpacing()
originT1 = readerT1.GetOutput().GetOrigin()

center = [0.0, 0.0, 0.0]
center[0] = originT2[0] + spacingT2[0] * 0.5 * (extentT2[0] + extentT2[1])#计算中心点的坐标
center[1] = originT2[1] + spacingT2[1] * 0.5 * (extentT2[2] + extentT2[3])
center[2] = originT2[2] + spacingT2[2] * 0.5 * (extentT2[4] + extentT2[5])

axialElements = [
    1, 0, 0, 0,
    0, 1, 0, 0,
    0, 0, 1, 0,
    0, 0, 0, 1
]

resliceAxesT2 = vtk.vtkMatrix4x4()
resliceAxesT2.DeepCopy(axialElements)
resliceAxesT2.SetElement(0, 3, center[0])
resliceAxesT2.SetElement(1, 3, center[1])
resliceAxesT2.SetElement(2, 3, center[2])

resliceAxesT1 = vtk.vtkMatrix4x4()
resliceAxesT1.DeepCopy(axialElements)
resliceAxesT1.SetElement(0, 3, center[0])
resliceAxesT1.SetElement(1, 3, center[1])
resliceAxesT1.SetElement(2, 3, center[2])

resliceT2 = vtk.vtkImageReslice()
resliceT2.SetInputConnection(readerT2.GetOutputPort())
resliceT2.SetOutputDimensionality(2)
resliceT2.SetResliceAxes(resliceAxesT2)
resliceT2.SetInterpolationModeToLinear()

resliceT1 = vtk.vtkImageReslice()
resliceT1.SetInputConnection(readerT1.GetOutputPort())
resliceT1.SetOutputDimensionality(2)
resliceT1.SetResliceAxes(resliceAxesT1)
resliceT1.SetInterpolationModeToLinear()

# Update reslice filters to compute the output images
resliceT2.Update()
resliceT1.Update()

# Get output images from the reslice filters
outputT2 = resliceT2.GetOutput()
outputT1 = resliceT1.GetOutput()

# Convert the output images to numpy arrays
arrayT2 = np.reshape(outputT2.GetPointData().GetScalars(),(outputT2.GetDimensions()[0],outputT2.GetDimensions()[1]))
arrayT1 = np.reshape(outputT1.GetPointData().GetScalars(),(outputT1.GetDimensions()[0],outputT1.GetDimensions()[1]))

# Normalize and convert to uint8
arrayT2 = (arrayT2 - arrayT2.min()) / (arrayT2.max() - arrayT2.min()) * 255 # 归一化
arrayT1 = (arrayT1 - arrayT1.min()) / (arrayT1.max() - arrayT1.min()) * 255

arrayT2 = arrayT2.astype('uint8')
arrayT1 = arrayT1.astype('uint8')



cv2.imshow("T1",arrayT1)
cv2.imshow("T2",arrayT2)
# cv2.imshow("rezie_t2",rezie_t2)

cv2.waitKey(0)