A robust non-local total-variation based image registration method under illumination changes in medical applications

Biomedical Signal Processing and Control 49 (2019) 96–112

Contents lists available at ScienceDirect

Biomedical Signal Processing and Control journal homepage: www.elsevier.com/locate/bspc

A robust non-local total-variation based image registration method under illumination changes in medical applications Khadijeh Aghajani a , Rohollah Yousefpour b , Mahbanou Zohrehvandi c a b c

Department of Computer Engineering, University of Mazandaran, Babolsar, Iran Department of Mathematical and Computer Sciences, University of Mazandaran, Babolsar, Iran Department of Electrical and Computer engineering, University of Malayer, Malayer, Iran

a r t i c l e

i n f o

Article history: Received 17 August 2017 Received in revised form 6 October 2018 Accepted 16 November 2018 Keywords: Nonrigid image registration Spatially-varying intensity distortion Weighted total variation

a b s t r a c t Independence of neighboring pixels and image stationarity are major concepts in conventional similarity metrics, used in image registration tasks. The accuracy of image registration decreases due to the presence of spatially varying intensity distortion in images. In this study, we hypothesized that changes in image illumination have limited total variation (TV). Accordingly, we introduced a similarity metric by reducing the weighted TV in the residual image. The primal dual method was then chosen to solve the proposed registration problem. The efficiency of the proposed method was compared to conventional methods, including the residual complexity (RC) method, the robust Huber similarity measure (RHSM), and the local linear reconstruction method (LLRM) which have been very successful in this field. The efficacy of the proposed method was confirmed by experimental findings on real-world and synthetic images. © 2018 Elsevier Ltd. All rights reserved.

1. Introduction Image registration is a process through which the coordinates of two images, which are possibly taken from the same object from different views, at different points of time, or by different sensors, are homogenized. It is known to have different medical, machine vision, and remote-sensing applications. Super resolution reconstruction [1], change assessment, robotic surgeries [2], and growth/development analysis [3,4] are among the image registration applications. So far, different methods of image registration, including intensity- and feature-based methods, have been introduced [4]. In the featured-based category, image registration includes feature extraction, finding correspondence between the features, and measurement of transformation function parameters [5]. In other words, this category elicits a subset of the corresponding points, versus considering the images’ overlapping pixels; the coordinates of those corresponding points are used for estimating the parameters of the transformation function. Some popular features in the stage of feature extraction include salient points (extracted by SIFT [6], SURF [7], and Harris [8]), geo-

E-mail addresses: [email protected] (K. Aghajani), [email protected] (R. Yousefpour), zohrehvandi [email protected] (M. Zohrehvandi). https://doi.org/10.1016/j.bspc.2018.11.001 1746-8094/© 2018 Elsevier Ltd. All rights reserved.

metric shapes (lines [5] and contours [9]), Gabor-type filters [10], and ˛-stable filters [11]. On the other hand, in the intensity-based category, parameters of transformation function are measured by optimizing the similarity function for pixels in both images. The process of image registration can be represented as follows [3,12]: T ∗ = argmaxT D(J, I(T )).

(1)

where J describes the reference image, I represents the moving image, and T is the geometric transformation function which aligns the images. The similarity metric, D(·, ·), is in fact a criterion to measure the similarity of the two images. By using optimization methods, the transformation function parameters can be obtained through maximizing the similarity metric. Consequently, introduction of a proper and robust similarity metric is a necessity in registration algorithms. Among the conventional similarity metrics, mutual information (MI), sum of squared differences (SSD), sum of absolute differences (SAD), and correlation coefficient (CC) are the most common [4]. These metrics are obtained by examining the difference in the intensities of the corresponding pixels, based on the assumption that intensity is spatially stationary and independent for every pixel [13]. The limitations of these methods include negligence of pixel dependence in real images and possible occurrence of spatially varying intensity distortion [14].

K. Aghajani et al. / Biomedical Signal Processing and Control 49 (2019) 96–112

97

Fig. 1. (a) Sample gray stripe; (b) the corrupted version using an additive spatially varying distortion field; (c) to (i) plots of similarity measures, including SAD, SSD, MI, CC, RC, RHSM, and TV relative to the translation of (b) over (a). Zero translation indicates precise image registration.

Conventional similarity metrics, including MI, CC, SSD, and SAD, are less robust in the event of complex intensity distortion and may fail in the registration of corrupted images. In order to illustrate this problem, an example is presented in Fig. (1). In this figure, image (b) is a corrupted version of image (a), which is contaminated using additive spatially varying intensity distortion. Images (c) to (i) show various similarity metrics, including total variation (TV) of residual images, for images created by horizontally relocation of image (b) over (a) with 14 pixel steps. Diagrams (c)–(f) do not present the optimal values with correct image alignment, while TV of the residual image (Diagram (i)) has a distinct optimum value with a wider range. Recently, more attention has been paid to image registration methods with illumination changes. The available methods can be divided into four categories. The first category includes methods with similarity metrics, such as local mutual information [15–18], in which spatially varying intensity distortion is supposed to be locally constant among neighboring pixels. Although the performance of local similarity measures is superior to global measures, local measures are both time- and space-consuming and show sensitivity to outliers and noise. The second category includes methods, which model pixel neighborhood interdependence on higher-order methods, such as MRF-based methods [19–21]. The third category is comprised of methods derived from structural image representation, which are applicable in multimodal image registration and intensity distortion of monomodal settings. Here, each pixel value is obtained depending on the local structure, not the intensity for encoding. Instead of using complex similarity metrics simpler met-

rics, such as SSD can be used to register the structural images. MIND [22], ALOST [23], entropy and Laplacian images [24], and LPCR [25] are some examples which belong to this category. The final category includes simultaneous intensity correction and registration. Here, solutions are obtained by considering the correction field and appropriate regularization terms for it. In the following two paragraphs, we consider some popular registration methods which belong to the fourth category. In the literature [13], the residual complexity (RC) metric has been used to solve the registration problem by assuming sparse spatially varying intensity distortion in the discrete cosine transform (DCT) domain. At first, the intensity correction field was removed from the objective function. This function is dependent on the residual image coding complexity in the DCT domain. In [26], an improvement to the RC method was performed by applying a multiplicative coefficient on the residual image. Ghaffari and Fatemizadeh defined a sparsity-based similarity metric, considering the residual image complexity in some transform domains, such as DCT or wavelet domains [27]. Moreover, Ghaffari and Fatemizadeh have been proposed another similarity metric called robust Huber similarity measure (RHSM) by using a maximum likelihood and M-estimator approaches [28]. They found that using both l1 and l2 norms of residual image in spatial transform domain can be an optimal similarity metric in the presence of Laplacian and Gaussian noises and spatially varying intensity distortion. So they utilized Huber norm of the residual which is a combination of l1 and l2 norms.

98

K. Aghajani et al. / Biomedical Signal Processing and Control 49 (2019) 96–112

In [29], it was assumed that illumination changes between two images in each region can be modeled by adding to and multiplying by the constant coefficients. Accordingly, an iterative method was proposed to detect the regions and estimate the coefficients in each region. Wang and Pan also proposed a dissimilarity measure, based on locally linear relationship between the two images [30]. The reference image in each local region was used to reconstruct the aligned floating image. In [31], a local linear reconstruction method (LLRM) is used to solve the registration problem. In this paper, it is assumed that the relationship between the reference and the moving image follows a linear relationship in the presence of spatially varying intensity distortion. Instead of an analytical solution, in [31], it has been assumed that multiplicative and additive coefficients, which are the same size of the image, have low TV. By applying this constraint and minimizing the sum of the squared of the reconstruction error, the registration parameters are estimated. In [32], after defining the multiplicative correction field, TV regularization term was used on the correction field. On the other hand, Keeling introduced an intensity-based similarity metric for image registration. To evade irregular scaling functions, Tikhonov regularization term was used on the scaling functions [33]. In [34], distortion is assumed to be an additive term; a ˇ-regularized ROF model was used as a regularization term. Here, it was hypothesized that changes in image illumination have low TV. The idea behind this assumption is that signals with perceivable details exhibit low total variation; that is, the integral of the gradient magnitude of a signal is low. This idea was first introduced by Rudin, Osher, and Rudin to reduce noise in images [38]. Overall, preservation of sharp discontinuities and noise elimination are the main objectives of TV-based methods. TV amount has been used as a regularization term on displacement fields in some image registration methods [35–37,39,40]. In this study, an intensity-based similarity metric was proposed to manage spatially varying intensity distortion. First, an intensity correction field was introduced, which brought the images into agreement. TV regularization term for the correction field was used, as the residual image is speculated to achieve minimum TV in correct alignment without outliers. In other words, we assumed that the residual image has few details. In addition, a weighting function on TV was used to diminish the smoothing effect on distortion across the edges. The weighting function was computed according to the similarity of the reference image’s neighboring pixels. It should be noted that, unlike the conventional TV-based image registration algorithm, we focused here on the use of TV regularization term on the intensity correction field. Because there was no need to solve the problem analytically, l1 norm was used as a data similarity term. It enjoys a higher resistance in the presence of the non Gaussian noise [35], and it is less sensitive to intensity variations [39]. Different approaches are available to resolve TV-based optimization problems. We used the primal-dual optimization for solving the proposed problem [41–43]. The difference between this paper and the method proposed in [31] is that, here, the distortion field is considered to be only additive. One of the disadvantages faced with the method presented in [31] was the spatial and time complexities of the proposed method. This is the result of considering the multiplicative and additive coefficients (along with the variables of the dual space). This makes it difficult to apply the algorithm on high-volume images (for example, remote sensing images and 3D medical images). Considering additive term (instead of multiplicative and additive terms) increased the speed compared to the previous version. However, there has not been any significant change in its performance, given the conducted experiments. Another difference between this research and [31] is considering the l1 norm of the changes in each direction (in TV regularization term) instead of Huber norm.

The rest of the paper has been organized as follows: in Section 2, the problem is generally defined and a numerical algorithm is proposed for resolving the problem. Section 3 focuses on the observations and experiments in the present study. Moreover, the proposed algorithm is compared with conventional algorithms in terms of performance, using real-world and simulated data. Section 4 presents the final conclusions. 2. Main idea and the proposed method 2.1. Definition of the problem Considering that the two images J :  ∈ R2 → R, I :  ∈ R2 → R are to be aligned, the following intensity relation between aligned images was assumed. I(T (x)) = J(x) + S(x) + (x)

(2)

in which I and J are the moving and the reference images, respectively.  = {x = (x1 , x2 )|1 ≤ x1 ≤ N, 1 ≤ x2 ≤ M} where N and M are the sizes of each dimension of the images. S is the additive correction field and T(x) is the transformation function that aligns I onto J. (x) is the reconstruction error which is modeled by a zero-mean Laplacian distribution. According to the MAP (maximum a posteriori) approach and Eq. (2), the registration problem can be described as follows: Problem: Given the reference image J and the moving image I, the transformation function, T, and the additive correction field, S, should be described to minimize the following energy function: E(S, T ) = D(I(T ), J + S) + ˛R(S) + Q (T )

(3)

where D(·, ·) denotes the data similarity term, R(S) and Q(T) are the regularization terms applied on S and T, respectively, D(·) describes dissimilarity between images I(T) and J + S. D(·, ·) is defined as follows: D(I, J + S) =



I(T (x)) − J(x) − S(x)1

(4)

x∈

in which ||·||p is the lp norm in Rn . As has been mentioned before, it was assumed that the additive correction field, S, has low TV amount. Thus R(S) was defined by the following equation: R(S) =

i=8 

wi (x)|(AS(x))i |

(5)

x ∈  i=1

Let us define 1 = (1, 0), 2 = (1, − 1), 3 = (0, − 1), 4 = (−1, − 1), 5 = (−1, 0), 6 = (−1, 1), 7 = (0, 1), 8 = (1, 1) and  = {i |i = 1, . . ., 8}. By these definitions we introduced operator A :  → R8 whose ith component (i = 1, 2, . . ., 8) at any given location x can be computed as:



i

(AS(x)) =

S(x + i ) − S(x)

if x + i ∈ , i = 1, 2, . . ., 8

0

otherwise

(6)

In simple term, this operator computes the intensity change at each direction i. wi (x) is a weighting function in direction i that is computed as Eq. (7). wi (x) = e−distance(x,x+i ,r)/

(7)

Here, distance(x, y, r) is measured by summing the squared differences in the corresponding pixels of two patches (size, (2r + 1)2 ), centered at x andy, respectively. distance(x, y, r) =



p ∈ Pr

G(p)(J(x + p) − J(y + p))2

(8)

K. Aghajani et al. / Biomedical Signal Processing and Control 49 (2019) 96–112

In which J is the reference image and Pr = {p = (p1 , p2 )||pi | ≤ r, i = 1,

99

2.3. Dual step

−p2 /((1/3)2 )

2 is a weighting function, used to reduce the 2}. G(p) = e effect of pixels, which are far from the center. Here, r was set to 2. The reference image was used to obtain the weighting function, so its computation time can be ignored, because the reference image is fixed during registration. Putting together the aforesaid notations and relations, the objective function can be rewritten as Eq. (9):

E(T, S) =



Apply gradient ascent method to the maximization problem (14). So the solution can be obtained by: [41,42] 

n+1

= (I + n ∂F ∗ )

The mapping (I

−1

n

( + n ˛−1 AS n ).

−1 +  ∂F ∗ )

wi (x)|(AS(x))i | + Q (T ),

(9)



x ∈  i=1

n+1

( i )

2.2. Numerical method

n+ 1

2

( −  2 ) = argmin { + F ∗ ()} 2n

i=8



n+1

[x]ba =

= ( i )

n+ 1 2

w i

minS,T max (, S, T ) = ˛−1 A∗ , S − F ∗ () + G(S) + Q (T ), 8NM

(10) F* (·)

⎧ ⎨ b if (x > b) ⎩

x

if (a ≤ x ≤ b)

a

if (x < a)

S

= (I + n ∂G)

−1

(S n − n ˛A∗ 

⎧ 1 ⎪ ⎨

we can rewrite Eq. (20) as

if  ∈ K

(12)

∞ o.w.

K = { = ( 1 , . . .,  8 ) : | i | ≤ wi ∀ i = 1 : 8}.

(13)

In two consecutive dual and primal steps, the primal-dual method can solve the problem. In the first step of nth iteration for a fixed Sn value, the following dual problem is solved: = argmax (, S n , T n )

(14)

In the second step, assuming that  n+1 is an approximate solution of the equation (14), Sn+1 and Tn+1 are computed by solving the following Primal problem. (S n+1 , T n+1 ) = argminS,T (

(20)

(21)

otherwise

(22)

where ∇ S  is the gradient descent direction of  with respect to S. ∇ S  can be found by

∇S =

∂ {˛S, A∗  + I(T ) − J − S1 } ∂S

(23)

in which the adjoint operator A* is defined as:

Here, K is given by

n+1

0

S n+1 = S n − n ∇ S ,

where IK denotes the so-called indicator function of the convex set K given by



).

if (10−8 < x)

F ∗ () = sup, AS − F(AS) = IK ()

(11)

n+1

If we approximate the derivative of |·| as

∂ |x| = −1 if (x < −10−8 ) ∂x ⎪ ⎩

S

(19)

Similar to the dual step the solution of (15) is obtained as n+1

where  · is the inner product, ∈ R is the dual of AS, and is the Legendre Fenchel conjugate of the function F(·) which can be computed as follows.

0

(18)

−wi

2.4. Primal step



wi (x)|(AS(x))i | and G(S) = x ∈  ||I(T( x)) − J( F(AS) = x∈ i=1 x) − S( x)||1 . Using the Legendre–Fenchel conjugate, the image registration problem can be represented as:



(17)

in which [ · ]ba is interpreted as:

In order to solve the minimization of Eq. (9), we used an iterative method called primal dual approach. This approach was proposed for solving problems in the form of minX F(AX) + G(X) [41,42]. In the original formulation, F : Y → R+ and G : X → R were convex functions, A was an operator from the space X to space Y, and Y was the dual space to X. In our problem, we defined the following convex functions:

IK () =

n

The solution can be determined as follows:

where ˛ and  are the trade-off parameters.



1

=  + n AS n , The solution of is denoted Prox (16) is determined through minimizing the following problem:

x∈

+˛

is called the “proximal map” of  ∂F* and

n+ F* . Let us define  2

I(T (x)) − J(x) − S(x)1

i=8 

(16)

n+1

, S, T )

(15)

We can stop the operation if it reaches a certain number of iterations or changes of the objective function becomes lower than a threshold. By considering  n and n as two step sizes (for dual and primal steps, respectively) in nth iteration, the following optimization algorithm updates the solutions as follows.

A∗ (x) =

i=8 

 i (x − i ) −  i (x).

(24)

i=1

The remaining point of this section discusses how we can update the transformation function. The updating was done in this step, according to the type of the transformation function. In this study, the cubic B-spline model was used to model the transformation function. Its parameters, the displacements of the control points (),were updated by computing ∇ T  and ∇  T. ∇ T  can be computed easily. Moreover, for estimating T a regularization term, Q(T), is considered. In order to prevent unnatural wraps, the Euclidean distance between all displacements of neighboring B-spline control points was penalized ( denotes the weight of the regularization term). According to the aforesaid explanations and equations presented in this section, the proposed iterative algorithm is summarized in Algorithm 1. Algorithm 1. Proposed iterative algorithm for image registration.

100

K. Aghajani et al. / Biomedical Signal Processing and Control 49 (2019) 96–112

2.5. Multiresolution implementation In order to improve the accuracy of registration methods and decrease computational complexity, a multi-resolution scheme can be used [12,13,30]. In this situation, by registering the lowerresolution images, initial values of the transformation function’s parameters for the higher resolution are obtained. The registration can be continued at higher resolution by considering the obtained initial transformation function. In the present study, we have also used this idea. The steps of multiresolution scheme are summarized in Algorithm 2. Algorithm 2. Multiresolution scheme for image registration algorithm.

3. Experimental results To model the transformation function T, free-form deformation (FFD) was used with hierarchical B-spline control points [44]. FFD is applicable for modeling non-rigid continuum motions. In this study, we considered T(x) = x + ˚(x), where the displacement fields ˚(x) can be obtained according to , which are the displacements of the control points on a B-spline control grid of size nx1 × nx2 . In the two-dimensional case, ˚(x) can be obtained using the following equation: ˚(x) = ˚(x1 , x2 ) =

3 3  

Bl (u)Bm (v)i+l,j+m

(25)

l=0 m=0 x1

where i = n  − 1, j = x1

2

nx x2

 − 1, u =

x1 nx1

1 − nx x1

, v =

x2 nx2

2 − nx x2

.

The Bl s represent the lth basis function of the B-spline as B0 (u) =

(1−u)3 , B1 (u) 6

3

2

3

2

points is an important parameter in this function. The transformation function becomes more global if the control points are far from each other. On the other hand, a shorter distance causes local deformation and leads to complex measurements. Therefore, determining the distance between the neighboring control points is a tradeoff between the model and computational complexities. Considering the extent of deformation in the experiments, a distance of 12–25 pixels seems appropriate. We evaluated the proposed method and compared it with other approaches in terms of image registration performance. In the quantitative analysis of the algorithms, one of the common problems is the absence of a ground truth in real-world images. Therefore, experiments have been carried out, based on the simu-

lated information in various studies. In the present study, non-rigid experiments were performed on the simulated and real data. Overall, the real-world images included clinical 4D CT scans, retinal images, neurosurgery images, and digital subtraction angiography (DSA) image pairs. The findings of the experiments are reported in two sections for the real-world and simulated images. Finally, the results were compared with SSD, MI, CC, RC, RHSM and LLRM methods, using the Medical Image Registration Toolbox (MIRT).1 In each section, setting of the parameters can be done according to Table (1) . 3.1. Simulated data Two-dimensional slices (216 × 180) of MR (Magnetic Resonance) Brain-web2 images were used in the experiments. Geometric and intensity distortion were created in the images as

3

+4 = 3u −6u , B2 (u) = −3u +3u6 +3u+1 , B3 (u) = u6 . 6 Obviously, changes in displacement of one control point influences displacement in the local neighborhood of the manipulated point. Mesh size or distance between the neighboring control

1 Medical Image Registration Toolbox for Matlab, https://sites.google.com/site/ myronenko/research/mirt. 2 http://www.bic.mni.mcgill.ca/brainweb/.

K. Aghajani et al. / Biomedical Signal Processing and Control 49 (2019) 96–112

101

Table 1 Setting of the required parameters.

Section 3.1 Section 3.2.1 Sections 3.2.2– 3.2.4

Level of multiresolution

grid size

MaxIter





˛ (PM)

3

12 × 12 15 × 15 × 15 20 × 20

300 150 300

10

0.1, 0.3 0.05 0.1–1

0.05

4

0.1

follows. First, one of the slices was considered to be the floating image I. Then, by applying a geometric transformation using ˆ true parameter to image I, the reference image J was a specific  created. To simulate the geometric transformation function, the cubic B-spline model was used. By disturbing a uniform grid of ˆ true was determined. Random perturbapoints (12 × 12 pixels),  tion was drawn in (−5.9:5.9) intervals from a uniform distribution. Finally, by registering the image I to J, the transformation function ˆ est , was estimated. parameter,  For quantitative evaluation of the presented similarity metˆ true and ric, two error parameters were used [13]. According to  ˆ est , the root mean squared error (RMSE) of transformation can be  obtained using the following equation:



transformationRMSE =

1  ˆ true ) − T (x,  ˆ est )22 . T (x,  ||

(26)

x∈

The intensity RMSE can be also calculated from the difference between the two images: the distortion free reference image and the registered image. The registered image can be obtained if the transformation function is applied to the original clean image (without intensity distortion). To evaluate the performance of the algorithm, spatially varying intensity distortion was created in addition to distortion of the geometric transformation function, and is considered to be both additive and multiplicative. In the experiments of this section, the ˛ parameter in Eq. (9) was set to 0.05 and the level of the image pyramid in the multiresolution schema was set to 3. 3.1.1. Selection of parameter  One of the parameters identified in image registration is the weight devoted to the regularization term in the transfer function, , as a tradeoff between transformation smoothness and image similarity. In this function, the regularization term can control local deformation. A higher value is considered if the transformation function mostly has a global nature, as well as low local deformation. Before carrying out experiments related to illumination changes, we briefly deal with establishing . Given that methods such as SSD, MI, and CC present ideal results without occuring illumination change, we first conducted an experiment on images by using the procedure explained earlier but this time without consideration illumination changes. The errors obtained using different methods at different values of  are plotted in Fig. (2). The  value was chosen from set = {0.01, 0.05, 0.1, 0.2, 0.5, 1, 2, 5, 10, 50}. The diagram shows that high values of  results in transfer function estimations that are characterized by very low levels of local deformation. The transfer function may also tend toward the identity function under increasing  value. However, very low values of  may generate unreal function estimations with high levels of unwanted warps because of the decreased effect of the regularization term. This can produce numerous errors in the presence of local minima. It has to be noted that, in the absence of distortion field, (the probability of this occurrence is very low in actual images), the results of the conventional algorithms are accurate. However, to be honest, the proposed method has presented, with a very slight difference, a little weaker performance than these methods because of the consideration of the term S, which is practically non-existent. Of course,

Fig. 2. Transformation error of various methods with respect to different parameter  settings.

it is expected that S will be completely zero, nevertheless, though small, its precision is slightly diminished in some points. Meanwhile, in this chart, the goal is examining the effect of changing the parameter  on registration error in each method. Having analyzed this experiment and all other experiments on illumination changes, we selected 0.1 as the  value for RC, RHSM, LLRM and the proposed method. However, because we observed numerous errors using SSD, MI, and CC methods, we used a higher value of  (i.e., 0.3) for these methods. Despite the high efficiency of the PM, it was more timeconsuming compared with some conventional methods such as SSD or CC. In Table 2, the computation times for the various methods applying on various image sizes are compared. The numbers reported in this table indicate that the time consumed by the proposed method for registering high-volume images has decreased compared to LLRM. 3.1.2. Additive intensity distortion The intensities of the floating and reference images were corrupted using Eq. (27). This equation is similar to the one used in [13,30,45]. −x − k  k=K 2 1 I(x) = I(x) + g(x, C) + exp 2( k ) K

(27)

k=1

rescale I to [0, 1] Here, the first term, I, represents the original image. The second term, g, represents a smoothly varying global intensity, while the last term shows a locally varying intensity field, including a combination of K Gaussian densities. Additionally, g(x, C) is the polynomial function of degree 3, with the random coefficient set C, where each element of C, ci , is randomly selected from the interval (−1 :1). K is a random variable belonging to the (1:8) interval, and k was selected randomly from the interval (25:35). k was also randomly selected from . This equation creates more complex intensity distortion than [13,30]. For example, in [13], the value of k was assumed to be 30, and the value of K was the same for both

102

K. Aghajani et al. / Biomedical Signal Processing and Control 49 (2019) 96–112

Table 2 Sizes of some image pairs used for evaluation and the average computation times (in minutes) of the various methods applying on them. Image size

180 × 250

280 × 330

500 × 500

650 × 800

256 × 256 × 94

com. time using SSD, CC com. time using MI com. time using RC, RHSM com. time using LLRM com. time using PM

0.18 0.19 0.26 0.42 0.23

0.25 0.30 0.51 0.82 0.48

1.1 1.8 1.1 2.8 1.6

2.1 3.2 2.2 5.3 2.4

11.1 22 16 56 23

Fig. 3. Synthetic experiment. (a) Original image. (b) Reference bias. (c) Moving bias. (d) The difference between moving and reference biases. (e) Moving image. (f) Original transformation. (g) Reference image. The results of registering (e) onto (g) using various methods are shown in Figs. 4 and (5).

images. K and k (k = 1 : K) can have different values for each image in each run. Fig. (3) is an example of how the experiments were conducted. The results of registering image (Fig. 3e) onto image (Fig. 3g) are presented in Figs. (4) and (5). According to the explanations presented regarding the establishment of  and its effect on output, we display the images created with different methods per two different values of  (0.1 and 1) in Figs. 4 and 5. The estimated transfer function and the produced image for each method are also presented. Fig. 4 shows that methods such as SSD, MI and CC produce erroneous results under illumination changes, thereby causing failed registration. Although the number of errors can be reduced by increasing the  value, this error reduction does not necessarily translate to registration of higher accuracy. This is confirmed by the fact that stringent restrictions are imposed on the transfer function. To more precisely examine and compare the efficiency of different methods, we displayed composite views through contour overlap after registration using different approaches in Fig. (6). To reach the error metric of the methods noted earlier, in each modality (PD, T1, and T2), fifteen two-dimensional slices were chosen randomly from the central slices of the three-dimensional data. For each slice, the explained evaluation process was repeated 20 times. During each run, the random intensity distortion and spatial distortions were reinitialized. The means and the standard deviations of the obtained RMSEs are reported in Table (2). The values reported in Table 2 reflect that SSD, MI, and CC methods are not applicable in the event of spatially varying intensity distortion. In the experiments, the RC and RHSM methods were sometimes observed to diverge. Errors made in these cases were ignored in the report.

3.1.3. Multiplicative distortion field On some occasions, changes in image intensity are multiplicative [13]. Although this hypothesis was not used to solve the problem, the accuracy of the proposed method was examined in the presence of multiplicative distortion fields. For modeling the intensity corruption of each image, the following equation was used:

I(x) = I(x).(0.2 + g(x, A) +

k=K 1

K

−x − k 22 exp

2( k )2

)

(28)

k=1

rescale I to [0, 1] Eq. (28) is also similar to those used in articles [13,30]. In the previous subsection, the terms of this equation were investigated. The registration methods were applied to the two-dimensional slices using various similarity metrics. Fifteen two-dimensional slices from each modality (T1, T2, and PD) were chosen randomly from the primary three-dimensional images. In each experiment, geometrical transformation and intensity distortion were reinitialized randomly. The mean (SD) RMSEs for each method are reported in Table 3. The values presented in this table indicate the inefficacy of SSD, MI, and CC methods in registering images in the event of spatially varying intensity distortion. In these experiments, when calculating the RMSEs, the high amounts of RC errors were not considered. The assumption of the multiplicative distortion field was not considered, and the results revealed that RC, RHSM and the proposed method are robust with this kind of distortion, while other methods, like the one using an additive distortion field, perform poorly. Moreover, LLRM, due to the consideration of the multiplicative and

K. Aghajani et al. / Biomedical Signal Processing and Control 49 (2019) 96–112

103

Fig. 4. Registration of the moving image (b) on the reference image (a). (c) True transformation. Results of registration (estimated transformation and the registered image) using SSD, MI, CC and RC methods per two different values of  (0.1 and 1). (d)–(g) SSD results. (h)–(k) MI results. (l)–(o) CC results. (p)–(s) RC results.

Table 3 Registration performance of SSD, MI, CC, RC, RHSM, LLRM, and proposed method (PM). Image corruption is due to additive distortion field (top, mean (SD) of transformation RMSEs; bottom, mean (SD) of intensity RMSEs). RMSE

SSD

MI

CC

RC

RHSM

LLRM

PM

Trans. Inten.

20.97(11.73) 0.27(0.041)

13.96(8.7) 0.29(0.03)

9.32(7.33) 0.26(0.03)

0.98(1.21) 0.091(0.017)

0.86(0.40) 0.090(0.0026)

0.86(0.70) 0.0822(0.008)

0.84(0.62) 0.0813(0.007)

104

K. Aghajani et al. / Biomedical Signal Processing and Control 49 (2019) 96–112

Fig. 5. Registration of the moving image (b) on the reference image (a). (c) True transformation. Results of registration (estimated transformation and the registered image) using RHSM, LLRM and PM methods per two different values of  (0.1 and 1). (d)–(g) RHSM results. (h)–(k) LLRM results. (l)–(o) PM results. (p) and (q) estimated S per two values of .

additive distortion field, is slightly more accurate than the proposed method. 3.1.4. Gaussian noise The efficiency of the proposed method in the presence of Gaussian noise was investigated. It is assumed that one of the images was corrupted by Gaussian noise. Fifteen slices from the T2 image

were used in these experiments. In Fig. (7), the mean and SD of transformation and intensity RMSEs are depicted with respect to changes in noise variance. According to the diagram of transformation RMSEs, RC, RHSM, LLRM and the proposed method were similarly efficient. In spite of considering l1 norm norm as a dissimilarity term in Eq. (9) (likewise for l2 norm in RC), the performances were similar. Although we sometime see a reduction in precision

K. Aghajani et al. / Biomedical Signal Processing and Control 49 (2019) 96–112

105

Fig. 6. Composite views through contour overlap after registration using various methods. (a) SSD result. (b) MI result. (c) CC result. (d) RC result. (e) RHSM result. (f) LLRM result. (g) PM result.

Fig. 7. The registration performance of the proposed and conventional methods. The reference image is contaminated by additive zero mean Gaussian noise ((a) and (b) represent the transformation and intensity RMSEs with respect to different levels of noise variance, respectively).

(even partial) of the proposed method compared to some methods, but generally the proposed method showed a slightly more effective performance in comparison with RC and RHSM methods.

3.1.5. Sensitivity to the parameter ˛ The sensitivity of the proposed method with respect to the amount of the ˛ parameter was examined. Fifty runs of both random geometrical distortion and additive intensity distortion (using Eq. (1)) were initialized and carried out for each ˛ ∈ {0.02, 0.035, 0.05, 0.07, 0.1, 0.15, 0.2, 0.5}. In Fig. (8), the means and standard deviations of transformation and intensity RMSEs obtained for various amounts of ˛ are plotted. As can be seen in this figure, ˛ values within the range of [0.035–0.15] produce similar acceptable results with limited advantage ˛ = 0.05.

3.1.6. Sensitivity to the weighting function The sensitivity of the proposed method to the weighting function was examined to investigate the effect of the weighting function on the performance of the proposed method by changing  in Eq. (7). A test similar to the one carried out in Section 3.1.1 was conducted using three slices from the T1 image and the additive distortion field (Eq. (27)). Fig. 9 shows the average error obtained with different values of . When  =∞, wi becomes 1, which means weight has no effect on Eq. (5). As the figure shows, the use of the weighting function slightly improves registration accuracy.

3.2. Real-world data Real-world images were used for the assessment, including neurosurgery images, retinal images, 4D CT scans, and DSA image pairs.

106

K. Aghajani et al. / Biomedical Signal Processing and Control 49 (2019) 96–112

Fig. 8. (a) Transformation and (b) intensity RMSEs in different parameter ˛ settings.

Fig. 9. (a) Transformation and (b) intensity RMSEs in different parameter  settings.

3.2.1. Registration of inhalation and exhalation images The 4D CT images were used to evaluate the proposed method, as well as SSD, MI, CC, RC, RHSM and LLRM methods. Ten CT scan pairs of maximum inhalation and exhalation phases were used (supplied by DIR-LAB at the University of Texas3 ). To highlight the rib cage, the images were cropped. Each case was identified according to the label. The image dimensions are presented in voxel, and voxel dimensions are indicated in millimeters for each case, as presented in Table (4). These images are three-dimensional. All the calculations for this mode are similar to those for the two-dimensional images. Here, a neighborhood was considered to be in three dimensions, but because of the high computational complexity, only six neighborhoods were considered for each voxel (two neighborhoods per dimension). In other words,  vectors were considered as follows:

3 Thoracic 4D CT images acquired as part of the standard planning process for the treatment of thoracic malignancies at The University of Texas, http://www.dir-lab. com/Downloads.html.

1 = (1, 0, 0), 2 = (−1, 0, 0), 3 = (0, 1, 0), 4 = (0, − 1, 0), 5 = (0, 0, 1), 6 = (0, 0, − 1), and  = {i |i = 1, . . ., 6}. Accordingly, the manner by which operator A works in this case can be easily generalized for the consideration of the third dimension. Some changes therefore occur in the dimensions of the dual variables. Calculating the weight function is simply a generalizable process. Such changes were considered in dealing with Eq. (7), where r is the radius of a circle. This circle was converted into a sphere with a radius of 2. The three-dimensional B-spline was also used as the transfer function. In this section, the level of multiresolution scheme is considered 3,  was set to 0.05, and the grid size was 15 × 15 × 15. Because of the large volume of images and the time-consuming nature of running different methods, the maximum number of iterations of Algorithm 1’s main loop was set to 150. In the present study, thoracic imaging experts identified the anatomical landmarks. The target registration error was determined by measuring the Euclidean distance between manually located landmarks in the reference image and their corresponding positions after algorithm application in the moving image. Fig.

K. Aghajani et al. / Biomedical Signal Processing and Control 49 (2019) 96–112

107

Fig. 10. Registration of a 4D CT scan image pair. This figure shows three image pairs of patient No. 8 in the transverse, coronal, and sagittal views (first, second, and third rows, respectively). The first column (a, d and g) presents the composite views through contour overlaps before image registration. The second column (b, e and h) indicates the composite views through contour overlaps following image registration with the proposed method. The third column (c, f and i) represents the estimated intensity distortion for the same slice of each orientation.

Fig. 11. The accuracy of 3D volume registration (mean and SD of errors using various methods).

108

K. Aghajani et al. / Biomedical Signal Processing and Control 49 (2019) 96–112

Table 4 Registration performance of SSD, MI, CC, RC, RHSM, LLRM, and proposed method (PM). Image corruption is due to multiplicative distortion field (top, mean (SD) of transformation RMSEs; bottom, mean (SD) of intensity RMSEs). RMSE

SSD

MI

CC

RC

RHSM

LLRM

PM

Trans. Inten.

13.49 (13.00) 0.06 (0.03)

8.75 (10.75) 0.042 (0.038)

8.35 (10.62) 0.037 (0.027)

1.01 (1.26) 0.012 (0.018)

1.09 (0.75) 0.012 (0.006)

0.80 (0.75) 0.008 (0.007)

0.84 (0.85) 0.009 (0.008)

Fig. 12. Registration of a retinal image pair: (a) reference image; (b) moving image; (c) composite view before image registration through contour overlaps; (d) to (j) composite view after registration using the SSD, MI, CC, RC, RHSM, LLRM, and the proposed method, respectively.; (k) estimated correction field; and (l) estimated transformation with the proposed method.

Table 5 The image volume dimensions are presented in voxel (voxel dimensions are represented in mm). The mean and SD of displacements in the corresponding landmarks are presented before registration. Case

Image dimension

Voxel size (mm)

Displacement before registration (mm)

1 2 3 4 5 6 7 8 9 10

256× 256× 94 256×256×112 256×256×104 256×256× 99 256×256×106 512× 512 × 128 512× 512 × 136 512× 512 × 128 512× 512 × 128 512× 512 × 120

0.97×0.97×2.5 1.16×1.16×2.5 1.15×1.15×2.5 1.13×1.13×2.5 1.10×1.10×2.5 0.97× 0.97 × 2.5 0.97× 0.97 × 2.5 0.97× 0.97 × 2.5 0.97× 0.97 × 2.5 0.97× 0.97 × 2.5

3.84 (2.78) 4.16 (3.95) 6.72 (4.10) 9.66 (4.93) 7.25 (5.55) 10.74 (6.94) 10.90 (7.48) 14.76 (9.05) 7.71 (4.01) 8.23 (6.86)

(10) presents the results of image registration using the proposed method for case No. 8. The average amount of errors and the corresponding SD for each case after registration with various methods are shown in Fig. (11). The mean and SD of errors before registration are reported in Table (5). By considering all 2775 landmarks (300 landmarks for cases 1–9 and 75 landmarks for case 10), the mean and SD of error in SSD, MI, CC, RC, RHSM, LLRM, and proposed method were 2.96 (2.748), 3.316 (2.714), 3.006 (2.795), 2.833 (2.671), 2.363 (2.097), 1.957 (1.84), and 1.942 (1.82), respectively.

Evidently, the proposed method exhibited a more effective performance, compared to others. 3.2.2. Retina images For evaluating the proposed method, two pairs of retinal images were used4 ; each pair was obtained within a 2-year interval. Overall, these images can be employed to diagnose diseases, such as glaucoma, diabetic retinopathy, and age-related macular degeneration [46]. Fig. (12) shows examples of these images. In such cases, image registration using intensity-based methods is quite challenging due to different intensity artifacts, such as blood vessels and non-uniform background [46,47]. As mentioned earlier, the proposed method considers the intensity distortion by S term; therefore, we expected it to show good performance on these images, as well. The problem was to register Fig. 12b onto Fig. 12a. Because the transformation function’s nature seemed to be global, a greater penalty coefficient was used for the regularization term on T. Moreover, the level of multiresolution scheme was considered 4 and the grid size was 20 × 20. In the experiments of this section and next sections, we have checked the registration results using different values of  and we have reported the best output obtained by each method. The corresponding  value of each result (using different methods) has

4

http://cmm.ensmp.fr/Anciens/zana/registration.html.

K. Aghajani et al. / Biomedical Signal Processing and Control 49 (2019) 96–112

109

Fig. 13. Registration of a retinal image pair: (a) reference image; (b) moving image; (c) composite view before image registration through contour overlaps; (d) composite view after registration using the RC method; (e) composite view after registration using the RHSM method; (f) composite view after registration using the LLRM method; (g) composite view after registration with the proposed method; (h) estimated correction field; and (i) estimated transformation with the proposed method.

Table 6 Registration of retina image pairs (RMSDs before and after registration obtained with different methods).

Image pair in Fig. 12 Image pair in Fig. 13

Before registration

SSD

MI

CC

RC

RHSM

LLRM

PM

16.5 9.2

15.6 11.8

15.5 8.1

16.2 9.1

5.5 4.1

6.5 1.3

5.8 1.3

3.2 1.3

been displayed above the obtained image. Fig. (12) presents the overlapping composites and contours of images obtained with the conventional and proposed methods. Evidently, the registration process failed when using SSD, MI, or CC method. It seems that RC and RHSM have nearly similar performances to LLRM and the proposed method, although the proposed method has a better performance. Another example of a retina image pair is shown in Fig. (13). The results of registering the images using RC, RHSM, LLRM and the proposed method are illustrated in this figure. Moreover, for better comparison, approximately 20 corresponding points are extracted manually for each pair. The registration error was obtained by computing the root mean squared distances (RMSD) between the locations of the corresponding points after registration. The obtained RMSDs for each image pair are reported in Table (6).

3.2.3. Digital subtraction angiography DSA is one of the most famous techniques in vascular imaging. In this approach, in order to highlight the blood vessels, a contrast agent is used. To investigate these vessels and to eliminate the background, assuming that the vessels have not moved, the image after injection is subtracted from the image taken before injection. In practice, this assumption is violated for some reasons such as breathing or movement by the patients, which leads to the appearance of unwanted artifacts in the subtraction of the two images. Examples of these images in DSA application are shown in Fig. 14a and b. In Fig. 14c, the obtained subtraction without any motion compensation is shown. Due to the nature of these images, using a registration algorithm as a preprocessing phase is required. As can be seen in Fig. 14a and b, this kind of registration problem has the challenge of a spatiallyvarying intensity distortion field. This leads to a decrease in the

110

K. Aghajani et al. / Biomedical Signal Processing and Control 49 (2019) 96–112

Fig. 14. Application of the image registration algorithm for coronary angiography images: (a) before the contrast agent injection; (b) after injection; (c) DSA image without motion compensation; (d)–(j) DSA images after compensation with SSD, MI, CC, RC, RHSM, LLRM, and proposed method, respectively; (k) estimated correction field; and (l) estimated transformation obtained with the proposed method.

efficiency of the conventional algorithms. However, it is expected that the RC method and the proposed method perform more desirably due to accounting for the intensity correction field. From Fig. 14d–j, the DSAs obtained from applying SSD, MI, CC, RC, RHSM, LLRM, and the proposed registration algorithms are displayed. As shown in these figures, there are fewer artifacts in DSA obtained by LLRM and the proposed registration algorithm.

3.2.4. Neurosurgery image pairs with outliers An important procedure in planning different neurosurgeries is imaging for purposes such as brain lesion resection. The proposed method was validated on a challenging image pair (i.e., MRI images of brain tumor resection). However, the registration methods could not manage the outliers robustly. It can be concluded that the obtained function is an oversmooth approximation of the real transformation function. Nevertheless, in continuum motions, it is not a very challenging problem, as it can be handled via changing the regularization term to, for instance, l1 norm [48]. The model can be considered efficient if the outlier space (volume and/or surface) is small in comparison to the whole image. Fig. (15) presents a preoperative T2 fluid-attenuated inversion recovery (FLAIR) MRI image, as well as an intraoperative image of the extensive resection of a T2 hyperintense lesion. In addition, this figure presents the application of conventional and proposed methods for registration of moving image on the reference image.

Here, the setting of the parameters were done similar to the values in Section 3.2.2. It is clear that, despite of presence of the outliers in these image pairs the proposed method and some conventional methods have efficient performances in registering the images. 4. Conclusion To register the mono-modal images, a robust similarity measure was presented for illumination changes assumed to have a low TV (total variation) between two images. Accordingly, an algorithm was presented to simultaneously register the images and correct intensities. The algorithm was effective in reducing the measure of weighted TV in the intensity correction field of the images. The weighting function of TV was used to decrease the smoothness effect across the edges of residual images. The weighting function was computed according to the similarity of the neighboring pixels of the reference image. The objective function presented here has two components: the l1 norm reconstruction error and the weighted TV regularizer. The primal dual method was used for numerical solving of the proposed optimization problem. To reduce computational complexity, a multiresolution schema was utilized. To evaluate the algorithm performance, both real-world and simulated data were used. Based on the simulated data, a more efficient model for illumination changes was proposed. It was observed that conventional methods, including SSD, MI, and CC, could not

K. Aghajani et al. / Biomedical Signal Processing and Control 49 (2019) 96–112

111

Fig. 15. (a) A pre-operative axial T2 FLAIR MR image containing a T2 hyper-intense mass lesion infiltrating the left temporal lobe and insula; (b) an intra-operative image showing extensive resection of the T2 hyper-intense lesion. (c) Composite view of contour overlap before registration. (d)–(j) Composite views of contour overlap after registration using SSD, MI, CC, RC, RHSM, LLRM, and the proposed method, respectively. (k) Estimated correction field and (l) estimated transformation using the proposed method.

register images, containing complex illumination changes. This behavior has rarely been observed for the RC and the RHSM methods. However, the behavior of the proposed algorithm was quite robust. Despite the elimination of the highly probable errors in the other algorithms, the means and standard deviations of reported errors in the simulated data were high. Experiments using the real data demonstrated the effectiveness of the proposed algorithm. In addition, according to the method presented in [31], the proposed method has a higher speed, while its accuracy has not changed significantly.

References [1] B. Ning, X. Gao, Multi-frame image super-resolution reconstruction using sparse co-occurrence prior and sub-pixel registration, Neurocomputing 117 (2013) 128–137. [2] F. Sauer, Image registration: enabling technology for image guided surgery and therapy., in: 27th Annual International Conference of the IEEE on Engineering in Medicine and Biology Society, IEEE-EMBS 2005, 2005, 2006, pp. 7242–7245. [3] F. Khalifa, G.M. Beache, G. Gimel’farb, J.S. Suri, A.S. El-Baz, State-of-the-art medical image registration methodologies: a survey, in: Multi Modality State-of-the-Art Medical Image Segmentation and Registration Methodologies, Springer US, 2011, pp. 235–280. [4] B. Zitova, J. Flusser, Image registration methods: a survey, Image Vision Comput. 21 (11) (2003) 977–1000. [5] A.A. Goshtasby, Image Registration: Principles, Tools and Methods, Springer Science & Business Media, 2012. [6] J. Ma, J.C. Chan, F. Canters, Fully automatic subpixel image registration of multiangle CHRIS/Proba data, IEEE Trans. Geosci. Remote Sens. 48 (7) (2010) 2829–2839. [7] R.J. Zhang, J.Q. Zhang, C. Yang, Image registration approach based on SURF], Infrared Laser Eng. 1 (2009) 041.

[8] D.J. Holtkamp, A.A. Goshtasby, Precision registration and mosaicking of multicamera images, IEEE Trans. Geosci. Remote Sens. 47 (10) (2009) 3446–3455. [9] C.A. Shah, Y. Sheng, L.C. Smith, Automated image registration based on pseudoinvariant metrics of dynamic land-surface features, IEEE Trans. Geosci. Remote Sens. 46 (11) (2008) 3908–3916. [10] D. Shen, C. Davatzikos, HAMMER: hierarchical attribute matching mechanism for elastic registration, IEEE Trans. Med. Imaging 21 (11) (2002) 1421–1439. [11] S. Liao, A. Chung, Feature based nonrigid brain MR image registration with symmetric alpha stable filters, IEEE Trans. Med. Imaging 29 (1) (2010) 106–119. [12] S. Klein, M. Staring, K. Murphy, M.A. Viergever, J.P.J.P. Pluim, Elastix: a toolbox for intensity-based medical image registration, IEEE Trans. Med. Imaging 29 (1) (2010) 196–205. [13] A. Myronenko, X. Song, Intensity-based image registration by minimizing residual complexity, IEEE Trans. Med. Imaging 29 (11) (2010) 1882–1891. [14] Q. Li, H. Ji, Multimodality image registration using local linear embedding and hybrid entropy, Neurocomputing 111 (2013) 34–42. [15] Z. Yi, S. Soatto, Nonrigid registration combining global and local statistics, in: IEEE Conference on Computer Vision and Pattern Recognition, CVPR 2009, IEEE, 2009, pp. 2200–2207. [16] G. Hermosillo, O. Faugeras, Dense image matching with global and local statistical criteria: a variational approach, in: Proceedings of the 2001 IEEE Computer Society Conference on Computer Vision and Pattern Recognition, CVPR 2001, vol. 1, I-73, IEEE, 2001. [17] D. Loeckx, P. Slagmolen, F. Maes, D. Vandermeulen, P. Suetens, Nonrigid image registration using conditional mutual information, IEEE Trans. Med. Imaging 29 (1) (2010) 19–29. [18] S. Klein, U.A. van der Heide, I.M. Lips, M. van Vulpen, M. Staring, J.P. Pluim, Automatic segmentation of the prostate in 3D MR images by atlas matching using localized mutual information, Med. Phys. 35 (4) (2008) 1407–1417. [19] P.P. Wyatt, J.A. Noble, MAP MRF joint segmentation and registration., in: Medical Image Computing and Computer-Assisted Intervention – MICCAI, Springer Berlin Heidelberg, 2002, pp. 580–587. [20] G.G. Zheng, X. Zhang, A unifying MAP-MRF framework for deriving new point similarity measures for intensity-based 2D-3D registration., in: 18th International Conference on Pattern Recognition, ICPR 2006,vol. 2, IEEE, 2006, pp. 1181–1185.

112

K. Aghajani et al. / Biomedical Signal Processing and Control 49 (2019) 96–112

[21] D. Mahapatra, Y. Sun, MRF-based intensity invariant elastic registration of cardiac perfusion images using saliency information, IEEE Trans. Biomed. Eng. 58 (4) (2011) 991–1000. [22] M.P. Heinrich, M. Jenkinson, M. Bhushan, T. Matin, F.V. Gleeson, M. Brady, J.A. Schnabel, MIND: Modality independent neighbourhood descriptor for multi-modal deformable registration, Med. Image Anal. 16 (7) (2012) 1423–1435. [23] Z. Li, D. Mahapatra, J.A. Tielbeek, J. Stoker, L.J. van Vliet, F.M. Vos, Image registration based on autocorrelation of local structure, IEEE Trans. Med. Imaging 35 (1) (2016) 63–75. [24] C. Wachinger, N. Navab, Entropy and Laplacian images: structural representations for multi-modal registration, Med. Image Anal. 16 (1) (2012) 1–17. [25] A. Wong, J. Orchard, Robust multimodal registration using local phase-coherence representations, J. signal Process. Syst. 54 (1-3) (2009) 89–100. [26] K. Aghajani, R. Yousefpour, M. Shirpour, M.T.M.T. Manzuri, Intensity based image registration by minimizing the complexity of weighted subtraction under illumination changes, Biomed. Signal Process. Control 25 (2016) 35–45. [27] A. Ghaffari, E. Fatemizadeh, Sparse-induced similarity measure: mono-modal image registration via sparse-induced similarity measure, Image Process. IET 8 (12) (2014) 728–741. [28] A. Ghaffari, E. Fatemizadeh, Robust Huber similarity measure for image registration in the presence of spatially-varying intensity distortion, Signal Process. 109 (2015) 54–68. [29] M.M. Fouad, R.M. Dansereau, A.D. Whitehead, Geometric image registration under arbitrarily-shaped locally variant illuminations, Signal Image Video Process. 6 (4) (2012) 521–532. [30] L. Wang, C. Pan, Nonrigid medical image registration with locally linear reconstruction, Neurocomputing 145 (2014) 303–315. [31] K. Aghajani, M.T. Manzuri, R. Yousefpour, A robust image registration method based on total variation regularization under complex illumination changes, Comput. Methods Programs Biomed. 134 (2016) 89–107. [32] J. Modersitzki, S. Wirtz, Combining homogenization and registration, in: Biomedical Image Registration, Springer Berlin Heidelberg, 2006, pp. 257–263. [33] S.L. Keeling, Image similarity based on intensity scaling, J. Math. Imaging Vision 29 (1) (2007) 21–34.

[34] M. Ebrahimi, A.A. Lausch, A.L. Martel, A Gauss-Newton approach to joint image registration and intensity correction, Comput. Methods Programs Biomed. 112 (3) (2013) 398–406. [35] M. Urschler, M. Werlberger, E. Scheurer, H. Bischof, Robust optical flow based deformable registration of thoracic CT images, Med. Image Anal. Clinic: A Grand Challenge (2010) 195–204. [36] C. Frohn-Schauf, S. Henn, K. Witsch, Multigrid based total variation image registration, Comput. Visual. Sci. 11 (2) (2008) 101–113. [37] W. Hu, Y. Xie, L. Li, W. Zhang, A total variation based nonrigid image registration by combining parametric and non-parametric transformation models, Neurocomputing 144 (2014) 222–237. [38] L.I. Rudin, S. Osher, E. Fatemi, Nonlinear total variation based noise removal algorithms, Physica D 60 (1) (1992) 259–268. [39] C. Hogea, C. Davatzikos, G. Biros, Medical Image Computing and Computer-Assisted Intervention – MICCAI, 2007. [40] M. Rajchl, J.S. Baxter, W. Qiu, A.R. Khan, A. Fenster, T.M. Peters, J. Yuan, Fast deformable image registration with non-smooth dual optimization, Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition Workshops (2016) 25–32. [41] T.F. Chan, G.H. Golub, P. Mulet, A nonlinear primal-dual method for total variation-based image restoration, SIAM J. Sci. Comput. 20 (6) (1999) 1964–1977. [42] M. Zhu, T. Chan, An efficient primal-dual hybrid gradient algorithm for total variation image restoration, UCLA CAM Report (2008) 08–34. [43] A. Chambolle, V. Caselles, D. Cremers, M. Novaga, T. Pock, An introduction to total variation for image analysis, Theoretical foundations and numerical methods for sparse recovery 9 (263-340) (2010) 227. [44] R. Szeliski, J. Coughlan, Spline-based image registration, Int. J. Comput. Vision 22 (3) (1997) 199–218. [45] M. Khader, A.B. Hamza, Nonrigid image registration using an entropic similarity, IEEE Trans. Inf. Technol. Biomed. 15 (5) (2011) 681–690. [46] J. Chen, J. Tian, N. Lee, J. Zheng, R.T. Smith, A.F. Laine, A partial intensity invariant feature descriptor for multimodal retinal image registration, IEEE Trans. Biomed. Eng. 57 (7) (2010) 1707. [47] L. Chen, X. Huang, J. Tian, Retinal image registration using topological vascular tree segmentation and bifurcation structures, Biomed. Signal Process. Control 16 (2015) 22–31. [48] W. Liu, E. Ribeiro, A survey on image-based continuum-body motion estimation, Image Vision Comput. 29 (8) (2011) 509–523.