VTK修炼之道26：图像基本操作_三维图像切片提取

1.三维图像切片提取

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

下面是切片提取的代码：

#include <vtkAutoInit.h>
VTK_MODULE_INIT(vtkRenderingOpenGL);

#include <vtkSmartPointer.h>
#include <vtkImageData.h>
#include <vtkMetaImageReader.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 <vtkInteractorStyleImage.h>

int main(int argc, char* argv[])
{
	vtkSmartPointer<vtkMetaImageReader> reader =
		vtkSmartPointer<vtkMetaImageReader>::New();
	reader->SetFileName("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];
	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]);
	//*****************************************************************//
	static double axialElements[16] = {
		1, 0, 0, 0,
		0, 1, 0, 0,
		0, 0, 1, 0,
		0, 0, 0, 1
	};

	vtkSmartPointer<vtkMatrix4x4> 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]);

	vtkSmartPointer<vtkImageReslice> reslice =
		vtkSmartPointer<vtkImageReslice>::New();
	reslice->SetInputConnection(reader->GetOutputPort());
	reslice->SetOutputDimensionality(2);
	reslice->SetResliceAxes(resliceAxes);
	reslice->SetInterpolationModeToLinear();
	//*****************************************************************//
	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();
	vtkSmartPointer<vtkImageMapToColors> colorMap =
		vtkSmartPointer<vtkImageMapToColors>::New();
	colorMap->SetLookupTable(colorTable);
	colorMap->SetInputConnection(reslice->GetOutputPort());
	//*****************************************************************//
	vtkSmartPointer<vtkImageActor> imgActor =
		vtkSmartPointer<vtkImageActor>::New();
	imgActor->SetInputData(colorMap->GetOutput());

	vtkSmartPointer<vtkRenderer> renderer =
		vtkSmartPointer<vtkRenderer>::New();
	renderer->AddActor(imgActor);
	renderer->SetBackground(1.0, 1.0, 1.0);

	vtkSmartPointer<vtkRenderWindow> renderWindow =
		vtkSmartPointer<vtkRenderWindow>::New();
	renderWindow->AddRenderer(renderer);
	renderWindow->Render();
	renderWindow->SetSize(640, 480);
	renderWindow->SetWindowName("Extract3Dslice");

	vtkSmartPointer<vtkRenderWindowInteractor> rwi =
		vtkSmartPointer<vtkRenderWindowInteractor>::New();
	vtkSmartPointer<vtkInteractorStyleImage> imagestyle =
		vtkSmartPointer<vtkInteractorStyleImage>::New();
	rwi->SetInteractorStyle(imagestyle);
	rwi->SetRenderWindow(renderWindow);
	rwi->Initialize();
	rwi->Start();

	return 0;
}

首先通过vtkMetaImageReader读取一张医学三维图像，并获取得到图像范围（extent），原点和像素间隔；由这三个参数可以计算图像的中心位置center；接下来定义了切面的变换矩阵axialElements，该矩阵的前三列分别表示x、y和z方向向量，第四列为中心点坐标；

代码中的axialElements表示切面变换矩阵与当前坐标系一致，且切面为过中心点center，并平行于XY平面的平面？？？当前，定义该切面时，也可以是其他平面，甚至是任意平面，但是必须要过图像内部点。

下面给出了一个常用的变换矩阵。

提取平行于XZ平面的切片：

static double coronalElements[16] = {
 1, 0, 0, 0,
 0, 0, 1, 0,
 0,-1, 0, 0,
 0, 0, 0, 1 }; 

提取平行于YZ平面的切片：

static double sagittalElements[16] = {
0, 0,-1, 0,
1, 0, 0, 0,
0,-1, 0, 0,
0, 0, 0, 1 }; 

提取斜切切片：

static double obliqueElements[16] = {
1, 0, 0, 0,
0, 0.866025, -0.5, 0,
0, 0.5, 0.866025, 0,
0, 0, 0, 1 }; 

注意使用这些变换矩阵的时候，需要将第四列替换为切片经过图像的一个点坐标，上例中将图像的中心添加到axialElements矩阵，并通过函数SetResliceAxes设置变换矩阵，SetOutputDimensionality(2)指定输出的图像为一个二维图像； 而函数SetInterpolationModeToLinear()则指定了切面提取中的差值方式为线性差值，另外该类中还提供了其他的插值方式：
SetInterpolationModeToNearestNeighbor()：最近邻方式
SetInterpolationModeToCubic()：三次线性差值
设置完毕后，执行Update()即可完成切面计算。

东灵提供的预想结果应该是：

然而，在实际运行中，却碰到了vtkMetaImageReader不能读取文件的问题，在之前32bit系统上，这都是正常的，暂时还不能肯定问题出现在哪里，需要进一步研究！！问题如下：

2.参看资料

1.《C++ primer》
2.《The VTK User’s Guide – 11thEdition》
3.《The Visualization Toolkit – AnObject-Oriented Approach To 3D Graphics (4th Edition)》
4.  张晓东, 罗火灵. VTK图形图像开发进阶[M]. 机械工业出版社, 2015.