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

2023-10-29

关注公众号：【小张Python】，为你准备了 50+ 本Python 精品电子书籍与 50G + 优质视频学习资料，后台回复关键字：1024 即可获取；如果对博文内容有什么疑问，后台添加作者【个人微信】，可与作者直接进行交流

序言，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 = 256 28=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() 设置计算等值面的法向量，提高渲染质量；

未命名文件 (2).png

上面这张图显示的就是 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 加速的光线投射算法。

以上就是本篇文章的全部内容，最后感谢阅读！

欢迎交流与联系，文章首发于公众号 (小张Python)

Reference:

https://blog.csdn.net/wp_veil/article/details/7047537;

https://blog.csdn.net/www_doling_net/article/details/44960713

本文内容由网友自发贡献，版权归原作者所有，本站不承担相应法律责任。如您发现有涉嫌抄袭侵权的内容，请联系:hwhale#tublm.com(使用前将#替换为@)