Research Article  Open Access
Retinal Fundus Image Registration via Vascular Structure Graph Matching
Abstract
Motivated by the observation that a retinal fundus image may contain some unique geometric structures within its vascular trees which can be utilized for feature matching, in this paper, we proposed a graphbased registration framework called GMICP to align pairwise retinal images. First, the retinal vessels are automatically detected and represented as vascular structure graphs. A graph matching is then performed to find global correspondences between vascular bifurcations. Finally, a revised ICP algorithm incorporating with quadratic transformation model is used at fine level to register vessel shape models. In order to eliminate the incorrect matches from global correspondence set obtained via graph matching, we proposed a structurebased sample consensus (STRUCTSAC) algorithm. The advantages of our approach are threefold: (1) global optimum solution can be achieved with graph matching; (2) our method is invariant to linear geometric transformations; and (3) heavy local feature descriptors are not required. The effectiveness of our method is demonstrated by the experiments with 48 pairs retinal images collected from clinical patients.
1. Introduction
In this paper, we proposed a novel graphbased registration framework, namely GMICP, to align pairwise retinal images. It consists of three independent steps. First, the retinal vessels are automatically detected and represented as vascular structure graphs; a graph matching (GM) is then performed to find global correspondences between vascular bifurcations; finally, the ICP algorithm is used at fine level to register vessel shape models.
The image registration techniques become increasingly important in clinical human retina diseases diagnosis and treatment. There are large demands for fast and fully automatic algorithms to register two or more retinal images which were taken at different times, from different views, and by different sensors. Accurate registration of intra and intermodality retinal images is required for information integration. Register a sequence of images of the same retina can give a complete view of the retinal fundus region which may lead to early diagnosis and treatment of AIDSCMV retinopathy [1]. In laser retinal surgery, the treatment of exudative AMD may cause highly 50% recurrence rate of CNV (Choroidal Neovascularization) due to inaccurate localization of the laser in fovea area. Realtime image registration method can be used to assist ophthalmologists with avoiding possible damages during the surgery [2].
Retinal images registration is a challenging task [3]. First, the retina is a curved surface; nonlinear deformation may occur when using a weakperspective uncalibrated camera. Second, image overlap may be small due to large viewpoint change between images. Third, retinal images may have large textureless regions with uneven illumination which make the extraction of retinal vessels very difficult. Fourth, two images taken long time apart or from diseased eyes can have physical changes in both structure and color of the retina.
In the past ten years, many techniques have been proposed to solve the retinal image registration problem in different ways. Generally, these methods can be classified into two categories: vesselbased and nonvesselbased methods.
The nonvesselbased methods refer to a class of techniques that do not make use of retinal vessels directly during image registration. These methods can be further subdivide into two classes: intensitybased and image descriptorbased methods. The intensitybased techniques usually rely on image intensities and gradients, the registration is carried out through optimizing a certain similarity function, such as meansquare error, mutual information [4], and crosscorrelation of images. Various of optimization method can be used for finding global optimum of the cost function, including downhill simplex, simulated annealing [4, 5], genetic algorithms [5], and so forth. In [6], a multiscale elastic registration scheme was proposed to take into account retinal intensity variants based on optical flows. The main limitations of intensitybased methods lie in two aspects: (1) these methods are highly reliant on the consistent intensities in two images and tend to fail due to the presence of nonuniform illumination and large textureless regions; (2) the optimization may have huge searching space so that the computational cost becomes a bottleneck for intensitybased methods to be applied in clinic.
The registration methods with image descriptors [7–9] are featurebased and become popular recently. Instead of optimizing a similarity function with whole image intensities, these methods extract local invariant descriptors as image features. And then, image registration is equivalent to finding correspondences between two feature sets. In [8], a generalized dualbootstrap iterative closest point (GDBICP) was proposed, in which SIFT local descriptor is used and the alignment process is driven by two types of keypoints: corner points and face points. In [7], to deal with multimodal registration problem, an edgedriven DBICP algorithm was developed by enriching the keypoint descriptor with shape context using edge points. In [9], a new salient feature region descriptor was proposed for low quality retinal images registration.
The vesselbased methods utilize retinal vascular features as a basis for image matching. Currently, most of the retinal image registration approaches are vesselbased [3, 10–15], because the retinal vessels are believed to be the most appropriate representation for retina. In general, the vesselbased methods consist of two independent steps: vessel segmentation and vascular featurebased registration. Once vessels were extracted, establishing point correspondences between segmented vessels is crucial. In [13], a dualbootstrap iterative closest point (DBICP) algorithm was introduced to match vascular centerlines. It started with some initial loworder estimations and iteratively refined the results by expending the bootstrap regions with higherorder transformation model. In order to make the matching robust to mismatches between feature points, hierarchical transformation models were employed in [3]. In [11], a hybrid retinal image registration approach is proposed by combining both intensitybased and vesselbased methods. To achieve both global and local alignments, an elastic matching scheme is used in [12] based on reconstructed vascular trees. Usually, the vesselbased method are supposed to be more reliable than nonvesselbased method because the vascular features are sufficient (but not too many) and accurate as landmarks and preferable spread over whole image region. The main drawback of vesselbased methods is the local convergence problem when initial misalignment is large and mass of segmentation noises exist.
In this paper, we study retinal image registration following vesselbased approaches. Inspired by the fact that a retinal image may contain some unique geometric structures within its vascular trees. We consider that the correspondence problem is not simply a onetoone mapping problem between point sets but rather a pairwise (structure) matching problem between two separate graphs. A brief overview of our approach is illustrated in Figure 1, in which the retinal feature extraction is preformed at first to generate a graph representation of the vessel networks. Afterward, a graph matchingbased registration is performed to achieve both globally and locally alignment of the retinal images.
To our knowledge, the graphbased approaches are relatively rare in the field of retinal image registration. Related work proposed recently were reported in [16, 17], where a graph transformation matching (GTM) algorithm was developed for retinal image registration and vascular characterization. However, their method is limited in three aspects: (1) a good initial guess of correspondences is crucial to GTM; (2) the number of graph nodes must be equal; and (3) the vessels are detected in a semiautomatic way.
There are three major contributions in this paper. Firstly, we demonstrated that the retinal fundus image registration problem can be efficiently solved within a graphmatching framework. Accordingly, we innovatively represent the retinal vascular features as an undirected graph model and a normalized vessel path distance measure is defined to characterize local similarities between two graph edges. Secondly, we proposed a hybrid registration framework that integrates graph matching and the ICP algorithm together to work with a hierarchy of transformation models. Thirdly, we proposed a structurebased sample consensus (STRUCTSAC) algorithm to eliminate incorrect correspondences from graph matching results.
The remainder of this paper is organized as follows. The feature extraction methods are described in Section 2. The GMICP registration framework is introduced in Section 3. Experiment results are shown in Section 4. Finally, we outline our conclusions in Section 5.
2. Feature Extraction
The key step of GMICP is the extraction of invariant image features. This section presents the feature extraction methods used in our registration framework, including vessel centerline detection and graphbased representation. A detailed flowchart of above two algorithms is illustrated in Figure 2.
(a)
(b)
2.1. Vessel Centerline Detection
Retinal vessels are the most unique and prominent anatomical structure in the retina and provide rich information for the localization of distinct features. A variety of techniques has been introduced for retinal vessel extraction, such as the Hessian method [14], the Gabor filter [18], matched filter [19], and parallel edge method [20]. In our implementation, a singlescale matched filter is employed to detect retinal vessels because it provides better responses to thin and lowcontrast vessels [19]. Although some multiscale techniques [18] can be used to achieve better vessel segmentation results. However, those methods are quite time consuming and not necessary for the image registration task.
Blood vessels in retinal images usually have poor local contrast so that the matched filter is designed to detect and enhance the local piecewise vessel segments through image convolution. To approximate the vessel profile, Gaussian function is commonly adopted as the kernel of matched filter [19]. However, in some cases, the intensity profile of some wide vessel segments is not Gaussian due to vessel central light reflection effect. Hence, we use a GaussianHermite mixture model described in [21] to approximate local vessel segment structures. The 1D secondorder GaussianHermite kernel is defined as where is a weight parameter, that when , the kernel is Gaussian. In order to detect a vessel oriented in any directions, a set rotated kernels is applied to a fundus image. Because of symmetry, only possible directions are required for computation. At each orientation, the twodimensional kernel can be represented as two onedimensional kernels applied in succession according to Cartesian separability rule. Let be the image intensity at position , a matchedfilter response (MFR) of the input image at orientation is defined as where is GaussianHermite kernel in normal direction of , is secondorder derivative of the GaussianHermite kernel along the direction of . To make the implementation easier, the image coordinate is rotated at each orientation instead of rotating convolution kernels. Finally, the matched filer is implemented by maximizing the responses over all orientations at each pixel, which can be expressed as
In our implementation, 12 matched filter orientations are used with a 15degree increment. In order to avoid the influence of image noise, a histogrammatch image filtering and a median image filtering was applied before MF vessel detection. After the main vascular structures were detected, we used a threshold probing technique proposed in [22] to distinguish between enhanced vessels and the image background. Finally, vessel thinning is performed to the binary vessel image so that the resulting patterns consist of lines and curves with one pixel wide only.
2.2. GraphBased Representation
Once the retinal vessel centerline (vessel shape) model is calculated, the transformation into a graph is straightforward. We construct a undirected vascular structure graph (VSG) , where each graph node represents a vascular bifurcation point or a vessel segment end point; each graph edge corresponds to a vascular segment between two vascular feature points. The attributes set carries some scalar numbers or vectors to characterize the local vessel segments. In our implementation, two distance measures were assigned to each graph edge: vessel path distance and Euclidean distance between two linked graph nodes. The vessel path distance measures the vessel path length between two vascular feature points in a binary vessel centerline image. To make the proposed method invariant to scale, we normalize in the following way: likewise for . In VSG, all the graph nodes are treated uniformly, we do not distinguish between the type of vertices on anatomical aspects.
The vascular structure graph model is built by following two steps. First, the bifurcations of vessels are extracted using a template matching technique described in [10]. Then, we performed a parallel breadth first search (BFS) on the monopixel binary retinal vessel network to build edges among those bifurcation points and vessel segment end points. The vessel path distance is calculated at the same time.
3. The GMICP Algorithm
In this section, we will introduce the proposed GMICP registration framework in details. As mentioned above, we solve retinal image registration problem through a featurebased approach. Hierarchical features were utilized in our method, including retinal vessel shape models, vascular bifurcations and the underlying vascular topological structures, which were extracted previously by using the methods described in section 2. The key problem to register two retinal images in different coordinate positions is to establish correspondences between two feature sets. Then, a closedform solution [23] can be applied to estimate the transformation parameters directly.
The simplest approach to the onetoone correspondence problem is to define some similarity measure between two feature points, for example, the Euclidean distance between SIFT descriptors. And then the matching is carried out in a highdimensional feature space following the nearest neighbor rule. This approach will fail when the ambiguities exist, such as nondiscriminative local appearance or repeated patterns. However, these situations are quite common to retinal images. Another way to match two sets of points is to employ some heuristic techniques, such as ICP [23] or EM [24] algorithm, to iteratively align two feature sets until it converges. These methods could be robust if the initial misalignment is small and the noise ratio is not high. However, the features extracted by using current automatic retinal image segmentation techniques will always contain a mass of noises and outliers. The transformation between two images could also be relatively large. As a result, these methods will often converge to local incorrect positions.
Our method was motivated by the observation that a retinal image may contain some unique geometric structures within its vascular trees. Indeed, we are trying to make use of those distinctive structures and to enforce some geometric consistency between pairs of feature correspondences. Furthermore, we find that the graphmatching approach is a natural and powerful way for doing this task. Thus, we want to obtain points correspondence by taking advantage of structure information and the edgetoedge comparison mechanism with graph matching. More local invariant and expressive features with relatively low computational costs can also benefit from this.
The flowchart of the proposed registration algorithm is given in Figure 3. In our implementation, the graduated assignment algorithm is employed to match vascular structure graphs extracted from two retinal images. In order to eliminate wrong matches from global correspondences set obtained via graph matching, a structurebased sample consensus (STRUCTSAC) algorithm is proposed. This algorithm is more efficient and reliable than RANSAC [25] by replacing the random sampling procedure with a sequential local structure sampling. Once the global transform is known, we use a kdtreebased ICP algorithm to register two vessel shape models to achieve better registration precision.
We use similarity transformation model for global correspondence calculation under the assumption which the apparent motion of the retina is generally rigid [3]. To correct some nonlinear motions in the boundary area with low overlapping images, a quadratic polynomial model is adopted for vessel shape model registration.
3.1. Problem Formulations
We consider two vessel structure graphs and , which were extracted from reference image and floating image , respectively. Let , . We do not assume that because there may be different numbers of features to be matched. The graphmatching problem is to find an optimal assignment between and that maximizes the sum of pairwise affinities between edges and . This is equivalent to seek a assignment matrix to maximize the objective function where is equal to 1 if node in corresponds to node in . is a compatibility matrix. carries potential matching weights corresponding to the pairs of points of and of . A rectangle rule may better illustrate this edgematching mechanism, which is shown in Figure 4.
The problem (5) can also be reformulated as an Integer Quadratic Programming (IQP) problem. Let us consider as a binary vector . Then, (5) may take the form as In order to achieve onetoone correspondences, an affine constraint is enforced to assignment matrix as well as to assignment vector
3.2. VSG Matching via Balanced Graduated Assignment
Due to its combinatorial nature, to solve the graph matching problem in an exact way could be NPHard and it is not applicable to our VSG matching problem. Therefore, we are seeking to find some approximate solutions. There are many approaches to inexact graph matching in the literature, such as the graduated assignment (GA) algorithm [26], spectral relaxation (SM) method [27], graph edit distancebased method [28], and shock graph grammarbased method [29]. Among these methods, the continuous relaxations approaches were thought to be the most successful ones [27].
In this paper, we employed a graduated assignment algorithm originally proposed by Gold and Rangarajan in [26] as our graphmatching framework. Our empirical evaluation shows that the method is efficient and reliable with VSG matching problem and even when outliers are present. The GA method relaxes the discrete IQP into a continuous nonconvex quadratic program (QP). It optimizes a quadratic cost function through minimizing its loworder Taylor expansion around previous solution. The deterministic annealing scheme was employed to find globally optimum solutions. Onetoone correspondence for graph nodes can be guaranteed via Sinkhorn twoway constraint procedure. Thus, (5) can be reformulated as in following expressions under the GA framework where is a doubly stochastic matrix, in which the maximum element in each row and column represents a correspondence between two graphs nodes. The term in (8) is a smoothing function which can push the minimum of the objective away from the local discrete points. The weighting parameter is introduced to avoid taking the logarithm of zero. The parameter controls the speed of annealing. We define the measure of compatibility between the links of two graphs as follows: where L is the maximum edge weights in both graphs. Two edge weights are involved for pairwise affinity calculation, where and as defined previously. is a weight parameter to balance the influence of two edge distance measures. Empirically, we set
From our experiments, we found that the ambiguities in the compatibility matrix may bring down the quality of graph matching in finding global correspondences. For example, one edge may have lots of potential matches in another graph. Such an edge is not discriminative. On the other hand, an edge with small group of potential matches may help to clarify the optimal matching. To this end, we employ a bistochastic normalization technique described in [27] to balance the compatibility matrix (To enhance correct correspondences and to suppress spurious correspondences). The normalization is carried out in three steps. First, the compatibility matrix W is represented as a edge similarity matrix S, where and . Then, the edge similarity matrix S is normalized iteratively by using following expressions until convergence After the normalization was finished, the edge similarity matrix S is converted back to the compatibility matrix W
3.3. STRUCTSAC
The global correspondences obtained by graph matching may contain a certain number of incorrect matches. In some extreme cases, the ratio of wrong matches can be greater than 80%. To eliminate the wrong matches from global correspondences set is crucial to our method.
The RANSAC algorithm [25] is a wellknown method to deal with outliers and has become a standard in the field of computer vision. It first fits model parameters with randomly selected subsets of input data and then evaluates the quality of the parameters on the whole input data. The process is performed repeatedly and terminated when the probability of finding a better model becomes lower than a userspecified threshold. The main drawback of RANSAC is that most of the verified model parameters are incorrect under arbitrary random sampling. This is inefficient and time consuming when dealing with large data sets.
In this paper, we proposed a local structurebased sample consensus (STRUCTSAC) algorithm, based on two simple observations: (i)correct matches between graph nodes are not isolated; (ii)most of the correctly matched graph edges are found next to each other.
The phenomenon described above can be interpreted by the property of graduated assignment that it trends to approximately find the maximum common subgraph (MCS) during graph matching. As a result, a good matched edge may have some good matched neighbors to support it. On the other hand, a noise or mismatched edge may influence its neighbors to be mismatched. Furthermore, we thought that the type of correctly matched structures in graphmatching results can be divided into two categories: stable structures and unstable structures (as illustrated in Figure 5). A graph node is stable if it has at least two correctly matched neighbors. An edge between two stable nodes is a stable edge. Otherwise, the structures are unstable. Thus, the goal of our algorithm is simplified to find some stable graph structures.
The STRUCTSAC algorithm was motivated by the hypothesizeandverify idea from RANSAC. However, instead of random sampling, we select subsets of input data in a deterministic way by utilizing the local structures of a sample point. A brief description of STRUCTSAC is as follows. Repeatedly, select a graph node and with its neighbors as the source subset. Find their corresponding (the correspondence was obtained via graph matching) nodes in other graph as the object subset. (As an option, the RANSAC strategy can be used here to further distinguish the subsets from noisy structures when the number of neighbors is large.) Compute the transformation parameters determined by the selected couples and verify it first with a local geometric consistency test using partial data and then with global model consistency test using whole input data points. The process is terminated when all the potential graph nodes were traversed and the best model parameter with lowest quality score is returned. In our implementation, the search is performed symmetrically and only on the nodes with 2 or 3 degrees. Since these nodes are stable in VSG and provide sufficient information for the estimation of model parameters. The effect of STRACSAC is demonstrated in Figure 6.
(a)
(b)
In local geometric consistency test, we utilize pairwise edge similarity and triplewise orientation similarity of graph nodes to distinguish local structure samples. In global model consistency test, we use sum of squared distance to measure the quality of model parameter between two point sets, as below: where T denotes the similarity transformation function, is a temporary correspondence set determined by transformation parameter . is a Huber kernel [30] to limit the influence of outliers
3.4. Vessel Shape Model Registration Based on HigherOrder Transformation Model
Since the retina is spherical and the retinal images captured with weakperspective cameras may have nonrigid distortion in localized regions, higherorder transformation model is needed for accurate registration. This is more important to the registration tasks with lowoverlap retinal images. Figure 7(a) shows two fused retinal images under linear registration where some misalignment vessels were indicated with arrows. In Figure 7(b), these vascular segments are well matched after a quadratic registration.
(a)
(b)
(c)
(d)
Our registration algorithms employ a hierarchy of transformation models [3], which is also known as the global to local strategy, to account for the nonlinearities described above. One advantage of hierarchical transformation is that it makes the algorithm robust to the mismatches caused by large interframe motions [3]. In our method, the global registration is achieved via graph matching and STRUCTSAC with similarity transformation model. At fine (local) level, we use a revised ICP algorithm to register vessel shape (centerline) models. A 12dimensional secondorder polynomial transformation model is employed in the local registration, that it maps a source point in floating vessel model onto the target point in reference vessel model as follows: where are transformation parameters. Once the correspondence is obtained at each iteration of ICP, the quadratic motion can be estimated in closedform by using a linear regression method described in [24].
In addition, several improvements are made to the standard ICP algorithm in order to achieve better performance.
(1) To accelerate the closest point indexing, a kdtree data structure is applied for underlying point sets modeling; (2) The registration is performed in a multiresolution scheme; (3) A point pair weighting strategy [31] is used to reject inconsistent matches.
4. Experiments
In this section, we will evaluate the performance of the proposed GMICP registration algorithm in two aspects: (1) saliency of the vascular structure graph, and (2) the overall performance of GMICP in retinal image registration. A total of 48 pairs of retinal images collected from clinical patients were involved for the algorithm evaluation. In order to make the evaluation consistent with images in different size, we resampled all test images to approximately using bicubic image interpolation.
Our retinal image registration program was designed with matlab and some core algorithms (such as graph matching and vascular bifurcations extraction) were implemented in CMEX. All the results reported in this paper were obtained on a PC with Intel Core 2 Duo 2.33 GHz CPU and 3 G RAM running under Windows 7 Professional.
4.1. Saliency of Vascular Structure Graph
The aim of our first experiment is to show the saliency and performance benefits of vascular structure graphs presented in this paper. Based on the vascular bifurcation points extracted with our methods, several spatial topological graph structures that commonly used in computer vision and graphmatching research were investigated in this study, including Delaunay triangulation graph (DT), minimum spanning tree of DT graph (DTMST), knearest neighbor graph (KNN) [17], and minimum spanning tree of KNN graph (KNNMST). Some examples for these graph structures are illustrated in Figure 8. The major difference between VSG and other graph representations is the anatomical properties of retinal vessels that VSG characterizes while the other graph representations do not provide.
(a)
(b)
(c)
(d)
The graph matching results with five different graph structures are shown in Table 1. Three evaluation criterions are presented in our results: the average number of correctly matched graph nodes, the recall rate and the success rate of registration. The recall value is defined as the number of correct matches with respect to the number of correspondences [32]. Our experiment results denotes that the VSG representation for retinal images substantially outperforms other topological graph structures in graph matching. In addition, we also plot recall rate against the percentage of structure (graph node) noise with gaussian curve fit of data in Figure 9. As can be seen, the VSG has obtained relatively high matching rate even with greater than 80% of outliers.

4.2. Overall Performance of GMICP Algorithm
Evaluating the accuracy of the retinal image registration is not an easy task because of lack of ground truth. In this section, we evaluate our algorithm with three assessment criteria:(i)normalized Correlation Coefficient (NCC), (ii)normalized Mutual Information (NMI), (iii)vessel Centerline Error Measure (CEM).
The first two terms measure appearance difference between reference and floating retinal images based on intensity statistics. The third term is a standard evaluation criteria that commonly adopted in the field of retinal image analysis [3, 13]. The CEM is defined as median distances over all corresponding points between two vessel centerline models, expressed as where and are two matched vessel points determined by correspondence set . denotes the final transformation function. In addition, the success rate and run time performance was investigated in our experiments as well. In our settings, a registration was considered successful if the CEM error is smaller than 3.0 pixels. A manual validation was also performed in order to ensure the correctness of our results.
In the second experiment, we compare the performance of GMICP with other two registration algorithms: the standard ICP and GDBICP. The standard ICP was implemented in Matlab and we try to use it to register twovessel shape models without initial alignment. The GDBICP [8] is a generalized version of DBICP which was originally proposed for retinal image registration. We had obtained a public copy of its binary program from the Internet and tested it under the quadratic transformation parameters setting. In order to assess the reliability of the graph matching algorithm, we also performed some image registration tests using the graph matching (GM) alone without ICP.
The results of comparative experiments on 48 retinal image pairs are listed in Table 2. As can be seen, the registration accuracy of GMICP (0.81 pixels in CEM) is very close to GDBICP (0.79 pixels in CEM) according to three error assessment criteria. Both of the two algorithms have achieved subpixel accuracy. In absence of fine level alignment, the registration with only graphmatching is still robust and obtained acceptable results (1.00 pixels in CEM). The GMICP successfully registered 95.83% of all pairs of test retinal images, two pairs failed due to poor quality of vessel segmentation from pathology images whereas GDBICP registered 97.92% and only one pair failed. It takes about 8 seconds on average to register a pair of images with GMICP, including the time cost in both feature extraction and registration, while GDBICP needs more than 20 seconds. In contrast with above three methods, the performance of standard ICP is not satisfactory in both registration accuracy and success rate. Besides, it is a little surprise to see that the NCC and NMI assessment of standard ICP is even better than the other three methods. It denotes that the intensitybased global similarity measures for retinal images are nonconvex and sometimes become unreliable as a result of the presence of largescale misalignment and inconsistent patterns.

In most of the cases, the performance of registration will be gradually degraded with less overlap of pair of images. We plot the vessel centerline error measure against the percentage of image overlap, as shown in Figure 10. It is observed that GMICP is accurate and robust with medium overlapped image pairs. Due to the fact that the graph matching scheme does not rely on direct distance metric between two feature sets, the GMICP registration framework is actually invariant to any linear geometric transformations, such as largescale initial translations and rotations. However, our method may fail with some low overlapped (under 40%) image pairs due to lack of enough consistent vascular segments for matching.
One difficulty to all vesselbased registration methods is that the reliability of registration is highly influenced by the vessel segmentation results. In our method, there are serval ways to deal with segmentation noise and structure inconsistency. First, small edges were removed during vascular graph representation according to some vessel network topological modification rules [15]. Second, a structurebased sample consensus algorithm was proposed to eliminate incorrect matches in global registration. Third, a point pair weighting strategy is used to reject inconsistent matches during fine level registration. The performance of GMICP for the images with high structure noise and pathologies is visually demonstrated in Figure 11.
5. Conclusion
In this paper, we proposed a graphbased algorithm framework for retinal image registration. Hierarchical retinal features were utilized in the framework, including vessel shape models, vascular bifurcations and the underlying vascular topological structures. The core idea behind our method is to obtain the point correspondences by taking into account structure information provided by human retina and retinal vessels. Thus, the image registration problem is transformed into a pairwise (edgetoedge) correspondence problem. Accordingly, we have employed a robust graph matching framework, the graduated assignment with normalization on compatibility matrix, to match two retinal vascular structure graphs. In order to make our method robust to structure noise, the STRUCTSAC algorithm is proposed to eliminate incorrect matches from graph matching results. The effectiveness of our method has been demonstrated experimentally with 48 pairs of real retinal fundus images.
The VSG representation for retinal vessel network is distinct, preferable spread and low computational cost. Additionally, it reflects the anatomical characterizations contained within a retinal image, which is prominent and almost unchangeable. Our experiments shows that the VSG model is robust even with high structure noise and substantially outperforms other regular topological structures for graph matching.
The most significant distinguishing property of our method is the structure matching and verification scheme for retinal image registration. Although graph matching have proven to be useful in many applications, such as character recognition, shape analysis [29], however, it is relatively rare for medical image registration. The reasons are twofold: (1) structure representation from image is not straightforward; (2) structure matching is not robust under high noise level. In our method, sufficient consistent structures are extracted from retinal images with proper selection of vessel segmentation and graph representation methods. The adopted graph matching framework incorporating with STRUCTSAC works fine with retinal vascular structure graphs. No initial guess is needed. Heavy local feature descriptors are not required (actually, only two simple distance measures are involved in the computation). Furthermore, it gives a convenient and practical way to distinguish between correct and incorrect structure correspondences.
The main limitation of the proposed method is the assumption that there are enough consistent vascular segments (edges) for matching. However, with some low overlapped or poor quality pathology image pairs, this assumption is not satisfied. As a result, the performance becomes poor (the number of stable graph nodes is less than 3) and even fails.
In this paper, we only addressed pairwise structure matching problem with retinal image registration. It would be interesting to cope with some higherorder graph matching methods [33] in our framework. As a further development, the registration accuracy can be improved by using more efficient and reliable vessel segmentation techniques. Also, our future work will focus on developing more robust graph representation for retinal features, which can accommodate both geometric and photometric invariants in an image.
Acknowledgments
This work is supported by the Program of the National Basic Research and Development Program of China (973) under Grant no. 2006CB705700, the Cheung Kong Scholars and Innovative Research Team in University (PCSIRT) under Grant no. IRT0645, the Chair Professors of Cheung Kong Scholars Program of Ministry of Education of China, CAS Hundred Talents Program, the National Natural Science Foundation of China under Grant no. 30873462, 30900334,60771068, the CAS Scientific Research Equipment Develop Program (YZ0642, YZ200766), and the Natural Science Basic Research Plan in Shaanxi Province of China under Grant no. 2009JQ8018.
References
 M. G. DeGrezia and M. Robinson, “Ophthalmic manifestations of HIV: an update,” The Journal of the Association of Nurses in AIDS Care, vol. 12, no. 3, pp. 22–32, 2001. View at: Google Scholar
 S. L. Fine, “Observations following laser treatment for choroidal neovascularization,” Archives of Ophthalmology, vol. 106, no. 11, pp. 1524–1525, 1988. View at: Google Scholar
 A. Can, C. V. Stewart, B. Roysam, and H. L. Tanenbaum, “A featurebased, robust, hierarchical algorithm for registering pairs of images of the curved human retina,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 24, no. 3, pp. 347–364, 2002. View at: Publisher Site  Google Scholar
 N. Ritter, R. Owens, J. Cooper, R. H. Eikelboom, and P. P. Van Saarloos, “Registration of stereo and temporal images of the retina,” IEEE Transactions on Medical Imaging, vol. 18, no. 5, pp. 404–418, 1999. View at: Publisher Site  Google Scholar
 G. K. Matsopoulos, N. A. Mouravliansky, K. K. Delibasis, and K. S. Nikita, “Automatic retinal image registration scheme using global optimization techniques,” IEEE Transactions on Information Technology in Biomedicine, vol. 3, no. 1, pp. 47–60, 1999. View at: Google Scholar
 J. C. Nunes, Y. Bouaoune, E. Delechelle, and PH. Bunel, “A multiscale elastic registration scheme for retinal angiograms,” Computer Vision and Image Understanding, vol. 95, no. 2, pp. 129–149, 2004. View at: Publisher Site  Google Scholar
 C.L. Tsai, C.Y. Li, G. Yang, and K.S. Lin, “The edgedriven dualbootstrap iterative closest point algorithm for registration of multimodal fluorescein angiogram sequence,” IEEE Transactions on Medical Imaging, vol. 29, no. 3, pp. 636–649, 2010. View at: Publisher Site  Google Scholar
 G. Yang, C. V. Stewart, M. Sofka, and C.L. Tsai, “Registration of challenging image pairs: initialization, estimation, and decision,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 29, no. 11, pp. 1973–1989, 2007. View at: Publisher Site  Google Scholar
 J. Zheng, J. Tian, Y. Dai, K. Deng, and J. Chen, “Retinal image registration based on salient feature regions,” in Proceedings of the 31st Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC '09), pp. 102–105, 2009. View at: Publisher Site  Google Scholar
 F. Zana and J. C. Klein, “A multimodal registration algorithm of eye fundus images using vessels detection and hough transform,” IEEE Transactions on Medical Imaging, vol. 18, no. 5, pp. 419–428, 1999. View at: Google Scholar
 T. Chanwimaluang, G. Fan, and S. R. Fransen, “Hybrid retinal image registration,” IEEE Transactions on Information Technology in Biomedicine, vol. 10, no. 1, pp. 129–142, 2006. View at: Publisher Site  Google Scholar
 B. Fang and Y. Y. Tang, “Elastic registration for retinal images based on reconstructed vascular trees,” IEEE Transactions on Biomedical Engineering, vol. 53, no. 6, Article ID 1634512, pp. 1183–1187, 2006. View at: Publisher Site  Google Scholar
 C. V. Stewart, C.L. Tsai, and B. Roysam, “The dualbootstrap iterative closest point algorithm with application to retinal image registration,” IEEE Transactions on Medical Imaging, vol. 22, no. 11, pp. 1379–1394, 2003. View at: Publisher Site  Google Scholar
 G. K. Matsopoulos, P. A. Asvestas, N. A. Mouravliansky, and K. K. Delibasis, “Multimodal registration of retinal images using self organizing maps,” IEEE Transactions on Medical Imaging, vol. 23, no. 12, pp. 1557–1562, 2004. View at: Publisher Site  Google Scholar
 F. Laliberté, L. Gagnon, and Y. Sheng, “Registration and fusion of retinal images—an evaluation study,” IEEE Transactions on Medical Imaging, vol. 22, no. 5, pp. 661–673, 2003. View at: Publisher Site  Google Scholar
 W. Aguilar, M.E. MartinezPerez, Y. Frauel, F. Escolano, M. A. Lozano, and A. EspinosaRomero, “Graphbased methods for retinal mosaicing and vascular characterization,” in Proceedings of the 6th IAPRTC15 International Workshop on GraphBased Representations in Pattern Recognition, vol. 4538 of Lecture Notes in Computer Science, pp. 25–36, 2007. View at: Google Scholar
 W. Aguilar, Y. Frauel, F. Escolano, M. E. MartinezPerez, A. EspinosaRomero, and M. A. Lozano, “A robust Graph Transformation Matching for nonrigid registration,” Image and Vision Computing, vol. 27, no. 7, pp. 897–910, 2009. View at: Publisher Site  Google Scholar
 J. V. B. Soares, J. J. G. Leandro, R. M. Cesar Jr., H. F. Jelinek, and M. J. Cree, “Retinal vessel segmentation using the 2D Gabor wavelet and supervised classification,” IEEE Transactions on Medical Imaging, vol. 25, no. 9, Article ID 1677727, pp. 1214–1222, 2006. View at: Publisher Site  Google Scholar
 M. Sofka and C. V. Stewart, “Retinal vessel centerline extraction using multiscale matched filters, confidence and edge measures,” IEEE Transactions on Medical Imaging, vol. 25, no. 12, pp. 1531–1546, 2006. View at: Publisher Site  Google Scholar
 A. H. Can, H. Shen, J. N. Turner, H. L. Tanenbaum, and B. Roysam, “Rapid automated tracing and feature extraction from retinal fundus images using direct exploratory algorithms,” IEEE Transactions on Information Technology in Biomedicine, vol. 3, no. 2, pp. 125–138, 1999. View at: Publisher Site  Google Scholar
 L. Wang, A. Bhalerao, and R. Wilson, “Analysis of retinal vasculature using a multiresolution hermite model,” IEEE Transactions on Medical Imaging, vol. 26, no. 2, pp. 137–152, 2007. View at: Publisher Site  Google Scholar
 A. Hoover, “Locating blood vessels in retinal images by piecewise threshold probing of a matched filter response,” IEEE Transactions on Medical Imaging, vol. 19, no. 3, pp. 203–210, 2000. View at: Publisher Site  Google Scholar
 P. J. Besl and N. D. McKay, “A method for registration of 3D shapes,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 14, no. 2, pp. 239–256, 1992. View at: Publisher Site  Google Scholar
 N. Ryan, C. Heneghan, and P. De Chazal, “Registration of digital retinal images using landmark correspondence by expectation maximization,” Image and Vision Computing, vol. 22, no. 11, pp. 883–898, 2004. View at: Publisher Site  Google Scholar
 Martin A. Fischler and Robert C. Bolles, “Random sample consensus: a paradigm for model fitting with applications to image analysis and automated cartography,” Communications of the ACM, vol. 24, no. 6, pp. 381–395, 1981. View at: Publisher Site  Google Scholar
 S. Gold and A. Rangarajan, “A graduated assignment algorithm for graph matching,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 18, no. 4, pp. 377–388, 1996. View at: Google Scholar
 T. Cour, P. Srinivasan, and J. Shi, “Balanced graph matching,” in Advanced in Neural Information Processing Systems, 2006. View at: Google Scholar
 R. Myers, R. C. Wilson, and E. R. Hancock, “Bayesian graph edit distance,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 22, no. 6, pp. 628–635, 2000. View at: Publisher Site  Google Scholar
 K. Siddiqi, A. Shokoufandeh, S. J. Dickinson, and S. W. Zucker, “Shock graphs and shape matching,” International Journal of Computer Vision, vol. 35, no. 1, pp. 13–32, 1999. View at: Publisher Site  Google Scholar
 G. D. Hager and P. N. Belhumeur, “Efficient region tracking with parametric models of geometry and illumination,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 20, no. 10, pp. 1025–1039, 1998. View at: Google Scholar
 S. Rusinkiewicz and M. Levoy, “Efficient variants of the ICP algorithm,” in Proceedings of the 3rd International Conference on 3D Digital Imaging and Modeling, pp. 145–152, 2001. View at: Google Scholar
 K. Mikolajczyk and C. Schmid, “A performance evaluation of local descriptors,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 27, no. 10, pp. 1615–1630, 2005. View at: Publisher Site  Google Scholar
 R. Zass and A. Shashua, “Probabilistic graph and hypergraph matching,” in Proceedings of the 26th IEEE Conference on Computer Vision and Pattern Recognition (CVPR '08), June 2008. View at: Publisher Site  Google Scholar
Copyright
Copyright © 2010 Kexin Deng et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.