C++/Python 三维图像切片

vtk9.2.6

以nii为例，其他的三位图像，用python可以参考https://blog.csdn.net/qq_39942341/article/details/122442949

C++

vtk9.2.6

main.cxx

#include<vtkNew.h>
#include<vtkNIFTIImageReader.h>
#include<vtkImageReader2.h>
#include<vtkMatrix4x4.h>
#include<vtkImageReslice.h>
#include<vtkImageData.h>
#include<vtkLookupTable.h>
#include<vtkImageMapToColors.h>
#include<vtkImageActor.h>
#include<vtkRenderer.h>
#include<vtkRenderWindow.h>
#include<vtkRenderWindowInteractor.h>
#include<vtkInteractorStyleImage.h>

void resliceByOrient(vtkImageReslice* reslice, vtkImageReader2* vtk_data, double center[3], int ori){
    static double axialElements[16]={
        1, 0, 0, 0,
        0, 1, 0, 0,
        0, 0, 1, 0,
        0, 0, 0, 1
    };
    static double sagittalElements[16]={
        0, 0, -1, 0,
        1, 0, 0, 0,
        0, -1, 0, 0,
        0, 0, 0, 1
    };
    static double coronalElements[16]={
        1, 0, 0, 0,
        0, 0, 1, 0,
        0, -1, 0, 0,
        0, 0, 0, 1
    };
    static double obliqueElements[16]={
        1, 0, 0, 0,
        0, 8.66025, -0.5, 0,
        0, 0.5, 0.866025, 0,
        0, 0, 0, 1
    };

    vtkNew<vtkMatrix4x4> resliceAxes;

    if(0 == ori){ // 横断面（transverse plane）A(z)
        resliceAxes->DeepCopy(axialElements);
    }
    else if(1 == ori){ // 矢状面（sagittal plane）S(x)
        resliceAxes->DeepCopy(sagittalElements);
    }
    else if(2 == ori){
        resliceAxes->DeepCopy(coronalElements);
    }
    else{
        resliceAxes->DeepCopy(obliqueElements);
    }
    for(int i = 0; i < 3; ++i){
        resliceAxes->SetElement(i, 3, center[i]);
    }

    //vtkNew<vtkImageReslice> reslice;
    reslice->SetInputConnection(vtk_data->GetOutputPort());
    reslice->SetOutputDimensionality(2);
    reslice->SetResliceAxes(resliceAxes);
    reslice->SetInterpolationModeToLinear();
    reslice->Update();
}



int main(){
    vtkNew<vtkNIFTIImageReader> reader;
    reader->SetFileName("/mnt/data/RibFrac101-image.nii.gz");
    reader->Update();

    double* origin = reader->GetDataOrigin();
    double* spacing = reader->GetDataSpacing();
    
    //257,414,177表示对应面的切片
    double center[3] = {
        origin[0] + spacing[0] * 257, //矢状面（sagittal plane）S
        origin[1] + spacing[1] * 414, //冠状面（coronal plane）C
        origin[2] + spacing[2] * 177 //横断面（transverse plane）A
    };

    vtkNew<vtkImageReslice> reslice;
    resliceByOrient(reslice, reader, center, 1);

    double* scale_range = reader->GetOutput()->GetScalarRange();
    vtkNew<vtkLookupTable> colorTable;
    colorTable->SetRange(scale_range[0], scale_range[1]);
    colorTable->SetValueRange(0.0, 1.0);
    colorTable->SetSaturationRange(0.0, 0.0);
    colorTable->SetRampToLinear();
    colorTable->Build();

    vtkNew<vtkImageMapToColors> colorMap;
    colorMap->SetLookupTable(colorTable);
    colorMap->SetInputConnection(reslice->GetOutputPort());
    colorMap->Update();

    vtkNew<vtkImageActor> imgActor;
    imgActor->SetInputData(colorMap->GetOutput());

    vtkNew<vtkRenderer> renderer;
    renderer->AddActor(imgActor);
    renderer->SetBackground(0.0, 0.0, 0.0);

    vtkNew<vtkRenderWindow> renderWindow;
    renderWindow->AddRenderer(renderer);
    renderWindow->SetSize(640, 400);

    vtkNew<vtkRenderWindowInteractor> rwi;
    vtkNew<vtkInteractorStyleImage> imageStyle;
    rwi->SetInteractorStyle(imageStyle);
    rwi->SetRenderWindow(renderWindow);
    rwi->Initialize();
    renderWindow->Render();
    rwi->Start();


    return 0;
}

CMakeLists.txt

cmake_minimum_required(VERSION 3.12 FATAL_ERROR)

project(main)

find_package(VTK COMPONENTS 
  CommonColor
  CommonCore
  CommonDataModel
  CommonMath
  CommonTransforms
  FiltersCore
  FiltersGeometry
  IOImage
  IOXML
  ImagingCore
  InteractionStyle
  RenderingContextOpenGL2
  RenderingCore
  RenderingFreeType
  RenderingGL2PSOpenGL2
  RenderingOpenGL2
)

if (NOT VTK_FOUND)
  message(FATAL_ERROR "main: Unable to find the VTK build folder.")
endif()

# Prevent a "command line is too long" failure in Windows.
set(CMAKE_NINJA_FORCE_RESPONSE_FILE "ON" CACHE BOOL "Force Ninja to use response files.")
add_executable(main MACOSX_BUNDLE main.cxx )
  target_link_libraries(main PRIVATE ${VTK_LIBRARIES}
)
# vtk_module_autoinit is needed
vtk_module_autoinit(
  TARGETS main
  MODULES ${VTK_LIBRARIES}
)

python

python3.8.10
vtk9.2.6

#!/usr/bin/env python
# _*_ coding:utf-8 _*_

# noinspection PyUnresolvedReferences
import vtkmodules.vtkInteractionStyle
# noinspection PyUnresolvedReferences
import vtkmodules.vtkRenderingOpenGL2
from vtkmodules.vtkCommonCore import (
    vtkLookupTable
)
from vtkmodules.vtkCommonMath import (
    vtkMatrix4x4
)
from vtkmodules.vtkIOImage import (
    vtkNIFTIImageReader
)
from vtkmodules.vtkImagingCore import (
    vtkImageReslice,
    vtkImageMapToColors
)
from vtkmodules.vtkInteractionStyle import (
    vtkInteractorStyleImage
)
from vtkmodules.vtkRenderingCore import (
    vtkImageActor,
    vtkRenderer,
    vtkRenderWindow,
    vtkRenderWindowInteractor
)


def resliceByOrient(vtk_data, center, ori):
    """
    重采样

    :param vtk_data: vtk的数据
    :param center: 采样中心(S,C,A)(x,y,z)
    :param ori: 1:S(x),2:C(y),0:A(z)
    :return: 重采样的数据
    """
    if 1 == ori:
        # 矢状面（sagittal plane）S(x)
        elements = [
            0, 0, -1, 0,
            1, 0, 0, 0,
            0, -1, 0, 0,
            0, 0, 0, 1
        ]
    elif 2 == ori:
        # 冠状面（coronal plane）C(y)
        elements = [
            1, 0, 0, 0,
            0, 0, 1, 0,
            0, -1, 0, 0,
            0, 0, 0, 1
        ]
    elif 0 == ori:
        # 横断面（transverse plane）A(z)
        elements = [
            1, 0, 0, 0,
            0, 1, 0, 0,
            0, 0, 1, 0,
            0, 0, 0, 1
        ]
    else:
        # 斜切切片
        elements = [
            1, 0, 0, 0,
            0, 0, 8.66025, -0.5, 0,
            0, 0.5, 0.866025, 0,
            0, 0, 0, 1
        ]
    resliceAxes = vtkMatrix4x4()
    resliceAxes.DeepCopy(elements)
    for i in range(3):
        resliceAxes.SetElement(i, 3, center[i])
        # resliceAxes.SetElement(i, 3, int(center[i]))
    reslice = vtkImageReslice()
    reslice.SetInputConnection(vtk_data.GetOutputPort())
    reslice.SetOutputDimensionality(2)
    reslice.SetResliceAxes(resliceAxes)
    reslice.SetInterpolationModeToLinear()
    reslice.Update()
    return reslice


def get_color_table(is_image, minimum=None, maximum=None, window=None, level=None):
    """
    获取颜色映射表

    :param is_image: 是否为图像
    :param minimum: 最小值
    :param maximum: 最大值
    :param window: 窗宽
    :param level: 窗位
    :return: 颜色映射表
    """
    colorTable = vtkLookupTable()
    if is_image:
        if minimum is None and maximum is None and window is not None and level is not None:
            minimum = (2 * level - window) / 2.0 + 0.5
            maximum = (2 * level + window) / 2.0 + 0.5
        elif minimum is None or maximum is None:
            raise Exception("can't get right min max")

        # 像素最小值到最大值，（可以根据窗宽窗位计算出像素范围）
        colorTable.SetRange(minimum, maximum)
        colorTable.SetValueRange(0.0, 1.0)
        colorTable.SetSaturationRange(0.0, 0.0)
        colorTable.SetRampToLinear()
    else:
        # TODO 改成正常的颜色
        colorTable.SetNumberOfTableValues(5)
        colorTable.SetRange(0.0, 4.0)
        colorTable.SetTableValue(0, 0, 0.0, 0.0, 0.0)
        colorTable.SetTableValue(1, 0.0, 1.0, 0.0, 0.4)
        colorTable.SetTableValue(2, 1.0, 0.0, 0.0, 0.4)
        colorTable.SetTableValue(3, 0.0, 0.0, 1.0, 0.4)
        colorTable.SetTableValue(4, 0.0, 0.0, 1.0, 0.4)
    colorTable.Build()
    return colorTable


if __name__ == '__main__':
    reader = vtkNIFTIImageReader()
    reader.SetFileName('RibFrac101-image.nii.gz')
    reader.Update()

    origin = reader.GetDataOrigin()
    spacing = reader.GetDataSpacing()
    direction = reader.GetDataDirection()
    extent = reader.GetDataExtent()
    # 257,414,177表示对应面的切片
    center = [
        origin[0] + spacing[0] * 257,  # 矢状面（sagittal plane）S
        origin[1] + spacing[1] * 414,  # 冠状面（coronal plane）C
        origin[2] + spacing[2] * 177,  # 横断面（transverse plane）A
    ]
    print(center)
    reslice = resliceByOrient(reader, center, 1)
    colorTable = get_color_table(True, *(reader.GetOutput().GetScalarRange()))

    colorMap = vtkImageMapToColors()
    colorMap.SetLookupTable(colorTable)
    colorMap.SetInputConnection(reslice.GetOutputPort())
    colorMap.Update()

    imgActor = vtkImageActor()
    imgActor.SetInputData(colorMap.GetOutput())

    renderer = vtkRenderer()
    renderer.AddActor(imgActor)
    renderer.SetBackground(0.0, 0.0, 0.0)

    renderWindow = vtkRenderWindow()
    renderWindow.AddRenderer(renderer)
    renderWindow.SetSize(640, 400)

    rwi = vtkRenderWindowInteractor()
    imageStyle = vtkInteractorStyleImage()
    rwi.SetInteractorStyle(imageStyle)
    rwi.SetRenderWindow(renderWindow)
    rwi.Initialize()
    renderWindow.Render()
    rwi.Start()

《VTK图形图像开发进阶》P97 三维图像切面提取