Retrospective Illumination Correction of Retinal Images

A method for correction of nonhomogenous illumination based on optimization of parameters of B-spline shading model with respect to Shannon's entropy is presented. The evaluation of Shannon's entropy is based on Parzen windowing method (Mangin, 2000) with the spline-based shading model. This allows us to express the derivatives of the entropy criterion analytically, which enables efficient use of gradient-based optimization algorithms. Seven different gradient- and nongradient-based optimization algorithms were initially tested on a set of 40 simulated retinal images, generated by a model of the respective image acquisition system. Among the tested optimizers, the gradient-based optimizer with varying step has shown to have the fastest convergence while providing the best precision. The final algorithm proved to be able of suppressing approximately 70% of the artificially introduced non-homogenous illumination. To assess the practical utility of the method, it was qualitatively tested on a set of 336 real retinal images; it proved the ability of eliminating the illumination inhomogeneity substantially in most of cases. The application field of this method is especially in preprocessing of retinal images, as preparation for reliable segmentation or registration.


Introduction
Improper scene illumination as well as nonideal acquisition conditions due to for example, misadjusted imaging system can introduce severe distortions into the resulting image. These distortions are usually perceived as smooth intensity variations across the image. According to the terminology commonly used in processing of magnetic resonance (MR) images, we call these systematic intensity level inhomogeneities as the bias field. With such unevenness, the subsequent image processing like image registration, segmentation, or pattern recognition may be substantially complicated; therefore, the correction of illumination inhomogeneities is highly desirable. Unfortunately, separating the bias field from the true underlying image is an underconstrained and ill-posed problem with fewer observations than free variables. For this reason, regularization of the problem is necessary.
Most existing bias correction methods assume that the bias field is multiplicative, slowly varying, and tissue independent. Many techniques ignore the noise and apply a log transform to make the bias field additive. The known illumination correction methods can be categorized in the following groups: filtering, segmentation based, surface fitting, and other methods.
Illumination inhomogeneities are generated during acquisition process in systems with different modalities. Here, the proposed illumination correction method will be applied on retinal images from confocal scanning laser ophthalmoscope (CSLO). This correction is an important preprocessing task in image segmentation and/or multimodal registration [1][2][3][4].
The earliest bias correction techniques were based on phantoms [5] with known structure; this approach enables estimation of the bias field, its inversion and removal. Provided the calibrated phantom is available, this method can be considered the ground truth but in real situations is not often applicable.
Linear filtering methods [6,7] try to estimate and eliminate the additive bias of the image using unsharp masking with defocus larger than any object in the image. These techniques can be extended, via nonlinear morphological approach to the multiplicative bias field estimation [6]. Homomorphic unsharp filtering methods [8][9][10] assume a separation of the low frequency bias field from the higher frequencies of the image structures. International Journal of Biomedical Imaging Some simple methods, as [11], require expert supervision, which is time-consuming and often too subjective. In this approach, the expert specifies some parameters related to the intensity models of the different tissues and forms a list of intensity values and locations related to the background or to object classes. Then, the bias field over the image may be obtained by least-squares fitting to the intensity values at the preselected points. These methods are closely connected to segmentation-based methods. Authors of [12] combine shading correction with adaptive segmentation. They use a fuzzy C-means algorithm to allow labeling of a pixel to be influenced by its immediate neighbors. Main benefits of this approach are high robustness to salt-and-pepper noise and computational efficiency of the algorithm.
The expectation maximization (EM) algorithm is proposed in [13] to compute iteratively the optimal smooth bias field corrupting the data based on classification into several tissue classes. Their formulation includes the bias distortion in the statistical model of the pixel distribution, that is, the bias field influences the distribution by locally changing its mean value. The algorithm iterates two steps, the E-step calculating the posterior tissue probabilities, and the M-step estimating the bias field. The Styner's method [14] is based on a simplified model of the imaging process, a parametric model of a tissue class statistics, and a polynomial model of the inhomogeneity field. The estimation of the parametric bias field is formulated as a non-linear energy minimization problem solved by an evolution strategy.
Segmentation-based approaches raise the problem of selecting the number of classes, which have to be explicitly modeled. Furthermore, these algorithms unfortunately tend to converge to a local nonoptimal minimum for some bias configurations, especially when more than two tissue classes are modeled [15]. The main problem of segmentation-based technique is that the determination of specific class-conditioned intensity is difficult when the image is incorrectly illuminated; however, the illumination correction is the main aim. Also, these methods usually assume that the intensity distribution of an image is normal and given by distributions of individual tissues, which may be often invalid. The above drawbacks are even more serious when correcting pathological image data. Therefore, Likar et al. [16] used the closed connection between intensity nonuniformity correction and segmentation, and proposed a method, which iteratively interleaves them, so that both of them gradually improve, until the final correction and segmentation is reached. The method is based on iterative minimization of the class square-error of intensity distributions caused by non-uniformity and on a nonparametric segmentation method. The method makes no assumption on the distribution of tissue intensity and does not require initialization.
To overcome the segmentation problems, bias correction methods not requiring the segmentation were designed based on a chosen image quality criterion. In [17], a presumed histogram matching method is presented. Authors of [18] suggested the iterative optimization method, which seeks the smooth multiplicative field that maximizes the higher frequency content of the distribution of tissue intensity. It requires a parametric model for the bias field but not a decomposition of the intensities. The quality criterion derived from the information theory has been also used in several applications. In [19], a polynomial multiplicative and additive shading model was introduced employing the cost-function based on image entropy and Powell optimizer. Mangin [20] uses image entropy combined with a measure of the field smoothness as the image quality cost function. Authors of [14] model the bias field using Legendre polynomials and use the energy function based on a multiclass model estimator and evolutional search algorithm.
Specific methods of illumination correction were proposed in the frame of retinal image processing and analysis. Simple and fast methods using large-kernel median filter to obtain a low-pass correction coefficients were used for CSLO image preprocessing in [21]. Authors of [22] model the bias field (background image) of a fundus (basic retinal) image as a white Gaussian random field and use Mahalanobis distance for background pixel classification. Contrast normalization using high-pass filtered image is used in [2] as one step of microaneurysm detection procedure. Additive model of nonuniform illumination is used in [3], together with adaptive histogram equalization.
Other approaches exist, for example, in applications including illumination correction for face recognition [23,24] and restoration of digitized documents [25,26]. The latter application frequently uses multiplicative model and specific properties of the processed images (e.g., illumination edges or geometrical distortion).
In this paper, we focus our attention on methods estimating the parametric illumination field using the quality criterion derived from the information theory. We use a multiplicative model of nonuniform illumination and parametric local bias model for formulation of criterion function and its derivatives (Section 2). The results are presented in Section 3 for different optimizers, which were tested for two image sets-artificially illuminated images and real CSLO images. The assessment of the results concludes the paper in Section 4.

Acquisition Model.
For the purpose of illumination correction process, we need a model of image creation. We assume, that each tissue class (vessels, optic disc, retinal surface) has a different mean value ρ(x) of the property measured by the imaging device. Moreover, every class of tissue has a characteristic texture, which can be modeled by the additive noise n tiss (x). The ideal output signal o(x) therefore consists of piecewise constant values plus additive noise. Due to finite size of the point spread function h(x) of the imaging device (different from the ideal Dirac impulse), this ideal signal is corrupted by convolution with h(x) and with the noise generated by the device n(x) (thermal or electronic noise and the noise coming from digitization during acquisition process). The overall equation of the observed image data can be formalized as follows (b(x) being the illumination), see Figure 1: (1) On the right-hand side of (1), we neglected the smoothing due to the imaging properties not playing a substantial role in our problem (the disturbance need not be taken into account in the retinal applications).
Thus, under the assumption that the bias involved in image creation process is multiplicative, and that it is the only important disturbance, we can reconstruct the original signal as where b(x) is the optimal illumination bias model controlled by parameters Φ i forming the vector Φ, while b −1 (x) is its reciprocal value. Example of a line input signal (measured line profile) can be seen on Figure 2. The choice of a proper bias model b(x) as well as of an appropriate criterion function, evaluating how well the non-homogenous illumination of the image is corrected for current values of the bias model parameters, is crucial.

Bias Model.
In this section, we extend the Likar's algorithm [19] and applied the modified algorithm for retrospective shading correction of retinal images. We derive an expression for the criterion of quality of illumination compensation based on Shannon's entropy and on Parzen windowing probability estimation. Further, as a novelty, we derive analytical expressions for derivatives of this criterion with respect to parameters of the used B-spline multiplicative illumination model. The analytical expressions ease substantially the optimization calculations compared to purely numerical evaluation of derivatives. Generally, in accordance with [14,19], the intensity transformation performed by the reciprocal bias model can be defined by a linear combination of K smooth basis functions where Φ i are the parameters of the transform defining the contribution of each basis. In order to regularize the problem of finding the bias field that would optimally correct the illumination inhomogeneity, we used the bases in the following form: where q i (x) are smooth polynomial basis functions, m c(x) is a mean correction coefficient, and n c(x) is a normalization coefficient. The correction coefficients are defined using constraints presented in [19], where the mean preserving condition preventing a global shift of the mean intensity of the resultant corrected image is formulated as Here, Θ is the size of the image domain Ω (or region of interest: ROI) in pixels. Further, as we use the multiplicative bias model, the final intensity is differently sensitive to values of different parameters, thus producing non-linear logarithmic scale (see [19] for detailed description). This fact is undesirable in optimization. Therefore, to assure that every parameter has the same influence on the intensity transform, a normalization constraint for every base function kernel We have found that the illumination distortion of retinal images shows a significant spatial variance, while both illumination models mentioned in [14,19] have rather a global character influencing the whole image despite possibly only locally defined distortions (e.g., when the image is corrupted in its upper left corner only, the polynomial model may significantly influence the opposite corner as well). Therefore, to model such local distortions, high-order models have to be used, which results in a high number of parameters to optimize, besides a danger of oscillatory character of the functions.
Hence, we decided to use the locally defined meancorrected and normalized reciprocal bias model based on Bsplines, formalized as follows: For the case of two dimensions, is a separable dth order B-spline kernel, P is a set of n x1 × n x2 control points of the transform regularly deployed over the extended image domain Ω with the spacing h=[h 1 , h 2 ] (h k possibly different but constant for each dimension, see Figure 3 for illustrative description). The division indicated in the expression for Δx in (7) should be understood element-wise; thus x i /h are the integer control-point indices. The Φ(xi) = φ i are scalar coefficients corresponding to the control point. These coefficients are the parameters of the bias model and express the weight of influence of each B-spline basis function. As we use the 3rd order B-spline, which is nonzero only for Δx k < −2, 2 > ∈, only 4 × 4 basis functions influence the computation of a current point x for two dimensions. For efficient implementation, the grid of control points has to be extended beyond the image borders. m c(xi) = m c i are coefficients ensuring the mean conservation condition (5) and n c(xi) = n ci are coefficients ensuring the normalization of the parameters (6); both these coefficients are precomputed ahead of the shading 4 International Journal of Biomedical Imaging correction. From the mean preserving condition (5) and from (7), we have which can be configured into The nontrivial solution which provides for each mean preserving constant m c i corresponding to a control point i with coordinates x i . Similarly, we can derive the normalization condition from (6) and (7) 1 Θ x∈Ω from where for each normalization constant n c i . Details of normalization coefficients computation in case of polynomial basis can be found in [19]; here we modified this approach for the B-spline representation. Finally, the restored image o(x) is provided using (2)

Criterion Function.
In order to find parameters of the bias model describing the undesirable illumination, we need to define a criterion, with respect to which the parameters would be optimized. In the works of Likar et al. [19] and Mangin [20], the image entropy is shown to be a suitable criterion. Their idea is that the illumination is additional information added to the information included in the original signal o(x) and because we would like to remove the illumination bias, the information content of the corrected image should be lower than that of the distorted image s(x). Therefore, when looking for parameters of the reciprocal bias model eliminating the non-homogenous illumination, the Shannon's entropy H(.) of the resulting image o(x), which is a measure of the information amount, should be minimized. Here, P(k) is the probability of intensity k appearing in any pixel of o(x). Although the brightness k is a discrete variable in reality (represented by 8 bits, i.e., 256 values), it may be well approximated by continuous variable κwith a probability density p(κ). Then, in order to analytically derive the derivatives of the criterion, we may assume that the amount of information can be described by the integral version of (14) for the continuous variable κ as However, the probability density p(κ) is not available explicitly and must be estimated from the image data transformed by the current parameters Φ, providing that the image has been generated by a homogeneous stochastic field. For this purpose, we used the Parzen windows (PW) technique [20], also known as kernel density estimator. In this scheme, the density is constructed by taking intensity samples o(x; Φ) from the transformed image and super-positioning kernel functions β (3) centered on the elements of Ω as illustrated in the Figure 4. More formally, is the probability density estimate at intensity value κ, which is obtained using PW, β (3) is 3rd order B-spline kernel, z is a constant parameter defining the width of the intensity class by controlling width of the super-positioned B-spline kernel (analogue to histogram bin size), and Θ is a normalization coefficient assuring that the integral of p(κ) is equal to one. This method can be looked at like at a convolution of impulses at different intensities with a spline kernel (see Figure 4). The advantage of this method is the possibility of expressing the probability derivatives analytically and also the possibility to speed-up the computation by taking into account only a subset of the whole sample set defined in the region of interest Ω.
Thanks to our formulation of the criterion (15) based on PW, both the probability estimation and the intensity transformation are defined continuously and therefore we can analytically derive partial derivatives of the criterion with respect to components of the parameter vector Φ. As H is a compound function of the form its derivative with respect to an element of Φ can be expressed as From here on, the ith control point (a linearly decoded control point of the rectangular grid) is characterized by a vector index k = [k 1 , k 2 ] with appropriate indices along both axes. If n x1 is the number of control points along the first dimension, then i = k 1 +n x1 k 2 .
By applying product rule, we obtain from (15) Calculation of this integral is approximated using the rectangle rule as Further from (16), The derivative of the B-spline kernel can be expressed using B-spline properties [27] as 6 International Journal of Biomedical Imaging The image o(x; Φ) derived from the observed image s(x) by application of the intensity transform with parameters Φ can be expressed using (2) and (7) as and for its partial derivatives we can write (in two dimensions): Finally, we have for the derivatives of the criterion where The grid controlling the compensation field is chosen adequately sparse as required by the smooth character of the to be compensated illumination unevenness. An example illustrating behaviour of the criterion function H(κ, Φ) as dependent on one element of the parameter vectorconcretely Φ [0,2] at i = 2, that is, for k = [0, 2]-topright image corner on 3 × 3 grid covering the image-is depicted on Figure 5(a) showing a distinct minimum. Below, also the course of partial derivative of H(κ, Φ) is shown, both as derived analytically by the described approach and as calculated numerically via finite differences. There is a good agreement between both versions namely, in the important area around the; however, important advantages of analytical derivatives are faster evaluation and smooth behavior in the parameter space. On Figure 5(b), the influence of this parameter on the illumination correction is illustrated for five particular values of Φ [0,2] influencing namely the correction in the upper right corner of the image; obviously, the best corrected image corresponds to the minimum of the criterion thus to the zero position of its derivative.

Optimization.
The overall illumination correction algorithm is based on the reciprocal illumination model b −1 with the parameters minimizing the entropy criterion H. The optimization aims at finding the optimum parameter vector The block diagram of the iterative algorithm is depicted on Figure 6. Various types of optimization techniques were studied in order to find the optimizer best suited for the particular properties of the used optimization criterion in the problem of non-homogenous illumination correction. The tested methods were namely, downhill simplex [28] (Amoeba), Powell's direction set [28], controlled random search (CRS) [29], gradient descent (GD) [28,30], conjugate gradient (CG) [28], and two versions of the limited memory Broyden, Fletcher, Goldfarb, Shannon (LBFGS) methods [31].
The comparison of results of the individual methods can be found in the next paragraph.

Experiments and Results
The proposed algorithm is supposed to be able to deal with non-homogenous illumination of retinal image data (384×384 pixel, 10 μm/pixel) obtained by means of the CSLO (Heidelberg Retinal Tomograph HRT II). In order to check the efficiency of the designed method on images with a known illumination inhomogeneity, a model of the HRT II image signal was designed, using segmentation of real HRT II images into five object classes with the estimated characteristic intensity values assigned to the objects as follows.
0-black background introduced by HRT II; 60-vessels; 120-vessel centers (brighter due to reflections); 90retinal tissue; 30-optic disc; 250-optic disc cup). Further, according to the image acquisition model (1), random tissue noise n tiss (uniformly distributed with standard deviation σ = 28) was added and the resulting image was convolved with the estimated point spread function of the real imaging system, approximated by 2D Gaussian kernel with standard deviation σ = 1. Finally, an artificial multiplicative illumination field controlled by a mesh of 3×3 nodes on the image area with random parameters was applied (see Figure 7).
A set of 40 simulated images was created via this model. Normalized parameters of the illumination field were uniformly distributed on interval [−15, 15]. Then, various optimizers were tested and evaluated with respect to achieved quality of the retrospective compensation of the illumination unevenness. The compensation quality was quantified using the fact that the known artificially introduced multiplicative illumination field and the final illumination correction field should be ideally inverse. Therefore, we can define the postcorrection illumination error ξ post as  International Journal of Biomedical Imaging  where b(x) is the introduced illumination field, b −1 (x) is the found correction of the illumination, and Θ is the size of the image domain Ω. The mean precorrection illumination error was set to ξ pre = 0.2875, that is, the intensity of the ideal image was corrupted on average by approximately 30%. The seven optimization algorithms mentioned in the previous section were tested and the average results obtained when using the simulated image set are summarized in Table 1. Here, the achieved after-correction total mean errors are presented, together with the computing times per image, numbers of iterations, and the achieved percentage suppression of the illumination unevenness. Examples of error fields resulting from these optimizers can be seen on Figure 8. The Amoeba, Powell and CRS methods need to evaluate the criterion value only; the derivatives of the criterion need not be computed. On the other hand, these three methods converge slowly and the total computational demands are higher compared to GD, CG, and LBFGS. Moreover, the worst problem of Amoeba, CRS, and especially of Powell algorithm was the unreliable convergence. The Powell optimizer was usually unable to find the optimum parameters of the correction illumination field (see Table 1). The false convergence to a side extreme is less frequent in case of smaller distortions introduced but even in this case the gradient-based algorithms LBFGS and GD outperform the nongradient algorithms. This can be explained by smoother description of the shape of the multidimensional cost function when using the analytical derivatives. Surprisingly, the simplest gradient descent optimizer with variable step proved to be the best with respect to the correction precision, suppressing the illumination errors by about 67 per cent on average (from ξ pre = 0.2875 to ξ post = 0.0926). The LBFGS optimizer was the fastest but the average achieved suppression of the simulated illumination error was only 51 per cent.
In frame of experiments, we tested also the influence of different number of image samples used for probability approximation. As can be seen in Table 2, even as little as 15% (about 20 000) of the available samples are statistically powerful enough to provide good approximation of probability density of image intensities. This results in faster computation of the criterion.
The next set of experiments concerned real retinal images with clearly visible illumination inhomogeneities that however were not known. The algorithm was tested on the set of 336 real (clinically obtained) retinal images. Obviously, the quality evaluation of non-homogenous illumination compensation could only be performed subjectively, due to lack of the golden standard (i.e., ideal illumination in each case). In majority of cases (about 95%), a substantial improvement was visible; in the remaining cases, the image remained practically unchanged, that is, no image was further distorted by the correction. On Figure 9, a result of the designed algorithm as applied to a real retinal HRT II image is illustrated; the intensity profiles clearly show the successful bias correction.

Conclusions
A method for efficient illumination correction was proposed, implemented, and verified: quantitatively on simulated images with known deterioration and qualitatively on an extensive set of real retinal images lacking this knowledge. It is based on estimating a B-spline polynomial shading model, the inversion of which provides correction of the input image, optimal in sense of the used informationbased criterion-Shannon's entropy. Previously published methods using a similar principle were modified namely, as to the type of the used correction function concerns, and as a novelty, the computation of the criterion is based on probability distribution estimates by Parzen windowing method, which enables us to derive analytical expressions for derivatives of the criterion. This consequently substantially speeds up the computations. Along with the efficient Bspline distortion model, it results in the possibility of efficient usage of gradient-based optimizers. We compared the efficiency and precision of three non-gradient and four gradient-based optimizers and found the classical gradient descent optimizer with variable step as the best for the purpose of illumination correction, formulated in the suggested way.
The quantitative tests were done on an image set artificially created with respect to characteristics of the imaging device (the method is primarily aimed at improving quality of retinal images taken by means of the HRT II CSLO); these tests have shown that the designed algorithm is capable of removing up to about 70% of artificially introduced illumination variability. Finally, the method was successfully qualitatively tested on a set containing 336 real CSLO clinically obtained retinal images.
According to results of some further tests, involving subsequent registration of multiple images, preprocessing of the retinal images by the proposed algorithm had a clearly positive influence on reliability of the registration of the images using the registration method developed by our group [1], when compared to registration results concerning the unprocessed images. A similar positive effect can be expected for even higher types of image analysis following the presented correction method.