Retinal Image Enhancement Using Robust Inverse Diffusion Equation and Self-Similarity Filtering

As a common ocular complication for diabetic patients, diabetic retinopathy has become an important public health problem in the world. Early diagnosis and early treatment with the help of fundus imaging technology is an effective control method. In this paper, a robust inverse diffusion equation combining a self-similarity filtering is presented to detect and evaluate diabetic retinopathy using retinal image enhancement. A flux corrected transport technique is used to control diffusion flux adaptively, which eliminates overshoots inherent in the Laplacian operation. Feature preserving denoising by the self-similarity filtering ensures a robust enhancement of noisy and blurry retinal images. Experimental results demonstrate that this algorithm can enhance important details of retinal image data effectively, affording an opportunity for better medical interpretation and subsequent processing.


Introduction
With the development of economy and the aging population, people's visual impairments has become a major public health problem all over the world. All kinds of ophthalmic diseases causing visual defects not only increase the burden of public health care system, more importantly, they also pose serious threats to social and economic activities [1,2]. Among them, as one of main blinding eye diseases, diabetic retinopathy is the most common ocular complication in diabetic patients, which includes a series of typical lesions with retinal microvascular and neuron damages caused by sugar metabolic abnormalities. It is a chronic and progressive blinding fundus disease which can be characterized by such clinical features as decreased vision, fundus bleeding and exudation, macular edema and hyperplastic lesions [3]. The overall prevalence of diabetic retinopathy was 34.6% according to a meta-analysis of 35 international renowned epidemiological studies (22896 cases of diabetes) of the world [4]. The fact of high prevalence, high blindness rate, high fashion trend, high social and economic burden, and low cognition rate makes things worse [2].
Early diagnosis and early treatment is an effective method for the control of diabetic retinopathy [5,6]. Fundus imaging by digital fundus camera is a standard diagnostic mode in ophthalmology, which captures the intensity of light reflected from the retinal surface in three different wavelength ranges [7,8]. By reason of imaging mechanism and system of fundus retina imaging itself, and the disturbance of various noise in image formation process, one often obtain noisy and blurry retinal image with nonuniform and distorted illumination, which is difficult to interpret medically and to process subsequently [7,[9][10][11]. Thus, it is indispensable to remove noise and disturbances, to improve signal-to-noise rate of image, to adjust image contrast and to enhance vessels and fine details of retinal image data [9,12,13]. By above image preprocessing, useful information in retinal image is highlighted, while useless one is weakened or removed, to make the result more suitable to clinical diagnosis and treatment [3,7,9,[14][15][16][17].
Important medical information of a retinal image lies in its retinal vessel network and local fine details. In order to accurately detect and evaluate diabetic retinopathy as soon as possible, it is very crucial to properly enhance possible retinal pathological features such as microaneurysm, bleeding and exudating spots. In this paper, a robust inverse diffusion equation is presented, which combines a powerful self-similarity filtering [27,28] for detail preserving image denoising. A flux corrected transport (FCT) technique [25,29] is used to control diffusion flux adaptively, which effectively eliminates overshoots inherent in the Laplacian operation. The proposed method extends our previous work [25] to enhance noisy images while avoiding annoying overshoots and noise magnification.
We organize this paper as follows. In Section II, related image enhancement methods by the Gamma transformation, the shock computing and the self-similarity filtering are introduced. In Section III, the robust inverse diffusion equation is built to enhance noisy and blurry retinal image data, where the flux corrected transport technique is elaborated in a subsequential process including three main steps. In Section IV, experiments on retinal images with typical diabetic retinopathy are carried out to verify the effectiveness of the proposed algorithm. Finally, conclusions and future work are included to end this paper in section V.

Related image enhancement methods (II)
Intensity transformation is the simplest technique in image enhancement through mapping a pixel value r into a pixel value z, among which the Gamma (power-law) transformation is one of basic transformation functions [18]. It is defined as where, b and γ are positive constants. The Gamma transformation with fractional values of γ can map a narrow range of input values into a wider range of output values. In a variety of devices the Gamma transformation is used to appropriately enhance image contrast and details for image capture, printing and display. In the past decades there has been an increasing research on partial differential equations (PDEs) in image enhancement [23,[30][31][32]. A great deal of successful applications of nonlinear evolving PDEs in image enhancement can mainly be attributed to their two basic characteristics: local operation and iterative processing. Osher and Rudin introduced a novel image sharpening technique, called the shock filter (SF) [30], which simulates the shock wave calculation in computational fluid mechanics: where sign is a sign function, r is a gradient operator, and u NN is the second directional derivative of image along local normal direction to isophote line. It detects an image edge using the zero-crossing of u NN , where a shock is formed at a speed of the gradient magnitude |ru|. Considering image noise in the estimation of edges, Alvarez and Mazorra added a smoothing kernel and coupled the anisotropic diffusion with the shock filter (ADSF) [31,32] for noise elimination and edge sharpening: where Δ is a Laplacian operator, G σ is a Gaussian kernel with standard deviation σ, u TT is the second directional derivative of image along local tangent direction, and c is a constant to balance the anisotropic diffusion and the shock filtering.
On the other hand, in order to effectively denoise images while preserving image details, a powerful non-local means algorithm was proposed by Buades et al. [27], which fully utilizes the big redundancy and the self-similarity of natural images in the photometric range. The discrete expression of the self-similarity filtering (SSF) algorithm is as follows. Let u be a noisy image defined in a discrete grid O & R 2 . The denoised intensity at the pixel (i, j) is expressed by where w ij (m, n) is an average weight which is determined by the similarity between the pixels (i, j) and (m, n), and is adopted as where N ij and N mn are similarity windows of size (2s + 1) × (2s + 1) centered at pixels (i, j) and (m, n), respectively. The term u(N ij ) denotes an image patch restricted in the similarity window N ij . The notation k Á k 2,a denotes a Gaussian weighted Euclidean distance between two image patches, where a is the standard deviation of the Gaussian function. The parameter h denotes a smoothing factor that controls the decay of the exponential function in the Eq (5). To reduce the computational burden and to improve the efficiency of the SSF filtering, the search window is always restricted to a proper local neighborhood (2f + 1) × (2f + 1) in O. The denominator in the Eq (4) is a normalizing factor.

Robust inverse diffusion equation (III)
As special inverse diffusion processing [23], although the shock computing [30][31][32][33] can effectively sharpen image edges and remove image noise, there are some inherent weaknesses for it to enhance retinal image. Firstly, for noisy and blurry retinal image data, it is difficult to estimate its local tangential and normal directions; and for finer details these directions are difficult to define and estimate. Secondly, in order to enhance tiny lesions important for the detection of diabetic retinopathy, where the value of image gradient is very small, it is improper for the shock computing at a speed of the gradient magnitude. Thirdly, the Gaussian and tangential smoothings in the shock computing easily erase important details when removing noise and smoothing the second directional derivative. Finally, unnatural artifacts may be produced around image edges where shocks are formed by the shock computing [23]. These defects will be compared and shown in following experiments of enhancing retinal images.
In order to overcome above difficulties, we present the following robust inverse diffusion (RID) equation: When solving numerically a nonlinear inverse diffusion equation like Eq (6) using a difference scheme, it must be discretized carefully because it is an instable process. Otherwise, numerical blowing up will appear inevitably. A strategy is to try to stop from numerical fluctuations before they appear, which is based on the Total Variation Diminishing (TVD) and nonlinear limiters [32,34]. The main idea of above flux corrected transport technique is to use a limiter function to control the change of the numerical solution by a nonlinear way, and the corresponding schemes satisfy the TVD condition and consequently eliminate above disadvantage effects.
An explicit Euler method with central difference scheme is used to approximate the Eq (6) except the gradient term. Below we detail a approach to it numerically. On the image grid O, the approximate solution is to satisfy: where l and Δt are spatial and temporal steps. Let l = 1, and d þ u k ij and d À u k ij are forward and backward difference schemes of u k ij , respectively. For example, along the x direction, ; the case is similar along the y direction. A limiter function M(p, q) is used to approximate the gradient term (see Fig 1): Here, λ is a constant to guarantee that tiny important details can also be enhanced effectively regardless of its small gradient magnitude. After the numerical discretization of Eq (6), it can be considered as an enhancement process by an iterative constrained Laplacian operation [25].
In order to improve the contrast of retinal image for the detection of tiny lesions, a Gamma transformation [18] is first used to enhance the image within proper gray levels. Then, the powerful self-similarity filtering [27] is employed to remove image noise, especially in the regions of interest (ROI). Finally, the proposed robust inverse diffusion is carried out to further sharpen important details of retinal image while avoiding overshoot artifacts. A whole flow chart of our proposed algorithm is shown in Fig 2.

Experimental results and analyses (IV)
Although retinal images can be represented in many color spaces (RGB, HSI, HSV, etc.), the selection of them highly depends on the application. In this paper, a retinal image enhancement algorithm is designed to help physicians in their task of early diagnose of retinopathy, and therefore the selected space must be as close as possible to human perception [35]. A well-established agreement is that the green channel in the RGB color space provides more blood vessel structural information and is less subject to non-uniform illumination, while the HSV color space does not preserve the fidelity of retinal images [35][36][37]. Because green light is absorbed by the blood and reflected by the retinal pigment epithelium, providing a good contrast for visualizing retinal vascular network, bleeding and exudation, we routinely extract and enhance the green channel (in the gray range of [0, 1]) from a RGB color fundus photograph [7].
In Fig 3, a retinal image with tiny microaneurysms of size 465 × 600 is enhanced for the early detection of diabetic retinopathy. Through sequential processings of the Gamma manipulation (γ = 0.6), the self-similarity filtering (h = 0.01, f = 5, s = 3) and the robust inverse diffusion (λ = 0.6, Δt = 0.15, k = 7), tiny microaneurysms and microvasculature are shown clearly due to image contrast improvement and noise removal. Moreover, our method produces fewer overshoot artifacts while avoiding noise magnification.
A further comparison is carried out with the histogram equalization [18], the ADSF filtering (c = 0.2, Δt = 0.5, n = 10) after the Gamma manipulation, and the Laplacian operation after the Gamma manipulation and the self-similarity filtering in Figs 4 and 5, where enhancement results and their zoomed local parts are shown when enhancing the macular area by these methods. The tiny microaneurysms are not too clear in the original image shown in Fig 3. The Gamma transformation improves image contrasts on the whole, but fine spots and details remain blurry without being highlighted compared with their surrounding regions. Although the histogram equalization can enhance image contrasts to some extent, it produces    nonuniform illumination distribution and noise magnification concealing some fine details. The noise magnification and artifacts (overshoots and halos) from the over-enhancing by the Laplacian operation and the ADSF filtering can be obviously observed. Moreover, the numerical blowing up will quickly come out for the multiple Laplacian operations [25], especially at bigger gradients around image edges. Only by the proposed method are tiny microaneurysms clearly shown while avoiding the artifacts and noise magnification.
In order to observe enhancement effects by these methods more clearly, local profiles (350th row, 250-300 columns) of different results are shown in Fig 6. One can see that, the Laplace operation produces overshoots and halos around two vessels. Because smaller u NN will also be enhanced indiscriminately for shock computing, the ADSF filtering over enhances image differences and leads to annoying artifacts and false edges in flat areas of image [23]. The proposed method does not produce annoying artifacts own to its proper constrained enhancement, providing a chance to early detect the diabetic retinopathy faithfully by image enhancement.
Next, in Fig 7, a retinal image with soft exudations of size 768 × 768 is enhanced to verify the proposed robust inverse diffusion for highlighting important medical features such as retinal vascular networks and local fine details. Through sequential processings of the Gamma manipulation (γ = 0.6), the self-similarity filtering (h = 0.03, f = 5, s = 3) and the robust inverse diffusion (λ = 0.5, Δt = 0.15, k = 7), one can see that, the proposed method removes noise effectively and preserves important image details. At the same time, vascular networks and exudative spots are shown more clearly while producing no artifacts.
In Fig 8, zoomed local parts of results by the proposed method are shown. Obviously, through sequential processings in three steps the degraded retinal image is greatly enhanced: image contrasts are improved, image noise is removed, and overshoots and halos are not produced, which further verify the advantages of the proposed method.
Local profiles (500th row, 510-560 columns) of different enhancement results are also shown in Fig 9. Similarly, one can see that, both the Laplace operation and the ADSF filtering produce overshoots and halos around two vessels. The proposed method does not produce annoying artifacts, ensuring a retinal image enhancement as faithful as possible.
Finally, both retinal images are enhanced by the SF filtering (Δt = 0.5, n = 10) after the Gamma manipulation in Fig 10. As discussed above, although the SF filtering enhances images by sharpening their edges, noise magnification and over-enhanced artifacts by the shock computing can be obviously observed. False over-enhanced details will make it difficult for a physician to identify abnormal pathologic changes correctly.
It is important to point out that the parameters in the proposed method will greatly affect the results of retinal image enhancement. For a specific system of fundus retina imaging, the parameters in the proposed method can be fixed by a certain amount of data simulations and tests.

Conclusions (V)
In retinal image data enhancement for early detection of diabetic retinopathy, it is crucial to highlight important pathological features such as microaneurysm, bleeding and exudating spots. In this paper, a robust inverse diffusion equation is presented by combining a powerful self-similarity filtering, where the flux corrected transport technique is used to eliminate overshoots inherent in the Laplacian operation. At the same time, the self-similarity filtering not only effectively removes image noise, but also avoids noise magnification common in image enhancement methods, resulting in a robust processing of noisy and blurry retinal image data.   Experimental results demonstrate that this algorithm can enhance important details of image data effectively without overshoots and noise magnification, affording an opportunity for better medical interpretation and subsequent processing.
For future research, we will further try to optimize the algorithm in the process of adaptive image enhancement according to the gray-level distribution of retinal lesions.