介紹如何在 ITK 的多階段(multi-stage)影像對準架構下,組合平移對準(translation registration)與仿射對準(affine registration)。
以下是使用多階段影像對準架構,將平移轉換(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