VTK学习笔记（十六）vtkImageReslice与vtkExtractVOI
1、关键代码：
先通过vtkMetaImageReader读取一副三维图像，获取图像范围、原点和像素间隔，由这三个参数可以计算图像的中心位置。
接下来定义了切面的变换矩阵axialElements，该矩阵的前三列分别表示X、Y和Z方向矢量，第四列为切面坐标系原点。通过修改切面坐标系原点，可以得到不同位置的切面图像。
然后将读取的图像作为vtkImageReslice的输入，通过函数SetResliceAxes()设置变换矩阵resliceAxis。
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();
2、完整代码：
#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 <vtkInteractorStyleImage.h>
#include <vtkCommand.h>
#include <vtkImageData.h>
#include <vtkMetaImageReader.h>
#include <vtkImageCast.h>
class vtkImageInteractionCallback : public vtkCommand
{
public:
static vtkImageInteractionCallback *New()
{
return new vtkImageInteractionCallback;
}
vtkImageInteractionCallback()
{
this->Slicing = 0;
this->ImageReslice = 0;
this->Interactor = 0;
}
void SetImageReslice(vtkImageReslice *reslice)
{
this->ImageReslice = reslice;
}
void SetImageMapToColors(vtkImageMapToColors *mapToColors)
{
this->mapToColors = mapToColors;
}
vtkImageReslice *GetImageReslice()
{
return this->ImageReslice;
}
void SetInteractor(vtkRenderWindowInteractor *interactor)
{
this->Interactor = interactor;
}
vtkRenderWindowInteractor *GetInteractor()
{
return this->Interactor;
}
virtual void Execute(vtkObject *, unsigned long event, void *)
{
vtkRenderWindowInteractor *interactor = this->GetInteractor();
int lastPos[2];
interactor->GetLastEventPosition(lastPos);
int currPos[2];
interactor->GetEventPosition(currPos);
if (event == vtkCommand::LeftButtonPressEvent)
{
this->Slicing = 1;
}
else if (event == vtkCommand::LeftButtonReleaseEvent)
{
this->Slicing = 0;
}
else if (event == vtkCommand::MouseMoveEvent)
{
if (this->Slicing)
{
vtkImageReslice *reslice = this->ImageReslice;
// Increment slice position by deltaY of mouse
int deltaY = lastPos[1] - currPos[1];
reslice->Update();
double sliceSpacing = reslice->GetOutput()->GetSpacing()[2];
vtkMatrix4x4 *matrix = reslice->GetResliceAxes();
// move the center point that we are slicing through
double point[4];
double center[4];
point[0] = 0.0;
point[1] = 0.0;
point[2] = sliceSpacing * 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]);
mapToColors->Update();
interactor->Render();
}
else
{
vtkInteractorStyle *style = vtkInteractorStyle::SafeDownCast(
interactor->GetInteractorStyle());
if (style)
{
style->OnMouseMove();
}
}
}
}
private:
int Slicing;
vtkImageReslice *ImageReslice;
vtkRenderWindowInteractor *Interactor;
vtkImageMapToColors *mapToColors;
};
int main()
{
vtkSmartPointer<vtkMetaImageReader> reader =
vtkSmartPointer<vtkMetaImageReader>::New();
reader->SetFileName ( "E:\\TestData\\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());
colorMap->Update();
vtkSmartPointer<vtkImageActor> imgActor =
vtkSmartPointer<vtkImageActor>::New();
imgActor->SetInputData(colorMap->GetOutput());
vtkSmartPointer<vtkRenderer> renderer =
vtkSmartPointer<vtkRenderer>::New();
renderer->AddActor(imgActor);
renderer->SetBackground(1, 1, 1);
vtkSmartPointer<vtkRenderWindow> renderWindow =
vtkSmartPointer<vtkRenderWindow>::New();
renderWindow->SetSize(500, 500);
renderWindow->AddRenderer(renderer);
vtkSmartPointer<vtkRenderWindowInteractor> renderWindowInteractor =
vtkSmartPointer<vtkRenderWindowInteractor>::New();
vtkSmartPointer<vtkInteractorStyleImage> imagestyle =
vtkSmartPointer<vtkInteractorStyleImage>::New();
renderWindowInteractor->SetInteractorStyle(imagestyle);
renderWindowInteractor->SetRenderWindow(renderWindow);
renderWindowInteractor->Initialize();
vtkSmartPointer<vtkImageInteractionCallback> callback =
vtkSmartPointer<vtkImageInteractionCallback>::New();
callback->SetImageReslice(reslice);
callback->SetInteractor(renderWindowInteractor);
callback->SetImageMapToColors(colorMap);
imagestyle->AddObserver(vtkCommand::MouseMoveEvent, callback);
imagestyle->AddObserver(vtkCommand::LeftButtonPressEvent, callback);
imagestyle->AddObserver(vtkCommand::LeftButtonReleaseEvent, callback);
renderWindowInteractor->Start();
return EXIT_SUCCESS;
}
3、CMakeLists.txt
cmake_minimum_required(VERSION 3.3 FATAL_ERROR)
project(testfilter)
set(VTK_DIR D:/ProgramFiles/vtk-9.0/lib/cmake/vtk-9.0)
find_package(VTK COMPONENTS
vtkCommonColor
vtkCommonCore
vtkCommonDataModel
vtkFiltersSources
vtkIOImage
vtkImagingStencil
vtkFiltersCore
vtkRenderingCore
vtkFiltersHybrid
vtkInteractionStyle
vtkRenderingFreeType
vtkRenderingOpenGL2
vtkRenderingVolumeOpenGL2
vtkImagingCore
vtkImagingColor
vtkImagingSources
vtkInteractionStyle
vtkRenderingContextOpenGL2
vtkRenderingGL2PSOpenGL2
vtkRenderingCore
QUIET
)
if (NOT VTK_FOUND)
message("Skipping testfilter: ${VTK_NOT_FOUND_MESSAGE}")
return()
endif()
message (STATUS "VTK_VERSION: ${VTK_VERSION}")
if (VTK_VERSION VERSION_LESS "8.90.0")
# old system
include(${VTK_USE_FILE})
endif()
file(GLOB_RECURSE mains ${CMAKE_CURRENT_SOURCE_DIR} *.cpp)
message(STATUS ${mains})
foreach(mainfile IN LISTS mains)
# Get file name without directory
get_filename_component(mainname ${mainfile} NAME_WE)
add_executable(${mainname} ${mainfile})
target_link_libraries (${mainname} ${VTK_LIBRARIES})
# vtk_module_autoinit is needed
if (NOT VTK_VERSION VERSION_LESS "8.90.0")
# vtk_module_autoinit is needed
vtk_module_autoinit(
TARGETS ${mainname}
MODULES ${VTK_LIBRARIES}
)
endif()
endforeach()
4、vtkImageReslice参数设置
设定的坐标点，vtkImageReslice生成的切面是一个无限大的二维图像，图像的中心点是设定的坐标点，在二维图像中是坐标原点(0,0)，根据需要可以设定二维图像左下角的坐标位置，也可以设定像素间距和宽高，然后从这个无限大的二维图像上截取一个矩形范围的二维图像作为输出。
vtkImageReslice获取的切面图像是一个有限范围的二维图像；可以想象从一个体数据中随意切一刀，会产生一个剖面，这个剖面的范围是不固定的，我们可以设置范围，也可以不设置范围。不设置范围，vtkImageReslice将使用默认值来框定二维图像的有效范围。设定范围后，vtkImageReslice就会从ResliceAxes变换矩阵的中心和向量信息计算出“视口”中的观察点位置（0,0）以及二维图像的左下角位置，一般就是相对于视口中心的偏移量，一般就是两个负值；
举例，如果要从当前ResliceAxes获取切面中获取一个256256的图像，XY方向上的像素间隔为1mm；如果是以当前视口的中心为二维图像的中心，则可以设置SetOutputOrigin的值为(spacing[0](extent[1]-extent[0])0.5-1, spacing[1]*(extent[3]-extent[2])0.5-1);
reslice->SetOutputExtent(0,255,0,255,0,0);
reslice->SetOutputSpacing(1,1,1);
reslice->SetOutputOrigin(-127.5,-127.5,0.0);
设置输出参数
SetOutputSpacing设置输出数据的体素间距。默认的输出间距是通过ResliceAxes排列的输入间距。
SetOutputOrigin设置输出数据的左下角开始的起点，默认的输出原点是通过ResliceAxes排列的输入原点。
SetOutputExtent设置输出数据的范围，默认的输出范围是通过ResliceAxes排列的输入范围。
需要注意的是：SetOutputOrigin设置的原点为XY相对位置，是像素*像素间距的值；
virtual void SetOutputSpacing(double x, double y, double z);
virtual void SetOutputSpacing(const double a[3]) { this->SetOutputSpacing(a[0], a[1], a[2]); }
vtkGetVector3Macro(OutputSpacing, double);
void SetOutputSpacingToDefault();
virtual void SetOutputOrigin(double x, double y, double z);
virtual void SetOutputOrigin(const double a[3]) { this->SetOutputOrigin(a[0], a[1], a[2]); }
vtkGetVector3Macro(OutputOrigin, double);
void SetOutputOriginToDefault();
virtual void SetOutputExtent(int a, int b, int c, int d, int e, int f);
virtual void SetOutputExtent(const int a[6])
{
this->SetOutputExtent(a[0], a[1], a[2], a[3], a[4], a[5]);
}
vtkGetVector6Macro(OutputExtent, int);
void SetOutputExtentToDefault();
vtkSetMacro(OutputDimensionality, int);
vtkGetMacro(OutputDimensionality, int);
参考：VTK笔记-计算MPR切面-vtkImageReslice输出视口设置
5、vtkExtractVOI
从一个大的体数据中抽取出一块小的体数据可以使用vtkExtractVOI实现。
常用接口：
vtkExtractVOI->SetInputData();
vtkExtractVOI->SetVOI(int,int,int,int,int,int);
vtkExtractVOI->SetSampleRate(int,int,int);
顾名思义，这三个接口分别是设置要抽取VOI的原始体数据，抽取VOI的范围，以及抽取数据后的采样率。
其中SetVOI()的6个参数值都是Int型的数字。而这个int型的数字是指抽取的范围在图像的Extent上的，而不是在图像的Bounds上的值。
当读入一组Dicom序列的图像时便形成了体数据，读取出该体数据的相关参数(origin, spacing, extent等)，在vtkBoxWidget的六个顶点bounds（XMin，XMax, YMin, YMax, ZMin, ZMax）形成的盒子区域，利用vtkExtractVOI抽取出这个盒子包围的区域。
vtkSmartPointer<vtkExtractVOI> extractVOI =
vtkSmartPointer<vtkExtractVOI>::New();
int extractXMin = (XMin-origin[0])/spacing[0];
int extractXMax = (XMax-origin[0])/spacing[0];
int extractYMin = (YMin-origin[1])/spacing[1];
int extractYMax = (YMax-origin[1])/spacing[1];
int extractZMin = (ZMin-origin[2])/spacing[2];
int extractZMax = (ZMax-origin[2])/spacing[2];
extractVOI ->SetInputData();
extractVOI->SetVOI(extractXMin ,extractXMax ,extractYMin ,extractYMax ,extractZMin ,extractZMax );
extractVOI->SetSampleRate(1,1,1);
extractVOI->Update();
参考：ExtractVOI
参考：使用vtkExtractVOI时未抽取出VOI的问题解决
6、vtkImageReslice实现crop
#include <vtkExtractVOI.h>
#include <vtkImageActor.h>
#include <vtkImageCast.h>
#include <vtkImageData.h>
#include <vtkImageMandelbrotSource.h>
#include <vtkImageMapper3D.h>
#include <vtkInteractorStyleImage.h>
#include <vtkNamedColors.h>
#include <vtkNew.h>
#include <vtkRenderWindow.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkRenderer.h>
#include <vtkNew.h>
#include <vtkNIFTIImageReader.h>
#include <vtkNIFTIImageWriter.h>
#include <vtkMatrix3x3.h>
#include <vtkImageReslice.h>
#include <string>
using namespace std;
int main(int, char*[])
{
string fileName = "D:/Dataset/Hip_Femur_Detection/Detect/BJ_02006.nii.gz";
vtkSmartPointer<vtkImageData> input;
vtkSmartPointer <vtkNIFTIImageReader> niftireader = vtkSmartPointer<vtkNIFTIImageReader>::New();
niftireader->SetFileName(fileName.c_str());
niftireader->Update();
input = niftireader->GetOutput();
int inputDims[3] = { 0 };
int extent[6] = {0};
double origin[3] = { 0 };
double spacing[3] = { 0 };
vtkNew<vtkMatrix3x3> mat;
input->GetDimensions(inputDims);
input->GetExtent(extent);
input->GetOrigin(origin);
input->GetSpacing(spacing);
auto matrix = input->GetDirectionMatrix();
string infileName = "input.nii.gz";
vtkSmartPointer<vtkNIFTIImageWriter> envImageWriter1 = vtkSmartPointer<vtkNIFTIImageWriter>::New();
envImageWriter1->SetInputData(input);
envImageWriter1->SetFileName(infileName.c_str());
envImageWriter1->Write();
std::cout << "Dims: "
<< " x: " << inputDims[0] << " y: " << inputDims[1]
<< " z: " << inputDims[2] << std::endl;
std::cout << "Number of points: " << input->GetNumberOfPoints()
<< std::endl;
std::cout << "Number of cells: " << input->GetNumberOfCells()
<< std::endl;
int* ext = input->GetExtent();
vtkSmartPointer<vtkImageReslice> pReslice = vtkSmartPointer<vtkImageReslice>::New();
pReslice->SetOutputOrigin(origin);
ext[5] = 0.8 * ext[5];
pReslice->SetOutputExtent(ext);
pReslice->SetAutoCropOutput(1);//！！！一定要加这句话
pReslice->SetInputData(input);
pReslice->Update();
vtkImageData* extracted = pReslice->GetOutput();
extracted->SetOrigin(origin);
int* extractedDims = extracted->GetDimensions();
std::cout << "Dims: "
<< " x: " << extractedDims[0] << " y: " << extractedDims[1]
<< " z: " << extractedDims[2] << std::endl;
std::cout << "Number of points: " << extracted->GetNumberOfPoints()
<< std::endl;
std::cout << "Number of cells: " << extracted->GetNumberOfCells()
<< std::endl;
string dstfileName = "resliced.nii.gz";
vtkSmartPointer<vtkNIFTIImageWriter> envImageWriter = vtkSmartPointer<vtkNIFTIImageWriter>::New();
envImageWriter->SetInputData(extracted);
envImageWriter->SetFileName(dstfileName.c_str());
envImageWriter->Write();
return EXIT_SUCCESS;
}
效果：
7、vtkExtractVOI实现crop
#include <vtkExtractVOI.h>
#include <vtkImageActor.h>
#include <vtkImageCast.h>
#include <vtkImageData.h>
#include <vtkImageMandelbrotSource.h>
#include <vtkImageMapper3D.h>
#include <vtkInteractorStyleImage.h>
#include <vtkNamedColors.h>
#include <vtkNew.h>
#include <vtkRenderWindow.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkRenderer.h>
#include <vtkNew.h>
#include <vtkNIFTIImageReader.h>
#include <vtkNIFTIImageWriter.h>
#include <vtkMatrix3x3.h>
#include <string>
using namespace std;
int main(int, char*[])
{
string fileName = "D:/Dataset/Hip_Femur_Detection/Detect/BJ_02006.nii.gz";
vtkSmartPointer<vtkImageData> input;
vtkSmartPointer <vtkNIFTIImageReader> niftireader = vtkSmartPointer<vtkNIFTIImageReader>::New();
niftireader->SetFileName(fileName.c_str());
niftireader->Update();
input = niftireader->GetOutput();
int inputDims[3] = { 0 };
int extent[6] = {0};
double origin[3] = { 0 };
double spacing[3] = { 0 };
vtkNew<vtkMatrix3x3> mat;
input->GetDimensions(inputDims);
input->GetExtent(extent);
input->GetOrigin(origin);
input->GetSpacing(spacing);
auto matrix = input->GetDirectionMatrix();
string infileName = "input.nii.gz";
vtkSmartPointer<vtkNIFTIImageWriter> envImageWriter1 = vtkSmartPointer<vtkNIFTIImageWriter>::New();
envImageWriter1->SetInputData(input);
envImageWriter1->SetFileName(infileName.c_str());
envImageWriter1->Write();
std::cout << "Dims: "
<< " x: " << inputDims[0] << " y: " << inputDims[1]
<< " z: " << inputDims[2] << std::endl;
std::cout << "Number of points: " << input->GetNumberOfPoints()
<< std::endl;
std::cout << "Number of cells: " << input->GetNumberOfCells()
<< std::endl;
vtkNew<vtkExtractVOI> extractVOI;
extractVOI->SetInputData(input);
extractVOI->SetVOI(0, inputDims[0]-1,
0, inputDims[1]-1, 0, inputDims[2] / 2.);
extractVOI->Update();
vtkImageData* extracted = extractVOI->GetOutput();
extracted->SetOrigin(origin);
int* extractedDims = extracted->GetDimensions();
std::cout << "Dims: "
<< " x: " << extractedDims[0] << " y: " << extractedDims[1]
<< " z: " << extractedDims[2] << std::endl;
std::cout << "Number of points: " << extracted->GetNumberOfPoints()
<< std::endl;
std::cout << "Number of cells: " << extracted->GetNumberOfCells()
<< std::endl;
string dstfileName = "extractVOI.nii.gz";
vtkSmartPointer<vtkNIFTIImageWriter> envImageWriter = vtkSmartPointer<vtkNIFTIImageWriter>::New();
envImageWriter->SetInputData(extracted);
envImageWriter->SetFileName(dstfileName.c_str());
envImageWriter->Write();
return EXIT_SUCCESS;
}
crop效果如下：
扫一扫
1204
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
