improper scene illumination as well as nonideal acquisition conditions due to for example , misadjusted imaging system can introduce severe distortions into the resulting 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 under - constrained and ill - posed problem with fewer observations than free variables . 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 pre - processing task in image segmentation and/or multimodal registration [ 14 ] . the earliest bias correction techniques were based on phantoms   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 . homomorphic unsharp filtering methods [ 810 ] assume a separation of the low frequency bias field from the higher frequencies of the image structures .   some simple methods , as , 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 . 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   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   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 . 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 . therefore , likar et al .   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 , a presumed histogram matching method is presented . authors of   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 , a polynomial multiplicative and additive shading model was introduced employing the cost - function based on image entropy and powell optimizer . mangin   uses image entropy combined with a measure of the field smoothness as the image quality cost function . authors of   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 . authors of   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   as one step of microaneurysm detection procedure . additive model of nonuniform illumination is used in , 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 . 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 ntiss(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)s(x)=[((x)+ntiss(x))b(x)]h(x )  + n(x)o(x)b(x)+n(x ) . 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 
 ( 2)o^(x)=s(x)b^(x    )n(x)b^(x    )s(x)b^(x    )=s(x)b^1(x     ) , 
							 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 . in this section , we extend the likar 's algorithm   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 ri(x ) . ( 3)b^1(x)=i=1kiri(x ) , 
							 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 : 
 ( 4)ri(x)=qi(x)mc(x)nc(x ) , 
							 where qi(x ) are smooth polynomial basis functions , c(x ) is a mean correction coefficient , and c(x ) is a   normalization coefficient . the correction coefficients are defined using constraints presented in , where the mean preserving condition preventing a global shift of the mean intensity of the resultant corrected image is formulated as 
 ( 5)1xs(x)=1xs(x)b^1(x ) . 
							 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   for detailed description ) . therefore , to assure that every parameter has the same influence on the intensity transform , a normalization constraint for every base function kernel ri(x ) is introduced as 
 ( 6)1x|s(x)ri(x)|=1 ,  i , 
  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 mean - corrected and normalized reciprocal bias model based on b - splines , formalized as follows : 
								 ( 7 )  b^1(x)=1+xip(xi)[((d)(xi)c(xi )  m )  nc(xi ) ] ,  where  xi = xxih . for the case of two dimensions , (x ) = (x1 )  (x2 ) is a separable dth order b - spline kernel , p is a set of nx1  nx2 control points of the transform regularly deployed over the extended image domain  with the spacing h=[h1 , h2 ] ( hk 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 xi / h are the integer control - point indices . the (xi ) 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 xk < 2,2> , only 4  4 basis functions influence the computation of a current point x   for two dimensions . for efficient implementation c(xi ) = ci are coefficients ensuring the mean conservation condition ( 5 ) and c(xi ) = ci 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 
 ( 8)1xs(x)=1x{s(x)[1+i=1ki((d)(xi)  mci )  nci ] } , 
							 which can be configured into 
 ( 9)1i=1k{ici  n[xs(x)((d)(xi)ci  m)]}=0 . the nontrivial solution ( i  0 ) of ( 9 ) is 
 ( 10)xs(x)((d)(xi)ci  m)=0 ,  i , 
							 which provides 
 ( 11)ci  m=x  s(x)(d)(xi)x  s(x ) 
							 for each mean preserving constant ci corresponding to a control point i with coordinates xi . similarly , we can derive the normalization condition from ( 6 ) and(7 ) 
 ( 12)1x|s(x)((d)(xi)ci  m )  nci|=1 ,  i 
							 from where 
 ( 13 )  nci=1x|s(x)((d)(xi)ci  m)| , 
							 for each normalization constant ci . details of normalization coefficients computation in case of polynomial basis can be found in ; here we modified this approach for the b - spline representation . finally , the restored image o^(x ) is provided using ( 2 ) by the intensity transform b^-1(x| ) with the found parameter vector . because this transformation may generally produce out - of - range 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 . 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 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 ) , 
								 ( 14)h=k  p(k)log  ( p(k ) ) , 
							 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 
								 ( 15)h=p()log  ( p())d. 
							 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 , 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  centered on the elements of  as illustrated in the figure 4 . more formally , 
								 ( 16)p(;)=1x(3)(o^(x;  ) z ) 
							 is the probability density estimate at intensity value  , which is obtained using pw ,  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 
								 ( 17)h(;)=f(p((o^( ) ) ) ) , 
							 its derivative with respect to an element of  can be expressed as 
								 ( 18)h(;)k=hppo^o^k .   from here on , the ith control point ( a linearly decoded control point of the rectangular grid ) is characterized by a vector index k = [ k1 , k2 ] with appropriate indices along both axes . if nx1 is the number of control points along the first dimension , then   i = k1+n x1k2 . by applying product rule , we obtain from ( 15 ) 
								 ( 19)h(;)p=(log  ( p(;))+1)d. calculation of this integral is approximated using the rectangle rule as 
								 ( 20)h(;)pk(log  ( p(k;))+1 ) . further from ( 16 ) , 
								 ( 21)p=x1. 
							 the derivative of the b - spline kernel can be expressed using b - spline properties   as 
								 ( 22)(3)()=(2)(+12)(2)(12 ) . 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 
								 ( 23)o^(x;)=s(x){1+xip(xi)[((n)(xi)c(xi )  m )  mc(xi ) ] } ,  where  xi = xxih 
							 and for its partial derivatives we can write ( in two dimensions ) : 
 ( 24)o^(x;)k=xip{s(x)1cin[(3)(x1xk1h1)cim ] 
													 1cin[(3)(x2xk2h2)cim ]  } . finally , we have for the derivatives of the criterion 
 ( 25)h(;)kp(;)k(log  ( p(;))+1 ) , 
							 where  
 ( 26)p(;)k=x1{xip(3)()|=(s(x)b1(x)/z)  s(x)1nci[(3)(x1xk1h1)mci]1nci  [(3)(x2xk2h2)mci]}. 
  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 i = 2 , that is , for k = [ 0 , 2]top - right 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  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 . the overall illumination correction algorithm is based on the reciprocal illumination model b with the parameters minimizing the entropy criterion h. the optimization aims at finding the optimum parameter vector 
 ( 27)opt = argmin  {h(s(x)b1(x))}. 
							 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   ( amoeba ) , powell 's direction set , controlled random search ( crs ) , gradient descent ( gd ) [ 28 , 30 ] , conjugate gradient ( cg ) , and two versions of the limited memory broyden , fletcher , goldfarb , shannon ( lbfgs ) methods .   the comparison of results of the individual methods can be found in the next paragraph . the proposed algorithm is supposed to be able to deal with non - homogenous illumination of retinal image data ( 384384  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 . 0black background introduced by hrt ii ; 60vessels ; 120vessel centers ( brighter due to reflections ) ; 90retinal tissue ; 30optic disc ; 250optic disc cup ) . further , according to the image acquisition model ( 1 ) , random tissue noise ntiss ( 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 33 nodes on the image area with random parameters was applied ( see figure 7 ) . 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 
						 ( 28)post=1x(1b(x)b^1(x))2 , 
					 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 defined as 
						 ( 29)pre=1x(1b(x))2 
					 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 . 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 . 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 information - based 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 b - spline 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 , 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 .
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 .