参考：https://itk.org/ITKSoftwareGuide/html/Book2/ITKSoftwareGuide-Book2ch3.html#x26-1740003.18
3.14 Demons Deformable Registration
针对单模可变形配准的问题，Insight Toolkit 提供了Thirion s demons算法的实现。在这个实现中，每个图像被视为一组等强度轮廓。其主要思想是，一个regular grid of forces deform通过将轮廓推向法向而使图像变形。由optical flow方程导出位移的方向和大小：D(X) · ∇f(X) = -(m(X) - f(X))
上式中，f(X)为fixed图像，m(X)为待配准的moving图像，D(X)为图像间的displacemnt或optical flow。在optical flow文献中，众所周知，公式3.33不足以局部指定D(X)，通常使用某种形式的正则化来确定。配准时，使用矢量在强度梯度方向上的投影：
但是，该方程对于较小的图像梯度值变得不稳定，从而导致较大的位移值。为了克服这个问题，Thirion重新对方程进行了归一化：
其中K是一个归一化因子，解释了强度和梯度之间的单位不平衡。这个因子被计算为像素间隔的平方平均值。K的加入保证了力的计算不受图像像素尺度的影响。
从初始变形场开始，demons算法使用公式3.35迭代更新field，以使第N次迭代的场由下式给出：
deformation的重建是一个病态问题，其中fixed和moving图像匹配有很多解决方案。例如，由于每个图像像素都可以自由地独立移动，所以m(X)中某个特定值的所有像素都有可能映射到f(X)中相同值的单个图像像素。由此产生的变形场对于实际应用可能是不现实的。解决该场的唯一选择是强制一个类似弹性elastic-like的行为，在迭代之间使用高斯滤波器平滑变形场。
在ITK中，demons算法被实现为有限差分求解器（FDS）框架的一部分，并在以下示例中演示了其用法。
3.14.1 Asymmetrical Demons Deformable Registration
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkImageRegionIterator.h"
#include "itkDemonsRegistrationFilter.h"//此示例演示了如何使用“demo”算法可变形地配准两个图像。第一步是包括头文件。
#include "itkHistogramMatchingImageFilter.h"
#include "itkCastImageFilter.h"
#include "itkResampleImageFilter.h"
#include "itkDisplacementFieldTransform.h"
// The following section of code implements a Command observer
// that will monitor the evolution of the registration process.
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::DemonsRegistrationFilter<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 (!(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>;
// Set up the file readers
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>;//demons算法依赖于这样的假设，即在一个物体上代表相同同源点的像素对待配准的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 :: DemonsRegistrationFilter中，deformation field表示为图像，其像素为浮点向量。
using DisplacementFieldType = itk::Image<VectorPixelType, Dimension>;
using RegistrationFilterType =
itk::DemonsRegistrationFilter<InternalImageType,
InternalImageType,
DisplacementFieldType>;
RegistrationFilterType::Pointer filter = RegistrationFilterType::New();
// Create the Command observer and register it with the registration filter.
CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
filter->AddObserver(itk::IterationEvent(), observer);
filter->SetFixedImage(fixedImageCaster->GetOutput());//输入的fixed图像只是fixed图像caseting滤波器的输出。输入的moving图像是直方图匹配滤波器的输出。
filter->SetMovingImage(matcher->GetOutput());
filter->SetNumberOfIterations(50);//demos配准过滤器具有两个参数：每次迭代后要执行的迭代次数和要应用于变形场的高斯平滑核的标准偏差。
filter->SetStandardDeviations(1.0);
filter->Update();//通过更新过滤器触发配准算法。过滤器输出是计算出的变形场。
using OutputPixelType = unsigned char;//itk::ResampleImageFilter可以用于使moving图像与输出deformation field一起变形。使用了itk::ResampleImageFilter的默认插值器，但需要指定输出图像的spacing和origin。
using OutputImageType = itk::Image<OutputPixelType, Dimension>;
using InterpolatorPrecisionType = double;
using WarperType = itk::ResampleImageFilter<MovingImageType,
OutputImageType,
InterpolatorPrecisionType,
float>;
WarperType::Pointer warper = WarperType::New();
warper->SetInput(movingImageReader->GetOutput());
warper->UseReferenceImageOn();
warper->SetReferenceImage(fixedImageReader->GetOutput());
using DisplacementFieldTransformType =
itk::DisplacementFieldTransform<InternalPixelType, Dimension>;//ResampleImageFilter需要变换，因此必须构造一个itk::DisplacementFieldTransform，然后将其设置为ResampleImageFilter的变换。按照先前的示例，将生成的变形或重新采样的图像写入文件。
auto displacementTransform = DisplacementFieldTransformType::New();
displacementTransform->SetDisplacementField(filter->GetOutput());
warper->SetTransform(displacementTransform);
// Write warped image out to file
using WriterType = itk::ImageFileWriter<OutputImageType>;
WriterType::Pointer writer = WriterType::New();
writer->SetFileName(argv[3]);
writer->SetInput(warper->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);
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;
}此外，还实现了deformed moving图像的梯度计算。这在PDE更新的一次迭代中提供了一定程度的force计算的对称性。本例中使用的方程为：
3.14.2 Symmetrical Demons Deformable Registration
这个例子演示了如何使用一个的demons算法变体来形变配准两幅图像。这个变体使用不同的公式来计算用于图像的forces，以便计算变形场。该方法利用fixed图像的梯度和deformed moving图像的梯度来计算forces。这种计算力的机制引入了相对于fixed和moving图像的选择的对称性。这种对称性只在PDE更新的一次迭代的计算过程中存在。这种机制不太可能在整个配准过程中实现完全对称。
demons算法依赖于这样的假设，即在一个物体上代表相同同源点的像素对待配准的fixed和moving图像具有相同的强度。在本例中，我们将对moving图像进行预处理，使用itk::HistogramMatchingImageFilter来匹配图像之间的强度。
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkSymmetricForcesDemonsRegistrationFilter.h"//头文件
#include "itkHistogramMatchingImageFilter.h"
#include "itkCastImageFilter.h"
#include "itkDisplacementFieldTransform.h"
#include "itkResampleImageFilter.h"
// The following section of code implements a Command observer
// that will monitor the evolution of the registration process.
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::SymmetricForcesDemonsRegistrationFilter<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 (!(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;
return EXIT_FAILURE;
}
constexpr unsigned int Dimension = 2;//其次，我们声明图像的类型。
using PixelType = unsigned short;
using FixedImageType = itk::Image<PixelType, Dimension>;
using MovingImageType = itk::Image<PixelType, Dimension>;
// Set up the file readers
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图像强度的重映射，我们声明了一个内部图像类型，该内部图像类型使用浮点像素类型，并将输入图像转换为内部图像类型。
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>;//基本思想是按照用户指定的分位数值数目匹配两个图像的直方图。为了鲁棒性，对直方图进行匹配，从而将背景像素从两个直方图中排除。对于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();//我们选择要表示直方图的bin数量以及要匹配直方图的点或分位数的数量。
using VectorPixelType = itk::Vector<float, Dimension>;//在itk::SymmetricForcesDemonsRegistrationFilter中，deformation field表示为像素为浮点向量的图像。
using DisplacementFieldType = itk::Image<VectorPixelType, Dimension>;
using RegistrationFilterType =
itk::SymmetricForcesDemonsRegistrationFilter<InternalImageType,
InternalImageType,
DisplacementFieldType>;
RegistrationFilterType::Pointer filter = RegistrationFilterType::New();
// Create the Command observer and register it with the registration filter.
CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
filter->AddObserver(itk::IterationEvent(), observer);
filter->SetFixedImage(fixedImageCaster->GetOutput());//输入的fixed图像只是fxied图像casting滤波器的输出。输入的moving图像是直方图匹配滤波器的输出
filter->SetMovingImage(matcher->GetOutput());
filter->SetNumberOfIterations(50);//demons配准过滤器具有两个参数：每次迭代后要执行的迭代次数和要应用于变形场的高斯平滑核的标准偏差。
filter->SetStandardDeviations(1.0);
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可以用输出变形场来warp moving图像。像itk::ResampleImageFilter一样，WarpImageFilter要求对输入图像进行重新采样，输入图像插值器，以及输出图像的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相对于矢量图像表示的变形场扭曲或变换输入图像。按照先前的示例，将生成的变形或重新采样的图像写入文件。 warper->SetDisplacementField( filter->GetOutput() );
// Write warped image out to file
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();
}
return EXIT_SUCCESS;
}
55
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
