一个读取CT图像序列，并进行表面重建的C++代码

不拱地的猪

已于 2024-10-09 10:54:26 修改

阅读量953

点赞数 19

文章标签： c++ vtk CatPhan503

于 2024-10-08 16:45:46 首次发布

本文链接：https://blog.csdn.net/kangdehua/article/details/142760667

版权

这篇文章中，介绍使用VTK进行读取CT图像（一个序列），然后进行表面重建。为什么不使用DCMTK呢？因为使用DCMTK需要一张一张读取，要自己写一个代码，还要创建一个容器来放读入的CT数据，比较复杂。在实际的工程中，我们需要寻找合适的工具来造出漂亮的项目。

1 为什么使用VTK

VTK自己就带有DICOM文件读取类，同时VTK中的imagedata 类可以存储一个图像矩阵，是一个很受欢迎的容器。VTK中还有表面重建算法marchingcubes，也很好用的。

1.1 vtkImageData的特点

数据结构：vtkImageData是基于结构化网格的数据表示形式，数据点被组织成一个规则的二维或三维数组。这些数据点称为像素或格点，每个像素包含一个或多个数值表示的属性信息。vtkImageData使用这种结构方便地表示二维图像和三维体数据。
数据类型：vtkImageData支持多种数据类型，如VTK_FLOAT、VTK_INT、VTK_SHORT、VTK_UNSIGNED_CHAR等，可以满足不同图像数据的存储需求。
维度信息：vtkImageData对象的维度描述了其包含像素的行数、列数和深度（对于三维图像）。这些维度信息决定了数据在空间中的范围和形状。
原点与间距：vtkImageData具有原点和间距属性，它们定义了空间中像素的位置和大小。原点表示数据中第一个像素的位置坐标，间距定义了相邻像素在各个方向上的间隔距离。
灵活性：vtkImageData提供了丰富的方法和属性，用于访问和修改单个像素的数值，进行图像处理以及数据分析。同时，它还支持与其他VTK类进行集成，实现更复杂的可视化和分析任务。

1.2 vtkMarchingCubes类的特点

等值面提取：
vtkMarchingCubes类主要用于从三维规则体数据中提取等值面。它通过分析体数据中的每个体素（小立方体），并根据给定的等值面值，构造出逼近等值面的三角面片。
算法效率：
Marching Cubes算法实际上是一个分而治之的方法，它将等值面的抽取分布于每一个体素中进行。这种分布式处理使得算法在处理大规模数据时具有较高的效率。
适用性广泛：
vtkMarchingCubes类不仅适用于医学图像重建，如CT和MRI数据，还广泛应用于地质勘探、气象预报、流体力学等领域中的三维数据可视化。
法向量计算：
vtkMarchingCubes类在提取等值面的同时，还可以计算每个三角面片的法向量。法向量的计算对于提高渲染质量和实现光照效果至关重要。
集成性强：
vtkMarchingCubes类与VTK库中的其他类具有良好的集成性，可以方便地与其他可视化工具结合使用，实现更复杂的可视化任务。

2 CPP代码实践

    string filePath="../data";  //定义文件路径（CT序列存放的位置）
    // 定义一个reader指针，指向vtkDICOMIMageReader类
    vtkSmartPointer<vtkDICOMImageReader> reader=vtkSmartPointer<vtkDICOMImageReader>::New();
    reader->SetDirectoryName(filePath.c_str());  //设置路径
    reader->Update();                            //读取
    float* ipp;                                  //定义一个指向ImagePositionPatient类的指针，该方法不安全
    ipp=reader->GetImagePositionPatient();
    // 定义一个指向ImageData容器类的指针，用来存储读入的CT图像。
    vtkSmartPointer<vtkImageData> originalData = reader->GetOutput();
    int dim[3];    // 维度
    double spa[3], ori[3];   //空间分辨率，原点
    originalData->GetDimensions(dim);
    originalData->GetSpacing(spa);
    originalData->GetOrigin(ori);
    //输出
    cout<<"Origin: "<<ori[0]<<","<<ori[1]<<","<<ori[2]<<endl;
    cout<<"ImagePositionPatient: "<<ipp[0]<<","<<ipp[1]<<","<<ipp[2]<<endl;
    cout<<"Dimension: "<<dim[0]<<","<<dim[1]<<","<<dim[2]<<endl;
    cout<<"Spacing: "<<spa[0]<<","<<spa[1]<<","<<spa[2]<<endl;

输出结果为：

Origin: 0,0,0
ImagePositionPatient: -449.414,-449.414,-80
Dimension: 384,384,64
Spacing: 1.17188,1.17188,2.5

我们可以看到原点在0,0,0,这个是VTK图像自定义的原点，在图像体积的左下角，如下图所示：

图1 VTK中默认的坐标系位置（红色到绿色到蓝色坐标系分别为X,Y,Z，右手系）

在实际的应用中，我们知道CT扫描的时候是有一个物理坐标系的。它的中心点（0,0,0）大致在CT孔颈的中心处，同时处在那一层的CT图像就是0层面。因为在实际的临床运用中，是需要坐标系的，比如远距离放射治疗，立体定向手术等。

使用vtkMarchingCubes类来获取表面重建

// 定义MarchingCubes指针
 vtkSmartPointer<vtkMarchingCubes> mc=vtkSmartPointer<vtkMarchingCubes>::New();
    mc->SetInputData(originalData);   //输入
    mc->SetValue(0,-100);             //阈值
    mc->Update();                     //数据更新

vtkSmartPointer<vtkPolyDataMapper> mapper=vtkSmartPointer<vtkPolyDataMapper>::New();
    mapper->SetInputConnection(mc->GetOutputPort());   //相当于演员的衣服

    vtkSmartPointer<vtkActor> actor=vtkSmartPointer<vtkActor>::New();  //演员
    actor->SetMapper(mapper);    //演员穿衣服
vtkSmartPointer<vtkAxesActor> axes=vtkSmartPointer<vtkAxesActor>::New();  //坐标轴
    axes->SetPosition(0,0,0);      //位置
    axes->AxisLabelsOff();         // 不显示X，Y，Z（很丑）
    axes->SetTotalLength(150,150,150);   // 坐标轴长度
    // 舞台
    vtkSmartPointer<vtkRenderer> renderer = vtkSmartPointer<vtkRenderer>::New();

    renderer->AddActor(axes);  //加入坐标轴
    renderer->AddActor(actor);  //加入演员

    vtkSmartPointer<vtkRenderWindow> renWin = vtkSmartPointer<vtkRenderWindow>::New();   //渲染窗口（可以调整大小）
    vtkSmartPointer<vtkRenderWindowInteractor> iren=vtkSmartPointer<vtkRenderWindowInteractor>::New();   //与窗口内容交互
    vtkSmartPointer<vtkInteractorStyleTrackballCamera> style=vtkSmartPointer<vtkInteractorStyleTrackballCamera>::New();  //交互模式（我比较喜欢）
    renWin->AddRenderer(renderer);
    renWin->SetSize(600,600);
    //renWin->Render();
    iren->SetInteractorStyle(style);
    iren->SetRenderWindow(renWin);
    iren->Initialize();
    iren->Start();

输出的结果如下：