ITK 有限元素法影像對準 FEM Registration 教學與範例

#include <itkImageFileReader.h>
#include <itkImageFileWriter.h>
#include <itkRescaleIntensityImageFilter.h>
#include <itkHistogramMatchingImageFilter.h>
#include <itkFEMRegistrationFilter.h>
#include <itkCheckerBoardImageFilter.h>

int main(int argc, char * argv[]) {
  const char *fixedImageName, *movingImageName;
  if (argc < 2) {
    std::cout << "Usage: " << argv[0] << " fixedImageFile movingImageFile"
              << std::endl;
    return EXIT_FAILURE;
  } else {
    fixedImageName = argv[1];
    movingImageName = argv[2];
  }

  // 定義適用於 2D FEM 影像對準的型別
  using DiskImageType = itk::Image<unsigned char, 2>;
  using ImageType = itk::Image<float, 2>;
  using ElementType = itk::fem::Element2DC0LinearQuadrilateralMembrane;
  using FEMObjectType = itk::fem::FEMObject<2>;
  using RegistrationType =
    itk::fem::FEMRegistrationFilter<ImageType, ImageType, FEMObjectType>;

  // 建立 FEMRegistrationFilter，並設定各種參數
  RegistrationType::Pointer registrationFilter = RegistrationType::New();
  registrationFilter->SetMaxLevel(1);
  registrationFilter->SetUseNormalizedGradient(true);
  registrationFilter->ChooseMetric(0);

  unsigned int maxiters = 20;
  float        E = 100;
  float        p = 1;
  registrationFilter->SetElasticity(E, 0);
  registrationFilter->SetRho(p, 0);
  registrationFilter->SetGamma(1., 0);
  registrationFilter->SetAlpha(1.);
  registrationFilter->SetMaximumIterations(maxiters, 0);
  registrationFilter->SetMeshPixelsPerElementAtEachResolution(4, 0);
  registrationFilter->SetWidthOfMetricRegion(1, 0);
  registrationFilter->SetNumberOfIntegrationPoints(2, 0);
  registrationFilter->SetDoLineSearchOnImageEnergy(0);
  registrationFilter->SetTimeStep(1.);
  registrationFilter->SetEmployRegridding(false);
  registrationFilter->SetUseLandmarks(false);

  // 讀取影像
  using FileSourceType = itk::ImageFileReader<DiskImageType>;

  FileSourceType::Pointer movingfilter = FileSourceType::New();
  movingfilter->SetFileName(movingImageName);
  FileSourceType::Pointer fixedfilter = FileSourceType::New();
  fixedfilter->SetFileName(fixedImageName);
  movingfilter->Update();
  fixedfilter->Update();

  // 將影像數值範圍調整為 0 至 255
  using FilterType =
    itk::RescaleIntensityImageFilter<DiskImageType, ImageType>;
  FilterType::Pointer movingrescalefilter = FilterType::New();
  FilterType::Pointer fixedrescalefilter = FilterType::New();

  movingrescalefilter->SetInput(movingfilter->GetOutput());
  fixedrescalefilter->SetInput(fixedfilter->GetOutput());

  constexpr double desiredMinimum = 0.0;
  constexpr double desiredMaximum = 255.0;

  movingrescalefilter->SetOutputMinimum(desiredMinimum);
  movingrescalefilter->SetOutputMaximum(desiredMaximum);
  movingrescalefilter->UpdateLargestPossibleRegion();
  fixedrescalefilter->SetOutputMinimum(desiredMinimum);
  fixedrescalefilter->SetOutputMaximum(desiredMaximum);
  fixedrescalefilter->UpdateLargestPossibleRegion();

  // 以直方圖分布的方式標準化影像值分布
  using HEFilterType =
    itk::HistogramMatchingImageFilter<ImageType, ImageType>;
  HEFilterType::Pointer IntensityEqualizeFilter = HEFilterType::New();

  IntensityEqualizeFilter->SetReferenceImage(fixedrescalefilter->GetOutput());
  IntensityEqualizeFilter->SetInput(movingrescalefilter->GetOutput());
  IntensityEqualizeFilter->SetNumberOfHistogramLevels(100);
  IntensityEqualizeFilter->SetNumberOfMatchPoints(15);
  IntensityEqualizeFilter->ThresholdAtMeanIntensityOn();
  IntensityEqualizeFilter->Update();

  // 設定 FEMRegistrationFilter 的固定影像與調動影像
  registrationFilter->SetFixedImage(fixedrescalefilter->GetOutput());
  registrationFilter->SetMovingImage(IntensityEqualizeFilter->GetOutput());

  //  In order to initialize the mesh of elements, we must first create
  //  ``dummy'' material and element objects and assign them to the
  //  registration filter.  These objects are subsequently used to
  //  either read a predefined mesh from a file or generate a mesh using
  //  the software.  The values assigned to the fields within the
  //  material object are arbitrary since they will be replaced with
  //  those specified earlier.  Similarly, the element
  //  object will be replaced with those from the desired mesh.

  // 建立材質性質
  itk::fem::MaterialLinearElasticity::Pointer m;
  m = itk::fem::MaterialLinearElasticity::New();
  m->SetGlobalNumber(0);
  m->SetYoungsModulus(registrationFilter->GetElasticity()); // 楊氏模數
  m->SetCrossSectionalArea(1.0);
  m->SetThickness(1.0);
  m->SetMomentOfInertia(1.0);
  m->SetPoissonsRatio(0.); // 蒲松比（不要設 1.0）
  m->SetDensityHeatProduct(1.0);

  // 建立 ElementType
  ElementType::Pointer e1 = ElementType::New();
  e1->SetMaterial(static_cast<itk::fem::Material::ConstPointer>(m));
  registrationFilter->SetElement(static_cast<itk::fem::Element::Pointer>(e1));
  registrationFilter->SetMaterial(m);

  // 進行影像對準
  registrationFilter->RunRegistration();

  // 儲存對準後的影像
  itk::ImageFileWriter<ImageType>::Pointer warpedImageWriter;
  warpedImageWriter = itk::ImageFileWriter<ImageType>::New();
  warpedImageWriter->SetInput(registrationFilter->GetWarpedImage());
  warpedImageWriter->SetFileName("warpedMovingImage.mha");
  warpedImageWriter->Update();

  // 儲存 displacement field
  using DispWriterType = itk::ImageFileWriter<RegistrationType::FieldType>;
  DispWriterType::Pointer dispWriter = DispWriterType::New();
  dispWriter->SetInput(registrationFilter->GetDisplacementField());
  dispWriter->SetFileName("displacement.mha");
  dispWriter->Update();

  // 將輸出影像轉型為指定類型
  using OutputImageType = itk::Image<unsigned char, 2>;
  using CastFilterType =
    itk::CastImageFilter<ImageType, OutputImageType>;
  using WriterType = itk::ImageFileWriter<OutputImageType>;

  WriterType::Pointer writer = WriterType::New();
  CastFilterType::Pointer caster = CastFilterType::New();

  // 產生棋盤式影像比較圖
  using CheckerBoardFilterType = itk::CheckerBoardImageFilter<ImageType>;
  CheckerBoardFilterType::Pointer checker = CheckerBoardFilterType::New();

  // 輸出影像對準前的棋盤式影像比較圖
  checker->SetInput1(registrationFilter->GetFixedImage());
  checker->SetInput2(registrationFilter->GetMovingImage());
  caster->SetInput(checker->GetOutput());
  writer->SetInput(caster->GetOutput());
  writer->SetFileName("before.png");
  writer->Update();

  // 輸出影像對準後的棋盤式影像比較圖
  checker->SetInput1(registrationFilter->GetFixedImage());
  checker->SetInput2(registrationFilter->GetWarpedImage());
  caster->SetInput(checker->GetOutput());
  writer->SetInput(caster->GetOutput());
  writer->SetFileName("after.png");
  writer->Update();

  return EXIT_SUCCESS;
}

cmake_minimum_required(VERSION 3.10.2)

# 設定專案名稱
project(FEMReg)

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

# 增加一個執行檔
add_executable(FEMReg FEMReg.cxx)

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

# 編譯程式
mkdir build
cd build
cmake ..

# 執行程式
./FEMReg RatLungSlice1.mha RatLungSlice2.mha

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

displacement = sitk.ReadImage("displacement.mha")
movingImage = sitk.ReadImage("moving.mha")

dispNDA = sitk.GetArrayViewFromImage(displacement)
movingNDA = sitk.GetArrayViewFromImage(movingImage)

ds = 4
X = np.arange(0, movingNDA.shape[0], ds)
Y = np.arange(0, movingNDA.shape[1], ds)
fig = plt.figure(figsize=(10, 10))
ax = fig.add_subplot()
ax.imshow(movingNDA, cmap='gray')
ax.quiver(X, Y, dispNDA[::ds,::ds,0], dispNDA[::ds,::ds,1], color='r', units='xy', scale=1)
plt.show()