Registration of retinal images from Public Health by minimising an error between vessels using an affine model with radial distortions

TL;DR

This paper presents a novel retinal image registration method that models affine transformations with radial distortions, using vessel overlap errors, achieving high accuracy on public datasets and outperforming existing methods.

Contribution

The paper introduces a new registration approach combining affine and radial distortion models, with vessel-based error minimization, improving accuracy over previous methods.

Findings

Successfully registers 96% of image pairs in a public dataset

Outperforms previous methods and state-of-the-art techniques

Achieves better registration results than reference methods on public datasets

Abstract

In order to estimate a registration model of eye fundus images made of an affinity and two radial distortions, we introduce an estimation criterion based on an error between the vessels. In [1], we estimated this model by minimising the error between characteristics points. In this paper, the detected vessels are selected using the circle and ellipse equations of the overlap area boundaries deduced from our model. Our method successfully registers 96 % of the 271 pairs in a Public Health dataset acquired mostly with different cameras. This is better than our previous method [1] and better than three other state-of-the-art methods. On a publicly available dataset, ours still better register the images than the reference method.

Tables

Table 1. Table 1 : Decreasing percentage of successful superimpositions for five methods in a PH dataset.

	Method	% correct
(i)	Current model (vessels)	96 %
(ii)	Current model (SIFT-points) [1]	92 %
(iii)	Lee et al. [10]	88 %
(v)	“REMPE” [3]	75 %
(iv)	“gdbicp” quadratic [11]	74 %

Table 2. Table 2 : AUC of the current method and “REMPE” for the categories 𝒮 𝒮 \mathcal{S} , 𝒫 𝒫 \mathcal{P} and 𝒜 𝒜 \mathcal{A} and the whole FIRE dataset.

	Method	$𝒮$	$𝒫$	$𝒜$	FIRE
(i)	Current method	0.942	0.632	0.768	0.810
(v)	“REMPE”	0.935	0.511	0.599	0.745

Table 3. Table 3 : Mean (and standard deviation) error of the current method and “REMPE” for the categories 𝒮 𝒮 \mathcal{S} , 𝒫 𝒫 \mathcal{P} and 𝒜 𝒜 \mathcal{A} and the whole FIRE dataset. Units are in pixels.

Method	$𝒮$	$𝒫$	$𝒜$	FIRE
Current	1.46 (1.12)	9.25 (10.00)	5.81 (7.21)	4.76 (7.47)
“REMPE”	1.63 (1.57)	12.64 (15.19)	14.05 (25.73)	6.96 (13.65)

Equations

$${\mathrm{u}}_{k_2}({\mathrm{p}}_2)=\mathrm{H}[{\mathrm{u}}_{k_1}({\mathrm{p}}_1)],$$

$$d^2(\mathrm{s},\mathcal{R})\approx \frac{d}{d-\rho }[(\mathrm{s}-{\mathrm{r}}^{\ast }).\mathrm{t}({\mathrm{r}}^{\ast }){]}^2+[(\mathrm{s}-{\mathrm{r}}^{\ast }).\mathrm{n}({\mathrm{r}}^{\ast }){]}^2$$

$${\mathrm{C}}_i=\left[\begin{array}{ccc}1& 0& -x_i\\ 0& 1& -y_i\\ -x_i& -y_i& -r_i^2+x_i^2+y_i^2\end{array}\right]$$

$${\mathrm{x}}^T{\mathrm{C}}_1^u\mathrm{x}={\mathrm{x}}^{\prime T}({\mathrm{H}}^{-1}{)}^T{\mathrm{C}}_1^u({\mathrm{H}}^{-1}){\mathrm{x}}^{\prime }={\mathrm{x}}^{\prime T}{\mathrm{C}}_1^r{\mathrm{x}}^{\prime },$$

$$\{\begin{array}{ccc}{\mathrm{x}}^{\prime T}{\mathrm{C}}_1^r{\mathrm{x}}^{\prime }& \le & 0\\ {\mathrm{x}}^{\prime T}{\mathrm{C}}_2^u{\mathrm{x}}^{\prime }& \le & 0\end{array}.$$

Full text

**Index Terms— ** eye fundus images, image registration, public health, radial distortion, vessel error

1 Introduction

The existence of diabetic retinopathy (DR) screening programmes has led to the creation of large Public Health (PH) databases of colour eye fundus images which allow to perform longitudinal (i.e. temporal) analysis. This analysis is facilitated by a perfect superimposition of the images. However, as images are often captured with different cameras with at least a year of interval, an appropriate method is necessary to correct [1]: (i) the different positions of the patient (rotation, translation, scaling), (ii) the change of the camera (scaling and radial distortion), (iii) the radial distortion caused by the projection of the retina (a spherical cap) onto the sensor plane, (iv) the radial distortion due to the camera optics and (v) the contrast changes between the images. For such a reason, we introduced in [1] a two-step method which consists of a pre-processing to correct the contrast variations and a registration model composed of an affinity and two radial distortion corrections. The model parameters are estimated with characteristic points extracted by the scale-invariant feature transform (SIFT) [2]. This estimation is generally sufficient in many images. However some may present noticeable differences on their external part, especially when the overlap area is small (less than 50 %). The aim of this paper is to address this issue by using the vessels to estimate the model in addition to the SIFT points. We will provide closed-form equations of the overlap area to efficiently select the vessel parts in this area. The paper is organised in two parts. Firstly, we will present our improved method. Secondly, we will compare it to a recent one “REMPE” [3]. We will use a PH dataset with 271 image pairs acquired mostly with different cameras [1]. We will recall the results we obtained in [1] in this dataset for three state-of-the-art methods. A second comparison will be performed in the publicly available dataset “FIRE” [4] associated with “REMPE” [3].

2 Method

A superimposition method requires a model of deformation and an error criterion to estimate its parameters. Let us remind the model and its estimation which were both presented in [1]. We will then present the new error criterion and its efficient computation by selecting the vessels in the overlap area.

2.1 The model and its estimation

Let $\mathrm{p}_{1},\mathrm{p}_{2}\in\mathbb{R}^{2}$ be two corresponding points in the initial images 1 and 2. Our model is based on an affine homography $\mathrm{H}$ and two radial distortion corrections $\mathrm{u}_{k_{1}}$ and $\mathrm{u}_{k_{2}}$:

[TABLE]

where $\mathrm{H}=\begin{bmatrix}\mathrm{A}&\mathrm{t}\\ O^{T}&1\end{bmatrix}$. $\mathrm{A}$ is a non-singular matrix of $\mathbb{R}^{2}$ representing the linear applications (rotation, translation, scaling, etc.). The vector $\mathrm{t}=[t_{x},t_{y}]^{T}\in\mathbb{R}^{2}$ is a translation. The radial distortion correction $\mathrm{u}_{k}$ is defined by the point $\mathrm{u}_{k}(\mathrm{p})=[\mathrm{p}^{u},1]^{T}\in\mathbb{R}^{3}$, in homogeneous coordinates [5]. The undistorted point $\mathrm{p}^{u}\in\mathbb{R}^{2}$ is given by $\mathrm{p}^{u}-\mathrm{c}=\overline{\mathrm{p}}/(1+k\left\|\overline{\mathrm{p}}\right\|^{2})$. $\overline{\mathrm{p}}=\mathrm{p}-\mathrm{c}$ is a point whose coordinates are centred on the image centre $\mathrm{c}$. If the same camera is used to capture both images, the real distortion parameters are equal to $k_{1}=k_{2}=k$.

The model (Eq. 1) is estimated by a several-stage approach (Fig. 1). (1) As many images in PH databases present a non-uniform brightness, a preprocessing corrects the colour contrast variations. (2) Characteristic points are extracted using SIFT algorithm [2] and matched between the images following the method presented in [1]. (3) The matched points serve to initialise the model and the number of distortions - 1 or 2 - is automatically selected. (4) An iterative estimation is performed on the parameters until the convergence of the error. Linear estimators are used to initialise the non-linear optimisers [1]. (5) A non-linear optimiser refines the model estimate [1]. In this paper, at the stages (4) and (5), we will replace the SIFT-point error we used in [1] by a criterion based on an error between the vessels.

2.2 An error criterion between the vessels

SIFT points are used for the initialisation stage (3) (Fig. 1). However, during the stages (4) and (5), an error between the closest vessels is minimised to ensure a better superimposition on the whole overlap area. In the plane $\mathbb{R}^{2}$, the error is measured between the closest vessel centrelines extracted by the method of [6]. Let $\mathcal{S}$ be the (centerline) curve of a source vessel and $\mathcal{R}$ the curve of the corresponding reference vessel. The curve-to-curve error $d(\mathcal{S},\mathcal{R})$ is the sum of the squared distances between all points $\mathrm{s}$ of $\mathcal{S}$ to the curve $\mathcal{R}$, $d(\mathcal{S},\mathcal{R})=\sum_{\mathrm{s}\in\mathcal{S}}d^{2}(\mathrm{s},\mathcal{R})$ [7]. The squared point-to-curve distance $d^{2}(\mathrm{s},\mathcal{R})$ is defined for every $\mathrm{s}\in\mathcal{S}$ as the squared Euclidean distance to its closest point $\mathrm{r}^{*}$ on $\mathcal{R}$, $d^{2}(\mathrm{s},\mathcal{R})=\min_{\mathrm{r}\in\mathcal{R}}\left\|\mathrm{r}-\mathrm{s}\right\|^{2}_{2}=\left\|\mathrm{r}^{*}-\mathrm{s}\right\|^{2}_{2}$. A second order approximation of the squared point-to-curve distance in discrete curves, was defined in [8] by

[TABLE]

where $\mathrm{t}(\mathrm{r}^{*})$ and $\mathrm{n}(\mathrm{r}^{*})$ are the unit tangent and normal vectors defined at $\mathrm{r}^{*}$. “.” is the scalar product. $\rho$ is the (signed) curvature radius at the point $\mathrm{r}^{*}$. $d$ is the signed distance to the closest point $\mathrm{r}^{*}$ defined by $d=\left\|\mathrm{s}-\mathrm{r}^{*}\right\|_{2}$ when $\mathrm{r}^{*}$ and $\mathrm{n}(\mathrm{r}^{*})$ lie on the same side of the curve and $d=-\left\|\mathrm{s}-\mathrm{r}^{*}\right\|_{2}$ otherwise. In practice, the curvature radii are computed once and for all iterations before the stage (4). $d^{2}(\mathrm{s},\mathcal{R})$ is estimated between the point $s$ and its closest vessel $\mathcal{R}$.

2.3 Equations of the overlap area to select the vessels

Knowing the equations of the overlap area allows to only select the vessels in this area. In each image $i\in\{1;2\}$, a circle is fitted on the boundaries of its field of view. Its centre corresponds to the image centre $\mathrm{c}_{i}=[x_{i},y_{i}]^{T}$ and its radius is denoted $r_{i}$. The circle equation of the image $i$ can be expressed by $(x-x_{i})^{2}+(y-y_{i})^{2}=r_{i}^{2}$ or by $\mathrm{x}^{T}\mathrm{C}_{i}\mathrm{x}=0$ in matrix form, where $\mathrm{x}=[x,y,1]^{T}$. The circle matrix $\mathrm{C}_{i}$ is given by:

[TABLE]

The radial distortion transforms the disk of radius $r_{i}$ into an undistorted disk of radius $r_{i}^{u}=\frac{r_{i}}{1+k_{i}r_{i}^{2}}$, whose equation is $\mathrm{x}^{T}\mathrm{C}_{i}^{u}\mathrm{x}=0$. $\mathrm{C}_{i}^{u}$ is the same matrix as $\mathrm{C}_{i}$ apart from $r_{i}$ which is replaced by $r_{i}^{u}$. Under the homography transformation $\mathrm{x}^{\prime}=\mathrm{H}\mathrm{x}$, the equation of the undistorted disk 1 becomes [5]:

[TABLE]

where $\mathrm{C}_{1}^{r}=(\mathrm{H}^{-1})^{T}\mathrm{C}_{1}^{u}(\mathrm{H}^{-1})$. $\mathrm{C}_{1}^{r}$ is the matrix of a conic [5] and in our case an ellipse. Indeed, the determinant of the $2\times 2$ top left hand block of the matrix is strictly positive $det(\mathrm{C}_{1}^{r}(1,2;1,2))=(det(A))^{2}>0$ [9, Chap. 7.5]. The overlap area between the circle $\mathrm{C}_{2}^{u}$ and the ellipse $\mathrm{C}_{1}^{r}$ is then defined by the points inside the circle and the ellipse:

[TABLE]

This equation system allows to only select the vessels inside the overlap area at each iteration of the model estimation.

3 Experiments and results

Experiments were made in order to compare the current method to others. Two datasets were used: a Public Health dataset [1] and a publicly available dataset.

3.1 Experiments in a Public Health dataset

The PH dataset is composed of 69 randomly selected patients coming from Diabetic Eye Screening Wales (DESW) programme in the United Kingdom. All patients had diabetes and different severity stages of retinopathy or maculopathy. Each of them had been screened annually for several years and 4 images were available per screening. We selected a series of 271 image pairs: (1) of sufficient quality, (2) with an approximate screening interval of one year between the examination events [1] and (3) captured when the screening service was renewing its eye fundus cameras. 63 % of the pairs were captured with different cameras - different resolutions and distortions. 10 pairs had a small overlap area of about 30 % of the superimposed image surface. All the retinal photographs were high quality and were captured according to a protocol including pupillary dilation. To assess the superimposition quality, all the registered overlap area were carefully checked by an expert according to the visual classification presented in [1]. Two categories were considered: (a) no noticeable difference (i.e. correct) and (b) noticeable difference (i.e. incorrect) with three subcategories: (b.1) differences of a small diameter vessel, (b.2) differences of the size of a large diameter vessel or (b.3) even larger. The three subcategories were grouped into a single one incorrect. Using this visual scale, we evaluated (i) the current method and we compared it to the results obtained in [1] in the same dataset for three other methods: (ii) the previous one [1] (iii) Lee et al.’s method [10] and (iv) “gdbicp” quadratic [11]. We added a comparison with a recent method (v) “REMPE” [3] based on a spherical eye assumption. Standard parameters were used.

3.2 Results in the Public Health dataset

In table 1, 96 % of the pairs are correctly superimposed (i) with the current method. This methods better registers the pairs than: (ii) the previous one, (iii) Lee et al. [10], (v) “REMPE” and (iv) “gdbicp” quadratic [11]. Figure 2 illustrates the current method (i) which successfully superimposes a pair with a small overlap, whereas the previous one (ii) fails. Our method is therefore better than the others (ii), (iii), (iv) and (v), in this PH dataset where 63 % of the pairs were captured with different cameras.

3.3 Experiments in a publicly available dataset: FIRE

We also comparee (i) the current method to (iv) “REMPE” using FIRE dataset which includes a ground truth [4]. It is composed of 129 images forming 134 image pairs and divided into 3 categories: (1) the category $\mathcal{S}$ which contains 71 pairs with high overlap and no anatomical differences; (2) the category $\mathcal{P}$ which includes 49 pairs with small overlap and no anatomical differences and (3) the category $\mathcal{A}$ which includes 14 pairs with high overlap and large anatomical changes. The ground truth was created by experts [4] who selected 10 corresponding points in each image. For both methods (i) and (v), we registered these points and we computed their mean Euclidean distances between each image of the pair.

3.4 Results in the publicly available dataset: FIRE

Registration error is given in the 2D plots of figure 3 according to the approach of [3, 4]. The horizontal axis is the mean error threshold value under which a registration is considered as successful. A step of 0.1 pixels is used for the error threshold. The vertical axis is the percentage of successful registrations. The Area Under Curve (AUC) is also computed for each curve. The higher the registration curve value (or the AUC) is, the better the registration is. The registration curves (Fig. 3), the AUC (table 2) and the mean error (table 3) are given for each category $\mathcal{S}$ , $\mathcal{P}$ or $\mathcal{A}$ and for the whole dataset. In each case, the registration curves and the AUC of (i) the current method are higher and greater than those of (v) “REMPE” (table 2). In each case, the mean error and standard deviation of (i) the current method are less than those of (v) “REMPE”. The current method is therefore better than “REMPE” in all the categories of FIRE dataset and especially for the categories $\mathcal{P}$ (with a small overlap) and $\mathcal{A}$ (with large anatomical differences). However, for these two categories the mean error remains high, 9.25 and 5.81 pixels respectively (table 3) for images of size $2912^{2}$ pixels. The superimposition can still be improved in these two categories especially for images with a small overlap (less than 50 %).

4 Conclusion and perpsectives

We have successfully achieved a new error criterion to estimate our affinity model with two radial distortions. It is based on an error between the vessels which are selected by the disk and the ellipse equations of the overlap area boundaries deduced from the model equation. Experiments have shown that our method successfully superimposes 96 % of the pairs from a PH dataset whose images are mostly acquired with different cameras. This is better than our previous method [1] and than three other state-of-the art methods [10, 3, 11]. In the publicly available dataset, FIRE [4], ours still better superimposes the images than the state-of-the-art method “REMPE” even if not all the pairs are perfectly superimposed. Nevertheless, the results show that our method is efficient for images of PH databases which are used for retinopathy screening and which present strong contrast variations and radial distortions.