Information Sciences 372 (2016) 492–504

Contents lists available at ScienceDirect

Information Sciences journal homepage: www.elsevier.com/locate/ins

Removing mismatches for retinal image registration via multi-attribute-driven regularized mixture model Gang Wang a,b,∗, Zhicheng Wang b,∗, Yufei Chen b, Qiangqiang Zhou b, Weidong Zhao b a b

School of Statistics and Management, Shanghai University of Finance and Economics, Shanghai 200433, China CAD Research Center, College of Electronic and Information Engineering, Tongji University, Shanghai 201804,China

a r t i c l e

i n f o

Article history: Received 8 September 2015 Revised 31 July 2016 Accepted 14 August 2016 Available online 16 August 2016 Keywords: Retinal image registration Point set registration Mixture model Kernel method

a b s t r a c t In order to address the problem of retinal image registration, this paper proposes and analyzes a novel and general matching algorithm called Multi-Attribute-Driven Regularized Mixture Model (MAD-RMM). Mismatches removal can play a key role in image registration, which refers to establish reliable matches between two point sets. Here the presented approach starts from multi-feature attributes which are used to guide the feature matching to identify inliers (correct matches) from outliers (incorrect matches), and then estimates the spatial transformation. In this paper, motivated by the problem of feature matching that the initial correspondence is always contaminated by outliers, thereby we formulate this issue as a probability deformable mixture model which consists of Gaussian components for inliers and uniform components for outliers. Moreover, the algorithm takes full advantage of using multiple attributes for better general matching performance. Here we are assuming all inliers are mapped into a high-dimensional feature space, namely reproducing kernel Hilbert space (RKHS), and the closed-form solution to the mapping function is given by the representation theorem with L2 norm regularization under the ExpectationMaximization (EM) algorithm. Finally, we evaluate the performance of the algorithm by applying it to retinal image registration on several datasets, where experimental results demonstrate that the MAD-RMM outperforms current state-of-the-art methods and shows the robustness to outliers on real retinal images. © 2016 Elsevier Inc. All rights reserved.

1. Introduction Image registration plays an important role in computer vision, pattern recognition, and medical diagnosis [43]. The task of this registration problem is to learn correspondences between two sets of features and to update the underlying transformation. More precisely, two or more images can be aligned together and mapped into the same scene with the optimal transformation estimated from the learned correspondence. As known widely, image registration is a fundamental problem and is still challenging in the field of vision-based research. From the perspective of the matching correspondence, image registration can be classified into two categories: areabased or feature-based methods. Intuitively, area-based methods use a specific similarity metric such as mutual information

∗

Corresponding authors. E-mail addresses: [email protected] (G. Wang), [email protected] (Z. Wang).

http://dx.doi.org/10.1016/j.ins.2016.08.041 0020-0255/© 2016 Elsevier Inc. All rights reserved.

G. Wang et al. / Information Sciences 372 (2016) 492–504

493

Fig. 1. Feature point matching of retinal image pairs. (a) (e) Retinal image pair (multi-modality). (b) (f) Feature matching, where yellow lines denote inliers and black lines denote outliers. (c) (g) Vector field. Inliers (blue arrows) and outliers (black arrows) are sampled from the matching result. Vector field interpolation [32] using line integral convolution (LIC) [7]. (d) (i) Registered images from (b) and (f) respectively, but (d) is with an obvious ghost. Best viewed in color. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

(MI) [27,34,40], cross correlation (CC) [12], entropy correlation coefficient (ECC) [9], and phase correlation [22], to get the correspondence between two images. Though these methods can get good registration results, they always suffer from intractable computational complexity, image distortion and illumination changes. By contrast, feature-based methods extract salient local keypoint including Harris corner [19], Scale Invariant Feature Transform (SIFT) [30], and Speed Up Robust Feature (SURF) [4], to match correspondences. Typically, sparse feature point information is used to describe the whole image, and it is easy to estimate the transformation with tractable computational complexity. In this paper, we mainly focus our study on the feature-based methods for retinal image registration. Feature-based image registration has in common with robust feature matching which needs to construct an initial alignment between local feature points. However, the existing uncertain mismatches make the matching task more difficult. In view of the transformation estimation needing a reliable feature correspondence, we can resort to a robust mismatches removal method. For instance, Fig. 1 shows the matching problem clearly. Initial feature matches are shown in Fig. 1(b), where mismatches make the registered image give rise to an obvious ghost in Fig. 1(d), while Fig. 1 (i) shows the perfect image registration after removing mismatches. Furthermore, vector field interpolation shows the smooth vector field learning from the feature matching. For retinal image registration, however, multi-modal image data, which is captured by different imaging sensors such as red-free and fluorescein angiography, is challenging to deal with. It is easy to extract salient feature points by SIFT or SURF descriptors for mono-modal image registration, while multi-modal image registration is difficult to construct a reliable descriptor. Though bifurcation-based methods are widely used for multi-modal retinal image registration, the vascular tree is hard to extract and segment in complex scenes such as poor-quality, occlusion by hemorrhage, and unhealthy area [1,10]. In order to address the above problems, then we propose a novel method called multi-attribute-driven regularized mixture model (MAD-RMM). The method can be used to learn the underlying correspondence efficiently when facing mismatches. More precisely, the initial matching correspondence is always contaminated by mismatches, let us define the correct matches and mismatches to be inliers and outliers, respectively. Subsequently, we use a Gaussian distribution and a uniform distribution to model inliers and outliers, respectively, then both of them can be absorbed into a mixture model. Under this mixture model, we make use of a novel multi-attribute-driven approach to identify inliers. In the MAD-RMM, a Maximum A Posteriori (MAP) solution of the unknown parameters of the transformation can be solved by the Expectation Maximum (EM) algorithm [13] where E-step computes the responsibilities and M-step updates the transformation, thus the underlying correspondence can be recovered by computing the maximum expectation after the algorithm convergence. Moreover, we add an L2 norm regularization term to constraint the non-linear transformation preserving smooth in a reproducing kernel Hilbert space (RKHS) [2,37] and avoiding the ill-posed problem in the optimization procedure. In many retinal image registration applications, though high-order transformation models show accurate results, we found that the linear transformation model such as affine [21], can fit most experimental cases, so we can apply the MAD-RMM algorithm with an affine model to register retinal images. Our contribution in this paper includes the following two aspects: (1) we propose a regression method to learn the correspondence for the feature-based matching. The initial correspondence is mapped to a high dimensional space by a specific kernel method, and the proposed learning approach can identify inliers from outliers in the special feature space (RKHS). (2) we introduce a multi-attribute-driven technique which can be efficient to guide the mixture model to recover the underlying inliers. Finally, in our previous work, we have presented the SURF-PIIFD [49] which improves the descriptor PIIFD (partial local intensity invariant feature descriptor) [10] for multi-model retinal image registration, and the SURF-PIIFD can be used to capture more reliable feature points to construct the initial correspondence.

494

G. Wang et al. / Information Sciences 372 (2016) 492–504

The remainder of the paper is organized as follows. Section 2 describes background material and related work. In Section 3, we present our multi-attribute-driven regularized mixture model. In Section 4, we apply the proposed method to register retinal images. Section 5 describes the implementation details. In Section 6, we illustrate the registration performance of the proposed method on various types of retinal image pairs with comparisons to other state-of-the-art methods. In Section 7, we give a brief discussion and conclusion. 2. Related work In this paper, we focus on digital retinal image registration which is widely used to diagnose varieties of diseases, including diabetic retinopathy, glaucoma, and age-related macular degeneration [42,58]. As such, numerous algorithms have been presented for image registration [1,3,8,10,26,44,47,49,54,55]. They aim to recover the underlying correspondence or find the robust estimation of transformation to align images. Here, we briefly overview these algorithms according to the two major categories: area-based and feature-based methods. 2.1. Area-based methods Area-based methods directly deal with the retinal image intensity values without attempting to detect salient structural information. The goal of these methods is to choose a certain similarity measure approach to compare the overlap area and then find an optimal transformation function to match and register an image pair. Mutual Information (MI) [27,34,40], Cross-Correlation (CC) [12], Entropy Correlation Coefficient (ECC) [9], and Phase Correlation (PC) [22] are the widely used methods. In this case, it will be feasible to formulate the registration problem and get good results in uncomplicated scenarios. However, due to using the dense intensity information, optimization algorithms (e.g., simulated annealing [23], genetic algorithms [17]) which are applied to minimize the measures of matching similarity, may fall into high computational complexity, where the huge searching space depends on the complexity of the selected transformation model. Note that the similarity measure is always disturbed by non-overlapping area caused by large camera view changes. Moreover, the performance of these area-based methods degrades when confronting with high illumination, content, background, or texture changes. 2.2. Feature-based methods In contrast to area-based methods, feature-based methods work directly with a sparse structural feature such as line or point-feature extracted from images, where the extracted feature is mainly used to represent the invariance between images. Image registration thereby will be reduced to the feature-based matching problem after extracting salient features. Thus the main objective of this registration is to learn the underlying correspondence from initial matching results which are contaminated by many outliers, and then find an optimal spatial transformation function between feature point sets under an iterative optimization. Briefly, the feature-based method consists of three components: design the invariant feature descriptor, identify the correct correspondence and fit the optimal transformation function. For instance, Liao et al. [29] have presented a feature guided deformable image registration framework which can be applied for face recognition. In order to analyze medical images, Wu et al. [53] have proposed a hierarchical attribute-guided registration method. Moreover, Kim et al. [25] have found that a sparse patch-based deformation estimation can be used for image registration. For retinal image registration, several features including bifurcation, fovea, optic disc, and corner [8,14,20,26,35,39,55,56] are widely used, where bifurcation (i.e., Y-shape feature), fovea and optic disc, from the view of the retinal structure, are invariant features to intensity, scale, rotation, and illumination changes. More precisely, assuming that the vascular tree has been detected, and bifurcations have been labeled with surrounding vessel orientations, then the probability for each matching bifurcation pair is given by this angle-based invariant. However, it is difficult to extract bifurcations, fovea, and optic disc in poor quality and unhealthy retinal images. To overcome this limitation, SIFT [30] and SURF [4] descriptors are applied to extract features, but they are not appropriate for multi-modal retinal image registration. Further, Chen et al. have designed the PIIFD [10] which is an SIFT-based variant, can be well applied for multi-modal retinal registration. In addition, Lee et al. [1] have exploited a low-dimensional step pattern analysis (LoSPA) method for unhealthy retinal images, but they need to pre-design many step patterns. Feature matching is used to fit correct matches and remove mismatches from the initial correspondence. Intuitively, function fitting or regression strategy will be feasible to formulate this mismatch removal problem. Many existing methods are as follows: Random Sample Consensus (RANSAC) algorithm [16] is the most popular strategy, it is easy to design and obtain well performance, but it always suffers from some limitations when dealing with a non-linear transformation or wide-baseline images. In order to address the limitation, many improved methods have been proposed: non-rigid RANSAC [46] has been designed to handle the non-linear transformation. Li and Hu [28] use support vector regression (SVR) to identify the correct correspondence. Furthermore, non-parameter-based methods use a non-linear mapping function to fit and identify inliers from outliers in high-dimensional RKHSs, then the identified inliers can be interpolated as a smooth coherent vector field. Examples of the methods consist of Vector Field Consensus (VFC) [32,33], Robust Point Matching via L2 -minimizing Estimation (RPM-L2 E), Mixture of Asymmetric Gaussians model (MoAG) [48,52], and Context-Aware Gaussian

G. Wang et al. / Information Sciences 372 (2016) 492–504

495

fields (CA-LapGF) [51] for mismatches removal. Generally, these methods need to construct the initial correspondence by a k-nearest neighbors (k-NN) method, such as the best bin first (BBF) [30]. Moreover, heuristic methods use an affinity matrix to recover point set correspondence without describing local features, where the location information or point coordinate is directly used to measure the similarity between point sets by distance-based or correlation-based strategies. Up to now, iterated closest point (ICP) algorithm [6] is widely used for image registration with a linear transformation, so as to exploit several ICP-based retinal registration methods, such as dualbootstrap ICP (DB-ICP) [44], generalized dual-bootstrap ICP (GDB-ICP) algorithm [54], and Edge-driven dual-bootstrap ICP (ED-DB-ICP) [47]. More precisely, each bootstrap region of the DB-ICP has three steps: refining the transformation estimation, expanding the bootstrap region, and estimating higher-order transformation. The ED-DB-ICP algorithm combines SIFT keypoints and vascular features to register multi-modal retinal images. Although, ICP-based methods need only one initial correct match to align image pairs successfully after the optimization algorithm reaching convergence, while the performance degrades when facing with poor quality or unhealthy retinal images. In addition, a robust point matching framework with thin plate spline (RPM-TPS) [11] is established for non-rigid image registration, where the soft-assignment technique is used to recover the correspondence from the affinity matrix, and the non-rigid transformation is modeled by TPS. Under this framework, Jian and Vemuri [24] use the Gaussian mixture model (GMM) to represent each feature point set, and then estimate the non-rigid geometrical transformation by L2 -norm distance, then they have presented GMM-TPS algorithm for point set registration, but it is sensitive to outliers. To improve the robustness to outliers, the coherence point drift (CPD) algorithm [38] formulates the registration problem as a mixture model with considering both inliers and outliers, and then it uses a Gaussian radial basis function instead of using the TPS to model the transformation. Similarly, kernel density estimation and fuzzy correspondence based method [50] is used to register the contaminated point sets. Another point-set representation method called shape context (SC) descriptor [5] is well used for point set registration. However, registration performance of these methods will degrade when facing a large number of outliers in the initial correspondence or in the affinity correlation matrix. 3. Removing mismatches 3.1. Notation We use the following notation in our method. • Image pair. Let I1 and I2 be the fixed target image and the moving image respectively. We aim to register the moving image onto the target image in the registration problem. M • Point set pair. Let XN×2 = {(xn , sn , pn )}N n=1 and YM×2 = { (ym , tm , qm )}m=1 be the feature point set extracted from I1 and I2 respectively, where location inform xn , ym ∈ R2 , shape context descriptor sn , tm ∈ RD1 , and local feature descriptor pn , qm ∈ RD2 , D1 and D2 denotes the dimension of feature descriptor. • Correspondence. Let S = {(xl , yl )}Ll=1 be the correspondence set that estimated by the multi-attribute-driven model iteratively.  by a mapping function T , i.e., Y  = T (Y, θ ), where θ denotes the transfor• Mapping function. Point set Y is moving to Y Niter mation parameter. Note that T = {Ti }i=1 is a function set, where Ti is the recovered transformation in the ith iteration and Niter denotes the maximum iteration number. • Others. ◦ denotes the Hadamard product, and  denotes the tensor product. 0 denotes a column vector of all zeroes, 1 denotes a column vector of all ones, and diag( · ) denotes the diagonal matrix. 3.2. Problem formulation Given an image pair I1 and I2 , the goal of registration problem is to align the moving image onto the target image. Due to both X and Y consist of multi-attribute information which is extracted from the given image pair, so feature pattern can be used to simplify this registration problem. Typically, location information is widely used to estimate the correspondence and transformation in point set registration. In order to make a robust registration method, a regularized Gaussian fields criterion [31] which combines shape context descriptor, can be used to register images. Intuitively, inlier points should be identified perfectly from the point set which is contaminated by outliers. To do this, according to the probabilistic deformable model, we assume the inlier satisfies a Gaussian distribution N (xn ; ym , θ , σ 2 ) with a range parameter σ , while the outlier satisfies a uniform distribution N1 . Motivated by the feature-guided strategy, multi-attribute information is joint to drive the probabilistic model and get more robust performance for registration, because it might appear that this strategy guides the weighting of the Gaussian component of this model. Thus, the registration problem can be formulated as a mixture model which takes the form as follows

p(xn |θ ) = γ1

M  m=1

Vnm N (xn ; ym , θ , σ 2 ) + γ

1 , N

(1)

where 0 ≤ γ 1 ≤ 1, 0 ≤ γ ≤ 1, γ1 + γ = 1, and Vnm = ωnm ◦ νnm denotes the joint weight of the inlier Gaussian component. The similarity measure between point pair with shape context [5] and a certain local feature descriptor are ωnm = exp(−C (sn , tm , α1 )) and νnm = exp(−H ( pn , qm , α2 )), respectively.

496

G. Wang et al. / Information Sciences 372 (2016) 492–504

Fig. 2. The flowchart of estimating correspondence via multi-attribute-driven mixture model.

The joint probability distribution (likelihood function L) for the whole data point set X is defined as follows

p(X |Y, θ , σ 2 ) =

N 

p( x n | θ ) = L ( θ | X ) ,

(2)

n=1

where the data point set satisfies i.i.d. (independent identically distributed). Based on the Gaussian mixture model (GMM) and Bayes’ rule, a MAP solution of θ can be given by minimizing the following negative log-likelihood function

E1 (θ ) = − log L(θ |X ) = −

N 

log

n=1

M 

p( x n | θ ) .

(3)

m=1

3.3. Multi-attribute-driven strategy 1 There are mainly two motivations: (1) in a traditional standard mixture model, such as CPD [38], where Vnm = M denotes that each point will have the equal opportunity to move towards its according target point. We thereby consider that each inlier Gaussian component should be assigned a large weight which can be calculated by multi-feature measuring. (2) for 2D retinal images, local intensity invariant feature can be easily described by PIIFD [10], and point positions can be described by SC [5]. Intuitively, the least square method is widely used to measure each point pair based on the location information (points coordinates) under L2 distance. By contrast, Gaussian distance between a pair of moving and target points is more robust, as demonstrated in [15,31]. Thus we can measure the global distance with location information as follows



 xn − T ( ym , θ )2 N (xn ; ym , θ , σ ) = α0 exp − , 2σ 2 2

(4)

where α0 = 2π1σ 2 is a positive value coefficient. From the view of shape similarity, it can be measured by SC [5] which denotes a feature descriptor for 2D point set. Given a point pair (xn , ym ), and its shape context is represented as (sn , tm ), then the cost of matching the point pair Cnm = C (sn , tm ) using the χ 2 test statistic is defined as follows

C (sn , tm , α1 ) = α1

K  (sn (k ) − tm (k ))2 , sn ( k ) + tm ( k )

(5)

k=1

where α1 = 0.5 is a constant. sn (k) and tm (k) denote the K-bin normalized histogram at xn and ym , respectively. Moreover, the local intensity invariant feature descriptor can be applied to construct the multi-attribute, such as SIFT [30] and SURF [4] for mono-modal images, while PIIFD [10] and SURF-PIIFD [49] for multi-modal retinal images. Then a bipartite matching method δ ( ·, ·, ·) (BBF [30]) is used to measure the feature similarity, so we can define it as follows

H ( pn , qm , α2 ) = δ ( pn , qm , α2 ),

(6)

where α2 = 1.5 is a positive value which denotes the matching threshold. For our proposed multi-attribute-driven probabilistic deformable model, as shown in Fig. 2, joint information thereby is used to guide the mixture model for estimating inliers iteratively.

G. Wang et al. / Information Sciences 372 (2016) 492–504

497

3.4. Transformation model and regularization In the deformable mixture model, the displacement function F is non-rigid for accurate point set registration. The moving points are mapped to the target point set in a reproducing kernel Hilbert space (RKHS). Radial basis function (RBF), including Gaussian RBF and Thin-plate splines (TPS), is used to construct the kernel in an RKHS H. In this paper, we use the y −y GRBF with ki j = exp(− 21  i β j 2 ) to construct a vector-valued kernel matrix K with element K (yi , y j ) = ki j . The non-rigid registration process is a regression problem which needs a regularization to avoid overfitting, then the Tikhonov regularization can be employed to mitigate this effect. According to the regularization theory [45], the new negative log-likelihood function with regularization term φ (F ) can be rewritten as

E ( θ ) = E1 ( θ ) +

λ 2

φ ( F ),

(7)

where λ ≥ 0 is a trade-off parameter and can be used to tune the empirical risk and the regularization penalty. Considering the smoothness constraint [38] on the displacement field, the regularization term comes from a prior as follows



p(F ) ∝ exp −

λ 2

 φ (F ) ,

(8)

where φ (F ) = F 2H , and  · 2H denotes the norm in the RKHS H. According to Bayes’ theorem, the posterior distribution p(θ |X, Y, σ 2 ) could be estimated by the given likelihood (2) and prior (8)

p(θ |X, Y, σ 2 ) ∝ p(X |Y, θ , σ 2 ) p(F ),

(9)

then the optimal parameter θ  can be estimated by the equation

θ  = arg max p(θ |X, Y, σ 2 ).

(10)

θ

3.5. The solution based on EM algorithm In order to seek an optimal solution of (10), in this paper, we use the expectation-maximization (EM) algorithm (a two-stage iterative optimization technique) [13] to solve the regularized mixture model. Considering the latent variants Z = {zi(n,m) = 1, zo(n,m) = 0}N,M in the complete data, in the E-step, we need to compute a posteriori probability distribution n=1,m=1

old of the mixture components, where P old (z ) for inliers, while P old (z ) for outliers. In the subsequent M-step, parameters Pnm nm i nm o F, σ 2 , and γ can be estimated and updated iteratively. We omit the terms independent of parameters  = {F, σ 2 , γ }, then the objective function can be rewritten as follows

Q(, old ) =

N M 1 

2σ 2

P old (zi )xn − (ym + F (ym ))2 +

n=1 m=1

 D log σ 2 P old (zi ) 2 N

M

n=1 m=1

− log (1 − γ )

N  M 

P

old

(zi ) − log γ

n=1 m=1

N  n=1

P

old

( zo ) +

λ 2

F 

(11)

2 H,

where old denotes the current estimated parameters which are used to evaluate the expectation and  denotes the new parameters that we optimize to increase Q. In the E-step, the posterior probabilities of mixture of Gaussian components can be calculated with the estimated values of parameters under the Bayes’ rule

ωnm νnm N (xn ; ym , θ , σ 2 ) , 2 k=1 ωnk νnk N (xn ; yk , θ , σ ) + c

Pnm = M

(12)

γM

where c = γ N is a constant. P (xn , ym ) = Pnm denotes the matching score value between xn and ym . Note that, in both spatial 1 and feature space, the closer distance between xn and ym , the larger value Pnm has. In the M-step, the optimal parameter group new = arg max Q(, old ), thus we can update parameters using the derivative of Q(, old ) with respect to . We obtain

γ1 = σ = 2

N

n=1

M

m=1 Pnm

N N

n=1

M

,



(13)

(

m=1 Pnm xn − ym +   2 Nn=1 M m=1 Pnm

F (ym ))2

.

(14)

498

G. Wang et al. / Information Sciences 372 (2016) 492–504

Due to the derivative of Q with respect to F is complex to solve [57]. Motivated by the literature [38], the functional form of displacement function F is solved using calculus of variation in the special functional space H. Following the vectorvalued representation theorem [36] in a Gaussian kernel based RKHS, we can obtain the optimal function

F (· ) =

M 

K (·, yk )Ok ,

(15)

k=1

where Ok is the coefficient needed to solve. Based on the form (15), the objective function (11) can be rewritten as

2 M  λ Q (F ) = Pnm xn − (ym + K (ym , yk )Ok + trace(OT KO ). 2 2 2σ n=1 m=1 k=1 N M 1 

(16)

Taking the derivative of (16) with respect to the coefficient row vector O = (O1 , . . . , OM ), we can get O by solving the following linear system

(K + λσ 2 diag(P)−1 )O = diag(P)−1 PX − Y,

(17)

= P where P N×M · 1M×1 , and the matrix P with elements P (n, m ) = Pnm . Then the mapping function can be represented as T (Y ) = Y + KO. According to the soft-assignment matching algorithm [41], we can get the correspondence set S = {n, m|Pnm ≥ τ } with a given threshold τ after the EM procedure converges. The proposed multi-attribute-driven regularized probabilistic mixture model (namely MAD-RMM) for robust feature matching is summarized in Algorithm 1. Algorithm 1: The MAD-RMM algorithm. Input: The target point set X, and the moving point set Y . , and the correspondence set S. Output: The transformed point set Y Initialize: σ 2 , β , λ, γ , τ , and O = 0. The kernel matrix K based on a Gaussian kernel with parameter β ; while iteration < termination condition do E-step: Compute the responsibility P by Eqs. (4), (5), (6), and (12); M-step: Update the parameters γ , σ 2 , and O using Eqs. (13), (14), and (17) respectively;  = Y + KO; Update the moving point set: Y end  by the optimal parameter O; Determine Y Determine S by the responsibility P and the threshold τ ; , and S; return Y

4. Application to retinal image registration In the MAD-RMM algorithm, a smooth mapping function from moving point set to the target point set is solved to recover the underlying correspondence and align them together. Generally, the proposed algorithm can be applied to 2D image matching with specific feature detectors and descriptors. In this paper, we focus on the retinal image registration, particularly the multi-modal data, and apply the MAD-RMM algorithm to retinal image registration. 4.1. SURF-PIIFD The PIIFD descriptor [10] is based on two assumptions: 1) the similar anatomical structure regions of one image would consist of similar outlines in the corresponding regions of another image, and 2) the gradient orientations at the corresponding control point locations would point to the same or opposite directions in the multi-modal retinal image. Thus, the descriptor is useful to extract repeatable, rotation, scale invariant, intensity, and affine partial invariant local features. Although the Harris corner detector is widely used for local feature matching, it is not scale-invariant and has low feature repeatability. Fortunately, the robust scale invariant feature, speed up robust feature (SURF) [4] has been designed, where the interest-point candidates can be detected by the maximum determinant of the Hessian matrix. Using the integral image, SURF approximates the different levels of scale space by adjusting the size of the box filters instead of the original image as used in SIFT. Finally, interest points can be found by non-maximum suppression in a 3 × 3 × 3 neighborhood around each sample point and the second-order Taylor expansion around the keypoint is used to improve the locations of the key points. As such, we have improved the PIIFD in extracting more reliable and repeatable keypoints.

G. Wang et al. / Information Sciences 372 (2016) 492–504

499

4.2. Transformation estimation It is worth noting that there are many parameter models are applied to retinal image registration, including translation, Euclidean, similarity, affine, bilinear, and quadratic as discussed in [44], here we choose lower-order transformations (affine or quadratic) for retinal images which have low resolution and high overlap rate. Given the correspondence set S = {(xl , yl )}Ll=1 estimated by the MAD-RMM, and the transformation

T () = yˆl , 0 where  = (φ φ 6

(18)

φ1 φ7

φ2 ), yˆl = yl for affine transformation, and φ8

 = (φφ0

6

φ1 φ7

φ2 φ8

φ3 φ9

φ4 φ10

φ5 ), yˆl = (1, yl0 , yl1 , yl0 φ11

◦yl1 , y2l0 , y2l1 )T for second-order polynomial transformation, thus the new objective function for retinal image registration can be rewritten as follows

Q(, σ 2 ) =

L 1 

2σ 2

xl − yˆl 2 + L log σ 2 .

(19)

l=1

We can directly take the partial derivatives of Q(, σ 2 ) with respect to each parameter, and then we can obtain the solution by solving the linear system of that equations. Removing mismatches for retinal image registration via the MAD-RMM algorithm is summarized in Algorithm 2. Algorithm 2: Retinal image registration. Input: The target image I1 , and the moving image I2 . Output: he transformed image  I2 , and the fusion image I f . Initialize: Construct feature point set X, Y using SC and SURF-PIIFD; while iteration < termination condition do E-step: Compute the responsibility P ; M-step: Update the parameters γ , σ 2 , and O;  = Y + KO; Update Y end Compute the transformation parameters by solving the linear systems after taking the partial derivatives of Eq. (19) with respect to σ 2 , ; Determine  I2 by  using the bilinear interpolation; Determine I f by aligning I1 and  I2 image pairs; return  I2 , and I f ;

5. Implementation details Scale parameter σ 2 controls the width of capture range for each Gaussian mixture model, and we initialize it as

σ2

N

M

x −y 2

n m = n=1 m2=1 . GRBF parameter β controls the width of the Gaussian kernel, and it produces locally smooth transMN formation and globally translation transformation with setting small and large values, respectively. In the experiments, we set β = 3. Regularization parameter λ trades off the empirical risk and the smoothness regularization term, here we set λ = 8. Outlier weight parameter γ is used to estimate the contamination degree of the initial correspondence, here we first estimate the value of γ , and then redo the EM algorithm using the fixed γ until convergence, so we initialize γ = 0.9. The maximum iteration number of EM algorithm is set to 100. It is worth noting that Harris-PIIFD [10] speeds up the feature detector process by reducing the scale of retinal images, while we detect features by SURF-PIIFD on the original images to keep more local feature information [49]. Moreover, inliers can be identified from outliers by the given threshold τ = 0.5.

6. Experiments and results In this section, we evaluate the performance of our proposed MAD-RMM algorithm against to four state-of-the-art methods: RANSAC [16], SVR [28], CPD [38], and GMM-TPS [24]. All the experiments were performed in Matlab 2015a on a PC with a 2.5 GHz Intel i5 Core system and 8GB RAM, where all the methods’ parameters are fixed throughout these experiments. 6.1. Data sets In this paper, we have employed six retinal image data sets for evaluating the proposed method. We have collected six data sets from the Internet and our local hospital. The first two data sets are mono-modal, while the other four data sets are multi-modal. More details are described as follows.

500

G. Wang et al. / Information Sciences 372 (2016) 492–504 Table 1 Min, max, median, and mean inlier ratios (%) of the six retinal image data sets are determined by ground-truth.

Min Max Median Mean

A

B

C

D

E

F

13.41% 77.71% 59.25% 54.85%

34.87% 70.00% 53.07% 56.27%

10.00% 81.74% 31.43% 35.57%

32.14% 85.71% 62.96% 59.90%

44.12% 94.44% 68.31% 71.35%

10.18% 59.46% 26.83% 30.22%

Data set A consists of a total of 19 pairs of color retinal images, which are collected from the public STARE1 (STructured Analysis of the REtina). Each image resolution is 700 × 605 pixels. Data set B, i.e., DRFDSIR2 (Digital Retinal Fundus Data Set for Image Registration) contains 8 retinal image pairs which are captured by Beijing Tongren Hospital. The photos were digitized to 1548 × 1260 pixels. Data set C consists of 9 pairs of intensity retinal images, the original images are collected from the test data of i2k retina software package3 . The resolution of each image is 1024 × 1024 pixels. Data set D contains 13 retinal image pairs which are collected from the Internet. Each image pair consists of one color fundus image and one FFA image. The image resolution of each image is from 445 × 332 pixels to 1068 × 829 pixels. Data set E consists of 10 low-quality retinal image pairs, which are collected from the Internet. The image resolution of each image is from 350 × 284 pixels to 500 × 514 pixels. Data set F consists of 75 retinal image pairs, obtained from Shanghai Xinshijie Eye Hospital. The images were required using a Heidelberg Retina Angiograph II (HRA2, Heidelberg Engineering, Germany) with a 30-degree field of view. The image resolution of each image is 770 × 770 pixels. 6.2. Evaluation metrics In this paper, we select four commonly used metrics to evaluate the performance of the comparison methods quantiT P , Precision: Pr = T P , Accuracy: Acc= T P+T N tatively, and they are Recall: Re= T P+ , and Root Mean Square Error (RMSE) of FN T P+F P T P+T N+F P+F N



registration: Er= K1 Kk=1 (xk −yk )2 , where TP denotes true positive (correctly identified inliers), FP denotes false positive (incorrectly identified inliers), TN denotes true negative (correctly identified outliers), and FN denotes false negative (incorrectly identified outliers). Accuracy and RMSE indicate the overall matching and registration performance, respectively. Note that we manually construct the matching ground-truth. More precisely, we use the SURF-PIIFD to extract feature points from retinal images, and then we carefully remove all mismatches (outliers) from the initial correspondence which is given by the BBF matching. Inlier ratios (min, max, median, and mean) of each dataset are shown in Table 1. 6.3. Results We evaluate the MAD-RMM against to two styles of matching methods: 1) feature matching methods with initial labels, such as RANSAC [16], and SVR [28], and 2) the methods without labels, such as CPD [38], and GMM-TPS [24]. These four state-of-the-art methods are widely used in feature matching and point set registration. Fig. 3 shows an illustration of the qualitative results on the six retinal image data sets using the MAD-RMM algorithm. The iterative procedure of the EM algorithm is used to illustrate the convergence speed of the proposed method. In the beginning, we obtain the initial correspondence (the first column) using the SURF-PIIFD [49] in Matlab, where the feature matching ratio α2 is set to 1.50. From top to bottom, the average outlier ratio for each image data set is 45.15%, 43.73%, 64.43%, 40.10%, 28.65%, and 69.78%, respectively. Almost inliers are identified correctly from outliers after several iterations (less than 10). Then we can use the identified correct matches for retinal image registration under the affine transformation model. Statistical analysis results are shown in Fig. 4. Acc is used to indicate the proportion of true results among the total number of examined cases. The Acc of the MAD-RMM for all data sets is statistically significantly higher than that of the other comparison methods in most scenarios. Table 2 shows the average Re, Pr, and Acc for all data sets, and it illustrates that the proposed method gives a higher accuracy and a better balance between Re and Pr. Although the max outlier ratio reaches 90.00% in data set C, the proposed method gives the highest Acc value (above 0.81). We observe that SVR and CPD can give a good performance, where SVR is a supervised method with the labeled data, while CPD is an unsupervised method based on the clustering assumption. Fig. 5 plots the accuracy-recall curves as the metric used in [48], where we can see that CPD and SVR perform better than RANSAC and GMM-TPS, while our MAD-RMM is far better than the others. 1 2 3

http://www.ces.clemson.edu/∼ahoover/stare/. http://www.mitk.net/download_3dmed.html. http://www.dualalign.com/retinal/image-registration-montage-software-overview.php.

G. Wang et al. / Information Sciences 372 (2016) 492–504

501

Fig. 3. Qualitative feature matching results on several typical retinal image pairs by the proposed MAD-RMM. From top to bottom, the retinal image pair example comes from data set A-F respectively. The columns show the matching process during the EM iteration. In the middle five columns, black arrows denote the identified outliers, blue arrows denote the identified inliers, green arrows denote incorrectly identified inliers, and cyan arrows denote unidentified inliers. The leftmost column and the rightmost column denote the initial correspondence and the final matches, respectively. Best viewed in color. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 4. Performance comparison on the six retinal image data sets using accuracy. RANSAC [16], SVR [28], CPD [38], GMM-TPS [24], and MAD-RMM are tested for comparing the feature point matching performance on retinal image pairs. Best viewed in color.

From the perspective of point set registration, a challenging problem is that the data are contaminated by a large number of outliers, thereby the registration performance degrades significantly. Table 3 shows the retinal image registration error results, GMM-TPS is sensitive to outliers because all points (including outliers) are modeled as a Gaussian component, CPD performs well because of its formulated mixture model, while MAD-RMM gives the lowest alignment errors. More precisely, the mean ± standard deviation (STD) of the registration error in Table 3 demonstrates the stability of the MAD-RMM. It can be observed that the significant performance improvement from CPD to MAD-RMM on whole data sets. Table 4 shows the comparison of average computation time. Because the average number of the initial correspondence for these data sets are 330, 10 0 0, 120, 68, 38, and 67, respectively, so the elapsed time of the proposed method becomes large when dealing with a large number of feature points. In our proposed mixture model, low-rank matrix approximation

502

G. Wang et al. / Information Sciences 372 (2016) 492–504

Table 2 Comparison of matching results on the six retinal image data sets using recall (Re), precision (Pr), and accuracy (Acc). (Ds denotes Data sets, Me denotes evaluation metrics). Ds

Me

RANSAC

SVR

CPD

GMM-TPS

MAD-RMM

A

Re Pr Acc Re Pr Acc Re Pr Acc Re Pr Acc Re Pr Acc Re Pr Acc

0.5398 0.9831 0.5311 0.4831 0.8396 0.4034 0.2307 0.5040 0.1056 0.4567 0.7383 0.3145 0.5649 0.7804 0.4164 0.8840 0.2775 0.2517

1.0 0 0 0 0.6756 0.7377 1.0 0 0 0 0.7015 0.7688 0.9756 0.4076 0.4744 0.9835 0.6904 0.7145 0.8975 0.7656 0.7272 0.3392 0.9940 0.4106

0.8175 0.6289 0.8119 0.7445 0.6874 0.7585 0.6618 0.6256 0.5303 0.9010 0.5972 0.6775 0.9500 0.6670 0.7296 0.2741 0.4878 0.6762

0.7339 0.1856 0.5037 0.8386 0.2975 0.6185 0.6111 0.4729 0.4341 0.6948 0.6433 0.4280 0.8550 0.3702 0.3565 0.7389 0.3249 0.4323

0.9366 0.8652 0.9251 0.9187 0.6917 0.8248 0.9067 0.6625 0.8192 0.9478 0.7806 0.8195 1.0 0 0 0 0.6529 0.7683 0.5917 0.7806 0.8459

B

C

D

E

F

Fig. 5. Performances of removing mismatches methods on all the data sets using accuracy-recall curves. Best viewed in color.

Table 3 Comparison of point set registration error (pixel) of CPD, GMM-TPS, and our MAD-RMM on the six retinal image data sets.

CPD GMM-TPS MAD-RMM

A

B

C

D

E

F

2.61 ± 4.36 8.92 ± 10.26 0.82 ± 1.41

3.12 ± 4.09 14.10 ± 13.26 1.14 ± 0.68

14.30 ± 17.06 24.18 ± 19.04 2.59 ± 1.85

5.40 ± 5.94 23.01 ± 21.75 1.97 ± 1.61

5.49 ± 8.45 14.14 ± 7.20 2.75 ± 2.50

20.09 ± 13.64 47.96 ± 15.26 4.07 ± 4.29

Table 4 Comparison of matching times (sec.) on the six retinal image data sets. A RANSAC SVR CPD GMM-TPS MAD-RMM

3.62 0.25 1.21 5.26 1.68

B ± ± ± ± ±

0.14 0.15 0.96 3.32 0.81

4.46 2.34 15.42 35.81 21.42

C ± ± ± ± ±

0.62 1.63 14.39 22.31 6.66

3.42 0.04 0.32 1.78 0.24

± ± ± ± ±

0.13 0.05 0.28 0.90 0.17

D

E

3.26 ± 0.04 0.03 ± 0.01 0.16 ± 0.12 0.87 ± 0.28 0.09 ± 0.12

3.24 0.02 0.08 0.56 0.03

F ± ± ± ± ±

0.04 0.00 0.05 0.13 0.02

2.38 0.01 0.09 0.69 1.05

± ± ± ± ±

0.06 0.01 0.04 0.18 1.01

G. Wang et al. / Information Sciences 372 (2016) 492–504

503

[48] is applied to reduce the computation time of solving the linear system (17) easily, and the fast Gaussian transformation  and PX. (FGT) [18] is used to compute the matrix-vector products P 6.4. Comparison with retinal image registration methods Many algorithms have been proposed for retinal image registration, as summarized in Section 2.2. As seen in Table 3 and discussed in Section 6.3, the MAD-RMM substantially outperforms ICP-based methods, such as GDB-ICP [54] (registration success rate: {(A, 100%), (B, 100%), (C, 67%), (D, 0%), (E, 50%), (F, 25%)}), and ED-DB-ICP [47] (registration success rate: {(A, 100%), (B, 50%), (C, 78%), (C, 30.76%), (E, 40%), (F, 25%)}), while the registration success rate of MAD-RMM is {(A, 100%), (B, 100%), (C, 100%), (C, 100%), (E, 100%), (F, 100%)}. 7. Discussion and conclusion Typically, the initial matches computed by local descriptors contain mismatches because of the ambiguity of the local appearance information around the keypoints. Then, an optimization procedure is applied that uses the global position information to refine the matches and eliminate outliers. Mismatch removal plays a key role in image matching and registration, and the proposed MAD-RMM is a hybrid method which combines feature matching and spatial mapping correspondence estimation. Note that the idea behind the algorithm is to start from an initial correspondence with many uncertain outliers which are needed to identify by the multi-attribute-driven regularization mixture model. Outliers are precisely modeled in the mixture model, and a smooth spatial mapping function can be estimated in a particularly feature space RKHS by the EM algorithm iteratively. Thus the proposed algorithm is successfully applied to retinal image registration. It is worth noting that the method is robust and general, which is able to handle non-linear transformations when facing a large number of outliers, and hence can be applied to various registration tasks in retinal images. Moreover, the multi-attribute-driven strategy plays a role of a guidance on the feature matching. It is infeasible to design a custom feature representation, while it is easy to identify correct matches under several certain attributes including locations, shape context, and a specific local descriptor. Qualitative and quantitative experimental evaluation shows better performance of the proposed MAD-RMM algorithm on several types of retinal image : two mono-modal, and four multi-modal data sets. Finally, the experimental results demonstrate that our method significantly outperforms four tested state-of-theart methods for estimating correspondence and removing mismatches. Acknowledgment This work was partially supported by National Natural Science Foundation of China (NSFC 61103070), and the Fundamental Research Funds for The Central Universities. References [1] J. Addison Lee, J. Cheng, B. Hai Lee, E. Ping Ong, G. Xu, D. Wing Kee Wong, J. Liu, A. Laude, T. Han Lim, A low-dimensional step pattern analysis algorithm with application to multimodal retinal image registration, in: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, 2015, pp. 1046–1053. [2] N. Aronszajn, Theory of reproducing kernels, Trans. AM. Math. Soc. (1950) 337–404. [3] A. Bardera, M. Feixas, I. Boada, M. Sbert, Image registration by compression, Inf. Sci. 180 (7) (2010) 1121–1133. [4] H. Bay, T. Tuytelaars, L. Van Gool, Surf: Speeded up robust features, in: Computer Vision-ECCV 2006, Springer, 2006, pp. 404–417. [5] S. Belongie, J. Malik, J. Puzicha, Shape matching and object recognition using shape contexts, IEEE Trans. Pattern Anal. 24 (4) (2002) 509–522. [6] P.J. Besl, N.D. McKay, A method for registration of 3-d shapes, IEEE Trans. Pattern Anal. 14 (2) (1992) 239–256. [7] B. Cabral, L.C. Leedom, Imaging vector fields using line integral convolution, in: Proceedings of the 20th annual conference on Computer Graphics and Interactive Techniques, ACM, 1993, pp. 263–270. [8] A. Can, C.V. Stewart, B. Roysam, H.L. Tanenbaum, A feature-based, robust, hierarchical algorithm for registering pairs of images of the curved human retina, IEEE Trans. Pattern Anal. 24 (3) (2002) 347–364. [9] T. Chanwimaluang, G. Fan, S.R. Fransen, Hybrid retinal image registration, IEEE Trans. Inf. Technol. B. 10 (1) (2006) 129–142. [10] J. Chen, J. Tian, N. Lee, J. Zheng, R. Smith, A.F. Laine, A partial intensity invariant feature descriptor for multimodal retinal image registration, IEEE Trans. Biomed. Eng. 57 (7) (2010) 1707–1718. [11] H. Chui, A. Rangarajan, A new point matching algorithm for non-rigid registration, Comput. Vis. Image Und. 89 (2) (2003) 114–141. [12] A.V. Cideciyan, Registration of ocular fundus images, IEEE Eng. Med. Biol. 14(1) 52–58. [13] A.P. Dempster, N.M. Laird, D.B. Rubin, Maximum likelihood from incomplete data via the em algorithm, J. R. Stat. Soc. B., (1977) (1977) 1–38. [14] B. Fang, Y.Y. Tang, Elastic registration for retinal images based on reconstructed vascular trees, IEEE Trans. Biomed. Eng. 53 (6) (2006) 1183–1187. [15] B.A.M.A. Faysal Boughorbel Andreas Koschan, Gaussian fields: a new criterion for 3d rigid registration, Pattern Recognit. 37 (7) (2004) 1567–1571. [16] M.A. Fischler, R.C. Bolles, Random sample consensus: a paradigm for model fitting with applications to image analysis and automated cartography, Commun. ACM 24 (6) (1981) 381–395. [17] D.E. Goldberg, J.H. Holland, Genetic algorithms and machine learning, Mach. Learn. 3 (2) (1988) 95–99. [18] L. Greengard, J. Strain, The fast gauss transform, SIAM J. Sci. Stat. Comput. 12 (1) (1991) 79–94. [19] C. Harris, M. Stephens, A combined corner and edge detector., in: Alvey vision conference, vol. 15, Manchester, UK, 1988, pp. 147–152. [20] W.E. Hart, M.H. Goldbaum, Registering retinal images using automatically selected control point pairs, in: IEEE International Conference on Image Processing, vol. 3, 1994, pp. 576–580. [21] R. Hartley, A. Zisserman, Multiple View Geometry in Computer Vision, Cambridge University, 2003. [22] J.-Z. Huang, T.-N. Tan, L. Ma, Y.-H. Wang, Phase correlation based iris image registration model, J. Comput. Sci. Technol. 20(3) 419–425. [23] L. Ingber, Simulated annealing: Practice versus theory, Math. Comput. Model. 18 (11) (1993) 29–57. [24] B. Jian, B.C. Vemuri, Robust point set registration using gaussian mixture models, IEEE Trans. Pattern Anal. 33 (8) (2011) 1633–1645.

504

G. Wang et al. / Information Sciences 372 (2016) 492–504

[25] M. Kim, G. Wu, Q. Wang, S.-W. Lee, D. Shen, Improved image registration by sparse patch-based deformation estimation, NeuroImage 105 (2015) 257–268. [26] F. Laliberte, L. Gagnon, Y. Sheng, Registration and fusion of retinal images-an evaluation study, IEEE Trans. Med. Imag. 22 (5) (2003) 661–673. [27] P.A. Legg, P.L. Rosin, D. Marshall, J.E. Morgan, Improving accuracy and efficiency of mutual information for multi-modal retinal image registration using adaptive probability density estimation, Comput. Med. Imag. Grap. 37 (7) (2013) 597–606. [28] X. Li, Z. Hu, Rejecting mismatches by correspondence function, Int. J. Comput. Vis. 89 (1) (2010) 1–17. [29] S. Liao, D. Shen, A. Chung, A markov random field groupwise registration framework for face recognition, IEEE Trans. Pattern Anal. 36 (4) (2014) 657–669. [30] D.G. Lowe, Distinctive image features from scale-invariant keypoints, Int. J. Comput. Vis. 60 (2) (2004) 91–110. [31] J. Ma, J. Zhao, Y. Ma, J. Tian, Non-rigid visible and infrared face registration via regularized gaussian fields criterion, Pattern Recognit. 48 (3) (2015) 772–784. [32] J. Ma, J. Zhao, J. Tian, X. Bai, Z. Tu, Regularized vector field learning with sparse approximation for mismatch removal, Pattern Recognit. 46 (12) (2013) 3519–3532. [33] J. Ma, J. Zhao, J. Tian, A.L. Yuille, Z. Tu, Robust point matching via vector field consensus, IEEE Trans. Image Process. 23 (4) (2014) 1706–1721. [34] F. Maes, A. Collignon, D. Vandermeulen, G. Marchal, P. Suetens, Multimodality image registration by maximization of mutual information, IEEE Trans. Med. Imag. 16 (2) (1997) 187–198. [35] G.K. Matsopoulos, P.A. Asvestas, N.A. Mouravliansky, K.K. Delibasis, Multimodal registration of retinal images using self organizing maps, IEEE Trans. Med. Imag. 23 (12) (2004) 1557–1563. [36] C.A. Micchelli, M. Pontil, On learning vector-valued functions, Neural Comput. 17 (1) (2005) 177–204. [37] H.Q. Minh, V. Sindhwani, Vector-valued manifold regularization, in: Proceedings of the 28th International Conference on Machine Learning, 2011, pp. 57–64. [38] A. Myronenko, X. Song, Point set registration: coherent point drift, IEEE Trans. Pattern Anal. 32 (12) (2010) 2262–2275. [39] J.C. Nunes, Y. Bouaoune, E. Delechelle, P. Bunel, A multiscale elastic registration scheme for retinal angiograms, Comput. Vis. Image Understand. 95 (2) (2004) 129–149. [40] J.P. Pluim, J.A. Maintz, M.A. Viergever, Mutual-information-based registration of medical images: a survey, IEEE Trans. Med. Imag. 22 (8) (2003) 986–1004. [41] A. Rangarajan, H. Chui, F.L. Bookstein, The softassign procrustes matching algorithm, in: Information Processing in Medical Imaging, Springer, 1997, pp. 29–42. [42] C. Sanchez-Galeana, C. Bowd, E.Z. Blumenthal, P.A. Gokhale, L.M. Zangwill, R.N. Weinreb, Using optical imaging summary data to detect glaucoma, Ophthalmology 108 (10) (2001) 1812–1818. [43] A. Sotiras, C. Davatzikos, N. Paragios, Deformable medical image registration: a survey, IEEE Trans. Med. Imaging 32 (7) (2013) 1153–1190. [44] C.V. Stewart, C.-L. Tsai, B. Roysam, The dual-bootstrap iterative closest point algorithm with application to retinal image registration, IEEE Trans. Med. Imag. 22 (11) (2003) 1379–1394. [45] A.N. Tikhonov, V.Y. Arsenin, Solutions of Ill-Posed Problems, Winston Washington, DC, 1977. [46] Q.-H. Tran, T.-J. Chin, G. Carneiro, M.S. Brown, D. Suter, In defence of RANSAC for outlier rejection in deformable registration, Springer, pp. 274–287. [47] C.-L. Tsai, C.-Y. Li, G. Yang, K.-S. Lin, The edge-driven dual-bootstrap iterative closest point algorithm for registration of multimodal fluorescein angiogram sequence, IEEE Trans. Med. Imag. 29 (3) (2010) 636–649. [48] G. Wang, Z. Wang, Y. Chen, W. Zhao, A robust non-rigid point set registration method based on asymmetric gaussian representation, Comput. Vis. Image Understand. 141 (2015) 67–80. [49] G. Wang, Z. Wang, Y. Chen, W. Zhao, Robust point matching method for multimodal retinal image registration, Biomed. Signal. Process. 19 (2015) 68–76. [50] G. Wang, Z. Wang, Y. Chen, W. Zhao, X. Liu, Fuzzy correspondences and kernel density estimation for contaminated point set registration, in: 2015 IEEE International Conference on Systems, Man, and Cybernetics, IEEE, 2015c, pp. 1936–1941. [51] G. Wang, Z. Wang, Y. Chen, Q. Zhou, W. Zhao, Context-aware gaussian field for non-rigid point set registration, in: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, IEEE, 2016, pp. 5811–5819. [52] G. Wang, Z. Wang, W. Zhao, Q. Zhou, Robust point matching using mixture of asymmetric gaussians for nonrigid transformation, in: Computer Vision–ACCV 2014, Springer, 2015d, pp. 433–444. [53] G. Wu, M. Kim, Q. Wang, D. Shen, S-hammer: Hierarchical attribute-guided, symmetric diffeomorphic registration for mr brain images, Human Brain Mapp. 35 (3) (2014) 1044–1060. [54] G. Yang, C.V. Stewart, M. Sofka, C.-L. Tsai, Registration of challenging image pairs: Initialization, estimation, and decision, IEEE Trans. Pattern Anal. 29 (11) (2007) 1973–1989. [55] F. Zana, J.-C. Klein, A multimodal registration algorithm of eye fundus images using vessels detection and hough transform, IEEE Trans. Med. Imag. 18 (5) (1999) 419–428. [56] B. Zhang, F. Karray, Q. Li, L. Zhang, Sparse representation classifier for microaneurysm detection and retinal blood vessel extraction, Inf. Sci. 200 (2012) 78–90. [57] J. Zhao, J. Ma, J. Tian, J. Ma, D. Zhang, A robust method for vector field learning with application to mismatch removing, in: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, IEEE, 2011, pp. 2977–2984. [58] L. Zhou, M.S. Rzeszotarski, L.J. Singerman, J.M. Chokreff, The detection and quantification of retinopathy using digital angiograms, IEEE Trans. Med. Imag. 13 (4) (1994) 619–626.