Image Registration Method Displacement 1

Overview

Code

using System;
using itk.simple;

namespace itk.simple.examples
{
    class IterationUpdate : Command
    {
        private ImageRegistrationMethod m_Method;

        public IterationUpdate(ImageRegistrationMethod m)
        {
            m_Method = m;
        }

        public override void Execute()
        {
            if (m_Method.GetOptimizerIteration() == 0)
            {
                Console.WriteLine(String.Format("\tLevel: {0,3}", m_Method.GetCurrentLevel()));
                VectorDouble scales = m_Method.GetOptimizerScales();
                Console.Write(String.Format("\tScales: [{0:F5}", scales[0]));
                for (int i = 1; i < scales.Count; i++)
                {
                    Console.Write(String.Format(", {0:F5}", scales[i]));
                }
                Console.WriteLine("]");
            }
            Console.WriteLine(String.Format("#{0}", m_Method.GetOptimizerIteration()));
            Console.WriteLine(String.Format("\tMetric Value: {0:F5}", m_Method.GetMetricValue()));
            Console.WriteLine(String.Format("\tLearning Rate: {0:F5}", m_Method.GetOptimizerLearningRate()));
            if (m_Method.GetOptimizerConvergenceValue() != Double.MaxValue)
            {
                Console.WriteLine(String.Format("\tConvergence Value: {0:E}", m_Method.GetOptimizerConvergenceValue()));
            }
        }
    }

    class MultiResolutionIterationUpdate : Command
    {
        private ImageRegistrationMethod m_Method;

        public MultiResolutionIterationUpdate(ImageRegistrationMethod m)
        {
            m_Method = m;
        }

        public override void Execute()
        {
            Console.WriteLine(String.Format("\tStop Condition: {0}", m_Method.GetOptimizerStopConditionDescription()));
            Console.WriteLine("============= Resolution Change =============");
        }
    }

    class ImageRegistrationMethodDisplacement1
    {
        static void Main(string[] args)
        {
            if (args.Length < 3)
            {
                Console.WriteLine("Usage: {0} <fixedImageFile> <movingImageFile> <outputTransformFile>",
                    System.AppDomain.CurrentDomain.FriendlyName);
                return;
            }

            Image fixedImage = SimpleITK.ReadImage(args[0], PixelIDValueEnum.sitkFloat32);
            Image movingImage = SimpleITK.ReadImage(args[1], PixelIDValueEnum.sitkFloat32);

            Transform initialTx = SimpleITK.CenteredTransformInitializer(fixedImage, movingImage,
                new AffineTransform(fixedImage.GetDimension()));

            ImageRegistrationMethod R = new ImageRegistrationMethod();

            // First registration stage with affine transform
            VectorUInt32 shrinkFactors = new VectorUInt32();
            shrinkFactors.Add(3); shrinkFactors.Add(2); shrinkFactors.Add(1);

            VectorDouble smoothingSigmas = new VectorDouble();
            smoothingSigmas.Add(2.0); smoothingSigmas.Add(1.0); smoothingSigmas.Add(1.0);

            R.SetShrinkFactorsPerLevel(shrinkFactors);
            R.SetSmoothingSigmasPerLevel(smoothingSigmas);

            R.SetMetricAsJointHistogramMutualInformation(20);
            R.MetricUseFixedImageGradientFilterOff();

            double learningRate = 1.0;
            uint numberOfIterations = 100;
            double convergenceMinimumValue = 1e-6;
            uint convergenceWindowSize = 10;
            ImageRegistrationMethod.EstimateLearningRateType estimateLearningRate =
                ImageRegistrationMethod.EstimateLearningRateType.EachIteration;
            R.SetOptimizerAsGradientDescent(learningRate, numberOfIterations,
                                           convergenceMinimumValue, convergenceWindowSize, estimateLearningRate);

            R.SetOptimizerScalesFromPhysicalShift();
            R.SetInitialTransform(initialTx);
            R.SetInterpolator(InterpolatorEnum.sitkLinear);

            IterationUpdate cmd = new IterationUpdate(R);
            R.AddCommand(EventEnum.sitkIterationEvent, cmd);

            MultiResolutionIterationUpdate cmd2 = new MultiResolutionIterationUpdate(R);
            R.AddCommand(EventEnum.sitkMultiResolutionIterationEvent, cmd2);

            Transform outTx1 = R.Execute(fixedImage, movingImage);

            Console.WriteLine("-------");
            Console.WriteLine(outTx1.ToString());
            Console.WriteLine(String.Format("Optimizer stop condition: {0}", R.GetOptimizerStopConditionDescription()));
            Console.WriteLine(String.Format(" Iteration: {0}", R.GetOptimizerIteration()));
            Console.WriteLine(String.Format(" Metric value: {0}", R.GetMetricValue()));

            // Second registration stage with displacement field
            Image displacementField = new Image(fixedImage.GetSize(), PixelIDValueEnum.sitkVectorFloat64);
            displacementField.CopyInformation(fixedImage);
            DisplacementFieldTransform displacementTx = new DisplacementFieldTransform(displacementField);
            double varianceForUpdateField = 0.0;
            double varianceForTotalField = 1.5;
            displacementTx.SetSmoothingGaussianOnUpdate(varianceForUpdateField, varianceForTotalField);

            R.SetMovingInitialTransform(outTx1);
            R.SetInitialTransform(displacementTx, true);

            R.SetMetricAsANTSNeighborhoodCorrelation(4);
            R.MetricUseFixedImageGradientFilterOff();

            // Reset shrink factors and smoothing sigmas for second stage
            shrinkFactors = new VectorUInt32();
            shrinkFactors.Add(3); shrinkFactors.Add(2); shrinkFactors.Add(1);

            smoothingSigmas = new VectorDouble();
            smoothingSigmas.Add(2.0); smoothingSigmas.Add(1.0); smoothingSigmas.Add(1.0);

            R.SetShrinkFactorsPerLevel(shrinkFactors);
            R.SetSmoothingSigmasPerLevel(smoothingSigmas);

            R.SetOptimizerScalesFromPhysicalShift();

            learningRate = 1.0;
            numberOfIterations = 300;
            convergenceMinimumValue = 1e-6;
            convergenceWindowSize = 10;
            estimateLearningRate = ImageRegistrationMethod.EstimateLearningRateType.EachIteration;
            R.SetOptimizerAsGradientDescent(learningRate, numberOfIterations,
                                           convergenceMinimumValue, convergenceWindowSize, estimateLearningRate);

            R.Execute(fixedImage, movingImage);

            Console.WriteLine("-------");
            Console.WriteLine(displacementTx.ToString());
            Console.WriteLine(String.Format("Optimizer stop condition: {0}", R.GetOptimizerStopConditionDescription()));
            Console.WriteLine(String.Format(" Iteration: {0}", R.GetOptimizerIteration()));
            Console.WriteLine(String.Format(" Metric value: {0}", R.GetMetricValue()));

            VectorOfTransform transforms = new VectorOfTransform();
            transforms.Add(outTx1);
            transforms.Add(displacementTx);
            CompositeTransform outTx = new CompositeTransform(transforms);
            SimpleITK.WriteTransform(outTx, args[2]);
        }
    }
}

#include <SimpleITK.h>
#include <iostream>
#include <stdlib.h>
#include <iomanip>

namespace sitk = itk::simple;


class IterationUpdate : public sitk::Command
{
public:
  IterationUpdate(const sitk::ImageRegistrationMethod & m)
    : m_Method(m)
  {}

  void
  Execute() override
  {
    // use sitk's output operator for std::vector etc..
    using sitk::operator<<;

    // stash the stream state
    std::ios state(NULL);
    state.copyfmt(std::cout);
    std::cout << std::fixed << std::setfill(' ') << std::setprecision(5);
    if (m_Method.GetOptimizerIteration() == 0)
    {
      std::cout << "\tLevel: " << std::setw(3) << m_Method.GetCurrentLevel() << std::endl;
      std::cout << "\tScales: " << m_Method.GetOptimizerScales() << std::endl;
    }
    std::cout << '#' << m_Method.GetOptimizerIteration() << std::endl;
    std::cout << "\tMetric Value: " << m_Method.GetMetricValue() << std::endl;
    std::cout << "\tLearning Rate: " << m_Method.GetOptimizerLearningRate() << std::endl;
    if (m_Method.GetOptimizerConvergenceValue() != std::numeric_limits<double>::max())
    {
      std::cout << "\tConvergence Value: " << std::scientific << m_Method.GetOptimizerConvergenceValue() << std::endl;
    }
    std::cout.copyfmt(state);
  }

private:
  const sitk::ImageRegistrationMethod & m_Method;
};

class MultiResolutionIterationUpdate : public sitk::Command
{
public:
  MultiResolutionIterationUpdate(const sitk::ImageRegistrationMethod & m)
    : m_Method(m)
  {}

  void
  Execute() override
  {
    // use sitk's output operator for std::vector etc..
    using sitk::operator<<;

    // stash the stream state
    std::ios state(NULL);
    state.copyfmt(std::cout);
    std::cout << std::fixed << std::setfill(' ') << std::setprecision(5);
    std::cout << "\tStop Condition: " << m_Method.GetOptimizerStopConditionDescription() << std::endl;
    std::cout << "============= Resolution Change =============" << std::endl;
    std::cout.copyfmt(state);
  }

private:
  const sitk::ImageRegistrationMethod & m_Method;
};


int
main(int argc, char * argv[])
{

  if (argc < 4)
  {
    std::cerr << "Usage: " << argv[0] << " <fixedImageFilter> <movingImageFile> <outputTransformFile>" << std::endl;
    return 1;
  }

  sitk::Image fixed = sitk::ReadImage(argv[1], sitk::sitkFloat32);

  sitk::Image moving = sitk::ReadImage(argv[2], sitk::sitkFloat32);

  sitk::Transform initialTx =
    sitk::CenteredTransformInitializer(fixed, moving, sitk::AffineTransform(fixed.GetDimension()));

  sitk::ImageRegistrationMethod R;

  {
    std::vector<unsigned int> shrinkFactors = { 3, 2, 1 };

    std::vector<double> smoothingSigmas = { 2.0, 1.0, 1.0 };

    R.SetShrinkFactorsPerLevel(shrinkFactors);
    R.SetSmoothingSigmasPerLevel(smoothingSigmas);
  }

  R.SetMetricAsJointHistogramMutualInformation(20);
  R.MetricUseFixedImageGradientFilterOff();

  {
    double                                                  learningRate = 1.0;
    unsigned int                                            numberOfIterations = 100;
    double                                                  convergenceMinimumValue = 1e-6;
    unsigned int                                            convergenceWindowSize = 10;
    sitk::ImageRegistrationMethod::EstimateLearningRateType estimateLearningRate = R.EachIteration;
    R.SetOptimizerAsGradientDescent(
      learningRate, numberOfIterations, convergenceMinimumValue, convergenceWindowSize, estimateLearningRate);
  }
  R.SetOptimizerScalesFromPhysicalShift();

  R.SetInitialTransform(initialTx);

  R.SetInterpolator(sitk::sitkLinear);

  IterationUpdate cmd(R);
  R.AddCommand(sitk::sitkIterationEvent, cmd);

  MultiResolutionIterationUpdate cmd2(R);
  R.AddCommand(sitk::sitkMultiResolutionIterationEvent, cmd2);

  sitk::Transform outTx1 = R.Execute(fixed, moving);

  std::cout << "-------" << std::endl;
  std::cout << outTx1.ToString() << std::endl;
  std::cout << "Optimizer stop condition: " << R.GetOptimizerStopConditionDescription() << std::endl;
  std::cout << " Iteration: " << R.GetOptimizerIteration() << std::endl;
  std::cout << " Metric value: " << R.GetMetricValue() << std::endl;


  sitk::Image displacementField = sitk::Image(fixed.GetSize(), sitk::sitkVectorFloat64);
  displacementField.CopyInformation(fixed);
  sitk::DisplacementFieldTransform displacementTx(displacementField);
  const double                     varianceForUpdateField = 0.0;
  const double                     varianceForTotalField = 1.5;
  displacementTx.SetSmoothingGaussianOnUpdate(varianceForUpdateField, varianceForTotalField);


  R.SetMovingInitialTransform(outTx1);
  R.SetInitialTransform(displacementTx, true);

  R.SetMetricAsANTSNeighborhoodCorrelation(4);
  R.MetricUseFixedImageGradientFilterOff();

  {
    std::vector<unsigned int> shrinkFactors = { 3, 2, 1 };

    std::vector<double> smoothingSigmas = { 2.0, 1.0, 1.0 };

    R.SetShrinkFactorsPerLevel(shrinkFactors);
    R.SetSmoothingSigmasPerLevel(smoothingSigmas);
  }

  R.SetOptimizerScalesFromPhysicalShift();

  {
    double                                                  learningRate = 1.0;
    unsigned int                                            numberOfIterations = 300;
    double                                                  convergenceMinimumValue = 1e-6;
    unsigned int                                            convergenceWindowSize = 10;
    sitk::ImageRegistrationMethod::EstimateLearningRateType estimateLearningRate = R.EachIteration;
    R.SetOptimizerAsGradientDescent(
      learningRate, numberOfIterations, convergenceMinimumValue, convergenceWindowSize, estimateLearningRate);
  }
  R.Execute(fixed, moving);


  std::cout << "-------" << std::endl;
  std::cout << displacementTx.ToString() << std::endl;
  std::cout << "Optimizer stop condition: " << R.GetOptimizerStopConditionDescription() << std::endl;
  std::cout << " Iteration: " << R.GetOptimizerIteration() << std::endl;
  std::cout << " Metric value: " << R.GetMetricValue() << std::endl;


  sitk::CompositeTransform outTx({ outTx1, displacementTx });
  sitk::WriteTransform(outTx, argv[3]);

  return 0;
}

import org.itk.simple.*;

class IterationUpdate extends Command {
    private ImageRegistrationMethod method;

    public IterationUpdate(ImageRegistrationMethod m) {
        method = m;
    }

    public void execute() {
        if (method.getOptimizerIteration() == 0) {
            System.out.printf("\tLevel: %3d%n", method.getCurrentLevel());
            VectorDouble scales = method.getOptimizerScales();
            System.out.printf("\tScales: [%.5f", scales.get(0));
            for (int i = 1; i < scales.size(); i++) {
                System.out.printf(", %.5f", scales.get(i));
            }
            System.out.println("]");
        }
        System.out.printf("#%d%n", method.getOptimizerIteration());
        System.out.printf("\tMetric Value: %.5f%n", method.getMetricValue());
        System.out.printf("\tLearning Rate: %.5f%n", method.getOptimizerLearningRate());
        if (method.getOptimizerConvergenceValue() != Double.MAX_VALUE) {
            System.out.printf("\tConvergence Value: %E%n", method.getOptimizerConvergenceValue());
        }
    }
}

class MultiResolutionIterationUpdate extends Command {
    private ImageRegistrationMethod method;

    public MultiResolutionIterationUpdate(ImageRegistrationMethod m) {
        method = m;
    }

    public void execute() {
        System.out.printf("\tStop Condition: %s%n", method.getOptimizerStopConditionDescription());
        System.out.println("============= Resolution Change =============");
    }
}

class ImageRegistrationMethodDisplacement1 {

    public static void main(String[] args) {

        if (args.length < 3) {
            System.out.println("Usage: ImageRegistrationMethodDisplacement1 <fixedImageFile> <movingImageFile> <outputTransformFile>");
            return;
        }

        Image fixedImage = SimpleITK.readImage(args[0], PixelIDValueEnum.sitkFloat32);
        Image movingImage = SimpleITK.readImage(args[1], PixelIDValueEnum.sitkFloat32);

        Transform initialTx = SimpleITK.centeredTransformInitializer(fixedImage, movingImage,
            new AffineTransform(fixedImage.getDimension()));

        ImageRegistrationMethod R = new ImageRegistrationMethod();

        // First registration stage with affine transform
        VectorUInt32 shrinkFactors = new VectorUInt32();
        shrinkFactors.add(3L); shrinkFactors.add(2L); shrinkFactors.add(1L);

        VectorDouble smoothingSigmas = new VectorDouble();
        smoothingSigmas.add(2.0); smoothingSigmas.add(1.0); smoothingSigmas.add(1.0);

        R.setShrinkFactorsPerLevel(shrinkFactors);
        R.setSmoothingSigmasPerLevel(smoothingSigmas);

        R.setMetricAsJointHistogramMutualInformation(20);
        R.metricUseFixedImageGradientFilterOff();

        double learningRate = 1.0;
        int numberOfIterations = 100;
        double convergenceMinimumValue = 1e-6;
        int convergenceWindowSize = 10;
        ImageRegistrationMethod.EstimateLearningRateType estimateLearningRate =
            ImageRegistrationMethod.EstimateLearningRateType.EachIteration;
        R.setOptimizerAsGradientDescent(learningRate, numberOfIterations,
                                       convergenceMinimumValue, convergenceWindowSize, estimateLearningRate);

        R.setOptimizerScalesFromPhysicalShift();
        R.setInitialTransform(initialTx);
        R.setInterpolator(InterpolatorEnum.sitkLinear);

        IterationUpdate cmd = new IterationUpdate(R);
        R.addCommand(EventEnum.sitkIterationEvent, cmd);

        MultiResolutionIterationUpdate cmd2 = new MultiResolutionIterationUpdate(R);
        R.addCommand(EventEnum.sitkMultiResolutionIterationEvent, cmd2);

        Transform outTx1 = R.execute(fixedImage, movingImage);

        System.out.println("-------");
        System.out.println(outTx1.toString());
        System.out.printf("Optimizer stop condition: %s%n", R.getOptimizerStopConditionDescription());
        System.out.printf(" Iteration: %d%n", R.getOptimizerIteration());
        System.out.printf(" Metric value: %f%n", R.getMetricValue());

        // Second registration stage with displacement field
        Image displacementField = new Image(fixedImage.getSize(), PixelIDValueEnum.sitkVectorFloat64);
        displacementField.copyInformation(fixedImage);
        DisplacementFieldTransform displacementTx = new DisplacementFieldTransform(displacementField);
        double varianceForUpdateField = 0.0;
        double varianceForTotalField = 1.5;
        displacementTx.setSmoothingGaussianOnUpdate(varianceForUpdateField, varianceForTotalField);

        R.setMovingInitialTransform(outTx1);
        R.setInitialTransform(displacementTx, true);

        R.setMetricAsANTSNeighborhoodCorrelation(4);
        R.metricUseFixedImageGradientFilterOff();

        // Reset shrink factors and smoothing sigmas for second stage
        shrinkFactors = new VectorUInt32();
        shrinkFactors.add(3L); shrinkFactors.add(2L); shrinkFactors.add(1L);

        smoothingSigmas = new VectorDouble();
        smoothingSigmas.add(2.0); smoothingSigmas.add(1.0); smoothingSigmas.add(1.0);

        R.setShrinkFactorsPerLevel(shrinkFactors);
        R.setSmoothingSigmasPerLevel(smoothingSigmas);

        R.setOptimizerScalesFromPhysicalShift();

        learningRate = 1.0;
        numberOfIterations = 300;
        convergenceMinimumValue = 1e-6;
        convergenceWindowSize = 10;
        estimateLearningRate = ImageRegistrationMethod.EstimateLearningRateType.EachIteration;
        R.setOptimizerAsGradientDescent(learningRate, numberOfIterations,
                                       convergenceMinimumValue, convergenceWindowSize, estimateLearningRate);

        R.execute(fixedImage, movingImage);

        System.out.println("-------");
        System.out.println(displacementTx.toString());
        System.out.printf("Optimizer stop condition: %s%n", R.getOptimizerStopConditionDescription());
        System.out.printf(" Iteration: %d%n", R.getOptimizerIteration());
        System.out.printf(" Metric value: %f%n", R.getMetricValue());

        VectorOfTransform transforms = new VectorOfTransform();
        transforms.add(outTx1);
        transforms.add(displacementTx);
        CompositeTransform outTx = new CompositeTransform(transforms);
        SimpleITK.writeTransform(outTx, args[2]);
    }
}

#!/usr/bin/env python

import sys
import os
import SimpleITK as sitk


def command_iteration(method):
    """ Callback invoked each iteration. """
    if method.GetOptimizerIteration() == 0:
        print(f"\tLevel: {method.GetCurrentLevel()}")
        print(f"\tScales: {method.GetOptimizerScales()}")
    print(f"#{method.GetOptimizerIteration()}")
    print(f"\tMetric Value: {method.GetMetricValue():10.5f}")
    print(f"\tLearningRate: {method.GetOptimizerLearningRate():10.5f}")
    if method.GetOptimizerConvergenceValue() != sys.float_info.max:
        print("\tConvergence Value: " + f"{method.GetOptimizerConvergenceValue():.5e}")


def command_multiresolution_iteration(method):
    """ Callback invoked at the end of each multi-resolution level. """
    print(f"\tStop Condition: {method.GetOptimizerStopConditionDescription()}")
    print("============= Resolution Change =============")


if len(sys.argv) < 4:
    print(
        "Usage:",
        sys.argv[0],
        "<fixedImageFilter> <movingImageFile>",
        "<outputTransformFile>",
    )
    sys.exit(1)

fixed = sitk.ReadImage(sys.argv[1], sitk.sitkFloat32)

moving = sitk.ReadImage(sys.argv[2], sitk.sitkFloat32)

initialTx = sitk.CenteredTransformInitializer(
    fixed, moving, sitk.AffineTransform(fixed.GetDimension())
)

R = sitk.ImageRegistrationMethod()

R.SetShrinkFactorsPerLevel([3, 2, 1])
R.SetSmoothingSigmasPerLevel([2, 1, 1])

R.SetMetricAsJointHistogramMutualInformation(20)
R.MetricUseFixedImageGradientFilterOff()

R.SetOptimizerAsGradientDescent(
    learningRate=1.0,
    numberOfIterations=100,
    estimateLearningRate=R.EachIteration,
)
R.SetOptimizerScalesFromPhysicalShift()

R.SetInitialTransform(initialTx)

R.SetInterpolator(sitk.sitkLinear)

R.AddCommand(sitk.sitkIterationEvent, lambda: command_iteration(R))
R.AddCommand(
    sitk.sitkMultiResolutionIterationEvent,
    lambda: command_multiresolution_iteration(R),
)

outTx1 = R.Execute(fixed, moving)

print("-------")
print(outTx1)
print(f"Optimizer stop condition: {R.GetOptimizerStopConditionDescription()}")
print(f" Iteration: {R.GetOptimizerIteration()}")
print(f" Metric value: {R.GetMetricValue()}")

displacementField = sitk.Image(fixed.GetSize(), sitk.sitkVectorFloat64)
displacementField.CopyInformation(fixed)
displacementTx = sitk.DisplacementFieldTransform(displacementField)
del displacementField
displacementTx.SetSmoothingGaussianOnUpdate(
    varianceForUpdateField=0.0, varianceForTotalField=1.5
)

R.SetMovingInitialTransform(outTx1)
R.SetInitialTransform(displacementTx, inPlace=True)

R.SetMetricAsANTSNeighborhoodCorrelation(4)
R.MetricUseFixedImageGradientFilterOff()

R.SetShrinkFactorsPerLevel([3, 2, 1])
R.SetSmoothingSigmasPerLevel([2, 1, 1])

R.SetOptimizerScalesFromPhysicalShift()
R.SetOptimizerAsGradientDescent(
    learningRate=1,
    numberOfIterations=300,
    estimateLearningRate=R.EachIteration,
)

R.Execute(fixed, moving)

print("-------")
print(displacementTx)
print(f"Optimizer stop condition: {R.GetOptimizerStopConditionDescription()}")
print(f" Iteration: {R.GetOptimizerIteration()}")
print(f" Metric value: {R.GetMetricValue()}")

compositeTx = sitk.CompositeTransform([outTx1, displacementTx])
sitk.WriteTransform(compositeTx, sys.argv[3])

if "SITK_NOSHOW" not in os.environ:
    sitk.Show(displacementTx.GetDisplacementField(), "Displacement Field")

    resampler = sitk.ResampleImageFilter()
    resampler.SetReferenceImage(fixed)
    resampler.SetInterpolator(sitk.sitkLinear)
    resampler.SetDefaultPixelValue(100)
    resampler.SetTransform(compositeTx)

    out = resampler.Execute(moving)
    simg1 = sitk.Cast(sitk.RescaleIntensity(fixed), sitk.sitkUInt8)
    simg2 = sitk.Cast(sitk.RescaleIntensity(out), sitk.sitkUInt8)
    cimg = sitk.Compose(simg1, simg2, simg1 // 2.0 + simg2 // 2.0)
    sitk.Show(cimg, "ImageRegistration1 Composition")

# Run with:
#
# Rscript --vanilla ImageRegistrationMethodDisplacement1.R fixedImageFilter movingImageFile outputTransformFile
#

library(SimpleITK)

commandIteration <- function(method)
{
    if (method$GetOptimizerIteration()==0) {
        cat("\tLevel:", method$GetCurrentLevel(),
            "\n\tScales:", method$GetOptimizerScales(), "\n")
    }
    msg <- paste("#", method$GetOptimizerIteration(),
                 "\n\tMetric Value: ", method$GetMetricValue(),
                 "\n\tLearning Rate: ", method$GetOptimizerLearningRate(),
                 '\n', sep="")
    cat(msg)
}

commandMultiIteration <- function(method)
{
    msg <- paste("\tStop Condition: ", method$GetOptimizerStopConditionDescription(),
                 "\n============= Resolution Change =============\n", sep="")
    cat(msg)
}

args <- commandArgs( TRUE )

if (length(args) != 3) {
    stop("3 arguments expected - fixedImageFilter, movingImageFile, outputTransformFile")
}

fixed <- ReadImage(args[[1]], 'sitkFloat32')
moving <- ReadImage(args[[2]], 'sitkFloat32')

initialTx <- CenteredTransformInitializer(fixed, moving, AffineTransform(fixed$GetDimension()))

R <- ImageRegistrationMethod()

R$SetShrinkFactorsPerLevel(c(3,2,1))
R$SetSmoothingSigmasPerLevel(c(2,1,1))

R$SetMetricAsJointHistogramMutualInformation(20)
R$MetricUseFixedImageGradientFilterOff()

R$SetOptimizerAsGradientDescent(learningRate=1.0,
                                numberOfIterations=100,
                                convergenceMinimumValue=1e-6,
                                convergenceWindowSize=10,
                                estimateLearningRate = 'EachIteration')
R$SetOptimizerScalesFromPhysicalShift()

R$SetInitialTransform(initialTx)

R$SetInterpolator('sitkLinear')

R$AddCommand( 'sitkIterationEvent', function() commandIteration(R) )
R$AddCommand( 'sitkMultiResolutionIterationEvent', function() commandMultiIteration(R) )

outTx1 <- R$Execute(fixed, moving)

cat("-------\n")
outTx1
cat("Optimizer stop condition:", R$GetOptimizerStopConditionDescription(), '\n')
cat("Iteration:", R$GetOptimizerIteration(), '\n')
cat("Metric value:", R$GetMetricValue(), '\n')


displacementField <- Image(fixed$GetSize(), 'sitkVectorFloat64')
displacementField$CopyInformation(fixed)
displacementTx <- DisplacementFieldTransform(displacementField)
rm(displacementField)
displacementTx$SetSmoothingGaussianOnUpdate(varianceForUpdateField=0.0,
                                            varianceForTotalField=1.5)

R$SetMovingInitialTransform(outTx1)
R$SetInitialTransform(displacementTx, inPlace=TRUE)

R$SetMetricAsANTSNeighborhoodCorrelation(4)
R$MetricUseFixedImageGradientFilterOff()

R$SetShrinkFactorsPerLevel(c(3,2,1))
R$SetSmoothingSigmasPerLevel(c(2,1,1))

R$SetOptimizerScalesFromPhysicalShift()
R$SetOptimizerAsGradientDescent(learningRate=1,
                                numberOfIterations=300,
                                convergenceMinimumValue=1e-6,
                                convergenceWindowSize=10,
                                estimateLearningRate = 'EachIteration')

R$Execute(fixed, moving)



cat("-------\n")
displacementTx
cat("Optimizer stop condition:", R$GetOptimizerStopConditionDescription(), '\n')
cat("Iteration:", R$GetOptimizerIteration(), '\n')
cat("Metric value:", R$GetMetricValue(), '\n')

outTx <- CompositeTransform(outTx1$GetDimension())
outTx$AddTransform(outTx1)
outTx$AddTransform(displacementTx)

WriteTransform(outTx,  args[[3]])