ITK 多階段 Multi-Stage 影像對準架構教學與範例

介紹如何在 ITK 的多階段（multi-stage）影像對準架構下，組合平移對準（translation registration）與仿射對準（affine registration）。

2D 多階段線性影像對準

以下是使用多階段影像對準架構，將平移轉換（translation transform）與仿射轉換（affine transform）串接之後，進行影像對準的 C++ 範例。

#include <itkImageRegistrationMethodv4.h>
#include <itkMattesMutualInformationImageToImageMetricv4.h>
#include <itkRegularStepGradientDescentOptimizerv4.h>
#include <itkConjugateGradientLineSearchOptimizerv4.h>
#include <itkTranslationTransform.h>
#include <itkAffineTransform.h>
#include <itkCompositeTransform.h>
#include <itkImageFileReader.h>
#include <itkImageFileWriter.h>
#include <itkImageMomentsCalculator.h>
#include <itkResampleImageFilter.h>
#include <itkCastImageFilter.h>
#include <itkCheckerBoardImageFilter.h>
#include <itkCommand.h>

// 觀察階段或層級改變的 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;

  // 可接受 const 物件的包裝函數
  void
  Execute(itk::Object * object, const itk::EventObject & event) override {
    Execute((const itk::Object *)object, event);
  }

  void
  Execute(const itk::Object * object, const itk::EventObject & event) override {
    if (!(itk::MultiResolutionIterationEvent().CheckEvent(&event))) {
      return;
    }

    std::cout << "\nObserving from class " << object->GetNameOfClass();
    if (!object->GetObjectName().empty()) {
      std::cout << " \"" << object->GetObjectName() << "\"" << std::endl;
    }

    const auto * registration = static_cast<const RegistrationType *>(object);
    if (registration == nullptr) {
      itkExceptionMacro(<< "Dynamic cast failed, object of type "
                        << object->GetNameOfClass());
    }

    unsigned int currentLevel = registration->GetCurrentLevel();
    typename RegistrationType::ShrinkFactorsPerDimensionContainerType
      shrinkFactors =
        registration->GetShrinkFactorsPerDimension(currentLevel);
    typename RegistrationType::SmoothingSigmasArrayType smoothingSigmas =
      registration->GetSmoothingSigmasPerLevel();

    std::cout << "-------------------------------------" << std::endl;
    std::cout << " Current multi-resolution level = " << currentLevel
              << std::endl;
    std::cout << "    shrink factor = " << shrinkFactors << std::endl;
    std::cout << "    smoothing sigma = " << smoothingSigmas[currentLevel]
              << std::endl;
    std::cout << std::endl;
  }
};

// 觀察影像對準疊代過程的 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::GradientDescentOptimizerv4Template<double>;
  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 (optimizer == nullptr) {
      return; // Optimizer 類型錯誤，不做任何處理
    }
    if (!(itk::IterationEvent().CheckEvent(&event))) {
      return;
    }
    std::cout << optimizer->GetCurrentIteration() << "   ";
    std::cout << optimizer->GetValue() << "   ";
    std::cout << optimizer->GetCurrentPosition() << "  "
              << m_CumulativeIterationIndex++ << std::endl;
  }

private:
  unsigned int m_CumulativeIterationIndex{ 0 };
};

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]" << std::endl;
    return EXIT_FAILURE;
  }

  constexpr unsigned int Dimension = 2;
  using PixelType = float;

  using FixedImageType = itk::Image<PixelType, Dimension>;
  using MovingImageType = itk::Image<PixelType, Dimension>;

  // 第一階段影像對準：Translation 轉換

  using TTransformType = itk::TranslationTransform<double, Dimension>;
  using TOptimizerType = itk::RegularStepGradientDescentOptimizerv4<double>;
  using MetricType =
    itk::MattesMutualInformationImageToImageMetricv4<FixedImageType,
                                                     MovingImageType>;

  // 這裡的 ImageRegistrationMethodv4 並沒有定義轉換的類型，
  // 預設會使用 Transform 作為轉換的輸出類型
  using TRegistrationType =
    itk::ImageRegistrationMethodv4<FixedImageType, MovingImageType>;

  TOptimizerType::Pointer    transOptimizer = TOptimizerType::New();
  MetricType::Pointer        transMetric = MetricType::New();
  TRegistrationType::Pointer transRegistration = TRegistrationType::New();

  transRegistration->SetOptimizer(transOptimizer);
  transRegistration->SetMetric(transMetric);

  // 設定轉換物件
  TTransformType::Pointer translationTx = TTransformType::New();
  transRegistration->SetInitialTransform(translationTx);

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

  // 讀取影像
  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]);

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

  // 設定 Registration 名稱
  transRegistration->SetObjectName("TranslationRegistration");

  // 僅使用單一層級、低解析度之模糊影像進行影像對準
  constexpr unsigned int numberOfLevels1 = 1;

  TRegistrationType::ShrinkFactorsArrayType shrinkFactorsPerLevel1;
  shrinkFactorsPerLevel1.SetSize(numberOfLevels1);
  shrinkFactorsPerLevel1[0] = 3;

  TRegistrationType::SmoothingSigmasArrayType smoothingSigmasPerLevel1;
  smoothingSigmasPerLevel1.SetSize(numberOfLevels1);
  smoothingSigmasPerLevel1[0] = 2;

  transRegistration->SetNumberOfLevels(numberOfLevels1);
  transRegistration->SetShrinkFactorsPerLevel(shrinkFactorsPerLevel1);
  transRegistration->SetSmoothingSigmasPerLevel(smoothingSigmasPerLevel1);

  transMetric->SetNumberOfHistogramBins(24);

  // RegularStepGradientDescentOptimizerv4 相關參數
  transOptimizer->SetNumberOfIterations(200);
  transOptimizer->SetRelaxationFactor(0.5);
  transOptimizer->SetLearningRate(16);
  transOptimizer->SetMinimumStepLength(1.5);

  // 觀察影像對準疊代過程的 Observer
  CommandIterationUpdate::Pointer observer1 = CommandIterationUpdate::New();
  transOptimizer->AddObserver(itk::IterationEvent(), observer1);

  // 觀察階段或層級改變的 Observer
  using TranslationCommandType =
    RegistrationInterfaceCommand<TRegistrationType>;
  TranslationCommandType::Pointer command1 = TranslationCommandType::New();
  transRegistration->AddObserver(itk::MultiResolutionIterationEvent(),
                                 command1);

  // 第二階段影像對準：Affine 轉換

  // 定義第二階段用的相關組件類型
  using ATransformType = itk::AffineTransform<double, Dimension>;
  using AOptimizerType =
    itk::ConjugateGradientLineSearchOptimizerv4Template<double>;

  // 預設使用 Transform 作為轉換的輸出類型，
  // 避免不同類型轉換在串接上形成型別錯誤
  using ARegistrationType =
    itk::ImageRegistrationMethodv4<FixedImageType, MovingImageType>;

  AOptimizerType::Pointer    affineOptimizer = AOptimizerType::New();
  MetricType::Pointer        affineMetric = MetricType::New();
  ARegistrationType::Pointer affineRegistration = ARegistrationType::New();

  affineRegistration->SetOptimizer(affineOptimizer);
  affineRegistration->SetMetric(affineMetric);

  affineMetric->SetNumberOfHistogramBins(24);

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

  // 計算固定影像的重心，作為旋轉中心點
  using FixedImageCalculatorType =
    itk::ImageMomentsCalculator<FixedImageType>;

  FixedImageCalculatorType::Pointer fixedCalculator =
    FixedImageCalculatorType::New();
  fixedCalculator->SetImage(fixedImage);
  fixedCalculator->Compute();

  FixedImageCalculatorType::VectorType fixedCenter =
    fixedCalculator->GetCenterOfGravity();

  // 建立 Affine 轉換
  ATransformType::Pointer affineTx = ATransformType::New();

  // 設定旋轉中心點
  const unsigned int numberOfFixedParameters =
    affineTx->GetFixedParameters().Size();
  ATransformType::ParametersType fixedParameters(numberOfFixedParameters);
  for (unsigned int i = 0; i < numberOfFixedParameters; ++i) {
    fixedParameters[i] = fixedCenter[i];
  }
  affineTx->SetFixedParameters(fixedParameters);

  // 設定轉換物件
  affineRegistration->SetInitialTransform(affineTx);

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

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

  // 設定 Registration 名稱
  affineRegistration->SetObjectName("AffineRegistration");

  // 前一個階段的輸出會以 DataObjectDecorator 包裝，
  // 此處以 SetMovingInitialTransformInput 設定為調動影像初始轉換
  affineRegistration->SetMovingInitialTransformInput(
    transRegistration->GetTransformOutput());

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

  // GradientDescentLineSearchOptimizerv4 相關參數
  affineOptimizer->SetScalesEstimator(scalesEstimator);
  affineOptimizer->SetDoEstimateLearningRateOnce(true);
  affineOptimizer->SetDoEstimateLearningRateAtEachIteration(false);
  affineOptimizer->SetLowerLimit(0);
  affineOptimizer->SetUpperLimit(2);
  affineOptimizer->SetEpsilon(0.2);
  affineOptimizer->SetNumberOfIterations(200);
  affineOptimizer->SetMinimumConvergenceValue(1e-6);
  affineOptimizer->SetConvergenceWindowSize(10);

  // 觀察影像對準疊代過程的 Observer
  CommandIterationUpdate::Pointer observer2 = CommandIterationUpdate::New();
  affineOptimizer->AddObserver(itk::IterationEvent(), observer2);

  // 使用兩個層級進行影像對準
  constexpr unsigned int numberOfLevels2 = 2;

  ARegistrationType::ShrinkFactorsArrayType shrinkFactorsPerLevel2;
  shrinkFactorsPerLevel2.SetSize(numberOfLevels2);
  shrinkFactorsPerLevel2[0] = 2;
  shrinkFactorsPerLevel2[1] = 1;

  ARegistrationType::SmoothingSigmasArrayType smoothingSigmasPerLevel2;
  smoothingSigmasPerLevel2.SetSize(numberOfLevels2);
  smoothingSigmasPerLevel2[0] = 1;
  smoothingSigmasPerLevel2[1] = 0;

  affineRegistration->SetNumberOfLevels(numberOfLevels2);
  affineRegistration->SetShrinkFactorsPerLevel(shrinkFactorsPerLevel2);
  affineRegistration->SetSmoothingSigmasPerLevel(smoothingSigmasPerLevel2);

  // 觀察階段或層級改變的 Observer
  using AffineCommandType = RegistrationInterfaceCommand<ARegistrationType>;
  AffineCommandType::Pointer command2 = AffineCommandType::New();
  affineRegistration->AddObserver(itk::MultiResolutionIterationEvent(),
                                  command2);

  // 實際進行兩階段的影像對準
  try {
    affineRegistration->Update();
    std::cout
      << "Optimizer stop condition: "
      << affineRegistration->GetOptimizer()->GetStopConditionDescription()
      << std::endl;
  } catch (itk::ExceptionObject & err) {
    std::cout << "ExceptionObject caught !" << std::endl;
    std::cout << err << std::endl;
    return EXIT_FAILURE;
  }

  // 輸出個階段影像對準結果
  std::cout << " Translation transform parameters after registration: "
            << std::endl
            << transOptimizer->GetCurrentPosition() << std::endl
            << " Last LearningRate: "
            << transOptimizer->GetCurrentStepLength() << std::endl;

  std::cout << " Affine transform parameters after registration: "
            << std::endl
            << affineOptimizer->GetCurrentPosition() << std::endl
            << " Last LearningRate: " << affineOptimizer->GetLearningRate()
            << std::endl;

  // 以 CompositeTransform 組合各階段影像對準結果
  using CompositeTransformType = itk::CompositeTransform<double, Dimension>;
  CompositeTransformType::Pointer compositeTransform =
    CompositeTransformType::New();
  compositeTransform->AddTransform(translationTx);
  compositeTransform->AddTransform(affineTx);

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

  // 設定影像轉換
  resample->SetTransform(compositeTransform);

  // 設定輸入影像
  resample->SetInput(movingImageReader->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>;
  using CastFilterType =
    itk::CastImageFilter<FixedImageType, 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<FixedImageType>;

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

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

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

  // 預設像素值
  resample->SetDefaultPixelValue(0);

  // 輸出影像對準前的棋盤式影像比較圖
  using TransformType = itk::IdentityTransform<double, Dimension>;
  TransformType::Pointer identityTransform;
  try {
    identityTransform = TransformType::New();
  } catch (const itk::ExceptionObject & err) {
    err.Print(std::cerr);
    return EXIT_FAILURE;
  }
  identityTransform->SetIdentity();
  resample->SetTransform(identityTransform);

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

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

  return EXIT_SUCCESS;
}

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

cmake_minimum_required(VERSION 3.10.2)

# 設定專案名稱
project(MultiStageReg)

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

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

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

編譯與執行的指令如下：

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

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

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

參考資料：MultiStageImageRegistration1.cxx、MultiStageImageRegistration2.cxx