VTK Learning Seven

医疗图像显示之切片提取

引言

切片是指三维图像中的一个切面对应的图像。切面可以是过图像内部一点且平行于XY、YZ、XZ平面的平面，也可以是任意的过三维图像内部一点任意方向的平面。通过提取切片可以方便的浏览和分析图像内部组织结构，是医学图像浏览软件中的一个重要的功能。在VTK中vtkImageReslice类实现图像切片提取功能。

Code

#include <vtkSmartPointer.h>
#include <vtkImageReader2.h>
#include <vtkMatrix4x4.h>
#include <vtkImageReslice.h>
#include <vtkLookupTable.h>
#include <vtkImageMapToColors.h>
#include <vtkImageActor.h>
#include <vtkRenderer.h>
#include <vtkRenderWindow.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkImageData.h>
#include <vtkMetaImageReader.h>
#include <vtkInteractorStyleImage.h>
#include <vector>

int main()
{
	vtkSmartPointer<vtkMetaImageReader> reader =
		vtkSmartPointer<vtkMetaImageReader>::New();
	reader->SetFileName("D:\\VTK\\demo\\ImageShow\\brain.mhd");
	reader->Update();

	int extent[6];
	double spacing[3];
	double origin[3];

	reader->GetOutput()->GetExtent(extent);
	reader->GetOutput()->GetSpacing(spacing);
	reader->GetOutput()->GetOrigin(origin);

	//计算中心点位置
	double center[3] = {0};
	center[0] = origin[0] + spacing[0] * 0.5 * (extent[0] + extent[1]);
	center[1] = origin[1] + spacing[1] * 0.5 * (extent[2] + extent[3]);
	center[2] = origin[2] + spacing[2] * 0.5 * (extent[4] + extent[5]);
	//横断面，平行于XY 平面
	static double axialElements[16] = {
		1, 0, 0, 0,
		0, 1, 0, 0,
		0, 0, 1, 0,
		0, 0, 0, 1
	};
	//由横断面绕x轴旋转90° ,冠状面，平行于XZ平面
	static double coronalElements[16] = {
		1, 0,0,0,
		0, 0,1,0,
		0,-1,0,0,
		0, 0,0,1
	};
	//由横断面绕x轴旋转90°，再绕y轴旋转90°，矢状面，平行于YZ平面
	static double sagittalElement[16] = {
		0, 0,-1,0,
		1, 0, 0,0,
		0,-1, 0,0,
		0, 0, 0,1
	};
	//由横断面绕x旋转30°
	static double obliqueElement[16] = {
		1,       0,       0,0,
		0,0.866025,    -0.5,0,
		0,	   0.5,0.866025,0,
		0,       0,       0,1
	};

	std::vector<vtkSmartPointer<vtkImageReslice>> reslices;

	auto resliceAxes = vtkSmartPointer<vtkMatrix4x4>::New();
	resliceAxes->DeepCopy(axialElements);
	resliceAxes->SetElement(0, 3, center[0]);
	resliceAxes->SetElement(1, 3, center[1]);
	resliceAxes->SetElement(2, 3, center[2]);

	auto reslice = vtkSmartPointer<vtkImageReslice>::New();
	reslice->SetInputConnection(reader->GetOutputPort());
	reslice->SetOutputDimensionality(2);
	reslice->SetResliceAxes(resliceAxes);
	reslice->SetInterpolationModeToLinear();
	reslices.push_back(reslice);

	auto resliceAxes1 = vtkSmartPointer<vtkMatrix4x4>::New();
	resliceAxes1->DeepCopy(coronalElements);
	resliceAxes1->SetElement(0, 3, center[0]);
	resliceAxes1->SetElement(1, 3, center[1]);
	resliceAxes1->SetElement(2, 3, center[2]);

	auto reslice1 = vtkSmartPointer<vtkImageReslice>::New();
	reslice1->SetInputConnection(reader->GetOutputPort());
	reslice1->SetOutputDimensionality(2);
	reslice1->SetResliceAxes(resliceAxes1);
	reslice1->SetInterpolationModeToLinear();
	reslices.push_back(reslice1);

	auto resliceAxes2 = vtkSmartPointer<vtkMatrix4x4>::New();
	resliceAxes2->DeepCopy(sagittalElement);
	resliceAxes2->SetElement(0, 3, center[0]);
	resliceAxes2->SetElement(1, 3, center[1]);
	resliceAxes2->SetElement(2, 3, center[2]);

	auto reslice2 = vtkSmartPointer<vtkImageReslice>::New();
	reslice2->SetInputConnection(reader->GetOutputPort());
	reslice2->SetOutputDimensionality(2);
	reslice2->SetResliceAxes(resliceAxes2);
	reslice2->SetInterpolationModeToLinear();
	reslices.push_back(reslice2);

	auto resliceAxes3 = vtkSmartPointer<vtkMatrix4x4>::New();
	resliceAxes3->DeepCopy(obliqueElement);
	resliceAxes3->SetElement(0, 3, center[0]);
	resliceAxes3->SetElement(1, 3, center[1]);
	resliceAxes3->SetElement(2, 3, center[2]);

	auto reslice3 = vtkSmartPointer<vtkImageReslice>::New();
	reslice3->SetInputConnection(reader->GetOutputPort());
	reslice3->SetOutputDimensionality(2);
	reslice3->SetResliceAxes(resliceAxes3);
	reslice3->SetInterpolationModeToLinear();
	reslices.push_back(reslice3);

	//生成颜色映射
	vtkSmartPointer<vtkLookupTable> colorTable =
		vtkSmartPointer<vtkLookupTable>::New();
	colorTable->SetRange(0, 1000);
	colorTable->SetValueRange(0.0, 1.0);
	colorTable->SetSaturationRange(0.0, 0.0);
	colorTable->SetRampToLinear();
	colorTable->Build();
	
	// Rows and columns
	unsigned int gridXDimensions = 2;
	unsigned int gridYDimensions = 2;

	// Need a renderer even if there is no actor
	std::vector<vtkSmartPointer<vtkRenderer>> renderers;
	for (size_t i = 0; i < gridXDimensions * gridYDimensions; i++)
	{
		renderers.push_back(vtkSmartPointer<vtkRenderer>::New());
		renderers[i]->SetBackground(1, 1, 1);
	}

	auto renderWindow = vtkSmartPointer<vtkRenderWindow>::New();
	int rendererSize = 256;
	renderWindow->SetSize(rendererSize * gridXDimensions,
		rendererSize * gridYDimensions);

	vtkSmartPointer<vtkRenderWindowInteractor> renderWindowInteractor =
		vtkSmartPointer<vtkRenderWindowInteractor>::New();

	// Create a pipeline for each function
	for (int row = 0; row < static_cast<int>(gridYDimensions); row++)
	{
		for (int col = 0; col < static_cast<int>(gridXDimensions); col++)
		{
			int index = row * gridXDimensions + col;

			// (xmin, ymin, xmax, ymax)
			double viewport[4] = { static_cast<double>(col) * rendererSize /
									  (gridXDimensions * rendererSize),

								  static_cast<double>(gridYDimensions - (row + 1)) *
									  rendererSize / (gridYDimensions * rendererSize),

								  static_cast<double>(col + 1) * rendererSize /
									  (gridXDimensions * rendererSize),

								  static_cast<double>(gridYDimensions - row) *
									  rendererSize /
									  (gridYDimensions * rendererSize) };

			renderWindow->AddRenderer(renderers[index]);
			renderers[index]->SetViewport(viewport);
			// Define the pipeline
			vtkSmartPointer<vtkImageMapToColors>mapcolors = vtkSmartPointer<vtkImageMapToColors>::New();
			mapcolors->SetLookupTable(colorTable);
			mapcolors->SetInputConnection(reslices[index]->GetOutputPort());
			mapcolors->Update();
			vtkSmartPointer<vtkImageActor> imgActor = vtkSmartPointer<vtkImageActor>::New();
			imgActor->SetInputData(mapcolors->GetOutput());
			renderers[index]->AddActor(imgActor);
		}
	}
	vtkSmartPointer<vtkInteractorStyleImage> imagestyle =
		vtkSmartPointer<vtkInteractorStyleImage>::New();
	renderWindowInteractor->SetInteractorStyle(imagestyle);
	renderWindowInteractor->SetRenderWindow(renderWindow);
	renderWindowInteractor->Initialize();
	renderWindowInteractor->Start();
	return EXIT_SUCCESS;
}

输出

横断面，冠状面，矢状面，斜切面