itk registration 15

3.15 Visualizing Deformation fields

3.15.1 Visualizing 2D deformation fields

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

 

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

#include "itkImageRegionIterator.h"

#include "itkDemonsRegistrationFilter.h"
#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;
  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>;
  MatchingFilterType::Pointer matcher = MatchingFilterType::New();
  
  matcher->SetInput(movingImageCaster->GetOutput());
  matcher->SetReferenceImage(fixedImageCaster->GetOutput());
 
  matcher->SetNumberOfHistogramLevels(1024);
  matcher->SetNumberOfMatchPoints(7);
  
  matcher->ThresholdAtMeanIntensityOn();
 
  using VectorPixelType = itk::Vector<float, Dimension>;
  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());
  filter->SetMovingImage(matcher->GetOutput());

  filter->SetNumberOfIterations(50);
  filter->SetStandardDeviations(1.0);
 
  filter->Update();
  
  using OutputPixelType = unsigned char;
  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>;
  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;
}

3.15.2 Visualizing 3D deformation fields

让我们创建一个三维形变场。我们将使用Thin Plate Splines来扭曲一个3D数据集，并创建一个形变场。我们将选择一组点landmark，并将它们平移landmarks对应关系的规范。请注意，出于说明目的随机抽取了landmark，但并非意图描绘真实的变形。landmark可以用几种方法来产生形变场。大多数技术会最小化一些表示变形场不规则性的正则化函数，这通常是该场的空间导数的某些函数。这里我们用薄板样条。薄板样条使正则化函数最小化：

下标表示f的偏导数。

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

#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"


// NOTE: the LBFGSOptimizerv4 does not invoke events

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

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

  std::cout << "Initial Parameters = " << std::endl;
  std::cout << transform->GetParameters() << std::endl;
  
  registration->SetInitialTransform(transform);
  registration->InPlaceOn();

  registration->SetFixedImage(fixedImage);
  registration->SetMovingImage(movingImageReader->GetOutput());
 
  using ScalesEstimatorType =
    itk::RegistrationParameterScalesFromPhysicalShift<MetricType>;
  ScalesEstimatorType::Pointer scalesEstimator = ScalesEstimatorType::New();
  scalesEstimator->SetMetric(metric);
  scalesEstimator->SetTransformForward(true);
  scalesEstimator->SetSmallParameterVariation(1.0);
  
  optimizer->SetGradientConvergenceTolerance(5e-2);
  optimizer->SetLineSearchAccuracy(1.2);
  optimizer->SetDefaultStepLength(1.5);
  optimizer->TraceOn();
  optimizer->SetMaximumNumberOfFunctionEvaluations(1000);
  optimizer->SetScalesEstimator(scalesEstimator);

  //  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(1);
  shrinkFactorsPerLevel[0] = 1;

  RegistrationType::SmoothingSigmasArrayType smoothingSigmasPerLevel;
  smoothingSigmasPerLevel.SetSize(1);
  smoothingSigmasPerLevel[0] = 0;

  registration->SetNumberOfLevels(numberOfLevels);
  registration->SetSmoothingSigmasPerLevel(smoothingSigmasPerLevel);
  registration->SetShrinkFactorsPerLevel(shrinkFactorsPerLevel);


  // Add time and memory probes
  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;
  }

  // Report the time and memory taken by the registration
  chronometer.Report(std::cout);
  memorymeter.Report(std::cout);
 
  OptimizerType::ParametersType finalParameters = transform->GetParameters();

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


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