RETINAL ANGIOGRAM REGISTRATION BY ESTIMATION OF DISTRIBUTION ALGORITHM

RETINAL ANGIOGRAM REGISTRATION BY ESTIMATION OF DISTRIBUTION ALGORITHM Johann Dréo ∗ Jean-Claude Nunes ∗∗ Pierre Truchetet ∗ Patrick Siarry ∗ Université Paris XII, Val de Marne, Laboratoire Image Signaux et Systèmes Intelligents (E.A. 3956), 61 avenue du Général de Gaulle, 94010 Créteil CEDEX ∗∗ Université Rennes 1, Laboratoire Traitement du Signal et de l'Image (UMR INSERM 642), Campus de Beaulieu, Bâtiment 22 35042 Rennes CEDEX ∗

Abstract: Retinal fundus photographs are employed as standard diagnostic tools in ophthalmology. We employ optimization techniques for registration of retinal angiograms, using non-linear pre-processing (Wiener ltering and morphological gradient) and computation of a similarity criterion. The present work makes a comparison between dierent optimization techniques, namely the optical ow minimization method, the Nelder-Mead local search, the CEDA and CHEDA metaheuristics. The impact of the resolution and median ltering of gradient image is studied and the robustness of the approaches is tested through experimental studies, performed on ICG angiographies. Our proposed method has shown interesting results, especially for high resolution registration problems. Copyright c 2006 IFAC ° Keywords: Image registration, optimization, meta-heuristic, image processing, biomedical images.

1. INTRODUCTION Registration is an important tool for solving many medical image analysis problems. Many common minimization strategies have been applied to image registration problems, such as exhaustive search, gradient descent, simplex method, simulated annealing, genetic algorithms and Powell's minimization (Ritter et al., 1999; Jenkinson and Smith, 2001). In most cases, the registration is performed in two steps: image processing and optimization of a similarity criterion. The image processing step aims at improving the image quality and extracting the relevant information required to perform

the optimization step. The optimization step, in term, must nd the optimal shift according to an objective function, describing the quality of the registration. The optimization step is often carried out by implementing mathematical optimization methods, that are generally only suitable for local optimization and may fail to nd a global optimum. As the image processing and the calculation of the objective function are time consuming, global optimization methods such as metaheuristics, that require more evaluations, are often avoided in favour of local optimization methods, that need less time to nd an optimum. But, when dealing with complex problems involving local optima,

the use of global optimization algorithms can be successful (Jenkinson and Smith, 2001). Combinations of temporal images or dierent image modalities are frequently used in order to help physicians in their diagnosis. During any angiographic sequence, there will inevitably be eye movements and it is essential that this phenomenon be corrected prior to the application of quantitative analysis. The use of the registration methods has become an important tool for computer-assisted diagnosis and therapy. The registration process permits the generation of deformation elds that reect the transformation of an image in a realistic way with respect to the given anatomy. The goal of the present work was to improve analysis quality in various applications of ophthalmology by improving angiogram registration. This paper presents registration of the retinal angiograms using a metaheuristic as a global optimization method and is composed of ve more sections. It rstly examines the registration problem and the image processing methods used (Section 2) and then introduces the optimization tools (Section 3), mainly an IDEA metaheuristic. Results are presented in Section 4 and an elaborated discussion is presented in Section 5. Finally, conclusion makes up the last section.

2.1 Retinal angiography images Normally, series of ocular fundus images are obtained through uorescein and/or ICG angiography. Both of these tests are useful for evaluating the retinal and choroidal circulation and in the diagnosis and treatment approaches for many retinal diseases such as Age Related Macular Degeneration (ARMD), Cyto-Megalo-Virus Retinitis (CMVR) and Diabetic Retinopathy (DR) (Richard et al., 1998). This technique consists of an injection of uorescein or ICG in the arm's cubital vein, followed by the observation of its distribution along retinal vessels at certain time instants. Retinal angiography (up to 36 frames) (Richard et al., 1998) is usually divided into three phases, early, mid and late. It thus provides the ability to visualize vascular or choroidal structures, and possibly detect existing diseases. The vessels are the only signicant visible structures (Richard et al., 1998) in all the images used in angiographic registration. However, variations in local intensity can give rise to many diculties, namely non-uniform background of the image, poor local contrast, eye movements, various kinds of noises, and the presence of blood vessels having non-connected endings, due to local intensity variations.

2.2 Image preprocessing 2. REGISTRATION OF RETINAL ANGIOGRAMS Rigourous experimentations have shown that our earlier proposed method in (Nunes et al., 2004), based on a local search is not suciently robust (only 30% in indocyanine green angiograms) in case of strong inter-image variations. In the majority of the previously published works on automatic registration of the retinal images (Berger et al., 1999; Can and Stewart, 1999; Hampson and Pesquet, 2000; Mukhopadhyay and Chanda, 2001; Pinz et al., 1998; Ritter et al., 1999; Simo and de Ves, 2001; Zana and Klein, 1999), it was assumed that the extraction of the features or the landmarks are known a priori. However, in late phase, the retinal angiogram images are characterized by poor local contrast. The proposed method can overcome the problems associated with detection in the vascular structures. This article presents an algorithm for the automatic iconic registration of uorescein and indocyanine green angiograms which is based on non-linear processing, summing of pixel-wise squared dierences and optimization techniques.

The ltering stage or noise reduction allows a more robust registration of the images. The rst stage namely Wiener ltering attempts to obtain edge-preserving ltering of the vascular structure, to permit angiographic registration. Since during angiography there are intensity variations of the retinal structures, we computed the registration transform from gradient images (morphological gradient) and not from original images.

2.3 Image registration The registration of the retinal images is required to recognize and quantify vascular retinopathies and choroïdopathies like in DR, CMVR and ARMD. Automated registration of the retinal images enables accurate comparisons between the images and equips to automate the calculation of the changes for both the lesions and the normal anatomic structures. Mostly classical search strategies have been used in motion analysis and image registration problems (Bangham et al., 1996; Hart and Goldbaum, 1994; Irani et al., 1994; Kim et al., 2001; Odobez and Bouthemy, 1994; Zhang and Blum, 2001).

One common approach in most of these strategies is the application of a multiresolution method (Odobez and Bouthemy, 1994; Zhang and Blum, 2001) based on optical ow, where the search is performed at increasingly higher resolutions. An important area of use of medical image registration is for retinal images (Berger et al., 1999). In order to perform the reliable registration of the retinal images, feature extraction is generally employed (Berger et al., 1999). Moreover, presence of low local contrast and noise makes the detection of vascular structure dicult. However, there exists no registration method for ICG angiograms. We proposed an iconic registration method based on the optical ow (Nunes et al., 2004), but it is not robust enough in the case of ICG images. This method can be suitably applied only in the case of the uorescein images.

2.4 Method proposed for the registration 2.4.1. Registration method An overview of several registration algorithms can be found in (Brown, 1992; Maintz and Viergerver, 1998; Hill et al., 2001). We propose to compute only the translation transform since it is more signicant and more delicate to be obtained. The proposed algorithm can be summarized in ve steps: Wiener ltering of the original image, morphological gradient computing, median ltering (optional step), computing the similarity measurement (summing of absolute intensity dierences) between the two images under the current transform, • and global optimization of similarity criteria.

• • • •

Thus, before transform computing, a non-linear pre-processing is performed. Over the last few years, intensity-based (or iconic) techniques have been applied to a number of registration problems. Their basic principle is to maximize a criterion measuring the intensity similarity of corresponding pixels. The similarity measurement is the calculation used to judge the closeness between two images under the current transformation. In (Roche et al., 1999), Roche et al. demonstrated what are the assumptions made corresponding to a number of popular similarity measurements in order to better understand their use, and nally how to choose a method which is the most appropriate one for a given class of problems. The sum of squared intensity dierences (SSD), the sum of absolute intensity dierences (SAD) (Yu et al., 1989), crosscorrelation (Cideciyan et al., 1992), entropy of the

dierence image, etc. are easy to compute and often aord simple minimization techniques. Many Computer Vision algorithms employ the sum of pixel-wise squared dierences between the pair of images as a statistical measurement of similarity. (Roche et al., 1999).

3. OPTIMIZATION TOOLS An optimization problem can be dened as a problem in which the aim is to nd the best of all possible solutions. More formally, nd a solution in the search space which has the minimum (or maximum) value of the objective function. In the case of the registration problem, the objective function must describe how good is a registration between two images. We will therefore consider the similarity function and the set of the feasible image shifts respectively as the objective function and the moves within the search space. Many dierent automatic registration methods have been proposed to date, and almost all of them share a common mathematical framework, i.e. optimizing a cost function. It has been demonstrated (Nunes et al., 2004) that in case of ICG angiograms the use of optical ow together with the standard multiresolution approach is not sufcient to nd the global minimum, reliably. One can potentially employ several optimization techniques that are available in practice. One of the most widely used categories is the category of the local search algorithms like gradient descent, iterated descent, simplex methods and so on. These methods can eciently solve simple optimization problems with only one optimum. When dealing with more dicult problems, with local optima, global search algorithms are often used, especially metaheuristics, like evolutionary algorithms, simulated annealing or ant colony algorithms. The advantage of local search algorithms is their rapidity, but their main drawback is that they can get easily trapped in local optima. In this paper we have tested three methods, a simplex local search algorithm called the NelderMead Search (NMS (Nelder and Mead, 1965)), an estimation of distribution algorithm (CEDA, (Bengoetxea et al., 2002)) and an hybrid of these two algorithms (CHEDA) (Dréo, 2004), for the problem under consideration. Estimation of distribution algorithms are populationnal metaheuristics where an explicit probability distribution is used for computing the transistion between two iterations.

3.1 Objective function and parameter setting The base of the objective function used to describe the problem is the similarity criterion, previously described in the section ??. The function is dened as f : N → R to be minimized. The CEDA and CHEDA algorithms are used with default parameter values (selection ratio: 0.5). The number of points is set to 100 and the allowed maximum number of iterations for the simplex, in the CHEDA algorithm, set to 10. The stopping criterion used is a limit of evaluation number. In this paper, the upper limit is xed to 2000 calls to the objective function. The NMS algorithm uses default parameter values. The algorithm stops if the current step length is smaller than 10−8 or if it has reached the limit of evaluation number. 4. RESULTS All angiographic images used in this paper were digitized from the video signal output at a resolution of 1024 × 1024 pixels and with a quantization of 8 bits per pixel. The algorithms have been tested on dierent resolutions from 10% to 100% of the original image size, and with or without additional median ltering of the gradient image.

4.1 Preliminary tests

remains almost unchanged for dierent resolutions.

4.2 Robustness In order to test the robustness of the methods, we have performed 20 tests on the same simple registration problem. Generally classical local optimization methods used in registration of angiograms are employed for low resolution problems. As high resolution images carry more informations, we have tested the impact of resolution on the eciency of the metaheuristic. Moreover, the additional median ltering is often used to improve the optimization test. However since it is the computation time which is of prime concern for this problem, we have tested its real impact when using the metaheuristics for the optimization step. The additional ltering is only relevant for high resolutions. The standard deviation is a good measure of the robustness of the algorithm, as it is a measure of the number of bad registrations, i.e. the greater the standard deviation is, more times the algorithm fails to nd a good registration. The NMS algorithm used alone can nd optimal positions, with standard deviation increasing with resolution, whereas the EDA algorithms show very little variations in their eciencies (Figure 2). Optima distribution

0

2

4

Frequency

6

8

10

In order to comprehend the behaviour of the optimization algorithm and to test its consistency for the problem, we used a simple registration problem where the global optimum is known (Figure 1).

0e+00

1e+07

2e+07

3e+07

4e+07

5e+07

Value

Fig. 2. Distribution of the optima for the CHEDA algorithm.

Fig. 1. Simple registration problem used.

4.3 Test case

Classical mathematical optimization techniques require less computation time when performing a registration at low resolutions. However, when dealing with a high resolution problem, the time taken can drastically increase. On the contrary, the computation time used by the metaheuristic

According to the results of the robustness tests, we have tested the metaheuristic for a high resolution (100% of the original image size) without additional median ltering. We have used a set of 4 angiograms, performing 20 optimizations on each image.

The metaheuristics could not determine the precise global optimum for all the runs. Indeed, if they always nd an approaching value for the ideal registration, they can lack in determining the accurate shift values. The CHEDA metaheuristic can achieve better robustness than CEDA method, which is also better than the NMS local search used alone. However, there can still be some problems in accurate registration, where a translational motion is not sufcient to represent the real transform. In spite of having a good dominant registration, there can be a residual error at the periphery of the angiogram.

5. DISCUSSION Our results demonstrate that using a metaheuristic for the optimization step of the registration of angiograms can be very useful when dealing with high resolution images. Indeed, the main advantage of the optimization methods presented in this paper lies in their constant evaluation cost. With classical methods like the optical ow, the time cost of the optimization drastically increases with the resolution. This is not the case for the methods presented above, as the number of evaluations remains constant with variations in the resolution. This feature of metaheuristics is due to their sampling operated on the objective function, which is not dependent on the way the problem is calculated but depends on the way the algorithm behaves. However, the computation time still remains a critical constraint and our result shows that for high resolutions, employing an additional median ltering is not an advantage as it decreases the robustness and increases the computation time. Our notion is that one should try to avoid too many additional ltering steps as this may lter out some of the relevant informations presented to the metaheuristic, and thus decrease its ability to nd the global optimum. The NMS method shows the same advantage as the metaheuristic from the point of view of computation time, as its performance remains almost unchanged with increase in resolution. However, our results show that the robustness of this local search decreases when optimizing high resolution problems. This is due to the existence of local optima that can trap the local search. Indeed, this problem is often encountered in classical local optimization algorithms used in registration and is solved by a multi-resolution approach. But such approaches are highly time consuming, and the use of a metaheuristic solves the problem posed by local optima without drastically increasing the computation time.

The problem of the peripheral error is less signicant compared to the dominant translation, which is most dicult to be obtained. This problem cannot be solved through a rigid approach. Indeed, an elastic registration method can be used to correct this residual transform, as shown in (Nunes et al., 2004). 6. CONCLUSION Registration is an important issue in many applications concerning medical image analysis, as motion correction and multimodal fusion. Because of eye movements during the angiographic acquisition, a registering stage of angiographic frames is necessary to improve the quantitative analysis of the retinal pathologies. The superposition of the images will allow direct comparison with the previous images, and this can be favourably utilized to judge the progress of the disease. However, angiographic data can vary considerably in intensity (e.g. the grey levels of the blood vessels vary from dark to bright in the angiographic phases) and intensity based registration algorithms are hence unsuitable for this purpose. The optimization techniques used are especially adapted for high resolution problems where more classical techniques cannot be favourably used due to the excessive time requirement. The metaheuristics are proved to have a better robustness than a local search algorithm and can achieve good qualitative registrations. Future evolves are intended to test these algorithms on more complex problems involving more parameters, such as ane registration, elastical registration or multi-resolution methods where the low time cost characteristic of metaheuristics should be a decisive factor. Furthermore, the methods can be better adapted to the problem, with the use of a tuned initial step or a more powerful similarity criterion, like mutual information. The authors wish to pursue these research directions in near future. ACKNOWLEDGEMENTS The authors also wish to thank Centre Hospitalier Intercommunal of Créteil (CHIC) for providing a large number of retinal angiographic images, and for partially supporting this work. REFERENCES J. A. Bangham, R. Harvey, and P. D. Ling. Morphological scale-space preserving transforms in many dimensions. J. Electronic Imaging, 5:283 299, 1996.

E. Bengoetxea, T. Miquélez, P. Larrañaga, and J. A. Lozano. Estimation of Distribution Algorithms, A New Tool for Evolutionary Computation, chapter Experimental Results in Function Optimization with EDAs in Continuous Domain, pages 181194. Genetic Algorithms and Evolutionary Computation. Kluwer Academic Publishers, 2002. J. W. Berger, M. E. Leventon, N. Hata, W. Wells, and R. Kinikis. Design considerations for a computer-vision-enabled ophtalmic augmented reality environment. Lectures Notes in Computer Science, 1205:399410, 1999. L. G. Brown. A survey of image registration techniques. ACM Comput. Surveys, 24:325 376, 1992. A. Can and C. V. Stewart. Robust hierarchical algorithm for constructing a mosaic from images of the curved human retina. In IEEE Conf. on Computer Vision and Pattern Recognition, volume 22, 1999. A. V. Cideciyan, S. G. Jacobson, C. M. Kemp, R. W. Knighton, and J. H. Nagel. Registration of high resolution images of the retina. In SPIE: Medical Imaging VI: Image Processing, volume 1652, pages 310322, 1992. J. Dréo. Adaptation de la méthode des colonies de fourmis pour l'optimisation en variables continues. Application en génie biomédical. PhD thesis, Université Paris 12 Val de Marne, 13 December 2004. F. J. Hampson and J. C. Pesquet. Motion estimation in the presence of illumination variations. Signal Processing: Image Communication, 16: 373381, 2000. W. E. Hart and M. H. Goldbaum. Registering retinal images using automatically selected control point pairs. In IEEE International Conference on Image Processing, 1994. D. L. G. Hill, P. G. Batchelor, M. Holden, and D. J. Hawkes. Medical Image Registration. Physics in Medicine and Biology, 46:145, 2001. M. Irani, B. Rousso, and S. Peleg. Computing occluding and transparent motion. IJCV, 12: 516, 1994. M. Jenkinson and S. Smith. A global optimisation method for robust ane registration of brain images. Medical Image Analysis, 5:143156, 2001. M. Kim, J. C. Jeon, J. S. Kwak, M. H. Lee, and C. Ahn. Moving object segmentation in video sequences by user interaction and automatic object tracking. Image and Vision Computing, 19:245260, 2001. J. B. A. Maintz and M. A. Viergerver. A survey of medical image registration. Med. Image Anal., 2:136, 1998. S. Mukhopadhyay and B. Chanda. Fusion of 2D grayscale images using multiscale morphology. Pattern Recognition, 34:606619, 2001.

J. A. Nelder and R. Mead. A simplex method for function minimization. Computer Journal, 7: 308313, 1965. J. C. Nunes, Y. Bouaoune, E. Deléchelle, and P. Bunel. A multiscale elastic registration scheme for retinal angiograms. Computer Vision and Images Understanding, 2004. J. M. Odobez and P. Bouthemy. Robust multiresolution estimation of parametric motion models applied to complex scenes. Technical Report 788, IRISA, 1994. A. Pinz, S. Bernögger, P. Datlinger, and A. Kruger. Mapping the human retina. IEEE Trans. on Med. Imag., 17:606619, 1998. G. Richard, G. Soubrane, and L. A. Yannuzzi. Fluorescein and ICG angiography. Thieme, 1998. N. Ritter, R. Owens, J. Cooper, R. H. Eikelboom, and P. P. V. Saarloos. Registration of stereo and temporal images of the retina. IEEE Trans. On Medical Imaging, 18:404418, 1999. A. Roche, G. Malandain, N. Ayache, and S. Prima. Towards a better comprehension of similarity measures used in medical image registration. In MICCAI, 1999. A. Simo and E. de Ves. Segmentation of macular uorescein angiographies, a statistical approach. Pattern Recognition, 34:795809, 2001. J. J.-H. Yu, B.-N. Hung, and C.-L. Liou. Fast algorithm for digital retinal image alignment. In IEEE Ann. Int. Conf. Engineering Medicine Biology Society, Images Twenty-First Century, volume 2, pages 374375, 1989. F. Zana and J. C. Klein. A multimodal registration algorithm of eye fundus images using vessels detection and Hough transform. IEEE Trans. on Medical Imaging, 18:419458, 1999. Z. Zhang and R. S. Blum. A hybrid image registration technique for a digital camera image fusion application. Information fusion, 2:135 149, 2001.