Research Article  Open Access
Retrospective Illumination Correction of Retinal Images
Abstract
A method for correction of nonhomogenous illumination based on optimization of parameters of Bspline 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 splinebased shading model. This allows us to express the derivatives of the entropy criterion analytically, which enables efficient use of gradientbased optimization algorithms. Seven different gradient and nongradientbased 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 gradientbased 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 nonhomogenous 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.
1. 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 illposed 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–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–10] assume a separation of the low frequency bias field from the higher frequencies of the image structures.
Some simple methods, as [11], require expert supervision, which is timeconsuming 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 leastsquares fitting to the intensity values at the preselected points. These methods are closely connected to segmentationbased methods. Authors of [12] combine shading correction with adaptive segmentation. They use a fuzzy Cmeans algorithm to allow labeling of a pixel to be influenced by its immediate neighbors. Main benefits of this approach are high robustness to saltandpepper 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 Estep calculating the posterior tissue probabilities, and the Mstep 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 nonlinear energy minimization problem solved by an evolution strategy.
Segmentationbased 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 segmentationbased technique is that the determination of specific classconditioned 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 squareerror of intensity distributions caused by nonuniformity 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 costfunction 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 largekernel median filter to obtain a lowpass 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 highpass 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.
2. Methods
2.1. 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 . 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:
On the righthand 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 is the optimal illumination bias model controlled by parameters forming the vector , while 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 as well as of an appropriate criterion function, evaluating how well the nonhomogenous illumination of the image is corrected for current values of the bias model parameters, is crucial.
(a)
(b)
2.2. 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 Bspline 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 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 (x) are smooth polynomial basis functions, (x) is a mean correction coefficient, and 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 nonlinear 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 r_{i}(x) is introduced as
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, highorder 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 th order Bspline kernel, 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 elementwise; thus x_{i}/h are the integer controlpoint indices. The 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 Bspline basis function. As we use the 3rd order Bspline, which is nonzero only for ∈, 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. are coefficients ensuring the mean conservation condition (5) and are coefficients ensuring the normalization of the parameters (6); both these coefficients are precomputed ahead of the shading correction. From the mean preserving condition (5) and from (7), we have which can be configured into The nontrivial solution (Φ_{i }≠ 0) of (9) is 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) from where for each normalization constant . Details of normalization coefficients computation in case of polynomial basis can be found in [19]; here we modified this approach for the Bspline representation. Finally, the restored image is provided using (2) by the intensity transform with the found parameter vector .
Because this transformation may generally produce outofrange intensity values, we use intermediate image representation with extended intensity range. The resultant image is finally computed using linear contrast transformation converting the intermediate image back into the original intensity range.
2.3. 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 nonhomogenous illumination, the Shannon’s entropy H(.) of the resulting image , 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 .
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 from the transformed image and superpositioning kernel functions 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, is 3rd order Bspline kernel, z is a constant parameter defining the width of the intensity class by controlling width of the superpositioned Bspline 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 speedup 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 = [, ] with appropriate indices along both axes. If n_{x1} is the number of control points along the first dimension, then i=+n_{x1}.
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 Bspline kernel can be expressed using Bspline properties [27] as
The image 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 vector—concretely at , that is, for —topright image corner on 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 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.
(a)
(b)
2.4. 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 nonhomogenous 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.
3. Experiments and Results
The proposed algorithm is supposed to be able to deal with nonhomogenous 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); 90—retinal tissue; 30—optic disc; 250—optic disc cup). Further, according to the image acquisition model (1), random tissue noise (uniformly distributed with standard deviation ) 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 . 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)
(b)
(c)
(d)
(e)
A set of 40 simulated images was created via this model. Normalized parameters of the illumination field were uniformly distributed on interval [, ]. 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 as where is the introduced illumination field, is the found correction of the illumination, and is the size of the image domain Ω. The mean precorrection illumination error defined as was set to , 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 aftercorrection 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 gradientbased 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 to ). The LBFGS optimizer was the fastest but the average achieved suppression of the simulated illumination error was only 51 per cent.
 
(a)
(b)
(c)
(d)
(e)
(f)
(g)
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 nonhomogenous 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.
(a)
(b)
(c)
(d)
4. 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 Bspline 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 gradientbased optimizers. We compared the efficiency and precision of three nongradient and four gradientbased 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.
Acknowledgments
This research has been supported by the National Research Centre DAR (sponsored by the Ministry of Education, Czech Republic) project no. 1M0572 and partly also by the research frame no. MSM 0021630513.
References
 L. Kubecka and J. Jan, “Registration of bimodal retinal images—improving modifications,” in Proceedings of the 26th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC '04), pp. 1695–1698, San Francisco, Calif, USA, September 2004. View at: Google Scholar
 A. D. Fleming, S. Philip, K. A. Goatman, J. A. Olson, and P. F. Sharp, “Automated microaneurysm detection using local contrast normalization and local vessel detection,” IEEE Transactions on Medical Imaging, vol. 25, no. 9, pp. 1223–1232, 2006. View at: Publisher Site  Google Scholar
 A. A.H. AbdelRazik Youssif, A. Z. Ghalwash, and A. A. S. AbdelRahman Ghoneim, “Optic disc detection from normalized digital fundus images by means of a vessels' direction matched filter,” IEEE Transactions on Medical Imaging, vol. 27, no. 1, pp. 11–18, 2008. View at: Publisher Site  Google Scholar
 R. Kolar, L. Kubecka, and J. Jan, “Registration and fusion of the autofluorescent and infrared retinal images,” International Journal of Biomedical Imaging, vol. 2008, no. 1, Article ID 513478, 11 pages, 2008. View at: Publisher Site  Google Scholar
 L. Axel, J. Costantini, and J. Listerud, “Intensity correction in surfacecoil MR imaging,” American Journal of Roentgenology, vol. 148, no. 2, pp. 418–420, 1987. View at: Google Scholar
 I. T. Young, “Shading correction: compensation for illumination and sensor inhomogeneities,” in Current Protocols in Cytometry, J. P. Robinson et al., Ed., vol. 1, pp. 2.11.1–2.11.12, John Wiley & Sons, New York, NY, USA, 2000. View at: Google Scholar
 R. Chrastek, G. Michelson, K. Donath, M. Wolf, and H. Niemann, “Vessel segmentation in retina scans,” in Proceedings of the 15th Biennial International EURASIP Conference Biosignal, pp. 252–254, Brno, Czech Republic, 2000. View at: Google Scholar
 B. H. Brinkmann, A. Manduca, and R. A. Robb, “Optimized homomorphic unsharp masking for MR grayscale inhomogeneity correction,” IEEE Transactions on Medical Imaging, vol. 17, no. 2, pp. 161–171, 1998. View at: Google Scholar
 R. B. Lufkin, T. Sharpless, B. Flannigan, and W. Hanafee, “Dynamicrange compression in surfacecoil MRI,” American Journal of Roentgenology, vol. 147, no. 2, pp. 379–382, 1986. View at: Google Scholar
 J. Haselgrove and M. Prammer, “An algorithm for compensation of surfacecoil images for sensitivity of the surface coil,” Magnetic Resonance Imaging, vol. 4, no. 6, pp. 469–472, 1986. View at: Google Scholar
 B. M. Dawant, A. P. Zijdenbos, and R. A. Margolin, “Correction of intensity variations in MR images for computeraided tissue classification,” IEEE Transactions on Medical Imaging, vol. 12, no. 4, pp. 770–781, 1993. View at: Publisher Site  Google Scholar
 M. N. Ahmed, S. M. Yamany, A. A. Farag, and T. Moriarty, “Bias field estimation and adaptive segmentation of MRI data using a modified fuzzy cmeans algorithm,” in Proceedings of the IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR '99), vol. 1, pp. 250–255, Fort Collins, Colo, USA, June 1999. View at: Google Scholar
 W. M. Wells III, W. E. L. Crimson, R. Kikinis, and F. A. Jolesz, “Adaptive segmentation of mri data,” IEEE Transactions on Medical Imaging, vol. 15, no. 4, pp. 429–442, 1996. View at: Google Scholar
 M. Styner, “Parametric estimate of intensity inhomogeneities applied to MRI,” IEEE Transactions on Medical Imaging, vol. 19, no. 3, pp. 153–165, 2000. View at: Publisher Site  Google Scholar
 R. Guillemaud and M. Brady, “Estimating the bias field of MR images,” IEEE Transactions on Medical Imaging, vol. 16, no. 3, pp. 238–251, 1997. View at: Google Scholar
 B. Likar, J. Derganc, and F. Pernuš, “Segmentationbased retrospective correction of intensity nonuniformity in multispectral MR images,” in Medical Imaging 2002: Image Processing, M. Sonka and J. M. Fitzpatrick, Eds., vol. 4684 of Proceedings of SPIE, pp. 1531–1540, San Diego, Calif, USA, February 2002. View at: Publisher Site  Google Scholar
 L. Wang, H.M. Lai, G. J. Barker, D. H. Miller, and P. S. Tofts, “Correction for variations in MRI scanner sensitivity in brain studies with histogram matching,” Magnetic Resonance in Medicine, vol. 39, no. 2, pp. 322–327, 1998. View at: Publisher Site  Google Scholar
 J. G. Sied, A. P. Zijdenbos, and A. C. Evans, “A nonparametric method for automatic correction of intensity nonuniformity in mri data,” IEEE Transactions on Medical Imaging, vol. 17, no. 1, pp. 87–97, 1998. View at: Google Scholar
 B. Likar, J. B. A. Maintz, M. A. Viergever, and F. Pernuš, “Retrospective shading correction based on entropy minimization,” Journal of Microscopy, vol. 197, no. 3, pp. 285–295, 2000. View at: Publisher Site  Google Scholar
 J.F. Mangin, “Entropy minimization for automatic correction of intensity nonuniformity,” in Proceedings of the IEEE Workshop on Mathematical Methods in Biomedical Image Analysis (MMBIA '00), pp. 162–169, IEEE Press, Hilton Head Island, SC, USA, June 2000. View at: Publisher Site  Google Scholar
 H. Niemann, R. Chrastek, B. Lausen et al., “Towards automated diagnostic evaluation of retina images,” Pattern Recognition and Image Analysis, vol. 16, no. 4, pp. 671–676, 2006. View at: Publisher Site  Google Scholar
 M. Foracchia, E. Grisan, and A. Ruggeri, “Luminosity and contrast normalization in retinal images,” Medical Image Analysis, vol. 9, no. 3, pp. 179–190, 2005. View at: Publisher Site  Google Scholar
 T. Chen, W. Yin, X. S. Zhou, D. Comaniciu, and T. S. Huang, “Total variation models for variable lighting face recognition,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 28, no. 9, pp. 1519–1524, 2006. View at: Publisher Site  Google Scholar
 J. Zhu, B. Liu, and S. C. Schwartz, “General illumination correction and its application to face normalization,” in Proceedings of the IEEE International Conference on Accoustics, Speech, and Signal Processing (ICASSP '03), vol. 3, pp. 133–136, April 2003. View at: Google Scholar
 M. S. Brown, M. Sun, R. Yang, Y. Lin, and W. B. Seales, “Restoring 2D content from distorted documents,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 29, no. 11, pp. 1904–1916, 2007. View at: Publisher Site  Google Scholar
 G. Meng, N. Zheng, S. Du, Y. Song, and Y. Zhang, “Shading extraction and correction for scanned book images,” IEEE Signal Processing Letters, vol. 15, pp. 849–852, 2008. View at: Publisher Site  Google Scholar
 M. Unser, A. Aldroubi, and M. Eden, “Bspline signal processing. Part I. Theory,” IEEE Transactions on Signal Processing, vol. 41, no. 2, pp. 821–833, 1993. View at: Publisher Site  Google Scholar
 W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes in C, Cambridge University Press, Cambridge, UK, 1995.
 M. M. Ali, A. Törn, and S. Viitanen, “A numerical comparison of some modified controlled random search algorithms,” Journal of Global Optimization, vol. 11, no. 4, pp. 377–385, 1997. View at: Google Scholar
 D. P. Mandic, “A generalized normalized gradient descent algorithm,” IEEE Signal Processing Letters, vol. 11, no. 2, pp. 115–118, 2004. View at: Publisher Site  Google Scholar
 D. C. Liu and J. Nocedal, “On the limited memory BFGS method for large scale optimization,” Mathematical Programming, vol. 45, no. 1–3, pp. 503–528, 1989. View at: Publisher Site  Google Scholar
Copyright
Copyright © 2010 Libor Kubecka et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.