基于VTK的CPR曲面重建

一、问题

        在医学成像中，管状结构（如：血管、支气管和结肠）的评估是一个备受关注的课题。CT和MRI提供了人体的三维容积数据集，其中包含这些感兴趣的对象。然而，从CT和MRI获得的数据包含许多不感兴趣的部份，使得无预处理的体绘制（VR、MIP、MinIP、SSD）通常不可能或不准确。此外，感兴趣的对象很难完全位于一个平面内。

        一种可视化小直径结构的常用方法是根据从中心线检测过程中获得的高级信息重新采样并可视化数据集，这个过程被称为CPR(Curved Planar Reformation)曲面重建。通过这种技术，管状结构的整个长度显示在一幅图像中，便于医生对血管异常（如：狭窄、闭塞、动脉瘤和血管壁钙化）进行观察。

        CPR本质上是一种曲面的多平面重建技术MPR（Multi-Planar Reformation / Reconstruction)，是常用的后处理技术。MPR将在某个平面（通常是轴位）获取的成像模式中的数据转换到另一个平面。

        下面的 Python 代码在CT序列数据上，输入自定义的折线点坐标，采用插值算法拟合成一条曲线，沿曲线生成曲面。

二、代码

# CPR曲面重建import vtk;# 由曲线生成曲面
def SweepLine(line, direction, distance, cols):rows = line.GetNumberOfPoints();    #待生成的surface的行数为曲线上点的个数spacing = distance / cols;          #列方向上的spacingsurface = vtk.vtkPolyData();# 生成点cols += 1;numberOfPoints = rows * cols;             #surface内点的个数numberOfPolys = (rows - 1) * (cols - 1);  #单元个数points = vtk.vtkPoints();points.Allocate(numberOfPoints);polys = vtk.vtkCellArray();polys.Allocate(numberOfPolys * 4);#生成每行每列的点坐标，x = [0.0, 0.0, 0.0];cnt = 0;for row in range(rows):for col in range(cols):p = [0.0, 0.0, 0.0];line.GetPoint(row, p);x[0] = p[0] + direction[0] * col * spacing;x[1] = p[1] + direction[1] * col * spacing;x[2] = p[2] + direction[2] * col * spacing;points.InsertPoint(cnt, x);cnt += 1;# 生成四边形pts = [0, 0, 0, 0];for row in range(rows-1):for col in range(cols-1):pts[0] = col + row * cols;pts[1] = pts[0] + 1;pts[2] = pts[0] + cols + 1;pts[3] = pts[0] + cols;polys.InsertNextCell(4, pts);   #每临近的四个点构成一个单元surface.SetPoints(points);surface.SetPolys(polys);return surface;if __name__ == '__main__':inputImage = 'E://data//01//CT1.25mm';# 读取CT序列imageReader = vtk.vtkDICOMImageReader();imageReader.SetDataByteOrderToLittleEndian();imageReader.SetDirectoryName(inputImage);imageReader.Update();resolution = 100;extent = [255, 255];thickness = 1.0;spacing = [1.0, 1.0];# 输入自定义折线顶点坐标points = vtk.vtkPoints()# (a, b, c, d)其中a表示点的序号，(b,c,d)表示点的三维坐标points.InsertPoint(0, 1.0, 0.0, 0.0)points.InsertPoint(1, 200, 100, 0.0)points.InsertPoint(2, 100, 200, 0.0)points.InsertPoint(3, 200, 300, 0.0)# 将点插值拟合成一条曲线spline = vtk.vtkParametricSpline()spline.SetPoints(points)splineFilter = vtk.vtkParametricFunctionSource()splineFilter.SetParametricFunction(spline)splineFilter.Update()# 由曲线生成曲面direction = [0.0, 0.0, 1.0];distance = 160;surface = SweepLine(splineFilter.GetOutput(), direction, distance, resolution);sampleVolume = vtk.vtkProbeFilter();sampleVolume.SetInputConnection(1, imageReader.GetOutputPort());sampleVolume.SetInputData(0, surface);# 根据数据源的数据范围设置映射表的窗宽窗位wlLut = vtk.vtkWindowLevelLookupTable();range = imageReader.GetOutput().GetScalarRange()[1] - imageReader.GetOutput().GetScalarRange()[0];level = (imageReader.GetOutput().GetScalarRange()[1] + imageReader.GetOutput().GetScalarRange()[0]) / 2.0;wlLut.SetWindow(range);wlLut.SetLevel(level);# 创建映射器和角色mapper = vtk.vtkDataSetMapper();mapper.SetInputConnection(sampleVolume.GetOutputPort());mapper.SetLookupTable(wlLut);mapper.SetScalarRange(0, 255);actor = vtk.vtkActor();actor.SetMapper(mapper);# 创建渲染器、渲染窗口和交互ren = vtk.vtkRenderer();renwin = vtk.vtkRenderWindow();renwin.AddRenderer(ren);renwin.SetWindowName("CurvedReformation");iren = vtk.vtkRenderWindowInteractor();iren.SetRenderWindow(renwin);# 将角色添加到场景ren.AddActor(actor);colors = vtk.vtkNamedColors();ren.SetBackground(colors.GetColor3d("DarkSlateGray"));# 设置相机以查看图像ren.GetActiveCamera().SetViewUp(0, 0, 1);ren.GetActiveCamera().SetPosition(0, 0, 0);ren.GetActiveCamera().SetFocalPoint(1, 0, 0);ren.ResetCamera();# 渲染和交互renwin.Render();iren.Start();

三、结果