INTRODUCTION
In the clinical practice, fundus examination is obligatory in diagnosing both ophthalmic diseases and circulation system diseases. Automatic macula localization plays an important role for ophthalmologists in facilitating the disease diagnosis and treatment. The computer-aided macula localization can ease ophthalmologists’ workload and improve the accuracy and efficiency.
Most of the existing macula localization methods are based on morphological features. The macula is approximately circular in shape, and dark red or reddish-brown in color. In the terms of the intensity, the macular areas are the darkest in fundus images. Accordingly, Lu and Lim[1]construct a linear operator based on the brightness and shape of the macula to localize the macula; Singh and Sivaswamy[2]propose to enhance the contrast of macula with the surrounding area in order to improve accuracy, but when macular morphological features are unobvious due to the lesion, the resulting localization becomes poor. An obvious characteristics of macula is that it hasn’t definite and clear boundary as optic disc (OD) and blood vessels (BV), which makes its localization problem more difficult than OD or BV.
Other morphological features related to OD and/or vascular structure are also used for macula localization, which are adopted in the methods[3-7]. These methods basically follow two steps, OD recognition and macula localization. Some methods are also combined with the macular morphological features, in order to improve accuracy[8]. Kamble et al[9]present an approach for fast and accurate localization of OD and fovea using one-dimensional scanned intensity profile analysis. Johny and Thomas[10]present the computer based approach for the diabetic macular edema detection. Issac et al[11]present a computer vision algorithm for automated and efficient localization of macula from low contrast and diabetic retinopathy affected fundus images. A statistical based model is used to detect macula in a specified region of fundus image which is designed using the geometric features of OD.Ramasubramanian and Selvaperumal[12]present a stand-alone application for the segmentation of OD and macula.
Macula localization using vascular structure information is based on the two facts, no BV in macular area, and the macular fovea located around the center of the maximum circlefitting into the vessel free area[13]. The macular area was preliminarily localized by seeking the darkest circular area, and further refined by selecting the vessel free area among the candidates.A vessel origin (VO) based localization scheme, with VO as the vertex of the parabola-like vasculature, is proposed to get fovea localization[14]. Deka et al[15]present an approach for detection of macula and fovea by investigating the structure of BV in the macular region localization of macula.
Some machine learning methods are also used for macula localization. Sedai et al[16]propose a two-stage deep learning framework for accurate segmentation of the fovea in retinal color fundus images. Tan et al[17]develop and train a convolutional neural network to automatically and simultaneously segment OD, fovea and BV.
The research of this paper attempts to present a novel method of macula localization with better accuracy and robustness,independent of morphological features. The core idea of our work is to identify the macular structure by using structure label transfer. Although there exist the individual differences for each fundus image, the physiological consistency guarantees that the individual fundus images are comparable to each other. Based on the fact, the location of the macula can be determined by compare a processed image with a reference image with a known macular label. When the similarity of two images is satisfied, the macular label can be transferred from the reference image to the processed one as a result of macula localization.
Our innovation point lies in structure label transfer, which is a process to match the processed image with the reference images with known structures and to determine the structures in the processed image by transferring the known structure label in the most similar reference image to the processed image. The transferred label represents a structure for the processed image, and macula localization in the processed image can be achieved in this way.
The rest of the paper is organized as follows: Section 2 presents our macula localization method, and depicts steps in detail. Section 3 compares the results of the proposed method with the existing methods in the experiments and discusses the performance of our method. Section 5 summarizes this paper.
MATERIALS AND METHODS
Our method uses the idea of structure label transfer, which can be depicted as a two-step process. Thefirst step is tofind the most similar images to the processed image from the reference images. We use the global image descriptor (GIST)to calculate the similarity between the processed image and the reference images, and select the most similar reference images as the candidate images. The second step is to determine the location of macular area. We use the local feature descriptor Partial Intensity Invariant Feature Descriptor (PIIFD) to form the correspondence between the processed image and the candidate images. Through feature matching, one can achieve structure label transfer from the candidate images to the processed image. As a result, the location of macular area in the processed image is determined according to the location of macular area in the reference images and their correspondence.The flowchart of the proposed method is presented in Figure 1 and will be described in the following subsections in detail.
Figure 1 The flowchart of the proposed method.
Structure Labelling Structure label transfer is premised on the basis of structure labelling. Structure labelling is used to build a set of reference maps as a priori knowledge by adding the labels by experts or ophthalmologists to the corresponding fundus images. The process is demonstrated in Figure 2, where yellow label corresponds to the center of a macular structure.Structure labelling is just applied to typical healthy fundus images. By structure labelling, the labelled maps corresponding to fundus images are constructed as a reference set for the subsequent macula localization. It should be emphasized that label maps are prepared in advance, and the number of label maps will in fluence the accuracy of macula localization.
Mask Extraction A mask is a binary image with white corresponding to the retinal area (view area) and black corresponding to the rest part. The retinal area is the only area we are concerned with, and the black area doesn't contain any useful information.
We use a simple and convenient mask extraction method proposed by Haar[18], which binarizes red channel of color fundus images with the threshold 0.137, followed by opening,closing, erosion operation with a 3×3 structural element and histogram equalization on the binary images.
Candidate Image Selection We use global image processing to select the most similar images to the processed image as the candidate images among all reference images. The candidate images will be used to establish an appropriate mapping to the processed image in order to transfer the macular label from the candidate images to the processed one. Because of the diversity of processed images, reference images must be adequate enough to meet a potential similarity matching between a processed image and a reference image. The more reference images are provided, the more possibly an arbitrary processed images can be matched by a similar reference image.
Figure 2 The illustration of structure labelling of the fundus image A: An original fundus image; B: The corresponding mask; C: The structure label map corresponding to the macula in A.
Figure 3 The demonstration of selection of candidate images A-F: Six most similar candidate images.
Global Image Descriptor Feature Extraction Taking into account computation complexity, feature similarity is more efficient than the pixel similarity for similarity matching. We adopt the GIST feature descriptor as a matching object to measure the similarity between the processed image and the reference images.The GIST is a general and abstract feature descriptor that describes globally the image with a set of indexes including degree of naturalness, openness, roughness, expansion, and ruggedness[19]. Its main applications lie in image classification.The GIST can be obtained by applying Gabor filter. Gaborfilter proposed by Gabor[20]is very sensitive to the image edge and texture, and has excellent characteristics of directionality and scalability. In practice, a multi-direction and multi-scale Gabor filter can be adaptive to edge features of arbitrary objects by adjusting its direction and scale parameters, in order to achieve the desired processing effect. The m×n Gaborfilters can be constructed by setting m direction parameters and n scale parameters. Assuming that the size of the processed image f (x, y) is w×h, the GIST of f (x,y) can be calculated as follows
where * is a convolution operator, cat is a concatenation operator, and the operation results is a m×n×w×h-dimensional vector.
In order to get more accurate feature representation for an image, a reasonable strategy is the block division of the processed image. The Gaborfilters are applied to each block respectively to produce a blockwise GIST feature. The GIST features of all blocks are finally merged to generate a set of block wise GIST features for a fundus image.
It should be noted that GIST features can be calculated in advance for the reference images. When the similarity of a processed image and a reference image is measured, it is only required to compute the GIST feature of the processed image, and then obtain the GIST feature distance between the processed image and the reference image. The GIST feature distance will serve as a measurement of their similarity. The less their GIST feature distance is, the more similar they are.
Global Image Matching Given a processed image, the similarity in terms of the GIST feature distance varies depending on the various reference images. Surely, not all reference images resemble the processed image. The reference images with the least feature distance is taken as the most similar images to the processed image. The most similar k reference images are selected as candidate images for the subsequent processing, and the other reference images would be omitted as dissimilar images, where the positive integer k can be set empirically.
Figure 3 illustrates 6 candidate images from 89 reference images in the DIARETDB1 dataset with regard to a processed image. The first one in the figure is a processed image, and the others are candidate images ranking according to the descending similarity (or the ascending feature distance).
Local Feature Matching The blockwise GIST feature descriptor guarantees that the similarity of selected candidate images to the processed image can retain the level of blocks.However, it is not enough to satisfy the requirement of label transfer. The structure label transfer requires that the matching should reach the level of features. Therefore, in order to transfer structure labels in candidate images to the processed image,finer feature matching is indispensable and necessary.Here, we will use the Harris-PIIFD descriptor[21]to perform feature matching.
Harris-PIIFD Descriptor PIIFD is a local feature descriptor invariant to image rotation, partially invariant to image intensity, affine transformation, and viewpoint/perspective change[21]. Harris-PIIFD feature, which serves fundus feature matching, uses feature corners instead of blood vessel branch points as control points, because there are enough corners in an image to be extracted, and the distribution of corners is uniform. PIIFDs are stable to the image rotation, and the PIIFD extraction is related to the main orientation of the control point.A large number of Harris feature corners are evenly distributed in a fundus image. Therefore, Harris detector is used to extract feature corners as control points. Harris detector uses Gaussian window for convolution to calculate the intensity changes in all directions and to represent them in the gradient. The feature corners aren’t directly used as feature points for feature matching, but rather provide candidate positions to calculate PIIFDs.
Before PIIFD extraction, the main orientation relative to the gradient is calculated and assigned to each candidate control point. The average rectangular gradient method is used to extract the main orientation. For each control point, its main orientation is expressed as
where
where hσis a Gaussian kernel, * represents convolution,
In order to extract Harris-PIIFDs, we sample the local neighborhood to determine the image gradient size and direction. A neighborhood of 16×16 pixel size centered at the current point is evenly divided into four 4×4 pixel regions. In each 4×4 pixel region of a neighborhood, one can calculate the gradient of each pixel, and accumulate the gradient direction in a uniform 16-bin direction histogram.
The image gradient is used to achieve partial intensity invariance. For this, the gradient intensity of the image needs to be normalized to reduce the influence of gradient change.By summing gradients in 2 opposite directions, the previous 16 directions are reduced to 8 directions. The extracted 128-dimensional feature vector is obtained as a Harris-PIIFD descriptor, which is used to calculate the similarity for feature matching.
Feature Matching After a set of feature points with feature descriptors are extracted, feature matching can be conducted,which can be divided into two steps, matching point selection and matching point filtering. The former is used to establish the mapping between the matching point pairs, and the latter is used to refine the selection and eliminate the wrong matching.In thefirst step, we adopt the nearest distance ratio method to speed up the matching. The nearest distance ratio method is a commonly used method of feature point matching. Its basic strategy is tofind the top two most similar points qj1, qj2in the candidate image for the feature point piin the processed image according to the Euclidean distance of the feature vectors. The Euclidean distance between piand qjkk=1, 2 can be expressed as follows,
where D (i, j1) and D (i, j2) represent the 1stand 2ndnearest distance, i.e. D (i, j1)≤D (i, j2), and ||·||2is an Euclidean norm operator. Here we define a ratio R=D (i, j1)/D (i, j2). Through different experiments, the correspondence between piand qjkk=1, 2 is considered as a unilateral match when their ratio is less than the threshold of τ=0.7. The above process is tofind the corresponding point pairs in the candidate images.
The process of matching point selection could yield multiple matching point pairs, and lead to one-to-many situation.This means that more than one point in the reference image correspond to a point in the processed image. Therefore,the resulting matching point pair needs to be refined, where random sample consensus (RANSAC) will be used in the second step. RANSAC is an iterative method to estimate parameters of a mathematical model from a set of observed data that contains outliers. The points within the estimated range are considered as “inliers”, otherwise “outliers”. After the removal of outliers, the parameters of the model will be re-estimated iteratively. After N iterations, the maximum point set as a support point set is determined. Then the stable transformation matrix H can be obtained as the final matrix according to the support point set.
Structure Label Transfer In the stage of structure label transfer, a weighted sum method is used to yield the structure label of the processed image. Its basic idea is that multiple candidate images can partially contribute their own label to the processed image according to their similarity. Assuming that f (x,y) is the resulting structure label at (x,y), n candidate images are selected from the reference images, and miis a structure label map corresponding to the ithtransformed candidate image, f (x,y) is given as follows,
Where ωiis the weight of the structure label transferring to the processed image, which can be obtained by calculating the feature distance of the PIIFD between the processed image and the ithcandidate image, defined as follows,
Figure 4 Experimental results of macula localization of fundus images.
where c is a normalized factor, P and Q are two corresponding feature point sets in the processed and the transformed candidate image, respectively. Obviously, the less the distance is, the more similar two images are, and the greater their weight is. Using Eq. (6), the location of macula in the processed image can be computed.
RESULTS
In the experiments, four datasets are used as test data, DRIVE,DIARETDB0, DIARETDB1 and STARE. The original images in the four datasets aren’t attached by macular labels.Therefore, the visible macular areas in reference images are manually labelled by experts and ophthalmologists, where“visible” refers to the state of macular area without occlusion or interference of serious lesions. These fundus images accompanied by corresponding macular labels in structure label maps are used as reference data.
In addition, the processed image is divided into blocks, and 32 Gaborfilters constructed by setting the direction and scale parameters to be 8 and 4 respectively are applied to each block to obtain the GIST features of the image.
Localization Performance We apply the proposed method to four datasets, and Figure 4 shows 6 examples of experimental results, including 1 healthy (the first) image and 5 lesion images (the others).
The columns from left to right are original images, preprocessed results, candidate images (3rdand 6thcolumns) and final localization results superimposed on the original images.The green and blue crosses represent the estimated and true macular location.
Performance Evaluation We use 3 quantitative indicators in order to evaluate and analyze the performance of the proposed method, including success rate, accuracy and robustness. We will compare the proposed method with the existing methods on the performance indicators.
1) Success rate: The success rate of macula localization is a percentage of successful macula localization among total test images. For convenience, one uses the OPTIC DISC DIAMETER (ODD) as a unit to measure the distance in the fundus image. In our experiments, the situation where the distance between the estimated and the true macular location is less than an ODD is taken as successful macula localization.
2) Accuracy: The error distance between the estimated and the true macular center is used to indicate accuracy. The smaller the distance is, the more accurate the localization results are.
3) Robustness: Robustness is an indicator used to describe the success rate contrast when dealing with lesion and healthy images, respectively, expressed as a ratio of success rate difference to the lowest success rate.
Comparison and Discussion of Performance Indicators We assess the success rate, accuracy and robustness of the macula localization in four datasets. The results are shown in Table 1 and Table 2.
For the sake of clarity, the accuracy of the proposed method is compared with the other methods. The accuracy of the macula localization is divided into three levels, and the proportion of each accuracy level of the localization results obtained using three different methods is listed in Table 3.
By analyzing Tables 1, 2 and 3, the performance of the proposed method can be summarized as follows. The average success rate in four datasets is 98.18%, where 100% success rate is achieved for DRIVE and DIARETDB1. For accuracy,the average error distance in four datasets is 0.151 ODD, and Table 3 shows that 88% of error distance is less than 0.35 ODD, which is obviously better than other methods. In four datasets, 49% of error distance is less than 0.1 ODD, 76% less than 0.2 ODD. Therefore, the high accuracy is a significant advantage of the proposed method. In terms of robustness,Table 2 shows that the proposed method has a slight difference in performance contrast when dealing with the healthy and the lesion images in four datasets. Even for severe lesion images, the proposed method can still maintain high success rate and high accuracy, e.g. 95.65% and 0.124 ODD in the case of STARE dataset, respectively, which indicates that the proposed method is highly robust and stable in the complicated situations.
In summary, the experimental results show that the proposed method has a high success rate, accuracy and robustness.
DISCUSSION
The research in this paper attempts to improve the accuracy and robustness of macula localization. We propose a method to identify the macular structure by using structure label transfer.Our method is outlined as the following steps, seeking the most similar images to the processed image in the dataset as the candidates, matching locally the processed image with the candidate images accompanied by the known structure labels at the level of features, transferring the structure label representing macular in the candidate image to the processed image, and obtaining the structure label in the processed image.
We perform the evaluation of localization results, respectively,and compare the proposed method with other existing methods.Through the comparison, we can see that the proposed method can avoid the interference of lesion to macular morphological features in macula localization, and outperforms the other methods in success rate, accuracy and robustness.
Table 1 Two performance indicators of macula localization applying the proposed method
Table 2 The robustness of macula localization applying the proposed method
Table 3 Comparison of the accuracy of the proposed method with other macula localization methods
ACKNOWLEDGEMENTS
Foundations:Supported by the National Key Research and Development Program of China (No.2016YFB0201503,No.2017YFC0602203); the 13th Five-Year Plan of the Science and Technology Research of the Education Department of Jilin Province (No.2016433); the National Natural Science Foundation of China (No.60905022); the PhD.Program Foundation of the Ministry of Education of China(No.20130061110054).
Conflicts of Interest:Guo XX, None; Li Q, None; Sun C,None; Lu YN, None.
REFERENCES
1 Lu S, Lim JH. Automatic macula detection from retinal images by a line operator. 2010 17th IEEE International Conference on Image Processing 2010:4073-4076.
2 Singh J, Sivaswamy J. Fundus foveal localization based on image relative subtraction-IReS approach, in: Proceedings of the 14th National Conference on Communications 2008.
3 Gagnon L, Lalonde M, Beaulieu M, Boucher MC. Procedure to detect anatomical structures in optical fundus images. Medical Imaging 2001:Image Processing. 2001.
4 Narasimha-Iyer H, Can A, Roysam B, Tanenbaum HL, Majerovics A.Integrated analysis of vascular and nonvascular changes from color retinal fundus image sequences. IEEE Trans Biomed Eng 2007;54(8):1436-1445.
5 Simandjuntak RA, Suksmono AB, Mengko TLR, Sovani I.Development of computer-aided diagnosis system for early diabetic retinopathy based on micro aneurysms detection from retinal images,Proceedings of 7th International Workshop on Enterprise networking and Computing in Healthcare Industry, 2005 on Enterprise networking and Computing in Healthcare Industry, 2005. HEALTHCOM 2005.
6 Medhi JP, Dandapat S. An effective fovea detection and automatic assessment of diabetic maculopathy in color fundus images. Comput Biol Med 2016;74:30-44.
7 Kao EF, Lin PC, Chou MC, Jaw TS, Liu GC. Automated detection of fovea in fundus images based on vessel-free zone and adaptive Gaussian template. Comput Methods Programs Biomed 2014;117(2):92-103.
8 Sinthanayothin C, Boyce JF, Cook HL, Williamson TH. Automated localisation of the optic disc, fovea, and retinal blood vessels from digital colour fundus images. Br J Ophthalmol 1999;83(8):902-910.
9 Kamble R, Kokare M, Deshmukh G, Hussin FA, Mériaudeau F.Localization of optic disc and fovea in retinal images using intensity based line scanning analysis. Comput Biol Med 2017;87:382-396.
10 Johny A, Thomas A. A novel approach for detection of diabetic macular edema. 2016 International Conference on Emerging Trends in Engineering, Technology and Science (ICETETS), 2016.
11 Issac A, Sengar N, Singh A, Dutta MK, Prinosil J, Riha K. An efficient imaging technique for automated macula localization from fundus images.2016 8th International Congress on Ultra Modern Telecommunications and Control Systems and Workshops (ICUMT), 2016.
12 Ramasubramanian B, Selvaperumal S. A stand-alone MATLAB application for the detection of optic disc and macula. 2016 International Conference on Advanced Communication Control and Computing Technologies (ICACCCT) 2016.
13 Li H, Chutatape O. Automated feature extraction in color retinal images by a model based approach. IEEE Trans Biomed Eng 2004;51(2):246-254.
14 Yu CY, Liu CC, Yu SS. A novel scheme for the fovea localization on retinal images. in: Computer, Consumer and Control (IS3C), 2014 International Symposium on, IEEE 2014;609-612.
15 Deka D, Medhi JP, Nirmala SR. Detection of macula and fovea for disease analysis in color fundus images. Recent Trends in Information Systems (ReTIS), 2015 IEEE 2nd International Conference on, IEEE,2015:231-236.
16 Sedai S, Tennakoon R, Roy P, Cao K, Garnavi R. Multi-stage segmentation of the fovea in retinal fundus images using fully convolutional neural networks. Biomedical Imaging (ISBI 2017), 2017 IEEE 14th International Symposium on, IEEE 2017:1083-1086.
17 Tan JH, Acharya UR, Bhandary SV, Chua KC, Sivaprasad S.Segmentation of optic disc, fovea and retinal vasculature using a single convolutional neural network. J Comput Sci 2017;20:70-79.
18 Frank ter Haar. Automatic localization of the optic disc in digital colour images of the human retina. Computer Science, Utrecht University(2005).
19 Oliva A, Torralba A. Modeling the shape of the scene: a holistic representation of the spatial envelope. Int J Comput Vision 2001;42(3):145-175.
20 Gabor D. Theory of communication. part 1: The analysis of information, Journal of the Institution of Electrical Engineers-Part III:Radio and Communication Engineering 1946;93(26):429-441.
21 Chen J, Tian J, Lee N, Zheng J, Smith RT, Laine AF. A partial intensity invariant feature descriptor for multimodal retinal image registration.IEEE Trans Biomed Eng 2010;57(7):1707-1718.
22 Liu Z, Opas S, Krishnan SM. Automatic image analysis of fundus photograph. Engineering in Medicine and Biology Society, 1997.Proceedings of the 19th Annual International Conference of the IEEE 1997;(2);524-525.
23 Sagar AV, Balasubramanian S, Chandrasekaran V. Automatic detection of anatomical structures in digital fundus retinal images. in: MVA,2007;483-486.