VTK笔记——医学图像的切片提取(vtkImageReslice)

医学图像的浏览和内部分析是很常见也很重要的功能，我们不仅可以从矢状面、冠状面和轴状面这样的切面去看，还可以从任意切面去看。在VTK术语中，切面（切片）就是图像数据。

vtkImageReslice

vtkImageReslice，功能非常强大，有“瑞士军刀”的美誉。它不仅可以提取切片，还有图像旋转，翻转，重采样，变形等功能，并且效率还很高。这里只用到它的提取功能。
切片提取
切片提取的关键是要给出ResliceAxes，实际上，它是一个4X4的矩阵。我们可以这样来理解：把这个矩阵看作是新坐标系在原坐标下的坐标表示，一个齐次坐标系。前三列分别表示新坐标系下x、y、z的向量坐标，第四列表示新坐标系下原点的坐标。例如：
0, 0, -1, 6
1, 0, 0, 0
0, -1, 0, 0
0, 0, 0, 1
这个矩阵表示了在新坐标系下，x方向向量为（0, 1, 0），y方向向量为（0, 0, -1），z方向向量为（-1, 0, 0）和原点坐标为（6, 0, 0）矩阵。需要注意到是x, y, z三个向量要符合右手法则。
下面两种写法，作用一样：
写法一：

 auto resliceAxes = vtkSmartPointer<vtkMatrix4x4>::New();auto ImageReslice = vtkSmartPointer<vtkImageReslice>::New();ImageReslice->SetResliceAxes(resliceAxes);

写法二：

 double x[3] = { 1, 0, 0 };double y[3] = { 0, 1, 0 };double z[3] = { 0, 0, 1 };double origin[3] = { 1, 0, 0 };ImageReslice->SetResliceAxesDirectionCosines(x, y, z);ImageReslice->SetResliceAxesOrigin(origin);

脑部CT图像

ResliceImage.cxx

#include "vtkSmartPointer.h"
#include "vtkMatrix4x4.h"
#include "vtkDICOMImageReader.h"
#include "vtkImageData.h"
#include "vtkImageReslice.h"
#include "vtkLookupTable.h"
#include "vtkImageMapToColors.h"
#include "vtkImageActor.h"
#include "vtkImageMapper3D.h"
#include "vtkRenderer.h"
#include "vtkRenderWindow.h"
#include "vtkRenderWindowInteractor.h"
#include "vtkInteractorStyleImage.h"
#include "vtkCommand.h"class vtkImageInteractionCallback : public vtkCommand
{public:static vtkImageInteractionCallback* New(){return new vtkImageInteractionCallback;}vtkImageInteractionCallback(){this->Slicing = 0;this->ImageReslice = nullptr;this->Interactor = nullptr;}void SetImageReslice(vtkImageReslice* reslice){this->ImageReslice = reslice;}void SetMapToColors(vtkImageMapToColors* colors){this->MapToColors = colors;}void SetInteractor(vtkRenderWindowInteractor* interactor){this->Interactor = interactor;}void SetRenderWindow(vtkRenderWindow* window){this->RenderWindow = window;}void Execute(vtkObject* caller, unsigned long eventId, void* callData) override{int lastPos[2], curPos[2];this->Interactor->GetLastEventPosition(lastPos);this->Interactor->GetEventPosition(curPos);if (eventId == vtkCommand::LeftButtonPressEvent){this->Slicing = 1;}else if (eventId == vtkCommand::LeftButtonReleaseEvent){this->Slicing = 0;}else if (eventId == vtkCommand::MouseMoveEvent){if (this->Slicing){int deltaY = lastPos[1] - curPos[1];this->ImageReslice->Update();double spacing = this->ImageReslice->GetOutput()->GetSpacing()[2];vtkMatrix4x4* matrix = this->ImageReslice->GetResliceAxes();double point[4], center[4];point[0] = 0.0;point[1] = 0.0;point[2] = spacing * deltaY;point[3] = 1.0;matrix->MultiplyPoint(point, center);matrix->SetElement(0, 3, center[0]);matrix->SetElement(1, 3, center[1]);matrix->SetElement(2, 3, center[2]);this->MapToColors->Update();this->OutputImageData = this->MapToColors->GetOutput();this->Interactor->Render();}else{vtkInteractorStyle* style = vtkInteractorStyle::SafeDownCast(this->Interactor->GetInteractorStyle());if (style)style->OnMouseMove();}}}private:int Slicing;vtkImageReslice* ImageReslice;vtkImageMapToColors* MapToColors;vtkRenderWindowInteractor* Interactor;vtkRenderWindow* RenderWindow;vtkImageData* OutputImageData;
};int main(int argc, char* argv[])
{if (argc < 2){std::cout << "Usage: " << argv[0] << " DicomDirectory" << std::endl;return EXIT_FAILURE;}auto reader = vtkSmartPointer<vtkDICOMImageReader>::New();reader->SetDirectoryName(argv[1]);reader->Update();int extent[6];double spacing[3];double origin[3];double center[3];reader->GetOutput()->GetExtent(extent);reader->GetOutput()->GetSpacing(spacing);reader->GetOutput()->GetOrigin(origin);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 sagittalElements[16] = {0, 0,-1, 0,1, 0, 0, 0,0,-1, 0, 0,0, 0, 0, 1 };static double cornalElements[16] = {1, 0, 0, 0,0, 0, 1, 0,0, -1,0, 0,0, 0, 0, 1 };static double axialElements[16] = {1, 0, 0, 0,0, 1, 0, 0,0, 0, 1, 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 };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]);double x[3] = { 1, 0, 0 };double y[3] = { 0, 1, 0 };double z[3] = { 0, 0, 1 };auto ImageReslice = vtkSmartPointer<vtkImageReslice>::New();ImageReslice->SetInputConnection(reader->GetOutputPort());ImageReslice->SetOutputDimensionality(2);
#if 0ImageReslice->SetResliceAxes(resliceAxes);
#elseImageReslice->SetResliceAxesDirectionCosines(x, y, z);ImageReslice->SetResliceAxesOrigin(center);
#endifImageReslice->SetInterpolationModeToLinear();auto lookupTable = vtkSmartPointer<vtkLookupTable>::New();lookupTable->SetRange(0, 2000);lookupTable->SetValueRange(0.0, 1.0);lookupTable->SetSaturationRange(0.0, 0.0);lookupTable->SetRampToLinear();lookupTable->Build();auto mapToColors = vtkSmartPointer<vtkImageMapToColors>::New();mapToColors->SetLookupTable(lookupTable);mapToColors->SetInputConnection(ImageReslice->GetOutputPort());mapToColors->Update();auto imageActor = vtkSmartPointer<vtkImageActor>::New();imageActor->GetMapper()->SetInputConnection(mapToColors->GetOutputPort());auto renderer = vtkSmartPointer<vtkRenderer>::New();renderer->AddActor(imageActor);auto renderWindow = vtkSmartPointer<vtkRenderWindow>::New();renderWindow->AddRenderer(renderer);auto imageStyle = vtkSmartPointer<vtkInteractorStyleImage>::New();auto interactor = vtkSmartPointer<vtkRenderWindowInteractor>::New();interactor->SetInteractorStyle(imageStyle);interactor->SetRenderWindow(renderWindow);interactor->Initialize();auto callback = vtkSmartPointer<vtkImageInteractionCallback>::New();callback->SetImageReslice(ImageReslice);callback->SetMapToColors(mapToColors);callback->SetInteractor(interactor);callback->SetRenderWindow(renderWindow);imageStyle->AddObserver(vtkCommand::LeftButtonPressEvent, callback);imageStyle->AddObserver(vtkCommand::LeftButtonReleaseEvent, callback);imageStyle->AddObserver(vtkCommand::MouseMoveEvent, callback);interactor->Start();return EXIT_SUCCESS;
}

Example Download

Ref

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