基于面绘制的MC算法以及基于体绘制的 Ray-casting 实现Dicom图像的三维重建(python实现)...

加入实验室后，经过张老师的介绍，有幸与某公司合共共同完成某个项目，在此项目中我主要负责的是三维 pdf 报告生成、Dicom图像上亮度、对比度调整以及 Dicom图像三维重建。今天主要介绍一下完成Dicom图像三维重建的过程以及自己的心得体会。实现Dicom三维图像重建最主要用的VTK(Visualization Toolkit，也就是可视化工具包)，由于今天的主题不是有关VTK，所以有关VTK的学习(包括VTK介绍、使用、实列)，可以参考此链接:https://blog.csdn.net/wishchin/article/details/12996693,个人建议：先把此教程中的前3个章节看完之后，在看此教程，这样能够更好的理解程序。接下来就让我们进入正题。

VTK将在可视化过程中经常遇到的细节屏蔽起来，并封装了一些常用的可视化算法，如将面绘制中常用的MC（MarchingCubes）算法和体绘制中常用的光线投射（Ray-Casting）算法封装成类的形式提供给使用者。这样在进行医学体数据的可视化时就可以直接使用VTK中已提供的相关类

整个项目的代码以及挂在GitHub上：https://github.com/tgpcai/Dicom_3D_Reconstruction，觉得哈不错的可以给楼主点一个start~

 

---

0.三维可视化的两种方式

（1）简单点说，三维可视化的目的就是让人们从屏幕上看到三维图像中有什么东西。众所周知，二维图像的显示是容易的，但是三维图像却不一样。过去由于技术的限制，得到了三维图像的数据，只能把它以二维图像的切片的形式展示给人看，这样显然不够直观。随着计算机科学的发展，使用计算机图形学技术为三维物体建模并实时渲染，实现场景漫游变成显示三维物体的主流方法，而过去切片显示的方式则逐渐被边缘化。

（2）由计算机图形学的知识我们可以知道，想显示三维图像中的内容，可以对这个“内容”的表面建立一个三角形网格模型。一旦得到了这个三角网格，那么渲染它就能够在屏幕上看到想要的内容，同时可以调节视角进行全方位的观察。所以第一类三维可视化方法就是基于这种思想：首先建立网格模型，之后渲染网格。这种方式被称为面绘制。

（3）还有一种叫做体绘制的方式，是直接将三维图像的体素点通过一定的透明度叠加计算后直接对屏幕上的像素点着色。这种方式的特点是能更加清楚的表现体数据内部细节，但是这种算法一般对计算机的压力也会比较大。

1.基于面绘制的MC算法

（0）首先基于MC的一系列算法需要明确一个“体元（Cell）”的概念。体元是在三维图像中由相邻的八个体素点组成的正方体方格，MarchingCubes算法的Cube的语义也可以指这个体元。注意区别体元和体素，体元是8个体素构成的方格，而每个体素（除了边界上的之外）都为8个体元所共享。

（1）面绘制：面绘制是采用分割技术对一系列的二维图像进行轮廓识别、提取等操作，最终还原出被检测物体的三维模型，并以表面的方式显示出来。

（2）面绘制实现三维重建。使用的是经典的 Marching Cubes 算法，也叫移动立方体法。

（3）采用面绘制，VTK中的数据流如下：source->filter(MC算法或者vtkContourFilter)->mapper->actor->render->renderwindow->interactor。

（4）MC算法简介：

• 首先，假定原始数据是离散的三维空间规则数据场，(断层扫描仪CT及核磁共振仪MRI产生的图像均属于这一类型)，读取这些数据，可得出这些数据的三个维度。
• 其次，以体元为单位来寻找三维图像中内容部分与背景部分的边界，在体元抽取三角片来拟合这个边界。
• 再者，遍历所有的体元，找出其中的三角片最后集合起来组成图像中实点表面的三角网格（Mesh）。
• 最后，建立好了三角形网格模型，对该模型进行渲染。

（5）VTK提供了两种提取等值面的类：vtkContourFilter滤波器和封装了MC(Marching Cubes)算法类vtkMarchingCubes。提取等值面之后的数据处理：通过vtkPolyDataNormals在等值面上产生法向量；通过vtkStripper在等值面上产生纹理或三角面片。

（6）利用MC算法提取等值面的代码实现：

 1 import vtk
 2 # source->filter(MC算法)->mapper->actor->render->renderwindow->interactor
 3 
 4 # 读取Dicom数据，对应source
 5 v16 = vtk.vtkDICOMImageReader()
 6 # v16.SetDirectoryName('D:/dicom_image/V')
 7 v16.SetDirectoryName('D:/dicom_image/vtkDicomRender-master/sample')
 8 
 9 # 利用封装好的MC算法抽取等值面，对应filter
10 marchingCubes = vtk.vtkMarchingCubes()
11 marchingCubes.SetInputConnection(v16.GetOutputPort())
12 marchingCubes.SetValue(0, 100)
13 
14 # 剔除旧的或废除的数据单元，提高绘制速度，对应filter
15 Stripper = vtk.vtkStripper()
16 Stripper.SetInputConnection(marchingCubes.GetOutputPort())
17 
18 # 建立映射，对应mapper
19 mapper = vtk.vtkPolyDataMapper()
20 # mapper.SetInputConnection(marchingCubes.GetOutputPort())
21 mapper.SetInputConnection(Stripper.GetOutputPort())
22 
23 # 建立角色以及属性的设置，对应actor
24 actor = vtk.vtkActor()
25 actor.SetMapper(mapper)
26 # 角色的颜色设置
27 actor.GetProperty().SetDiffuseColor(1, .94, .25)
28 # 设置高光照明系数
29 actor.GetProperty().SetSpecular(.1)
30 # 设置高光能量
31 actor.GetProperty().SetSpecularPower(100)
32 
33 # 定义舞台，也就是渲染器，对应render
34 renderer = vtk.vtkRenderer()
35 
36 # 定义舞台上的相机，对应render
37 aCamera = vtk.vtkCamera()
38 aCamera.SetViewUp(0, 0, -1)
39 aCamera.SetPosition(0, 1, 0)
40 aCamera.SetFocalPoint(0, 0, 0)
41 aCamera.ComputeViewPlaneNormal()
42 
43 # 定义整个剧院(应用窗口)，对应renderwindow
44 rewin =