itk registration 5

最新推荐文章于 2024-05-08 15:43:05 发布

shchojj

最新推荐文章于 2024-05-08 15:43:05 发布

阅读量279

点赞数

分类专栏： rigid registration

原文链接：https://itk.org/ITKSoftwareGuide/html/Book2/ITKSoftwareGuide-Book2ch3.html#x26-1740003.18

版权

rigid registration 专栏收录该内容

22 篇文章 8 订阅

订阅专栏

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

3.5多模配准

基于互信息评价的度量方法能够很好地克服多模配准的困难。

3.5.1 Mattes Mutual Information

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

//包含头文件
#include "itkImageRegistrationMethodv4.h"
#include "itkTranslationTransform.h"
#include "itkMattesMutualInformationImageToImageMetricv4.h"
#include "itkRegularStepGradientDescentOptimizerv4.h"

#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"

#include "itkResampleImageFilter.h"
#include "itkCastImageFilter.h"
#include "itkCheckerBoardImageFilter.h"
#include "itkMersenneTwisterRandomVariateGenerator.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::RegularStepGradientDescentOptimizerv4<double>;
  using OptimizerPointer = const OptimizerType *;

  void
  Execute(itk::Object * caller, const itk::EventObject & event) override
  {
    Execute((const itk::Object *)caller, event);
  }

  void
  Execute(const itk::Object * object, const itk::EventObject & event) override
  {
    auto optimizer = static_cast<OptimizerPointer>(object);
    if (!itk::IterationEvent().CheckEvent(&event))
    {
      return;
    }
    std::cout << optimizer->GetCurrentIteration() << "   ";
    std::cout << optimizer->GetValue() << "   ";
    std::cout << optimizer->GetCurrentPosition() << 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 [defaultPixelValue]" << std::endl;
    std::cerr << "[checkerBoardAfter] [checkerBoardBefore]" << std::endl;
    std::cerr << "[numberOfBins] [numberOfSamples]";
    std::cerr << "[useExplicitPDFderivatives ] " << std::endl;
    return EXIT_FAILURE;
  }

  constexpr unsigned int Dimension = 2;
  using PixelType = float;

  using FixedImageType = itk::Image<PixelType, Dimension>;
  using MovingImageType = itk::Image<PixelType, Dimension>;

  using TransformType = itk::TranslationTransform<double, Dimension>;
  using OptimizerType = itk::RegularStepGradientDescentOptimizerv4<double>;
  using RegistrationType = itk::
    ImageRegistrationMethodv4<FixedImageType, MovingImageType, TransformType>;

  using MetricType =
    itk::MattesMutualInformationImageToImageMetricv4<FixedImageType,
                                                     MovingImageType>;

  OptimizerType::Pointer    optimizer = OptimizerType::New();
  RegistrationType::Pointer registration = RegistrationType::New();

  registration->SetOptimizer(optimizer);

  MetricType::Pointer metric = MetricType::New();//创建度量方法
  registration->SetMetric(metric);//连接到配准对象

  unsigned int numberOfBins = 24;

  if (argc > 7)
  {
    numberOfBins = std::stoi(argv[7]);
  }


  metric->SetNumberOfHistogramBins(numberOfBins);
 //为了计算图像梯度，使用了基于ImageFunction的图像梯度计算器来代替图像梯度滤波器。图像渐变方法在超类ImageToImageMetricv4中定义
  metric->SetUseMovingImageGradientFilter(false);
  metric->SetUseFixedImageGradientFilter(false);

  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]);

  registration->SetFixedImage(fixedImageReader->GetOutput());
  registration->SetMovingImage(movingImageReader->GetOutput());
 
//注意，较大的学习率值将使优化器不稳定。另一方面，小的值可能导致优化器需要太多的迭代来遍历成本函数的极值。优化这个参数的简单方法是从很小的值开始，可能在{1.0,5.0}的范围内。一旦调整了其他配准参数以产生收敛性，您可能希望重新访问学习率并开始增加它的值，直到您观察到优化变得不稳定为止。该参数的理想值是在优化的参数空间上仍然保持稳定路径的情况下迭代次数最少的参数。记住，这个参数是一个应用在度量梯度上的乘法因子。因此，它对优化器步长的影响与度量值本身成比例。具有较大值的度量将要求您为学习率使用较小的值，以便维护类似的优化器行为。
  optimizer->SetLearningRate(8.00);//请注意，在ITKv4注册框架中，优化器总是试图最小化代价函数，而指标总是返回一个参数和改进优化的导数结果，因此这个指标计算负的相互信息。
  optimizer->SetMinimumStepLength(0.001);
  optimizer->SetNumberOfIterations(200);
  optimizer->ReturnBestParametersAndValueOn();
 //每当regular step gradient descent optimizer在参数空间的移动方向发生变化时，它都会减小步长的大小。减小步长的速率由松弛因子控制。该因子的默认值是0.5。然而，这个值对于有噪声的度量可能是不够的，因为它们往往会在优化器上引起不稳定的移动，从而导致许多方向的变化。在这些条件下，优化器将迅速缩小步长，但它仍然离代价函数的极值位置太远。在本例中，我们将松弛因子设置为一个高于默认值的数字，以防止步长过早收缩。
  optimizer->SetRelaxationFactor(0.8);
//与其使用整个虚拟域(通常是固定的图像域)进行配准，我们可以使用一个空间采样点集，通过提供一个任意的点列表来评估metric。期望点列表在固定的图像域内，根据需要在内部将点转换为虚拟域。用户可以通过SetFixedSampledPointSet()定义点集，通过调用SetUsedFixedSampledPointSet()来使用点集。
 
  CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
  optimizer->AddObserver(itk::IterationEvent(), observer);
 
  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);

//此外，用户可以在每个等级为配准框架定义采样百分比和采样策略，而不是直接处理度量。在这种情况下，配准过滤器根据输入strategy(REGULAR、RANDOM)管理fixed图像空间上的采样操作，并在内部将采样点集传递给度量。
  RegistrationType::MetricSamplingStrategyEnum samplingStrategy =
    RegistrationType::MetricSamplingStrategyEnum::RANDOM;

//所使用的空间样本数量取决于图像的内容。如果图像是平滑的，不包含很多细节，空间样本的数量通常可以低至固定图像中像素总数的1%。另一方面，如果图像很详细，则可能需要使用更高的比例，比如20%到50%。增加样本的数量可以提高度量的平滑性，因此，当这个度量与依赖度量值连续性的优化器一起使用时，会有所帮助。当然，这样做的代价是，在每次度量的评估中，大量的样本会导致更长的计算时间。
 //使度量达到极限的一种机制是禁用采样并使用FixedImageRegion中的所有像素。这可以通过SetUseSampledPointSet(false)方法完成。您可能想要尝试这个选项，只有当您正在微调所有其他参数的注册。在当前的例子中我们不使用这个方法。

  double samplingPercentage = 0.20;//实验证明，样本数不是配准过程的关键参数。当您开始微调您自己的配准过程时，您应该开始使用较高的样本数量，例如在固定图像中像素数量的20%到50%的范围内。一旦你成功地配准了你的图像，你就可以逐步减少样本的数量，直到你找到一个很好的折中方法来计算一个度量的评价。请注意，如果值中的噪声导致优化器需要更多的迭代来收敛，那么对度量进行非常快的评估是没有用处的。然后，您必须在迭代过程中研究度量值的行为。
 //在ITKv4中，单个虚拟域或空间样本点集用于注册过程的所有迭代。单个样本集的使用会产生平稳的成本函数，可以改善优化程序的功能。空间点集是伪随机生成的。为了获得可重现的结果，应设置一个整数种子。
  registration->SetMetricSamplingStrategy(samplingStrategy);
  registration->SetMetricSamplingPercentage(samplingPercentage);

  registration->MetricSamplingReinitializeSeed(121213);

  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;
  }

  TransformType::ParametersType finalParameters =
    registration->GetOutput()->Get()->GetParameters();

  double TranslationAlongX = finalParameters[0];
  double TranslationAlongY = finalParameters[1];
 
  unsigned int numberOfIterations = optimizer->GetCurrentIteration();

  double bestValue = optimizer->GetValue();

  std::cout << std::endl;
  std::cout << "Result = " << std::endl;
  std::cout << " Translation X = " << TranslationAlongX << std::endl;
  std::cout << " Translation Y = " << TranslationAlongY << std::endl;
  std::cout << " Iterations    = " << numberOfIterations << std::endl;
  std::cout << " Metric value  = " << bestValue << std::endl;
  std::cout << " Stop Condition  = "
            << optimizer->GetStopConditionDescription() << std::endl;

  using ResampleFilterType =
    itk::ResampleImageFilter<MovingImageType, FixedImageType>;

  ResampleFilterType::Pointer resample = ResampleFilterType::New();

  resample->SetTransform(registration->GetTransform());
  resample->SetInput(movingImageReader->GetOutput());

  FixedImageType::Pointer fixedImage = fixedImageReader->GetOutput();

  PixelType defaultPixelValue = 100;

  if (argc > 4)
  {
    defaultPixelValue = std::stoi(argv[4]);
  }

  resample->SetSize(fixedImage->GetLargestPossibleRegion().GetSize());
  resample->SetOutputOrigin(fixedImage->GetOrigin());
  resample->SetOutputSpacing(fixedImage->GetSpacing());
  resample->SetOutputDirection(fixedImage->GetDirection());
  resample->SetDefaultPixelValue(defaultPixelValue);


  using OutputPixelType = unsigned char;

  using OutputImageType = itk::Image<OutputPixelType, Dimension>;

  using CastFilterType =
    itk::CastImageFilter<FixedImageType, OutputImageType>;

  using WriterType = itk::ImageFileWriter<OutputImageType>;

  WriterType::Pointer     writer = WriterType::New();
  CastFilterType::Pointer caster = CastFilterType::New();

  writer->SetFileName(argv[3]);

  caster->SetInput(resample->GetOutput());
  writer->SetInput(caster->GetOutput());
  writer->Update();

  using CheckerBoardFilterType = itk::CheckerBoardImageFilter<FixedImageType>;

  CheckerBoardFilterType::Pointer checker = CheckerBoardFilterType::New();

  checker->SetInput1(fixedImage);
  checker->SetInput2(resample->GetOutput());

  caster->SetInput(checker->GetOutput());
  writer->SetInput(caster->GetOutput());

  resample->SetDefaultPixelValue(0);

  TransformType::Pointer identityTransform = TransformType::New();
  identityTransform->SetIdentity();
  resample->SetTransform(identityTransform);

  if (argc > 5)
  {
    writer->SetFileName(argv[5]);
    writer->Update();
  }


  resample->SetTransform(registration->GetTransform());
  if (argc > 6)
  {
    writer->SetFileName(argv[6]);
    writer->Update();
  }

  return EXIT_SUCCESS;
}

图3.12:优化器每次迭代的转换和度量值序列。

图3.12(左上)显示了优化器搜索参数空间时紧随其后的平移序列。右上角的图显示了优化器最后一次迭代的收敛池。同一图的底部显示了优化器搜索参数空间时计算的度量值序列。

但是，您必须注意，在优化参数的微调过程中涉及到许多重要问题。例如，在估计相互信息中使用的桶的数量对优化器的性能有显著的影响。为了演示这一效果，在执行相同的示例时，使用了从10到30个不同的箱子数量值。如果您重复这个实验，您将注意到，根据使用的容器数量，优化器的路径可能会在局部极小值的早期被困住。图3.13显示了多个路径。