Image Registration Method BSpline 3

If you are not familiar with the SimpleITK registration framework we recommend that you read the registration overview before continuing with the example.

Overview

This example performs registration of multi-modality images with a multi-resolution BSpline approach.

A BSplineTransform usually has a large number of parameters which increases the complexity and duration of optimizing the deformation. The multi-resolution BSpline approach initially performs the registration at a lower resolution with fewer parameters at the first level and then adapts or resamples the BSpline control points to a higher resolution at subsequent levels. This transformation adaption is done concurrently with the ImageRegistrationMethod’s multi-level feature. This enables the setting of the shrink factor, smoothing sigmas, sampling percentage and BSpline resolution to vary per level to efficiently solve a diverse set of registration problems.

The initial transform is the low resolution BSpline. The scaling factor of each BSpline used per level is specified with the ImageRegistration::SetInitialBSpline method’s third argument as an integer array. The first value is usually 1 and it is reasonable to double the resolution at each registration level. For this example the last resolution is 5, so that the result is comparable to the transformation from the previous example.

It can be important to monitor and observe the transform at each level or iteration. When the “inplace” option is enabled, the transform passed as the initial transform will be up to date during the registration process which enables it to be used in event commands.

See also: Image Registration Method BSpline 1, Image Registration Method BSpline 2.

Code

C#C++JavaPythonR

using itk.simple;

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

        public IterationUpdate(ImageRegistrationMethod m, Transform bsplineTransform)
        {
            m_Method = m;
            m_BSplineTransform = bsplineTransform;
        }

        public override void Execute()
        {
            if (m_Method.GetOptimizerIteration() == 0)
            {
                // The BSpline is resized before the first optimizer
                // iteration is completed per level. Print the transform object
                // to show the adapted BSpline transform.
                Console.WriteLine(m_BSplineTransform.ToString());
            }

            Console.WriteLine(String.Format("{0,3} = {1,10:F5}",
                             m_Method.GetOptimizerIteration(),
                             m_Method.GetMetricValue()));
        }
    }

    class MultiResolutionIterationUpdate : Command
    {
        private ImageRegistrationMethod m_Method;

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

        public override void Execute()
        {
            if (m_Method.GetCurrentLevel() > 0)
            {
                Console.WriteLine(String.Format("Optimizer stop condition: {0}",
                                 m_Method.GetOptimizerStopConditionDescription()));
                Console.WriteLine(String.Format(" Iteration: {0}", m_Method.GetOptimizerIteration()));
                Console.WriteLine(String.Format(" Metric value: {0}", m_Method.GetMetricValue()));
            }

            Console.WriteLine("--------- Resolution Changing ---------");
        }
    }

    class ImageRegistrationMethodBSpline3
    {
        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);

            VectorUInt32 transformDomainMeshSize = new VectorUInt32();
            for (uint i = 0; i < fixedImage.GetDimension(); i++)
            {
                transformDomainMeshSize.Add(2);
            }

            BSplineTransform tx = SimpleITK.BSplineTransformInitializer(fixedImage, transformDomainMeshSize);

            Console.WriteLine(String.Format("Initial Number of Parameters: {0}", tx.GetNumberOfParameters()));

            ImageRegistrationMethod R = new ImageRegistrationMethod();
            R.SetMetricAsJointHistogramMutualInformation();

            double learningRate = 5.0;
            uint numberOfIterations = 100;
            double convergenceMinimumValue = 1e-4;
            uint convergenceWindowSize = 5;
            R.SetOptimizerAsGradientDescentLineSearch(learningRate, numberOfIterations,
                                                     convergenceMinimumValue, convergenceWindowSize);

            R.SetInterpolator(InterpolatorEnum.sitkLinear);

            VectorUInt32 scaleFactors = new VectorUInt32();
            scaleFactors.Add(1); scaleFactors.Add(2); scaleFactors.Add(5);
            R.SetInitialTransformAsBSpline(tx, true, scaleFactors);

            VectorUInt32 shrinkFactors = new VectorUInt32();
            shrinkFactors.Add(4); shrinkFactors.Add(2); shrinkFactors.Add(1);
            R.SetShrinkFactorsPerLevel(shrinkFactors);

            VectorDouble smoothingSigmas = new VectorDouble();
            smoothingSigmas.Add(4); smoothingSigmas.Add(2); smoothingSigmas.Add(1);
            R.SetSmoothingSigmasPerLevel(smoothingSigmas);

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

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

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

            Console.WriteLine("-------");
            Console.WriteLine(tx.ToString());
            Console.WriteLine(outTx.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()));

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

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

                Image outputImage = resampler.Execute(movingImage);
                Image simg1 = SimpleITK.Cast(SimpleITK.RescaleIntensity(fixedImage), PixelIDValueEnum.sitkUInt8);
                Image simg2 = SimpleITK.Cast(SimpleITK.RescaleIntensity(outputImage), PixelIDValueEnum.sitkUInt8);
                Image cimg = SimpleITK.Compose(simg1, simg2, SimpleITK.Divide(SimpleITK.Add(simg1, simg2), 2.0));

                SimpleITK.Show(cimg, "Image Registration Composition");
            }
        }
    }
}