1. Introduction
Image registration refers to the matching of coordinates of different images, which are obtained in different time periods or by different imaging techniques. As one of the most challenging tasks in image processing, image registration has been widely used in remote sensing data analysis, computer vision, medical imaging, and other fields [
1,
2,
3,
4,
5,
6,
7,
8]. Many image registration models have been proposed during recent decades, for example, the Demons model [
9], vectorial total variation model (VTV model) [
10], and fractional-order model [
11]. For more details, one can refer to [
12]. Mathematically, image registration is stated in the following way. Suppose that
is an open bounded set and
are two real functions defined on
. In image registration,
T and
D are viewed as two images, the source image and target image, respectively. The goal of registration is to find a geometric deformation
such that the deformed image
looks like
as much as possible. Concerning the problem of how to characterize the similarity between
and
, many metrics are proposed, such as the normalized cross correlation (NCC) [
13,
14], the normalized mutual information (NMI) [
14], and the sum of squared distance (SSD) [
15]. We mainly focus on mono-modality image registration. Therefore, SSD is used in this paper. That is, we aim to find the minimum of the SSD, which is defined by
where
and in small deformable image registration,
is approximately divided into trivial identity part
and displacement part
, i.e.,
. With these notations, the variational framework of image registration is formulated as:
where, in here and in what follows,
,
is an admissible solution set.
Equation (
2) is ill-posed [
16,
17]. A common way to overcome the ill-posedness is to add some regularization
[
11,
12,
13,
15,
16,
17,
18,
19,
20,
21,
22]. Based on this idea, the image registration model is reformulated as follows:
where
is a constant to balance
and
.
For different applications, different regularizations are introduced. For instance, in order to induce a diffusion type model,
is adopted in [
23]. In addition, for the purpose of fitting the fractional-order smoothness of image texture, some fractional-order regularizations [
24] are also introduced. For all these models, the diffeomorphic image registration models have become the focus of several groups [
11,
15,
19,
20,
24,
25,
26,
27,
28,
29]. In diffeomorphic registration, one is expected to produce a deformation
, which is a continuous bijection (1-to-1 and onto). By constraining the deformation into this kind of mapping, the topology of the object tissue is preserved. This is an advantage of diffeomorphic image registration. However, in some applications, i.e., the image pairs
T and
D with strong noises or pathological tissues, the registration deformation may be discontinuous. For example, in respiratory movement, the diaphragm is severely deformed but the chest remains almost rigid [
30], which leads to discontinuous deformation. This discontinuous deformation cannot be well simulated by diffeomorphic image registration models. To simulate discontinuous deformation, some other image registration models were also proposed [
11,
24]. In [
30], a registration model was introduced by letting
. This is called total variation (TV) registration. As we all know, TV preserves the edge of tissue well. However, serious staircase effect occurs in the TV model. To overcome the staircase effect of the TV model, some image registration models were also proposed for improvement, for instance, the bounded deformation (BD) model [
12]. In [
12], by viewing the image process as an elastic deformation, the authors propose a bounded deformation model, which addresses the discontinuous deformation well. However, the algorithm in [
12] is easily trapped into a local minimum of object functional
because it is non-convex on
. In fact, most of the published registration algorithms are all easily trapped into local minimum, which is caused by the non-convexity of the object functional. To overcome this problem, Han, Wang, and Zhang [
25] proposed a multiscale image registration approach for 2D diffeomorphic image registration. Theoretical results and numerical tests are also presented in [
25] to validate the fact that the multiscale approach can effectively overcome the local minimum of the object functional for image registration. Note that the multiscale approach in [
25] is only suitable for diffeomorphic image registration. As an extension, we propose a multiscale approach for BD model, which can overcome the local minimum of the object functional and simulate the discontinuous deformation well. The proposed multiscale approach is essentially an iterative process. The convergence of the iterative process is shown, and several numerical tests are performed to show the good performance of the proposed multiscale approach.
The rest of the paper is organized as follows. In
Section 2, we propose a multiscale image registration approach and show the convergence of the proposed approach. In
Section 3, an image registration algorithm of the proposed multiscale approach is presented. In
Section 4, several numerical tests are performed to show the good performance of the proposed multiscale approach.
2. Proposed Model
In [
12], authors proposed a BD image registration model, which is formulated as follows:
where
,
,
,
and
,
, there
. Here, the
norm is defined as follows:
.
The existence of the solution of (
4) is proved in [
12]. Moreover, an explicit scheme for (
4) is also presented. However, the proposed algorithm in [
12] is easily trapped into a local minimum of the object functional (
4), because the object functional in (
4) is non-convex on
. To overcome the local minimum of the object functional, we extend the diffeomorphism-based multiscale approach in [
25] into the BD model (
4) by setting an increasing parameter sequence
in each iteration. The proposed multiscale approach for BD model is listed as follows:
Step 1: Solving the following variational problem to got an initial value:
where,
and
. Based on (
5), we defined
.
Step 2: Based on the deformation
in Step 1, we updated the source image by
and solved the following variational problem to get a new displacement filed
:
where
. Based on the result of (
6), we defined
and
.
Step n: We updated the source image by
and solved the following variational problem:
where
. Based on (
7), we defined
and
. We repeated this process until it reached the desired accuracy.
Remark 1. (5)–(7) is a process of function composition, which we can write uniformly as follows:where , , Note that we set to be an increasing sequence in (8) for the following two reasons: It increases the ratio of the similarity term in the object functional. Benefiting from this setting, plays a dominant role in object functional as n is large enough. This improves the accuracy of image registration.
An increasing helps the algorithm get out of the local minimum of the object functional in former iterative steps.
For the selection of in the actual calculation, we set it as an exponential function, i.e., , where is the maximum scale number. The registration result is optimized by adjusting the values of parameters a and .
Theorem 1. Define and set to be an upper bounded sequence with as n is large enough, then the deformation induced by (5)–(7) converges to the global minimizer of on . Proof. By the fact , we know that for any .
This implies,
is a decreasing sequence with lower bound. By the monotone bounded convergence theorem [
31], we know that the sequence
is convergent.
On the other hand, by (
9) we obtain that
By (
11) and (
12), we get that
.
That is,
converges to the global minimizer of
on
. Moreover, we know that
is an approximate minimizer of
on
because
Note that here we use the condition as n is large enough. □
Remark 2. One can also see that the similarity sequence is independent of these parameters in (5)–(7). This implies that the final registration result does not rely on parameters. This ensures the robustness of the proposed multiscale approach. 3. Numerical Implementation of the Proposed Multiscale Approach
In this section, we discuss the numerical implementation of the proposed multiscale approach (
5)–(
7). In the discussion in
Section 2, one can notice that the iterative process in (
5)–(
7) is uniformly rewritten as the variational problem (
8). That is, it is expected to solve the variational problem (
8) for each
n in the proposed iterative process. Now, we focus on the numerical implementation of the variational problem (
8).
Firstly, we hope to transform the variational problem (
8) into some partial differential equations (PDE). For this purpose, we choose a perturbation
along the minimizer
. This induces a new function
. Then, we obtain that
By fixing and , is then a function on t, which is denoted by .
Note that
is a minimizer of (
8); therefore, one can obtain that
. That is
for arbitrary
, then
. Based on this notation, we know that
Due to
being arbitrary, and the variational principle in [
32], we can obtain the Euler-Lagrange equation:
Equation (
16) is a nonlinear elliptic PDE. It is difficult to solve (
16) directly. Motivated by [
18], here, we use the gradient descent flow method to solve (
16). By introducing a new variable
, the gradient descent flow equation for (
16) is formulated as follows:
Using the similar method of [
18], one can prove that (
16) is the steady state of (
17). That is, (
16) and (
17) are equivalent as
t is large enough.
For numerical implementation, we discretize the region as follows: is assumed to be a square, i.e., , for , so we define , , for .
By using the finite difference method,
and
can be approximated by the following formulas:
, where is the time step.
In a similar way, we obtain that
where
With these approximations, (
17) is approximated by the following equations:
where
and
For simplicity, we set
in (
23). The numerical implementation for variational problems (
8) is summarized as Algorithm 1.
Algorithm 1 Gradient descent algorithm |
Initialization: the source image T, target image D, , , , and maximum iteration times N; Whiledo 1. Calculate and by ( 18), ( 24), ( 25); 2. Update the displacement filed as: by ( 23); 3. Calculate using interpolation; 4. Set and repeat 2–4 until the maximum iteration times is reached. Output: for some . |
Based on Algorithm 1, the numerical implementation for (
5)–(
8) is stated as Algorithm 2.
Algorithm 2 Multiscale gradient descent algorithm of bounded deformation function |
Initialization: Given maximum scale number and ; For n do 1. Let and update ; 2. Use Algorithm 1 to obtain for ; 3. Compute , where and . Output: . |
4. Numerical Tests
In this section, we use several numerical tests to demonstrate the good performance of the proposed multiscale approach. All the numerical tests were implemented under Windows 10 and Matlab R2019a with 11th Gen Intel(R) Core(TM) i5-1135G7 @ 2.40 GHz and 16.0 GB memory. For this purpose, we compared the proposed multiscale approach with the BD model in [
12], LDDMM-demons algorithm in [
33], and TV model in [
30] to validate the efficiency of the multiscale approach. In particular, the BD model can be viewed as a single scale approach without iteration (scale = 1). The used data sets included synthetic images, medical images, and underwater distorted images. The comparison between the BD model and the proposed multiscale approach helps to validate the fact that the proposed approach has advantages in getting out of the local minimum of the object functional. The size of all images used for registration is
. In order to evaluate the similarity between the deformed source image and target image, we used the relative sum of squared difference to quantitatively evaluate the registration results. The relative sum of squared difference is defined as follows:
where
.
Test 1 (Test for synthetic images) In this test, we presented two numerical tests by comparing these four algorithms: proposed multiscale approach (Algorithm 2), BD model algorithm in [
12], the LDDMM-demons algorithm in [
33], and TV model in [
30]. These two tests were performed on synthetic image pairs, Rectangle and C-E, respectively. By using these four different algorithms to register the above two image pairs, respectively, we list the registration results in
Figure 1,
Figure 2,
Figure 3 and
Figure 4 and
Table 1. Moreover, sensitivity tests for parameter
are also listed in
Figure 1. One can notice that from
Figure 1f, the final registration result is not sensitive to the parameter, while
and
only affects the convergence speed. This coincides with the result in [
25]. In
Figure 2 and
Figure 4, we can see that
gradually decreases with respect to the scale number
n. In addition, the final
achieves 1.3% and 0.27%, respectively. This indicates
, which implies that the proposed multiscale approach leads to an accurate result. With the quantitative comparison results in
Table 1, we obtain the following two conclusions: (I) The proposed multiscale approach has advantages in overcoming the local minimum of the object functional. In fact, in the comparison between the BD model and the proposed algorithm, we see that the proposed Algorithm 2 obviously improves the registration result. Moreover, one can notice that the final
for the BD model on the Rectangle image pair is 71%. This shows that the algorithm in [
12] is trapped into a local minimum. Compared with
led by the local minimum, the proposed Algorithm 2 achieves a very good result with 1.3%. This validates the fact that the proposed Algorithm 2 has advantages on overcoming the local minimum of the object functional. This is also the main motivation for this work. In addition, the final
of image pair C-E(
) shows that the proposed Algorithm 2 can also address the smooth registration well (note that C-E induces a continuous deformation). At last, with
Figure 1e and
Figure 3e, we can conclude the sequence
converges as
. This validates that the proposed Algorithm 2 induces a convergence similarity sequence. This coincides with the conclusion at the end of
Section 2. (II) The proposed Algorithm 2 addresses the discontinuous registration well. One can notice that the image pair Rectangle induces a discontinuous deformation (this is the main reason why this image pair is selected for numerical comparison). It follows from
Table 1 that the three discontinuous methods (Algorithm 2, BD model, TV model) obviously perform much better than the diffeomorphic model (LDDMM). Moreover, one can also see that the proposed Algorithm 2 achieves the smallest
. This shows that the proposed Algorithm 2 is competitive to some state-of-the-art registration algorithms.
Test 2 (Test for medical images) Similar to Test 1, we used Algorithm 2, the BD model algorithm in [
12], the LDDMM-demons algorithm in [
33], and the TV model in [
30] to register three different kinds of medical image pairs: lung, hand, and liver. The final registration results are listed in
Figure 5,
Figure 6,
Figure 7,
Figure 8,
Figure 9 and
Figure 10 and
Table 2. With
Figure 5e,
Figure 7e and
Figure 9e, we conclude that the proposed multiscale approach plays an important role in overcoming the local minimum. On the other hand, in
Figure 5d,
Figure 7d and
Figure 9d, we see that the final registration result
looks nearly the same with the target image
. This shows the high accuracy of the proposed Algorithm 2 in view of computer vision. Concerning the quantitative comparison between the BD model and the proposed multiscale approach, it follows from
Table 2 that the proposed multiscale approach obviously improves the result of the BD model, though it costs much more CPU time. Note that here we care more about the registration result rather than the CPU time. Particularly, the final result (
) of the BD model demonstrates that the BD model gets into a local minimum of the object functional. This is mainly caused by the nonconvexity of the object functional. A sharp contrast to this result, the final registration result (
) for the proposed Algorithm 2 shows that the proposed approach can effectively improve the final registration results. Concerning the process for the details of improvement, one can see
Figure 6,
Figure 8 and
Figure 10. In addition, through the quantitative comparison results in
Table 2, we see that the proposed Algorithm 2 is competitive to the other three registration algorithms. Moreover, one can see that the BD model achieves the worst result among these four algorithms, while the BD-based multiscale approach achieves the best image registration result. This further validates the fact that the proposed approach has advantages in overcoming the local minimum of the object functional.
Test 3 (Test for underwater images) At the end of this section, we test the proposed Algorithm 2 on very special natural image pairs: underwater images. In underwater images, serious distortion occurs, which makes it difficult for image registration. In this test, we selected two image pairs (Square and Coin) as experimental data. We used Algorithm 2, the BD model algorithm in [
12], the LDDMM-demons algorithm in [
33], and the TV model in [
30] to register these two image pairs, and the final registration results are listed in
Figure 11,
Figure 12,
Figure 13 and
Figure 14 and
Table 3. In
Figure 11e and
Figure 13e, we see that the
induced by proposed the multiscale approach converges to the global minimum of
on
. Additionally, with
Figure 11d and
Figure 13d, we know that the deformed image
induced by the proposed Algorithm 2 looks the same with target image
. At last, one can see from
Table 3 that the proposed multiscale approach achieves the best result among the four comparison results. This shows the high accuracy of the proposed algorithm. Moreover, via comparison through data Coin, we see that the LDDMM and TV model are trapped in the local minimum for some default parameters, while the local minimum does not occur in all the numerical tests performed by the proposed Algorithm 2. This validates the fact that the proposed Algorithm 2 is robust and independent of parameters.
Analysing the registration results of these three tests, we conclude that the proposed Algorithm 2 has the advantage of improving the registration results for different kinds of image pairs (i.e., natural iamge and medical image). Moreover, via the relationship between and scale number listed in these three tests, we also conclude that the proposed Algorithm 2 can help to overcome the local minimum of the object functional and make the registration more accurate.