ITK 使用 TriangleMeshToBinaryImageFilter 將 Surface 轉為二元影像教學與範例

介紹如何使用 ITK 的 TriangleMeshToBinaryImageFilter 將網格（mesh）包圍的範圍轉為二元影像（binary image）。

Python 語言版本

在使用 ITK 的 TriangleMeshToBinaryImageFilter 將網格轉為二元影像的時候，需要指定網格與影像的類型，可用的類型可以用 GetTypes 函數查詢：

import SimpleITK as sitk
import numpy as np
import itk
import matplotlib.pyplot as plt

# 查詢 TriangleMeshToBinaryImageFilter 支援的類型
itk.TriangleMeshToBinaryImageFilter.GetTypes()

<itkTemplate itk::TriangleMeshToBinaryImageFilter>
Options:
  [<class 'itk.itkMeshBasePython.itkMeshD3'>, <class 'itk.itkImagePython.itkImageD3'>]
  [<class 'itk.itkMeshBasePython.itkMeshF3'>, <class 'itk.itkImagePython.itkImageF3'>]
  [<class 'itk.itkMeshBasePython.itkMeshSS3'>, <class 'itk.itkImagePython.itkImageSS3'>]
  [<class 'itk.itkMeshBasePython.itkMeshUC3'>, <class 'itk.itkImagePython.itkImageUC3'>]
  [<class 'itk.itkMeshBasePython.itkMeshUS3'>, <class 'itk.itkImagePython.itkImageUS3'>]

這裡我們設定網格與影像的類型皆為無號短整數：

# 設定網格與影像的類型
meshType = itk.Mesh[itk.US, 3]
imageType = itk.Image[itk.US, 3]

接著讀取網格與影像檔案：

# 讀取 Mesh 網格檔案
meshReader = itk.MeshFileReader[meshType].New()
meshIO = itk.OBJMeshIO.New()
meshReader.SetMeshIO(meshIO)
meshReader.SetFileName("annotation/ccf_2017/structure_meshes/477.obj")
meshReader.Update()
mesh = meshReader.GetOutput()

# 讀取影像檔案
image = itk.imread("average_template/average_template_25.nrrd")

若要查看網格的資料，可以使用 itkwidgets 來顯示：

# 顯示網格資料
import itkwidgets
itkwidgets.view(geometries = mesh)

透過 TriangleMeshToBinaryImageFilter 將網格轉為二元影像（遮罩影像）：

# 建立 TriangleMeshToBinaryImageFilter
mesh2binFilter = itk.TriangleMeshToBinaryImageFilter[meshType, imageType].New()

# 設定輸入網格
mesh2binFilter.SetInput(mesh)

# 根據影像設定影像資訊（spacing、size、origin）
mesh2binFilter.SetInfoImage(image)

# 亦可自行設定影像資訊
#mesh2binFilter.SetSpacing(spacing)
#mesh2binFilter.SetSize(size)
#mesh2binFilter.SetOrigin(origin)

# 設定值區域內外值
mesh2binFilter.SetInsideValue(255)
mesh2binFilter.SetOutsideValue(0)

# 實際執行
mesh2binFilter.Update()
mask = mesh2binFilter.GetOutput()

這裡產生的 mask 就是遮罩影像，網格內部的值會是 255，外部則會是 0。

# 顯示影像切面
nda = itk.GetArrayViewFromImage(mask)
fig, axs = plt.subplots(1, 3, figsize=(15, 5))
axs[0].imshow(nda[nda.shape[0]//2,:,:], cmap = 'gray')
axs[1].imshow(nda[:,nda.shape[1]//2,:], cmap = 'gray')
axs[2].imshow(nda[:,:,nda.shape[2]//2], cmap = 'gray')
plt.show()

C++ 語言版本

這是要 C++ 語言中使用 TriangleMeshToBinaryImageFilter 的範例，基本概念都相同，只是以不同程式語言實作而已。

#include <itkMesh.h>
#include <itkMeshFileReader.h>
#include <itkImage.h>
#include <itkImageFileReader.h>
#include <itkImageFileWriter.h>
#include <itkCastImageFilter.h>
#include <itkTriangleMeshToBinaryImageFilter.h>

int main() {
  // 設定 Mesh 資料型態（浮點數、維度 3）
  constexpr unsigned int Dimension = 3;
  using MeshPixelType = double;
  using MeshType = itk::Mesh<MeshPixelType, Dimension>;

  // 建立 Mesh Reader
  using MeshReaderType = itk::MeshFileReader<MeshType>;
  MeshReaderType::Pointer meshReader = MeshReaderType::New();

  // 指定 Mesh 檔案名稱
  meshReader->SetFileName("input_mesh.obj");

  // 設定影像資料型態
  using PixelType = unsigned char;
  using ImageType = itk::Image<PixelType, Dimension>;

  // 建立 Image Reader
  using ImageReaderType = itk::ImageFileReader<ImageType>;
  ImageReaderType::Pointer imageReader = ImageReaderType::New();

  // 指定影像檔案名稱
  imageReader->SetFileName("input_image.nrrd");

  // 建立 TriangleMeshToBinaryImageFilter
  using FilterType = itk::TriangleMeshToBinaryImageFilter<MeshType, ImageType>;
  FilterType::Pointer filter = FilterType::New();

  // 指定輸入 Mesh
  filter->SetInput(meshReader->GetOutput());

  // 設定影像資訊
  filter->SetInfoImage(imageReader->GetOutput());

  // 設定 Mesh 內部數值
  filter->SetInsideValue(itk::NumericTraits<PixelType>::max());

  // 實際進行轉換
  try {
    filter->Update();
  } catch (itk::ExceptionObject & error) {
    std::cerr << "Error: " << error << std::endl;
    return EXIT_FAILURE;
  }

  // 建立 Image Writer，輸出轉換結果影像
  using WriterType = itk::ImageFileWriter<ImageType>;
  WriterType::Pointer writer = WriterType::New();
  writer->SetFileName("output_image.nrrd");
  writer->SetInput(filter->GetOutput());

  // 實際寫入檔案
  try {
    writer->Update();
  } catch (itk::ExceptionObject & error) {
    std::cerr << "Error: " << error << std::endl;
    return EXIT_FAILURE;
  }

  return EXIT_SUCCESS;
}

搭配以下 CMakeLists.txt 以 CMake 編譯：

cmake_minimum_required(VERSION 3.10.2)

# 設定專案名稱
project(ConvMesh)

# 尋找並引入 ITK 函式庫
find_package(ITK REQUIRED)
include(${ITK_USE_FILE})

# 增加一個執行檔
add_executable(ConvMesh ConvMesh.cpp)

# 定義執行檔連結方式
target_link_libraries(ConvMesh ${ITK_LIBRARIES})

參考資料：ITK Examples、ITKPythonPackage