加入实验室后, 经过张老师的介绍, 有幸与某公司合共共同完成某个项目, 在此项目中我主要负责的是三维 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~
1. 基于面绘制的 MC 算法
(1)面绘制: 面绘制是采用分割技术对一系列的二维图像进行轮廓识别, 提取等操作, 最终还原出被检测物体的三维模型, 并以表面的方式显示出来.
(2)面绘制实现三维重建. 使用的是经典的 Marching Cubes 算法, 也叫移动立方体法.
(3)采用面绘制, VTK 中的数据流如下: source->filter(MC 算法或者 vtkContourFilter)->mapper->actor->render->renderwindow->interactor.
(4)MC 算法简介:
首先, 假定原始数据是离散的三维空间规则数据场,(断层扫描仪 CT 及核磁共振仪 MRI 产生的图像均属于这一类型);
其次, 给出所求等值面的值(为了在这一数据场中构造等值面);
最后, 找出等值面经过的体元位置, 求出该体元内的等值面并计算出相关参数(以便由常用的图形软件包或图形硬件提供的面绘制功能绘制等值
(5)VTK 提供了两种提取等值面的类: vtkContourFilter 滤波器和封装了 MC(Marching Cubes)算法类 vtkMarchingCubes. 提取等值面之后的数据处理: 通过 vtkPolyDataNormals 在等值面上产生法向量; 通过 vtkStripper 在等值面上产生纹理或三角面片.
(6)利用 MC 算法提取等值面的代码实现:
- import vtk
- # source->filter(MC 算法)->mapper->actor->render->renderwindow->interactor
- # 读取 Dicom 数据, 对应 source
- v16 = vtk.vtkDICOMImageReader()
- # v16.SetDirectoryName('D:/dicom_image/V')
- v16.SetDirectoryName('D:/dicom_image/vtkDicomRender-master/sample')
- # 利用封装好的 MC 算法抽取等值面, 对应 filter
- marchingCubes = vtk.vtkMarchingCubes()
- marchingCubes.SetInputConnection(v16.GetOutputPort())
- marchingCubes.SetValue(0, 100)
- # 剔除旧的或废除的数据单元, 提高绘制速度, 对应 filter
- Stripper = vtk.vtkStripper()
- Stripper.SetInputConnection(marchingCubes.GetOutputPort())
- # 建立映射, 对应 mapper
- mapper = vtk.vtkPolyDataMapper()
- # mapper.SetInputConnection(marchingCubes.GetOutputPort())
- mapper.SetInputConnection(Stripper.GetOutputPort())
- # 建立角色以及属性的设置, 对应 actor
- actor = vtk.vtkActor()
- actor.SetMapper(mapper)
- # 角色的颜色设置
- actor.GetProperty().SetDiffuseColor(1, .94, .25)
- # 设置高光照明系数
- actor.GetProperty().SetSpecular(.1)
- # 设置高光能量
- actor.GetProperty().SetSpecularPower(100)
- # 定义舞台, 也就是渲染器, 对应 render
- renderer = vtk.vtkRenderer()
- # 定义舞台上的相机, 对应 render
- aCamera = vtk.vtkCamera()
- aCamera.SetViewUp(0, 0, -1)
- aCamera.SetPosition(0, 1, 0)
- aCamera.SetFocalPoint(0, 0, 0)
- aCamera.ComputeViewPlaneNormal()
- # 定义整个剧院(应用窗口), 对应 renderwindow
- rewin = vtk.vtkRenderWindow()
- # 定义与 actor 之间的交互, 对应 interactor
- interactor = vtk.vtkRenderWindowInteractor()
- # 将相机添加到舞台 renderer
- renderer.SetActiveCamera(aCamera)
- aCamera.Dolly(1.5)
- # 设置交互方式
- style = vtk.vtkInteractorStyleTrackballCamera()
- interactor.SetInteractorStyle(style)
- # 将舞台添加到剧院中
- rewin.AddRenderer(renderer)
- interactor.SetRenderWindow(rewin)
- # 将角色添加到舞台中
- renderer.AddActor(actor)
- # 将相机的焦点移动至中央, The camera will reposition itself to view the center point of the actors,
- # and move along its initial view plane normal
- renderer.ResetCamera()
- interactor.Initialize()
- interactor.Start()
结果如下:
(7)利用 vtkContourFilter 滤波器提取等值面的代码实现:
- # 抽取轮廓 (等值面) 的操作对象是标量数据.
- # 其思想是: 将数据集中标量值等于某一指定恒量值的部分提取出来. 对于 3D 的数据集而言, 产生的是一个等值面; 对于 2D 的数据集而言, 产生的是一个等值线.
- # 其典型的应用有气象图中的等温线, 地形图中的等高线. 对于医学数据而言, 不同的标量值代表的是人体的不同部分, 因而可以分别提取出人的皮肤或骨头.
- # 抽取轮廓的功能是由一个过滤器实现的, 如 vtkContourFilter,vtkMarchingCubes.vtkContourFilter 可以接受任意数据集类型作为输入, 因而具有 一般性.
- # 使用 vtkContourFilter 时, 除了需要设置输入数据集外, 还需要指定一个或多个用于抽取的标量值. 可用如下两种方法进行设置.
- #
- # 使用方法 SetValue()逐个设置抽取值. 该方法有个两个参数: 第一个参数是抽取值的索引号, 表示第几个 抽取值. 索引号从 0 开始计数; 第二个参数就是指定的抽取值.
- # 使用方法 GenerateValues()自动产生一系列抽取值. 该方法有三个参数: 第一个参数是抽取值的个数, 后面两个参数是抽取值的取值范围.
- # coding=utf-8
- import vtk
- # source-filter--mapper--actor--render--renderwindow--interactor
- aRenderer = vtk.vtkRenderer() # 渲染器
- renWin = vtk.vtkRenderWindow() # 渲染窗口, 创建窗口
- renWin.AddRenderer(aRenderer) # 渲染窗口
- # renWin.Render()
- iren = vtk.vtkRenderWindowInteractor() # 窗口交互
- iren.SetRenderWindow(renWin)
- # The following reader is used to read a series of 2D slices(images)
- # that compose the volume.Theslicedimensions are set, and the
- # pixel spacing.The data Endianness must also be specified.The reader
- # uses the FilePrefix in combination with the slice number to construct
- # filenames using the format FilePrefix. % d.(In this case the FilePrefix
- # is the root name of the file.
- v16 = vtk.vtkDICOMImageReader()
- # v16.SetDirectoryName('D:/dicom_image/V')
- v16.SetDirectoryName('D:/dicom_image/vtkDicomRender-master/sample')
- # An isosurface, or contour value of 500 is known to correspond to the
- # skin of the patient.Once generated, a vtkPolyDataNormals filter is
- # used to create normals for smooth surface shading during rendering.
- skinExtractor = vtk.vtkContourFilter()
- skinExtractor.SetInputConnection(v16.GetOutputPort())
- skinExtractor.SetValue(0, -10)
- # skinExtractor.GenerateValues(2, 100, 110)
- skinNormals = vtk.vtkPolyDataNormals()
- skinNormals.SetInputConnection(skinExtractor.GetOutputPort())
- skinNormals.SetFeatureAngle(60.0)
- skinMapper = vtk.vtkPolyDataMapper() # 映射器
- skinMapper.SetInputConnection(skinNormals.GetOutputPort())
- skinMapper.ScalarVisibilityOff()
- skin = vtk.vtkActor()
- # 设置颜色 RGB 颜色系统就是由三个颜色分量: 红色 (R), 绿色(G) 和蓝色 (B) 的组合表示,
- # 在 VTK 里这三个分量的取值都是从 0 到 1,(0, 0, 0)表示黑色,(1, 1, 1)表示白色.
- # vtkProperty::SetColor(r,g, b)采用的就是 RGB 颜色系统设置颜色属性值.
- #skin.GetProperty().SetColor(0, 0, 1)
- skin.SetMapper(skinMapper)
- skin.GetProperty().SetDiffuseColor(1, .49, .25)
- skin.GetProperty().SetSpecular(.5)
- skin.GetProperty().SetSpecularPower(20)
- # skin.GetProperty().SetRepresentationToSurface()
- # 构建图形的方框
- outlineData = vtk.vtkOutlineFilter()
- outlineData.SetInputConnection(v16.GetOutputPort())
- mapOutline = vtk.vtkPolyDataMapper()
- mapOutline.SetInputConnection(outlineData.GetOutputPort())
- outline = vtk.vtkActor()
- outline.SetMapper(mapOutline)
- outline.GetProperty().SetColor(0, 0, 0)
- # 构建舞台的相机
- aCamera = vtk.vtkCamera()
- aCamera.SetViewUp(0, 0, -1)
- aCamera.SetPosition(0, 1, 0)
- aCamera.SetFocalPoint(0, 0, 0)
- aCamera.ComputeViewPlaneNormal()
- # Actors are added to the renderer.An initial camera view is created.
- # The Dolly() method moves the camera towards the Focal Point,
- # thereby enlarging the image.
- aRenderer.AddActor(outline)
- aRenderer.AddActor(skin)
- aRenderer.SetActiveCamera(aCamera)
- # 将相机的焦点移动至中央, The camera will reposition itself to view the center point of the actors,
- # and move along its initial view plane normal
- aRenderer.ResetCamera()
- # aCamera.Dolly(1.5)
- # aCamera.Roll(180)
- # aCamera.Yaw(60)
- aRenderer.SetBackground(250, 250, 250)
- # renWin.SetSize(640, 480)
- # 该方法是从 vtkRenderWindow 的父类 vtkWindow 继承过来的, 用于设置窗口的大小, 以像素为单位.
- renWin.SetSize(500, 500)
- aRenderer.ResetCameraClippingRange()
- style = vtk.vtkInteractorStyleTrackballCamera()
- iren.SetInteractorStyle(style)
- iren.Initialize()
- iren.Start()
结果如下:
(8)与可视化窗口的交互方式: 可以使用鼠标与三维图形交互, 比如用鼠标滚轮可以对三维图形放大, 缩小; 按下鼠标左键不放, 然后移动鼠标, 可以转动三维图形; 按下鼠标左键, 同时按下 Shift 键, 移动鼠标, 可以移动整个三维图形, 等等. 其他的功能你也可以试着摸索一下, 比如按下 Ctrl 键时再按鼠标左键; 鼠标停留在柱体上, 然后按下 P 键; 按一下字母 E 将关闭窗口.
(9)整个过程的简要总结:
(a)读取数据以及数据处理: 首先, 读取切片数据, 并将其转换为我们的开发工具 VTK 所支持的一种数据表达形式(vtkImageData). 我们给 CT 数据建立的是比较抽象的等值面模型, 最后将物理组件与抽象的模型结合在一起来建立对 CT 数据的可视化, 以帮助用户正确理解数据. 利用 VTK 中的 vtkDICOMImageReader 我们可以很方便的读取切片数据.
(b)提取等值面: 接着我们就可以用算法对所读取的数据进行处理了. 比如采用的经典 MC 的面绘制方法, 首先利用 vtkMarchingCubes 类来提取出某一 CT 值的等值面 (利用 vtksetValue() 来设置需要提取的值), 但这时的等值面其实仍只是一些三角面片, 还必须由 vtkStripper 类将其拼接起来形成连续的等值面. 这样就把读取的原始数据经过处理转换为应用数据, 也即由原始的点阵数据转换为多边形数据然后由 vtkPolyDataMapper 将其映射为几何数据, 并将其属性赋给窗口中代表它的演员, 将结果显示出来. 在实际应用中 Visualization Toolkit 支持多表面重建. 我们可以设置多个参数值, 提取出多个等值面并同时显示出来, 如何设置多个参数值呢? 可以通过 VTK 自带的 GenerateValues()函数. 常见的比如人体皮肤所对应的 value 值为 500, 人体骨骼所对应的 value 值为 1150.
(c)显示结果: 通过前面这些工作, 我们基本上已经完成了对数据的读取处理映射等步骤, 下面我们就要对数据进行显示了. 通常这些步骤也叫做渲染引擎. 可以通过调整 value 值和 actor 的相应属性达到重建三维图形的不同效果.
最主要就是设置相机, 设置 actor 的相关属性(颜色, 亮度, 透明度等等).
2. 基于体绘制的 Ray-casting 算法
(1)体绘制: 体绘制是将三维空间的离散数据直接转换为最后的立体, 图像而不必生成中间几何图元(面绘制需要), 其中心思想是为每一个体素指定一个不透明度, 并考虑每一个体素对光线的透射, 发射和反射作用.
(2)体绘制达到的效果: 体绘制的目标是在一副图片上展示空间体细节. 举例而言, 你面前有一间房子, 房子中有家具, 家电, 站在房子外面只能看到外部形状(类似于面绘制的效果), 无法观察到房子的布局或者房子中的物体; 假设房子和房子中的物体都是半透明的, 这样你就可以同时查看到所有的细节. 这就是体绘制所要达到的效果.
(3)体绘制常用的算法: 光线投射算法 ( Ray-casting ), 错切 - 变形算法( Shear-warp ), 频域体绘制算法( Frequency Domain ) 和抛雪球算法( Splatting ). 其中又以光线投射算法最为重要和通用.
(4)光线投射算法 ( Ray-casting ) 原理: 从图像平面的每个像素都沿着视线方向发出一条射线, 此射线穿过体数据集, 按一定步长进行采样, 由内插计算每个采样点的颜色值和不透明度, 然后由前向后或由后向前逐点计算累计的颜色值和不透明度值, 直至光线完全被吸收或穿过物体. 该方法能很好地反映物质边界的变化, 使用 Phong 模型, 引入镜面反射, 漫反射和环境反射能得到很好的光照效果, 在医学上可将各组织器官的性质属性, 形状特征及相互之间的层次关系表现出来, 从而丰富了图像的信息.(借鉴百度百科)
(5)体绘制的原理和面绘制完全不相同. 面绘制需要生成中间图元, 而体绘制则是直接在原图上进行绘制, 内容需求较面绘制小. 每切换一个视角需要重新对所有的像素点进行颜色和透明度计算, 需要时间比面绘制长.
(6)代码实现基于体绘制的 Ray-casting 算法
- # This example reads a volume dataset and displays it via volume rendering(体绘制).
- import vtk
- from vtk.util.misc import vtkGetDataRoot
- # Create the renderer, the render Windows, and the interactor. The renderer
- # draws into the render Windows, the interactor enables mouse- and
- # keyboard-based interaction with the scene.
- ren = vtk.vtkRenderer()
- renWin = vtk.vtkRenderWindow()
- renWin.AddRenderer(ren)
- iren = vtk.vtkRenderWindowInteractor()
- iren.SetRenderWindow(renWin)
- # The following reader is used to read a series of 2D slices (images)
- # that compose the volume. The slice dimensions are set, and the
- # pixel spacing. The data Endianness must also be specified. The reader
- # usese the FilePrefix in combination with the slice number to construct
- # filenames using the format FilePrefix.%d. (In this case the FilePrefix
- # is the root name of the file: quarter.)
- # v16 = vtk.vtkVolume16Reader()
- # v16.SetDataDimensions(64, 64)
- # v16.SetImageRange(1, 93)
- # v16.SetDataByteOrderToLittleEndian()
- # v16.SetFilePrefix("D:/dicom_image/headsq/quarter")
- # v16.SetDataSpacing(3.2, 3.2, 1.5)
- v16 = vtk.vtkDICOMImageReader()
- # v16.SetDirectoryName('D:/dicom_image/vtkDicomRender-master/sample')
- v16.SetDirectoryName('D:/dicom_image/V')
- # The volume will be displayed by ray-cast alpha compositing.
- # A ray-cast mapper is needed to do the ray-casting, and a
- # compositing function is needed to do the compositing along the ray.
- volumeMapper = vtk.vtkGPUVolumeRayCastMapper()
- volumeMapper.SetInputConnection(v16.GetOutputPort())
- volumeMapper.SetBlendModeToComposite()
- # The color transfer function maps voxel intensities to colors.
- # It is modality-specific, and often anatomy-specific as well.
- # The goal is to one color for flesh (between 500 and 1000)
- # and another color for bone (1150 and over).
- volumeColor = vtk.vtkColorTransferFunction()
- volumeColor.AddRGBPoint(0, 0.0, 0.0, 0.0)
- volumeColor.AddRGBPoint(500, 1.0, 0.5, 0.3)
- volumeColor.AddRGBPoint(1000, 1.0, 0.5, 0.3)
- volumeColor.AddRGBPoint(1150, 1.0, 1.0, 0.9)
- # The opacity transfer function is used to control the opacity
- # of different tissue types.
- volumeScalarOpacity = vtk.vtkPiecewiseFunction()
- volumeScalarOpacity.AddPoint(0, 0.00)
- volumeScalarOpacity.AddPoint(500, 0.15)
- volumeScalarOpacity.AddPoint(1000, 0.15)
- volumeScalarOpacity.AddPoint(1150, 0.85)
- # The gradient opacity function is used to decrease the opacity
- # in the "flat" regions of the volume while maintaining the opacity
- # at the boundaries between tissue types. The gradient is measured
- # as the amount by which the intensity changes over unit distance.
- # For most medical data, the unit distance is 1mm.
- volumeGradientOpacity = vtk.vtkPiecewiseFunction()
- volumeGradientOpacity.AddPoint(0, 0.0)
- volumeGradientOpacity.AddPoint(90, 0.5)
- volumeGradientOpacity.AddPoint(100, 1.0)
- # The VolumeProperty attaches the color and opacity functions to the
- # volume, and sets other volume properties. The interpolation should
- # be set to linear to do a high-quality rendering. The ShadeOn option
- # turns on directional lighting, which will usually enhance the
- # appearance of the volume and make it look more "3D". However,
- # the quality of the shading depends on how accurately the gradient
- # of the volume can be calculated, and for noisy data the gradient
- # estimation will be very poor. The impact of the shading can be
- # decreased by increasing the Ambient coefficient while decreasing
- # the Diffuse and Specular coefficient. To increase the impact
- # of shading, decrease the Ambient and increase the Diffuse and Specular.
- volumeProperty = vtk.vtkVolumeProperty()
- volumeProperty.SetColor(volumeColor)
- volumeProperty.SetScalarOpacity(volumeScalarOpacity)
- # volumeProperty.SetGradientOpacity(volumeGradientOpacity)
- volumeProperty.SetInterpolationTypeToLinear()
- volumeProperty.ShadeOn()
- volumeProperty.SetAmbient(0.9)
- volumeProperty.SetDiffuse(0.9)
- volumeProperty.SetSpecular(0.9)
- # The vtkVolume is a vtkProp3D (like a vtkActor) and controls the position
- # and orientation of the volume in world coordinates.
- volume = vtk.vtkVolume()
- volume.SetMapper(volumeMapper)
- volume.SetProperty(volumeProperty)
- # Finally, add the volume to the renderer
- ren.AddViewProp(volume)
- # Set up an initial view of the volume. The focal point will be the
- # center of the volume, and the camera position will be 400mm to the
- # patient's left (which is our right).
- camera = ren.GetActiveCamera()
- c = volume.GetCenter()
- camera.SetFocalPoint(c[0], c[1], c[2])
- camera.SetPosition(c[0] + 400, c[1], c[2])
- camera.SetViewUp(0, 0, -1)
- # Increase the size of the render Windows
- renWin.SetSize(640, 480)
- # Interact with the data.
- iren.Initialize()
- renWin.Render()
- iren.Start()
结果如下:
(7)体绘制的整个过程包括 VTK 数据量都与面绘制类似. 同样可以通过调整 actor 的相应属性达到重建三维图形的不同效果.
3. 心得体会与总结
(1)其实无论是面绘制还是体绘制都需要一定的 VTK 知识, 所以先了解 VTK 的一些基础知识才能帮助你更好的掌握这些方法.
(2)有关 VTK 整个数据流的过程可以用一下的例子进行类比, 方便理解(虽然这个类比不是非常形象):
当我们去看舞台剧的时候, 我们坐在台下, 展现在我们面前的是一个舞台, 舞台上有各式的灯光, 各样的演员. 演员出场的时候肯定是会先化妆, 有些演员可能会打扮成高富帅, 有些演员可能会化妆成白富美. 观众有时还会与台上的演员有一定的互动.
整个剧院就好比 VTK 程序的渲染窗口(vtkRenderWindow); 舞台就相当于渲染场景(vtkRenderer); 而那些高富帅, 白富美就是我们程序中的 Actor(有些文献翻译成 "演员", 有些翻译成 "角色", 这里我们不作翻译); 台上的演员与台下观众的互动可以看成是程序的交互(vtkRenderWindowInteractor); 演员与观众的互动方式有很多种, 现场的观众可以直接上台跟演员们握手拥抱, 电视机前的可以发短信, 电脑, 移动终端用户等可以微博关注, 加粉等等, 这就好比我们程序里的交互器样式(vtkInteractorStyle);
舞台上的演员我们都能一一分辨出来, 不会把高富帅弄混淆, 是因为他们化的妆, 穿的服饰都不一样, 这就相当于我们程序里 vtkActor 的不同属性(vtkProperty); 台下观众的眼睛可以看作是 vtkCamera, 前排的观众因为离得近, 看台上演员会显得比较高大, 而后排的观众看到的会显得小点, 每个观众看到的东西在他的世界里都是唯一的, 所以渲染场景 Renderer 里的 vtkCamera 对象也是只有一个; 舞台上的灯光可以有多个, 所以渲染场景里的 vtkLight 也存在多个
如下图所示, 可以加深理解
参考文献: https://vtk.org/doc/nightly/html/annotated.html
以上就是本次学习的内容, 欢迎交流与讨论
来源: https://www.cnblogs.com/XDU-Lakers/p/10822840.html