Optimal Joint Decoding/deblurring Method for Optical Images

We present a novel joint nonuniform illumination estimation and deblurring method for bar code signals based on a penalized nonlinear squares objective function. The objective function is based on the proper parameterization of a bar code signal and nonuniform illumination as well as a regularization on the illumination using a smoothness penalty. By the minimization of the objective function, the proposed method simultaneously estimates the bar code signal and illumination in the spatial domain. In simulations and experiments, the proposed method showed improved performance compared with two conventional bar code decoding methods without deblurring or nonuniform illumination correction. In a few iterations, the proposed method was able to decode test bar code signals that were not decodable due to blurring or nonuniform illumination.


Introduction
Decoding information from an acquired bar code image requires accurate estimation of the true lengths of bar and space patterns in the bar code image since information is encoded in the lengths of the patterns.To do this, one of the most widely used methods involves detecting edges in the acquired image using an edge detector and subsequently estimating the lengths by computing the differences between the detected edges [1].This local edge based method often does not perform satisfactorily if the image is blurred due to defocusing since the true starting locations of the patterns do not correspond to the detected edges in the acquired image.This has been a concern in bar code signal processing.Researchers have therefore investigated various methods for the deblurring of a bar code signal including improved local methods such as a modified edge detector method [1] and a peak detection method [2].
Recently, drawbacks of the local methods have been noted, such as sensitivity to noise and limited accuracy in detecting the starting locations of patterns [3].Several researchers have investigated the statistical properties of speckle noise and its effect on the accuracy of the detected edge locations for laser based scanning systems [4,5].In addition, a thorough investigation on edge localization errors due to additive noise has been reported [6].On the other hand, to overcome the problems, researchers have also investigated alternative global approaches, such as a penalized least squares method [3] and a nonlinear least squares method based on the parameterization of a bar code signal using the locations of edges [7].The deblurring problem, which requires the restoration of a bi-level bar code signal and a convolutional blurring kernel from an acquired signal, is an underdetermined ill-posed inverse problem.Thus, the penalized least squares method incorporated a regularizing penalty function based on the total variation of a bar code signal [3].Alternatively, the nonlinear least squares method parameterized a bar code signal using edge locations under the assumption that the number of edges is known.Since the bar code signal was parameterized using a small number of parameters, the problem was considered well-posed although it lacks rigorous mathematical proof [7].Compared with the total variation method, the parametric method has the advantage that only a small number of parameters need to be estimated without requiring regularization, thereby improving speed without causing difficulty in tuning the regularization parameter [7].Although the approach is limited for bar code signals whose number of edges is known, there exists a class of bar code signals whose number of edges is predetermined, such as UPC and PDF417 [8].Moreover, in practical situations, the number of edges can be estimated by proper selection of edge candidates during edge detection.In fact, the assumption that the number of edges is known has been often used for modeling the interaction of edges [1,9].
Although the global approaches showed satisfactory performance for deblurring of bar code signals in simulations [3,7], the performance for real bar code images often suffers since the images are usually distorted by nonuniform illumination.To alleviate the problem, attempts have been made to improve lighting devices for bar code scanners [10,11].However, even with such additional hardware devices, the nonuniformity cannot be removed completely.Moreover, demand is increasing for a general purpose hand-held device with imaging capability, such as a mobile phone with a digital camera, having the capability of decoding bar codes [12].Images acquired by such a device may be more severely distorted by nonuniform illumination.The effect of nonuniform illumination on the intensities of an acquired image is a spatially varying multiplicative field acting on the intensities of an intrinsic image, since the acquired image is the product of the illumination and the reflectance of the intrinsic image [13,14].In the presence of such nonuniform illumination on an image, the performance of image analysis such as segmentation and registration is often degraded [15,16].We have observed a similar adverse effect of nonuniform illumination in the deblurring of a bar code signal.Therefore, for correct decoding of a bar code signal acquired under nonuniform illumination, both nonuniform illumination compensation and deblurring are essential.
Compared with the deblurring problem, the nonuniform illumination compensation problem has rarely been studied for bar code images.To our knowledge, no investigation on retrospective nonuniform illumination compensation has been reported, although several inventions on improving hardware devices for more uniform illumination have been published [10,11].The analogous problem, however, has been intensively studied in microscopy image processing [15,17], magnetic resonance (MR) image corrections [16,18,19] and document image processing [20].The nature of the nonuniform illumination compensation (often called shading correction) is the estimation of the unknown multiplicative field (nonuniform illumination for microscopy and document images, inhomogeneous magnetic field for MR images) and unknown intrinsic image from an observed image.Simultaneously estimating the nonuniform illumination and the intrinsic image in the spatial domain is a well known underdetermined ill-posed problem [14].Consequently, conventional methods often separate the illumination and the intrinsic image in the frequency domain or in the distribution domain (such as a histogram domain).For example, based on the assumption that illumination is slowly varying in the spatial domain, methods that separate illumination based on homomorphic filtering (lowpass filtering of the observed signal after taking the logarithm) have been investigated [15,18].The performance of these approaches is limited since the intrinsic image also has low frequency components.Other methods rely on the assumption that the distribution of the observed image is the blurred distribution of the intrinsic image due to nonuniform illumination.Such methods seek a multiplicative field that minimizes the entropy of the estimated intrinsic image [17] or that maximizes the frequency contents of the distribution of the estimated intrinsic image [19].
For nonuniform illumination compensation and deblurring for a bar code signal, we believe that the existing methods are not the most appropriate since they do not utilize the fact that a bar code signal is a bi-level waveform.In addition, since the existing methods only attempt to compensate for nonuniform illumination effects, deblurring and edge detection must be performed separately in order to decode the bar code.
In this study, we propose a joint nonuniform illumination estimation and deblurring method that simultaneously estimates nonuniform illumination and a bar code signal in the spatial domain.We parameterize a bar code signal using the locations of edges and model nonuniform illumination using a sum of B-spline functions that has been shown to be effective in approximating a continuous function [21].Since the minimization of square errors between an observed signal and a signal model can be ill-posed even though the number of parameters in our model is much smaller than the number of observations, we incorporate a regularizing penalty based on the a priori information that illumination is slowly varying in the spatial domain.The proposed method utilizes the fact that the true intrinsic image is a blurred bi-level image based on accurate modeling of a bar code signal.Moreover, since nonuniform illumination and the bar code signal are jointly estimated, the proposed method does not suffer from performance degradation in deblurring due to nonuniform illumination effects.

Model
We model an observed bar code signal under unknown illumination by the product of a blurred bar code signal and a spatially varying illumination field with additive noise as follows: where t denotes one-dimensional (1D) spatial location, * denotes the convolution operator, I(t), G(t) and x(t) denote unknown illumination, blurring kernel, and clean bar code signal, respectively, c > 0 denotes unknown constant background and n(t) denotes additive noise.Estimation of I(t), G(t), x(t) and c from y(t) is an ill-posed problem since the problem is underdetermined, and thus its solution is not unique [14].To overcome this difficulty, we parameterize the illumination, the blurring kernel, and the bar code signal using a small number of parameters.
Although the problem still may be ill-posed despite the parameterization, we choose this model since we can design an effective regularization method for the model, as will be discussed in the next section.We assume that the number of edges, K, in a bar code signal is known.Based on the assumption, we parameterize a clean bar code signal as follows: where e = [e 1 , e 2 , ..., e K ] denotes the locations of edges such that e i+1 > e i , and u(t) is the unit step function.We also parameterize the blurring kernel using the Gaussian function with unknown standard deviation σ > 0 as follows: Although the true blurring kernel is neither space invariant nor an exact Gaussian function [22], the simple, space invariant Gaussian model has been effectively used for modeling blurring kernels in many other investigations and applications including modeling the line spread function in flying spot laser based bar code scanners [1][2][3]6,7].We model spatially varying illumination using third order B-spline functions that are uniformly located in the spatial domain as follows: where β 3 (t) is a third order B-spline function centered at zero, L is the number of the Bspline functions, b = [b 1 , ..., b L ] denotes the coefficients for the B-spline function and l j , j = 1, ..., L denote uniformly distributed center locations of the B-spline functions.We predetermine the locations of the B-spline functions.In addition, we impose positivity constraints on the coefficients for the B-Spline functions to ensure that illumination is positive everywhere.The Bspline based model in Eq. ( 4) has been shown effective for approximating a continuous function by proper selection of the coefficients [21].With the models in Eq. (1)-Eq.( 4), the goal of joint nonuniform illumination estimation and deblurring of a bar code signal is the accurate estimation of e, b, c, and σ from y(t).Although accurate estimation of e is sufficient for correct decoding of a bar code signal, we think that accurate estimation of b, c, and σ is also necessary for the accurate estimation of e.If the parameters b, c, and σ are not accurately estimated, estimation of e with the inaccurately estimated parameters corresponds to estimation based on an incorrect signal model.It has been reported that the performance of an estimator based on an incorrect model deteriorates, i.e., the estimator may have larger bias and/or variance than that of an estimator based on a correct model [23].

Regularization
Based on the models defined in Eq. (1)-Eq.( 4), we consider the problem of estimating the parameter vector θ , where the number of total parameters N = K + 2 + L, from the samples of the observed signal y(t i ), i = 1, ..., M subject to the constraints such as e i < e i+1 , σ > 0, c > 0 and b j > 0. To solve the problem, one effective approach is determining a parameter vector θ that minimizes square errors between an observed signal and a model subject to the constraints.This leads to a constrained optimization problem.Although constraints do exist, we focus on unconstrained optimization since, for most real bar code signals, the constraints are not violated during optimization provided the initial estimate is close to the optimal parameter vector.Within the unconstrained optimization framework, the most natural approach is determining an estimator θ that minimizes the sum of nonlinear squares: where the objective function F(θ ) is defined as where f i (θ ) is the sampled residual defined as where t i denotes the sample locations, and [G(t i ; σ ) * x(t i ; e)] denotes [G(t; σ ) * x(t; e)] t=t i for notational simplicity.The nonlinear least squares problem in Eq. ( 5) can be ill-posed, even though the number of variables N is usually smaller than the number of observations M, since the minimizer may not be unique [24].Note that the minimizer is determined by zeroing the gradient of the objective function F(θ ) as follows: where T , and the Jacobian matrix For θ in a neighborhood of θ , if J( θ ) (J(θ )| θ = J( θ ) for notational simplicity) is rank-deficient, i.e., the rank of J( θ ) is smaller than N, the nonlinear least squares problem is ill-posed since the solution of Eq. ( 8) is not unique [24].Equivalently, the Hessian matrix ∇ 2 θ F(θ ) ∈ R N×N of the objective function F(θ ) can be singular at θ since, in a neighborhood of θ , the Hessian matrix is well approximated by the Jacobian as follows [25]: where (i, j) th element of the Hessian matrix is defined as follows: Therefore, it is difficult to apply the Gauss-Newton method or its variations to solve the nonlinear least squares problem since the approximated Hessian is not invertible [24,25].For an ill-posed problem, in order to determine unique solution and stabilize a numerical procedure for the solution, regularization based on a priori information about the solution is required [24,25].
The regularization is usually accomplished by adding a regularizing penalty function to the objective function in such a way that the approximated Hessian matrix becomes positive definite [26].
A situation causing the nonlinear least squares method with nonuniform illumination model to be ill-posed is when illumination is indistinguishable from the blurred bar code signal.For example, suppose that there are two different illuminations I(t i ; b 1 ) and I(t i ; b 2 ) that have the same shapes as the blurred bar code signal [G(t i ; σ 2 ) * x(t i ; e 2 ) + c 2 )] and [G(t; σ 1 ) * x(t i ; e 1 ) + c 1 )], respectively, such that: By the mean value theorem [27], there exists a nonzero vector h ∈ R 1×N and a parameter vector θ on the line connecting in the parameter space such that: where, 0 denotes the 1× M zero vector.Equation ( 13) implies that J( θ ) is rank deficient.To prevent singular Hessian matrices during optimization, we regularize I(t; b) such that illumination would not have a waveform similar to the blurred bar code model.Before designing a regularization method, we first consider the nonlinear least squares problem in Eq. ( 5)-Eq.( 7) using the uniform illumination model, i.e., the problem of estimating , where b > 0. For this model, all the B-spline coefficients b j , j = 1, .., L in Eq. ( 4) are set by b.The nonlinear least squares problem with uniform illumination model is equivalent to estimating σ , c, b, and a bi-level signal x(t i ) ∈ {0, 1}, i = 1, ..., M that has the total variation of K (number of edges).Note that the set of bi-level signals having total variation of K is the same as the set of bi-level signals composed of K edges.It has been shown that the estimation of x(t i ) having a bounded total variation from y(t i ) is well-posed with some constraints such that σ ≥ 0 and b ∈ P, where P is a compact subset of positive real numbers [3].One can also show that two different parameter vectors do not satisfy Eq. ( 12) by extending the uniqueness proof in [3].Based on the result, we assume that the nonlinear least squares with uniform illumination model is well-posed and the unique minimizer is determined by zeroing the gradient of the objective function.For most real bar code signals, this assumption is true since the constraints are not violated during optimization.Based on the assumption, the approximated Hessian at some θs in a neighborhood of the minimizer is positive definite, i.e., for all θ s ∈ R K+3 , the following inequality holds: where Next, we consider a regularization method for the nonlinear least squares method with nonuniform illumination model.From inequality (14), it is straightforward to show that where, is not a function of b j , it is straightforward to show the following equality: where θ = [ẽ, σ , c, b] with a b vector that does not necessarily have uniform B-spline coefficients.Based on inequality (14), Eq. ( 15) and Eq. ( 16), the approximated Hessian of the objective function F(θ ) satisfies the following inequality: Note that inequality (17) does not guarantee the positive definiteness of the approximated Hessian since θ T J( θ ) T J( θ )θ > 0 is required.To overcome the problem, we design a new objective function that includes a regularizing penalty function: where λ > 0 is a regularization parameter and the roughness penalty function R(•) is defined as follows: where R is the matrix representation of the roughness penalty function [28].The penalty function in Eq. ( 19) penalizes the changes of neighboring B-spline coefficients, thereby encouraging the estimated illumination to be smooth.It has been shown that if the change of neighboring Bspline coefficients is bounded by a constant α > 0, then the first order derivatives of nonuniform illumination in Eq. ( 4) are bounded by α/L where L is the distance between two neighboring B-spline coefficients [29].The first term in the objective function in Eq. ( 18) encourages errors in fitting data to the observations to be small, while the second term encourages estimated illumination to be smooth.The regularization parameter λ controls the trade-off between the fidelity to the observed data and the smoothness of the illumination.For θ in a neighborhood of the minimizer, the Hessian matrix of the objective function Φ(θ ) in Eq. ( 18) is well approximated as follows: We claim that the approximated Hessian is positive definite, satisfying the following inequality: The claim is true since for θ that has uniform B-spline coefficients, θ T J( θ ) T J( θ )θ > 0 by Eq. ( 17) while λ θ T Rθ = 0, and for θ whose B-spline coefficients are not uniform, λ θ T Rθ > 0 even if θ T J( θ ) T J( θ )θ = 0. Since the approximated Hessian in Eq. ( 20) is positive definite in a neighborhood of the minimizer, one may apply the Gauss-Newton algorithm to minimize the objective function in Eq. (18).We minimize the objective function using the standard Gauss-Newton algorithm that constructs a sequence of parameters {θ k } as follows: where Due to the constraints such as σ > 0, e i+1 > e i , i = 1, .., K, c > 1 and b j > 0, j = 1, .., L, one must solve the constrained nonlinear least squares problem.To do this, we have incorporated active set strategy [30] with the Gauss-Newton algorithm.However, in practice, no constraint was violated during iterations for most bar code signals that were correctly decoded.For a severely blurred and noisy image, we observed that certain constraints were violated during iterations.For these cases, the bar code was not correctly decoded.

Initial estimates
As in many optimization problems, accurate initial estimation is very important for an effective minimization of the objective function in Eq. (18).We estimate initial edge locations using local maxima and minima of the signal acquired after wavelet-based denoising of an observed signal, as described in the previous investigation [7].We also initialize σ using the distance from the last local maxima to the last local minima in the spatial domain, as described in the previous investigations [1,7].We estimate the magnitude of illumination using the difference between the maximum and minimum of the observed signal.Based on the assumption that the illumination is close to uniform, we initialize all the B-spline coefficients by the estimated magnitude of illumination.If true illumination were uniform and true σ were very small, the initial illumination estimate would be close to the true illumination.We also initialize the constant background c by dividing the minimum of the observed signal by the estimated magnitude of illumination.

Simulations
To evaluate the performance of the proposed joint nonuniform illumination estimation and deblurring method, we conducted simulations using UPC bar code signals.Bars and spaces in a true UPC bar code image can have one of four different lengths, each length being an integer multiple of a minimum length up to four times the minimum.For simplicity of notation, we denote the minimum length by length 1 and other lengths by length 2 through length 4 based on the ratios of their lengths to the minimum.Figure 1(a) shows a two-dimensional (2D) UPC bar code image representing code '012345678905'.Figure 1(b) shows a clean 1D signal of the same code with a constant background of unity (without loss of generality, we assume that the interval of interest starts at t = 0 and ends at t = 1).If Fig. 1(a) were an ideal image, Fig. 1(b) would correspond to a horizontal line scan of Fig. 1(a).In Fig. 1(b), the length of the bar code signal is 0.9406 (excluding the background) and the length 1 corresponds to 0.0099.
We synthesized the 1D blurred signal shown in Fig. 2(a) by convolving the signal shown in Fig. 1(b) with the Gaussian convolutional kernel whose σ is 0.007.We also generated two  synthetic nonuniform illuminations (shown in Fig. 2(b)) having Gaussian and linear functional forms, respectively.We multiplied each synthetic nonuniform illumination with the signal in Fig. 2(a), and added white Gaussian noise for generating synthetic observations.Figure 2(c) and Fig. 2(d) show two synthetic observations with the Gaussian and the linear illuminations, respectively.For comparison, we have included a synthetic observation with uniform illumination whose magnitude is unity in Fig. 2(c) and Fig. 2(d).
To decode the signal shown in Fig. 2(c), we applied three different methods: (1) the edge based method that attempts to decode the signal using initial estimates of edge locations, (2) the nonlinear least squares with uniform illumination model (NLS-UI) method and (3) the proposed method (the nonlinear least squares with nonuniform illumination model).For the proposed method, we placed 128 B-spline functions uniformly in the spatial domain.
Figure 3(a) shows the fitted bar code signal using the edge based method, which may be considered as the initial fitting for the proposed method (to simplify the illustration, the true c value was used instead of the initial one).Figure 3(b) and Fig. 3(c) show the fitted signals by using the NLS-UI method and the proposed method, respectively.Figure 3(d) shows the estimated illumination using the proposed method with the regularization parameter λ = 50.As shown in Fig. 3(a), Fig. 3(b) and Fig. 3(c), the fitted signal using the proposed method is closer to the observed signal than that using the other methods.Note that the fitted signal in Fig. 3(b) has larger fitting errors around peaks than that using the proposed method.In addition, the estimated nonuniform illumination in Fig. 3(d) was close to the true illumination that was used to synthesize the observation.
We also applied the three methods to the observation shown in Fig. 2(d), although the figures are not shown here.The results were similar to the results shown in Fig. 3: the proposed method gave the best match of the fitted signal to the observed signal.
The goodness of fit of the estimated signal to the observed signal is only one measure of the performance of each method.The goal of bar code deblurring is the accurate estimation of edge locations, and accurate estimation of the length of each bar and space pattern is the key element to correctly decode a bar code signal.Thus we compute the estimated lengths of patterns, Δ = [ Δ1 , Δ2 , ..., ΔK−1 ], whose elements are defined as follows: where ti is the location of each estimated edge.If the elements of Δ satisfy the following rela-  tionship, Δi1 < Δi2 < Δi3 < Δi4 , where i1, i2, i3, and i4 denote the indices for length 1, length 2, length 3, and length 4 in the true bar code signal, one may decode the bar code correctly by proper grouping of the estimated lengths. Figure 4(a), Fig. 4(b) and Fig. 4(c) show Δ for the signal shown in Fig. 2(c) using the edge based method, the NLS-UI method and the proposed method, respectively.For comparison, we have included the true length of patterns in the figures.As shown in Fig. 4(a), the edge based method was not able to decode the bar code signal correctly since Δ37 and Δ39 (length 1 in the true signal) are larger than Δ40 and Δ42 (length 2).In addition, Δ22 (length 3) is almost the same as Δ36 (length 4).The NLS-UI method was likewise not able to decode the bar code.In Fig. 4(b), Δ25 (length 2 in the true signal) was larger than Δ26 (length 3).In addition, Δ54 (length 2) was larger than Δ55 (length 3).Compared with the conventional methods, the proposed method showed significantly improved results.As shown in Fig. 4(c), the estimated lengths Δ by the proposed method are very close to the true lengths.Clearly, the bar code can be correctly decoded by grouping the estimated lengths using the threshold 1 through threshold 3 shown in Fig. 4(c).
We repeated the same simulations for the signal in Fig. 2(d) and obtained similar results.Neither the edge based method nor the NLS-UI method decoded the signal correctly, while the proposed method decoded correctly.For brevity, these figures are not shown here.To quantify the performance of each method in terms of correct estimation of edge locations, we computed the maximum error (ME) of the estimated edge locations as follows: where t true i denotes the true locations of edges.We use the ME for evaluating the performance of each method instead of usual mean square error since a large estimation error of a single edge location can result in the failure of decoding entire bar code [6].To evaluate the statistical properties of each method, we applied each of the three methods 100 times for the signal shown in Fig. 2(c) and Fig. 2(d) with different noise realizations, and computed the average of the ME.We repeated these simulations for different size blurring kernels and for signals with different signal-to-noise ratios (SNR).Table 1 and Table 2 summarize the results for the signals shown in Fig. 2(c) and Fig. 2(d), respectively.In the tables, we represent the level of blur by 2σ /T (T is the length 1 in the signal shown in Fig. 1(b)) since the edge localization error for a stripe (bar or space) due to a blurring kernel is approximately inversely proportional to the width of the stripe.As shown in the tables, the proposed method exhibited smaller ME than the existing methods in all simulations.Since the ultimate goal of bar code decoding is to decode a bar code signal correctly, we also evaluated the performance of each method by a decode rate that was computed by the percentage of correct decoding for the signal shown in Fig. 2(c) in 100 different noise realizations.We repeated the computation of the decode rate for different size blurring kernels and for signals with different SNR.Table 3 summarizes the results.As shown in the table, the proposed method showed significantly improved decode rate compared with the other methods in all simulations.An interesting observation was that, for blur level 1.21, the decode rate of the NLS-UI method for higher SNR signals (20dB and 25dB) was worse than that for lower SNR signal (15dB).This is due to the fact that the NLS-UI method was not able to decode noise free signal since the method was based on an incorrect model (uniform illumination model while true illumination was nonuniform).For low SNR signal, we suspect that large variability of estimated edge locations due to noise happened to result in correct decoding in some realizations.When modeling error was small, i.e., true illumination was close to unform illumination, we observed that the decode rate of the NLS-UI method improved as SNR increased although results for this are not included for brevity.
In the preceding simulations, the regularization parameter λ was determined manually.To investigate the effect of λ , we applied the proposed method to the signals shown in Fig. 2(c) and Fig. 2(d) while changing the value of λ .Figure 5(a) and Fig. 5(b) show the estimated nonuniform illuminations with different λ values for the signals shown in Fig. 2(c) and Fig. 2(d), respectively.As shown in the figures, small λ yielded rapidly changing illumination that is not close to the true illumination.As a result, a very small regularization parameter such as λ = 0.1 prevented the proposed method from accurately decoding the bar code.Using a very large λ value, such as λ = 1000 in Fig. 5(a) and λ = 10000 in Fig. 5(b), also led to incorrect decoding, yielding almost uniform illumination.In other words, if λ is too large, the proposed method yields results very similar to those of the NLS-UI method.However, for a wide range of λ such as 0.1 < λ < 1000, the proposed method was able to correctly decode the signals.
Figure 6(a) and Fig. 6(b) show the changes of the objective function and the ME of the estimated edge locations during optimization using the proposed method with different λ values for the signals shown in Fig. 2(c) and Fig. 2(d).As shown in Fig. 6(a), the proposed method converged to the minima within less than 10 iterations for every λ value.Note that the converged values of the objective function using smaller λ values are smaller than that using larger λ .This result is expected since a smaller regularization parameter implies more fidelity to the measured data.However, as shown in Fig. 6(b), the ME of the estimated edge locations with a small λ such as 0.1 or 1 is larger than that with λ = 50.This result implies that simply minimizing data fitting errors may yield physically unrealistic results, since rapidly changing illumination can better fit a noisy observation.Therefore, incorporating a priori information about the solution, such as the smoothness of illumination, is important for accurate decoding of a bar code signal.Figure 6(c) and Fig. 6(d) show the changes of the objective function and the ME of the estimated edge locations for the signal shown in Fig. 2(d).The results are similar to those in Fig. 6(a) and Fig. 6(b) in the sense that too small or too large λ yields larger ME than a wellchosen λ value such as 500.Although the performance of the proposed method depends on λ , it is worth while to note that the method was able to decode the bar code signals correctly over a wide range of λ .

Experiments
We acquired real bar code images using a digital camera (Canon EOS-350D).Figure 7(a) and Fig. 7(b) show two bar code images acquired under nonuniform illumination.In Fig. 7(a), illumination around the central portion of the image was brighter than around the outer parts.In Fig. 7(b), illumination was brighter on the left portion of the image than on the right portion.We scanned in horizontal lines across each image for decoding.Figure 7(c) and Fig. 7(d) show the scanned 1D signals from Fig. 7(a) and Fig. 7(b), respectively.Unlike the synthesized 1D signals used for simulations, space patterns in these images have higher intensity values than the bar patterns since the intensity levels of the brighter parts are higher in the acquired image.
We applied the two existing methods as well as the proposed method to the real bar code signals shown in Fig. 7(c).Figure 8(a), Fig. 8(b) and Fig. 8(c) show the fitted signals using the three methods.In addition, Fig. 8(d) shows the estimated illumination using the proposed method with λ = 100.Figure 8(a) through Fig. 8(c) show that the proposed method gives a fitted signal more closely matching the observed signal than the other methods do.This result is consistent with the results in the simulations.
Figure 9(a), Fig. 9(b) and Fig. 9(c) show the fitted signals for the observed signal shown in Fig. 7(d) using the three methods.In addition, Fig. 9(d) shows the estimated illumination using the proposed method.The figures demonstrate that the proposed method yields smaller fitting errors than the others.
To evaluate the performance of the three methods in terms of decodability, we computed the lengths of the estimated patterns, Δ, by each method.Figure 10( show the estimated Δ for the real bar code signal shown in Fig. 7(c) using the edge method, the NLS-UI method and the proposed method, respectively.Since the true length of a length 1 pattern is not known in a real image, we estimated the true length by dividing the sum of estimated lengths by 95 that is the total length of the true patterns when length 1 is unity (since the estimated length 1 is based on an estimated results, it can be different for different estimations).Using the estimated length 1, we generated correct Δ based on the ratios of true patterns to the estimated length 1, and included such correct Δ for comparison.
As shown in Fig. 10(a), the edge based method did not decode the bar code correctly, giving Δ17 (length 4 in the actual bar code) is smaller than Δ25 (length 3) and Δ12 (length 2) smaller than Δ16 (length 1).Similarly, the NLS-UI method did not decode the signal correctly.In Fig. 10(b), Δ1 (length 1 in the actual bar code) is larger than Δ8 and Δ10 (length 2).In addition, Δ4 (length 3) is smaller than Δ5 (length 2).Compared with the existing methods, the proposed method was able to estimate the lengths of the patterns more accurately as shown in Fig. 10(c).It is clear that the bar code can be correctly decoded using the thresholds shown in Fig. 10(c).
Figure 11(a), Fig. 11(b) and Fig. 11(c) show the lengths of patterns estimated by the three methods for the real signal shown in Fig. 7(d).Again, the edge method was not able to decode the signal, giving Δ11 and Δ13 (length 2 in the actual bar code) smaller than Δ19 (length 1).As shown in Fig. 11(b), the NLS-UI method also gave a number of errors, especially for the right hand part of the signal.For example, Δ47 (length 2) is smaller than Δ48 (length 1).In addition, Δ59 (length 1) is even larger than Δ36 (length 4).In contrast, Fig. 11(c) shows that the proposed method was able to correctly decode the signal.We also studied the effect of the regularization parameter in decoding the real bar code images.Figure 12(a) and Fig. 12(b) show the estimated nonuniform illuminations with different λ values for the real signals shown in Fig. 7(c) and Fig. 7(d), respectively.As expected, Fig. 12 shows that very large λ yielded almost uniform illumination, thereby suffering from the same problem as the NLS-UI method.This is consistent with the simulation results.On the other hand, with very small λ such as 10 −2 , the estimated illumination is not physically realistic, changing rapidly in the spatial domain.This is due to the fact that the penalized nonlinear least squares problem becomes increasingly ill-conditioned as λ becomes smaller [24], implying that the estimation with smaller λ is more sensitive to noise.As in the simulations, although very small or very large λ values give incorrect decoding, the proposed method correctly decodes bar codes for a broad range of λ values, such as 1 < λ < 1000.minima within 10 iterations for every case except for very small λ such as 10 −2 .For λ = 10 −2 , the numerical optimization procedure showed oscillatory behavior, likely due to the fact the approximated Hessian becomes close to singular [24].

Discussion
In the presence of blurring and nonuniform illumination in acquired bar code images, the local edge method and the NLS-UI method, which may be considered the same method described in [7] except for the slight extension of including background intensity, did not perform satisfactorily.The edge based method was not able to decode most blurred bar code signals even in simulations with uniform illumination.The NLS-UI method also failed to decode most bar code images acquired under nonuniform illumination.Compared with the two existing methods, the proposed joint nonuniform illumination estimation and deblurring method showed significantly improved results in both simulations and experiments.
One may conceive of a two-step method that first applies one of the existing nonuniform illumination correction methods [15][16][17][18]20] to an observed bar code signal, then applies a deblurring method such as the one described in [7].We believe that our proposed method has advantages over such a hypothetical method.Firstly, unlike the proposed method, the existing  nonuniform illumination correction methods do not incorporate the fact that a bar code signal is a bi-level waveform, although a comparative study is deferred to future research.Secondly, since the proposed method jointly performs the estimation of illumination and deblurring of a bar code, it overcomes the problem that errors in nonuniform illumination estimation lead to in errors in deblurring.The hypothetical two-step method involves an optimization of parameters in illumination model followed by a separate optimization of parameters involved in deblurring.Apparently, iterations of these two stages are required for convergence [31].This requires significantly more computation.Moreover, for many cases, it has been reported that such alternating optimization is less effective than simultaneous optimization [31].
We have designed a regularization method based on the assumption that illumination is slowly varying in the spatial domain.By incorporating a regularizing penalty function that quadratically penalizes the changes of neighboring B-spline coefficients, we have shown that the penalized nonlinear least squares problem for the proposed method is well-posed.The argument is based on the assumption that the objective function of the nonlinear least squares with uniform illumination model is locally convex, so the Hessian of the objective function is locally positive definite.This assumption is based on the well-posedness proof of the constrained nonlinear least squares with uniform illumination model in the previous investigation [3].One may conceive of an illumination model using fewer numbers of parameters, such as a third order polynomial.For such a model, regularization might not be necessary.Investigating the advantages and disadvantages of such a model compared with the proposed B-spline based model is deferred to future research.As shown in the simulations and experiments, the design of the regularization parameter λ is important for accurate estimation of illumination and the locations of edges.Using a very large λ in the proposed method yields estimation results similar to the NLS-UI method.Using a very small λ yields unrealistic estimation of illumination and inaccurate estimation of edges due to noise and numerical instability.However, for a broad range of λ values, the proposed method was able to correctly decode bar code signals that were not decodable by conventional methods, suggesting that the tuning of λ may not be very difficult in practice.In our simulations and experiments, we determined λ manually.We plan to investigate an automatic method for tuning of the regularization parameter.
We applied the Gauss-Newton method with active set strategy for solving the constrained penalized nonlinear least squares in our method.Although we designed a constrained optimization method to accommodate constraints such as the positivity of σ , the positivity of illumination, the positivity of background intensity and the orders of edge locations, no constraint was violated during iterations for most bar code images that were successfully decoded.It may be possible to solve the nonlinear least squares problem using other optimization methods such as the variable metric method, the Levenberg-Marquadt method, etc [31].However, it is unlikely that  there will be great improvement using an alternative method since the Gauss-Newton method converged quickly to the minimum.By virtue of the fast convergence, we believe our proposed method is very promising for the real time decoding of bar codes.Moreover, for many bar code images, we think that it might be possible to decode bar codes correctly after only two or three iterations.

Conclusion
We have proposed a novel joint nonuniform illumination estimation and deblurring method for bar code signals based on a penalized nonlinear squares objective function.Based on the proper parameterization of a bar code signal and illumination as well as a regularization on the illumination using a smoothness penalty, the proposed method is able to simultaneously estimate the bar code signal and illumination in the spatial domain.In the simulations and experiments, the proposed method showed improved performance compared with the edge based method without deblurring and the NLS-UI method without nonuniform illumination correction.In particular, the proposed method was able to decode test bar code images in a few iterations that were not decodable by the two methods.Since blurring and nonuniform illumination are unavoidable in many situations, we believe that the proposed method will be useful for decoding bar code images.

Fig. 2 .
Fig. 2. Synthetic blurred bar code signals with nonuniform illuminations: (a) synthetic blurred bar code signal; (b) synthetic nonuniform illuminations; (c) synthetic observation with the Gaussian illumination; (d) synthetic observation with the linear illumination.

Fig. 3 .
Fig. 3. Fitted signals for the observation in Fig. 2(c) using the three methods, plus the estimated illumination by the proposed method: (a) fitted bar code signal using the edge method; (b) fitted bar code signal using the NLS-UI method; (c) fitted bar code signal using the proposed method; (d) estimated illumination by the proposed method.

Fig. 4 .
Fig. 4. Estimated lengths of patterns for the signal shown in Fig. 2(c) using the three methods: (a) the edge method; (b) the NLS-UI method; (c) the proposed method.

Fig. 5 .
Fig. 5.Estimated nonuniform illuminations for the signals shown in Fig. 2(c) and Fig. 2(d) with different λ values: (a) for the signal shown in Fig. 2(c); (b) for the signal shown in Fig. 2(d).

Fig. 6 .
Fig. 6.Change of the objective functions and ME during iterations for the signals shown in Fig. 2(c) and Fig. 2(d): (a) the objective function for the signal shown in Fig. 2(c); (b) the ME for the signal shown in Fig. 2(c); (c) the objective function for the signal shown in Fig. 2(d); (d) the ME for the signal shown in Fig. 2(d).

Fig. 7 .
Fig. 7. Bar code images under nonuniform illumination: (a) bar code image whose central part is brighter due to nonuniform illumination; (b) bar code image whose left part is brighter due to nonuniform illumination; (c) horizontal line scan of the image in (a); (d) horizontal line scan of the image in (b).

Figure 13 (
a) and Fig. 13(b) show the changes of the objective function for the real signals of Fig. 7(c) and Fig. 7(d) with different λ values.The optimization procedure converged to the

Fig. 8 .
Fig. 8. Fitted signals for the observation in Fig. 7(c)) using the three methods, plus the estimated illumination using the proposed method: (a) the edge method; (b) the NLS-UI method; (c) the proposed method; (d) estimated illumination by the proposed method with λ = 10.

Fig. 9 .
Fig. 9. Fitted signals for the observation in Fig. 7(d) by the three methods, plus the estimated illumination by the proposed method: (a) the edge method; (b) the NLS-UI method; (c) the proposed method; (d) estimated illumination by the proposed method with λ = 10.

Fig. 10 .
Fig. 10.Estimated lengths of patterns by the three methods for the signal in Fig. 7(c): (a) the edge method; (b) the NLS-UI method; (c) the proposed method.

Fig. 11 .
Fig. 11.Estimated lengths of patterns by the three methods for the signal in Fig. 7(d): (a) the edge method; (b) the NLS-UI method; (c) the proposed method.

Table 1 .
Average ME values of the three methods for the signal shown in Fig.2(c) with different size blurring kernels and SNR values (the unit for ME value is 10 −3 ).

Table 2 .
Average ME values of the three methods for the signal shown in Fig.2(d) with different size blurring kernels and SNR values (the unit for ME value is 10 −3 ).

Table 3 .
Decode rate values of the three methods for the signal shown in Fig.2(c) with different size blurring kernels and SNR values (the unit for decode rate is %.)