Evaluation of Free Form Deformation and Demons Registration with Discontinuities

Medical image registration plays an important part in most today’s clinical procedures. Registration goal is to find transformation which warps one image into the space of another. Registration of moving organs in human body has a significant part in therapy planning. This task is harder in cases when one organ (tissue) slides along another, i.e. in a case of discontinuities in the motion field. Discontinuities introduce unwanted transformations which often lead to poor or unsatisfied registration results. In this paper we evaluate one form of discontinuities for two well-known and used registration algorithms namely Free Form Deformation and Demons.


INTRODUCTION
Image registration is the process of overlaying two or more images of the same scene taken at different times, from different viewpoints, and/or by different sensors.The goal of registration is to find spatial alignment of images [1].Establishing the correspondences of spatial information in medical images and equivalent structures in the body is fundamental to medical image interpretation and analysis.The images might be acquired with different sensors or the same sensor at different times.Radiologist with subjective judgment of relative size, shape, and spatial relationships of visible structures and physiology inferred from intensity distributions, uses this information for developing a diagnosis, planning a therapy, and monitoring the disease progression or response.During the previous twenty years a number of image registration algorithms were proposed.
In Thirion [2] was proposed a method for non/rigid image registration.Rueckert [3] proposed the use of Free Form Deformations for non/rigid image registration.During the time, several improvements regarding the precision and registration speed were developed.In Muyan-Ozcelik [4], GPU implementation of Demons algorithm was developed.Image registration is a high time consuming operation.In this paper, Demons registration for 3D CT lung images was developed on a Graphics Processing Units (GPU), wherein speed improvements of 55 times were reported.Better description of Demons algorithm can be found in Penec [5].In Gu [6], five variants of demons algorithm were developed on a GPU.In Zhao [7] a framework for capturing very large deformations was developed using demons algorithm.
In Shi [8] was proposed a sparse Free Form Deformation for which fine local details such as motion discontinuities could be captured.In Yang [9], demons algorithm with locally adaptive regularization was proposed.
In this paper, two widely used image registration algorithms are analysed for the case of motion discontinuities, and those are Free Form Deformations and Demons algorithm.

METHODS
Image registration process can be described as in Fig. 1.In general, image registration needs two images, here called moving and fixed.As an output of registration algorithm we have registered image.During registration, an algorithm iteratively transforms moving image to match fixed image.Final transformed image should be "similar" to fixed image (Fig. 2).Similarity is calculated using the criterion which usually compares corresponding pixels from input images, and also terminates the registration process.Image registration has applications in many fields; one addressed in this paper deals with medical imaging, specifically radiological imaging.In this section we will give a review for the two most used algorithms for non-rigid image registration.

Figure 1 Three components of the image registration algorithm
Figure 2 The transformation maps pixels from the coordinate system of first image to the coordinate system of another one.

Free Form Deformations
Free Form Deformations algorithm uses Thin Plate Splines [10,11,12].A thin-plate spline (TPS) belongs to the family of splines based on radial basis functions.In medical image registration, thin-plate splines are widely used, for example, as an approximation solution where the degree of approximation depends on the confidence of landmark localization.They have been formed by Duchon and Meinguet for the surface interpolation of scattered data [10,11].Radial basis function (RBF) spline can be defined as a linear combination of n radial basis functions ( ) s θ .
( ) The transformation can be defined as three separate thin plate spline functions T = (t 1 , t 2 , t 3 ) T (Eq.( 1)), which yield mapping between images where the coefficients correspond to a fine part of mapping of the spline-based transformation.
The goal of free-form deformations is to provide the means for modelling arbitrary deformations applied to objects.The general idea is to deform an image by manipulating a regular grid of control points that are distributed across the image.The control points can be moved, and the position of the individual pixels between control points is computed from positions of surrounding points.Compared to other type of splines (thin-plate spline, elastic-body spline), B-splines are locally controlled.This fact leads to computational efficiency for large number of points.In [3] was presented the framework for FFD registering breast images, which consists of combined local and global transformation (Eq.( 2)).
Combined transform is defined as the following ( , , ) ( , , ) ( , , ) FFD is deformed an object by manipulating an underlying mesh of control points.On the image volume, we define the mesh of control points with uniform spacing.FFD can be defined as a 3D tensor product of 1-D cubic B-splines: where: Here l B represents l th basis function of the B-spline (Fig. 3).
An advantage of FFD is a number of degrees of freedom, which has greater significance in non-rigid registration cases.Spacing of control points in the algorithm defines the number of degrees of freedom and computational complexity.Multi-resolution approach as described in Forsey [13] can be applied.Here we have hierarchy of control point meshes at different resolutions, whereas the resolution increases at each resolution level.This way, the sum of local transformations gives the local transformation.To obtain smooth transformations, this spline-based FFD is constrained with the regularizer, as defined in Wahba [14].The part of cost function which models TPS is defined with the following Therefore, the cost function is defined as the following 0 ( , ) ( ( ), ( ( ))) ( ) The smoothness part in Eq. ( 4) is the binding energy of thin-plate of metal.Minimizing the previous cost function, we find an "optimal" transformation.In the previous equation, the term   is the image similarity defined with mutual information, and the second term  ℎ corresponds to the cost associated with the smoothness of transformation.Finding optimum of Eq. ( 5) can be done with some of the well-known optimization methods, for example Gradient descent method.

DEMONS ALGORITHM
Thirion [3] has proposed the images to non-rigidly register by means of the diffusion process.In that paper he introduced the image entities (called demons) that push according to local characteristics of images in the similar way Maxwell did for solving the Gibbs paradox in thermodynamics.Basically, the task for demons is to sort the particles.This process moves particles outwards and inwards, which depends on the relation of scene and model.A good description of the Demons algorithm is given in Pennec [5].We can illustrate the idea for twodimensional case with two images S and M (Fig. 4 is differential force of the interaction between the static and moving image which is an 'external' force.Thirion has proposed iterative calculation of Eq. ( 6).The displacement is regularized in each iteration by Gaussian filter with a variance of 2 σ .
Here the regularization plays an important role in preserving the geometric continuity and suppressing the noise.From Eq. ( 6) it can be seen that information is driven only from static image by using its gradient.In [15] authors propose to expand the original method for calculating demons with another 'active' part or force m f  , which is based on the gradient information from the moving image By combining Eq. ( 6) and ( 7) we have the total force expressed as Here the influences of 'passive' and 'active' forces are mixed, whereas 'active' force is used to emphasize the contribution of the gradient information from moving image.Normalization factor α enables to adjust the force strength in each iteration, and smaller α can be used for relatively large deformations and then reduced when the algorithm approaches convergence.Then Eq.( 8) is given as follows (m ) (s ) (s ) Thirion's Demons method can be described with the flowchart (Fig. 5).The algorithm is iterative and runs for specified number of iterations.The problem with big deformations can be solved with multi-resolution approach, which drives the registration process from coarse to fine fashion.Depending on the deformation, four resolution levels are usually enough.Thirion's Demons uses the optical flow constant intensity assumption which is satisfied only for small deformations.This can be a problem when used in settings which include large deformations, for example in the case of adaptive radiation therapy.Big displacements break the intensity constancy assumption and lead to poor registration results.Multi-resolution approach can be applied for handling this problem, and leads the optimization problem into a global optimum.

RESULTS
The common starting point in the registration evaluation is to use synthetic images.Synthetic images have ground-truth information about deformation which is not available in many practical uses.This data can give us the overview about the registration quality, especially if we can simulate deformation which is physically valid for clinical use.
In this part we will show the performance of two algorithms in a case of motion discontinuities.In a case of non-rigid image registration, the validation of results is a challenging problem.Since registration is an ill-posed problem, there should be infinitely many mappings that can be a solution to the problem.Our interest is only in cases where it is physically consistent according to type of tissue involved.For this evaluation we start with a well-known method for the image registration: FFD B-spline algorithm as described in Rueckert [3].For all consequent simulations in this part, we use MATLAB environment for easy prototyping and comparisons between different algorithms.
For the quantitative analysis of the obtained deformation field was used Average Angular Error (AAE) [16].AAE can be calculated with the following: where   = ( 1  ,  2  ) denotes the ground truth deformation field.
The first pair of simulated images consists of three rectangular parts with two grey levels as fixed image (Fig. 6).In this image, the middle part is shifted for 5 pixels downwards and this is the floating image which serves us as ground-truth data.This experiment has a sliding part in the middle of image, where abruptly the change of image intensity in these regions causes discontinuous motion.6 shows the result of registration.In the first row, the pair of images to be registered is shown, and the ground truth motion field for each pair.Registration results for FFD algorithm is given in Fig. 6 in second row, deformed image and motion field, respectively.For each registration, Average Angle Error (AAE) was calculated and results are shown in Fig. 6 in the third row.
For the Demons algorithm, the following filter parameter was used: 2, 5, 10, 15, 25, 30, 35, 40 and 45 given in pixels.For each simulation other parameters were kept unchanged.Registration results are given in the fourth row of Fig. 6, while Average Angle Error (AAE) was calculated and results are shown in the sixth row in the same figure .From the results for block image it can be seen that both algorithms have problems with recovering the true motion.By comparing the registration results and the motion field we can see that FFD algorithm cannot reproduce the true motion in the middle rectangle.Demons can give better visual results, but the motion around the edges is not sharp and discontinuities are not well preserved.This is visible from the AAE, i.e.Demons gives lower value than FFD.For the other two experiments, discontinuity preserving is better for both algorithms.From the visual inspection and from the AAE it can be seen that Demons for this type of images produce better registrations.
It is visible from all experiments that increasing the resolution of control points corresponds to lower AAE, but not necessarily improved registration results in the region where discontinuities lie, and also introduce some errant registrations in the right part of image on the skull edge.Also lower control point resolutions in this example cannot truly align the right part of image to its correct position, and it can be seen that, for some pixels, it is lower than true position.Better deformation field pushes some regions in their correct position, but result is unacceptable, because the deformation on the edges and also wrong results are propagated on the left side where deformation should be zero.From color rendering of the motion field onto images one can see that the difference on the left side of images and in the region where discontinuity is located, is small but noticeable.
As we can see from the results for vertical wedge-like motion, the registered image is pretty well matched except for the left and right part of the image.The result for Demons is better, only the region with wedge undergoes the motion and small motion is visible near the discontinuities.

CONCLUSION
Image registration plays an important role in today's clinical practice.Recovering transformation which maps one image into the space of the second one is generally a tedious task.The need for a clinical procedure which includes the organ motion has lead to developing the registration algorithms which take into account when one organ slides along each other.This case needs special algorithms which incorporate discontinuities into the cost function.
Free Form Deformation and Demons algorithm were evaluated with synthetically deformed images where deformation is previously known.Average Angle Error was used as the evaluation measure.Results in this work show the performance of FFD and Demons algorithm for synthetic images with discontinuities.
From the AAE measure and also from visual inspection one can see that Demons algorithm performs better for this type of motion.One can see that beside good local support, FFD algorithm gives worse results than Demons for tested pairs of MR images, as well as synthetic images.On the other side, Demons has lower computational requirements than FFD and better performance in our experiments.Therefore, it is recommended to use Demons algorithm due to its speed and simplicity and also good results in most cases.

Figure 3
Figure 3 TPS neighbourhoods and four basis functions ).A demon d is located at spatial position d and ( ) 0 ≠ R d .According to the gradient R(d) and the image difference ( ) ( ) S d M d − , the demons produce the pushing force p. Template is pushed according to S(d), if ( ) ( ) S d M d < and -S(d), if ( ) ( ) S d M d > .Using gradient informationfor expressing the demons forces may not be efficient when gradient values are close to 0. For the estimation of demons forces, optical flow equation is used.For a given point P, let s be the intensity in static image S and m the intensity of the moving image M. Then the estimated displacement u required for the point is ∇ is a gradient of the static image.Here S ∇ represents 'internal' force originating from the static image.

Figure 4
Figure 4 An example of Demons 'forces'

Figure 5
Figure 5 Flowchart of Demons registration algorithm

Figure 6
Figure 6 Evaluation results for FFD and Demons registration.First row: image to be registered and corresponding ground truth flow; second row: registration with FFD algorithm, third row: change of AAE for FFD experiments; fourth row: registrations with Demons algorithm, fifth row: change of AAE for Demons experiments Evaluation of FFD registration will be carried for synthetic blocks and synthetic example of brain MR images.There are three synthetic examples: blocks, vertical motion and vertical shear motion.Moving image is obtained by shifting the corresponding part of fixed image 5 pixels downwards for each experiment.Then this image was registered with the original one using the FFD with different control points spacing, and Demons algorithm with different sigma values.For the evaluation we use the following control point neighbourhoods given in pixels 6×6, 9×9, 15×5, 26×26, 48×48, 92×92, for vertical shear motion and vertical motion, and additional for blocks 132×132 and 260×260, wherein other parameters are kept unchanged.

Fig.
Fig.6shows the result of registration.In the first row, the pair of images to be registered is shown, and the ground truth motion field for each pair.Registration results for FFD algorithm is given in Fig.6in second row, deformed image and motion field, respectively.For each registration, Average Angle Error (AAE) was calculated and results are shown in Fig.6in the third row.For the Demons algorithm, the following filter parameter was used: 2, 5, 10, 15, 25, 30, 35, 40 and 45 given in pixels.For each simulation other parameters were kept unchanged.Registration results are given in the fourth row of Fig.6, while Average Angle Error (AAE) was calculated and results are shown in the sixth row in the same figure.