分别用 VTK 体绘制和面绘制来实现医学图像三维重建

简介: 序言,VTK介绍:VTK 全称为 The Visualization Toolkit (可视化工具),是一个开源、跨平台、自由获取、支持并行计算的图形应用函数;拥有3D 渲染的最新工具、提供3D交互模式以及2D绘图等。

序言,VTK介绍:


VTK 全称为 The Visualization Toolkit (可视化工具),是一个开源、跨平台、自由获取、支持并行计算的图形应用函数;拥有3D 渲染的最新工具、提供3D交互模式以及2D绘图等。

VTK 包含一个C++类库,目前提供了众多语言接口,例如 Java、Python、TCL;在三维函数库OpenGL 的基础上采用面向对象设计方法发展起来

图形学基本概念和数据结构,是VTK的核心,VTK是通过 Pipline的形式来输送数据,实现预览效果。


三维重建


在 VTK 中,提供了两种重建方式:体绘制和面绘制 (一般来说用VTK做重建,医学图像领域较多,如 Dicom、mha、mhd;当然 VTK 也实现点云重建)

面绘制

利用面绘值用到VTK封装到的 Marching Cube 算法,简称 MC算法,MC 算法的实现主要分为三部分:

1,确定包含等值面的体元

首先介绍一下 体元的概念,体元是三维图像中由相邻的八个体素点组成的正方体方格,英语也叫 Cube,体元中角点函数值分为两种情况,一种是大于等于给定等值面的值 C0 ,则将角点设为 1 称该角点在等值面内部,否则设为0,在等值面之外,


一般来说,会出现一个角点在内,一个角点在外,则角点之间的连线(也就是体元的边)必然与等值面相交,根据这个原理就能判断等值面与哪些体元相交。

体元内每个角点(顶点)有两种情况:0和1,一共8个角点即分为256种( 2 8 = 256 2^8 = 2562

8

=256 ),根据平面对称性、中心对称性,256种最终降到15种


2,确定等值面与体元边界的交点

找到含有等值面的体元之后,接下来就是确定等值面与体元边界的交点,体元间的数值都是呈线性变化,求交点时一般采用的是线性插值,如 Case0 中等值面的两个端点 一个在外为( 标记0) ,一个在内 ( 标记为1 ) 则交点为0.5;


3,求等值面的法向量

以上步骤 1,2,3 为实现 MC 算法步骤流程,但利用 VTK ,不需要这么繁琐,主要算法步骤都已经封装到 vtkMarchingCube 类中,使用 vtkMarchingCube 时,需要设置三个参数:


SetValue(int i,double value) 设置第i 个等值面的值b,(提醒一下,医学图像中的灰度值范围不是 0-256 而是0-65326,但大部分取值范围都在0-1000)。


SetNumberofContours(int number),设置等值面的个数


ComputerNormalsOn() 设置计算等值面的法向量,提高渲染质量;

上面这张图显示的就是 vtk 呈像的基本流程,下面是仿照官网写的用面绘制来对图像重建的代码部分:

#include<vtkRenderWindow.h>
#include<vtkRenderWindowInteractor.h>
#include<vtkDICOMImageReader.h>
#include<vtkMarchingCubes.h>
#include<vtkPolyDataMapper.h>
#include<vtkStripper.h>
#include<vtkActor.h>
#include<vtkProperty.h>
#include<vtkCamera.h>
#include<vtkOutlineFilter.h>
#include<vtkOBJExporter.h>
#include<vtkRenderer.h>
#include<vtkMetaImageReader.h>
#include<vtkInteractorStyleTrackballCamera.h>
#include<iostream>
#include<string.h>
//需要进行初始化,否则会报错
#include <vtkAutoInit.h> 
#include<vtkRenderingVolumeOpenGL2ObjectFactory.h>
#include<vtkRenderingOpenGL2ObjectFactory.h>
using namespace std;
int main()
{
  ///Marching Cube; 
  vtkObjectFactory::RegisterFactory(vtkRenderingOpenGL2ObjectFactory::New());
  vtkObjectFactory::RegisterFactory(vtkRenderingVolumeOpenGL2ObjectFactory::New());
  vtkSmartPointer<vtkRenderer> ren = vtkSmartPointer<vtkRenderer>::New();
  vtkSmartPointer<vtkRenderWindow> renWin = vtkSmartPointer<vtkRenderWindow>::New();//WINDOW;
  renWin->AddRenderer(ren);
  vtkSmartPointer<vtkRenderWindowInteractor> iren = vtkSmartPointer<vtkRenderWindowInteractor>::New();//wininteratcor;
  iren->SetRenderWindow(renWin);
  vtkSmartPointer<vtkDICOMImageReader> reader = vtkSmartPointer<vtkDICOMImageReader>::New();
  reader->SetDirectoryName("E:/DIcom_Data/DICOM");
  reader->SetDataByteOrderToLittleEndian();
  reader->Update();
  /*vtkDICOMImageReader *reader = vtkDICOMImageReader::New();
  reader->SetDirectoryName("E:/Coding Pra/VTK/VTK_Examples_StandardFormats_Input_DicomTestImages/DICOM");
  reader->SetDataByteOrderToLittleEndian();
  reader->Update();*/
  cout << "读取数据完毕" << endl;
  cout << "The width is" << reader->GetWidth() << endl;
  cout << "The height is" << reader->GetHeight() << endl;
  cout << "The depth is" << reader->GetPixelSpacing() << endl;
  cout << "The Output port is" << reader->GetOutputPort() << endl;
  vtkSmartPointer<vtkMarchingCubes> marchingcube = vtkSmartPointer<vtkMarchingCubes>::New();
  marchingcube->SetInputConnection(reader->GetOutputPort());//获得读取的数据的点集;
  marchingcube->SetValue(0, 200);//Setting the threshold;
  marchingcube->ComputeNormalsOn();//计算表面法向量;
  vtkSmartPointer<vtkStripper> Stripper = vtkSmartPointer<vtkStripper>::New();
  Stripper->SetInputConnection(marchingcube->GetOutputPort());//获取三角片
  vtkSmartPointer<vtkPolyDataMapper> Mapper = vtkSmartPointer<vtkPolyDataMapper>::New();//将三角片映射为几何数据;
  Mapper->SetInputConnection(Stripper->GetOutputPort());
  Mapper->ScalarVisibilityOff();//
  vtkSmartPointer<vtkActor> actor = vtkSmartPointer<vtkActor>::New();//Created a actor;
  actor->SetMapper(Mapper);//获得皮肤几何数据
  actor->GetProperty()->SetDiffuseColor(1, .49, .25);//设置皮肤颜色;
  actor->GetProperty()->SetSpecular(0.3);//反射率;
  actor->GetProperty()->SetOpacity(1.0);//透明度;
  actor->GetProperty()->SetSpecularPower(20);//反射光强度;
  actor->GetProperty()->SetColor(1, 0, 0);//设置角的颜色;
  actor->GetProperty()->SetRepresentationToWireframe();//线框;
  //vtkSmartPointer<vtkCamera> camera = vtkSmartPointer<vtkCamera>::New();//Setting the Camera;
  //camera->SetViewUp(0, 0, -1);//设置相机向上方向;
  //camera->SetPosition(0, 1, 0);//位置:世界坐标系,相机位置;
  //camera->SetFocalPoint(0, 0, 0);//焦点,世界坐标系,控制相机方向;
  //camera->ComputeViewPlaneNormal();//重置视平面方向,基于当前的位置和焦点;
  vtkSmartPointer<vtkOutlineFilter> outfilterline = vtkSmartPointer<vtkOutlineFilter>::New();
  outfilterline->SetInputConnection(reader->GetOutputPort());
  vtkSmartPointer<vtkPolyDataMapper> outmapper = vtkSmartPointer<vtkPolyDataMapper>::New();
  outmapper->SetInputConnection(outfilterline->GetOutputPort());
  vtkSmartPointer<vtkActor> OutlineActor = vtkSmartPointer<vtkActor>::New();
  OutlineActor->SetMapper(outmapper);
  OutlineActor->GetProperty()->SetColor(0, 0, 0);//线框颜色
  ren->AddActor(actor);
  ren->AddActor(OutlineActor);
  //ren->SetActiveCamera(camera);//设置渲染器的相机;
  ren->ResetCamera();
  ren->ResetCameraClippingRange();
  //camera->Dolly(1.5);//使用Dolly()方法延伸着视平面法向移动相机;
  ren->SetBackground(1, 1, 1);//设置背景颜色;
  renWin->SetSize(1000, 600);
  vtkInteractorStyleTrackballCamera *style = vtkInteractorStyleTrackballCamera::New();
  iren->SetInteractorStyle(style);
  renWin->Render();
  iren->Initialize();
  iren->Start();
  vtkSmartPointer<vtkOBJExporter> porter = vtkSmartPointer<vtkOBJExporter>::New();
  porter->SetFilePrefix("E:/ceshi/aaa/regist_after/polywrite.obj");//重建图像输出
  porter->SetInput(renWin);
  porter->Write();
  return EXIT_SUCCESS;
}


上面就是 VTK 基于 Marching Cube算法实现的重建效果:

体绘制重建

体绘制时分为两部分

1,定义 vtkVoluemRayCastMapper 对象

体绘制中最常用的方法 ;vtkVolumeRayCastMapper() 光线投影,体绘制时,首先定义一个Mapper 然后接受两个输入:


SetInput(vtkImageDate *) 用于设置输入图像数据;

SetVolumeRayCastFunction(vtkVolumeRayCastFunction *) 用于设置光线投影函数类型;


2,利用 vtkVolumeProperty 定义体绘制属性;

SetScalarOpacity() 设置灰度不透明函数;

SetColor() 颜色传输函数;


3, 定义 vtkVolume 对象接收 Mapper对象和 Property 对象

SetMapper()接受 Mapper 对象;

SetProperty() 接受 Property 对象;

vtk 中体绘制 核心就是改变 Mapper 和 vtkVolumeRayCastFunction() ,上面中vtkColumeRayCastMapper 只是 VolumeMapper 其中的一种,且投影函数类 vtkVolumeRayCastFunction 一共有三个子类:


vtkVolumeRayCastCompositeFunction

vtkVolumeRayCasMIPFunction、

vtkVolumeRayCastIsosurfaceFunction,

因此,其细分的话vtk中的体绘制也不止一种

而下面这个是最常用到的(``vtkVolumeRayCastMapper+vtkVolumeRayCastCompositeFunction`)

//体绘制
#include<vtkRenderWindowInteractor.h>
#include<vtkDICOMImageReader.h>
#include<vtkCamera.h>
#include<vtkActor.h>
#include<vtkRenderer.h>
#include<vtkVolumeProperty.h>
#include<vtkProperty.h>
#include<vtkPolyDataNormals.h>
#include<vtkImageShiftScale.h>
#include "vtkVolumeRayCastMapper.h"
#include<vtkPiecewiseFunction.h>
#include<vtkColorTransferFunction.h>
#include<vtkVolumeRayCastCompositeFunction.h>
#include<vtkRenderWindow.h>
#include<vtkImageCast.h>
#include<vtkVolumeRayCastCompositeFunction.h>
#include<vtkOBJExporter.h>
#include<vtkOutlineFilter.h>
#include<vtkPolyDataMapper.h>
#include<vtkInteractorStyleTrackballCamera.h>
#include<vtkRenderingVolumeOpenGL2ObjectFactory.h>
#include<vtkRenderingOpenGL2ObjectFactory.h>
#include<vtkMetaImageReader.h>
#include<vtkLODProp3D.h>
//体绘制加速
//Gpu光照映射
#include<vtkGPUVolumeRayCastMapper.h>
#include<iostream>
int main()
{
  vtkObjectFactory::RegisterFactory(vtkRenderingOpenGL2ObjectFactory::New());
  vtkObjectFactory::RegisterFactory(vtkRenderingVolumeOpenGL2ObjectFactory::New());
  //定义绘制器;
  vtkRenderer *aRenderer = vtkRenderer::New();//指向指针;
  vtkSmartPointer<vtkRenderWindow> renWin = vtkSmartPointer<vtkRenderWindow>::New();
  renWin->AddRenderer(aRenderer);
  vtkRenderWindowInteractor *iren = vtkRenderWindowInteractor::New();
  iren->SetRenderWindow(renWin);
  //读取数据;
  /*vtkDICOMImageReader *reader = vtkDICOMImageReader::New();
  reader->SetDirectoryName("E:/Coding Pra/VTK/VTK_Examples_StandardFormats_Input_DicomTestImages/DICOM");
  reader->SetDataByteOrderToLittleEndian();*/
  vtkSmartPointer<vtkDICOMImageReader> reader = vtkSmartPointer<vtkDICOMImageReader>::New();
  reader->SetDirectoryName("E:/DIcom_Data/DICOM");
  reader->SetDataByteOrderToLittleEndian();
  //图像数据预处理,类型转换:通过 vtkimageCast 将不同类型数据集转化为 vtk 可以处理的数据集;
  vtkImageCast *cast_file = vtkImageCast::New();
  cast_file->SetInputConnection(reader->GetOutputPort());
  cast_file->SetOutputScalarTypeToUnsignedShort();
  cast_file->Update();
  //透明度映射函数定义;
  vtkPiecewiseFunction *opacityTransform = vtkPiecewiseFunction::New();
  opacityTransform->AddPoint(0, 0.0);
  opacityTransform->AddPoint(20, 0.0);
  opacityTransform->AddPoint(200, 1.0);
  opacityTransform->AddPoint(300, 1.0);
  //颜色映射函数定义,梯度上升的
  vtkColorTransferFunction *colorTransformFunction = vtkColorTransferFunction::New();
  colorTransformFunction->AddRGBPoint(0.0, 0.0, 0.0, 0.0);
  colorTransformFunction->AddRGBPoint(64.0, 0.0, 0.0, 0.0);
  colorTransformFunction->AddRGBPoint(128.0, 1.0, 0.0, 0.0);
  colorTransformFunction->AddRGBPoint(192.0, 1.0, 0.0, 0.0);
  colorTransformFunction->AddRGBPoint(255.0, 1.0, 0.0, 0.0);
  vtkPiecewiseFunction *gradientTransform = vtkPiecewiseFunction::New();
  gradientTransform->AddPoint(0, 0.0);
  gradientTransform->AddPoint(20, 2.0);
  gradientTransform->AddPoint(200, 0.1);
  gradientTransform->AddPoint(300, 0.1);
  //体数据属性;
  vtkVolumeProperty *volumeProperty = vtkVolumeProperty::New();
  volumeProperty->SetColor(colorTransformFunction);
  volumeProperty->SetScalarOpacity(opacityTransform);
  volumeProperty->SetGradientOpacity(gradientTransform);
  volumeProperty->ShadeOn();//应用
  volumeProperty->SetInterpolationTypeToLinear();//直线间样条插值;
  volumeProperty->SetAmbient(0.4);//环境光系数;
  volumeProperty->SetDiffuse(0.6);//漫反射;
  volumeProperty->SetSpecular(0.2);
  volumeProperty->SetSpecularPower(10);//高光强度;
  计算光照效应;利用 vtkBolumeRayCaseMapper进行计算;
  //vtkVolumeRayCastMapper *volunemapper = vtkVolumeRayCastMapper::New();
  //vtkVolumeRayCastCompositeFunction *compositeFunction = vtkVolumeRayCastCompositeFunction::New();
  //光纤映射类型定义:
  vtkSmartPointer<vtkVolumeRayCastCompositeFunction> compositecast =
    vtkSmartPointer<vtkVolumeRayCastCompositeFunction>::New();
  //Mapper定义,
  vtkSmartPointer<vtkVolumeRayCastMapper> hiresMapper = 
    vtkSmartPointer<vtkVolumeRayCastMapper>::New();
  hiresMapper->SetInputData(cast_file->GetOutput());
  hiresMapper->SetVolumeRayCastFunction(compositecast);
  vtkSmartPointer<vtkLODProp3D> prop = vtkSmartPointer<vtkLODProp3D>::New();
  prop->AddLOD(hiresMapper,volumeProperty,0.0);
  //
  //volunemapper->SetVolumeRayCastFunction(compositeFunction);//载入体绘制方法;
  //volunemapper->SetInputConnection(cast_file->GetOutputPort());
  //vtkFixedPointVolumeRayCastMapper *fixedPointVolumeMapper = vtkFixedPointVolumeRayCastMapper::New()
  //fixedPointVolumeMapper->SetInput()
  vtkVolume *volume = vtkVolume::New();
  volume->SetMapper(hiresMapper);
  volume->SetProperty(volumeProperty);//设置体属性;
  double volumeView[4] = { 0,0,0.5,1 };
  vtkOutlineFilter *outlineData = vtkOutlineFilter::New();//线框;
  outlineData->SetInputConnection(reader->GetOutputPort());
  vtkPolyDataMapper *mapOutline = vtkPolyDataMapper::New();
  mapOutline->SetInputConnection(outlineData->GetOutputPort());
  vtkActor *outline = vtkActor::New();
  outline->SetMapper(mapOutline);
  outline->GetProperty()->SetColor(0, 0, 0);//背景纯黑色;
  aRenderer->AddVolume(volume);
  aRenderer->AddActor(outline);
  aRenderer->SetBackground(1, 1, 1);
  aRenderer->ResetCamera();
  //重设相机的剪切范围;
  aRenderer->ResetCameraClippingRange();
  renWin->SetSize(800, 800);
  renWin->SetWindowName("测试");
  vtkRenderWindowInteractor *iren2 = vtkRenderWindowInteractor::New();
  iren2->SetRenderWindow(renWin);
  //设置相机跟踪模式
  vtkInteractorStyleTrackballCamera *style = vtkInteractorStyleTrackballCamera::New();
  iren2->SetInteractorStyle(style);
  renWin->Render();
  iren2->Initialize();
  iren2->Start();
  vtkOBJExporter *porter = vtkOBJExporter::New();
  porter->SetFilePrefix("E:/ceshi/aaa/regist_after/esho.obj");
  porter->SetInput(renWin);
  porter->Write();
  porter->Update();
  return EXIT_SUCCESS;
}



上面是体绘制的结果,相对来说体绘制需要计算资源更大些, vtk 在这方面有所考虑,提供了vtKGPUVolumeRayCastMapper GUP 加速的光线投射算法。


相关文章
|
6天前
|
存储 数据采集 数据可视化
【C++】PACS医学图像存储和传输系统源码带三维重建
【C++】PACS医学图像存储和传输系统源码带三维重建
39 0
|
存储 编解码 搜索推荐
PACS系统源码(医学图像存储与传输)支持三维重建与还原
影像阅片 影像阅片是PACS最核心的部分,主要用来给医生提供调阅影像和影像处理,基础功能一般厂商都有,比如序列、旋转、放大缩小、标注、窗宽调整、四角信息设置、定位线、比例尺、测量、裁剪、伪彩等等,三维重建是一个亮点功能,很多厂商目前由于技术瓶颈尚未实现。这套PACS系统源码是带三维重建和还原的,是符合市场需求的PACS系统。
180 0
PACS系统源码(医学图像存储与传输)支持三维重建与还原
|
6天前
|
存储 编解码 安全
带三维重建和还原的PACS源码 医学影像PACS系统源码
带三维重建和还原的PACS源码 医学影像PACS系统源码 PACS及影像存取与传输系统”( Picture Archiving and Communication System),为以实现医学影像数字化存储、诊断为核心任务,从医学影像设备(如CT、CR、DR、MR、DSA、RF等)获取影像,集中存储、综合管理医学影像及病人相关信息,建立数字化工作流程。系统可实现检查预约、病人信息登记、计算机阅片、电子报告书写、胶片打印、数据备份等一系列满足影像科室日常工作的功能,并且由于影像数字化存储,用户可利用影像处理与测量技术辅助诊断、方便快捷地查找资料或利用网络将资料传输至临床科室,还可与医院HIS、L
22 0
|
6天前
|
存储 数据采集 固态存储
带三维重建和还原功能的医学影像管理系统(pacs)源码
带三维重建和还原功能的医学影像管理系统(pacs)源码
59 0
|
6天前
|
存储 数据可视化 vr&ar
突破传统 重新定义:3D医学影像PACS系统源码(包含RIS放射信息) 实现三维重建与还原
突破传统,重新定义PACS/RIS服务,洞察用户需求,关注应用场景,新一代PACS/RIS系统,系统顶层设计采用集中+分布式架构,满足医院影像全流程业务运行,同时各模块均可独立部署,满足医院未来影像信息化扩展新需求、感受新时代影像服务便捷性、易用性!系统基于平台化设计,与第三方服务自然接入无压力,从功能多样化到调阅速度快;覆盖(放射、超声、内镜、病理、核医学、心血管、临床科室等,是以影像采集、传输、存储、诊断、报告书写和科室管理)为核心应用的模块化PACS/RIS系统,实现了全院级影像信息的合理共享与应用。
23 0
突破传统 重新定义:3D医学影像PACS系统源码(包含RIS放射信息) 实现三维重建与还原
|
6天前
|
存储 数据采集 编解码
【PACS】医学影像管理系统源码带三维重建后处理技术
【PACS】医学影像管理系统源码带三维重建后处理技术
61 0
|
6天前
|
C++
【C++医学影像PACS】CT检查中的三维重建是什么检查?
【C++医学影像PACS】CT检查中的三维重建是什么检查?
70 0
|
6天前
|
数据采集 存储 数据可视化
医院影像PACS系统三维重建技术(获取数据、预处理、重建)
开放式体系结构,完全符合DICOM3.0标准,提供HL7标准接口,可实现与提供相应标准接口的HIS系统以及其他医学信息系统间的数据通信。
38 3
|
6天前
|
存储 编解码 监控
【C++】医学影像PACS三维重建后处理系统源码
系统完全符合国际标准的DICOM3.0标准
34 2
|
6天前
|
存储
医院PACS系统全套源码 强大的三维重建功能
对非DICOM影像,如超声、病理、心电图等进行了集成,做到了可以同时处理DICOM标准图像和非DICOM图像。
24 1