Abstract
Dirichlet, anisotropic, and Huber regularization terms are presented for efficient registration of deformable images. Image registration, an ill-posed optimization problem, is solved using a gradient-descent-based method and some fundamental theorems in calculus of variations. Euler-Lagrange equations with homogeneous Neumann boundary conditions are obtained. These equations are discretized by multigrid and finite difference numerical techniques. The method is applied to the registration of brain MR images of size . Computational results indicate that the presented method is quite fast and efficient in the registration of deformable medical images.
1. Introduction
The purpose of image registration is to match two or more images of the same scene or object to each other, which can be obtained at different times, perspectives, and sensors such as MRI, X-ray, CT, PET, and tomography. In the present day, correct matching of similar images are quite useful and still a challenging problem. Although there exist some efficient image registration methods (for instance, [1–4]) in the literature, finding reliable and efficient image registration techniques are still significantly important and active research area among image processing problems.
Mathematically, variational techniques are based on minimizing an energy or cost function that usually consist of the sum of a data term and a regularization term (e.g., see [5]). In this paper, we express the image registration algorithm as a variational optimization problem, which consists of a sum of the similarity measure and a regularization term. Our method incorporates -norm sense similarity measures with several different regularization terms such as Dirichlet, anisotropic, and Huber.
Organization of the paper is as follows. In the first section, we present optimal control formulation of the image registration problem. We use the sum of squared differences in the -norm sense as the similarity measure and introduce three different regularization terms. We are mostly motivated by the efficient performance of these regularization terms in the statistical approach of image restoration problems (for instance, [6, 7] apply these regularization terms to image denoising problems). It is plausible to explore the utility of these statistically derived regularization terms on the deformation field. In the same section, we solve the existing optimization problem in a systematic way by the techniques in variational calculus such as gradient-descent-based methods and numerical methods such as finite differences. Computational results given in the last section indicate that new regularization terms provide fast, stable, and efficient image registration techniques.
2. Optimal Control Formulation of Image Registration Problem
The state of the art of the image registration problem can be expressed in the following way. Assume that both the template and reference images are defined on the same domain . Then, the image registration problem can be formulated as the optimization problem for the functional where denotes a similarity measure between the template image and the reference image , is the deformation field, is displacement field, is the set of all possible admissible transformations, is a regularization term, and is a regularization constant. Because reference and template images are obtained from different distances, angles, times, and sometimes even by different individuals, a deformation field may occur between these images. A deformation field is a vector field that maps pixels (or coordinates) of reference image to the corresponding ones of the template image. One of the major goals of this paper is to compute the deformation field in a systematic way.
We choose the -norm type similarity measure defined as This similarity measure is often referred to as the “sum of squared differences” (SSD) measure. Note that other similarity measures can be selected depending on the problem. We choose the similarity measure (2.3) due to its well-known effectiveness (e.g., see [8]), for the convenience in computations and for easily adapting the regularization terms in numerical solutions.
Without the regularizing term in functional (2.2), the image registration problem (2.1) is ill-posed; furthermore, imaging data usually is not smooth due to edges, folding, or other unwanted deformations. Ill-posed problems are widely used in PDE-based image processing problems and inverse problems. For instance, in [9] the authors solve ill-posed equations with transformed argument as an inverse problem. An optimization problem is said to be well-posed if the solution of the problem uniquely exists and the solution depends continuously on the data of the problem. If one of these three conditions are not satisfied, it is called an ill-posed problem. Image registration is an ill-posed optimal control problem. In order to overcome the ill-posedness of the optimization problem (2.1) and to assure smooth solutions, we introduce the additional regularization terms. The main idea behind adding a regularization term is to smoothen the problem with respect to both the functional and the solution so that well-posedness is assured and efficient computational methods can be defined to determine minimizers. Typical regularization terms associated with image registration problems include curvature, diffusion, elasticity, and fluid. Details about each of these regularization approaches can be seen, for example, in [1].
In this paper, we shall introduce several different types of regularization terms, , mostly in the form of where is a smooth function, is the gradient vector of , and is the -norm of . For the convenience of the computations, we define . Using formulas (2.3) and (2.4) in functional (2.2), we can express (2.2) as Next we review some fundamental mathematical lemmas and theorems that characterize the existence of a unique solution of a minimization problem to which we will fit our optimization problem.
Lemma 2.1. If a Gateux differentiable function, say , defined on an open set has a local extremum at , then .
The solutions of are known as Euler-Lagrange equations of the problem. If is convex, then a solution of is a solution of a minimization problem. This Lemma and both of the following Theorems are well known in calculus of variations (e.g., see [10] or [11]).
Theorem 2.2. Assume is a nonempty, closed, convex subset of a Hilbert space . Let be bilinear, symmetric, bounded, elliptic function, and with . Then, there exists a unique such that
It is clear from (2.6) that we can express the optimization problem given by (2.1) as a general regularization problem as: minimize the cost functional with
over a domain (). Next, we present a theorem for the solution of this problem.
Theorem 2.3. The solution of the general variational problem (2.7) is given as the solution of the Euler-Lagrange equations under the conditions stated in Theorem 2.2. denotes the (partial) derivative of with respect to . This solution can be interpreted as the steady-state solution to the following nonlinear elliptic PDE, called gradient descent flow: In both cases, appropriate boundary conditions must be applied.
Now, we solve the optimization problem (2.1) for Dirichlet, anisotropic, and Huber regularization terms. We start with Dirichlet variational integral For this regularization term, we can write functional (2.5) as Considering the functional (2.12) in the general functional (2.2) and by using Theorem 2.3, the solution of the optimization problem (2.1) is given by Discretizing (2.13) with a finite-difference method, we get where is computed with any interpolation method. Defining we obtain We solve this equation with a multigrid method. Now we express the implementation of this numerical solution: (1)set ;, (2)for to max, (i)compute ,(ii)update , ,(iii)compute ,(iv)solve ,(3)end. We implemented this algorithm in Matlab programming language where we chose the , , . Let us notice that because similar implementation procedures can be written for anisotropic and Huber-type regularization terms, we will omit their detail matlab programming implementations for the sake of brevity.
Using Theorem 2.3, the minimizer of (2.5) can be interpreted as the steady state solution to the following nonlinear elliptic PDE, called gradient descent flow: assuming homogeneous Neumann boundary conditions.
It was proved in [12] that the anisotropic diffusion PDE of Perona and Malik ([13]) given by is the gradient descent flow for the variational integral where is given by where . Hence, for these variational integrals we get the same equation with (2.17) and defined by (2.20).
Motivated by the robustness of the Huber M-filter in the probabilistic approach of image denoising ([6, 7]), we define the Huber variational integral as where
One may check that the Huber variational integral is well defined, convex, and coercive. It follows from Theorem 2.3 that the minimization problem has a unique solution when . Using the Euler-Lagrange variational principle, it follows that the Huber gradient descent flow is given by where is the Huber M-estimator weight function with homogeneous boundary conditions. Hence, for the Huber variational integral we can express (2.17) as where is defined in (2.25).
3. Computational Results
In this example, we demonstrate the registration of two-dimensional noise-free brain MR images in the size of . Computational examples indicate that all regularization terms produce similarly good registration quality but that the cost associated with Huber model approach is, on the average, less than that for others. The SSD registers images in an inverse way with its rate. In other words, the smaller SSD gives a better registration from the quality and efficiency point of views. Table 1 shows the changes in SSD for each model used versus number of iterations. Duration of the registration with the Huber model is almost 1 minute, and the one with anisotropic and Dirichlet is almost 2 minutes. The template, reference, and registered images with Huber regularization term are shown in Figure 1. Because the registered image is similar to the other regularization terms, we omit them inhere for the sake of brevity. We applied the presented method to some other brain MR images and obtained the similar results.
(a)
(b)
(c)
4. Concluding Remarks
Deformable image registration is a significant branch of image processing. It has a broad application in medical and nonmedical imaging. In this paper, we presented a number of deformable image registration techniques. Our method incorporates sum of squared differences similarity measure with several different regularization terms such as Dirichlet, anisotropic, and Huber. By variational calculus methods and gradient-descent-based optimization techniques, we solve the existing optimization problem and obtain the corresponding optimality system in a systematic manner. We were mostly motivated by the efficient performance of these regularization terms in the statistical approach of image denoising problems, (for instance, see [6, 7]).
In a future work, we will investigate the applications of these image registration techniques to the registration of noisy and blurred images. Furthermore, we plan to compare the strength of these registration techniques with some well-known image registration methods from speed, quality, and effectiveness point of view in detail.