Guillaume Noyel
email: guillaume.noyel@i-pri.org

Rebecca Thomas
email: r.l.thomas@swansea.ac.uk

Gavin Bhakta
email: gavin.bhakta@drssw.wales.nhs.uk

Andrew Crowder
email: andrew.crowder@drssw.wales.nhs.uk

David Owens
email: owensdr@cardiff.ac.uk

Peter Boyle
email: peter.boyle@i-pri.org

Superimposition of eye fundus images for longitudinal analysis from large public health databases

Keywords: eye fundus images, image registration, invertible model, longitudinal analysis, public health databases, radial distortion, diabetic retinopathy Mathematics Subject Classification (MSC): 68U10 Image processing 65D19 Computational issues in computer and robotic vision

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Diabetic Retinopathy (DR) as one of the major causes of visual impairment in the world represents a major public health challenge. It is a complication of both types of diabetes mellitus, which affects the light perception part of the eye (retina). DR may lead to the development of sight threatening lesions and without adequate and timely treatment the patients could lose their sight and eventually become blind (International Diabetes Federation and The Fred Hollows Foundation, 2015;[START_REF] Scanlon | A practical manual of diabetic retinopathy management[END_REF]. DR is often asymptomatic until an advanced stage, thereby screening to detect sight threatening DR at an early stage is essential which has resulted in the introduction of DR Screening services in many countries such as UK [START_REF] Harding | Grading and disease management in national screening for diabetic retinopathy in England and Wales[END_REF], USA, the Netherlands, France, etc. The commonest screening method involves acquiring eye fundus images on an annual or biennial basis.

As these DR screening programs have been in existence over several years, performing longitudinal analysis of the eye fundus images of the same patient is now possible. However, in order to accurately compare the evolution of DR over time, the images must be perfectly super-imposed.

The direct superimposition of two images of the same patient never gives good results (see Figure 1). Indeed, for two separate photographic-eye examinations the patient is never in exactly in same position and also the camera may differ. Therefore, the super-imposition method has to take into account the different causes of the deformation such as: i) the position of the patient: by taking into account rotation, translation and scaling, ii) the change of the camera: by using scaling and radial correction, iii) the projection of a 3D scene assimilated to a sphere (the retina of the eye) onto the 2D plane of the sensor of the camera: by using radial correction, iv) the radial deformation due to the optics of the camera: by using radial correction and v) the colour variability between images due to the light intensity and sensor. To perform a superimpositionalso named registrationtwo stages are regarded as necessary: a model of deformation and a matching criterion to fit the model. Registration of medical images has been a very active field for research during the last decades especially in the field of radiology. For an overview of medical image registration methods see e.g. [START_REF] Maintz | A survey of medical image registration[END_REF][START_REF] Oliveira | Medical image registration: a review[END_REF][START_REF] Viergever | A survey of medical image registration -under review[END_REF]. [START_REF] Pluim | Mutual-information-based registration of medical images: a survey[END_REF] have also written a survey about mutual-information-based registration of medical images and [START_REF] Sotiras | Deformable medical image registration: a survey[END_REF] about deformable registration. In this paper, we will focus on the registration of eye fundus images.

Several models of super-imposition between pairs of eye fundus images are in existence. The earliest methods relied on fluorescein images and based on a composition of translation, rotation and scaling -i.e. an affine transformation model (Zana and Klein, 1999a, b). The bifurcations of the vessels were used to match the points and fit the model. Another matching criterion consists in the minimisation of image intensity differences [START_REF] Cideciyan | Registration of ocular fundus images: an algorithm using cross-correlation of triple invariant image descriptors[END_REF][START_REF] Matsopoulos | Automatic Retinal Image Registration Scheme Using Global Optimization Techniques[END_REF][START_REF] Ritter | Registration of Stereo and Temporal Images of the Retina[END_REF][START_REF] Adal | A hierarchical coarse-to fine approach for fundus image registration Biomedical image registration[END_REF]. Matching by salient points using SURF -Speeded Up Robust Features -detector, was used in [START_REF] Cattin | Retina Mosaicing Using Local Features[END_REF].

Other methods are based on a similarity (i.e. a rotation and a translation) [START_REF] Matsopoulos | Automatic Retinal Image Registration Scheme Using Global Optimization Techniques[END_REF][START_REF] Lloret | th International Conference on Pattern Recognition[END_REF] and an elastic model of deformation (You et al., 2005;[START_REF] Fang | Elastic registration for retinal images based on reconstructed vascular trees[END_REF]. In [START_REF] Jian | A Partial Intensity Invariant Feature Descriptor for Multimodal Retinal Image Registration[END_REF]; [START_REF] Ghassabi | An efficient approach for robust multimodal retinal image registration based on UR-SIFT features and PIIFD descriptors[END_REF], new descriptors PIIFD (Partial Intensity Invariant Feature Descriptor) have been introduced for multimodal image superimposition between auto-fluorescence, infrared and red-free images (i.e. the green component of a colour eye fundus images) using an affine and a quadratic model of deformation. This method, powerful for multimodal image registration, has not been designed and tested for the superimposition of colour eye fundus images.

More recently, it has been shown that a quadratic model gives better registration results [START_REF] Adal | A hierarchical coarse-to fine approach for fundus image registration Biomedical image registration[END_REF][START_REF] Can | A feature-based, robust, hierarchical algorithm for registering pairs of images of the curved human retina[END_REF][START_REF] Stewart | The dual-bootstrap iterative closest point algorithm with application to retinal image registration[END_REF][START_REF] Chanwimaluang | Hybrid retinal image registration[END_REF][START_REF] Ryan | Registration of digital retinal images using landmark correspondence by expectation maximization Image[END_REF]. The difficulty inherent in these models is to estimate their parameters. To overcome such a limitation, a radial distortion model has been introduced by [START_REF] Lee | Feature-based pairwise retinal image registration by radial distortion correction[END_REF] and compared to previous methods in [START_REF] Lee | Objective and expert-independent validation of retinal image registration algorithms by a projective imaging distortion model[END_REF]. It consists of adding a radial model to the affine transformation in order to correct the effects of radial distortion due to the geometry of the camera and of the eye. Superimposition of eye fundus images have also been performed in the tri-dimensional space using a model of near planar surface [START_REF] Yuping | Retinal image registration from 2D to 3D[END_REF] or an ellipsoid model [START_REF] Hernandez-Matas | IEEE 38th Annual International Conference of the Engineering in Medicine and Biology Society (EMBC)[END_REF].

However, for all these methods the superimposition of eye fundus images is performed with colour images acquired by the same camera during the same examination. Superimposing colour images acquired by different cameras with at least a year of interval and stored in large screening databases still remains a challenging problem. Two new issues are appearing: i) the radial distortion must be corrected for the two different lenses, because the camera has changed and ii) the colour of the eye may have changed due to a different source of light and to an evolution of the anatomy of the eye (e.g. cataract removal, evolution of the retinopathy to a more severe stage, etc.)

Currently, for colour eye fundus images, no existing method has been designed to address the problem of the change of camera lenses, camera and light for retinal examination used for DR screening (i.e. at a year of interval). In addition, in the existing methods, even when one radial distortion is corrected by a radial distortion correction or by a second order model, the registration criterion requires a similar and uniform intensity (i.e. colour) of the images, particularly if the extraction of anatomical features such as the vessels is needed. These methods are working in laboratory conditions (several acquisitions on the same day with the same equipment) but failed when they are used on public health databases constituted by routine screening across years of patients with Diabetes.

In this paper, our contribution has been to address this challenge by presenting a robust superimposition method designed for longitudinal screening of large public health image based databases. After having taken into account the variation of colour, we introduce a new model of image registration based on two radial distortions and a homography. Importantly, this model is invertible which is useful for image deformation (i.e. image warping) and makes it suitable for fast computation, a key factor for an efficient analysis of large databases. We also introduce linear equations to estimate the parameters of the model in a fast and simple way. Finally, we provide a complete framework for robust image superimposition. The registration criterion is based on a standard salient point extraction, which does not rely on the extraction of anatomical features.

The paper is organised as follows: after presenting our method to superimpose pairs of images, it will be validated on a simulated montage and using different patient databases. It will be compared with two well-established methods of image registration. Then, the results will be interpreted and discussed before the conclusion.

Methods

During a photographic eye examination, eye fundus (retina) two images of both eyes are acquired in two positions i.e. a 45 degree "nasal" and "macular" field (Figure 2). As patients with diabetes are annually screened this make series of images available for longitudinal analysis. This is a different problem from large mosaics of eye fundus images acquired during the same examination with the same camera [START_REF] Chanwimaluang | Hybrid retinal image registration[END_REF]. Consequently, the aim of this study was to develop a superimposition method of images in the same positions while being captured during two different exams and often with different cameras and resolutions.

A schematic description of the study is represented in Figure 3. The different stages have been designed to provide efficient solutions to the superimposition of images acquired for practical screening. Between two examinations, the camera and lighting might have been changed producing differences in colour requiring a pre-processing stage consisting of (1) normalizing the colour of the eye fundus image [START_REF] Noyel | Contrast enhancement of eye fundus images[END_REF]. In addition, the position, resolution and radial distortion between images might be different requiring an image transformation estimated thanks to the (2) extraction of characteristic points in pairs of images (3) a matching procedure, (4) the use of a model correcting radial distortion of both images and (5) the estimation of the parameters of the model by a robust optimisation. Eventually, the method is validated (6) on a simulated montage and the superimpositions of the image of the database are verified (7) and compared to state of the art methods.

In this paper, the image pairs of the same eye in the same position, "nasal" or "macular", are registered independently. The brightness of eye fundus images is non-uniform due to various reasons: disease such as cataract, motion of the patient, acquisition conditions and differences in absorptions of the light in the eye [START_REF] Walter | Application de la morphologie mathématique au diagnostic de la rétinopathie diabétique à partir d'images couleur[END_REF][START_REF] Walter | Automatic Analysis of Color Fundus Photographs and Its Application to the Diagnosis of Diabetic Retinopathy[END_REF]. Some parts of the images appear as bright while others are dark. Moreover, the possible change of the eye fundus camera between two separate examinations may contribute to a change in the colour between two images of the same eye (Figure 4). We have used a method [START_REF] Noyel | Contrast enhancement of eye fundus images[END_REF] to correct the variations of colour contrast between the images. The method is based on a colour model consistent with the physical principles of image formation. The contrast of the dark or the bright elements are adjusted in a way that provides a similar colour aspect to lesions such as micro-aneurysms or to anatomical structures such as vessels or veins. Results can be seen in Figure 4. Our colour methods goes further than previous existing grey-level methods using polynomial adjustment [START_REF] Walter | Automatic detection of microaneurysms in color fundus images[END_REF], mathematical morphology [START_REF] Zhang | Exudate detection in color retinal images for mass screening of diabetic retinopathy[END_REF] or histogram equalisation [START_REF] Zuiderveld | Contrast limited adaptive histogram equalization[END_REF].

Extraction of characteristic points

The pre-processing is followed by the extraction of several salient points (Figure 5) using the Scale-Invariant Feature Transform -SIFT -algorithm [START_REF] Lowe | Distinctive Image Features from Scale-Invariant Keypoints[END_REF][START_REF] Vedaldi | VLFEat: An open and portable library of computor vision algorithms[END_REF]. The SIFT algorithm has been designed to be robust to the variation of observation angle and to some variations in lighting. Briefly SIFT consists of extracting key points based on a multiscale analysis. Then, series of descriptors are computed for each salient point. These descriptors are used for point matching.

Point matching with a preliminary estimate of the image deformation

Point matching is a necessary step before the estimation of the model of deformation. With the SIFT method to extract characteristic points, [START_REF] Lowe | Distinctive Image Features from Scale-Invariant Keypoints[END_REF] has proposed a matching method. However, it was not robust enough to estimate the model of deformation in the database. Therefore, a three-step procedure of matching has been proposed:

(a) A first matching by Lowe's method followed by a refined selection of the correspondence vectors according to their size and orientation (see hereinafter). (b) A preliminary estimate of the image deformation by a homography. The position of the characteristic points is transformed using a homography (see hereinafter). (c) Step (a) is run a second time using the transformed points.

Selection of the correspondence vectors

First of all, the region of interest of both images are extracted and joint together after resizing and padding to have the same vertical size (Figure 5). As the image deformation is closed to a similarity (i.e. rotation, translation and scaling), the length of ratios and angles are invariant by similarity [START_REF] Hartley | Multiple view geometry in computor vision[END_REF]. This property is used to remove the correspondence vectors 𝑣 between inconsistent matched points. With the two images joint together, a twostep selection is performed on the lengths 𝑙 and orientations 𝜃 of the correspondence vectors:

(i) only the vectors 𝑣 whose length 𝑙 𝑣 and orientation 𝜃 𝑣 is in the interval {|𝑙 𝑣 -𝐸{𝑙 𝑣 }| ≤ 𝜎{𝑙 𝑣 } 𝑎𝑛𝑑 |𝜃 𝑣 -𝐸{𝜃 𝑣 }| ≤ 5° } are kept. 𝐸{} is the mean and 𝜎{} the standard deviation of a variable. (ii) Among the selected vectors 𝑣 ̃, only the vectors whose length 𝑙 𝑣 ̃ and orientation 𝜃 𝑣 ̃ is in the interval {|𝑙 𝑣 ̃-𝐸{𝑙 𝑣 ̃}| ≤ max(3𝜎{𝑙 𝑣 ̃}, 5% × 𝑦𝑠𝑖𝑧𝑒) 𝑎𝑛𝑑 |𝜃 𝑣 ̃-𝐸{𝜃 𝑣 ̃}| ≤ max (5°, 𝜎{𝜃 𝑣 ̃})} are kept. 𝑦𝑠𝑖𝑧𝑒 is the number of lines in the image.

Notice: The value 5° has been empirically selected and perfectly works for all the processed images.

Estimation of a homography

The affine homography (or affine transformation) 𝐻 is defined as: The estimation of a homography is performed using a robust maximum likelihood estimation with the RANdom SAmple Consensus (RANSAC) algorithm applied on the "linear normalised DLT (Direct Linear Transform)" algorithm [START_REF] Hartley | Multiple view geometry in computor vision[END_REF]). This is the initialisation step of the "gold standard algorithm" of [START_REF] Hartley | Multiple view geometry in computor vision[END_REF]. This initialisation algorithm estimates around 100 homographies on randomly selected 4 pairs of points. Then, the homography with the minimum error, when transforming the matched points, is selected. This allows a removal of the incorrect matchings. However, as there is also a radial distortion in the image, the deformation is not entirely modelled by a homography. Therefore, some homographies must be discarded. In particular, those with a scaling factor on the 𝑥 and 𝑦 axis with a relative difference greater than 1%, because the sensor resolution is almost the same in 𝑥 and in 𝑦. For this purpose, several estimates (until 50) using the initialisation algorithm are performed if the relative difference between the scaling factors is greater than 1%. If the value of 1% is never reached, then the homography with the smallest relative difference between the scaling factors is kept. Finally, a non-linear estimation is performed with the algorithm of Levenberg-Marquardt [START_REF] More | The Levenberg-Marquardt algorithm: implementation and theory[END_REF][START_REF] Bonnans | Numerical optimization[END_REF]. It serves as a good initialisation of the estimation of the complete deformation model.

𝐻 = [ 𝐴 𝑇 𝑂 𝑇 1 ] = [

Model of deformation

The model of deformation ensures a correct superimposition between the images. Several deformations are taken into account: (i) the difference in terms of positions of the eye between a pair of images will be corrected by an affine transformation (i.e. and homography) and (ii) the radial deformations due to the projection of the eye into the camera and due to the optics of the camera [START_REF] Hartley | Multiple view geometry in computor vision[END_REF]) will be corrected using a radial transformation. For this purpose, [START_REF] Lee | Feature-based pairwise retinal image registration by radial distortion correction[END_REF]; [START_REF] Lee | Retinal image mosaicing using the radial distortion correction model[END_REF], have proposed a model coupling a unique radial transformation for both images and a homography. [START_REF] Lee | Objective and expert-independent validation of retinal image registration algorithms by a projective imaging distortion model[END_REF] have made the comparison with two other second order models.

Their approach was extended by defining a model with one homography 𝐻 and two radial distortions, i.e. one for each image. Indeed, the camera may have changed between screening exams on a large number of patients.

The radial distortion due the background of the sphere surface of the eye and of the radial distortion of the camera was modelled by a division model [START_REF] Fitzgibbon | Simultaneous linear estimation of multiple view geometry and lens distortion[END_REF] in the following way:

𝑃 ̅ 𝑑 = (1 + 𝑘(𝑟 𝑑 ) 2 ) . 𝑃 ̅ 𝑢 (2)
𝑃 𝑑 ∈ ℝ 2 are the distorted coordinates in the original (i.e. distorted) image. 𝑃 𝑢 ∈ ℝ 2 are the undistorted coordinates in the undistorted image. 𝑃 ̅ 𝑑 ∈ ℝ 2 are the distorted coordinates centred on the image centre 𝐶: 𝑃 ̅ 𝑑 = 𝑃 𝑑 -𝐶. 𝑃 ̅ 𝑢 ∈ ℝ 2 are the undistorted coordinates centred on the image centre 𝐶: 𝑃 ̅ 𝑢 = 𝑃 𝑢 -𝐶. 𝑟 𝑑 = ‖𝑃 𝑑 -𝐶‖ = ‖𝑃 ̅ 𝑑 ‖ ∈ ℝ, is the distance of the deformed coordinates 𝑃 𝑑 from the optic centre 𝐶 (i.e. the image centre). 𝑘 is a real parameter of distortion. Its normalised version, 𝑘 ̃= 𝑘(1 + ‖𝐶‖ 2 ), is in the interval [-0.2 ; 0.2].

The model was named division model because the distorted coordinates were divided by the radial distortion 𝑃 ̅ 𝑑 (1 + 𝑘(𝑟 𝑑 ) 2 ) ⁄ = 𝑃 ̅ 𝑢 . The distorted image corresponds to the original image and the undistorted image is the image after the correction of radial distortion.

Given 𝑃 1 𝑑 and 𝑃 2 𝑑 the coordinates of the points in the original (i.e. deformed) images 1 and 2, 𝑘 1 and 𝑘 2 the distortion parameters, 𝐶 1 and 𝐶 2 the image centres, the model mapping image 1 into image 2 is defined as follows:

𝑃 ̅ 2 𝑑 (1 + 𝑘 2 (𝑟 2 𝑑 ) 2 ) + 𝐶 2 = 𝐻 [ 𝑃 ̅ 1 𝑑 (1 + 𝑘 1 (𝑟 1 𝑑 ) 2 ) + 𝐶 1 ] (3)
If the same camera is used to acquire both images, the distortion parameters are equal, 𝑘 1 = 𝑘 2 = 𝑘, and the model corresponds to the model of [START_REF] Lee | Feature-based pairwise retinal image registration by radial distortion correction[END_REF]. However, in our approach, the model is estimated after having extracted and matched the points in the pair of original images (target and reference). Therefore, the radial distortion correction is performed after the detection of the feature correspondence points in the original (i.e. distorted) image. It has been programmed in the subsequent way.

After point matching, the model is estimated by minimising the deformation error 𝑒𝑟𝑟 between the undistorted points in the image 2 estimated by the model,

(𝑃 ̅ 2 𝑢 ) 𝑒𝑠𝑡 = 𝐻 [ 𝑃 ̅ 1 𝑑 (1+𝑘 1 (𝑟 1 𝑑 ) 2 )
+ 𝐶 1 ] -𝐶 2 , and the undistorted points extracted in the image 2, 𝑃 ̅ 2 𝑢 = 𝑃 ̅ 2 𝑑 / (1 + 𝑘 2 (𝑟 2 𝑑 ) 2 ). The deformation error is 𝑒𝑟𝑟 = ‖(𝑃 ̅ 2 𝑢 ) 𝑒𝑠𝑡 -𝑃 ̅ 2 𝑢 ‖.

Estimation of the model parameters

The model parameters are estimated in (Figure 6) by a different method of [START_REF] Lee | Feature-based pairwise retinal image registration by radial distortion correction[END_REF][START_REF] Abramoff | Optimal registration of multiple deformed images using a physical model of the imaging distortion[END_REF]. In particular, during the initialisation step, the type of the model is selected: either a single distortion model, if the images are acquired by the same camera, or a two distortions model if two cameras are used. In addition, the radial distortion is estimated after the homography without needing any preliminary initialisation by a calibration of the camera [START_REF] Hartley | Multiple view geometry in computor vision[END_REF]. Moreover, we have introduced linear estimators at each step of the process.

Figure 6. Flowchart of the superimposition method based on a homography and a radial distortion process.

The estimation process is as follows (Figure 6). The homography and the radial distortions are initialised using linear estimators. Then, an iterative estimation is performed independently for the homography and the radial distortions until the deformation error has converged. In order to ensure the error has converged, at each step a linear estimator is used as an initialiser of a non-linear (iterative) optimiser such as Levenberg-Marquardt [START_REF] More | The Levenberg-Marquardt algorithm: implementation and theory[END_REF][START_REF] Bonnans | Numerical optimization[END_REF]. Indeed, without initialisation, the non-linear estimator may converge towards a local minimum of the error [START_REF] Hartley | Multiple view geometry in computor vision[END_REF]) (p. 110). The convergence criterion is defined by: 𝑒𝑟𝑟 is the mean error at iteration 𝑛, 𝜀 = 0.01 is the tolerance on the error, 𝑡𝑜𝑙 = 0.01 is a tolerance on the relative error between iterations and 𝑀𝑎𝑥𝐼𝑡𝑒𝑟 = 100 is the maximum number of iterations. For each iteration, the normalised radial distortion parameters 𝑘 ̃ must be in the interval [-0.2 ; 0.2]. If not, the algorithm stops and the model estimate with the smallest error is selected.

𝑐𝑜𝑛𝑣𝑒𝑟𝑔𝑒𝑛𝑐𝑒 = [𝑒𝑟𝑟 𝑛 < 𝜀] 𝑎𝑛𝑑 [ 𝑒𝑟𝑟 𝑛 -𝑒𝑟𝑟 𝑛-1 𝑒𝑟𝑟 𝑛-1 < 𝑡𝑜𝑙] 𝑎𝑛𝑑 [𝑘 ̃ 𝑜𝑟 {𝑘 ̃1𝑎𝑛𝑑 𝑘 ̃2} ∈ [-0.2 ; 0.2]] 𝑎𝑛𝑑 [𝑛 < 𝑀𝑎𝑥𝐼𝑡𝑒𝑟]
The final optimisation is performed jointly on the homography and the radial distortion with a linear estimator followed by a non-linear optimiser such as Levenberg-Marquardt [START_REF] More | The Levenberg-Marquardt algorithm: implementation and theory[END_REF][START_REF] Bonnans | Numerical optimization[END_REF] for a single distortion or for two distortions the "trust region method" [START_REF] Moré | Recent Developments in Algorithms and Software for Trust Region Methods[END_REF][START_REF] Bonnans | Numerical optimization[END_REF] with the bounds [-0.2 ; 0.2] for the radial distortion parameters.

Notice: Even, if linear equations are introduced to estimate the model parameters, a final non-linear estimator is used for two reasons. Firstly, an estimate of the model is performed using the (non-linear) parametrisation of the matrix 𝐴 into a set of two angles and two scaling factors [START_REF] Hartley | Multiple view geometry in computor vision[END_REF]) (p. 40), which have a direct physical interpretation. Secondly, the location of the correspondence points may be re-estimated jointly with the model parameters using the "Sampson error" or the "gold standard error" (Hartley and Zisserman, 2004) (p. 114). In this paper, the position of the correspondence points has not been re-estimated. It will be studied in a future paper.

Estimation of the model with two radial distortions.

When a different camera is used for each image of a pair, two radial distortions must be estimated, 𝑘 1 and 𝑘 2 . Let us define linear estimators to initialise the non-linear estimators of the model parameters.

Linear estimator of the radial distortion parameters 𝑘 1 and 𝑘 2

Equations ( 1) and (3) imply that:

𝑘 1 𝑘 2 [(𝑟 1 𝑑 𝑟 2 𝑑 ) 2 . 𝐷] + 𝑘 1 [𝑟 1 𝑑 2 . 𝐷 + 𝑟 1 𝑑 2 . 𝑃 ̅ 2 𝑑 ] + 𝑘 2 [𝑟 2 𝑑 2 . 𝐷 -𝑟 2 𝑑 2 . 𝐴𝑃 ̅ 1 𝑑 ] = -[𝑃 ̅ 2 𝑑 + 𝐷 -𝐴𝑃 ̅ 1 𝑑 ] (4)
with 𝐷 = 𝐶 2 -𝐴𝐶 1 -𝑇. Equation ( 4) is a quadratic equation with its unknowns 𝑘 1 and 𝑘 2 :

𝑘 1 𝑘 2 [𝑀 12 ] + 𝑘 1 [𝑀 1 ] + 𝑘 2 [𝑀 2 ] + [𝑀 0 ] = 0.
The variables 𝑀 𝑖 , 𝑖 = {0,1,2,12}, have two components 𝑀 𝑖 = [𝑀 𝑖 𝑥 , 𝑀 𝑖 𝑦 ] 𝑡 , 𝑀 𝑖 𝑥 , 𝑀 𝑖 𝑦 ∈ ℝ 𝑁 . For each component, the equation is verified for each of the 𝑁 pairs of matched points. This leads to a set of two systems of quadratic equations. The set of equations corresponds to the intersection of two conics which as a unique pair of solutions according to Bernstein's theorem for bivariate polynomials (Sturmfels) (p. 37). This result can be explained in a more intuitive way; between a pair of matched points, there is one pair of solutions per equation (i.e. two pairs).

As the solution must be valid for both systems of equations, there is only a unique solution pair (𝑘 1 , 𝑘 2 ) valid for all the matched points. The other solutions are only valid for a single pair of matched points. The unique solution pair (𝑘 1 , 𝑘 2 ) satisfies both systems of equations together and can be obtained by elimination of the second order term:

𝑘 1 𝑘 2 + 𝑀 1 𝑥 𝑀 12 𝑥 𝑘 1 + 𝑀 2 𝑥 𝑀 12 𝑥 𝑘 2 + 𝑀 0 𝑥 𝑀 12 𝑥 = 0 𝑘 1 𝑘 2 + 𝑀 1 𝑦 𝑀 12 𝑦 𝑘 1 + 𝑀 2 𝑦 𝑀 12 𝑦 𝑘 2 + 𝑀 0 𝑦 𝑀 12 𝑦 = 0 } ⇒ 𝑘 1 [ 𝑀 1 𝑥 𝑀 12 𝑥 - 𝑀 1 𝑦 𝑀 12 𝑦 ] + 𝑘 2 [ 𝑀 2 𝑥 𝑀 12 𝑥 - 𝑀 2 𝑦 𝑀 12 𝑦 ] = -[ 𝑀 0 𝑥 𝑀 12 𝑥 - 𝑀 0 𝑦 𝑀 12 𝑦 ]. (5)
The unique solution pair (𝑘 1 , 𝑘 2 ) of the equation ( 5) is determined using the method of least squares. When the camera is the same for both images, only one radial distortion, 𝑘 = 𝑘 1 = 𝑘 2 , needs to be estimated in the equation (5).

Linear estimator of the radial distortion parameters 𝑘 1 and 𝑘 2 and the homography

From equations (1) and (3), we obtain:

𝑘 1 𝑘 2 𝐷 [(𝑟 1 𝑑 𝑟 2 𝑑 ) 2 ] + 𝑘 1 𝐷 [𝑟 1 𝑑 2 ] + 𝑘 1 [𝑟 1 𝑑 2 . 𝑃 ̅ 2 𝑑 ] + 𝑘 2 𝐷 [𝑟 2 𝑑 2 ] -𝑘 2 𝐴 [𝑟 2 𝑑 2 . 𝑃 ̅ 1 𝑑 ] -𝐴[𝑃 ̅ 1 𝑑 ] + 𝐷[1] = -[𝑃 ̅ 2 𝑑 ] (6)
Equation ( 6) is a linear equation with its variables in the bracket. It is valid for each component (𝑥 and y) of the vectors

𝑃 ̅ 1 𝑑 = [𝑃 ̅ 1 𝑑,𝑥 , 𝑃 ̅ 1 𝑑,𝑦 ] 𝑡 , 𝑃 ̅ 2 𝑑 = [𝑃 ̅ 2 𝑑,𝑥 , 𝑃 ̅ 2 𝑑,𝑦 ]
𝑡 and 𝐷 = [𝐷 𝑥 , 𝐷 𝑦 ] 𝑡 , leading to a system of two equations:

{ 𝑘 1 𝑘 2 𝐷 𝑥 [(𝑟 1 𝑑 𝑟 2 𝑑 ) 2 ] + 𝑘 1 𝐷 𝑥 [𝑟 1 𝑑 2 ] + 𝑘 1 [𝑟 1 𝑑 2 𝑃 ̅ 2 𝑑,𝑥 ] + 𝑘 2 𝐷 𝑥 [𝑟 2 𝑑 2 ] -𝑘 2 𝑎 11 [𝑟 2 𝑑 2 𝑃 ̅ 1 𝑑,𝑥 ] -𝑘 2 𝑎 12 [𝑟 2 𝑑 2 𝑃 ̅ 1 𝑑,𝑦 ] -𝑎 11 [𝑃 ̅ 1 𝑑,𝑥 ] -𝑎 12 [𝑃 ̅ 1 𝑑,𝑦 ] + 𝐷 𝑥 [1] = -[𝑃 ̅ 2 𝑑,𝑥 ] 𝑘 1 𝑘 2 𝐷 𝑦 [(𝑟 1 𝑑 𝑟 2 𝑑 ) 2 ] + 𝑘 1 𝐷 𝑦 [𝑟 1 𝑑 2 ] + 𝑘 1 [𝑟 1 𝑑 2 𝑃 ̅ 2 𝑑,𝑦 ] + 𝑘 2 𝐷 𝑦 [𝑟 2 𝑑 2 ] -𝑘 2 𝑎 21 [𝑟 2 𝑑 2 𝑃 ̅ 1 𝑑,𝑥 ] -𝑘 2 𝑎 22 [𝑟 2 𝑑 2 𝑃 ̅ 1 𝑑,𝑦 ] -𝑎 21 [𝑃 ̅ 1 𝑑,𝑥 ] -𝑎 22 [𝑃 ̅ 1 𝑑,𝑦 ] + 𝐷 𝑦 [1] = -[𝑃 ̅ 2 𝑑,𝑦 ] (7)
The 20 variables -in the brackets -correspond to the data. The 18 parameters of the set of both equations are estimated using the method of least squares applied to the set of both equations. As there are multiple relations between the 18 estimator parameters, only the 8 parameters of the model -𝑘 1 , 𝑘 2 and 𝐻 (through the intermediary of 𝐴 and 𝐷) -are unknown. An admissible solution for 𝑘 1 , 𝑘 2 and 𝐻 is deduced from the 18 parameters of the linear estimator. The value of 𝑘 1 is directly obtained (3 rd term of the 1 st equation) as well as those of 𝐷 𝑥 and 𝐷 𝑦 (9 th terms of both equations) and 𝑎 𝑖𝑗 (7 th and 8 th terms of both equations). The value of 𝑘 2 is deduced, e.g. from the 4 th term of the 1 st equation. Then, this admissible solution is introduced as the initialisation of the non-linear estimation.

Image warping

In order to analyse a large database, a fast algorithm of image warping is needed. Forward warping is time consuming and so we therefore use inverse warping. However, the registration model needs to be invertible [START_REF] Wolberg | Digital image warping[END_REF]. The radial distortion is modelled in equation ( 3) by a division model [START_REF] Fitzgibbon | Simultaneous linear estimation of multiple view geometry and lens distortion[END_REF]. [START_REF] Wonpil | An embedded camera lens distortion correction method for mobile computing applications[END_REF] and [START_REF] Park | Lens distortion correction using ideal image coordinates[END_REF] have computed an approximate transformation for a standard distortion method. Here, the exact inversion of the division model is computed.

Given 𝑟 𝑢 = ‖𝑃 𝑢 -𝐶‖ = ‖𝑃 ̅ 𝑢 ‖ ∈ ℝ, the distance of the undistorted coordinates 𝑃 𝑢 from the optic centre 𝐶, using equation ( 2), gives:

𝑟 𝑑 = (1 + 𝑘(𝑟 𝑑 ) 2 ) 𝑟 𝑢 (8)
Equations ( 2) and ( 8), implies that:

𝑃 ̅ 𝑢 = 𝑟 𝑢 𝑟 𝑑 𝑃 ̅ 𝑑 = 𝑊(𝑃 ̅ 𝑑 ) (9)
In order to use invert warping, it is necessary to determine 𝑊 -1 transforming the undistorted points 𝑃 𝑢 into the distorted points 𝑃 𝑑 . From equation ( 9), it is equivalent to determine 𝑟 𝑢 knowing 𝑟 𝑑 .

Equation ( 8) is equivalent to: 𝑘𝑟 𝑢 𝑟 𝑑 2 -𝑟 𝑑 + 𝑟 𝑢 = 0, a second order equation in 𝑟 𝑑 . In the demonstration hereinafter, we show that its discriminant, 𝛥 = 1 -4𝑘𝑟 𝑢 2 , is strictly positive and that the inverse transformation 𝑊 -1 corresponds to the root:

𝑟 𝑑 = 𝑊 -1 (𝑟 𝑢 ) = 1 -√1 -4𝑘𝑟 𝑢 2 2𝑘𝑟 𝑢 (10)
Demonstration: let us demonstrate (i) the positivity of the discriminant and (ii) the choice of the selected root.

We have:

𝑟 𝑢 = ‖𝑃 𝑢 -𝐶‖ ≤ ‖𝐶‖ |𝑘| = |𝑘 ̃| 1+‖𝐶‖ 2 ≤ 2.10 -1 1+‖𝐶‖ 2 } ⇒ |𝑘|𝑟 𝑢 = |𝑘 ̃| 1+‖𝐶‖ 2 ‖𝑃 𝑢 -𝐶‖ ≤ 2. 10 -1 ‖𝐶‖ 1+‖𝐶‖ 2 ≤ 2.10 -1 ‖𝐶‖ which implies that 1 2|𝑘|𝑟 𝑢 ≥ ‖𝐶‖ 4.
10 -1 and 4|𝑘|𝑟 𝑢 2 ≤ 8.10 -1 < 1 .

Therefore, we have: 𝛥 = 1 -4𝑘𝑟 𝑢 2 > 0 (i) q.e.d.

When the discriminant is positive, there exists two roots: a) 1+√Δ 2𝑘𝑟 𝑢 ≥ ‖𝐶‖ 4.10 -1 = 2.5‖𝐶‖, which corresponds to a solution where the pixel locations are rejected outside a given radius, 2.5‖𝐶‖, This is an aberrant solution. b) 1-√Δ 2𝑘𝑟 𝑢 ≤ ‖𝐶‖ 4.10 -1 = 2.5‖𝐶‖, which corresponds to a solution where the pixel locations are inside a given radius, 2.5‖𝐶‖. This is the correct solution. (ii) q.e.d. Therefore, the transformation used is invertible (Equation ( 10)). An invertible image warping method compared to a non-invertible method dramatically reduces the computation time [START_REF] Wolberg | Digital image warping[END_REF]. We experimented a reduction from about 10 minutes to a few seconds (about 5 seconds) on a standard computer using Matlab (16GB RAM, processor Intel i7-4702HQ, 2.20GHz).

Experimental validation

Validation by a simulated montage

A simulated montage was created by superimposing two eye fundus images and deforming them according to the method presented by [START_REF] Lee | Objective and expert-independent validation of retinal image registration algorithms by a projective imaging distortion model[END_REF]. Real registered images were selected with an overlap percentage of 80% corresponding to the most common case in longitudinal databases. No modification of colour was done to the images. After adding equally spaced landmarks (i.e. the ground truth), the images were cut and deformed according to the model of [START_REF] Lee | Objective and expert-independent validation of retinal image registration algorithms by a projective imaging distortion model[END_REF]. An affine transformation has been used, with rotation scaling and shearing. Then, the image was modified by a projective distortion. The radius of the eye ball is equal to the ratio between the radius of the disk of the image divided by the observation angle of the camera (45 degrees). Then, the images were registered and the error was measured between the landmarks after registration and their true position (i.e the ground truth).

The quality of image superimposition was evaluated through the simulated montage for three approaches: (i) our method with a single radial distortion, (ii) our method with two radial distortions and for (iii) a state-of the art method using a quadratic model "gdbicp" [START_REF] Yang | Registration of challenging image pairs: initialization, estimation, and decision[END_REF]. The ground truth of the simulated montage allows a comparison between our model and the quadratic model "gdbicp". The duration of each superimposition method is also measured.

Validation with a public health database

Our method needs to be visually evaluated and compared to other state-of-the-art methods in public health database used for retinopathy screening. In order to assess the evolution of Diabetic Retinopathy several screening programs in the world are in existence. Trials have been performed in a database of 69 randomly selected patients coming from the Diabetic Eye Screening Wales (DESW) program in the United Kingdom [START_REF] Thomas R L, Dunstan | Incidence of diabetic retinopathy in people with type 2 diabetes mellitus attending the Diabetic Retinopathy Screening Service for Wales: retrospective analysis[END_REF]. All patients had diabetes and different severity stages, including no signs, of retinopathy or maculopathy. All the selected patients have been screened annually for several years, seven years in average. For each of them we have kept two series of two examinations with an approximate screening interval of one year between the examination events. For each event exam, four images are available, two positions (nasal and macular) per eye. There were two series of images acquired for different years, with the first series are made up of 271 pairs of sufficient image quality (63% are acquired with different cameras) and the second series included of 268 pairs (9% are acquired with different cameras). In the first series, 11 pairs have a small overlapping area (about 30 % of the surface of the registered pair). For each position, we have performed the superimposition of the images between the two different examinations. All the retinal photographs were high quality images made according to a protocol including pupillary dilation.

In addition, in order to show that our method is not only useful on high quality images acquired following pupillary dilation but also on low quality images without dilation of the pupils prior to photography, a second validation test is performed to superimpose the images of low/poor quality for 5 patients with diabetes. The quality of image superimposition is checked for the 20 pairs of images (i.e. 4 per patient). The acquisitions conditions were significantly harsher compared to the high quality images and the quality of images was quite heterogeneous in part due to the lack of pupillary dilatation prior to photography. In this database, the patients with diabetes have two retinal examination with at least a year of interval using the same camera.

As there is no reference registration for the databases (i.e. no ground truth), for all pairs of images the superimposition has been visually checked by the same expert. The complete overlapping area of each pair has been carefully checked. The classification has been done in two categories: 1) no noticeable difference (i.e. correct) and 2) noticeable difference (i.e. incorrect). The category "incorrect" includes three sub categories: i) differences of the size of a small diameter vessel, ii) differences of the size of a large diameter vessel or iii) even larger. In this paper, the three subcategories have been grouped into a single category "incorrect". For each set of about 270 images, the duration of one check for one model is longer than an hour and a half. For each method and each pair of each set a first check has been done and a second check has been performed for verification by the same person. For almost all the pairs the results have been in accordance between both checks. In the rare cases of discrepancy, an adjudication has always been taken in the most restrictive way (i.e. "incorrect"). Therefore, the total duration of the check process was longer than thirty hours for the five methods tested hereinafter.

A visual evaluation has been performed on the two series of high quality images and on the series of low quality images in order to compare our approach to two other well established state-of-the-art whose code is publicly available: "gdbicp" [START_REF] Stewart | The dual-bootstrap iterative closest point algorithm with application to retinal image registration[END_REF][START_REF] Yang | Registration of challenging image pairs: initialization, estimation, and decision[END_REF] or that we have reprogrammed [START_REF] Lee | Feature-based pairwise retinal image registration by radial distortion correction[END_REF]. In particular, we have tested five methods (see Table 3): (i) our method on the normalised images with an automatic choice of one or two radial distortions according to the number of cameras used, (ii) our model with a single radial distortion for all the original image pairs without colour stabilisation, in order to be similar to the state-of-the-art-method of [START_REF] Lee | Feature-based pairwise retinal image registration by radial distortion correction[END_REF], (iii) the state-of-the-art-method "gdbicp" with a quadratic model [START_REF] Stewart | The dual-bootstrap iterative closest point algorithm with application to retinal image registration[END_REF] on the original images, (iv) "gdbicp" with a radial model on the original images [START_REF] Yang | Registration of challenging image pairs: initialization, estimation, and decision[END_REF] and "gdbicp" with a quadratic model [START_REF] Stewart | The dual-bootstrap iterative closest point algorithm with application to retinal image registration[END_REF] on the normalised images.

Notice: To facilitate the comparison between (i) our model and (ii) the model similar to [START_REF] Lee | Feature-based pairwise retinal image registration by radial distortion correction[END_REF], we have performed the model fitting with our method for both models, using the characteristics points extracted by SIFT, whereas [START_REF] Lee | Feature-based pairwise retinal image registration by radial distortion correction[END_REF] use the centreline of the vessels and the branch centres for their model.

Results

Validation by a simulated montage

Using our method, with a single distortion, the mean registration error is of 0.86 pixels (standard deviation 1.75 pixels) in images of size 1568 x 2352 pixels (Figure 7) with vessels of maximum diameter of 30 pixels. The relative error respective to the image diagonal is 0.03%, and respective to the vessels is 2.9 %. With two distortions, the mean registration error is of 0.92 pixels (standard deviation 1.94 pixels) and the relative error is 0.03 % respective to the image and 3.1 % respective to the vessels (see Table 1). In the Figure 7, the error is mainly located on the external part of the superimposed pair of images. Our method is compared to another state-of-the-art method "gdbicp" which superimposes eye fundus images with a quadratic model [START_REF] Stewart | The dual-bootstrap iterative closest point algorithm with application to retinal image registration[END_REF]. Using the simulated montage, the mean registration error is of 2.44 pixels (standard deviation 1.64 pixels). As a reminder, the model of deformation proposed in this paper with a single radial distortion corresponds to the model of [START_REF] Lee | Feature-based pairwise retinal image registration by radial distortion correction[END_REF], however, with a different method to estimate the parameters.

As a speed test, the duration of the different superimposition methods is measured with the images of the simulated montage (see Table 2). The experiments have been done on a recent computer (16GB RAM, processor Intel i7-4702HQ, 2.20GHz). The inverse warping method with one or two distortions lasts approximately 5.5 s using Matlab. Homography and 1 radial distortion (normalised images) (our method) 64.5

Homography and 2 radial distortions (normalised images) (our method) 71.7

"gdbicp" quadratic (initial images) 169

Validation with a public health database and comparison to other state-of-the-art methods

A visual evaluation is essential on the public health databases for the five methods previously described. Preliminary results have been presented in [START_REF] Noyel | American Diabetes Association 77th Scientific Sessions -ADA2017[END_REF]. (i) First of all, our superimposition method is applied on standardised images. If the same camera acquired both images of a pair, our model is used with a single radial distortion. If both cameras were used, the model with two radial distortions is applied. In Figure 1 In the first series of 271 pairs, 2 pairs have small differences in the external part of the superimposition. These differences are of the size of the diameter of a vessel (Figure 10). When the percentage of overlapping surface is low (around 30 %) compared to the surface of the superimposed image, the superimpositions are correct in 92 % of the pairs and 96% if we consider the pairs in the same position. In the second series of 268 pairs of images acquired a few years later, the superimposition was successful for 98 % of the pairs. Within the pairs incorrectly superimposed, the external part was presenting differences whereas the central part was perfectly superimposed. For both series, the pairs with a sufficient overlap (greater than 50 %), no noticeable difference has been perceived between them. Incorrect superimposition with noticeable differences (blurred vessel) using the method similar to [START_REF] Lee | Feature-based pairwise retinal image registration by radial distortion correction[END_REF] on the initial images (𝑘 = -0.0334). The differences might be more visible in the electronic version.

(ii) Our approach is compared to the one similar to [START_REF] Lee | Feature-based pairwise retinal image registration by radial distortion correction[END_REF] (Table 3). In the first series with 271 pairs (respectively second series with 268 pairs), the method similar to [START_REF] Lee | Feature-based pairwise retinal image registration by radial distortion correction[END_REF] finds a superimposition solution for 88% of pairs (respectively 95%) while ours gives a correct result for 92% of the pairs (respectively 98%). In figure 9, with their approach, a noticeable difference appears near a vessel located on the external part of the image, where the effect of the radial deformation is the most important, while there is no noticeable difference using our method. (iii) Our method is compared to the state of the art method "gdbicp" [START_REF] Stewart | The dual-bootstrap iterative closest point algorithm with application to retinal image registration[END_REF][START_REF] Yang | Registration of challenging image pairs: initialization, estimation, and decision[END_REF] which can make superimposition of images with a quadratic model (Table 3). In the first series with 271 pairs (respectively second series with 268 pairs), their method has found a superimposition solution for 74% of pairs (resp. 58%) while our method gives a correct result for 92% of the pairs (resp. 98%). (iv) In addition, when using with "gdbicp" a homography and two radial distortions using a multiplicative model -which is therefore not invertible -less pairs are registered (48% in the first series and 33% in the second series). The method of "gdbicp" with a quadratic model gives better results for eye fundus images than "gdbicp" with radial distortions because, "gdbicp" with a quadratic model was specifically designed to superimpose eye fundus images in [START_REF] Stewart | The dual-bootstrap iterative closest point algorithm with application to retinal image registration[END_REF].

In order to illustrate the robustness of our method, the superimposition is performed on the database of 5 patients with 20 pairs of low quality images. In each case, there was no noticeable difference in the superimposition (i) with our method, even on the border of the images. (ii) The method similar to the one of [START_REF] Lee | Feature-based pairwise retinal image registration by radial distortion correction[END_REF] gives a correct superimposition for 95% of the pairs, (iii) "gdbicp" with a quadratic model 40% and (iv) "gdbicp" with a homography and two radial distortions 25%. In this second validation, our method gives better result than the others.

Interpretation & discussion

The accuracy of the superimposition methods has been validated on a simulated montage (Table 1). The superimposition error is smaller for our model including one or two radial distortions than with a quadratic model of superimposition such as the one used in "gdbcip" quadratic [START_REF] Stewart | The dual-bootstrap iterative closest point algorithm with application to retinal image registration[END_REF][START_REF] Yang | Registration of challenging image pairs: initialization, estimation, and decision[END_REF][START_REF] Adal | A hierarchical coarse-to fine approach for fundus image registration Biomedical image registration[END_REF][START_REF] Jian | A Partial Intensity Invariant Feature Descriptor for Multimodal Retinal Image Registration[END_REF]. This confirms the results of [START_REF] Lee | Objective and expert-independent validation of retinal image registration algorithms by a projective imaging distortion model[END_REF] who previously showed that quadratic models have a greater error than a homography and a radial distortion model. Such results demonstrate that our method gives a superimposition without noticeable difference, by a human observer. Therefore, this approach is suitable to perform an analysis in a large public health database.

In terms of speed, our method is 2.4 times faster than "gdbicp" quadratic (Table 2), even using a slower computer language (Matlab versus a compiled language). The duration results show that our method is computationally efficient to analyse large databases. The good speed results are due to the linear estimators of the model parameters and the use of an invertible model of deformation. In future, we will look for an even more efficient implementation using a compiled language, which is faster.

In order to assess the efficiency for public health purposes, our superimposition method has been validated on public health databases with high quality images of 69 patients (two series including 271 pairs and 268 different pairs) and 5 patients (20 pairs) with low quality images. In the first series, the majority of the pairs were acquired with two cameras while in the second series the majority were acquired with the same camera. Whatever is the number of cameras, there is no noticeable difference in the superimposed images if the overlap is sufficient (more than 50% about) (Table 3-i).

For the pairs with a small overlapping area, we have developed another algorithm using, in addition to the matched points, the distance between the retinal vessels. This algorithm is similar to those described by [START_REF] Can | A feature-based, robust, hierarchical algorithm for registering pairs of images of the curved human retina[END_REF] and [START_REF] Lee | Objective and expert-independent validation of retinal image registration algorithms by a projective imaging distortion model[END_REF]. These findings will be presented in a future paper. However, the interest of the superimposition is to compare the evolutions in public health databases over many years, which is only useful when the image overlap is large enough. This means that our method, efficient when the image overlap is sufficient, is suitable to be applied to analyse large databases.

(v) Finally, in order to show the influence of the colour stabilisation on the superimposition, the "gdbicp" method with a quadratic model is applied on images with colour stabilisation (i.e. normalised). The results obtained with colour stabilisation are much better than without: in the first series with 271 pairs, 88% vs. 48%, in the second series with 268 pairs, 87% vs. 33% and in the third series with 20 low quality pairs of 100 % vs. 25%.

Compared to the other state-of-the-art method of [START_REF] Lee | Feature-based pairwise retinal image registration by radial distortion correction[END_REF] made of a homography and a single distortion (ii), for all the three series of images, our method (i) systematically improves the obtained results (Table 3). This result is due to the use of two radial distortions, instead of one, when images of different size are acquired by different cameras. This shows the importance of using an adapted model to the type of cameras used to acquire a pair of images. Therefore, our model goes further than the previous one [START_REF] Lee | Feature-based pairwise retinal image registration by radial distortion correction[END_REF] that was not taking into account the change of camera. We underline that in [START_REF] Lee | Feature-based pairwise retinal image registration by radial distortion correction[END_REF][START_REF] Lee | Objective and expert-independent validation of retinal image registration algorithms by a projective imaging distortion model[END_REF] the authors were interested by large mosaicking with images acquired by the same camera during the same patient examination and therefore, they were using a corresponding database (i.e. images of similar colour and same distortion parameter). Our database is made of images of the patients acquires with at least a year of interval by different cameras and they may be of low quality.

Compared to the method "gdbicp" quadratic (Yang et al., 2007) (iii), which uses a quadratic model, our method (i) better registers the eye fundus images. Our method still better superimposes the images than the method "gdbicp" quadratic with normalised images (v). However, when the image overlapping area is small (i.e. less than 50 %), "gdbicp" quadratic (v) better registers than our method. This is due to the detection of vessels landmarks, which makes the model fitting more efficient with a small overlapping. However, in this image configuration, the superimposition solution given by "gdbicp" quadratic (iii, v) often presents a strong shearing deformation that has no physical justification. This means that the quadratic model is not well adapted when the overlapping area is small. This confirms the results obtained by [START_REF] Lee | Objective and expert-independent validation of retinal image registration algorithms by a projective imaging distortion model[END_REF] who have shown that quadratic models introduces more errors than models based on a homography and a radial distortion. Indeed, with 12 parameters to be fitted (vs 8 parameters for our model -6 for the homography and 2 for the radial distortions) quadratic models have more degree of freedom that is a source of additional errors. Therefore, our model goes further than the previous approaches with quadratic models [START_REF] Can | A feature-based, robust, hierarchical algorithm for registering pairs of images of the curved human retina[END_REF][START_REF] Lee | Objective and expert-independent validation of retinal image registration algorithms by a projective imaging distortion model[END_REF][START_REF] Jian | A Partial Intensity Invariant Feature Descriptor for Multimodal Retinal Image Registration[END_REF][START_REF] Adal | A hierarchical coarse-to fine approach for fundus image registration Biomedical image registration[END_REF]. When the overlapping area is small (less than 50 %), the method "gdbicp" quadratic with normalised images (v) is complementary to our method (i) due to the detection of vessel landmarks.

To summarise the comparison, our method (i) better superimposes the images of our public health databases (Table 3). The second model is the one of [START_REF] Lee | Feature-based pairwise retinal image registration by radial distortion correction[END_REF] (ii), whose results are systematically improved by our method. The third method is "gdbicp" quadratic [START_REF] Yang | Registration of challenging image pairs: initialization, estimation, and decision[END_REF] (iii) which is improved when normalising the images (v). The method gdbicp" quadratic, and also the method described in [START_REF] Lee | Feature-based pairwise retinal image registration by radial distortion correction[END_REF], is based on the detection of vessel landmarks. Therefore, when they are used on images with a small overlapping, the vessel detection makes them more efficient and complementary to our method. We will present such improvement for our method in a future paper.

The obtained results with the validation montage and with the public health databases demonstrate that our model better corrects the errors coming from the different positions of the patient during image acquisition, the change in the camera employed (resolution and optical lenses) and the projection of a 3D scene onto a plane and the variability of colour between images.

In addition, the colour stabilisation step is useful when images have strong contrast variations. The results obtained with "gdbcip" quadratic with the colour stabilisation (v), which are better than and without (iii), show the importance of using normalised images for colour in public health databases (Table 3). In our method (i), the SIFT points are extracted and matched on the normalised images for colour contrast before estimating the registration model. Experiments have shown that using normalised images gives more robust results when extracting and matching the SIFT points, even if these points are known to be robust to colour variations [START_REF] Lowe | Distinctive Image Features from Scale-Invariant Keypoints[END_REF][START_REF] Vedaldi | VLFEat: An open and portable library of computor vision algorithms[END_REF]. Importantly, even when using normalised images with "gdbicp", our method remains better.

Moreover, our method has been tested for images acquired with a field of view (FOV) of 45° which are used for DR screening by single-field fundus photography [START_REF] Williams | Single-field fundus photography for diabetic retinopathy screening: a report by the[END_REF]. In addition, our approach could be useful for automatic detection of referral patients due to Diabetic Retinopathy [START_REF] Fleming | Automated grading for diabetic retinopathy: a large-scale audit using arbitration by clinical experts[END_REF][START_REF] Decencière | TeleOphta: Machine learning and image processing methods for teleophthalmology[END_REF][START_REF] Abramoff | Automated analysis of retinal images for detection of referable diabetic retinopathy[END_REF][START_REF] Quellec | Automatic detection of referral patients due to retinal pathologies through data mining[END_REF].

Despite the good results obtained, our method does have a limitation. For images with a smaller overlap (e.g. 30% of the surface of the mosaic image), the superimposition may present small differences on the external part. To address this issue, we have developed another algorithm using in addition to the matched points, the distance between the vessels. It will be presented in a future paper.

Conclusions

We have therefore successfully achieved a new method to answer the main challenges to superimpose eye fundus images coming from large public health databases. In addition to the previously existing methods, ours has been designed to deal with changes in terms of camera, lens, image resolution and strong colour variations between two retinal exams. The presented method consists of fitting a registration model composed of a homography and one or two radial distortions on salient points extracted in images after colour stabilisation. The choice between the number of radial distortions is made automatically depending of the number of cameras used to acquire the pair of images. Then, the image warping is performed using a division model, which is invertible, which makes it fast to compute. Linear equations are introduced to estimate the parameters of the model in a fast way, followed by non-linear estimators. Importantly, the method is easy to use and does not require to extract intrinsic characteristics such as the vessels or their branch points.

The efficiency of our method has been validated on a simulated montage with a superimposition error, which is less than the error of a quadratic model. Its efficiency has also been demonstrated on public health databases with high quality images and images of lower quality, due to differences in the conditions of acquisition, and pairs of images acquired with the same or different cameras. Nevertheless, the results show that there is no noticeable difference between the images from two examinations with the eye in the same position (nasal or macular). The superimposition is correct in more than 92% of the cases (96% with an overlapping area greater than 50 %). In these databases, our method better superimposes the images than the two state-of-the-art methods of [START_REF] Lee | Feature-based pairwise retinal image registration by radial distortion correction[END_REF] and of "gdbicp" [START_REF] Stewart | The dual-bootstrap iterative closest point algorithm with application to retinal image registration[END_REF][START_REF] Yang | Registration of challenging image pairs: initialization, estimation, and decision[END_REF].

Figure 1 .

 1 Figure 1. Superimposition of eye fundus images acquired with different cameras and image resolution. (a) Naive superimposition of images. (b) Perfect superimposition of images by using a model. For display purpose, in each subfigure (a) and (b), one of the image of the pair is in false colours (i.e. bluish colours).

  Figure 2. (a) and (b) macular eye fundus images acquired during two exams with 1 year of interval on a patient. (c) and (d) nasal eye fundus images acquired during two exams with 1 year of interval on another patient. For both patients, different cameras, resolution and lighting conditions were used.

Figure 3 .

 3 Figure3. Framework of eye fundus superimposition of images for longitudinal screening. After colour stabilisation (pre-processing) and extraction of characteristic points, the points are matched and the model is estimated. The differences are corrected in term of radial distortions and resolutions between the images.

  Figure 4. (a), (b) Eye fundus images of the same patient acquired during two exams with 1.5 years of interval with the same camera. These low quality images have been acquired in harsh conditions. (c), (d) Colour stabilisation of eye fundus images (a) and (b).

Figure 5 .

 5 Figure 5. Point extraction and matching between images. (a) Initial matching (green arrows). (b) Matching after simplification (red arrows).

Figure 7 .

 7 Figure 7. Validation of the superimposition model by registering a pair of images previously deformed. The green stars corresponds to the points of the reference image and the blue dots to the points of the current image. The arrows represent the registration errors between the two images.

Figure 8 .

 8 Figure 8. Superimposition of 3 pairs of eye fundus images with correction of one or two radial distortions. (a) Superimposition of the pair in the Figure 2 (a) and (b), with two radial distortions 𝑘 ̃1 = -0.0546, 𝑘 ̃2 = -0.0646. (b) Superimposition of the pair in the Figure 2 (c) and (d), with two radial distortions 𝑘 ̃1 = -0.0695, 𝑘 ̃2 = -0.0536. (c) Superimposition of the pair in the Figure 4 (a) and (b), with a single radial distortion 𝑘 ̃= -0.0819.

  and in Figure 8, four examples of superimpositions with the radial distortion model are shown. One can notice the good quality of the superimposition. In the Figure 8, (a) and (b), two radial distortions have been used to correct the deformation of two different cameras whereas in the Figure 8 (c) a single radial distortion corrects the deformation of a single camera. The validation results for the complete database are presented in table 3.

  Figure 9. (a) Correct superimposition with our method based on a homography and two radial distortions with normalised images (𝑘 1 = -0.0616 and 𝑘 2 = -0.0685).(b) Incorrect superimposition with noticeable differences (blurred vessel) using the method similar to[START_REF] Lee | Feature-based pairwise retinal image registration by radial distortion correction[END_REF] on the initial images (𝑘 = -0.0334). The differences might be more visible in the electronic version.

Figure 10 .

 10 Figure10. Examples of superimposition with our method. (a) Good superimposition. (b) Superimposition with small differences (in the white circle) and its enlargement (c). (d) Superimposition with small differences (in the white circle) and its enlargement (e). There is only a small overlap in the pair of (d) and (e).

  ∀𝑖, 𝑗 ∈ {1,2}, 𝑎 𝑖𝑗 ∈ ℝ, is a 2× 2 non-singular matrix (i.e. with a nonzero determinant) and represents the linear applications. 𝑇 = [ 𝑡 𝑥 𝑡 𝑦 ] , 𝑡 𝑥 , 𝑡 𝑦 ∈ ℝ, is a translation.

			𝑎 11 𝑎 12 𝑡 𝑥
			𝑎 21 𝑎 22 𝑡 𝑦	]	(1)
			0	0	1
	𝐴 = [	𝑎 11 𝑎 12 𝑎 21 𝑎 22	] ,

Table 1 .

 1 Superimposition error for different methods of superimposition with the simulated montage

	Model	Mean (pixels)	Standard Deviation (pixels)	Mean relative error respective to the image (%)	Mean relative error respective to the vessels (%)
	Homography and 1 radial
	distortion	0.86	1.75	0.03 %	2.9 %
	(normalised images) (our method)
	Homography and 2 radial
	distortions	0.92	1.94	0.03 %	3.1 %
	(normalised images) (our method)
	"gdbicp" quadratic (initial images)	2.44	1.64	0.08 %	8.1 %

Table 2 .

 2 Duration of each method of superimposition

	Model

Table 3 .

 3 Percentage of correct superimpositions of image pairs of eye fundus in public health databases for different

	methods and models

Funding

This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

Conflict of interest

No conflict of interest.