ITK 多層級 Multi-Resolution 影像對準架構教學與範例

介紹如何在 ITK 的多層級（multi-resolution）影像對準架構下，進行線性仿射影像對準（affine registration）。

2D 多層級線性影像對準

以下是使用多層級影像對準架構，進行線性仿射影像對準的 C++ 範例。

#include <itkAffineTransform.h>
#include <itkCenteredTransformInitializer.h>
#include <itkMultiResolutionImageRegistrationMethod.h>
#include <itkMattesMutualInformationImageToImageMetric.h>
#include <itkRegularStepGradientDescentOptimizer.h>
#include <itkRecursiveMultiResolutionPyramidImageFilter.h>
#include <itkImage.h>
#include <itkImageFileReader.h>
#include <itkImageFileWriter.h>
#include <itkResampleImageFilter.h>
#include <itkCastImageFilter.h>
#include <itkCheckerBoardImageFilter.h>
#include <itkCommand.h>

// 定義與實作觀察影像對準過程的 Observer
class CommandIterationUpdate : public itk::Command {
public:
  using Self = CommandIterationUpdate;
  using Superclass = itk::Command;
  using Pointer = itk::SmartPointer<Self>;
  itkNewMacro(Self);

protected:
  CommandIterationUpdate() = default;

public:
  using OptimizerType = itk::RegularStepGradientDescentOptimizer;
  using OptimizerPointer = const OptimizerType *;

  void Execute(itk::Object * caller, const itk::EventObject & event) override {
    Execute((const itk::Object *)caller, event);
  }

  void Execute(const itk::Object * object, const itk::EventObject & event) override {
    auto optimizer = static_cast<OptimizerPointer>(object);
    if (!(itk::IterationEvent().CheckEvent(&event))) {
      return;
    }
    std::cout << std::fixed;
    std::cout << optimizer->GetCurrentIteration() << "   ";
    std::cout << optimizer->GetValue() << "   ";
    std::cout << optimizer->GetCurrentPosition() << "  "
              << m_CumulativeIterationIndex++ << std::endl;
  }

private:
  unsigned int m_CumulativeIterationIndex{ 0 };
};

// 定義與實作在解析度變動時控制影像對準參數的 Observer
template <typename TRegistration>
class RegistrationInterfaceCommand : public itk::Command {
public:
  using Self = RegistrationInterfaceCommand;
  using Superclass = itk::Command;
  using Pointer = itk::SmartPointer<Self>;
  itkNewMacro(Self);

protected:
  RegistrationInterfaceCommand() = default;

public:
  using RegistrationType = TRegistration;
  using RegistrationPointer = RegistrationType *;
  using OptimizerType = itk::RegularStepGradientDescentOptimizer;
  using OptimizerPointer = OptimizerType *;
  void Execute(itk::Object * object, const itk::EventObject & event) override {
    if (!(itk::IterationEvent().CheckEvent(&event))) {
      return;
    }
    auto registration = static_cast<RegistrationPointer>(object);
    auto optimizer =
      static_cast<OptimizerPointer>(registration->GetModifiableOptimizer());

    std::cout << "-------------------------------------" << std::endl;
    std::cout << "MultiResolution Level : " << registration->GetCurrentLevel()
              << std::endl;
    std::cout << std::endl;

    // 自訂步長規則
    if (registration->GetCurrentLevel() == 0) {
      optimizer->SetMaximumStepLength(16.00);
      optimizer->SetMinimumStepLength(0.01);
    } else {
      optimizer->SetMaximumStepLength(optimizer->GetMaximumStepLength() / 4.0);
      optimizer->SetMinimumStepLength(optimizer->GetMinimumStepLength() / 10.0);
    }
  }
  void Execute(const itk::Object *, const itk::EventObject &) override {
    return;
  }
};

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

  // 影像維度
  constexpr unsigned int Dimension = 2;

  // 影像像素資料型態
  using PixelType = float;

  // 影像型態
  using ImageType = itk::Image<PixelType, Dimension>;

  // 設定影像轉換類型為線性的 AffineTransform（參數型態可為 float 或 double）
  using TransformType = itk::AffineTransform<double, Dimension>;

  // Optimizer
  using OptimizerType = itk::RegularStepGradientDescentOptimizer;

  // 使用線性內插
  using InterpolatorType =
    itk::LinearInterpolateImageFunction<ImageType, double>;

  // Metric
  using MetricType =
    itk::MattesMutualInformationImageToImageMetric<ImageType, ImageType>;

  using OptimizerScalesType = OptimizerType::ScalesType;

  // 多層級影像對準
  using RegistrationType =
    itk::MultiResolutionImageRegistrationMethod<ImageType, ImageType>;

  // 建立 Optimizer、Interpolator、Registration、Metric、Transform 物件
  OptimizerType::Pointer    optimizer = OptimizerType::New();
  InterpolatorType::Pointer interpolator = InterpolatorType::New();
  RegistrationType::Pointer registration = RegistrationType::New();
  MetricType::Pointer       metric = MetricType::New();
  TransformType::Pointer    transform = TransformType::New();

  registration->SetOptimizer(optimizer);
  registration->SetInterpolator(interpolator);
  registration->SetMetric(metric);
  registration->SetTransform(transform);

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

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

  fixedImageReader->SetFileName(argv[1]);  // 讀取固定影像
  movingImageReader->SetFileName(argv[2]); // 讀取調動影像

  // 設定影像來源
  registration->SetFixedImage(fixedImageReader->GetOutput());
  registration->SetMovingImage(movingImageReader->GetOutput());

  fixedImageReader->Update();

  // 設定固定影像作用範圍
  registration->SetFixedImageRegion(
    fixedImageReader->GetOutput()->GetBufferedRegion());

  // 使用 CenteredTransformInitializer 初始化轉換參數
  using TransformInitializerType =
    itk::CenteredTransformInitializer<TransformType, ImageType, ImageType>;
  TransformInitializerType::Pointer initializer =
    TransformInitializerType::New();
  initializer->SetTransform(transform);
  initializer->SetFixedImage(fixedImageReader->GetOutput());
  initializer->SetMovingImage(movingImageReader->GetOutput());
  initializer->MomentsOn();
  initializer->InitializeTransform();
  registration->SetInitialTransformParameters(transform->GetParameters());

  // 設定各參數尺度，前 NxN 個參數為線性轉換矩陣，後 N 個參數為平移轉換參數
  OptimizerScalesType optimizerScales(transform->GetNumberOfParameters());

  optimizerScales[0] = 1.0; // 線性轉換矩陣 M11
  optimizerScales[1] = 1.0; // 線性轉換矩陣 M12
  optimizerScales[2] = 1.0; // 線性轉換矩陣 M21
  optimizerScales[3] = 1.0; // 線性轉換矩陣 M22
  optimizerScales[4] = 1.0 / 1e7; // 平移轉換參數 X
  optimizerScales[5] = 1.0 / 1e7; // 平移轉換參數 Y

  optimizer->SetScales(optimizerScales);

  // MattesMutualInformationImageToImageMetric 相關參數
  metric->SetNumberOfHistogramBins(128);
  metric->SetNumberOfSpatialSamples(50000);

  // 設定隨機取樣用亂數種子
  metric->ReinitializeSeed(76926294);

  // 設定 ExplicitPDFDerivatives
  // UseExplicitPDFDerivatives = True
  //   計算速度快，耗費記憶體，適用於參數量較少的狀況，例如剛性轉換
  // UseExplicitPDFDerivatives = False
  //   計算速度慢，節省記憶體，適用於參數量較少的狀況，例如 BSpline 轉換
  metric->SetUseExplicitPDFDerivatives(true);

  // 設定疊代次數上限與 Relaxation Factor
  optimizer->SetNumberOfIterations(200);
  optimizer->SetRelaxationFactor(0.8);

  // 設定觀察影像對準過程的 Observer
  CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
  optimizer->AddObserver(itk::IterationEvent(), observer);

  // 設定在解析度變動時控制影像對準參數的 Observer
  using CommandType = RegistrationInterfaceCommand<RegistrationType>;
  CommandType::Pointer command = CommandType::New();
  registration->AddObserver(itk::IterationEvent(), command);

  // 設定層級數
  registration->SetNumberOfLevels(3);

  // 進行影像對準
  try {
    registration->Update();
    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;
  }

  std::cout << "Optimizer Stopping Condition = "
            << optimizer->GetStopCondition() << std::endl;

  using ParametersType = RegistrationType::ParametersType;
  ParametersType finalParameters = registration->GetLastTransformParameters();

  double MatrixValue11 = finalParameters[0];
  double MatrixValue12 = finalParameters[1];
  double MatrixValue21 = finalParameters[2];
  double MatrixValue22 = finalParameters[3];

  double TranslationAlongX = finalParameters[4];
  double TranslationAlongY = finalParameters[5];

  unsigned int numberOfIterations = optimizer->GetCurrentIteration();

  double bestValue = optimizer->GetValue();

  // 輸出影像對準結果
  std::cout << std::fixed;
  std::cout << "Result = " << std::endl;
  std::cout << " Translation X = " << TranslationAlongX << std::endl;
  std::cout << " Translation Y = " << TranslationAlongY << std::endl;
  std::cout << " Matrix        = " << MatrixValue11 << " " << MatrixValue12 << std::endl;
  std::cout << "                 " << MatrixValue21 << " " << MatrixValue22 << std::endl;
  std::cout << " Iterations    = " << numberOfIterations << std::endl;
  std::cout << " Metric value  = " << bestValue << std::endl;

  // 建立 ResampleFilter 套用影像轉換
  using ResampleFilterType =
    itk::ResampleImageFilter<ImageType, ImageType>;

  TransformType::Pointer finalTransform = TransformType::New();

  finalTransform->SetParameters(finalParameters);
  finalTransform->SetFixedParameters(transform->GetFixedParameters());

  ResampleFilterType::Pointer resample = ResampleFilterType::New();

  resample->SetTransform(finalTransform);
  resample->SetInput(movingImageReader->GetOutput());

  ImageType::Pointer fixedImage = fixedImageReader->GetOutput();

  // 預設背景像素值
  PixelType backgroundGrayLevel = 100;

  resample->SetSize(fixedImage->GetLargestPossibleRegion().GetSize());
  resample->SetOutputOrigin(fixedImage->GetOrigin());
  resample->SetOutputSpacing(fixedImage->GetSpacing());
  resample->SetOutputDirection(fixedImage->GetDirection());
  resample->SetDefaultPixelValue(backgroundGrayLevel);

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

  // 影像轉型 CastImageFilter
  using CastFilterType =
    itk::CastImageFilter<ImageType, OutputImageType>;
  using WriterType = itk::ImageFileWriter<OutputImageType>;

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

  // 設定出檔案
  writer->SetFileName(argv[3]);

  caster->SetInput(resample->GetOutput());
  writer->SetInput(caster->GetOutput());
  writer->Update();

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

  CheckerBoardFilterType::Pointer checker = CheckerBoardFilterType::New();

  checker->SetInput1(fixedImage);
  checker->SetInput2(resample->GetOutput());

  caster->SetInput(checker->GetOutput());
  writer->SetInput(caster->GetOutput());

  resample->SetDefaultPixelValue(0);

  // 輸出影像對準前的棋盤式影像比較圖
  TransformType::Pointer identityTransform = TransformType::New();
  identityTransform->SetIdentity();
  resample->SetTransform(identityTransform);

  if (argc > 4) {
    writer->SetFileName(argv[4]);
    writer->Update();
  }

  // 輸出影像對準後的棋盤式影像比較圖
  resample->SetTransform(finalTransform);
  if (argc > 5) {
    writer->SetFileName(argv[5]);
    writer->Update();
  }

  return EXIT_SUCCESS;
}

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

cmake_minimum_required(VERSION 3.10.2)

# 設定專案名稱
project(MultiResReg)

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

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

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

編譯以執行的指令如下：

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

# 執行程式
./MultiResReg BrainT1SliceBorder20.png \
  BrainProtonDensitySliceShifted13x17y.png \
  output.png checkerboard1.png checkerboard2.png

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

3D 多層級線性影像對準

3D 版的線性影像對準基本上跟 2D 版本幾乎相同，只是更改影像維度而已，另外稍微修改一下參數尺度。

#include <itkAffineTransform.h>
#include <itkCenteredTransformInitializer.h>
#include <itkMultiResolutionImageRegistrationMethod.h>
#include <itkMattesMutualInformationImageToImageMetric.h>
#include <itkRegularStepGradientDescentOptimizer.h>
#include <itkRecursiveMultiResolutionPyramidImageFilter.h>
#include <itkImage.h>
#include <itkImageFileReader.h>
#include <itkImageFileWriter.h>
#include <itkResampleImageFilter.h>
#include <itkCastImageFilter.h>
#include <itkCheckerBoardImageFilter.h>
#include <itkCommand.h>

// 定義與實作觀察影像對準過程的 Observer
class CommandIterationUpdate : public itk::Command {
public:
  using Self = CommandIterationUpdate;
  using Superclass = itk::Command;
  using Pointer = itk::SmartPointer<Self>;
  itkNewMacro(Self);

protected:
  CommandIterationUpdate() = default;

public:
  using OptimizerType = itk::RegularStepGradientDescentOptimizer;
  using OptimizerPointer = const OptimizerType *;

  void Execute(itk::Object * caller, const itk::EventObject & event) override {
    Execute((const itk::Object *)caller, event);
  }

  void Execute(const itk::Object * object, const itk::EventObject & event) override {
    auto optimizer = static_cast<OptimizerPointer>(object);
    if (!(itk::IterationEvent().CheckEvent(&event))) {
      return;
    }
    std::cout << std::fixed;
    std::cout << optimizer->GetCurrentIteration() << "   ";
    std::cout << optimizer->GetValue() << "   ";
    std::cout << optimizer->GetCurrentPosition() << "  "
              << m_CumulativeIterationIndex++ << std::endl;
  }

private:
  unsigned int m_CumulativeIterationIndex{ 0 };
};

// 定義與實作在解析度變動時控制影像對準參數的 Observer
template <typename TRegistration>
class RegistrationInterfaceCommand : public itk::Command {
public:
  using Self = RegistrationInterfaceCommand;
  using Superclass = itk::Command;
  using Pointer = itk::SmartPointer<Self>;
  itkNewMacro(Self);

protected:
  RegistrationInterfaceCommand() = default;

public:
  using RegistrationType = TRegistration;
  using RegistrationPointer = RegistrationType *;
  using OptimizerType = itk::RegularStepGradientDescentOptimizer;
  using OptimizerPointer = OptimizerType *;
  void Execute(itk::Object * object, const itk::EventObject & event) override {
    if (!(itk::IterationEvent().CheckEvent(&event))) {
      return;
    }
    auto registration = static_cast<RegistrationPointer>(object);
    auto optimizer =
      static_cast<OptimizerPointer>(registration->GetModifiableOptimizer());

    std::cout << "-------------------------------------" << std::endl;
    std::cout << "MultiResolution Level : " << registration->GetCurrentLevel()
              << std::endl;
    std::cout << std::endl;

    // 自訂步長規則
    if (registration->GetCurrentLevel() == 0) {
      optimizer->SetMaximumStepLength(32.00);
      optimizer->SetMinimumStepLength(0.01);
    } else {
      optimizer->SetMaximumStepLength(optimizer->GetMaximumStepLength() / 4.0);
      optimizer->SetMinimumStepLength(optimizer->GetMinimumStepLength() / 10.0);
    }
  }
  void Execute(const itk::Object *, const itk::EventObject &) override {
    return;
  }
};

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

  // 影像維度
  constexpr unsigned int Dimension = 3;

  // 影像像素資料型態
  using PixelType = float;

  // 影像型態
  using ImageType = itk::Image<PixelType, Dimension>;

  // 設定影像轉換類型為線性的 AffineTransform（參數型態可為 float 或 double）
  using TransformType = itk::AffineTransform<double, Dimension>;

  // Optimizer
  using OptimizerType = itk::RegularStepGradientDescentOptimizer;

  // 使用線性內插
  using InterpolatorType =
    itk::LinearInterpolateImageFunction<ImageType, double>;

  // Metric
  using MetricType =
    itk::MattesMutualInformationImageToImageMetric<ImageType, ImageType>;

  using OptimizerScalesType = OptimizerType::ScalesType;

  // 多層級影像對準
  using RegistrationType =
    itk::MultiResolutionImageRegistrationMethod<ImageType, ImageType>;

  // 建立 Optimizer、Interpolator、Registration、Metric、Transform 物件
  OptimizerType::Pointer    optimizer = OptimizerType::New();
  InterpolatorType::Pointer interpolator = InterpolatorType::New();
  RegistrationType::Pointer registration = RegistrationType::New();
  MetricType::Pointer       metric = MetricType::New();
  TransformType::Pointer    transform = TransformType::New();

  registration->SetOptimizer(optimizer);
  registration->SetInterpolator(interpolator);
  registration->SetMetric(metric);
  registration->SetTransform(transform);

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

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

  fixedImageReader->SetFileName(argv[1]);  // 讀取固定影像
  movingImageReader->SetFileName(argv[2]); // 讀取調動影像

  // 設定影像來源
  registration->SetFixedImage(fixedImageReader->GetOutput());
  registration->SetMovingImage(movingImageReader->GetOutput());

  fixedImageReader->Update();

  // 設定固定影像作用範圍
  registration->SetFixedImageRegion(
    fixedImageReader->GetOutput()->GetBufferedRegion());

  // 使用 CenteredTransformInitializer 初始化轉換參數
  using TransformInitializerType =
    itk::CenteredTransformInitializer<TransformType, ImageType, ImageType>;
  TransformInitializerType::Pointer initializer =
    TransformInitializerType::New();
  initializer->SetTransform(transform);
  initializer->SetFixedImage(fixedImageReader->GetOutput());
  initializer->SetMovingImage(movingImageReader->GetOutput());
  initializer->MomentsOn();
  initializer->InitializeTransform();
  registration->SetInitialTransformParameters(transform->GetParameters());

  // 設定各參數尺度，前 NxN 個參數為線性轉換矩陣，後 N 個參數為平移轉換參數
  OptimizerScalesType optimizerScales(transform->GetNumberOfParameters());

  optimizerScales[0] = 1.0; // 線性轉換矩陣 M11
  optimizerScales[1] = 1.0; // 線性轉換矩陣 M12
  optimizerScales[2] = 1.0; // 線性轉換矩陣 M13
  optimizerScales[3] = 1.0; // 線性轉換矩陣 M21
  optimizerScales[4] = 1.0; // 線性轉換矩陣 M22
  optimizerScales[5] = 1.0; // 線性轉換矩陣 M23
  optimizerScales[6] = 1.0; // 線性轉換矩陣 M31
  optimizerScales[7] = 1.0; // 線性轉換矩陣 M32
  optimizerScales[8] = 1.0; // 線性轉換矩陣 M33

  optimizerScales[9]  = 1.0 / 1e7; // 平移轉換參數 X
  optimizerScales[10] = 1.0 / 1e7; // 平移轉換參數 Y
  optimizerScales[11] = 1.0 / 1e7; // 平移轉換參數 Z

  optimizer->SetScales(optimizerScales);

  // MattesMutualInformationImageToImageMetric 相關參數
  metric->SetNumberOfHistogramBins(128);
  metric->SetNumberOfSpatialSamples(50000);

  // 設定隨機取樣用亂數種子
  metric->ReinitializeSeed(76926294);

  // 設定 ExplicitPDFDerivatives
  // UseExplicitPDFDerivatives = True
  //   計算速度快，耗費記憶體，適用於參數量較少的狀況，例如剛性轉換
  // UseExplicitPDFDerivatives = False
  //   計算速度慢，節省記憶體，適用於參數量較少的狀況，例如 BSpline 轉換
  metric->SetUseExplicitPDFDerivatives(true);

  // 設定疊代次數上限與 Relaxation Factor
  optimizer->SetNumberOfIterations(200);
  optimizer->SetRelaxationFactor(0.8);

  // 設定觀察影像對準過程的 Observer
  CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
  optimizer->AddObserver(itk::IterationEvent(), observer);

  // 設定在解析度變動時控制影像對準參數的 Observer
  using CommandType = RegistrationInterfaceCommand<RegistrationType>;
  CommandType::Pointer command = CommandType::New();
  registration->AddObserver(itk::IterationEvent(), command);

  // 設定層級數
  registration->SetNumberOfLevels(3);

  // 進行影像對準
  try {
    registration->Update();
    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;
  }

  std::cout << "Optimizer Stopping Condition = "
            << optimizer->GetStopCondition() << std::endl;

  using ParametersType = RegistrationType::ParametersType;
  ParametersType finalParameters = registration->GetLastTransformParameters();

  double MatrixValue11 = finalParameters[0];
  double MatrixValue12 = finalParameters[1];
  double MatrixValue13 = finalParameters[2];
  double MatrixValue21 = finalParameters[3];
  double MatrixValue22 = finalParameters[4];
  double MatrixValue23 = finalParameters[5];
  double MatrixValue31 = finalParameters[6];
  double MatrixValue32 = finalParameters[7];
  double MatrixValue33 = finalParameters[8];

  double TranslationAlongX = finalParameters[9];
  double TranslationAlongY = finalParameters[10];
  double TranslationAlongZ = finalParameters[11];

  unsigned int numberOfIterations = optimizer->GetCurrentIteration();

  double bestValue = optimizer->GetValue();

  // 輸出影像對準結果
  std::cout << std::fixed;
  std::cout << "Result = " << std::endl;
  std::cout << " Translation X = " << TranslationAlongX << std::endl;
  std::cout << " Translation Y = " << TranslationAlongY << std::endl;
  std::cout << " Translation Z = " << TranslationAlongZ << std::endl;
  std::cout << " Matrix        = " << MatrixValue11 << " " << MatrixValue12 << " " << MatrixValue13 << std::endl;
  std::cout << "                 " << MatrixValue21 << " " << MatrixValue22 << " " << MatrixValue23 << std::endl;
  std::cout << "                 " << MatrixValue31 << " " << MatrixValue32 << " " << MatrixValue33 << std::endl;
  std::cout << " Iterations    = " << numberOfIterations << std::endl;
  std::cout << " Metric value  = " << bestValue << std::endl;

  // 建立 ResampleFilter 套用影像轉換
  using ResampleFilterType =
    itk::ResampleImageFilter<ImageType, ImageType>;

  TransformType::Pointer finalTransform = TransformType::New();

  finalTransform->SetParameters(finalParameters);
  finalTransform->SetFixedParameters(transform->GetFixedParameters());

  ResampleFilterType::Pointer resample = ResampleFilterType::New();

  resample->SetTransform(finalTransform);
  resample->SetInput(movingImageReader->GetOutput());

  ImageType::Pointer fixedImage = fixedImageReader->GetOutput();

  // 預設背景像素值
  PixelType backgroundGrayLevel = 100;

  resample->SetSize(fixedImage->GetLargestPossibleRegion().GetSize());
  resample->SetOutputOrigin(fixedImage->GetOrigin());
  resample->SetOutputSpacing(fixedImage->GetSpacing());
  resample->SetOutputDirection(fixedImage->GetDirection());
  resample->SetDefaultPixelValue(backgroundGrayLevel);

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

  // 影像轉型 CastImageFilter
  using CastFilterType =
    itk::CastImageFilter<ImageType, OutputImageType>;
  using WriterType = itk::ImageFileWriter<OutputImageType>;

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

  // 設定出檔案
  writer->SetFileName(argv[3]);

  caster->SetInput(resample->GetOutput());
  writer->SetInput(caster->GetOutput());
  writer->Update();

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

  CheckerBoardFilterType::Pointer checker = CheckerBoardFilterType::New();

  checker->SetInput1(fixedImage);
  checker->SetInput2(resample->GetOutput());

  caster->SetInput(checker->GetOutput());
  writer->SetInput(caster->GetOutput());

  resample->SetDefaultPixelValue(0);

  // 輸出影像對準前的棋盤式影像比較圖
  TransformType::Pointer identityTransform = TransformType::New();
  identityTransform->SetIdentity();
  resample->SetTransform(identityTransform);

  if (argc > 4) {
    writer->SetFileName(argv[4]);
    writer->Update();
  }

  // 輸出影像對準後的棋盤式影像比較圖
  resample->SetTransform(finalTransform);
  if (argc > 5) {
    writer->SetFileName(argv[5]);
    writer->Update();
  }

  return EXIT_SUCCESS;
}

編譯方式完全相同，而執行時則輸入 3D 的影像：

# 執行程式
./MultiResReg brainweb165a10f17.mha \
  brainweb1e1a10f20Rot10Tx15.mha \
  output.mha checkerboard1.mha checkerboard2.mha

這裡所使用的測試影像可從 Kitware 網站上取得。

參考資料：MultiResImageRegistration1.cxx、MultiResImageRegistration2.cxx、MultiResImageRegistration3.cxx