医疗图像显示之切片提取
引言
切片是指三维图像中的一个切面对应的图像。切面可以是过图像内部一点且平行于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;
}
输出
横断面,冠状面,矢状面,斜切面