Python-医学图像处理之三维重建（进行切片级重建）

       对于从事医学图像处理的小伙伴而言，医学图像三维重建并不是一个陌生的东西啦~ 例如，在对图像分割结果进行展示或者验证时，我们常常通过对分割结果进行三维可视化的方式进行展示和说明。

       那废话不多说，今天就来教大家如何根据自己的分割结果进行三维重建。

       这里呢，我用现在正在做的韧带分割进行说明。

       首先，通过深度学习或者传统方法对医学图像进行分割，得到二值化的分割结果（如图1所示）。将分割结果放置在一个文件夹里。

图1 分割结果

       这里要注意，我们的分割结果的图片名称，需要按医学图像切片的顺序命名，不然重建时会出现错乱！！！！

       然后，读入数据。这里呢，我们使用的是vtk库。

reader = vtk.vtkPNGReader()
reader.SetDataScalarTypeToUnsignedChar()
reader.SetFileDimensionality(3)

path = r"./sag/" #读取 图像分割结果文件夹 路径

reader.SetFilePrefix(path)  # 文件夹路径
reader.SetFilePattern("%s%d.png")

       如果是JPG或JPEG格式的图片，可以使用vtk.vtkJPEGReader()函数进行替换。当然png对应的地方也不要忘记进行替换哦~

       接下来，是一些参数的设置。这些参数将决定我们重建结果的正确性！！！

reader.SetDataExtent(3, 384, 3, 384, 1, 50) #384是图像大小， 倒数两个参数（1,50）：重建的切片范围

       这里，需要修改图像大小和重建范围。这里对重建范围进行一下说明，我这里的1,50也就是说，从我分割结果的第1张到第50张进行重建，小伙伴们可以根据自己的需求进行修改。

reader.SetDataSpacing(0.5625, 0.5625, 1.76) #前两个是 像素距离，最后一个是 层厚

       这里需要修改的是像素距离和层厚，这个地方，小伙伴们根据自己的医学图像数据进行对应的修改就行。

      之后呢，就是vtk三维重建的常规套路啦，建立actor、mapper等等，这里就不过多赘述啦~

      这里给大家展示一下我重建的结果（如图2所示）。

图2 重建结果

       最后附上代码，希望对各位小伙伴有所帮助鸭~

import vtk

reader = vtk.vtkPNGReader()
reader.SetDataScalarTypeToUnsignedChar()
reader.SetFileDimensionality(3)

path = r"./sag/" #读取文件路径
path2 = "./FX_sag_patella_lig.ply" #保存模型路径

reader.SetFilePrefix(path)  # 文件夹路径
reader.SetFilePattern("%s%d.png")

reader.SetDataExtent(3, 384, 3, 384, 1, 50) #384是图像大小， 倒数两个参数（1,50）：重建的切片范围
reader.SetDataSpacing(0.5625, 0.5625, 1.76) #前两个是 像素距离，最后一个是 层厚

reader.SetDataOrigin(0, 0, 0)
reader.Update()

skinExtractor = vtk.vtkContourFilter()
skinExtractor.SetInputConnection(reader.GetOutputPort())
skinExtractor.SetValue(0, 255)
skinExtractor.ComputeGradientsOn()
skinExtractor.ComputeScalarsOn()

smooth = vtk.vtkSmoothPolyDataFilter()
smooth.SetInputConnection(skinExtractor.GetOutputPort())
smooth.SetNumberOfIterations(300)

skinNormals = vtk.vtkPolyDataNormals()
skinNormals.SetInputConnection(smooth.GetOutputPort())
skinNormals.SetFeatureAngle(100)
skinNormals.Update()

skinMapper = vtk.vtkPolyDataMapper()
skinMapper.SetInputConnection(skinNormals.GetOutputPort())
skinMapper.ScalarVisibilityOff()

skin = vtk.vtkActor()
skin.SetMapper(skinMapper)
skin.SetOrigin(skin.GetCenter())

arender = vtk.vtkRenderer()
renWin = vtk.vtkRenderWindow()
renWin.AddRenderer(arender)
iren = vtk.vtkRenderWindowInteractor()
iren.SetRenderWindow(renWin)

arender.AddActor(skin)

arender.ResetCamera()
arender.SetBackground(.2, .3, .4)
arender.ResetCameraClippingRange()

renWin.SetSize(500, 500)
style = vtk.vtkInteractorStyleTrackballCamera()
iren.SetInteractorStyle(style)

vtkWriter = vtk.vtkPLYWriter()
vtkWriter.SetInputData(skinNormals.GetOutput())
vtkWriter.SetFileName(path2)
vtkWriter.Write()

renWin.Render()
iren.Initialize()
iren.Start()