itk registration 13
参考：https://itk.org/ITKSoftwareGuide/html/Book2/ITKSoftwareGuide-Book2ch3.html#x26-1740003.183.13Deformable Registration3.13.1FEM-Based Image Registrationhttps://github.com/InsightSoftwareConsortium/IT
·
参考：https://itk.org/ITKSoftwareGuide/html/Book2/ITKSoftwareGuide-Book2ch3.html#x26-1740003.18
3.13 Deformable Registration
3.13.1 FEM-Based Image Registration
#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 properties
m->SetGlobalNumber(0);
m->SetYoungsModulus(registrationFilter->GetElasticity());// Young's modulus of the membrane
m->SetCrossSectionalArea(1.0); // Cross-sectional area
m->SetThickness(1.0); // Thickness
m->SetMomentOfInertia(1.0); // Moment of inertia
m->SetPoissonsRatio(0.); // Poisson's ratio -- DONT CHOOSE 1.0!!
m->SetDensityHeatProduct(1.0); // Density-Heat capacity product
ElementType::Pointer e1 = ElementType::New();// Create the element type
e1->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图像配准
这个例子演示了如何使用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
#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:
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
{
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
这个示例演示了在多分辨率方案中使用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 *;
void
Execute(itk::Object * caller, const itk::EventObject & event) override
{
Execute((const itk::Object *)caller, event);
}
void
Execute(const itk::Object * object, const itk::EventObject & event) override
{
auto optimizer = static_cast<OptimizerPointer>(object);
if (!itk::IterationEvent().CheckEvent(&event))
{
return;
}
std::cout << 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 identity
const 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 dimensions
for (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-splines
for (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 level
for (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
这个示例演示了如何使用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 *;
void
Execute(itk::Object * caller, const itk::EventObject & event) override
{
Execute((const itk::Object *)caller, event);
}
void
Execute(const itk::Object * object, const itk::EventObject & event) override
{
auto optimizer = static_cast<OptimizerPointer>(object);
if (!(itk::IterationEvent().CheckEvent(&event)))
{
return;
}
std::cout << 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 transform
using 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 identity
using 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
这个例子演示了如何使用一个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
本示例说明了如何使用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:
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
{
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;
}
开放原子开发者工作坊旨在鼓励更多人参与开源活动，与志同道合的开发者们相互交流开发经验、分享开发心得、获取前沿技术趋势。工作坊有多种形式的开发者活动，如meetup、训练营等，主打技术交流，干货满满，真诚地邀请各位开发者共同参与！
更多推荐
1
0- 0
已为社区贡献3条内容
所有评论(0)