VTK cultivation 26: basic operation of image_ Three dimensional image slice extraction

1. 3D image slice extraction

A slice is an image corresponding to a section of a three-dimensional image. The section can be a plane that passes through an inner point of the image and is parallel to the XY, YZ and XZ planes, or it can be any plane that passes through an inner point of the three-dimensional image in any direction. It is an important function of medical image browsing software to browse and analyze the internal tissue structure of images easily by extracting slices. In VTK, VTKImageReslice class realizes the function of image slice extraction.
Here is the code for slice extraction:

#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;
}

Firstly, a medical 3D image is read through VTK Metaimagereader, and the extent, origin and pixel interval of the image are obtained. From these three parameters, the center position of the image can be calculated. Next, we define the transformation matrix AxialElements. The first three columns of the matrix represent the x, y, and z vectors, and the fourth column represents the center coordinates.
AxialElements in the code represents a plane whose plane transformation matrix is consistent with the current coordinate system, whose plane is parallel to the XY plane, and whose plane is past the center point. Now, when you define that slice, it could be any plane, any plane, but you have to go through the interior of the image.
A common transformation matrix is given below.
Extraction of slices parallel to the XZ plane:

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

Extract slices parallel to the YZ plane:

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

Extract oblique slice:

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

Note that when using these transformation matrices, you need to replace the fourth column with a point coordinate of the image that the slice passes through. In the example above, add the center of the image to the AXIALElements matrix and set the transformation matrix with the function SetResliceAxes. SetOutputDimensionality(2) specifies that the output image is a two-dimensional image.
And function SetInterpolationModeToLinear () specifies the difference of edge extraction method for the linear difference, the class also provides other interpolation methods:

SetInterpolationModeToNearestNeighbor () : nearest neighbor approach

SetInterpolationModeToCubic () : three linear difference

After the setup is complete, execute Update() to complete the section calculation.

The expected results provided by Dongling should be:

However, in the actual operation, we encountered the problem that VTKMetaimagereader could not read the file, which was normal in the previous 32bit system, and we are not sure where the problem appears temporarily, which needs further study!! Here are the questions:


2. See the data

1. C++ Primer

2. The VTK User’s Guide — 11 The Edition

3. The Visualization Toolkit — AnObject-Oriented Approach To 3D Graphics (4th Edition)

4. Zhang Xiaodong, Luo Huoling. Advance of VTK Graphics Image Development [M]. Machinery Industry Press, 2015.


Read More: