转载：VTK笔记-裁剪分割-三维曲线或几何切割体数据（黑山老妖）

https://blog.csdn.net/liushao1031177/article/details/118419221

捕鲸叉

1603人浏览 · 2022-06-26 13:58:04

捕鲸叉 · 2022-06-26 13:58:04 发布

学习整理了一下黑山老妖的博客《VTK笔记-裁剪分割-三维曲线或几何切割体数据》，使用阈值填充的方法获得切割体数据的结果。原文URL:https://blog.csdn.net/liushao1031177/article/details/118419221

本人测试的图形显示如下：

整理后的源代码如下：

#pragma once
#include <vtkNIFTIImageReader.h>
#include <vtkInformation.h>
#include <vtkVertexGlyphFilter.h>
#include <vtkSelectEnclosedPoints.h>
#include <vtkFloatArray.h>
#include <vtkShortArray.h>
#include <vtkActor.h>
#include <vtkDICOMImageReader.h>

#include <vtkPlane.h>
#include <vtkPolygon.h>
#include <vtkPolyDataMapper.h>
#include <vtkGPUVolumeRayCastMapper.h>
#include <vtkPolyDataToImageStencil.h>
#include <vtkImageData.h>
#include <vtkImageStencilToImage.h>
#include <vtkLineSource.h>
#include <vtkAutoInit.h>
#include <vtkDataSet.h>
#include <vtkPointData.h>
#include <vtkVolumeProperty.h>
#include <vtkPiecewiseFunction.h>
#include <vtkVolume.h>
#include <vtkColorTransferFunction.h>
#include <vtkRenderer.h>
#include <vtkRenderWindow.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkPlaneSource.h>
#include <vtkCamera.h>
#include <vtkProperty.h>

VTK_MODULE_INIT(vtkRenderingOpenGL2);
VTK_MODULE_INIT(vtkRenderingVolumeOpenGL2);
VTK_MODULE_INIT(vtkRenderingFreeType);
VTK_MODULE_INIT(vtkInteractionStyle);

double* spacing;
double* origin;
int* dm;

class ClipVolumeByMask
{
public:
	static void Test()
	{
		vtkNew<vtkRenderer> ren;
		vtkNew<vtkDICOMImageReader> reader;
		reader->SetDirectoryName("D:\\VTK-DATA\\Dicom-Samples\\xuebaolin");
		reader->Update();
		vtkImageData* pImageData = reader->GetOutput();
		dm = pImageData->GetDimensions();
		spacing = reader->GetDataSpacing();
		origin = reader->GetDataOrigin();

		int extent[6] = { 0, dm[0], 0, dm[1], 0, dm[2] - 1 };
		double camera_position[3] = { origin[0] + spacing[0] * dm[0] / 2, origin[1] - spacing[1] * dm[1] * 2 , origin[2] + spacing[2] * dm[2] / 2 };
		double top_left[5][3] =
		{
			{ camera_position[0] - 30,camera_position[1],camera_position[2] - 20 },
			{ camera_position[0] + 30,camera_position[1],camera_position[2] - 20 },
			{ camera_position[0] + 30,camera_position[1],camera_position[2] + 20 },
			{ camera_position[0] - 30,camera_position[1],camera_position[2] + 20 },
			{ 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> cylinderCutterBody = GenerateFivePrismatic(top_left, another_polygon);
		vtkSmartPointer<vtkPolyDataMapper> cutterMapper = vtkSmartPointer<vtkPolyDataMapper>::New();
		cutterMapper->SetInputData(cylinderCutterBody);
		vtkSmartPointer<vtkActor> cutterActor = vtkSmartPointer<vtkActor>::New();
		cutterActor->SetMapper(cutterMapper);
#pragma endregion

#pragma region 裁剪
		vtkNew<vtkImageData> whiteImage;

		whiteImage->SetSpacing(spacing);
		whiteImage->SetExtent(extent);
		whiteImage->SetOrigin(origin);
		whiteImage->AllocateScalars(VTK_UNSIGNED_CHAR, 1);

		vtkNew<vtkPolyDataToImageStencil> polyToStencil; // <--
		polyToStencil->SetInputData(cylinderCutterBody);
		polyToStencil->SetOutputOrigin(origin);
		polyToStencil->SetOutputSpacing(spacing);
		polyToStencil->SetOutputWholeExtent(whiteImage->GetExtent());
		polyToStencil->Update();

		const unsigned char inval = 0; //纯黑
		const unsigned char outval = 255; //纯白

		vtkNew<vtkImageStencilToImage> stencilToImage; // <--
		stencilToImage->SetInputConnection(polyToStencil->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(cutterActor);
		ren->AddActor(volumeActor);

		planeActor->SetVisibility(false);
		cutterActor->SetVisibility(false);

		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 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] = { spacing[0] * dm[0],0,0 };
		double dis2 = plane->DistanceToPlane(point_100);
		double point_010[3] = { 0,spacing[0] * dm[0],0 };
		double dis3 = plane->DistanceToPlane(point_010);
		double point_110[3] = { spacing[0] * dm[0],spacing[0] * dm[0],0 };
		double dis4 = plane->DistanceToPlane(point_110);
		double point_001[3] = { 0,0,(dm[2] - 1) };
		double dis5 = plane->DistanceToPlane(point_001);
		double point_101[3] = { spacing[0] * dm[0],0,(dm[2] - 1) };
		double dis6 = plane->DistanceToPlane(point_101);
		double point_011[3] = { 0,spacing[0] * dm[0],(dm[2] - 1) };
		double dis7 = plane->DistanceToPlane(point_011);
		double point_111[3] = { spacing[0] * dm[0],spacing[0] * dm[0],(dm[2] - 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(spacing[0] * dm[0], 0, 0);
		points_volume->InsertNextPoint(spacing[0] * dm[0], spacing[0] * dm[0], 0);
		points_volume->InsertNextPoint(0, spacing[0] * dm[0], 0);
		points_volume->InsertNextPoint(0, 0, (dm[2] - 1));
		points_volume->InsertNextPoint(spacing[0] * dm[0], 0, (dm[2] - 1));
		points_volume->InsertNextPoint(spacing[0] * dm[0], spacing[0] * dm[0], (dm[2] - 1));
		points_volume->InsertNextPoint(0, spacing[0] * dm[0], (dm[2] - 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;
	}
};

int main(int, char* [])
{
	ClipVolumeByMask::Test();
	return 0;
}

下图是我本人采用上述算法思想做的体绘制裁剪，效果基本出来了：