Image Registration Method Exhaustive

Overview

This script demonstrates the use of the Exhaustive optimizer in the ImageRegistrationMethod to estimate a good initial rotation position.

Because gradient descent base optimization can get stuck in local minima, a good initial transform is critical for reasonable results. Search a reasonable space on a grid with brute force may be a reliable way to get a starting location for further optimization.

The initial translation and center of rotation for the transform is initialized based on the first principle moments of the intensities of the image. Then in either 2D or 3D a Euler transform is used to exhaustively search a grid of the rotation space at a certain step size. The resulting transform is a reasonable guess where to start further registration.

Code

C#C++JavaPythonR

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)
            {
                VectorDouble scales = m_Method.GetOptimizerScales();
                Console.Write("Scales: [" + scales[0]);
                for (int i = 1; i < scales.Count; i++)
                {
                    Console.Write(", " + scales[i]);
                }
                Console.WriteLine("]");
            }

            VectorDouble pos = m_Method.GetOptimizerPosition();
            Console.Write("{0,3} = {1,7:F5} : [{2:F5}",
                         m_Method.GetOptimizerIteration(),
                         m_Method.GetMetricValue(),
                         pos[0]);
            for (int i = 1; i < pos.Count; i++)
            {
                Console.Write(", {0:F5}", pos[i]);
            }
            Console.WriteLine("]");
        }
    }

    class ImageRegistrationMethodExhaustive
    {
        static void Main(string[] args)
        {
            if (args.Length < 3)
            {
                Console.WriteLine("Usage: ImageRegistrationMethodExhaustive <fixedImageFile> <movingImageFile> <outputTransformFile>");
                return;
            }

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

            ImageRegistrationMethod R = new ImageRegistrationMethod();

            R.SetMetricAsMattesMutualInformation(50); // numberOfHistogramBins

            int samplePerAxis = 12;
            Transform tx = null;

            if (fixedImage.GetDimension() == 2)
            {
                tx = new Euler2DTransform();
                // Set the number of samples (radius) in each dimension, with a default step size of 1.0
                VectorUInt32 exhaustiveSteps = new VectorUInt32(new uint[] { (uint)(samplePerAxis / 2), 0, 0 });
                R.SetOptimizerAsExhaustive(exhaustiveSteps);
                // Utilize the scale to set the step size for each dimension
                VectorDouble scales = new VectorDouble(new double[] { 2.0 * Math.PI / samplePerAxis, 1.0, 1.0 });
                R.SetOptimizerScales(scales);
            }
            else if (fixedImage.GetDimension() == 3)
            {
                tx = new Euler3DTransform();
                VectorUInt32 exhaustiveSteps = new VectorUInt32(new uint[] {
                    (uint)(samplePerAxis / 2),
                    (uint)(samplePerAxis / 2),
                    (uint)(samplePerAxis / 4),
                    0, 0, 0
                });
                R.SetOptimizerAsExhaustive(exhaustiveSteps);
                VectorDouble scales = new VectorDouble(new double[] {
                    2.0 * Math.PI / samplePerAxis,
                    2.0 * Math.PI / samplePerAxis,
                    2.0 * Math.PI / samplePerAxis,
                    1.0, 1.0, 1.0
                });
                R.SetOptimizerScales(scales);
            }

            // Initialize the transform with a translation and the center of rotation from the moments of intensity.
            tx = SimpleITK.CenteredTransformInitializer(fixedImage, movingImage, tx);

            R.SetInitialTransform(tx);
            R.SetInterpolator(InterpolatorEnum.sitkLinear);

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

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

            Console.WriteLine("-------");
            Console.WriteLine(outTx.ToString());
            Console.WriteLine("Optimizer stop condition: " + R.GetOptimizerStopConditionDescription());
            Console.WriteLine(" Iteration: " + R.GetOptimizerIteration());
            Console.WriteLine(" Metric value: " + R.GetMetricValue());

            SimpleITK.WriteTransform(outTx, args[2]);

            if (Environment.GetEnvironmentVariable("SITK_NOSHOW") == null)
            {
                ResampleImageFilter resampler = new ResampleImageFilter();
                resampler.SetReferenceImage(fixedImage);
                resampler.SetInterpolator(InterpolatorEnum.sitkLinear);
                resampler.SetDefaultPixelValue(1);
                resampler.SetTransform(outTx);

                Image output = resampler.Execute(movingImage);

                Image simg1 = SimpleITK.Cast(SimpleITK.RescaleIntensity(fixedImage), PixelIDValueEnum.sitkUInt8);
                Image simg2 = SimpleITK.Cast(SimpleITK.RescaleIntensity(output), PixelIDValueEnum.sitkUInt8);
                Image cimg = SimpleITK.Compose(simg1, simg2, SimpleITK.Divide(SimpleITK.Add(simg1, simg2), 2));
                SimpleITK.Show(cimg, "ImageRegistrationExhaustive Composition");
            }
        }
    }
}