Registration of fluorescein angiography and optical coherence tomography images of curved retina via scanning laser ophthalmoscopy photographs

: Accurate and automatic registration of multimodal retinal images such as ﬂuorescein angiography (FA) and optical coherence tomography (OCT) enables utilization of supplementary information. FA is a gold standard imaging modality that depicts neurovascular structure of retina and is used for diagnosing neurovascular-related diseases such as diabetic retinopathy (DR). Unlike FA, OCT is non-invasive retinal imaging modality that provides cross-sectional data of retina. Due to diﬀerences in contrast, resolution and brightness of multimodal retinal images, the images resulted from vessel extraction of image pairs are not exactly the same. Also, prevalent feature detection, extraction and matching schemes do not result in perfect matches. In addition, the relationships between retinal image pairs are usually modeled by aﬃne transformation, which cannot generate accurate alignments due to the non-planar retina surface. In this paper, a precise registration scheme is proposed to align FA and OCT images via scanning laser ophthalmoscopy (SLO) photographs as intermediate images. For this purpose, ﬁrst a retinal vessel segmentation is applied to extract main blood vessels from the FA and SLO images. Next, a novel global registration is proposed based on the Gaussian model for curved surface of retina. For doing so, ﬁrst a global rigid transformation is applied to FA vessel-map image using a new feature-based method to align it with SLO vessel-map photograph, in a way that outlier matched features resulted from not-perfect vessel segmentation are completely eliminated. After that, the transformed image is globally registered again considering Gaussian model for curved surface of retina to improve the precision of the previous step. Eventually a local non-rigid transformation is exploited to register two images perfectly. The experimental results indicate the presented scheme is more precise compared to other registration methods.


Introduction
Nowadays, many researches have been done in medical image processing to present systems that can detect and analyze various diseases automatically. Minimizing expenditures, increasing detection speed, reducing subjective mistakes and enhancing accuracy are the chief goals of these systems. Particularly in ophthalmology, a lot of researches have been conducted to propose and apply systems for diagnosis and analysis of several retinal diseases including Diabetic Retinopathy (DR) and Age-related Macular Degeneration (AMD) [1].
Image registration is the procedure of aligning two or more images into a unique common coordinate system. Applying registration on multimodal images leads to integration of information, thus benefitting from complementary nature of images. There are various ocular imaging modalities such as Fluorescein Angiography (FA) and Optical Coherence Tomography (OCT), each of which shows some features of retina better than the others. FA is an invasive imaging technique that clearly illustrates retina vasculature through injection of fluorescein dye. FA is known as the gold standard ocular imaging for presenting neurovascular structure and assessing vasculature diseases such as DR and Diabetic Macular Edema (DME) [2]. On the other hand, OCT is non-invasive, prevalent ocular imaging that works based on the Michelson interferometry rule using near-infrared light instead of sound in ultrasound and presents a cross-sectional view of retinal structure. OCT can be used to identify macular edema, macular cysts, vitreomacular traction, sub-retinal fluid, pigment epithelial detachment, choroidal neovascularization and measuring retinal thickness [3].
As mentioned before, FA is used to diagnose DR as a common retinal disease with hyperfluorescent spots as hallmarks of Microaneurysm (MA) [4]. The greater number of MA spots indicates the DR progression. Albeit FA depicts neurovascular structure of retina exactly, it utilizes intravenous injection that is harmful for human body and requires an expert photographer with the patient remaining strictly still. Unlike FA, OCT is non-invasive with retinal in-depth data. Therefore, clinicians usually use information of both imaging modalities for better diagnosis and evaluation of retinal diseases. Great number of papers have evaluated MA visual properties via intuitive comparison of FA and OCT images simultaneously [4][5][6][7][8] whereas information of both images can be added together using automatic and precise retinal image registration. Furthermore, due to reduction in subject-dependent errors, this can be more helpful for clinical purposes.
So far, many retinal image registration techniques have been proposed which align images of various modalities into a unique coordinate system. From a technical point of view, these techniques can be classified into two categories: segmentation-based methods and non-segmentation-based ones [9]. In the former, the retinal image is partitioned into meaningful segments such as vessels and optic disc making the analysis of the image simpler [10]. In the latter, the full image contents are processed using image processing techniques. A comprehensive review, analysis, and classification of the retinal vessel segmentation algorithms and methodologies is provided by [11], also a survey of state-of-the-art techniques for extracting retinal vessels is presented in [12]. In [13], a fast and accurate registration algorithm is proposed based on Salient Feature Regions (SFR) followed by local rigid transformation. The alignment process is then augmented with a global second order polynomial transformation. This method is only applicable to mono-modal fundus images. In [14], a registration algorithm for poor quality multimodal retinal images using Harris corner detector and Partially Intensity Invariant Feature Descriptor (PIIFD) is presented. A multimodal registration of Spectral Domain OCT (SD-OCT) volumes and fundus photographs is suggested in [15] which uses Features from Accelerated Segment Test (FAST) feature detector and Histogram of Oriented Gradients (HOG) descriptor. In this method, the affine transformation is computed using RANdom SAmple Consensus (RANSAC) [16]. A multi-resolution difference of Gaussian pyramid with Saddle detector (D-Saddle) is proposed in [17], to detect feature points on the low-quality region that consists of vessels with varying contrasts and sizes. In [18] a two-step framework is proposed to register multimodal retinal images. In the first step, HOG descriptor matching is applied to the mean phase image for transformation computation and in the second step, the deformation registration method based on Modality Independent Neighborhood Descriptor (MIND) is used to improve the registration accuracy. In [19], reliable feature point sets are extracted based on SURF detector and PIIFD descriptor from model and target retinal images. Then Mixture Feature and Structure Preservation (MFSP) is used to map feature point sets for finding exact point correspondences. An Adaptive RANSAC (A-RANSAC) method is proposed in [20] for multimodal retinal image registration. In this method, the threshold value is chosen so that the Root Mean Square Error (RMSE) and the number of removed matches are optimized simultaneously. It should be noted that all of the aforementioned methods and techniques are included in none segmentation-based category. The following methods are from the segmentation-based one.
In [21] a registration of OCT projection images with color fundus photographs is proposed based on a hierarchical local feature matching that utilizes the search of local maximization of the similarity function. Reference [22] registers pairs of FA frames for the sake of segmenting fluorescein leakage in patients with DME where two-step registration method performs global and local registration in turn. In global registration step, after vessel extraction using exploratory dijkstra forest algorithm, Speeded Up Robust Features (SURF) features are detected, extracted, and matched, then M-estimator Sample Consensus (MSAC) [23] algorithm is applied for transformation computation. Thereafter in local registration step, intensity multi-resolution registration is performed on local patches to obtain optimal results. In [24], a joint vessel segmentation and deformable registration model is proposed based on Convolutional Neural Network (CNN). In vessel segmentation, a style loss guides the model to generate segmentation maps and helps transform images of different modalities into consistent representations. In deformable registration, a content loss helps finding dense correspondence for multi-modal images. An image registration algorithm is proposed in [25] to trace changes in the retina structure across modalities using vessel segmentation and automatic landmark detection. The segmentation of the vessels is done using a U-Net and the detection of the vessel junctions is achieved with Mask R-CNN. In [26] a new method is presented to register SD-OCT and FA images. For this purpose SLO images captured by Heidelberg Spectrally HRA2/OCT device are used as intermediate images. To the best of our knowledge, only [26] registers OCT images with FA photographs via corresponding SLO images.
As a result of differences in multimodal retinal images, when common feature-based registration methods are applied in RANSAC or MSAC, the outlier matches are not removed completely. State-of-the-art mismatch removal techniques are proposed in [27][28][29][30]. An efficient approach is designed in [27] based on maintaining the local neighborhood structures of the potential true matches. The problem is formulated into a mathematical model, and a closed-form solution with linearithmic time and linear space complexities is derived. In [28], the mismatch removal is converted to a two-class classification problem, and learning a general classifier to determine the correctness of an arbitrary putative match. Method in [29] adaptively clusters the putative matches into several motion consistent clusters together with an outlier/mismatch cluster. The classic Density Based Spatial Clustering method of Applications with Noise (DBSCAN) is customized in the context of feature matching to implement spatial clustering. An iterative clustering strategy is designed to promote the matching performance in case of severely degraded data. In [30], a feature guided Gaussian Mixture Model (GMM) is proposed for the non-rigid registration of retinal images. The problem is formulated as an estimation of a feature guided mixture of densities: a GMM is fitted to one point set in which the centers of the Gaussian densities characterized by spatial positions associated with local appearance descriptors are constrained to coincide with the other point set.
Due to differences in contrast, resolution and brightness of the multimodal retinal images, prevalent feature detection, extraction and matching schemes do not result in perfect matches. Also, the images resulted from vessel extraction of image pairs are not exactly the same and the accuracy of retinal image registration methods that use vessel segmentation is largely dependent on the accuracy of vessel segmentation. In addition, the relationships between retinal image pairs are usually modeled by affine transformation, which cannot generate accurate alignments due to the non-planar retina surface.
In this paper, a novel method is proposed for registering FA and SD-OCT images using SLO photographs. For doing so, after applying some preprocessing techniques on FA and SLO images, the vessels are extracted in both images and the gray level images are converted to new black and white vessel-map ones. Next, a new feature detector-descriptor method is applied on these vessel-map photographs to remove outlier matches completely. For this goal, first, SURF features are detected and extracted from both images. Likewise, FAST-HOG framework is applied on both images to detect and extract new feature set. Applying only one of these feature detector-descriptors separately in RANSAC or MSAC, does not result in perfect matches. So, SURF feature descriptors from two images are matched using approximate nearest neighbors [31]. The same procedure is iterated with HOG feature descriptors. Since this matching method may include some erroneous matches, the MSAC method is applied separately on matched SURF and HOG features to refine them, but because of vessel segmentation errors and different modalities, some outliers still remained. After that the two vectors of matched featuresF resulted from the previous steps are concatenated to enhance matched features. After applying MSAC again on concatenated features and obtaining perfect matches, global rigid transformation is calculated and applied to the vessel-map and the original images. Combining these features and applying the MSAC for three times, guarantee complete elimination of outlier matches. In addition, to consider non-planar surface of retina, the transformed images are globally aligned again applying a Gaussian model for the surface of retina which increases the accuracy of the previous step. This reduces the amount of deformation effects resulted from the subsequent local registration step. High image deformations in multimodal image registration, could affect the primary image, especially in parts of image that contain no vessels. Eventually displacement field between transformed and reference original images are estimated using method presented in [32] to register two images perfectly. Our dataset contains 36 image pairs of 21 subjects with diabetic retinopathy. Using accurate registration of FA and SD-OCT images, enables us to exactly detect location and morphological appearance of diabetic retinopathy symptoms such as MAs and leakages in OCT B-scans.
The proposed methods are explained in section 2. section 3 presents experimental results and finally, this paper is concluded in section 4.

Methods
Because of the different nature of two modalities, direct registration of OCT images and FA photographs is not possible; one presenting a depth view of retina and the other showing the surface of it. On the other hand, Heidelberg Spectralis HRA2/OCT device provides the capability of simultaneous SLO/OCT imaging by means of a single light source that guarantees exact correspondence between SLO image of retina surface and cross sectional B-scans of OCT. In this paper precise automatic registration of FA images with OCT photographs is proposed using registration of FA images with SLO photographs. Comprehensive block diagram of the proposed four-step method is demonstrated in Fig. 1. As shown, overall process is divided into four main sections: data acquisition, preprocessing, global registration, and local registration. In the following, the steps are explained in detail.

Data acquisition
In this research, dataset is comprised of 36 pairs of FA and SLO images of 21 subjects with diabetic retinopathy, where SLO image pixels are perfectly in correspondence with OCT B-scans ( Fig. 2 and Fig. 3). The FA, OCT, and SLO images are captured via Heidelberg Spectralis HRA2/OCT device in Didavaran Eye Clinic, Isfahan. Moreover the FA and SLO images are the same size as 768×768 pixels and FA images were captured with two different fields of views (30 and 55 degrees). Here, the FA and SLO images are considered as moving (source) and fixed (target) images, respectively.

Preprocessing
In this step, first the contrast of FA and SLO images are enhanced utilizing the Contrast Limited Adaptive Histogram Equalization (CLAHE) method [33]. Then FA image pixels intensities are complemented to become similar to SLO image. Next, the vessels of both moving and fixed photographs are extracted. In order to do so, the method presented in [34] is applied to precisely segment the main vessel structure of the smoothed image from original one, then binarization using Otsu's threshold selection method [35] and elimination of small objects from binary image are performed. Due to differences in contrast, resolution and brightness of the multimodal retinal images, the binary images resulted from vessel extraction of FA and SLO images are not exactly the same (Fig. 4).

Global registration
Global registration process consists of bringing FA and SLO images into spatial alignment. Here, FA image is mapped into the space of SLO photograph. For the registration scheme, it should be taken into account that the retina surface is not flat [1,36,37]. The registration schemes in [36,37] use the quadratic model for the surface of retina. Since these methods consider only the small shifts in eye and then model them as translations and rotations, they cannot be applied directly to multimodal and multi-scale images. Also, quadratic model does not adequately describe the retina surface for field of view of our image data set and leads to inaccurate results. On the other hand, studies in [1] describes the retina as a Gaussian surface, but they did not use this theory in related applications. Therefore, a new two-step global registration is proposed applying Gaussian model to improve the global registration accuracy. In our proposed method, the first step aligns two images globally using a new feature-based method and making them of the same scale, and then Gaussian model is applied to the system to calculate displacements to decrease the effects of deformations in the consequent local registration step. Rigid global registration is obtained using the following steps:

Control points detection
From another point of view, image registration techniques can be classified into two categories: intensity based techniques and feature based ones. In the former, whole image or sub-images are registered whereas the latter uses distinct image features as control points and matches corresponding features of two images for registration. Control points are well-defined image pixels that can be detected robustly and are distributed all over the image. In the proposed method, features are detected, extracted and matched to be used as input of transformation computation stage. Therefore, after vessel extraction, the control points should be detected in vessel-map images. As mentioned before, due to inevitable errors in vessel segmentation algorithms, using a feature detector-descriptor and MSAC alone does not lead to perfect matches. To address this problem, two feature detection algorithms are utilized: Features from Accelerated Segment Test (FAST) corner points [38] and Speeded Up Robust Features (SURF) [39].
FAST corner detector FAST corner detector utilizes a Bresenham circle of radius 3 containing 16 pixels around the candidate pixel ( Fig. 6(a)). If the pixel satisfies either of the two following conditions, it can be determined as the corner point: where S and I x denote a set of N contiguous pixels of Bresenham circle and the intensity of pixel x, respectively. Also I p stands for the intensity of candidate pixel and T for the threshold. In these conditions the parameters T and N are selected in a manner that there is a compromise between number of corner points and computational complexity; here N can be 9, 10, 11 and 12. Also a high speed test is used to reject non-corner pixels quickly in a way that if at least three pixels out of four specific pixels of Bresenham circle, namely pixel 1, 5, 9, and 13 are not less than I p − T or greater than I p + T, the evaluated pixel is not a corner, otherwise all 16 pixels should be examined [38] (Fig. 5(a)).
Here, L xx (x, σ) denotes the convolution of Gaussian second order derivatives with the image point x at scale σ, similar to L xy (x, σ) and L yy (x, σ). For speeding up, the approximations of second order derivatives are used instead. The point could be considered as blob where det(H approx ) of Eq. (2) becomes maximum in that point.

Feature extraction
When features have been detected, the local patch around each feature can be described as a compact representation named feature vector. In proposed method, the SURF (fast-hessian) and FAST detected features are extracted using SURF and HOG feature descriptors. SURF feature descriptor To achieve rotation-invariant descriptor, main orientation is assigned to the previously-detected SURF features. For doing so, Haar wavelet responses are calculated in vertical and horizontal directions in circular neighborhood of radius 6×s around detected feature, where s denotes the image scale of the feature. After that, appropriate Gaussian weights are applied. The main orientation is appraised via computing the sum of all wavelet responses within a sliding orientation window of 60 degrees angle. Then for a (4 × 4) sub-region aligned to the estimated orientation around detected feature, horizontal and vertical wavelet responses are taken and a vector v = ( d x , d y , |d x |, |d y |) is formed for each pixel leading to 64-dimension (4 × 4 × 4) feature vector [39].

HOG feature descriptor
Here HOG descriptor is used to describe FAST features [41]. Individual gradients may be affected by noise. This is why histogram of gradients is utilized to represent image patches that makes descriptor robust against noise. For this aim, first image is divided to blocks and each block is partitioned into cells. Thereafter gradient value and orientation are calculated for each pixel within the cell. N bins are considered, each of which represents an angle starting from 0 degree, ending by 180 − ( 180 / N ) degree with steps of 180 / N degree for unsigned gradients. Within the cell, the gradient values corresponding to the same orientation are summed leading to N-dimension vector where each dimension contains aggregate of corresponding gradient values. Taking into account that if the gradient orientation is between two bins, the corresponding value should be split proportionally into both adjacent bins. Finally, feature vector is normalized to be not affected by lighting variations. In the proposed method, N=9, cell size = [8,8] pixels, and block size = [2,2] cells, so HOG vector has 36 dimensions (9×2×2). Figure 7 visualizes HOG descriptors on vessel map images.

Feature matching
After FAST and SURF features are detected and extracted utilizing HOG and SURF descriptors, HOG descriptors of both binary images containing vessel maps should be matched. This procedure is applied for SURF features in the same way. To do this, the method in [31] is used which presents an automatic algorithm for searching approximate nearest neighbor considering data sets and desired precision to speed up the matching process. According to [31], depending on the data set and required precision, two methods namely searching hierarchical k-means trees with a priority search order [31] and multiple randomized kd-trees [42] can perform perfectly. Afterwards, since the abovementioned matching methods may include some erroneous matches, the MSAC method is applied separately on matched SURF and HOG features to refine them and obtain better matches. This results in more precise inlier set of matched features, but as can be seen in Fig. 8(a) some outliers still remained. MSAC method is a variant of RANSAC. For more details, the thresholding in RANSAC can be understood as a top-hat cost function whereas MSAC uses truncated quadratic cost function [43]. The steps of MSAC are as follows: a. Select a subset of two pairs from matched features (may contain incorrect matches) and compute the similarity transformation from this subset.
b. Apply the transformation to the remaining matched features from moving image and calculate the distances between transformed feature points and corresponding feature points from fixed image. The set of transformed features of moving image with their distances from corresponding features of fixed image are less than the previously-mentioned threshold are considered as inliers and denoted by S i .
c. Repeat the steps a and b for N times (N=1000) and determine the largest number of members for S i as inliers. In this step, transformation computation is used only to refine inlier matches and overall transformation estimation is done in the next step.

Transformation computation
To compute the similarity transformation which is registration objective, the matched FAST and SURF features of fixed and moving images are joined. Thereafter to calculate transformation, the MSAC method is utilized again and final transformation is re-computed utilizing all inliers associated with the largest number of S i . Using SURF and FAST features together in this manner has two main advantages: first, it benefits from properties of both feature detection-description frameworks and secondly, refining matched features using the MSAC for three times, causes the erroneous matches to be removed completely. Figure 8 shows how final MSAC step refines the previously-matched features and perfects them. Then the computed transformation is applied to the original images. The proposed feature detection, description and matching scheme is flexible and can be generalized. In fact, other feature detection, extraction and matching techniques can be added to the scheme to increase the number of inlier matches under the penalty of longer execution time. In general, for k number of separate feature detectors, m-sac should be applied for k + 1 times.

Applying the Gaussian model for retina surface
As mentioned before, the retina surface is not flat and in the proposed method, Gaussian model is applied to improve the global registration accuracy. Figure 9 illustrates a point on retina surface that is expressed in two different 3D coordinate systems related to two different cameras namely C p and C q . These points are represented as P T = (X, Y, Z) and Q T = (X , Y , Z ) in coordinate systems of C p and C q . Assume the transformed FA and SLO images are shown by I p and I q corresponding to P and Q coordinate systems. The Gaussian equation which models the retina surface and relates X, Y and Z in the coordinate system of C p is: where A 1 , A 2 , A 3 , . . . , A 5 are unknown parameters. Here, the aim is to derive the transformation that maps the image coordinates p to q, using the transformation that relates P and Q. Because of the global registration of previous step, the rigid transformation relates P and Q (Eq. (4)): where Q, P, R and t are C q and C p coordinate systems, orthonormal rotation and translation matrices. Here, weak perspective projection model adequately projects the 3D retina surface to 2D image plane [36,37]. The projections of P and Q in images I p and I q considering homogenous coordinates is as follows: where M p and M q are weak perspective camera projection matrices: In Eq. (7), α x , α y , α x , and α y are pixel dimension parameters, c x , c y , c x , and c y denote the centers of projection in C p and C q , while s and s are scaling parameters of the weak perspective projection. From Eq. (6): Applying Eq. (8) to Eq. (3) results in: where a 1 = A 1 , a 2 = c x +(α x ×A 2 ) s and so on. Now from Eq. (5), Eq. (8), and Eq. (9), the X and Y , are calculated in terms of x and y. Finally x and y can be obtained from Eq. (8) as follows: The unknown parameters d 11 , d 12 , . . . , d 18 , d 21 , . . . , d 28 are calculated using robust least squares method. For doing so, the first step of global registration is applied again to the transformed images and the matched features are utilized as the points to find parameters d 11 through d 28 .
As shown in Fig. 10 (also in Fig. 13 and Table 1) the proposed Gaussian model significantly improves the result of first step compared to quadratic model and decreases the deformations effects of the consequent local registration step.

Local registration
Since global registration aligns images approximately ( Fig. 11(a), Fig. 11(b), Fig. 11(d) and Fig. 11(e)), local registration is required to estimate displacement field and register images accurately (Fig. 11(c) and Fig. 11(f)). For local registration, the method proposed in [32] is applied. To achieve this, one image called diffusing model should be locally deformed to become similar to the other fixed one, named scene. In this approach, diffusing model (transformed FA image) is assumed as a grid model that must be deformed. This model contains vertices that are particles. On the other hand, contours of the scene objects (SLO image objects) are assumed to be membrane with effectors named demons. These demons are placed throughout the membrane and the presumptive vectors which are perpendicular to this contour, oriented from inside of image objects to the outside. In grid model, vertices polarity are labeled as either inside or outside and accordingly in the other image, demons locally apply forces to push the model inside the other image objects or vice versa [44]. For example, the demon corresponding to vertex with inside polarity, pushes that vertex of model toward inside of the image object. This method applies transformation iteratively. Magnitudes of forces applied from demons are constant within an iteration, but decreases between iterations for convergence. On the other words, local registration using diffusing model can be divided into two steps:

Pre-computation of demons set
The demons should be extracted from the scene (SLO). They can be elicited from all image pixels, image contours or edges. Each demon includes information of its spatial position, direction from the inside to the outside, current displacement from the corresponding point of the model.

Iterative process
After applying an initial transformation to the model, the iterative process is started. Each iteration is composed of two steps: a. Computing associated elementary force for each demon which depends on the demon direction and the polarity of grid model at the point that corresponds to the demon.
b. Computing next transformation using the previous one and all elementary demon forces computed in step a.
In the proposed method, to increase robustness and speed of local registration process, the number of iterations in each multi-resolution pyramid level equals 100 and three pyramid levels is considered. Figure 12 indicates data set images after local registration which are zoomed around Region of Interest (ROI).

Results
In this paper, three standard and prevalent parameters have been considered to assess the performance of proposed registration method: Root-Mean-Square (RMS) error, success rate and speed of registration.

RMS error
In mono-modal image registration, the pixel intensity values of two images can be utilized to determine distance between corresponding pixels. However in this case, for multimodal images, distances are calculated according to Eq. (11) where n and ||p i || stand for the number of landmarks and Euclidean distances of i-th corresponding landmarks.
Since vessel branches and crossovers are easy to locate even in multimodal retinal images, these points have been determined as landmarks in each image. Here n equals 15. The mean and standard deviation of RMS error of proposed method after global and local registration steps is indicated in Table 1, also the box-plot of RMS errors of global and local registration steps are demonstrated in Fig. 13. The process of measuring RMS errors is done by two observers, and the differences between mean and standard deviation of RMS errors measured by them is negligible. Nevertheless, the maximum values are reported in Table 1. Table 2 compares mean RMS error of proposed method with other methods. It can be seen from this table that due to combining SURF and FAST features, applying Gaussian model for retina surface, and final local registration, in the proposed method the RMS error have been dramatically reduced.

Success rate
Success rate determines what percentage of images are aligned correctly. Registration success or failure, is directly observable. The success rate of different methods is compared in Table 3. Figure 14 indicates qualitative evaluation of the proposed method.

Speed
Running time of the proposed method is compared with other methods in Table 4. It is worth mentioning that, our experiments were done using a PC with Windows 7 64-bit operating system, 4 GB RAM, and Intel Core i5 CPU 2.5 GHz whereas other reported results were obtained utilizing a PC with Windows 7 64-bit operating system, 64 GB RAM, and Intel Xeon i5 CPU 3.7 GHz.

Discussions and conclusions
In this paper, a novel method is proposed to register OCT images with FA photographs in an accurate manner. This method uses SLO images as intermediate, and includes four main steps namely image acquisition, preprocessing, global registration, and local registration. Because of using SURF and FAST features together and refining the matches for three times, the proposed method has inherited the advantages of both methods resulting in enhanced success rate. On the other hand, applying Gaussian model for retina surface and local registration after global registration leads to dramatic reduction in RMS error. Precise registration of FA and OCT images via SLO photographs enables us to exactly detect location and morphological appearance of diabetic retinopathy symptoms such as MAs and leakages in OCT B-scans. In this research, after registration of FA and OCT images using proposed method, morphological appearance and location of MAs in the OCT b-scans are estimated exploiting the same properties in the corresponding FA images. Using the fact that MAs in the FA images are indicated by hyperfluorescent dots, and utilizing the position of these dots in FA photographs, we analyzed location and morphological appearance of MAs in OCT b-scans (Fig. 15). For this purpose, 38 MAs which are located in 20 images are analyzed. In 8 out of 38 cases, MAs had complete or incomplete ring sign, whereas remaining MAs have no ring sign. Some leakages in OCT B-scans were observed as blob signs whereas other leakages did not show any visual sign. Figure 16 illustrates blob sign of a leakage in FA and corresponding B-scan of OCT image. Also, many researches have been conducted to analyze appearance of MAs in OCT B-scans [5,46].

Funding
Office of Vice Chancellor for Research and Technology, University of Isfahan (1390734).