ITK 樣條曲線 BSpline 非剛性影像對準教學與範例

介紹如何使用 ITK 的樣條曲線（BSpline）轉換對 2D 影像進行非剛性影像對準。

這個範例中使用 BSplineTransform 轉換對 2D 的影像進行非剛性影像對準（non-rigid image registration），由於樣條曲線（BSpline）轉換的參數數量很龐大，所以這個例子中採用 LBFGSOptimizerv4 而不用普通的 RegularStepGradientDescentOptimizer 或 ConjugateGradientLineSearchOptimizer。

由於 BSplineTransform 轉換是由樣條曲線網格所組成的複雜變形，龐大的參數量可容許非常多種變形方式，但也會需要比較大量的計算。

#include <itkImageRegistrationMethodv4.h>
#include <itkCorrelationImageToImageMetricv4.h>
#include <itkTimeProbesCollectorBase.h>
#include <itkMemoryProbesCollectorBase.h>
#include <itkBSplineTransform.h>
#include <itkLBFGSOptimizerv4.h>
#include <itkImageFileReader.h>
#include <itkImageFileWriter.h>
#include <itkResampleImageFilter.h>
#include <itkCastImageFilter.h>
#include <itkSquaredDifferenceImageFilter.h>
#include <itkBSplineTransformInitializer.h>
#include <itkTransformToDisplacementFieldFilter.h>
#include <itkCheckerBoardImageFilter.h>

// LBFGSOptimizerv4 不會觸發事件，因此不用建立 Observer

int main(int argc, char * argv[]) {
  if (argc < 3) {
    std::cerr << "Missing Parameters " << std::endl;
    std::cerr << "Usage: " << argv[0];
    std::cerr << " fixedImageFile  movingImageFile" << std::endl;
    return EXIT_FAILURE;
  }

  // 影像維度
  constexpr unsigned int ImageDimension = 2;
  using PixelType = float;

  // 影像維度
  using FixedImageType = itk::Image<PixelType, ImageDimension>;
  using MovingImageType = itk::Image<PixelType, ImageDimension>;

  // 轉換空間維度
  const unsigned int SpaceDimension = ImageDimension;

  // Spline 曲線階數
  constexpr unsigned int SplineOrder = 3;

  // 座標資料類型
  using CoordinateRepType = double;

  // 建立 BSpline 轉換
  using TransformType = itk::BSplineTransform<CoordinateRepType, SpaceDimension, SplineOrder>;
  TransformType::Pointer transform = TransformType::New();

  // 建立 Optimizer
  using OptimizerType = itk::LBFGSOptimizerv4;
  OptimizerType::Pointer optimizer = OptimizerType::New();

  // 建立 Metric
  using MetricType = itk::CorrelationImageToImageMetricv4<FixedImageType, MovingImageType>;
  MetricType::Pointer metric = MetricType::New();

  // 建立 Registration
  using RegistrationType = itk::ImageRegistrationMethodv4<FixedImageType, MovingImageType>;
  RegistrationType::Pointer registration = RegistrationType::New();

  // 設定 Registration 的 Metric 與 Optimizer
  registration->SetMetric(metric);
  registration->SetOptimizer(optimizer);

  // 讀取影像
  using FixedImageReaderType = itk::ImageFileReader<FixedImageType>;
  using MovingImageReaderType = itk::ImageFileReader<MovingImageType>;

  FixedImageReaderType::Pointer fixedImageReader = FixedImageReaderType::New();
  MovingImageReaderType::Pointer movingImageReader = MovingImageReaderType::New();

  fixedImageReader->SetFileName(argv[1]);
  movingImageReader->SetFileName(argv[2]);

  fixedImageReader->Update();
  FixedImageType::ConstPointer fixedImage = fixedImageReader->GetOutput();

  // 以 BSplineTransformInitializer 設定 BSpline 轉換用的固定參數
  using InitializerType = itk::BSplineTransformInitializer<TransformType, FixedImageType>;
  InitializerType::Pointer transformInitializer = InitializerType::New();

  unsigned int numberOfGridNodesInOneDimension = 8;
  TransformType::MeshSizeType meshSize;
  meshSize.Fill(numberOfGridNodesInOneDimension - SplineOrder);

  transformInitializer->SetTransform(transform);
  transformInitializer->SetImage(fixedImage);
  transformInitializer->SetTransformDomainMeshSize(meshSize);
  transformInitializer->InitializeTransform();

  // 初始化轉換參數
  transform->SetIdentity();

  std::cout << "Initial Parameters = " << std::endl;
  std::cout << transform->GetParameters() << std::endl;

  // 設定 Registration 的轉換
  registration->SetInitialTransform(transform);

  // 將輸入的轉換物件直接用於輸出的轉換物件
  registration->InPlaceOn();

  // 設定輸入影像
  registration->SetFixedImage(fixedImage);
  registration->SetMovingImage(movingImageReader->GetOutput());

  // 自動估計參數尺度
  using ScalesEstimatorType = itk::RegistrationParameterScalesFromPhysicalShift<MetricType>;
  ScalesEstimatorType::Pointer scalesEstimator = ScalesEstimatorType::New();
  scalesEstimator->SetMetric(metric);
  scalesEstimator->SetTransformForward(true);
  scalesEstimator->SetSmallParameterVariation(1.0);

  optimizer->SetScalesEstimator(scalesEstimator);

  // 設定 Optimizer 參數
  optimizer->SetGradientConvergenceTolerance(5e-2);
  optimizer->SetLineSearchAccuracy(1.2);
  optimizer->SetDefaultStepLength(1.5);
  optimizer->TraceOn();
  optimizer->SetMaximumNumberOfFunctionEvaluations(1000);

  // 使用單一層級 Registration
  constexpr unsigned int numberOfLevels = 1;

  RegistrationType::ShrinkFactorsArrayType shrinkFactorsPerLevel;
  shrinkFactorsPerLevel.SetSize(1);
  shrinkFactorsPerLevel[0] = 1;

  RegistrationType::SmoothingSigmasArrayType smoothingSigmasPerLevel;
  smoothingSigmasPerLevel.SetSize(1);
  smoothingSigmasPerLevel[0] = 0;

  registration->SetNumberOfLevels(numberOfLevels);
  registration->SetSmoothingSigmasPerLevel(smoothingSigmasPerLevel);
  registration->SetShrinkFactorsPerLevel(shrinkFactorsPerLevel);


  // 測量時間與記憶體用量的探針
  itk::TimeProbesCollectorBase chronometer;
  itk::MemoryProbesCollectorBase memorymeter;

  std::cout << std::endl << "Starting Registration" << std::endl;

  try {
    memorymeter.Start("Registration");
    chronometer.Start("Registration");

    registration->Update();

    chronometer.Stop("Registration");
    memorymeter.Stop("Registration");

    std::cout << "Optimizer stop condition = "
              << registration->GetOptimizer()->GetStopConditionDescription()
              << std::endl;
  } catch (itk::ExceptionObject & err) {
    std::cerr << "ExceptionObject caught !" << std::endl;
    std::cerr << err << std::endl;
    return EXIT_FAILURE;
  }

  // 輸出時間與記憶體用量
  chronometer.Report(std::cout);
  memorymeter.Report(std::cout);

  // 取得對準結果轉換參數
  OptimizerType::ParametersType finalParameters = transform->GetParameters();

  std::cout << "Last Transform Parameters" << std::endl;
  std::cout << finalParameters << std::endl;

  // 建立 ResampleFilter 套用影像轉換
  using ResampleFilterType = itk::ResampleImageFilter<MovingImageType, FixedImageType>;
  ResampleFilterType::Pointer resample = ResampleFilterType::New();

  // 設定轉換
  resample->SetTransform(transform);

  // 設定輸入影像
  resample->SetInput(movingImageReader->GetOutput());

  // 設定輸出影像參數
  resample->SetSize(fixedImage->GetLargestPossibleRegion().GetSize());
  resample->SetOutputOrigin(fixedImage->GetOrigin());
  resample->SetOutputSpacing(fixedImage->GetSpacing());
  resample->SetOutputDirection(fixedImage->GetDirection());
  resample->SetDefaultPixelValue(100);

  // 輸出影像類型
  using OutputPixelType = unsigned char;
  using OutputImageType = itk::Image<OutputPixelType, ImageDimension>;

  // 影像轉型 Filter
  using CastFilterType = itk::CastImageFilter<FixedImageType, OutputImageType>;
  CastFilterType::Pointer caster = CastFilterType::New();
  caster->SetInput(resample->GetOutput());

  // 寫入影像
  using WriterType = itk::ImageFileWriter<OutputImageType>;
  WriterType::Pointer writer = WriterType::New();
  writer->SetFileName("output.png");
  writer->SetInput(caster->GetOutput());
  writer->Update();

  // 輸出影像對準的 Deformation Field
  using VectorPixelType = itk::Vector<float, ImageDimension>;
  using DisplacementFieldImageType = itk::Image<VectorPixelType, ImageDimension>;
  using DisplacementFieldGeneratorType =
    itk::TransformToDisplacementFieldFilter<DisplacementFieldImageType,
                                            CoordinateRepType>;

  // 建立並設定 Displacement Field Generator
  DisplacementFieldGeneratorType::Pointer dispfieldGenerator =
    DisplacementFieldGeneratorType::New();
  dispfieldGenerator->UseReferenceImageOn();
  dispfieldGenerator->SetReferenceImage(fixedImage);
  dispfieldGenerator->SetTransform(transform);
  dispfieldGenerator->Update();

  // 將 Displacement Field 寫入檔案
  using FieldWriterType = itk::ImageFileWriter<DisplacementFieldImageType>;
  FieldWriterType::Pointer fieldWriter = FieldWriterType::New();
  fieldWriter->SetInput(dispfieldGenerator->GetOutput());
  fieldWriter->SetFileName("displacement.mha");
  fieldWriter->Update();

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

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

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

  return EXIT_SUCCESS;
}

將這個 C++ 程式儲存為 BSplineReg.cxx 之後，搭配以下 CMakeLists.txt 設定檔以 CMake 編譯。

cmake_minimum_required(VERSION 3.10.2)

# 設定專案名稱
project(BSplineReg)

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

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

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

編譯與執行的指令如下：

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

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

這裡的 RatLungSlice1.mha 與 RatLungSlice2.mha 圖檔可以從 ITK 的 GitHub 網站取得，執行後除了可以產生對準後的圖形之外，也可以產生影像對準前後的棋盤式影像比較圖以及 displacement field。

這裡產生的 displacement.mha 可以使用以下 Python 程式繪製圖形：

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

displacement = sitk.ReadImage("displacement.mha")
movingImage = sitk.ReadImage("RatLungSlice2.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()

參考資料：DeformableRegistration4.cxx