itk registration 13

参考：https://itk.org/ITKSoftwareGuide/html/Book2/ITKSoftwareGuide-Book2ch3.html#x26-1740003.18

3.13 Deformable Registration

3.13.1 FEM-Based Image Registration

https://github.com/InsightSoftwareConsortium/ITK/blob/master/Examples/RegistrationITKv4/DeformableRegistration1.cxx

#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"#include "itkRescaleIntensityImageFilter.h"
#include "itkHistogramMatchingImageFilter.h"#include "itkFEMRegistrationFilter.h"using DiskImageType = itk::Image<unsigned char, 2>;//接下来，我们using类型别名实例化所有必要的类。我们定义了计划用于解决二维配准问题的图像和元素类型。我们定义了多个元素类型，以便在不重新编译代码的情况下使用它们。
using ImageType = itk::Image<float, 2>;
using ElementType = itk::fem::Element2DC0LinearQuadrilateralMembrane;
using ElementType2 = itk::fem::Element2DC0LinearTriangularMembrane;
using FEMObjectType = itk::fem::FEMObject<2>;using FileImage3DType = itk::Image<unsigned char, 3>;//注意，为了解决三维配准问题，我们将简单地定义3D图像和元素类型来代替上述类型。以下声明可用于3D问题：
using Image3DType = itk::Image<float, 3>;
using Element3DType = itk::fem::Element3DC0LinearHexahedronMembrane;
using Element3DType2 = itk::fem::Element3DC0LinearTetrahedronMembrane;
using FEMObject3DType = itk::fem::FEMObject<3>;using RegistrationType =itk::fem::FEMRegistrationFilter<ImageType, ImageType, FEMObjectType>;//一旦实例化了所有必需的组件，就可以实例化itk :: FEMRegistrationFilter，这取决于图像的输入和输出类型。int
main(int argc, char * argv[])
{const char *fixedImageName, *movingImageName;if (argc < 2){std::cout << "Image file names missing" << std::endl;std::cout << "Usage: " << argv[0] << " fixedImageFile movingImageFile"<< std::endl;return EXIT_FAILURE;}else{fixedImageName = argv[1];movingImageName = argv[2];}RegistrationType::Pointer registrationFilter = RegistrationType::New();//为了开始配准，我们声明FEMRegistrationFilter的实例并设置其参数。为简单起见，我们将其称为registrationFilter。registrationFilter->SetMaxLevel(1);registrationFilter->SetUseNormalizedGradient(true);registrationFilter->ChooseMetric(0);unsigned int maxiters = 20;float        E = 100;float        p = 1;registrationFilter->SetElasticity(E, 0);registrationFilter->SetRho(p, 0);registrationFilter->SetGamma(1., 0);registrationFilter->SetAlpha(1.);registrationFilter->SetMaximumIterations(maxiters, 0);registrationFilter->SetMeshPixelsPerElementAtEachResolution(4, 0);registrationFilter->SetWidthOfMetricRegion(1, 0);registrationFilter->SetNumberOfIntegrationPoints(2, 0);registrationFilter->SetDoLineSearchOnImageEnergy(0);registrationFilter->SetTimeStep(1.);registrationFilter->SetEmployRegridding(false);registrationFilter->SetUseLandmarks(false);using FileSourceType = itk::ImageFileReader<DiskImageType>;FileSourceType::Pointer movingfilter = FileSourceType::New();movingfilter->SetFileName(movingImageName);FileSourceType::Pointer fixedfilter = FileSourceType::New();fixedfilter->SetFileName(fixedImageName);std::cout << " reading moving " << movingImageName << std::endl;std::cout << " reading fixed " << fixedImageName << std::endl;try{movingfilter->Update();}catch (const itk::ExceptionObject & e){std::cerr << "Exception caught during reference file reading "<< std::endl;std::cerr << e << std::endl;return EXIT_FAILURE;}try{fixedfilter->Update();}catch (const itk::ExceptionObject & e){std::cerr << "Exception caught during target file reading " << std::endl;std::cerr << e << std::endl;return EXIT_FAILURE;}using FilterType =itk::RescaleIntensityImageFilter<DiskImageType, ImageType>;FilterType::Pointer movingrescalefilter = FilterType::New();FilterType::Pointer fixedrescalefilter = FilterType::New();movingrescalefilter->SetInput(movingfilter->GetOutput());fixedrescalefilter->SetInput(fixedfilter->GetOutput());constexpr double desiredMinimum = 0.0;constexpr double desiredMaximum = 255.0;movingrescalefilter->SetOutputMinimum(desiredMinimum);movingrescalefilter->SetOutputMaximum(desiredMaximum);movingrescalefilter->UpdateLargestPossibleRegion();fixedrescalefilter->SetOutputMinimum(desiredMinimum);fixedrescalefilter->SetOutputMaximum(desiredMaximum);fixedrescalefilter->UpdateLargestPossibleRegion();using HEFilterType =itk::HistogramMatchingImageFilter<ImageType, ImageType>;HEFilterType::Pointer IntensityEqualizeFilter = HEFilterType::New();IntensityEqualizeFilter->SetReferenceImage(fixedrescalefilter->GetOutput());IntensityEqualizeFilter->SetInput(movingrescalefilter->GetOutput());IntensityEqualizeFilter->SetNumberOfHistogramLevels(100);IntensityEqualizeFilter->SetNumberOfMatchPoints(15);IntensityEqualizeFilter->ThresholdAtMeanIntensityOn();IntensityEqualizeFilter->Update();registrationFilter->SetFixedImage(fixedrescalefilter->GetOutput());registrationFilter->SetMovingImage(IntensityEqualizeFilter->GetOutput());itk::ImageFileWriter<ImageType>::Pointer writer;writer = itk::ImageFileWriter<ImageType>::New();std::string ofn = "fixed.mha";writer->SetFileName(ofn.c_str());writer->SetInput(registrationFilter->GetFixedImage());try{writer->Write();}catch (const itk::ExceptionObject & excp){std::cerr << excp << std::endl;return EXIT_FAILURE;}ofn = "moving.mha";itk::ImageFileWriter<ImageType>::Pointer writer2;writer2 = itk::ImageFileWriter<ImageType>::New();writer2->SetFileName(ofn.c_str());writer2->SetInput(registrationFilter->GetMovingImage());try{writer2->Write();}catch (const itk::ExceptionObject & excp){std::cerr << excp << std::endl;return EXIT_FAILURE;}itk::fem::MaterialLinearElasticity::Pointer m;//为了初始化网格元素，我们必须首先创建“dummy”material和元素对象，并将它们分配给配准过滤器。这些对象随后用于从文件中读取预定义的网格或使用软件生成网格。分配给material对象中的字段的值是任意的，因为它们将被前面指定的值替换。类似地，元素对象将被替换为来自所需网格的对象。m = itk::fem::MaterialLinearElasticity::New();    // Create the material propertiesm->SetGlobalNumber(0);m->SetYoungsModulus(registrationFilter->GetElasticity());// Young's modulus of the membranem->SetCrossSectionalArea(1.0); // Cross-sectional aream->SetThickness(1.0);          // Thicknessm->SetMomentOfInertia(1.0);    // Moment of inertiam->SetPoissonsRatio(0.);       // Poisson's ratio -- DONT CHOOSE 1.0!!m->SetDensityHeatProduct(1.0); // Density-Heat capacity productElementType::Pointer e1 = ElementType::New();// Create the element typee1->SetMaterial(m);registrationFilter->SetElement(e1);registrationFilter->SetMaterial(m);registrationFilter->RunRegistration();//运行配准：itk::ImageFileWriter<ImageType>::Pointer warpedImageWriter;//要输出配准产生的图像，我们可以调用GetWarpedImage（）。图像以浮点格式写入。warpedImageWriter = itk::ImageFileWriter<ImageType>::New();warpedImageWriter->SetInput(registrationFilter->GetWarpedImage());warpedImageWriter->SetFileName("warpedMovingImage.mha");try{warpedImageWriter->Update();}catch (const itk::ExceptionObject & excp){std::cerr << excp << std::endl;return EXIT_FAILURE;}using DispWriterType = itk::ImageFileWriter<RegistrationType::FieldType>;//可以输出配准产生的位移场；我们可以调用GetDisplacementField（）来获取多分量图像。DispWriterType::Pointer dispWriter = DispWriterType::New();dispWriter->SetInput(registrationFilter->GetDisplacementField());dispWriter->SetFileName("displacement.mha");try{dispWriter->Update();}catch (const itk::ExceptionObject & excp){std::cerr << excp << std::endl;return EXIT_FAILURE;}return EXIT_SUCCESS;
}

图3.46：基于FEM的可变形配准之前和之后的棋盘格比较。

图3.46显示了对活鼠数据集的两个时间间隔切片应用基于FEM的形变配准的结果。两幅图像在配准前(左)和配准后(右)的棋盘图比较。两张图像均取自同一只活鼠，第一张是在吸入空气后，第二张是在呼气后。膈肌和肋间肌的松弛导致肺组织变形，这两种肌肉都对肺组织施加压力，导致空气被排出。

下面是一个文档化的示例参数文件，可用于此可变形注册示例。这个示例演示了不使用多分辨率策略的基本注册问题的设置。因此，参数只需要在(每个元素像素的#)和(最大迭代)之间设置一个值。为了使用多分辨率策略，您必须在金字塔的每一层为这些参数指定值。

3.13.2 BSplines图像配准

https://github.com/InsightSoftwareConsortium/ITK/blob/master/Examples/RegistrationITKv4/DeformableRegistration4.cxx

这个例子演示了如何使用itk::BSplineTransform类在ITKv4注册框架中对两个2D图像进行注册。由于b样条变换的参数很多，我们将使用itk::LBFGSOptimizerv4，而不是简单的最速下降或共轭梯度下降优化器。

BSplineTransform的参数空间由与BSpline网格节点相关的所有形变集合组成。大量的参数使得它可以表示各种各样的变形，但需要大量的计算时间。

#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"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 << " [differenceOutputfile] [differenceBeforeRegistration] ";std::cerr << " [deformationField] ";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;//现在，我们使用坐标表示的类型，空间的尺寸以及BSpline的顺序作为模板参数，实例化BSplineTransform的类型。constexpr unsigned int SplineOrder = 3;using CoordinateRepType = double;using TransformType =itk::BSplineTransform<CoordinateRepType, SpaceDimension, SplineOrder>;using OptimizerType = itk::LBFGSOptimizerv4;using MetricType =itk::CorrelationImageToImageMetricv4<FixedImageType, MovingImageType>;using RegistrationType =itk::ImageRegistrationMethodv4<FixedImageType, MovingImageType>;MetricType::Pointer       metric = MetricType::New();OptimizerType::Pointer    optimizer = OptimizerType::New();RegistrationType::Pointer registration = RegistrationType::New();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();TransformType::Pointer transform = TransformType::New();//变换对象构造。using InitializerType =itk::BSplineTransformInitializer<TransformType, FixedImageType>;//bBSpline变换的fixed参数应该在配准之前定义。这些参数定义了变换网格的origin、dimension、direction和mesh size，并根据fixed的图像空间网格的规格进行设置。我们可以使用itk::BSplineTransformInitializer来初始化一个bBSpline变换的fixed参数。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();//在设置了变换的fixed参数之后，我们将初始变换设置为恒等变换。就像在创建的参数空间中将所有变换参数设置为零一样。std::cout << "Initial Parameters = " << std::endl;std::cout << transform->GetParameters() << std::endl;registration->SetInitialTransform(transform);registration->InPlaceOn();//调用InPlaceOn()意味着当前初始化的变换将直接优化并嫁接到输出，因此可以将其视为输出变换对象。否则，初始变换将被复制或“克隆”到输出变换对象，并在配准过程中对复制的对象进行优化。registration->SetFixedImage(fixedImage);registration->SetMovingImage(movingImageReader->GetOutput());using ScalesEstimatorType =itk::RegistrationParameterScalesFromPhysicalShift<MetricType>;//itk :: RegistrationParameterScalesFromPhysicalShift类用于在设置优化程序之前估算参数scales。ScalesEstimatorType::Pointer scalesEstimator = ScalesEstimatorType::New();scalesEstimator->SetMetric(metric);scalesEstimator->SetTransformForward(true);scalesEstimator->SetSmallParameterVariation(1.0);optimizer->SetGradientConvergenceTolerance(5e-2);//现在，将scalesEstimator传递给itk :: LBFGSOptimizerv4，我们还设置了优化器的其他参数。optimizer->SetLineSearchAccuracy(1.2);optimizer->SetDefaultStepLength(1.5);optimizer->TraceOn();optimizer->SetMaximumNumberOfFunctionEvaluations(1000);optimizer->SetScalesEstimator(scalesEstimator);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 (const 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();//变换对象在配准过程中更新，并传递给重采样器，以将moving图像空间映射到fixed图像空间。std::cout << "Last Transform Parameters" << std::endl;std::cout << finalParameters << std::endl;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>;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());try{writer->Update();}catch (const itk::ExceptionObject & err){std::cerr << "ExceptionObject caught !" << std::endl;std::cerr << err << std::endl;return EXIT_FAILURE;}using DifferenceFilterType =itk::SquaredDifferenceImageFilter<FixedImageType,FixedImageType,OutputImageType>;DifferenceFilterType::Pointer difference = DifferenceFilterType::New();WriterType::Pointer writer2 = WriterType::New();writer2->SetInput(difference->GetOutput());if (argc >= 5){difference->SetInput1(fixedImageReader->GetOutput());difference->SetInput2(resample->GetOutput());writer2->SetFileName(argv[4]);try{writer2->Update();}catch (const itk::ExceptionObject & err){std::cerr << "ExceptionObject caught !" << std::endl;std::cerr << err << std::endl;return EXIT_FAILURE;}}if (argc >= 6){writer2->SetFileName(argv[5]);difference->SetInput1(fixedImageReader->GetOutput());difference->SetInput2(movingImageReader->GetOutput());try{writer2->Update();}catch (const itk::ExceptionObject & err){std::cerr << "ExceptionObject caught !" << std::endl;std::cerr << err << std::endl;return EXIT_FAILURE;}}using VectorPixelType = itk::Vector<float, ImageDimension>;using DisplacementFieldImageType =itk::Image<VectorPixelType, ImageDimension>;using DisplacementFieldGeneratorType =itk::TransformToDisplacementFieldFilter<DisplacementFieldImageType,CoordinateRepType>;DisplacementFieldGeneratorType::Pointer dispfieldGenerator =DisplacementFieldGeneratorType::New();dispfieldGenerator->UseReferenceImageOn();dispfieldGenerator->SetReferenceImage(fixedImage);dispfieldGenerator->SetTransform(transform);try{dispfieldGenerator->Update();}catch (const itk::ExceptionObject & err){std::cerr << "Exception detected while generating deformation field";std::cerr << " : " << err << std::endl;return EXIT_FAILURE;}using FieldWriterType = itk::ImageFileWriter<DisplacementFieldImageType>;FieldWriterType::Pointer fieldWriter = FieldWriterType::New();fieldWriter->SetInput(dispfieldGenerator->GetOutput());if (argc >= 7){fieldWriter->SetFileName(argv[6]);try{fieldWriter->Update();}catch (const itk::ExceptionObject & excp){std::cerr << "Exception thrown " << std::endl;std::cerr << excp << std::endl;return EXIT_FAILURE;}}return EXIT_SUCCESS;
}

3.13.3 Level Set Motion for Deformable Registration

https://github.com/InsightSoftwareConsortium/ITK/blob/master/Examples/RegistrationITKv4/DeformableRegistration5.cxx

#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"#include "itkImageRegionIterator.h"#include "itkLevelSetMotionRegistrationFilter.h"//本示例演示了如何使用level set motion来可变形地配准两个图像。第一步是包含头文件。
#include "itkHistogramMatchingImageFilter.h"
#include "itkCastImageFilter.h"
#include "itkDisplacementFieldTransform.h"
#include "itkResampleImageFilter.h"class CommandIterationUpdate : public itk::Command
{
public:using Self = CommandIterationUpdate;using Superclass = itk::Command;using Pointer = itk::SmartPointer<CommandIterationUpdate>;itkNewMacro(CommandIterationUpdate);protected:CommandIterationUpdate() = default;using InternalImageType = itk::Image<float, 2>;using VectorPixelType = itk::Vector<float, 2>;using DisplacementFieldType = itk::Image<VectorPixelType, 2>;using RegistrationFilterType =itk::LevelSetMotionRegistrationFilter<InternalImageType,InternalImageType,DisplacementFieldType>;public:voidExecute(itk::Object * caller, const itk::EventObject & event) override{Execute((const itk::Object *)caller, event);}voidExecute(const itk::Object * object, const itk::EventObject & event) override{const auto * filter = static_cast<const RegistrationFilterType *>(object);if (filter == nullptr){return;}if (!(itk::IterationEvent().CheckEvent(&event))){return;}std::cout << filter->GetMetric() << std::endl;}
};int
main(int argc, char * argv[])
{if (argc < 4){std::cerr << "Missing Parameters " << std::endl;std::cerr << "Usage: " << argv[0];std::cerr << " fixedImageFile movingImageFile ";std::cerr << " outputImageFile " << std::endl;std::cerr << " [outputDisplacementFieldFile] " << std::endl;return EXIT_FAILURE;}constexpr unsigned int Dimension = 2;//其次，我们声明图像的类型。using PixelType = unsigned short;using FixedImageType = itk::Image<PixelType, Dimension>;using MovingImageType = itk::Image<PixelType, Dimension>;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]);using InternalPixelType = float;//图像文件阅读器的设置与前面的示例类似。为了支持moving图像强度的re-mapping，我们声明了一个内部图像类型，该内部图像类型使用浮点像素类型，并将输入图像转换为内部图像类型。using InternalImageType = itk::Image<InternalPixelType, Dimension>;using FixedImageCasterType =itk::CastImageFilter<FixedImageType, InternalImageType>;using MovingImageCasterType =itk::CastImageFilter<MovingImageType, InternalImageType>;FixedImageCasterType::Pointer fixedImageCaster =FixedImageCasterType::New();MovingImageCasterType::Pointer movingImageCaster =MovingImageCasterType::New();fixedImageCaster->SetInput(fixedImageReader->GetOutput());movingImageCaster->SetInput(movingImageReader->GetOutput());using MatchingFilterType =itk::HistogramMatchingImageFilter<InternalImageType, InternalImageType>;//level set motion算法依赖于这样的假设，即表示物体上相同同源点的像素对待配准的fixed图像和moving图像具有相同的强度。在本例中，我们将对moving图像进行预处理，使用itk::HistogramMatchingImageFilter来匹配图像之间的强度。
//基本思想是按照用户指定的分位数值数目匹配两个图像的直方图。为了鲁棒性，对直方图进行匹配，从而将背景像素从两个直方图中排除。对于MR图像，一个简单的步骤是排除所有小于图像平均灰度值的灰度值。MatchingFilterType::Pointer matcher = MatchingFilterType::New();matcher->SetInput(movingImageCaster->GetOutput());//对于此示例，我们将moving图像设置为source图像或input图像，将fixed图像设置为reference图像。matcher->SetReferenceImage(fixedImageCaster->GetOutput());matcher->SetNumberOfHistogramLevels(1024);//然后，我们选择要表示直方图的bin数量以及要匹配直方图的点或分位数的数量。matcher->SetNumberOfMatchPoints(7);matcher->ThresholdAtMeanIntensityOn();//简单的背景提取通过在平均强度处进行阈值处理来完成。using VectorPixelType = itk::Vector<float, Dimension>;//在itk :: LevelSetMotionRegistrationFilter中，deformation表示为图像，其像素为浮点向量。using DisplacementFieldType = itk::Image<VectorPixelType, Dimension>;using RegistrationFilterType =itk::LevelSetMotionRegistrationFilter<InternalImageType,InternalImageType,DisplacementFieldType>;RegistrationFilterType::Pointer filter = RegistrationFilterType::New();CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();filter->AddObserver(itk::IterationEvent(), observer);filter->SetFixedImage(fixedImageCaster->GetOutput());//输入的fixed图像只是fixed图像casting滤波器的输出。输入的moving图像是直方图匹配滤波器的输出。filter->SetMovingImage(matcher->GetOutput());filter->SetNumberOfIterations(50);//level set motion配准滤波器具有两个参数：要执行的迭代次数和在计算梯度之前要应用于图像的高斯平滑核的标准偏差。filter->SetGradientSmoothingStandardDeviations(4);filter->Update();//通过更新过滤器触发配准算法。过滤器输出是计算出的变形场。using InterpolatorPrecisionType = double;using TransformPrecisionType = float;using WarperType = itk::ResampleImageFilter<MovingImageType,MovingImageType,InterpolatorPrecisionType,TransformPrecisionType>;using InterpolatorType =itk::LinearInterpolateImageFunction<MovingImageType,InterpolatorPrecisionType>;WarperType::Pointer       warper = WarperType::New();//WarpImageFilter可以用输出deformation来warp moving图像。像itk::ResampleImageFilter一样，WarpImageFilter要求对输入图像进行重新采样，输入图像interpolator，以及输出图像的spacing和origin。InterpolatorType::Pointer interpolator = InterpolatorType::New();FixedImageType::Pointer   fixedImage = fixedImageReader->GetOutput();warper->SetInput(movingImageReader->GetOutput());warper->SetInterpolator(interpolator);warper->UseReferenceImageOn();warper->SetReferenceImage(fixedImage);using DisplacementFieldTransformType =itk::DisplacementFieldTransform<TransformPrecisionType, Dimension>;auto displacementTransform = DisplacementFieldTransformType::New();displacementTransform->SetDisplacementField(filter->GetOutput());warper->SetTransform(displacementTransform);//与ResampleImageFilter不同，WarpImageFilter根据矢量图像所表示的变形场来扭曲或变换输入图像。所产生的warped或重新采样的图像将按照前面的示例写入文件。warper->SetDisplacementField( filter->GetOutput() );using OutputPixelType = unsigned char;using OutputImageType = itk::Image<OutputPixelType, Dimension>;using CastFilterType =itk::CastImageFilter<MovingImageType, OutputImageType>;using WriterType = itk::ImageFileWriter<OutputImageType>;WriterType::Pointer     writer = WriterType::New();CastFilterType::Pointer caster = CastFilterType::New();writer->SetFileName(argv[3]);caster->SetInput(warper->GetOutput());writer->SetInput(caster->GetOutput());writer->Update();if (argc > 4) // if a fourth line argument has been provided...{using FieldWriterType = itk::ImageFileWriter<DisplacementFieldType>;//将变形场写为矢量图像。可以使用以下代码完成。FieldWriterType::Pointer fieldWriter = FieldWriterType::New();fieldWriter->SetFileName(argv[4]);fieldWriter->SetInput(filter->GetOutput());fieldWriter->Update();}if (argc > 5) // if a fifth line argument has been provided...{using VectorImage2DType = DisplacementFieldType;using Vector2DType = DisplacementFieldType::PixelType;VectorImage2DType::ConstPointer vectorImage2D = filter->GetOutput();VectorImage2DType::RegionType region2D =vectorImage2D->GetBufferedRegion();VectorImage2DType::IndexType index2D = region2D.GetIndex();VectorImage2DType::SizeType  size2D = region2D.GetSize();using Vector3DType = itk::Vector<float, 3>;using VectorImage3DType = itk::Image<Vector3DType, 3>;using VectorImage3DWriterType = itk::ImageFileWriter<VectorImage3DType>;VectorImage3DWriterType::Pointer writer3D =VectorImage3DWriterType::New();VectorImage3DType::Pointer vectorImage3D = VectorImage3DType::New();VectorImage3DType::RegionType region3D;VectorImage3DType::IndexType  index3D;VectorImage3DType::SizeType   size3D;index3D[0] = index2D[0];index3D[1] = index2D[1];index3D[2] = 0;size3D[0] = size2D[0];size3D[1] = size2D[1];size3D[2] = 1;region3D.SetSize(size3D);region3D.SetIndex(index3D);VectorImage2DType::SpacingType spacing2D = vectorImage2D->GetSpacing();VectorImage3DType::SpacingType spacing3D;spacing3D[0] = spacing2D[0];spacing3D[1] = spacing2D[1];spacing3D[2] = 1.0;vectorImage3D->SetSpacing(spacing3D);vectorImage3D->SetRegions(region3D);vectorImage3D->Allocate();using Iterator2DType = itk::ImageRegionConstIterator<VectorImage2DType>;using Iterator3DType = itk::ImageRegionIterator<VectorImage3DType>;Iterator2DType it2(vectorImage2D, region2D);Iterator3DType it3(vectorImage3D, region3D);it2.GoToBegin();it3.GoToBegin();Vector2DType vector2D;Vector3DType vector3D;vector3D[2] = 0; // set Z component to zero.while (!it2.IsAtEnd()){vector2D = it2.Get();vector3D[0] = vector2D[0];vector3D[1] = vector2D[1];it3.Set(vector3D);++it2;++it3;}writer3D->SetInput(vectorImage3D);writer3D->SetFileName(argv[5]);try{writer3D->Update();}catch (const itk::ExceptionObject & excp){std::cerr << excp << std::endl;return EXIT_FAILURE;}}return EXIT_SUCCESS;
}

3.13.4 BSplines Multi-Grid Image Registration

https://github.com/InsightSoftwareConsortium/ITK/blob/master/Examples/RegistrationITKv4/DeformableRegistration6.cxx

这个示例演示了在多分辨率方案中使用itk::BSplineTransform类。这里我们运行了3个层次的分辨率。第一级配准是用低分辨率的样条网格进行的。然后，一个常见的做法是在每一层增加b样条网格的分辨率(或者，类似地，控制点网格的大小)。

为此，我们介绍了transform适配器的概念。每个阶段的每个级别都由transform适配器定义，它描述了如何通过增加前一级别的分辨率来将变换适应到当前级别。在这里，我们使用itk::BSplineTransformParametersAdaptor类在每个分辨率级别调整BSpline变换参数。注意，对于许多变换，例如仿射，适配器的概念可能是荒谬的，因为变换参数的数量在分辨率级别之间不会改变。

这个示例使用itk::LBFGS2Optimizerv4，它是quasi-Newtown无界limited-memory Broyden Fletcher Goldfarb Shannon (LBFGS)优化器的新实现。无边界版本不需要指定参数空间的边界，因为参数的数量在每次B-Spline解析时都会改变，所以这种实现是首选的。

由于这个示例与前面使用BSplineTransform的示例非常相似，因此我们省略了已经讨论过的大部分细节，而将重点放在与多分辨率方法相关的方面。

#include "itkImageRegistrationMethodv4.h"
#include "itkMeanSquaresImageToImageMetricv4.h"#include "itkBSplineTransform.h"//头文件
#include "itkLBFGS2Optimizerv4.h"
#include "itkBSplineTransformParametersAdaptor.h"#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"#include "itkResampleImageFilter.h"
#include "itkCastImageFilter.h"
#include "itkSquaredDifferenceImageFilter.h"#include "itkIdentityTransform.h"#include "itkBSplineTransformInitializer.h"
#include "itkTransformToDisplacementFieldFilter.h"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::LBFGS2Optimizerv4;using OptimizerPointer = const OptimizerType *;voidExecute(itk::Object * caller, const itk::EventObject & event) override{Execute((const itk::Object *)caller, event);}voidExecute(const itk::Object * object, const itk::EventObject & event) override{auto optimizer = static_cast<OptimizerPointer>(object);if (!itk::IterationEvent().CheckEvent(&event)){return;}std::cout << optimizer->GetCurrentIteration() << " = ";std::cout << optimizer->GetValue() << std::endl;std::cout << "\t" << optimizer->GetCurrentGradientNorm() << " "<< optimizer->GetCurrentParameterNorm() << " "<< optimizer->GetCurrentStepSize() << std::endl;}
};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 << " [differenceOutputfile] [differenceBeforeRegistration] ";std::cerr << " [deformationField] ";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;//我们使用坐标表示的类型，空间的尺寸以及BSpline的顺序作为模板参数来实例化BSplineTransform的类型。constexpr unsigned int SplineOrder = 3;using CoordinateRepType = double;using TransformType =itk::BSplineTransform<CoordinateRepType, SpaceDimension, SplineOrder>;using OptimizerType = itk::LBFGS2Optimizerv4;using MetricType =itk::MeanSquaresImageToImageMetricv4<FixedImageType, MovingImageType>;using RegistrationType =itk::ImageRegistrationMethodv4<FixedImageType, MovingImageType>;MetricType::Pointer       metric = MetricType::New();OptimizerType::Pointer    optimizer = OptimizerType::New();RegistrationType::Pointer registration = RegistrationType::New();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]);FixedImageType::ConstPointer fixedImage = fixedImageReader->GetOutput();registration->SetFixedImage(fixedImage);registration->SetMovingImage(movingImageReader->GetOutput());fixedImageReader->Update();TransformType::Pointer outputBSplineTransform = TransformType::New();//我们transform转换对象，初始化其参数并将其连接到配准对象。using InitializerType =itk::BSplineTransformInitializer<TransformType, FixedImageType>;// Initialize the fixed parameters of transform (grid size, etc).InitializerType::Pointer transformInitializer = InitializerType::New();unsigned int numberOfGridNodesInOneDimension = 8;TransformType::MeshSizeType meshSize;meshSize.Fill(numberOfGridNodesInOneDimension - SplineOrder);transformInitializer->SetTransform(outputBSplineTransform);transformInitializer->SetImage(fixedImage);transformInitializer->SetTransformDomainMeshSize(meshSize);transformInitializer->InitializeTransform();using ParametersType = TransformType::ParametersType;// Set transform to identityconst unsigned int numberOfParameters =outputBSplineTransform->GetNumberOfParameters();ParametersType parameters(numberOfParameters);parameters.Fill(0.0);outputBSplineTransform->SetParameters(parameters);registration->SetInitialTransform(outputBSplineTransform);registration->InPlaceOn();constexpr unsigned int numberOfLevels = 3;//配准过程分为三个级别。为每个级别设置收缩因子和平滑sigma。RegistrationType::ShrinkFactorsArrayType shrinkFactorsPerLevel;shrinkFactorsPerLevel.SetSize(numberOfLevels);shrinkFactorsPerLevel[0] = 3;shrinkFactorsPerLevel[1] = 2;shrinkFactorsPerLevel[2] = 1;RegistrationType::SmoothingSigmasArrayType smoothingSigmasPerLevel;smoothingSigmasPerLevel.SetSize(numberOfLevels);smoothingSigmasPerLevel[0] = 2;smoothingSigmasPerLevel[1] = 1;smoothingSigmasPerLevel[2] = 0;registration->SetNumberOfLevels(numberOfLevels);registration->SetSmoothingSigmasPerLevel(smoothingSigmasPerLevel);registration->SetShrinkFactorsPerLevel(shrinkFactorsPerLevel);RegistrationType::TransformParametersAdaptorsContainerType adaptors;//创建变换适配器以针对此多分辨率方案的每个级别修改可变形变换的灵活性。TransformType::PhysicalDimensionsType fixedPhysicalDimensions;// First, get fixed image physical dimensionsfor (unsigned int i = 0; i < SpaceDimension; i++){fixedPhysicalDimensions[i] =fixedImage->GetSpacing()[i] *static_cast<double>(fixedImage->GetLargestPossibleRegion().GetSize()[i] - 1);}
// Create the transform adaptors specific to B-splinesfor (unsigned int level = 0; level < numberOfLevels; level++){using ShrinkFilterType =itk::ShrinkImageFilter<FixedImageType, FixedImageType>;ShrinkFilterType::Pointer shrinkFilter = ShrinkFilterType::New();shrinkFilter->SetShrinkFactors(shrinkFactorsPerLevel[level]);shrinkFilter->SetInput(fixedImage);shrinkFilter->Update();TransformType::MeshSizeType requiredMeshSize;// A good heuristic is to double the b-spline mesh resolution at each levelfor (unsigned int d = 0; d < ImageDimension; d++){requiredMeshSize[d] = meshSize[d] << level;}using BSplineAdaptorType =itk::BSplineTransformParametersAdaptor<TransformType>;BSplineAdaptorType::Pointer bsplineAdaptor = BSplineAdaptorType::New();bsplineAdaptor->SetTransform(outputBSplineTransform);bsplineAdaptor->SetRequiredTransformDomainMeshSize(requiredMeshSize);bsplineAdaptor->SetRequiredTransformDomainOrigin(shrinkFilter->GetOutput()->GetOrigin());bsplineAdaptor->SetRequiredTransformDomainDirection(shrinkFilter->GetOutput()->GetDirection());bsplineAdaptor->SetRequiredTransformDomainPhysicalDimensions(fixedPhysicalDimensions);adaptors.push_back(bsplineAdaptor);}registration->SetTransformParametersAdaptorsPerLevel(adaptors);using ScalesEstimatorType =itk::RegistrationParameterScalesFromPhysicalShift<MetricType>;ScalesEstimatorType::Pointer scalesEstimator = ScalesEstimatorType::New();scalesEstimator->SetMetric(metric);scalesEstimator->SetTransformForward(true);scalesEstimator->SetSmallParameterVariation(1.0);optimizer->SetScalesEstimator(scalesEstimator);optimizer->SetSolutionAccuracy(1e-4);optimizer->SetHessianApproximationAccuracy(5);optimizer->SetMaximumIterations(100);optimizer->SetMaximumLineSearchEvaluations(10);CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();optimizer->AddObserver(itk::IterationEvent(), observer);std::cout << "Starting Registration " << std::endl;try{registration->Update();std::cout << "Optimizer stop condition = "<< registration->GetOptimizer()->GetStopConditionDescription()<< std::endl;}catch (const itk::ExceptionObject & err){std::cerr << "ExceptionObject caught !" << std::endl;std::cerr << err << std::endl;return EXIT_FAILURE;}using ResampleFilterType =itk::ResampleImageFilter<MovingImageType, FixedImageType>;ResampleFilterType::Pointer resample = ResampleFilterType::New();resample->SetTransform(outputBSplineTransform);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>;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());try{writer->Update();}catch (const itk::ExceptionObject & err){std::cerr << "ExceptionObject caught !" << std::endl;std::cerr << err << std::endl;return EXIT_FAILURE;}using DifferenceFilterType =itk::SquaredDifferenceImageFilter<FixedImageType,FixedImageType,OutputImageType>;DifferenceFilterType::Pointer difference = DifferenceFilterType::New();WriterType::Pointer writer2 = WriterType::New();writer2->SetInput(difference->GetOutput());if (argc >= 5){difference->SetInput1(fixedImageReader->GetOutput());difference->SetInput2(resample->GetOutput());writer2->SetFileName(argv[4]);try{writer2->Update();}catch (const itk::ExceptionObject & err){std::cerr << "ExceptionObject caught !" << std::endl;std::cerr << err << std::endl;return EXIT_FAILURE;}}if (argc >= 6){writer2->SetFileName(argv[5]);difference->SetInput1(fixedImageReader->GetOutput());difference->SetInput2(movingImageReader->GetOutput());try{writer2->Update();}catch (const itk::ExceptionObject & err){std::cerr << "ExceptionObject caught !" << std::endl;std::cerr << err << std::endl;return EXIT_FAILURE;}}using VectorPixelType = itk::Vector<float, ImageDimension>;using DisplacementFieldImageType =itk::Image<VectorPixelType, ImageDimension>;using DisplacementFieldGeneratorType =itk::TransformToDisplacementFieldFilter<DisplacementFieldImageType,CoordinateRepType>;DisplacementFieldGeneratorType::Pointer dispfieldGenerator =DisplacementFieldGeneratorType::New();dispfieldGenerator->UseReferenceImageOn();dispfieldGenerator->SetReferenceImage(fixedImage);dispfieldGenerator->SetTransform(outputBSplineTransform);try{dispfieldGenerator->Update();}catch (const itk::ExceptionObject & err){std::cerr << "Exception detected while generating deformation field";std::cerr << " : " << err << std::endl;return EXIT_FAILURE;}using FieldWriterType = itk::ImageFileWriter<DisplacementFieldImageType>;FieldWriterType::Pointer fieldWriter = FieldWriterType::New();fieldWriter->SetInput(dispfieldGenerator->GetOutput());if (argc >= 7){fieldWriter->SetFileName(argv[6]);try{fieldWriter->Update();}catch (const itk::ExceptionObject & excp){std::cerr << "Exception thrown " << std::endl;std::cerr << excp << std::endl;return EXIT_FAILURE;}}return EXIT_SUCCESS;
}

3.13.5 BSplines Multi-Grid Image Registration in 3D

https://github.com/InsightSoftwareConsortium/ITK/blob/master/Examples/RegistrationITKv4/DeformableRegistration7.cxx

这个示例演示了如何使用itk::BSplineTransform类来执行两个3D图像的配准。示例代码在很大程度上与3.13.4节中给出的代码相同。主要区别在于，在本例中，我们将映像维设置为3，并将itk::LBFGSOptimizerv4优化器替换为itk::LBFGSBOptimizerv4。我们进行修改是因为我们发现LBFGS在初始位置处于或接近最佳时表现不佳;因此我们在无约束模式下使用LBFGSB。

BSplineTransform的参数空间由与BSpline网格节点相关的所有形变集合组成。大量的参数使它能够表示各种各样的形变，代价是需要大量的计算时间。

#include "itkImageRegistrationMethodv4.h"
#include "itkMeanSquaresImageToImageMetricv4.h"#include "itkBSplineTransform.h"//头文件
#include "itkLBFGSBOptimizerv4.h"#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"#include "itkResampleImageFilter.h"
#include "itkCastImageFilter.h"
#include "itkSquaredDifferenceImageFilter.h"#include "itkIdentityTransform.h"#include "itkBSplineTransformInitializer.h"
#include "itkTransformToDisplacementFieldFilter.h"#include "itkCommand.h"
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::LBFGSBOptimizerv4;using OptimizerPointer = const OptimizerType *;voidExecute(itk::Object * caller, const itk::EventObject & event) override{Execute((const itk::Object *)caller, event);}voidExecute(const itk::Object * object, const itk::EventObject & event) override{auto optimizer = static_cast<OptimizerPointer>(object);if (!(itk::IterationEvent().CheckEvent(&event))){return;}std::cout << optimizer->GetCurrentIteration() << "   ";std::cout << optimizer->GetCurrentMetricValue() << "   ";std::cout << optimizer->GetInfinityNormOfProjectedGradient() << std::endl;}
};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 << " [differenceOutputfile] [differenceBeforeRegistration] ";std::cerr << " [deformationField] ";return EXIT_FAILURE;}constexpr unsigned int ImageDimension = 3;using PixelType = float;using FixedImageType = itk::Image<PixelType, ImageDimension>;using MovingImageType = itk::Image<PixelType, ImageDimension>;const unsigned int     SpaceDimension = ImageDimension;//我们使用坐标表示的类型，空间的尺寸以及BSpline的顺序作为模板参数来实例化BSplineTransform的类型constexpr unsigned int SplineOrder = 3;using CoordinateRepType = double;using TransformType =itk::BSplineTransform<CoordinateRepType, SpaceDimension, SplineOrder>;using OptimizerType = itk::LBFGSBOptimizerv4;using MetricType =itk::MeanSquaresImageToImageMetricv4<FixedImageType, MovingImageType>;using RegistrationType =itk::ImageRegistrationMethodv4<FixedImageType, MovingImageType>;MetricType::Pointer       metric = MetricType::New();OptimizerType::Pointer    optimizer = OptimizerType::New();RegistrationType::Pointer registration = RegistrationType::New();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]);FixedImageType::ConstPointer fixedImage = fixedImageReader->GetOutput();registration->SetFixedImage(fixedImage);registration->SetMovingImage(movingImageReader->GetOutput());fixedImageReader->Update();TransformType::Pointer outputBSplineTransform = TransformType::New();// Initialize the transformusing InitializerType =itk::BSplineTransformInitializer<TransformType, FixedImageType>;InitializerType::Pointer transformInitializer = InitializerType::New();unsigned int numberOfGridNodesInOneDimension = 8;TransformType::MeshSizeType meshSize;meshSize.Fill(numberOfGridNodesInOneDimension - SplineOrder);transformInitializer->SetTransform(outputBSplineTransform);transformInitializer->SetImage(fixedImage);transformInitializer->SetTransformDomainMeshSize(meshSize);transformInitializer->InitializeTransform();// Set transform to identityusing ParametersType = TransformType::ParametersType;const unsigned int numberOfParameters =outputBSplineTransform->GetNumberOfParameters();ParametersType parameters(numberOfParameters);parameters.Fill(0.0);outputBSplineTransform->SetParameters(parameters);registration->SetInitialTransform(outputBSplineTransform);//构造，初始化对象并将其传递给配准方法。registration->InPlaceOn();//  A single level registration process is run using//  the shrink factor 1 and smoothing sigma 0.constexpr unsigned int numberOfLevels = 1;RegistrationType::ShrinkFactorsArrayType shrinkFactorsPerLevel;shrinkFactorsPerLevel.SetSize(numberOfLevels);shrinkFactorsPerLevel[0] = 1;RegistrationType::SmoothingSigmasArrayType smoothingSigmasPerLevel;smoothingSigmasPerLevel.SetSize(numberOfLevels);smoothingSigmasPerLevel[0] = 0;registration->SetNumberOfLevels(numberOfLevels);registration->SetSmoothingSigmasPerLevel(smoothingSigmasPerLevel);registration->SetShrinkFactorsPerLevel(shrinkFactorsPerLevel);const unsigned int numParameters =outputBSplineTransform->GetNumberOfParameters();//接下来，我们设置LBFGSB优化器的参数。注意，此优化器不支持scale estimator，并将所有参数设置为1。同样，我们应该为每个变量设置边界条件，其中boundSelect[i]可以设置为:UNBOUNDED、LOWERBOUNDED、BOTHBOUNDED、UPPERBOUNDED。OptimizerType::BoundSelectionType boundSelect(numParameters);OptimizerType::BoundValueType     upperBound(numParameters);OptimizerType::BoundValueType     lowerBound(numParameters);boundSelect.Fill(OptimizerType::UNBOUNDED);upperBound.Fill(0.0);lowerBound.Fill(0.0);optimizer->SetBoundSelection(boundSelect);optimizer->SetUpperBound(upperBound);optimizer->SetLowerBound(lowerBound);optimizer->SetCostFunctionConvergenceFactor(1e+12);optimizer->SetGradientConvergenceTolerance(1.0e-35);optimizer->SetNumberOfIterations(500);optimizer->SetMaximumNumberOfFunctionEvaluations(500);optimizer->SetMaximumNumberOfCorrections(5);CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();optimizer->AddObserver(itk::IterationEvent(), observer);std::cout << "Starting Registration " << std::endl;try{registration->Update();std::cout << "Optimizer stop condition = "<< registration->GetOptimizer()->GetStopConditionDescription()<< std::endl;}catch (const itk::ExceptionObject & err){std::cerr << "ExceptionObject caught !" << std::endl;std::cerr << err << std::endl;return EXIT_FAILURE;}OptimizerType::ParametersType finalParameters =outputBSplineTransform->GetParameters();std::cout << "Last Transform Parameters" << std::endl;std::cout << finalParameters << std::endl;// Finally we use the last transform in order to resample the image.using ResampleFilterType =itk::ResampleImageFilter<MovingImageType, FixedImageType>;ResampleFilterType::Pointer resample = ResampleFilterType::New();resample->SetTransform(outputBSplineTransform);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>;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());try{writer->Update();}catch (const itk::ExceptionObject & err){std::cerr << "ExceptionObject caught !" << std::endl;std::cerr << err << std::endl;return EXIT_FAILURE;}using DifferenceFilterType =itk::SquaredDifferenceImageFilter<FixedImageType,FixedImageType,OutputImageType>;DifferenceFilterType::Pointer difference = DifferenceFilterType::New();WriterType::Pointer writer2 = WriterType::New();writer2->SetInput(difference->GetOutput());// Compute the difference image between the fixed and resampled moving image.if (argc >= 5){difference->SetInput1(fixedImageReader->GetOutput());difference->SetInput2(resample->GetOutput());writer2->SetFileName(argv[4]);try{writer2->Update();}catch (const itk::ExceptionObject & err){std::cerr << "ExceptionObject caught !" << std::endl;std::cerr << err << std::endl;return EXIT_FAILURE;}}// Compute the difference image between the fixed and moving image before registration.if (argc >= 6){writer2->SetFileName(argv[5]);difference->SetInput1(fixedImageReader->GetOutput());difference->SetInput2(movingImageReader->GetOutput());try{writer2->Update();}catch (const itk::ExceptionObject & err){std::cerr << "ExceptionObject caught !" << std::endl;std::cerr << err << std::endl;return EXIT_FAILURE;}}// Generate the explicit deformation field resulting from the registration.using VectorPixelType = itk::Vector<float, ImageDimension>;using DisplacementFieldImageType =itk::Image<VectorPixelType, ImageDimension>;using DisplacementFieldGeneratorType =itk::TransformToDisplacementFieldFilter<DisplacementFieldImageType,CoordinateRepType>;/** Create an setup displacement field generator. */DisplacementFieldGeneratorType::Pointer dispfieldGenerator =DisplacementFieldGeneratorType::New();dispfieldGenerator->UseReferenceImageOn();dispfieldGenerator->SetReferenceImage(fixedImage);dispfieldGenerator->SetTransform(outputBSplineTransform);try{dispfieldGenerator->Update();}catch (const itk::ExceptionObject & err){std::cerr << "Exception detected while generating deformation field";std::cerr << " : " << err << std::endl;return EXIT_FAILURE;}using FieldWriterType = itk::ImageFileWriter<DisplacementFieldImageType>;FieldWriterType::Pointer fieldWriter = FieldWriterType::New();fieldWriter->SetInput(dispfieldGenerator->GetOutput());if (argc >= 7){fieldWriter->SetFileName(argv[6]);try{fieldWriter->Update();}catch (const itk::ExceptionObject & excp){std::cerr << "Exception thrown " << std::endl;std::cerr << excp << std::endl;return EXIT_FAILURE;}}return EXIT_SUCCESS;
}

3.13.6 Image Warping with Kernel Splines

https://github.com/InsightSoftwareConsortium/ITK/blob/master/Examples/RegistrationITKv4/LandmarkWarping2.cxx

这个例子演示了如何使用一个KernelBase样条和两组landmarks来deform图像。

#include "itkImage.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkDisplacementFieldTransform.h"
#include "itkResampleImageFilter.h"
#include "itkVector.h"//包含头文件
#include "itkLandmarkDisplacementFieldSource.h"
#include <fstream>int
main(int argc, char * argv[])
{if (argc < 3){std::cerr << "Missing Parameters " << std::endl;std::cerr << "Usage: " << argv[0];std::cerr << " landmarksFile fixedImage ";std::cerr << "movingImage deformedMovingImage" << std::endl;return EXIT_FAILURE;}constexpr unsigned int Dimension = 2;using VectorComponentType = float;using VectorType = itk::Vector<VectorComponentType, Dimension>;using DisplacementFieldType = itk::Image<VectorType, Dimension>;using PixelType = unsigned char;using FixedImageType = itk::Image<PixelType, Dimension>;using MovingImageType = itk::Image<PixelType, Dimension>;using FixedReaderType = itk::ImageFileReader<FixedImageType>;using MovingReaderType = itk::ImageFileReader<MovingImageType>;using MovingWriterType = itk::ImageFileWriter<MovingImageType>;FixedReaderType::Pointer fixedReader = FixedReaderType::New();fixedReader->SetFileName(argv[2]);try{fixedReader->Update();}catch (const itk::ExceptionObject & excp){std::cerr << "Exception thrown " << std::endl;std::cerr << excp << std::endl;return EXIT_FAILURE;}MovingReaderType::Pointer movingReader = MovingReaderType::New();MovingWriterType::Pointer movingWriter = MovingWriterType::New();movingReader->SetFileName(argv[3]);movingWriter->SetFileName(argv[4]);FixedImageType::ConstPointer fixedImage = fixedReader->GetOutput();using DisplacementSourceType =itk::LandmarkDisplacementFieldSource<DisplacementFieldType>;//读取fixed图像和moving图像后，将从itk::LandmarkDisplacementFieldSource类实例化变形器对象，并设置图像空间和方向的参数。DisplacementSourceType::Pointer deformer = DisplacementSourceType::New();deformer->SetOutputSpacing(fixedImage->GetSpacing());deformer->SetOutputOrigin(fixedImage->GetOrigin());deformer->SetOutputRegion(fixedImage->GetLargestPossibleRegion());deformer->SetOutputDirection(fixedImage->GetDirection());using LandmarkContainerType = DisplacementSourceType::LandmarkContainer;//然后创建Source和target landmarks，并从文件流中读取点本身。using LandmarkPointType = DisplacementSourceType::LandmarkPointType;LandmarkContainerType::Pointer sourceLandmarks =LandmarkContainerType::New();LandmarkContainerType::Pointer targetLandmarks =LandmarkContainerType::New();LandmarkPointType sourcePoint;LandmarkPointType targetPoint;std::ifstream pointsFile;pointsFile.open(argv[1]);unsigned int pointId = 0;pointsFile >> sourcePoint;pointsFile >> targetPoint;while (!pointsFile.fail()){sourceLandmarks->InsertElement(pointId, sourcePoint);targetLandmarks->InsertElement(pointId, targetPoint);++pointId;pointsFile >> sourcePoint;pointsFile >> targetPoint;}pointsFile.close();deformer->SetSourceLandmarks(sourceLandmarks);//然后将source地标和target地标对象分配给变形器。deformer->SetTargetLandmarks(targetLandmarks);try{deformer->UpdateLargestPossibleRegion();//在deformer上调用UpdateLargestPossibleRegion（）之后，可以通过GetOutput（）方法获取displacement field。}catch (const itk::ExceptionObject & excp){std::cerr << "Exception thrown " << std::endl;std::cerr << excp << std::endl;return EXIT_FAILURE;}DisplacementFieldType::Pointer displacementField = deformer->GetOutput();using InterpolatorPrecisionType = double;using TransformPrecisionType = float;using FilterType = itk::ResampleImageFilter<MovingImageType,MovingImageType,InterpolatorPrecisionType,TransformPrecisionType>;FilterType::Pointer warper = FilterType::New();using InterpolatorType =itk::LinearInterpolateImageFunction<MovingImageType,InterpolatorPrecisionType>;InterpolatorType::Pointer interpolator = InterpolatorType::New();warper->SetInterpolator(interpolator);using DisplacementFieldTransformType =itk::DisplacementFieldTransform<TransformPrecisionType, Dimension>;auto displacementTransform = DisplacementFieldTransformType::New();displacementTransform->SetDisplacementField(displacementField);warper->SetTransform(displacementTransform);warper->SetOutputSpacing(displacementField->GetSpacing());warper->SetOutputOrigin(displacementField->GetOrigin());warper->SetOutputDirection(displacementField->GetDirection());warper->SetSize(displacementField->GetLargestPossibleRegion().GetSize());warper->SetInput(movingReader->GetOutput());movingWriter->SetInput(warper->GetOutput());try{movingWriter->Update();}catch (const itk::ExceptionObject & excp){std::cerr << "Exception thrown " << std::endl;std::cerr << excp << std::endl;return EXIT_FAILURE;}return EXIT_SUCCESS;
}

3.13.7 Image Warping with BSplines

https://github.com/InsightSoftwareConsortium/ITK/blob/master/Examples/RegistrationITKv4/BSplineWarping1.cxx

本示例说明了如何使用BSplineTransform变形2D图像。

b样条网格现在应该在每个节点上提供coefficients。由于这是一个二维网格，每个节点应该接收两个coefficients。每个coefficient对代表这个节点上的一个位移向量。coefficients可以以数组的形式传递给b-spline，其中第一组元素是所有节点位移的第一个分量，第二组元素由所有节点位移的第二个分量组成。

在这个例子中，我们从一个文件中读取这样的displacements，但是为了方便起见，我们使用每个节点的(x,y)位移对来编写这个文件。因此，当将从文件中读取的元素分配给数组的元素时，应该重新组织这些元素。我们通过将文件中所有奇数元素存储在数组的第一个块中，将文件中所有偶数元素存储在数组的第二个块中来实现这一点。最后，使用SetParameters()方法将数组传递给b -spline变换。

#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"#include "itkResampleImageFilter.h"#include "itkBSplineTransform.h"
#include "itkTransformFileWriter.h"#include <fstream>//  The following section of code implements a Command observer used to monitor the evolution of the registration process.#include "itkCommand.h"
class CommandProgressUpdate : public itk::Command
{
public:using Self = CommandProgressUpdate;using Superclass = itk::Command;using Pointer = itk::SmartPointer<Self>;itkNewMacro(Self);protected:CommandProgressUpdate() = default;public:voidExecute(itk::Object * caller, const itk::EventObject & event) override{Execute((const itk::Object *)caller, event);}voidExecute(const itk::Object * object, const itk::EventObject & event) override{const auto * filter = static_cast<const itk::ProcessObject *>(object);if (!itk::ProgressEvent().CheckEvent(&event)){return;}std::cout << filter->GetProgress() << std::endl;}
};int
main(int argc, char * argv[])
{if (argc < 5){std::cerr << "Missing Parameters " << std::endl;std::cerr << "Usage: " << argv[0];std::cerr << " coefficientsFile fixedImage ";std::cerr << "movingImage deformedMovingImage" << std::endl;std::cerr << "[deformationField] [transformFile]" << std::endl;return EXIT_FAILURE;}constexpr unsigned int ImageDimension = 2;//首先，我们为fixed和moving图像以及图像读取器定义必要的类型。using PixelType = unsigned char;using FixedImageType = itk::Image<PixelType, ImageDimension>;using MovingImageType = itk::Image<PixelType, ImageDimension>;using FixedReaderType = itk::ImageFileReader<FixedImageType>;using MovingReaderType = itk::ImageFileReader<MovingImageType>;using MovingWriterType = itk::ImageFileWriter<MovingImageType>;FixedReaderType::Pointer fixedReader = FixedReaderType::New();fixedReader->SetFileName(argv[2]);try{fixedReader->Update();}catch (const itk::ExceptionObject & excp){std::cerr << "Exception thrown " << std::endl;std::cerr << excp << std::endl;return EXIT_FAILURE;}MovingReaderType::Pointer movingReader = MovingReaderType::New();MovingWriterType::Pointer movingWriter = MovingWriterType::New();movingReader->SetFileName(argv[3]);movingWriter->SetFileName(argv[4]);FixedImageType::ConstPointer fixedImage = fixedReader->GetOutput();using FilterType =itk::ResampleImageFilter<MovingImageType, FixedImageType>;FilterType::Pointer resampler = FilterType::New();using InterpolatorType =itk::LinearInterpolateImageFunction<MovingImageType, double>;InterpolatorType::Pointer interpolator = InterpolatorType::New();resampler->SetInterpolator(interpolator);FixedImageType::SpacingType   fixedSpacing = fixedImage->GetSpacing();//使用fixed图像中的值在resampler中设置相应的值。FixedImageType::PointType     fixedOrigin = fixedImage->GetOrigin();FixedImageType::DirectionType fixedDirection = fixedImage->GetDirection();resampler->SetOutputSpacing(fixedSpacing);resampler->SetOutputOrigin(fixedOrigin);resampler->SetOutputDirection(fixedDirection);FixedImageType::RegionType fixedRegion = fixedImage->GetBufferedRegion();FixedImageType::SizeType   fixedSize = fixedRegion.GetSize();resampler->SetSize(fixedSize);resampler->SetOutputStartIndex(fixedRegion.GetIndex());resampler->SetInput(movingReader->GetOutput());movingWriter->SetInput(resampler->GetOutput());const unsigned int     SpaceDimension = ImageDimension;//现在，我们使用坐标表示的类型，空间的尺寸以及B-spline的顺序作为模板参数，实例化BSplineTransform的类型。constexpr unsigned int SplineOrder = 3;using CoordinateRepType = double;using TransformType =itk::BSplineTransform<CoordinateRepType, SpaceDimension, SplineOrder>;TransformType::Pointer bsplineTransform = TransformType::New();constexpr unsigned int numberOfGridNodes = 7;//使用来自fixed图像和mesh的值填充B-spline变换的参数。TransformType::PhysicalDimensionsType fixedPhysicalDimensions;TransformType::MeshSizeType           meshSize;for (unsigned int i = 0; i < SpaceDimension; i++){fixedPhysicalDimensions[i] =fixedSpacing[i] * static_cast<double>(fixedSize[i] - 1);}meshSize.Fill(numberOfGridNodes - SplineOrder);bsplineTransform->SetTransformDomainOrigin(fixedOrigin);bsplineTransform->SetTransformDomainPhysicalDimensions(fixedPhysicalDimensions);bsplineTransform->SetTransformDomainMeshSize(meshSize);bsplineTransform->SetTransformDomainDirection(fixedDirection);using ParametersType = TransformType::ParametersType;const unsigned int numberOfParameters =bsplineTransform->GetNumberOfParameters();const unsigned int numberOfNodes = numberOfParameters / SpaceDimension;ParametersType parameters(numberOfParameters);std::ifstream infile;infile.open(argv[1]);for (unsigned int n = 0; n < numberOfNodes; ++n){infile >> parameters[n];infile >> parameters[n + numberOfNodes];}infile.close();bsplineTransform->SetParameters(parameters);//最后，使用SetParameters()将数组传递给B-Spline变换。CommandProgressUpdate::Pointer observer = CommandProgressUpdate::New();resampler->AddObserver(itk::ProgressEvent(), observer);resampler->SetTransform(bsplineTransform);//至此，我们准备将transform用作resample过滤器的一部分。我们通过在管道的最后一个过滤器（在本例中为writer）上调用Update（）来触发管道的执行。try{movingWriter->Update();}catch (const itk::ExceptionObject & excp){std::cerr << "Exception thrown " << std::endl;std::cerr << excp << std::endl;return EXIT_FAILURE;}using VectorType = itk::Vector<float, ImageDimension>;using DisplacementFieldType = itk::Image<VectorType, ImageDimension>;DisplacementFieldType::Pointer field = DisplacementFieldType::New();field->SetRegions(fixedRegion);field->SetOrigin(fixedOrigin);field->SetSpacing(fixedSpacing);field->SetDirection(fixedDirection);field->Allocate();using FieldIterator = itk::ImageRegionIterator<DisplacementFieldType>;FieldIterator fi(field, fixedRegion);fi.GoToBegin();TransformType::InputPointType    fixedPoint;TransformType::OutputPointType   movingPoint;DisplacementFieldType::IndexType index;VectorType displacement;while (!fi.IsAtEnd()){index = fi.GetIndex();field->TransformIndexToPhysicalPoint(index, fixedPoint);movingPoint = bsplineTransform->TransformPoint(fixedPoint);displacement = movingPoint - fixedPoint;fi.Set(displacement);++fi;}using FieldWriterType = itk::ImageFileWriter<DisplacementFieldType>;FieldWriterType::Pointer fieldWriter = FieldWriterType::New();fieldWriter->SetInput(field);if (argc >= 6){fieldWriter->SetFileName(argv[5]);try{fieldWriter->Update();}catch (const itk::ExceptionObject & excp){std::cerr << "Exception thrown " << std::endl;std::cerr << excp << std::endl;return EXIT_FAILURE;}}if (argc >= 7){fieldWriter->SetFileName(argv[6]);try{using TransformWriterType = itk::TransformFileWriter;TransformWriterType::Pointer transformWriter =TransformWriterType::New();transformWriter->AddTransform(bsplineTransform);transformWriter->SetFileName(argv[6]);transformWriter->Update();}catch (const itk::ExceptionObject & excp){std::cerr << "Exception thrown " << std::endl;std::cerr << excp << std::endl;return EXIT_FAILURE;}}return EXIT_SUCCESS;
}

itk registration 13相关推荐

itk registration 14
参考:https://itk.org/ITKSoftwareGuide/html/Book2/ITKSoftwareGuide-Book2ch3.html#x26-1740003.18 3.14 De ...
itk registration 9
参考:https://itk.org/ITKSoftwareGuide/html/Book2/ITKSoftwareGuide-Book2ch3.html#x26-1740003.18 3.9 Tra ...
手把手教你Linux下基于ITK开发基础（1）
手把手教你Linux下基于ITK开发基础(1) *requirement: requirement version git 2.17.1 g++ Ubuntu 7.5.0-3 ubuntu1~18.0 ...
Mongo北京大会3月3号召开！报名抢注火爆进行中！（免费）
Mongo Beijing 简介 Mongo北京开发者聚会将讨论如何用开源的非关系型数据库MongoDB构建多彩的应用.由MongoDB工程师和经验丰富的MongoDB用户做演讲.演讲分为中文和英文, ...
龙芯俱乐部邀请参加“2008开源中国开源世界”高峰论坛
<script type="text/javascript" src="http://www.loongsonclub.cn/bbs/include/javascr ...
mahout安装测试
Mahout 是 Apache Software Foundation(ASF) 旗下的一个开源项目,提供一些可扩展的机器学习领域经典算法的实现,旨在帮助开发人员更加方便快捷地创建智能应用程序.Apa ...
ITK:Watch注册Registration
ITK:Watch注册Registration 内容提要输出结果 C++实现代码内容提要观看使用VTK进行的注册迭代输出结果输出: Optimizer stop condition: Reg ...
VS2019(x64)+ITK4.13.2手动新建项目配置ITK(不使用CMakeLists.txt构建)
完整博客见:Win10-64位(或Win7)+VS2019+ITK4.13.2+VTK8.2+CMake3.16环境配置目录: 3.6 手动新建项目配置ITK(不使用CMakeLists.txt构 ...
ITK入门教程（13）点集之得到点中的存储数据
目录 1.概述 2.过程讲解 3.代码 4.结果 1.概述 Itk::PointSet 类是与图像类交互的.因此点集中的点可以很方便地保存从图像计算得到的值.与点相关的值被指定为PixelType像素 ...

itk registration 13

3.13 Deformable Registration

itk registration 13相关推荐

最新文章

热门文章