介紹如何在 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 網站取得，執行後除了可以產生對準後的圖形之外，也可以產生影像對準前後的棋盤式影像比較圖。
