Next Article in Journal
Optimal PID Controllers for AVR Systems Using Hybrid Simulated Annealing and Gorilla Troops Optimization
Next Article in Special Issue
Characteristic Analysis and Circuit Implementation of a Novel Fractional-Order Memristor-Based Clamping Voltage Drift
Previous Article in Journal
Fractal Perturbation of the Nadaraya–Watson Estimator
Previous Article in Special Issue
Dynamical Analysis of a Novel Fractional-Order Chaotic System Based on Memcapacitor and Meminductor
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Multiscale Approach for Bounded Deformation Image Registration

Department of Mathematics, Wuhan University of Technology, Wuhan 430070, China
*
Author to whom correspondence should be addressed.
Fractal Fract. 2022, 6(11), 681; https://doi.org/10.3390/fractalfract6110681
Submission received: 8 August 2022 / Revised: 15 November 2022 / Accepted: 16 November 2022 / Published: 18 November 2022

Abstract

:
Deformable image registration is a very important topic in the field of image processing. It is widely used in image fusion and shape analysis. Generally speaking, image registration models can be divided into two categories: smooth registration and non-smooth registration. During the last decades, many smooth registration models (i.e., diffeomorphic registration) were proposed. However, image with strong noise may lead to discontinuous deformation, which cannot be modelled by smooth registration. To simulate this kind of deformation, some non-smooth registration models were also proposed. However, numerical algorithms for these models are easily trapped into a local minimum because of the nonconvexity of the object functional. To overcome the local minimum of the object functional, we propose a multiscale approach for a non-smooth registration model: the bounded deformation (BD) model. The convergence of the approach is shown, and numerical tests are also performed to show the good performance of the proposed multiscale approach.

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 Ω R 2 is an open bounded set and T , D : Ω R 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 h : Ω Ω such that the deformed image T h ( · ) looks like D ( · ) as much as possible. Concerning the problem of how to characterize the similarity between T h ( · ) and D ( · ) , 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
Ω [ T ( h ( x ) ) D ( x ) ] 2 d x .
where x = x 1 , x 2 and in small deformable image registration, h ( x ) is approximately divided into trivial identity part x and displacement part u ( x ) = u 1 ( x ) , u 2 ( x ) , i.e., h ( x ) x + u ( x ) . With these notations, the variational framework of image registration is formulated as:
u = arg min u K S ( u ) ,
where, in here and in what follows, S ( u ) = Ω [ T ( x + u ( x ) ) D ( x ) ] 2 d x , K 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 R ( u )  [11,12,13,15,16,17,18,19,20,21,22]. Based on this idea, the image registration model is reformulated as follows:
u = arg min u K S ( u ) + λ R ( u ) ,
where λ > 0 is a constant to balance R ( u ) and S ( u ) .
For different applications, different regularizations are introduced. For instance, in order to induce a diffusion type model, R ( u ) = u L 2 ( Ω ) 2 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 h : Ω Ω , 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 R ( u ) = u L 1 ( Ω ) . 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 S ( u ) + λ R ( u ) because it is non-convex on u . 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:
u = arg min u K Θ S ( u ) + λ R ( u ) ,
where Θ , λ > 0 , R ( u ) = Ω u + T u d x , K = u = ( u 1 , u 2 ) B D ( Ω ) | u | Ω = 0 , Ω = ( a 1 , b 1 ) × ( a 1 , b 1 ) and B D ( Ω ) = u L 1 Ω ; R 2 | ε ( u ) ( Ω ) < , ε ( u ) ( Ω ) = Ω u + T u d x , there L 1 ( Ω ; R 2 ) = u : Ω R 2 | Ω u ( x ) d x < . Here, the L 1 norm is defined as follows: u 1 = Ω u d x = Ω i = 1 2 u i 2 d x .
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 u . 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 Θ n 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:
    u 0 = arg min u K Θ 0 Ω T h ( x ) D ( x ) 2 d x + λ R ( u ) ,
    where, h ( x ) = x + u ( x ) and Θ 0 > 0 . Based on (5), we defined h ^ 0 ( x ) = x + u 0 ( x ) .
  • Step 2: Based on the deformation h ^ 0 ( x ) in Step 1, we updated the source image by T h ^ 0 ( · ) and solved the following variational problem to get a new displacement filed u 1 :
    u 1 = arg min u K Θ 1 Ω T h ^ 0 ( x + u ( x ) ) D ( x ) 2 d x + λ R ( u ) ,
    where Θ 1 > 0 . Based on the result of (6), we defined h 1 ( x ) = x + u 1 ( x ) and h ^ 1 ( x ) = h ^ 0 h 1 ( x ) .
  • Step n: We updated the source image by T h ^ n 1 ( · ) and solved the following variational problem:
    u n = arg min u K Θ n Ω T h ^ n 1 ( x + u ( x ) ) D ( x ) 2 d x + λ R ( u ) ,
    where Θ n > 0 . Based on (7), we defined h n ( x ) = x + u n ( x ) and h ^ n ( x ) = h ^ n 1 h n ( x ) . 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:
u n = arg min u K E n ( u ) n N .
where T n ( x ) = T h ^ n 1 ( x ) , E n ( u ) = S n ( u ) + λ R ( u ) , S n ( u ) = Θ n Ω T n ( x + u ( x ) ) D ( x ) 2 d x .
Note that we set Θ n to be an increasing sequence in (8) for the following two reasons:
  • It increases the ratio of the similarity term S n ( u ) in the object functional. Benefiting from this setting, S n ( u ) plays a dominant role in object functional as n is large enough. This improves the accuracy of image registration.
  • An increasing Θ n helps the algorithm get out of the local minimum of the object functional in former iterative steps.
For the selection of Θ n n = 1 K M in the actual calculation, we set it as an exponential function, i.e., Θ n = a n × Θ 0 ( n N + ) , where K M is the maximum scale number. The registration result is optimized by adjusting the values of parameters a and Θ 0 .
Theorem 1.
Define K ¯ = h = h 0 h 1 h n h i ( x ) = x + u i ( x ) a n d u i K and set Θ n to be an upper bounded sequence with Θ n λ as n is large enough, then the deformation h ^ n · induced by (5)–(7) converges to the global minimizer of S ( u ) on K ¯ .
Proof. 
By the fact 0 K , we know that E n ( u n ) E n ( 0 ) for any n N + .
That is,
E n ( u n ) = Θ n T h ^ n ( · ) D ( · ) L 2 ( Ω ) 2 + λ R ( u n ) Θ n T h ^ n 1 ( · ) D ( · ) L 2 ( Ω ) 2 = E n ( 0 ) .
Therefore,
T h ^ n ( · ) D ( · ) L 2 ( Ω ) 2 T h ^ n 1 ( · ) D ( · ) L 2 ( Ω ) 2 n N + .
This implies, T h ^ n ( · ) D ( · ) L 2 ( Ω ) 2 is a decreasing sequence with lower bound. By the monotone bounded convergence theorem [31], we know that the sequence T h ^ n ( · ) D ( · ) L 2 ( Ω ) 2 is convergent.
Define
lim n T h ^ n ( · ) D ( · ) L 2 ( Ω ) 2 = δ R .
On the other hand, by (9) we obtain that
T h ^ n ( · ) D ( · ) L 2 ( Ω ) 2 + λ Θ n R ( u n ) T h ^ n 1 ( · ) D ( · ) L 2 ( Ω ) 2 .
By (11) and (12), we get that R ( u n ) n 0 .
That is, h ^ n ( · ) converges to the global minimizer of E n ( u ) on K ¯ . Moreover, we know that h ^ ( · ) is an approximate minimizer of S ( u ) on K ¯ because
min u K ¯ 1 Θ n E n ( u ) min u K ¯ S ( u ) .
Note that here we use the condition Θ n λ as n is large enough. □
Remark 2.
One can also see that the similarity sequence T h ^ n ( · ) D ( · ) L 2 ( Ω ) 2 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 η C 0 ( Ω ; R 2 ) along the minimizer u n . This induces a new function u ¯ = u n + t η K . Then, we obtain that
E ( u n + t η ) = Θ n Ω T n ( x + u n ( x ) + t η ) D ( x ) 2 d x + λ Ω u n + t η + T u n + t η d x .
By fixing u n and η , E ( u n + t η ) is then a function on t, which is denoted by φ ( t ) = E ( u n + t η ) .
Note that u n is a minimizer of (8); therefore, one can obtain that E ( u n ) E ( u ¯ ) . That is φ ( 0 ) φ ( t ) for arbitrary t R , then d φ ( t ) d t | t = 0 = 0 . Based on this notation, we know that
d E ( u n + t η ) d t | t = 0 = 2 Θ n Ω ( T n ( x + u n ( x ) ) D ( x ) ) T n ( x + u n ( x ) ) η d x + λ Ω u n + T u n u n + T u n η + T η d x = 2 Θ n Ω ( T n ( x + u n ( x ) ) D ( x ) ) T n ( x + u n ( x ) ) η d x 2 λ Ω div u n + T u n u n + T u n η d x = Ω [ 2 Θ n ( T n ( x + u n ( x ) ) D ( x ) ) T n ( x + u n ( x ) ) 2 λ div u n + T u n u n + T u n ] η d x = 0 .
Due to η C 0 ( Ω ; R 2 ) being arbitrary, and the variational principle in [32], we can obtain the Euler-Lagrange equation:
0 = Θ n ( T n ( x + u n ( x ) ) D ( x ) ) T n ( x + u n ( x ) ) λ div u n + T u n u n + T u n .
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 t > 0 , the gradient descent flow equation for (16) is formulated as follows:
u n t = Θ n ( T n ( x + u n ( x ) ) D ( x ) ) T n ( x + u n ( x ) ) + λ div u n + T u n u n + T u n .
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., Ω ( a 1 , b 1 ) × ( a 1 , b 1 ) , for N N + , so we define h = b 1 a 1 N , x i , j = ( x 1 , i , x 2 , j ) , x i x 1 , i = i h , x j x 2 , j = j h for i , j = 0 , 1 , 2 , , N .
By using the finite difference method, u n and u n t can be approximated by the following formulas:
( u n ) i , j = u 1 n x 1 u 1 n x 2 u 2 n x 1 u 2 n x 2 i , j 1 h ( u 1 n ) i + 1 , j ( u 1 n ) i , j ( u 1 n ) i , j + 1 ( u 1 n ) i , j ( u 2 n ) i + 1 , j ( u 2 n ) i , j ( u 2 n ) i , j + 1 ( u 2 n ) i , j ,
u n t ( u n ) t + 1 ( u n ) t δ t , where δ t is the time step.
In a similar way, we obtain that
( u n + T u n ) i , j 2 u 1 n x 1 u 1 n x 2 + u 2 n x 1 u 2 n x 1 + u 1 n x 2 2 u 2 n x 2 i , j ,
m i , j u n + T u n i , j = 4 u 1 n x 1 i , j 2 + 2 u 1 n x 2 + u 2 n x 1 i , j 2 + 4 u 2 n x 2 i , j 2 ,
u n + T u n u n + T u n i , j = ( u n + T u n ) i , j u n + T u n i , j , i f u n + T u n 0 0 , i f u n + T u n = 0 ,
div u n + T u n u n + T u n i , j = ( ( div 1 ) i , j , ( div 2 ) i , j ) T ,
where ( div 1 ) i , j = ϕ 11 x 1 i , j + ϕ 12 x 2 i , j ( ϕ 11 ) i + 1 , j ( ϕ 11 ) i , j h + ( ϕ 12 ) i , j + 1 ( ϕ 12 ) i , j h , ( div 2 ) i , j = ϕ 21 x 1 i , j + ϕ 22 x 2 i , j ( ϕ 21 ) i + 1 , j ( ϕ 21 ) i , j h + ( ϕ 22 ) i , j + 1 ( ϕ 22 ) i , j h ;   ( ϕ 11 ) i , j 2 m i , j u 1 n x 1 i , j , ( ϕ 12 ) i , j = ( ϕ 21 ) i , j ( u 1 n x 2 ) i , j + ( u 2 n x 1 ) i , j m i , j , ( ϕ 22 ) i , j 2 m i , j u 2 n x 2 i , j .
With these approximations, (17) is approximated by the following equations:
u n i , j ( t + 1 ) = u n i , j ( t ) + δ u n i , j ( t ) ,
where
δ u n i , j ( t ) = δ t [ λ div u n + T u n u n + T u n i , j ( t ) Θ n T n ( x + u n ( x ) ) i , j ( t ) D i , j T n ( x + u n ( x ) ) i , j ( t ) ] ,
and
T ( t ) i , j = T x ( t ) , T y ( t ) i , j T = ( T ( t ) ) i + 1 , j ( T ( t ) ) i , j , ( T ( t ) ) i , j + 1 ( T ( t ) ) i , j T .
For simplicity, we set h = 1 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, Θ n , λ , δ t , and maximum iteration times N;
While t N do
1. Calculate ( u n ) ( t ) , δ ( u n ) ( t ) and T ( t ) by (18), (24), (25);
2. Update the displacement filed as: ( u n ) ( t + 1 ) = ( u n ) ( t ) + δ ( u n ) ( t ) by (23);
3. Calculate T ( x + ( u n ) ( t ) ( x ) ) using interpolation;
4. Set t = t + 1 and repeat 2–4 until the maximum iteration times is reached.
Output: T x + ( u n ) ( t ) ( x ) for some t N .
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 K M and λ ;
For n = 0 : K M do
1. Let u = 0 and update Θ n ;
2. Use Algorithm 1 to obtain ( u n ) i , j for i , j = 1 , 2 , , N 1 ;
3. Compute T n ( · ) = T h ^ n ( · ) , where h ^ n = h ^ n 1 h n and ( h n ) i , j = x i , j + ( u n ) i , j .
Output: T h ^ K M ( · ) .

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 129 × 129 . 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:
R e S S D ( T , D , u ) = S S D ( T ( x + u ) , D ) S S D ( T , D ) ,
where S S D ( T , D ) = 1 2 i , j ( T i , j D i , j ) 2 .
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 Θ n are also listed in Figure 1. One can notice that from Figure 1f, the final registration result is not sensitive to the parameter, while Θ n and Θ n only affects the convergence speed. This coincides with the result in [25]. In Figure 2 and Figure 4, we can see that S n ( u ) gradually decreases with respect to the scale number n. In addition, the final R e S S D achieves 1.3% and 0.27%, respectively. This indicates T h ^ K M ( · ) D ( · ) , 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 R e S S D 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 R e S S D = 71 % 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 R e S S D of image pair C-E( R e S S D = 0.27 % ) 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 S n ( u ) converges as n . 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 R e S S D . 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 T h ^ K M ( · ) looks nearly the same with the target image D ( · ) . 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 ( R e S S D = 16.85 % ) 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 ( R e S S D = 5.88 % ) 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 S n ( u ) induced by proposed the multiscale approach converges to the global minimum of S ( u ) on K . Additionally, with Figure 11d and Figure 13d, we know that the deformed image T h ^ K M ( · ) induced by the proposed Algorithm 2 looks the same with target image D ( · ) . 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 R e S S D 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.

5. Conclusions and Discussion

In this paper, we propose a multiscale approach for the BD model, which addresses the discontinuous registration well. By setting an increasing sequence Θ n in the multiscale approach, we show that the similarity sequence S n ( u ) converges to the global minimum of the original similarity. A numerical algorithm is also presented to simulate the multiscale approach. Based on the proposed algorithm, several numerical tests are performed to validate the main motivation for this paper: (I) Overcoming the local minimum of the similarity; (II) Addressing the discontinuous registration caused by strong noise or pathological tissue.
For future research, one could try to extend the 2D image registration into 3D multiscale image registration. The study of 3D image registration is a very important topic. This makes it possible for clinical application and diagnosis.

Author Contributions

H.H. contributed to the conceptualization and methodology of the study. Y.D. and H.H. performed the formal analysis and investigation. Y.D. wrote the manuscript. Y.D. and H.H. performed the software and revising. All authors have read and agreed to the published version of the manuscript.

Funding

This work was supported by National Key Research and Development Program of China (No. 2020YFA0714200) and National Natural Science Foundation of China (No. 11771127, No. 11871386, No. 11871387, No. 11931012, and No. 11901443).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data are available from the corresponding author upon request.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Brown, L.G. A survey of image registration techniques. ACM Comput. Surv. 1992, 24, 325–376. [Google Scholar] [CrossRef]
  2. Lester, H.; Arridge, S.R. A survey of hierarchical non-linear medical image registration. Pattern Recognit. 1999, 32, 129–149. [Google Scholar] [CrossRef]
  3. Maintz, J.B.A.; Viergever, M.A. A survey of medical image registration. Med. Image Anal. 1998, 2, 1–36. [Google Scholar] [CrossRef]
  4. Mohamed, A.; Zacharaki, E.I.; Shen, D.G.; Davatzikos, C. Deformable registration of brain tumor images via a statistical model of tumor-induced deformation. Med. Image Anal. 2006, 10, 752–763. [Google Scholar] [CrossRef]
  5. Sotiras, A.; Davatzikos, C.; Paragios, N. Deformable medical image registration: A survey. IEEE Trans. Med. Imaging 2013, 32, 1153–1190. [Google Scholar] [CrossRef] [Green Version]
  6. Zitova, B.; Flusser, J. Image registration methods: A survey. Image Vis. Comput. 2003, 21, 977–1000. [Google Scholar] [CrossRef] [Green Version]
  7. Quan, D.; Wang, S.; Gu, Y.; Gu, Y.; Lei, R.Q. Deep feature correlation learning for multi-modal remote sensing image registration. IEEE Trans. Geosci. Remote Sens. 2022, 60, 1–16. [Google Scholar] [CrossRef]
  8. Zheng, Z.Y.; Cao, W.M.; Duan, Y.; Cao, G.T.; Lian, D.L. Multi-strategy mutual learning network for deformable medical image registration. Neurocomputing 2022, 501, 102–112. [Google Scholar] [CrossRef]
  9. Thirion, J.P. Image matching as a diffusion process: An analogy with Maxwell’s demons. Med. Image Anal. 1998, 2, 243–260. [Google Scholar] [CrossRef] [Green Version]
  10. Chumchob, N. Vectorial total variation-based regularization for variational image registration. IEEE Trans. Image Process. 2013, 22, 4551–4559. [Google Scholar] [CrossRef] [PubMed]
  11. Zhang, J.P.; Chen, K. Variational image registration by a total fractional-order variation model. J. Comput. Phys. 2015, 293, 442–461. [Google Scholar] [CrossRef] [Green Version]
  12. Nie, Z.W.; Yang, X.P. Deformable image registration using functions of bounded deformation. IEEE Trans. Med. Imaging 2019, 38, 1488–1500. [Google Scholar] [CrossRef]
  13. Modersitzki, J. FAIR: Flexible Algorithms for Image Registration; SIAM: Philadelphia, PA, USA, 2009; pp. 171–191. [Google Scholar]
  14. Sengupta, D.; Gupta, P.; Biswas, A. A survey on mutual information based medical image registration algorithms. Neurocomputing 2022, 486, 174–188. [Google Scholar] [CrossRef]
  15. Zhang, D.P.; Chen, K. A novel diffeomorphic model for image registration and its algorithm. J. Math. Imaging Vis. 2018, 60, 1261–1283. [Google Scholar] [CrossRef] [Green Version]
  16. Amit, Y. A nonlinear variational problem for image matching. SIAM J. Sci. Comput. 1994, 15, 207–224. [Google Scholar] [CrossRef] [Green Version]
  17. Christensen, G.E.; Rabbitt, R.D.; Miller, M.I. Deformable templates using large deformation kinematics. IEEE Trans. Image Process. 1996, 5, 1435–1447. [Google Scholar] [CrossRef]
  18. Han, H.; Wang, A. A fast multi grid algorithm for 2D diffeomorphic image registration model. J. Comput. Appl. Math. 2021, 394, 113576. [Google Scholar] [CrossRef]
  19. Lui, L.M.; Ng, T.C. A splitting method for diffeomorphism optimization problem using Beltrami coefficients. J. Sci. Comput. 2014, 63, 573–611. [Google Scholar] [CrossRef]
  20. Lui, L.M.; Wen, C. Geometric registration of high-genus surfaces. SIAM J. Imaging Sci. 2013, 7, 337–365. [Google Scholar] [CrossRef] [Green Version]
  21. Vercauteren, T.; Pennec, X.; Perchant, A. Non-parametric diffeomorphic image registration with the demons algorithm. In Proceedings of the International Conference on Medical Image Computing and Computer-Assisted Intervention, Brisbane, Australia, 29 October–2 November 2007; pp. 319–326. [Google Scholar]
  22. Zhang, D.P.; Chen, K. 3D orientation-preserving variational models for accurate image registration. SIAM J. Imaging Sci. 2020, 13, 1653–1691. [Google Scholar] [CrossRef]
  23. Chumchob, N.; Chen, K. A robust multigrid approach for variational image registration models. J. Comput. Appl. Math. 2011, 236, 653–674. [Google Scholar] [CrossRef] [Green Version]
  24. Han, H.; Wang, Z.P. A diffeomorphic image registration model with fractional-order regularization and Cauchy–Riemann constraint. SIAM J. Imaging Sci. 2020, 13, 1240–1271. [Google Scholar] [CrossRef]
  25. Han, H.; Wang, Z.P.; Zhang, Y.M. Multiscale approach for two-dimensional diffeomorphic image registration. Multiscale Model. Simul. 2021, 19, 1538–1572. [Google Scholar] [CrossRef]
  26. Modin, K.; Nachman, A.; Rondi, L. A multiscale theory for image registration and nonlinear inverse problems. Adv. Math. 2019, 346, 1009–1066. [Google Scholar] [CrossRef] [Green Version]
  27. Han, H.; Wang, Z.P.; Zhang, Y.M. Multiscale approach for three-dimensional conformal image registration. SIAM J. Imaging Sci. 2022, 15, 1431–1468. [Google Scholar] [CrossRef]
  28. Zhang, J.P.; Li, Y.Y. Diffeomorphic image registration with an optimal control relaxation and its implementation. SIAM J. Imaging Sci. 2021, 14, 1890–1931. [Google Scholar] [CrossRef]
  29. Cai, C.W.; Wang, L.; Ying, S.H. Symmetric diffeomorphic image registration with multi-label segmentation masks. Mathematics 2022, 10, 1946. [Google Scholar] [CrossRef]
  30. Pock, T.; Urschler, M.; Zach, C. A duality based algorithm for TV-L1-optical-flow image registration. In Proceedings of the International Conference on Medical Image Computing and Computer-Assisted Intervention, Brisbane, Australia, 29 October–2 November 2007; pp. 511–518. [Google Scholar]
  31. Zorich, V.A.; Paniagua, O. Mathematical Analysis II; Springer: Berlin/Heidelberg, Germany, 2016. [Google Scholar]
  32. Evans, L.C. Partial Differential Equations, 2nd ed.; AMS: Providence, RI, USA, 2010. [Google Scholar]
  33. Li, J.; Shi, Y.; Tran, G.; Dinov, I.; Wang, D.; Toga, A. Fast local trust region technique for diffusion tensor registration using exact reorientation and regularization. IEEE Trans. Med. Imaging 2013, 33, 1005–1022. [Google Scholar] [CrossRef]
Figure 1. Registration result of Algorithm 2 on image pairs Rectangle. (a) T ( · ) . (b) D ( · ) . (c) | T D | . (d) T h ^ K M ( · ) . (e) R e S S D change (Rectangle) with scale number. (f) Θ n / R e S S D / s c a l e .
Figure 1. Registration result of Algorithm 2 on image pairs Rectangle. (a) T ( · ) . (b) D ( · ) . (c) | T D | . (d) T h ^ K M ( · ) . (e) R e S S D change (Rectangle) with scale number. (f) Θ n / R e S S D / s c a l e .
Fractalfract 06 00681 g001
Figure 2. Registration result for different scales on image pairs Rectangle. (a) T ( · ) . (b) R e S S D = 73.29 % . (c) R e S S D = 43.84 % . (d) R e S S D = 14.75 % . (e) R e S S D = 1.65 % . (f) R e S S D = 1.3 % .
Figure 2. Registration result for different scales on image pairs Rectangle. (a) T ( · ) . (b) R e S S D = 73.29 % . (c) R e S S D = 43.84 % . (d) R e S S D = 14.75 % . (e) R e S S D = 1.65 % . (f) R e S S D = 1.3 % .
Fractalfract 06 00681 g002
Figure 3. Registration result of Algorithm 2 on image pairs C-E. (a) T ( · ) . (b) D ( · ) . (c) | T D | . (d) T h ^ K M ( · ) . (e) R e S S D change (C-E) with scale number.
Figure 3. Registration result of Algorithm 2 on image pairs C-E. (a) T ( · ) . (b) D ( · ) . (c) | T D | . (d) T h ^ K M ( · ) . (e) R e S S D change (C-E) with scale number.
Fractalfract 06 00681 g003
Figure 4. Registration result for different scales on image pairs C-E. (a) T ( · ) . (b) R e S S D = 37.35 % . (c) R e S S D = 1.83 % . (d) R e S S D = 0.97 % . (e) R e S S D = 0.46 % . (f) R e S S D = 0.27 % .
Figure 4. Registration result for different scales on image pairs C-E. (a) T ( · ) . (b) R e S S D = 37.35 % . (c) R e S S D = 1.83 % . (d) R e S S D = 0.97 % . (e) R e S S D = 0.46 % . (f) R e S S D = 0.27 % .
Fractalfract 06 00681 g004
Figure 5. Registration result of Algorithm 2 on image pairs Lung. (a) T ( · ) . (b) D ( · ) . (c) | T D | . (d) T h ^ K M ( · ) . (e) R e S S D change (Lung) with scale number.
Figure 5. Registration result of Algorithm 2 on image pairs Lung. (a) T ( · ) . (b) D ( · ) . (c) | T D | . (d) T h ^ K M ( · ) . (e) R e S S D change (Lung) with scale number.
Fractalfract 06 00681 g005
Figure 6. Registration result for different scales on image pairs Lung. (a) T ( · ) . (b) R e S S D = 6 % . (c) R e S S D = 3.08 % . (d) R e S S D = 2.98 % . (e) R e S S D = 2.91 % . (f) R e S S D = 2.88 % .
Figure 6. Registration result for different scales on image pairs Lung. (a) T ( · ) . (b) R e S S D = 6 % . (c) R e S S D = 3.08 % . (d) R e S S D = 2.98 % . (e) R e S S D = 2.91 % . (f) R e S S D = 2.88 % .
Fractalfract 06 00681 g006
Figure 7. Registration result of Algorithm 2 on image pairs Hand. (a) T ( · ) . (b) D ( · ) . (c) | T D | . (d) T h ^ K M ( · ) . (e) R e S S D change (Hand) with scale number.
Figure 7. Registration result of Algorithm 2 on image pairs Hand. (a) T ( · ) . (b) D ( · ) . (c) | T D | . (d) T h ^ K M ( · ) . (e) R e S S D change (Hand) with scale number.
Fractalfract 06 00681 g007
Figure 8. Registration result for different scales on image pairs Hand. (a) T ( · ) . (b) R e S S D = 12 % . (c) R e S S D = 6.18 % . (d) R e S S D = 4.96 % . (e) R e S S D = 4.4 % . (f) R e S S D = 4.04 % .
Figure 8. Registration result for different scales on image pairs Hand. (a) T ( · ) . (b) R e S S D = 12 % . (c) R e S S D = 6.18 % . (d) R e S S D = 4.96 % . (e) R e S S D = 4.4 % . (f) R e S S D = 4.04 % .
Fractalfract 06 00681 g008
Figure 9. Registration result of Algorithm 2 on image pairs Liver. (a) T ( · ) . (b) D ( · ) . (c) | T D | . (d) T h ^ K M ( · ) . (e) R e S S D change (Liver) with scale number.
Figure 9. Registration result of Algorithm 2 on image pairs Liver. (a) T ( · ) . (b) D ( · ) . (c) | T D | . (d) T h ^ K M ( · ) . (e) R e S S D change (Liver) with scale number.
Fractalfract 06 00681 g009
Figure 10. Registration result for different scales on image pairs Liver. (a) T ( · ) . (b) R e S S D = 25.59 % . (c) R e S S D = 10.53 % . (d) R e S S D = 6.63 % . (e) R e S S D = 5.99 % . (f) R e S S D = 5.88 % .
Figure 10. Registration result for different scales on image pairs Liver. (a) T ( · ) . (b) R e S S D = 25.59 % . (c) R e S S D = 10.53 % . (d) R e S S D = 6.63 % . (e) R e S S D = 5.99 % . (f) R e S S D = 5.88 % .
Fractalfract 06 00681 g010
Figure 11. Registration result of Algorithm 2 on image pairs Underwater square. (a) T ( · ) . (b) D ( · ) . (c) | T D | . (d) T h ^ K M ( · ) . (e) R e S S D change (Underwater square) with scale number.
Figure 11. Registration result of Algorithm 2 on image pairs Underwater square. (a) T ( · ) . (b) D ( · ) . (c) | T D | . (d) T h ^ K M ( · ) . (e) R e S S D change (Underwater square) with scale number.
Fractalfract 06 00681 g011
Figure 12. Registration result for different scales on image pairs of Underwater square. (a) T ( · ) . (b) R e S S D = 19.32 % . (c) R e S S D = 3.79 % . (d) R e S S D = 2.16 % . (e) R e S S D = 1.69 % . (f) R e S S D = 1.43 % .
Figure 12. Registration result for different scales on image pairs of Underwater square. (a) T ( · ) . (b) R e S S D = 19.32 % . (c) R e S S D = 3.79 % . (d) R e S S D = 2.16 % . (e) R e S S D = 1.69 % . (f) R e S S D = 1.43 % .
Fractalfract 06 00681 g012
Figure 13. Registration result of Algorithm 2 on image pairs of Underwater coin. (a) T ( · ) . (b) D ( · ) . (c) | T D | . (d) T h ^ K M ( · ) . (e) R e S S D change (Underwater coin) with scale number.
Figure 13. Registration result of Algorithm 2 on image pairs of Underwater coin. (a) T ( · ) . (b) D ( · ) . (c) | T D | . (d) T h ^ K M ( · ) . (e) R e S S D change (Underwater coin) with scale number.
Fractalfract 06 00681 g013
Figure 14. Registration result for different scales on image pairs of Underwater coin. (a) T ( · ) . (b) R e S S D = 12.48 % . (c) R e S S D = 2.95 % . (d) R e S S D = 1.85 % . (e) R e S S D = 1.55 % . (f) R e S S D = 1.42 % .
Figure 14. Registration result for different scales on image pairs of Underwater coin. (a) T ( · ) . (b) R e S S D = 12.48 % . (c) R e S S D = 2.95 % . (d) R e S S D = 1.85 % . (e) R e S S D = 1.55 % . (f) R e S S D = 1.42 % .
Fractalfract 06 00681 g014
Table 1. Quantitative comparison of registration results of four different algorithms (test for synthetic image pairs).
Table 1. Quantitative comparison of registration results of four different algorithms (test for synthetic image pairs).
DataAlgorithmRe-SSD (%)CPU/s
RectangleProposed1.3213
BD7116.8
LDDMM18.74143
TV31.56249.3
C-EProposed0.2783.7
BD0.7820.7
LDDMM3.15101.8
TV4.29210.2
Table 2. Quantitative comparison of registration results of four different algorithms (test for medical image pairs).
Table 2. Quantitative comparison of registration results of four different algorithms (test for medical image pairs).
DataAlgorithmRe-SSD (%)CPU/s
LungProposed2.88102.8
BD6.7315.9
LDDMM5.55123.2
TV4.6398.3
HandProposed4.04195.7
BD8.319.8
LDDMM5.99114.9
TV8.83112.2
LiverProposed5.88119.75
BD16.858.6
LDDMM15.14118.5
TV19.6630.3
Table 3. Quantitative comparison of registration results of four different algorithms (test for underwater distorted image pairs).
Table 3. Quantitative comparison of registration results of four different algorithms (test for underwater distorted image pairs).
DataAlgorithmRe-SSD (%)CPU/s
SquareProposed1.43242.4
BD8.288.5
LDDMM8.8620.6
TV2.99244.7
CoinProposed1.42272.65
BD6.9522.2
LDDMM58.99248
TV97.83213.14
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Du, Y.; Han, H. Multiscale Approach for Bounded Deformation Image Registration. Fractal Fract. 2022, 6, 681. https://doi.org/10.3390/fractalfract6110681

AMA Style

Du Y, Han H. Multiscale Approach for Bounded Deformation Image Registration. Fractal and Fractional. 2022; 6(11):681. https://doi.org/10.3390/fractalfract6110681

Chicago/Turabian Style

Du, Yunfeng, and Huan Han. 2022. "Multiscale Approach for Bounded Deformation Image Registration" Fractal and Fractional 6, no. 11: 681. https://doi.org/10.3390/fractalfract6110681

APA Style

Du, Y., & Han, H. (2022). Multiscale Approach for Bounded Deformation Image Registration. Fractal and Fractional, 6(11), 681. https://doi.org/10.3390/fractalfract6110681

Article Metrics

Back to TopTop