VTK笔记-裁剪分割-三维曲线或几何切割体数据

  先说重点，本文是按照开源代码MorphoDig实现的，大家可以直接去下载代码学习；

前言

  最近，在研究图像窗口上圈选范围切割三维体数据。网上大部分的博客都只列出了裁剪效果或者是大概的切割计算逻辑，没有给出具体的逻辑和实现方式。VTK资料并不多，而且网上大多数资料都是在重复，实在不能苟同这种“国之利器不可示于人”的想法。大家还是应该多交流多分享，让VTK广泛传播，给后人留下丰富的学习资料和例子，当从我做起。

医学图像重建

  医学图像三维重建技术是根据成像设备上的二维图像集合，利用计算机重建并合成病灶区或者人体器官的三维图像，可选择不同的重建算法使人体不同组织内部空间结构和物理属性被显示出来，直观展示出人体病灶及其周围组织器官的形状、空间分布等信息，方便医生对病变部位内部结构的可视化观察和诊断分析。随着医学图像三维重建技术的不断发展，医生可依据病人实际病情的具体诊断需要，重点三维显示可能的病变区域细节，可大大提高确诊率。医学图像三维重建技术已成为现阶段可视化领域的研究热点。

裁剪研究整理

  医学图像分割是可视化的基础，需要将某个器官组织区域分割出来，为进行可视化做预处理，而且分割结果的质量会直接影响三维可视化的应用。医学图像分割和可视化是两个不可分离的功能，对于医学诊断，手术模拟及人体解剖等方面具有重要的应用价值，是近年来医学领域的研究热点。这里的分割不同与裁剪操作，分割更注重根据图像中各个不同脏器组织在CT扫描后生成的不同数值进行一系列的算法计算，分割出针对某一脏器的分割算法。例如，有针对血管分割提取的算法，有针对骨骼的提取算法，在医学图像处理中有重要意义；
  交互式数据裁剪是在数据可视化后对数据进行交互式裁剪操作，是在图像分割后的补充；利用VTK等可视化工具中的裁剪工具切除体数据中选中的或者选中外的数据，保留想要观察诊断的数据。VTK中现在可以使用的裁剪方法有使用定义规则几何体（平面、圆柱体等多面体）来切割几何体，有使用不规则几何体和不规则几何体切割三维重建图像数据。然后衍生出其他从自身需求出发，通过上面两种途径实现裁剪体数据的目的。比如，使用在图像窗口绘制曲线对数据裁剪，使用圈选范围裁剪数据，使用规划规则几何体的移动路线裁剪体数据。
  对医学图像体数据及重构几何模型进行虚拟剖切，可以方便地看到内部的组织，便于观察和诊断，可用于医疗放射治疗规划。王卫红，秦绪佳，郑红波的《医学图像3维重建模型的虚拟剖切算法》提出了针对医学图像重建的表面几何模型，提出了对模型进行平面剖切，立体开窗及任意交互切割的算法。平面剖切和开窗是用剖切面或剖切体对重建模型施以剖切，在剖切面上生成边序列及顶点序列；由此边序列和顶点序列生成封闭的边界轮廓，确定各轮廓的包含关系；对封闭轮廓包围的截面区域进行Delaunay三角剖分，得到完整的剖切后的表面模型。任意交互切割过程是交互生成切割路径，确定切割边界，并沿切割边界对表面模型进行切割。
  李楚雅，杨勇，史天才提出了一种三维模型曲线切割方法，也算是一种在平面投影来判断是否在圈选范围内的算法；实际的操作如下：第一步,将三维体数据模型各像素依据当前场景投影为二维模型,计算出当前场景下,该三维模型中各像素在屏幕上的坐标：将三维体数据模型中各像素坐标,按照当前显示场景,投影到显示屏幕上,并得到相应屏幕坐标,所有像素的投影结果即作为二维投影模型；按照三维体数据中各像素点顺序将投影结果记录；第二步,采用种子生长法,以曲线为区域边界,扩展出封闭曲线内部和外部区域,并作不同标记：以曲线为边界,计算出曲线内部一点作为扩展种子点,基于该点,采用种子生长法,以曲线为区域边界,扩展出封闭曲线内部,并作标记,曲线上像素做另一标记,并将曲线各区域像素点以当前场景给予坐标值,形成一个二维区域图；第三步,将三维投影模型与封闭曲线区域进行对比,识别出在曲线内部和外部的三维投影模型部分：将三维投影模型中每一像素,依据投影后的坐标,在曲线区域图中,查找相应坐标像素状态,得出该像素点是在曲线内部还是外部或在曲线上,并为每一像素标记,标记其在曲线区域的哪一部分；第四步,将三维模型各像素投影所在区域应用于对应三维模型各像素中：依据投影顺序,将投影模型像素判断结果一一对应地反馈给原始三维模型。这篇论文的重点可能在于种子生长法来计算曲线内部的坐标结合，不过我没有获取到这篇论文，只能从摘要中分析一下；
  2015有篇论文《基于2D轮廓线的体数据切割算法》，作者李燕提出了一种基于2D轮廓曲线的交互式体数据切割算法，并结合可视化工具包和布尔操作实现单平面和包围盒的体数据切割。其中轮廓线体数据切割算法根据操作者在图像窗口任意勾画的2D轮廓曲线构建切割几何体，并将轮廓曲线和体数据通过坐标变换映射到掩膜图像内，在2D空间上实现三维图像数据的切割。论文的实验结果表明轮廓线切割算法可以处理任意形状的几何体。论文着重在归纳了当时的三维图像的裁剪方法介绍说明以及如何使用，他使用计算三维体数据投影在2D平面上的掩膜区域内来判断是否进行裁剪去除，自己的创新点在于如何生成掩膜图像，他采用了一种注入式反向填充算法将闭合不自交的曲线填充成为二值图像的掩膜图像，这里不详细讲述他这篇论文的实现，以后专门写篇一篇这个论文的阅读笔记；我自己按照他的论文实现了一下，论文的短板在于计算每个点在平面上的投影计算较为耗时，遍历整个三维体数据耗时更长，有博客按照这个论文实现后，耗时在30s以上，而我自己的实现耗时更长。
  博主Beyond欣实现的投影法进行裁剪VTK：利用qt实现体绘制剪裁，有需要的可以看一下他的实现代码。
  在本笔记之前，我写的几篇笔记都是在关于裁剪分割的。有的笔记是关于裁剪多边形数据，有的笔记是关于计算点在平面上的投影计算，还有的笔记是多个多边形形成闭合多边形数据，有的笔记是从从多边形数据转换为图像数据，其实这些笔记都是为了一步一步的走向使用闭合多边形数据制作图像掩膜数据对体数据进行分割的目的。

图像掩膜数据分割体数据

  下图是物体在不同平面上的投影，之前说到的平面投影法，就是体数据中的像素点坐标投影到平面上，判断是否在圈选范围内；
  如果我们在屏幕窗口上绘制一个四边形abcd，获取体数据8角坐标与屏幕屏幕上的投影距离得到最远的距离dis，通过平面的法向量可以距离屏幕屏幕距离dis的平面上的四边形ABCD的坐标，四边形abcd和四边形ABCD组成了一个四棱柱，由四边形abcd和四边ABCD组成四棱柱的上面和下面，中间棱柱有aAbBcCdDaA组成的三角带构成。之前笔记有验证过怎么从多个多边形组成闭合多边形数据见VTK笔记-几何数据-两个平面上的多边形形成闭合多边形。

  之前笔记VTK笔记-多边形数据转换图像数据-vtkPolyData转换为vtkImageData验证了如何将多边形数据转换为图像数据，在这里的功能就是将多边形形成的闭合多边形，通过设置起始坐标和三轴方向上的长度，获取到一个二值的图像数据，可以用作掩膜数据；之前笔记VTK笔记-裁剪分割-vtkClipVolume-vtkSelectPolyData类讲到了使用vtkClipVolume类裁剪体数据，不过我使用的VTK8.2版本中设置隐函数的尺寸较大时，就会报出内存分配失败的错误。有朋友推荐了一个GitHub的开源代码，我按照他的思路简单验证了一下，目测是可以实现，平面勾画2D轮廓线，通过创建Mask图像结合体渲染中的SetMask方法对数据进行裁剪；如果要重复对数据进行裁剪，那么就需要重复对Mask进行取交集，更新Mask数值，这个功能可以放在后面实现。现在就是实现设定一个闭合多边形，并从构建的闭合多形性所在平面的法向量计算出另外一个平行面上的多边形坐标，构成穿透体数据的闭合多边形。
  如下图中，在平面上绘制一个闭合的五边形，计算体数据的8个顶点投影在平面上的距离，计算出另外五个点的坐标，形成五棱的柱体，该柱体穿透了体数据，与体数据交集部分即是我们要裁剪的区域；



  注意：这里北部的五个点形成了一个四边形，其实应该是一个五边形，之前有笔记《VTK笔记-数据集-绘制凹多边形-多边形三角网格化-vtkTriangleFilter类》解决了该问题，有需要的可以跳过去看看。这里不再赘述。
  这里使用了vtkGPUVolumeRayCastMapper类提供的函数：
  SetMaskTypeToBinary和SetMaskTypeToLabelMap用来设置掩膜计算的类型，分为了Binary二值类型和LabelMap类型；
  如果掩膜是二进制MASK类型，则体绘制仅限于二进制遮罩内的区域。假设二进制掩码的数据类型为UCHAR，值为255（内部）和0（外部）。
  掩膜也可以是Label Map。标签映射的数据类型必须为UCHAR，即最多可以有256个标签。标签0保留为特殊标签。在标签值为0的体素中，使用vtkVolumeProperty提供的默认传递函数。对于标签值大于0的体素，将使用使用vtkVolumeProperty的标签API提供的颜色传递函数。对于标签值大于0的体素，颜色传递函数与默认颜色传递函数混合，混合权重由MaskBlendFactor确定。

void SetMaskInput(vtkImageData* mask);
vtkGetObjectMacro(MaskInput, vtkImageData);
enum
{
	BinaryMaskType = 0,
	LabelMapMaskType
};
vtkSetMacro(MaskType, int);
vtkGetMacro(MaskType, int);
void SetMaskTypeToBinary();
void SetMaskTypeToLabelMap();

  通过接口SetMaskInput和GetMaskInput来设置或者获取掩膜图像数据；

void SetMaskInput(vtkImageData* mask);
vtkGetObjectMacro(MaskInput, vtkImageData);

  如果不使用vtkGPUVolumeRayCastMapper类，那么也可以使用Mask掩膜数据对原始体数据进行处理，或者是进行掩膜计算，使用vtkImageMask，可以参考之前的笔记《VTK笔记-图像融合-vtkImageBlend-vtkImageMask》；

vtkSmartPointer<vtkImageMask> mask_filter = vtkSmartPointer<vtkImageMask>::New();
mask_filter->SetInputData(0, pImageData);
mask_filter->SetInputData(1, stencilToImage->GetOutput());
mask_filter->SetMaskedOutputValue(-0, -0, -0);		
mask_filter->Update();

  mask_filter->GetOutput()就是掩膜处理后的图像数据；

示例

  由于只是做一个验证性的demo，所以使用了一个半身DICOM序列，数据的长宽高都是固定给出的，在视平面上固定五个点坐标和法向量是之前根据相机的自动调整获取到的，绘制了视平面、五棱柱体和体数据的外边框；
  计算步骤如下：
    1.加载DICOM图像，DICOM信息，包括了XYZ三方向上的尺寸和像素间距等信息；生成vtkImageData体数据；
    2.根据五个顶点坐标和法向量，构建五边形和视平面plane；
    3.根据体数据的8个顶点，计算出距离plane最远的距离；根据该距离划定另外一个平面上的五个点的坐标；
    4.构建闭合的五棱柱体；并从vtkPolyData转换为vtkImageData掩膜数据；
    5.使用vtkGPUVolumeRayCastMapper设置SetMaskData掩膜数据，设置掩膜类型为二值类型；
    6.体渲染流程，绘制其他线，点，面等几何数据；
完整代码如下：

#pragma once
#include "vtkNIFTIImageReader.h"
#include "vtkInformation.h"
#include <vtkVertexGlyphFilter.h>
#include <vtkSelectEnclosedPoints.h>
#include <vtkFloatArray.h>
#include "dcmtk\config\osconfig.h"
#include "dcmtk\dcmdata\dctk.h"
#include <dcmtk/dcmjpeg/djdecode.h>  /* for dcmjpeg decoders */
#include <dcmtk/dcmjpeg/djencode.h>
#include <dcmtk/dcmjpls/djdecode.h>		//for JPEG-LS decode
#include <dcmtk/dcmjpls/djencode.h>		//for JPEG-LS encode
#include <vtkPlane.h>
#include <vtkPolygon.h>
#include <vtkGPUVolumeRayCastMapper.h>
#include <vtkPolyDataToImageStencil.h>
#include <vtkImageStencilToImage.h>
#include <vtkLineSource.h>
#include "vtkAutoInit.h" 
VTK_MODULE_INIT(vtkRenderingOpenGL2);
VTK_MODULE_INIT(vtkRenderingVolumeOpenGL2);
VTK_MODULE_INIT(vtkRenderingFreeType);
VTK_MODULE_INIT(vtkInteractionStyle);
class ClipVolumeByMask {
public:
	static void Test() {
		vtkNew<vtkRenderer> ren;

		DJDecoderRegistration::registerCodecs();
		DJLSEncoderRegistration::registerCodecs();		//JPEG-LS encoder registerCodecs
		DJLSDecoderRegistration::registerCodecs();		//JPEG-LS decoder registerCodecs

		int slice_num = 438;
		int single_img = 512 * 512;
		int size_all = slice_num * single_img;
		int extent[6] = { 0,511,0,511,0,slice_num - 1 };
		double spacing[3] = { 0.586, 0.586, 1.0 };
		double origin[3] = { 0 };
#pragma region 1.从DICOM文件中获取像素，生成vtkImageData
		short* pBuf = new short[size_all] {0};
		Read_Dicom_File(pBuf, slice_num, single_img);

		vtkNew<vtkShortArray> dataArray;
		dataArray->SetArray(pBuf, size_all, 1);
		vtkNew<vtkImageData> pImageData;

		pImageData->AllocateScalars(VTK_SHORT, 1);
		pImageData->SetDimensions(512, 512, slice_num);
		pImageData->GetPointData()->SetScalars(dataArray);
		pImageData->SetSpacing(spacing);
		pImageData->SetOrigin(origin);
#pragma endregion

#pragma region 闭合多边形
		double camera_position[3] = { 149.72299999999998, -1025.8604209057046, 218.50000000000000 };
		double top_left[5][3] = {
			{ camera_position[0] - 20,camera_position[1],camera_position[2] - 10 },
			{ camera_position[0] + 20,camera_position[1],camera_position[2] - 10 },
			{ camera_position[0] + 20,camera_position[1],camera_position[2] + 10 }, 
			{ camera_position[0] - 20,camera_position[1],camera_position[2] + 10 },
			{ camera_position[0],camera_position[1],camera_position[2]} };

		double normal_camera[3] = { 0.0000000000000000, -1.0000000000000000, 0.0000000000000000 };

		double another_polygon[5][3] = { 0 };
		Get_Another_Plane_Polygon(normal_camera, top_left, another_polygon);

#pragma region 生成五棱柱体
		vtkSmartPointer<vtkPolyData> geometry = GenerateFivePrismatic(top_left, another_polygon);		
		vtkSmartPointer<vtkPolyDataMapper> geometryMapper = vtkSmartPointer<vtkPolyDataMapper>::New();
		geometryMapper->SetInputData(geometry);
		vtkSmartPointer<vtkActor> geometryActor = vtkSmartPointer<vtkActor>::New();
		geometryActor->SetMapper(geometryMapper);
#pragma endregion
		
#pragma region 裁剪
		vtkNew<vtkImageData> whiteImage;

		whiteImage->SetSpacing(spacing);
		whiteImage->SetExtent(extent);
		whiteImage->SetOrigin(origin);
		whiteImage->AllocateScalars(VTK_UNSIGNED_CHAR, 1);
		vtkNew<vtkPolyDataToImageStencil> pol2stenc;
		pol2stenc->SetInputData(geometry);
		pol2stenc->SetOutputOrigin(origin);
		pol2stenc->SetOutputSpacing(spacing);
		pol2stenc->SetOutputWholeExtent(whiteImage->GetExtent());
		pol2stenc->Update();
		const unsigned char inval = 0;
		const unsigned char outval = 255;
		vtkNew<vtkImageStencilToImage> stencilToImage;
		stencilToImage->SetInputConnection(pol2stenc->GetOutputPort());
		stencilToImage->SetInsideValue(inval);
		stencilToImage->SetOutsideValue(outval);
		
		stencilToImage->SetOutputScalarType(VTK_UNSIGNED_CHAR);
		stencilToImage->Update();
#pragma endregion
		// 绘制两个平面上的五个点，已经前面的五边形
		vtkSmartPointer<vtkActor> pointsActor = DrawPolyFivePoint(camera_position, top_left, another_polygon);
		vtkSmartPointer<vtkActor> linesActor = DrawPolyFiveLine(camera_position, top_left, another_polygon);
		// 绘制视平面
		vtkSmartPointer<vtkActor> planeActor = DrawViewPlane(normal_camera, camera_position);
		vtkSmartPointer<vtkActor> volumeActor = DrawOutline();
#pragma endregion

#pragma region 使用vtkGPUVolumeRayCastMapper体渲染方式
		vtkNew<vtkGPUVolumeRayCastMapper> gpuMapper;
		gpuMapper->SetInputData(pImageData);
		gpuMapper->SetMaskTypeToBinary();
		gpuMapper->SetMaskInput(stencilToImage->GetOutput());
		
		vtkNew<vtkVolumeProperty> volumeProperty;
		volumeProperty->SetInterpolationTypeToLinear();
		volumeProperty->SetAmbient(0.4);
		volumeProperty->SetDiffuse(1.6);  //漫反射
		volumeProperty->SetSpecular(0.2); //镜面反射
		volumeProperty->ShadeOn();  //打开或者关闭阴影测试
		//设置不透明度
		vtkNew<vtkPiecewiseFunction> compositeOpacity;
		compositeOpacity->AddPoint(220, 0.00);
		compositeOpacity->AddPoint(255, 1.0);
		volumeProperty->SetScalarOpacity(compositeOpacity); //设置不透明度传输函数

		vtkNew<vtkColorTransferFunction> color;
		color->AddRGBPoint(0.000, 0.00, 0.00, 0.00);
		color->AddRGBPoint(60, 1.00, 0.52, 0.30);
		color->AddRGBPoint(150, 1.00, 1.00, 1.00);
		color->AddRGBPoint(200, .5, .3, .2);
		volumeProperty->SetColor(color);

		vtkNew<vtkVolume> volume;
		volume->SetMapper(gpuMapper);
		volume->SetProperty(volumeProperty);
#pragma endregion
#pragma region 渲染管线		
		ren->SetBackground(0.3, 0.4, 0.3);
		ren->AddVolume(volume);
		ren->GetActiveCamera()->SetViewUp(0, 0, 1);
		ren->GetActiveCamera()->SetFocalPoint(0, 0, 0);
		ren->GetActiveCamera()->SetPosition(0, -1, 0);
		ren->ResetCamera();
		
		ren->AddActor(planeActor);
		ren->AddActor(linesActor);
		ren->AddActor(pointsActor);
		ren->AddActor(geometryActor);
		ren->AddActor(volumeActor);
		vtkNew<vtkRenderWindow> rw;
		rw->AddRenderer(ren);
		rw->SetSize(640, 480);
		rw->Render();
		rw->SetWindowName("VolumeRendering PipeLine");

		vtkNew<vtkRenderWindowInteractor> rwi;
		rwi->SetRenderWindow(rw);

		rw->Render();
		rwi->Start();
#pragma endregion
	}
private:
	static void Read_Dicom_File(short* pBuf, int slice_num, int szie_single_img) {
		for (int i = 1; i <= slice_num; i++){
			std::string fileName = "G:\\Data\\" + std::to_string(i) + ".DCM";

			DcmFileFormat *myFileFormat = new DcmFileFormat;
			OFCondition cond = myFileFormat->loadFile(fileName.c_str());

			DcmElement* pElement = nullptr;
			myFileFormat->getDataset()->chooseRepresentation(EXS_DeflatedLittleEndianExplicit, NULL);
			myFileFormat->getDataset()->findAndGetElement(DCM_PixelData, pElement);
			if (pElement != NULL) {
				int size = pElement->getLength() / 2;
				OFCondition cond = pElement->loadAllDataIntoMemory();
				if (cond.good()) {
					Uint16 * ptr;
					cond = pElement->getUint16Array(ptr);
					if (cond.good()) {
						memcpy(pBuf + (i - 1)*szie_single_img, ptr, size * sizeof(short));
					}
				}
			}
			delete myFileFormat;
		}
	}
	static void Get_Another_Plane_Polygon(double normal[3], double polygon[5][3], double another_polygon[5][3])	{
		double max_distance = Cal_Max_Distance(normal, polygon[0]);
		for (int i = 0; i < 5; i++){
			another_polygon[i][0] = polygon[i][0] - max_distance * normal[0];
			another_polygon[i][1] = polygon[i][1] - max_distance * normal[1];
			another_polygon[i][2] = polygon[i][2] - max_distance * normal[2];
		}
	}

	static double Cal_Max_Distance(double normal[3], double origin[3])	{
		vtkNew<vtkPlane> plane;
		plane->SetNormal(normal);
		plane->SetOrigin(origin);
		double point_000[3] = { 0,0,0 };
		double dis1 = plane->DistanceToPlane(point_000);
		double point_100[3] = { 0.586 * 511,0,0 };
		double dis2 = plane->DistanceToPlane(point_100);
		double point_010[3] = { 0,0.586 * 511,0 };
		double dis3 = plane->DistanceToPlane(point_010);
		double point_110[3] = { 0.586 * 511,0.586 * 511,0 };
		double dis4 = plane->DistanceToPlane(point_110);
		double point_001[3] = { 0,0,(438 - 1) };
		double dis5 = plane->DistanceToPlane(point_001);
		double point_101[3] = { 0.586 * 511,0,(438 - 1) };
		double dis6 = plane->DistanceToPlane(point_101);
		double point_011[3] = { 0,0.586 * 511,(438 - 1) };
		double dis7 = plane->DistanceToPlane(point_011);
		double point_111[3] = { 0.586 * 511,0.586 * 511,(438 - 1) };
		double dis8 = plane->DistanceToPlane(point_111);
		return dis8;
	}
	static vtkSmartPointer<vtkActor> DrawViewPlane(double normal[3], double center[3])	{
		vtkSmartPointer<vtkPlaneSource> plane = vtkSmartPointer<vtkPlaneSource>::New();
		plane->SetNormal(normal);
		plane->SetCenter(center);
		plane->SetOrigin(center[0] - 100, center[1], center[2] - 100);
		plane->SetPoint1(center[0] + 100, center[1], center[2] - 100);
		plane->SetPoint2(center[0] - 100, center[1], center[2] + 100);
		vtkSmartPointer<vtkPolyDataMapper> planeMapper = vtkSmartPointer<vtkPolyDataMapper>::New();
		planeMapper->SetInputConnection(plane->GetOutputPort());
		vtkSmartPointer<vtkActor> planeActor = vtkSmartPointer<vtkActor>::New();
		planeActor->SetMapper(planeMapper);
		return planeActor;
	}
	static vtkSmartPointer<vtkActor> DrawPolyFivePoint(double center[3],double top_left[4][3],double another_polygon[4][3]) {
		vtkSmartPointer<vtkPoints> points1 = vtkSmartPointer<vtkPoints>::New();
		points1->InsertNextPoint(center);
		points1->InsertNextPoint(top_left[0]);
		points1->InsertNextPoint(top_left[1]);
		points1->InsertNextPoint(top_left[2]);
		points1->InsertNextPoint(top_left[3]);
		points1->InsertNextPoint(another_polygon[0]);
		points1->InsertNextPoint(another_polygon[1]);
		points1->InsertNextPoint(another_polygon[2]);
		points1->InsertNextPoint(another_polygon[3]);

		vtkSmartPointer<vtkPolyData> pointsPolydata = vtkSmartPointer<vtkPolyData>::New();
		pointsPolydata->SetPoints(points1);

		vtkSmartPointer<vtkVertexGlyphFilter> vertexGlyphFilter = vtkSmartPointer<vtkVertexGlyphFilter>::New();
		vertexGlyphFilter->AddInputData(pointsPolydata);
		vertexGlyphFilter->Update();

		vtkSmartPointer<vtkPolyDataMapper> pointsMapper = vtkSmartPointer<vtkPolyDataMapper>::New();
		pointsMapper->SetInputConnection(vertexGlyphFilter->GetOutputPort());

		vtkSmartPointer<vtkActor> pointsActor = vtkSmartPointer<vtkActor>::New();

		pointsActor->SetMapper(pointsMapper);
		pointsActor->GetProperty()->SetPointSize(10);//定义点的尺寸大小，这样点才能在画布上显示出来
		pointsActor->GetProperty()->SetColor(0.0, 0.0, 1);

		return pointsActor;
	}

	static vtkSmartPointer<vtkActor> DrawPolyFiveLine(double center[3], double top_left[4][3], double another_polygon[4][3]) {
		vtkSmartPointer<vtkPoints> points1 = vtkSmartPointer<vtkPoints>::New();
		points1->InsertNextPoint(center);
		points1->InsertNextPoint(top_left[0]);
		points1->InsertNextPoint(top_left[1]);
		points1->InsertNextPoint(top_left[2]);
		points1->InsertNextPoint(top_left[3]);
		vtkSmartPointer<vtkCellArray> line_polys = vtkSmartPointer<vtkCellArray>::New();
		line_polys->InsertNextCell(6);

		for (size_t i = 0; i < 5; i++)
		{
			line_polys->InsertCellPoint(i);
		}
		line_polys->InsertCellPoint(0);
		vtkSmartPointer<vtkPolyData> geometry_line = vtkSmartPointer<vtkPolyData>::New();
		geometry_line->SetPoints(points1);
		geometry_line->SetLines(line_polys);
		vtkSmartPointer<vtkPolyDataMapper> linesMapper = vtkSmartPointer<vtkPolyDataMapper>::New();
		linesMapper->SetInputData(geometry_line);

		vtkSmartPointer<vtkActor> linesActor = vtkSmartPointer<vtkActor>::New();

		linesActor->SetMapper(linesMapper);
		linesActor->GetProperty()->SetColor(255, 0, 0);
		return linesActor;
	}
	static vtkSmartPointer<vtkPolyData> GenerateFivePrismatic(double top_left[4][3], double another_polygon[4][3]) {
		vtkSmartPointer<vtkPolyData> geometry = vtkSmartPointer<vtkPolyData>::New();
		vtkSmartPointer<vtkPoints> points = vtkSmartPointer<vtkPoints>::New();
		for (size_t i = 0; i < 5; i++) {
			points->InsertPoint(2 * i, top_left[i]);
			points->InsertPoint(2 * i + 1, another_polygon[i]);
		}
		vtkSmartPointer<vtkCellArray> polys = vtkSmartPointer<vtkCellArray>::New();
		vtkSmartPointer<vtkCellArray> strips = vtkSmartPointer<vtkCellArray>::New();
		// 前平面
		polys->InsertNextCell(5);
		for (int i = 0; i < 5; i++) {
			polys->InsertCellPoint(2 * i);
		}
		// 后平面
		polys->InsertNextCell(5);
		for (int i = 0; i < 5; i++) {
			polys->InsertCellPoint(2 * i + 1);
		}
		// 中间柱面
		strips->InsertNextCell(10 + 2);
		for (int i = 0; i < 10; i++) {
			strips->InsertCellPoint(i);
		}
		strips->InsertCellPoint(0);
		strips->InsertCellPoint(1);
		geometry->SetPoints(points);
		geometry->SetPolys(polys);
		geometry->SetStrips(strips);
		return geometry;
	}		
	static vtkSmartPointer<vtkActor> DrawOutline() {
		vtkSmartPointer<vtkPoints> points_volume = vtkSmartPointer<vtkPoints>::New();
		points_volume->InsertNextPoint(0, 0, 0);
		points_volume->InsertNextPoint(0.586 * 511, 0, 0);
		points_volume->InsertNextPoint(0.586 * 511, 0.586 * 511, 0);
		points_volume->InsertNextPoint(0, 0.586 * 511, 0);
		points_volume->InsertNextPoint(0, 0, (438 - 1));
		points_volume->InsertNextPoint(0.586 * 511, 0, (438 - 1));
		points_volume->InsertNextPoint(0.586 * 511, 0.586 * 511, (438 - 1));
		points_volume->InsertNextPoint(0, 0.586 * 511, (438 - 1));

		vtkSmartPointer<vtkCellArray> line_volume = vtkSmartPointer<vtkCellArray>::New();
		line_volume->InsertNextCell(5);
		for (int i = 0; i < 4; i++) {
			line_volume->InsertCellPoint(i);
		}
		line_volume->InsertCellPoint(0);
		line_volume->InsertNextCell(5);
		for (int i = 4; i < 8; i++) {
			line_volume->InsertCellPoint(i);
		}
		line_volume->InsertCellPoint(4);

		for (int i = 0; i < 4; i++) {
			line_volume->InsertNextCell(2);
			line_volume->InsertCellPoint(0 + i);
			line_volume->InsertCellPoint(4 + i);
		}

		vtkSmartPointer<vtkPolyData> geometry_volume = vtkSmartPointer<vtkPolyData>::New();
		geometry_volume->SetPoints(points_volume);
		geometry_volume->SetLines(line_volume);
		vtkSmartPointer<vtkPolyDataMapper> volumeMapper = vtkSmartPointer<vtkPolyDataMapper>::New();
		volumeMapper->SetInputData(geometry_volume);

		vtkSmartPointer<vtkActor> volumeActor = vtkSmartPointer<vtkActor>::New();
		volumeActor->SetMapper(volumeMapper);
		return volumeActor;
	}
};

  最后的显示效果为，左图为裁剪内部数据留下外部数据，右图为裁剪外部数据留下内部数据；

遗留一个问题是：裁剪体数据后，边缘有锯齿，如何解决还在寻找方法。

参考资料

1.vtk的几种分割算法总结
2.利用VTK手动切割体数据，曲面切割
3.vtk体绘制的任意切割
4.vtkGPUVolumeRayCastMapper -> SetMaskInput : vtkVolumePicker picks masked cells