基于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;          #列方向上的spacing
    surface = 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();

三、结果