Retinex by Two Bilateral Filters - Semantic Scholar
Retinex by Two Bilateral Filters Michael Elad The Computer Science Department The Technion - Israel Institute of Technology, Haifa 32000 Israel. Email: [email protected] Abstract. Retinex theory deals with the removal of unfavorable illumination effects from images. This ill-posed inverse problem is typically regularized by forcing spatial smoothness on the recoverable illumination. Recent work in this field suggested exploiting the knowledge that the illumination image bounds the image from above, and the fact that the reflectance is also expected to be smooth. In this paper we show how the above model can be improved to provide a non-iterative retinex algorithm that handles better edges in the illumination, and suppresses noise in dark areas. This algorithm uses two specially tailored bilateral filters – the first evaluates the illumination and the other is used for the computation of the reflectance. This result stands as a theoretic justification and refinement for the recently proposed heuristic use of the bilateral filter for retinex by Durand and Dorsey. In line with their appealing way of speeding up the bilateral filter, we show that similar speedup methods apply to our algorithm.
1
Introduction
Retinex theory deals with the removal of unfavorable illumination effects from a given image. A commonly assumed model suggests that any given image S is the pixel-wise multiplication of two images, the reflectance R and the illumination L, i.e., S = R · L. A look-up-table log operation transfers this multiplication into an addition, resulting with s = log(S) = log(L) + log(R) = ` + r. Clearly, the recovery of ` from s is an ill-posed inverse problem. Solving it is typically done by introducing a regularization that forces a spatial smoothness on the recoverable illumination. Thus, early heuristic and successful retinex methods, such as the homomorphic filtering algorithm [1] and many others (e.g., [2–5]), proposed a low-pass filter on s, or an algorithm that amounts to this effect, to obtain a rough estimate of `. In this paper we refrain from reviewing this literature, and limit our approach to this topic by building upon a recent study presented in [6]. The work in [6] describes several improvements to the classical retinex models. One improvement refers to the passivity of the reflectance, assumed to satisfy 0 ≤ R ≤ 1. As a direct consequence we have that L ≥ S, implying that the illumination image should be an envelope image bounding S from above. Due to the monotonicity of the log operation we have ` ≥ s. Merging the above with the desire to get spatially smooth ` may lead to the trivial and meaningless result of a constant image, ` = max(s). The remedy, as proposed in [6], is to assume that
2
Michael Elad
k` − sk2 should be small, implying that ` should upper envelope s while being close to it. Based on these modifications, the reconstruction of the illumination can be posed as the following quadratic programming (QP) problem o n 2 2 2 (1) min λ k` − sk2 + kDx `k2 + kDy `k2 . `≥s
The operators Dx and Dy represent horizontal and vertical discrete derivatives, forcing this way spatial smoothness. A second ingredient introduced in [6, 7] is a smoothness penalty forced also on the reflectance image r = s − `. This added to (1) gives n o n o 2 2 2 2 2 min λ k` − sk2 + kDx `k2 + kDy `k2 + α kDx (s − `)k2 + kDy (s − `)k2 .(2) `≥s
Note that, since s = ` + r is enforced, the new term contradicts the illumination smoothness, as r and ` cannot be jointly smooth. Thus, the effect is to gain some smoothness in r at the expense of losing some of it in `. The justification for this is the desire to lead r to be “nice-looking”, as natural images should be. Based on the above model, an efficient multi-scale algorithm has been proposed in [6] to estimate ` and thus r. The work in [7] used the same model to propose a simplified estimate solvers based on known implementation constraints. More recently, [8] further simplified the computation of ` by introducing a spatial recursive smoothing filter. While the above model is general enough and covers the correct forces to be used in the solution of the retinex problem, it has several flaws: – Hallows: A commonly encountered artifact with retinex algorithms is the existence of hallows. This is a direct consequence of the smoothness assumption discussed above. When passing from a strongly illuminated region to a dark zone (e.g., on a border of a shadowed area), the smoothness forces the illumination to remain high in the dark region near this edge and smoothly descend to grasp the illumination within the dark region. Thus, when removed, the dark regions near such edges remain dark, resulting with these hallow effects. Such effects can be also obtained in the bright areas near such illumination edges, if the constraint ` ≥ s is not practiced. Then those bright areas become further brighter. – Noise: In dark regions of the image these retinex algorithms are expected to yield a contrast stretching, very much similar to the effect caused by standard Gamma correction. The stretching causes a magnification of the noise, and this becomes evident especially in low-quality images, or ones with noticeable compression. The constraint s = ` + r implies that the noise migrates as a whole to the two ingredients, rather than being suppressed. – Iterative Solution: The above model formulation leads naturally to the need for an iterative solver. The work in [7] and [8] bypassed this limitation, but with a price on the final outcome’s quality. In this paper we propose an alternative model for retinex, and a numerical algorithm that builds on it. The new model is similar to the one in (2), in the
Lecture Notes in Computer Science
3
sense that all the presented forces are included. However, this new model is enhanced to solve the above mentioned shortcomings. More specifically, smoothness of the illumination and the reflectance are forced using a robust statistics method, and hallows are avoided. The smoothness terms used are very much in the spirit of the bilateral filter [9, 10], having a wide stencil effect that enables avoiding the need for an iterative or multi-scale solver. We use different smoothing formulation for the reflectance and the illumination to handle them differently, and absorb the constraint ` ≥ s in a natural way. This leads to a two-stage algorithm that applies two variations of the bilateral filter, first estimating the illumination, and then the reflectance. The new model suppresses noise by allowing ` + r to deviate from s, implying that the residual should be the additive noise we want to discard of. The new model and accompanying algorithm stand as a theoretic justification and refinement for the recently proposed heuristic use of the bilateral filter for retinex as appeared in [11]. In line with their appealing way of speeding up the bilateral filter, we show similar speedup methods for our two bilateral filter variations. This paper is organized as follows: Section 2 presents the bilateral filter that this work is building on. Section 3 then turn to describe the new model for retinex, and the algorithm that emerges from it. Speedup methods are discussed in Section 4. Section 5 presents some results, and Section 6 concludes this paper.
2
Denoising by the Bilateral Filter
In this Section we present the bilateral filter, designed for the removal of additive noise from images [9]. We also discuss its origins as described in [10–13]. These will serve us as we turn later to consider the retinex problem. Consider an image s contaminated by additive noise. Our goal is to develop an edge-preserving smoothing algorithm that effectively removes most of the noise while preserving the image details. A maximum a-posteriori probability (MAP) formulation of this problem as presented in [10] yields 2
min λ kˆ s − sk2 + sˆ
P X
P X
T
(Cm,n sˆ − sˆ) W[m,n] (s) (Cm,n sˆ − sˆ) .
(3)
m=−P n=−P
The operators Cm,n are shift operators, moving the image sˆ by m pixels horizontally and n pixels vertically. The matrices W[m,n] are diagonal matrices that down-weight large edge entries in s so as not to smooth over edges of the image. The choice W[m,n] (s) = I ∀ m, n leads to the non-robust option that makes the overall problem QP as in (2). Choosing these weights to be inversely proportional to |Cm,n s − s| leads to the ability to handle edges in the image better. Note that using weighting here parallels the use of robust statistics - more on this relationship and can be found in [10, 12]. The fact that smoothness is forced in a wide neighborhood implies that even a simple iteration to minimize this functional will be very effective. Indeed, the work in [10] established that the bilateral filter as presented by [9] is an
4
Michael Elad
approximate solver of this programming task. More specifically, it was shown that the bilateral filter amounts to a single Jacobi iteration over this penalty term. Here we briefly show this property and its meaning. The Jacobi step is constructed using the gradient and the diagonal of the Hessian of the penalty function in (3). The gradient is given by P X ∂F {ˆ s} = 2λ(ˆ s − s) + 2 ∂ˆ s
P X
T
(Cm,n − I) W[m,n] (s) (Cm,n − I) sˆ. (4)
m=−P n=−P
The Hessian of F is given by P X ∂ 2 F {ˆ s} = 2λI + 2 ∂ˆ s2
P X
T
(Cm,n − I) W[m,n] (s) (Cm,n − I) .
(5)
m=−P n=−P
Denoting the main diagonal of the Hessian as the matrix 0.5M(s)1 , and assuming an initialization sˆ = s, the first Jacobi iteration to minimize F gives ( )−1 ¯ ¯ ∂ 2 F {ˆ s} ¯¯ ∂F {ˆ s} ¯¯ sˆ1 = sˆ0 − diag · (6) ∂ˆ s2 ¯sˆ0 =s ∂ˆ s ¯sˆ=s " # P P X X T −1 = I − M(s) · (Cm,n − I) W[m,n] (s) (Cm,n − I) s. m=−P n=−P
The above represents an operator that multiplies the image s. This operator applies a weighted sum of the input pixels in a stencil of (2P + 1)-by-(2P + 1) pixels to compute the output, and these weights are dependent on W[m,n] (s) and the local differences between the center pixel s[k, j] and its neighbors s[k − m, j − n]. Thus, this is a spatially adaptive FIR filter of some sort. In [10] it was shown that if the [m, n] weight at the pixel [k, j] is chosen as W[m,n] (k, j) =
ρ0 {s[k, j] − s[k − m, j − n]} · V [m, n], s[k, j] − s[k − m, j − n]
(7)
then we obtain the very filter that Tomasi and Manduchi proposed in [9]. For this equivalence we have to choose λ = 1, ρ(x) = 1−exp(−x2 /2σ 2 ), and V [m, n] being a Gaussian kernel. Still, we can consider many other robust functions and weights V [m, n] that give a filter very much in line with the spirit of the bilateral filter. Interestingly, this filter is a discrete version of the short-time effective kernel of the Beltrami flow as discussed in [14, 15]. This implies that this algorithm has deep roots in the geometric understanding of images as manifolds. While the above analysis is helpful in understanding the origins of the bilateral filter, it is hard to understand how it is applied in practice. As shown in [9], the effective filter computes every output pixel sˆ1 [k, j] by sˆ1 [k, j] =
P X
P X
a[m, n, k, j]s[k − m, j − n],
m=−P n=−P 1
The additional 0.5 comes to null the factor 2 in the gradient term.
(8)
Lecture Notes in Computer Science
where
a[m, n, k, j] =
³ 2 2 exp − m2σ+n − 2 s
(s[k,j]−s[k−m,j−n])2 2σr2
Z[k, j]
5
´ .
(9)
The term Z[k, j] normalizes these weights to sum to one. This filter assigns per every neighbor a weight inversely proportional to its Euclidean distance (m2 +n2 ) and inversely proportional to its distance in gray-value from the center pixel. The parameters σr and σs governs the behavior of the filters - more on those can be found in [9, 10].
3
Retinex by Two Bilateral Filters
In this section we present the new model for the retinex problem that uses the bilateral smoothness term. We use this model to develop the two bilateral filters that compose our novel retinex algorithm. Hallows in the retinex result could be avoided by allowing ` to be piece-wise 2 smooth. This could be easily accomplished by replacing the terms kDx `k2 + 2 kDy `k2 with kDx `k1 + kDy `k1 , TV [16], or any other robust statistics based penalty, and there are numerous options of the like. However, adopting such local terms implies a need for many iterations in the numerical solution. Thus, we consider instead the bilateral smoothness. For brevity of notations, we denote hereafter BW,P {x} =
P X
P X
T
(Cm,n x − x) W[m,n] (s) (Cm,n x − x) .
(10)
m=−P n=−P
Starting from the quadratic programming problem posed in (2), we propose the following alternative model for retinex n o n o 2 2 min λ` k` − sk2 + BW` ,P` {`} + α λr kr − s + `k2 + BWr ,Pr {r} .(11) `, r: `≥s
The first part handles the smoothness of the illumination ` and its proximity to s, while bounding it from above. The second part introduces the smoothness of the reflectance r, and requires it to be close to the residual image s − `. Thus, noise can be discarded by becoming the residual s − ` − r. Note also that our notations hints to the fact that we will consider different weights and parameters in the smoothness terms for ` and r. The formulation given in (11) leads to a decomposition that seeks both ` and r as unknown, and one does not imply the other as before. Instead of optimizing with respect to both ` and r in parallel (which is an option we have not explored in our work, but one that can certainly be addressed based on the model we have posed), we adopt a two stage process, first estimating `, based on the first part in (11), and then given `, we evaluate r. Starting with the quest for `, let us attempt to evaluate it such that it addresses only the first term in (11). Thus, we seek a solution to the problem min
`: `≥s
2
λ` · k` − sk2 + BW` ,P` {`} .
(12)
6
Michael Elad
Clearly, without the constraint ` ≥ s, the above is equivalent to the problem posed in (3), and as such, the bilateral filter is an excellent solver candidate. Thus, the natural question we should pose here is how the constraint should be accommodated, in a way that preserves the convenience of the bilateral filter. We propose to introduce a special choice of weights W` that handle the constraint implicitly. These new weights are based on Equation (7), but using a one-sided robust function ρ, ½ 1 − exp(−x2 /2σr2 ) x ≤ 0 ρ(x) = . (13) ∞ x>0 This alternative choice of weights introduces a simple modification to the bilateral filter, where, among the (2P + 1)-by-(2P + 1) neighbors per each pixel, we consider only those that satisfy s[k, j] ≤ s[k − m, j − n]. This way, the local averaging is done with non-negative normalized weights, while combining only pixels that have higher gray values than the center pixel, resulting with a final outcome that must satisfy `1 [k, j] ≥ s[k, j]. Thus, this new filter will necessarily achieve both a satisfaction of the constraint (by virtue of the weights), while reducing the newly defined penalty term that still considers smoothness as we desire. We refer hereafter to this filter as the envelope-bilateral filter. In practice, the above implies that the bilateral filter as presented in section 2 is slightly changed. Parallel to (8) and (9), in the envelope-bilateral filter every output pixel `1 [k, j] is evaluated by `1 [k, j] =
P X
P X
a[m, n, k, j]s[k − m, j − n],
(14)
m=−P n=−P
where a[m, n, k, j] =
³ 2 2 exp − m2σ+n − 2 s
(s[k,j]−s[k−m,j−n])2 2σr2
´ · µ{s[k − m, j − n] − s[k, j]}
Z[k, j]
The notation µ{x} stands for the step-function, being 1 for non-negative x and zero elsewhere. The term Z[k, j] normalizes these weights to sum to one, as before. Note that from the above description it is clear that if s[k, j] is the peak of its (2P + 1)2 neighborhood, then its filtering amounts to `1 [k, j] = s[k, j], since in this case all weights are zero and only a[0, 0, k, j] = 1. Assuming that the above stage has been completed, we have an estimate of ` and we now turn to evaluate r. We consider the second term in (11), solving 2
min λr · kr − (s − `)k2 + BWr ,Pr {r} . r
(15)
Since the image s − ` is given, this is the very bilateral filter formulation in (2). Thus, an application of the bilateral filter on the image s − ` should lead to the desired r. However, due to the transform to the log-domain, the noise that should be discarded from the reflectance image resides mostly in the regions where s
.
Lecture Notes in Computer Science
7
is low. Thus, we can better direct the above bilateral filter by using σr to be inversely proportional to s to reflect this matter. A choice of the form σr [k, j] = (C1 · s[k, j]p + C2 )−1 could be used to this effect2 . Nothing in the definition or the implementation of the bilateral filter prevents having such spatially adaptive parameter. This will ensure that r is hardly smoothed in regions where s is bright, while it is being smoothed in darker regions.
4
Speeding Up the Retinex Algorithm
In their paper, Durand and Dorsey proposed a wonderful speedup algorithm for the bilateral filter, and this algorithm can be applied directly to both our two bilateral filter versions. Here we outline the basic ideas of this speedup, starting from Equation (14), although everything said applies just as well to the second bilateral filter. Referring to s[k, j] in these equations as a constant c, we can re-write these equations as ¶ µ P P X X 1 m2 + n2 · · (16) exp − Z[k, j] 2σs2 m=−P n=−P · µ ¶ ¸ (c − s[k − m, j − n])2 · exp − · µ{s[k − m, j − n] − c}s[k − m, j − n] 2σr2 ¶ µ P P X X 1 m2 + n2 = · · g[k − m, j − n]. exp − Z[k, j] 2σs2
`1 [k, j] =
m=−P n=−P
This expression is a convolution between the image g[k, j], being · ½ ¾ ¸ (c − s[k, j])2 g[k, j] = exp − · µ{s[k, j] − c}s[k, j] , 2σr2
(17)
and the Gaussian blur. Thus, we could apply a sequence of such convolutions, scanning the values of s[k, j] in the range [0, loge 255], and then merging the results, choosing the proper values from each output, based on the s[k, j] values. Note that the value of Z[k, j] is given by a similar expression Z[k, j] =
P X
P X
m=−P n=−P
µ ¶ m 2 + n2 exp − u[k − m, j − n], 2σs2
(18)
where ½
(c − s[k, j])2 u[k, j] = exp − 2σr2 2
¾ · µ{s[k, j] − c}.
(19)
Recall that S[k, j] ∈ [1, 255] as we shift by 1 to avoid singularities, and we have also 0 ≤ s[k, j] ≤ 5.54.
8
Michael Elad
Thus, its computation can also be done using a sequence of similar convolutions. Durand and Dorsey proposed two ways to further speed-up the evaluation of `1 : (i) piece-wise linear approximation; and (ii) multi-scale implementation. The first idea is to scan the values of s[k, j] in the range [0, loge 255] with jumps, and interpolate in between. Practically speaking, using 30−50 equispaced jumps in the range [0, loge 255] are found to induce almost no change to the outcome. Since our weights include a step-function discontinuity, the interpolation should be done as a one-sided operation, always preferring to adopt the larger c to avoid a violation of the ` ≥ s constraint. This causes the interpolation to loose some of of its accuracy, but our experiments show that this lose is mild and unnoticeable. As to the multi-scale option, since images are convolved in the above expressions with wide-range Gaussian smoothers, a pre down-scale and post up-scale yield a substantial gain in run-time with almost no change in the outcome. The gain is especially noticed for wide supports (P À 1, and σs À 1). On top of these two ideas, note that the required convolutions required are all separable. Furthermore, when σs is large enough, the effective convolving kernel is the square step function. In such a case further speedup can be obtained using the computation of the integral image [17].
5
Results
An interesting idea reported in [6] is to return some of the illumination to the reflectance when presenting the final output image. Thus, the output image is computed as Out = R[k, j] · L[k, j]1/γ = S[k, j] · L[k, j]1/γ−1 . Reflectance images are typically unrealistic looking, and with a modest and reduced effect of illumination returned to it, the final image enjoys both the desired brightness and the natural appearance. We have made use of this idea in the following presented results. The illumination is returned to the original image by applying Gamma-correction on it using γ = 3, and multiplying it back by the estimated reflectance. Figures 1-2 present two pairs of original images3 and their retinex results. In these two cases, the use of spatially varying σr has very little effect because the images are of high quality, and thus we do not show it. Figure 3 presents results for a third image, where the dark region is noisy, and thus the two versions are shown side-by-side for comparison (with parameters p = 8, C1 = 5e − 3, C2 = 0.3). We should note that in processing color images we apply the retinex algorithm to the luminance (V) layer in the HSV color representation, and leave the chromatic layers unchanged. In all three cases we used the following setup parameters: the envelope bilateral filter parameters used are P` = 15, σr = 0.3, and σs = 100. The second bilateral filter used Pr = 4, σr = 0.3 or an adaptive method as described earlier, and σs = 100. The speedup algorithm was used with scale down factor of 2 : 1, and grey-value steps of 0.1. Figures 4 and 5 return to the first example, presenting several accompanying results. Figure 4 shows the obtained reflectance and illumination results (gray3
These images and the one in Figure 3 are from the NASA retinex web-page.
Lecture Notes in Computer Science
9
Fig. 1. Example 1 - An original image (left) and its retinex result (right).
Fig. 2. Example 2 - An original image (left) and its retinex result (right).
Fig. 3. Example 3 - An original image (left) and its retinex results using a regular bilateral filter for computing r (middle) and using the spatially adaptive σr (right).
10
Michael Elad
value images referring to the V-layer). As can be seen, what we call ‘reflectance’ is far from being satisfactory to describe the image, and indeed there is room to return of illumination. This Figure indicates that our separation is not perfect and there is a leakage between r and `. In fact, our ` are r stand for large scale intensity components and small scale corrections, respectively, both estimated with preservation of discontinuities. Still, the final outcome is satisfactory because of the illumination return. Figure 5 describes the Out-to-In correspondence of the overall retinex algorithm, showing that while the retinex process generally resembles a Gammacorrection effect, is has a different effect of varying Gamma. This idea is further expanded, showing in Figure 5 the image γeffective = log In/ log Out as a function of the location. This gives the effective Gamma correction that should be applied in every pixel to reproduce the obtained result.
Fig. 4. The reflectance and the illumination images in Example 1.
3
2.8
2.6
2.4
2.2
2
1.8
1.6
1.4
1.2
1
Fig. 5. Example 1 - Right: The Input-Output mapping, and overlayed on it is the Gamma-correction that corresponds to γ = 3; Left: The effective Gamma correction value per pixel.
The overall improvement in speed introduced by the speedup algorithm depends on many of the parameters that are mentioned above, and on implemen-
Lecture Notes in Computer Science
11
tation issues. We compared two efficient implementations of the bilateral filter - both implemented with Matlab. The first sweeps through the support of the filter, applying operations on complete images, and the other being the speedup algorithm mentioned above. For the parameters used here we obtained a factor of 5 − 10 shorter run time with the speedup algorithm. Figure 6 presents a comparison between the new algorithm and the one reported in [6] on a severely degraded image4 . For this comparison we changed the color space to YCbCr, and choose γ = 2.3 in the illumination return, both done to match with the alternative algorithm. The results show strong hallows in the previous method, while those are fully suppressed by our algorithm.
Fig. 6. An original image (top), the new retinex algorithm (bottom left), and the one reported in[6] (bottom right).
6
Conclusion
In this paper we have presented a new model for the retinex problem – removal of undesired illumination effects from an image. The new model enables a better handling of edges in the illumination that causes hallow effects, and it enables the 4
Curtesy of Eyal Gordon, The CS department - The Technion.
12
Michael Elad
suppression of noise in dark areas. An algorithm based on this model has been developed, leading to two specially tailored bilateral filters, the first evaluates the illumination and the second is used for the computation of the reflectance. Our work stands as a theoretic justification and refinement for the recently proposed heuristic use of the bilateral filter for retinex by Durand and Dorsey. We have used their way of speeding up the bilateral to propose a similar speedup methods for our filters.
References 1. Faugeras, O.D. (1979) “Digital image color processing within the framework of a human visual system”, IEEE Trans.on ASSP, 27: 380–393. 2. Land, E.H. (1983) “Recent advances in the Retinex theory and some implications for cortical computations: Color vision and the natural image”, Proc. Nat. Acad. Sci. USA, 80: 5163–5169. 3. Blake, A. (1985)“Boundary conditions of lightness computation in mondrian world”, Computer Vision Graphics and Image Processing, 32: 314–327. 4. Terzopoulos, D. (1986) “Image analysis using multigrid relaxation methods”, IEEE Trans. on PAMI, 8: 129–139. 5. Funt, B., Ciurea, F., and McCann, J. (2004) Retinex in Matlab, Journal of the Electronic Imaging, 13(1):48-57. Terzopoulos, D. (1986) “Image analysis using multigrid relaxation methods”, IEEE Trans. on PAMI, 8: 129–139. 6. Kimmel, R., Elad, M., Shaked, D., Keshet, R., and Sobel. I. (2003) “A variational framework for retinex”, International Journal of Computer Vision 52(1): 7–23. 7. M. Elad, M., Kimmel, R., Shaked, D., and Keshet, R. (2003) “Reduced complexity retinex algorithm via the variational approach”, Journal on Visual Communication and Image Representation, 14(4): 369–388. 8. Shaked, D. and Keshet, R. (2002) “Robust recursive envelope operators for fast retinex”, Hewlett-Packard Research Laboratories Technical Report, HPL-2002-74R1. 9. Tomasi, C. and Manduchi, R. (1998) “Bilateral filtering for gray and color images”, Proc. 6th Int. Conf. Computer Vision, New Delhi,India, 839–846. 10. Elad, M. (2002) “On the bilateral filter and ways to improve it”, IEEE Trans. On Image Processing, 11(10): 1141–1151. 11. Durand, F. and Dorsey, J. (2002) “Fast bilateral filtering for the display of highdynamic-range images”, SIGGRAPH 2002: 257–266. 12. Barash, D. (2001) “A fundamental relationship between bilateral filtering, adaptive smoothing, and the nonlinear diffusion equation”, IEEE Trans. Pattern Anal. Mach. Intell., 24(6): 844–847. 13. Mrazek, P., Weickert, J., and Bruhn, A. (2005) On robust estimation and smoothing with spatial and tonal kernels. In Geometric Properties from Incomplete Data. Kluwer, Dordrecht, 2005, to appear. 14. Sochen, N., Kimmel, R., and Malladi, R. (1998) “A geometrical framework for low level vision”, IEEE Trans. Image Processing, 7: 310–318. 15. Sochen, N., Kimmel, R. and Bruckstein, A.M. (2001) “Diffusions and confusions in signal and image processing”, J. Math. Imag. Vis., 14(3): 195–209. 16. Chan, T.F., Osher, S., and Shen, J. (2001) “The digital TV filter and nonlinear denoising”, IEEE Trans. Image Processing, 10: 231–241. 17. Viola, P. and Jones, M. (2001) “Rapid object detection using a boosted cascade of simple features”, Conference on Computer Vision and Pattern recognition.
1
Introduction
Retinex theory deals with the removal of unfavorable illumination effects from a given image. A commonly assumed model suggests that any given image S is the pixel-wise multiplication of two images, the reflectance R and the illumination L, i.e., S = R · L. A look-up-table log operation transfers this multiplication into an addition, resulting with s = log(S) = log(L) + log(R) = ` + r. Clearly, the recovery of ` from s is an ill-posed inverse problem. Solving it is typically done by introducing a regularization that forces a spatial smoothness on the recoverable illumination. Thus, early heuristic and successful retinex methods, such as the homomorphic filtering algorithm [1] and many others (e.g., [2–5]), proposed a low-pass filter on s, or an algorithm that amounts to this effect, to obtain a rough estimate of `. In this paper we refrain from reviewing this literature, and limit our approach to this topic by building upon a recent study presented in [6]. The work in [6] describes several improvements to the classical retinex models. One improvement refers to the passivity of the reflectance, assumed to satisfy 0 ≤ R ≤ 1. As a direct consequence we have that L ≥ S, implying that the illumination image should be an envelope image bounding S from above. Due to the monotonicity of the log operation we have ` ≥ s. Merging the above with the desire to get spatially smooth ` may lead to the trivial and meaningless result of a constant image, ` = max(s). The remedy, as proposed in [6], is to assume that
2
Michael Elad
k` − sk2 should be small, implying that ` should upper envelope s while being close to it. Based on these modifications, the reconstruction of the illumination can be posed as the following quadratic programming (QP) problem o n 2 2 2 (1) min λ k` − sk2 + kDx `k2 + kDy `k2 . `≥s
The operators Dx and Dy represent horizontal and vertical discrete derivatives, forcing this way spatial smoothness. A second ingredient introduced in [6, 7] is a smoothness penalty forced also on the reflectance image r = s − `. This added to (1) gives n o n o 2 2 2 2 2 min λ k` − sk2 + kDx `k2 + kDy `k2 + α kDx (s − `)k2 + kDy (s − `)k2 .(2) `≥s
Note that, since s = ` + r is enforced, the new term contradicts the illumination smoothness, as r and ` cannot be jointly smooth. Thus, the effect is to gain some smoothness in r at the expense of losing some of it in `. The justification for this is the desire to lead r to be “nice-looking”, as natural images should be. Based on the above model, an efficient multi-scale algorithm has been proposed in [6] to estimate ` and thus r. The work in [7] used the same model to propose a simplified estimate solvers based on known implementation constraints. More recently, [8] further simplified the computation of ` by introducing a spatial recursive smoothing filter. While the above model is general enough and covers the correct forces to be used in the solution of the retinex problem, it has several flaws: – Hallows: A commonly encountered artifact with retinex algorithms is the existence of hallows. This is a direct consequence of the smoothness assumption discussed above. When passing from a strongly illuminated region to a dark zone (e.g., on a border of a shadowed area), the smoothness forces the illumination to remain high in the dark region near this edge and smoothly descend to grasp the illumination within the dark region. Thus, when removed, the dark regions near such edges remain dark, resulting with these hallow effects. Such effects can be also obtained in the bright areas near such illumination edges, if the constraint ` ≥ s is not practiced. Then those bright areas become further brighter. – Noise: In dark regions of the image these retinex algorithms are expected to yield a contrast stretching, very much similar to the effect caused by standard Gamma correction. The stretching causes a magnification of the noise, and this becomes evident especially in low-quality images, or ones with noticeable compression. The constraint s = ` + r implies that the noise migrates as a whole to the two ingredients, rather than being suppressed. – Iterative Solution: The above model formulation leads naturally to the need for an iterative solver. The work in [7] and [8] bypassed this limitation, but with a price on the final outcome’s quality. In this paper we propose an alternative model for retinex, and a numerical algorithm that builds on it. The new model is similar to the one in (2), in the
Lecture Notes in Computer Science
3
sense that all the presented forces are included. However, this new model is enhanced to solve the above mentioned shortcomings. More specifically, smoothness of the illumination and the reflectance are forced using a robust statistics method, and hallows are avoided. The smoothness terms used are very much in the spirit of the bilateral filter [9, 10], having a wide stencil effect that enables avoiding the need for an iterative or multi-scale solver. We use different smoothing formulation for the reflectance and the illumination to handle them differently, and absorb the constraint ` ≥ s in a natural way. This leads to a two-stage algorithm that applies two variations of the bilateral filter, first estimating the illumination, and then the reflectance. The new model suppresses noise by allowing ` + r to deviate from s, implying that the residual should be the additive noise we want to discard of. The new model and accompanying algorithm stand as a theoretic justification and refinement for the recently proposed heuristic use of the bilateral filter for retinex as appeared in [11]. In line with their appealing way of speeding up the bilateral filter, we show similar speedup methods for our two bilateral filter variations. This paper is organized as follows: Section 2 presents the bilateral filter that this work is building on. Section 3 then turn to describe the new model for retinex, and the algorithm that emerges from it. Speedup methods are discussed in Section 4. Section 5 presents some results, and Section 6 concludes this paper.
2
Denoising by the Bilateral Filter
In this Section we present the bilateral filter, designed for the removal of additive noise from images [9]. We also discuss its origins as described in [10–13]. These will serve us as we turn later to consider the retinex problem. Consider an image s contaminated by additive noise. Our goal is to develop an edge-preserving smoothing algorithm that effectively removes most of the noise while preserving the image details. A maximum a-posteriori probability (MAP) formulation of this problem as presented in [10] yields 2
min λ kˆ s − sk2 + sˆ
P X
P X
T
(Cm,n sˆ − sˆ) W[m,n] (s) (Cm,n sˆ − sˆ) .
(3)
m=−P n=−P
The operators Cm,n are shift operators, moving the image sˆ by m pixels horizontally and n pixels vertically. The matrices W[m,n] are diagonal matrices that down-weight large edge entries in s so as not to smooth over edges of the image. The choice W[m,n] (s) = I ∀ m, n leads to the non-robust option that makes the overall problem QP as in (2). Choosing these weights to be inversely proportional to |Cm,n s − s| leads to the ability to handle edges in the image better. Note that using weighting here parallels the use of robust statistics - more on this relationship and can be found in [10, 12]. The fact that smoothness is forced in a wide neighborhood implies that even a simple iteration to minimize this functional will be very effective. Indeed, the work in [10] established that the bilateral filter as presented by [9] is an
4
Michael Elad
approximate solver of this programming task. More specifically, it was shown that the bilateral filter amounts to a single Jacobi iteration over this penalty term. Here we briefly show this property and its meaning. The Jacobi step is constructed using the gradient and the diagonal of the Hessian of the penalty function in (3). The gradient is given by P X ∂F {ˆ s} = 2λ(ˆ s − s) + 2 ∂ˆ s
P X
T
(Cm,n − I) W[m,n] (s) (Cm,n − I) sˆ. (4)
m=−P n=−P
The Hessian of F is given by P X ∂ 2 F {ˆ s} = 2λI + 2 ∂ˆ s2
P X
T
(Cm,n − I) W[m,n] (s) (Cm,n − I) .
(5)
m=−P n=−P
Denoting the main diagonal of the Hessian as the matrix 0.5M(s)1 , and assuming an initialization sˆ = s, the first Jacobi iteration to minimize F gives ( )−1 ¯ ¯ ∂ 2 F {ˆ s} ¯¯ ∂F {ˆ s} ¯¯ sˆ1 = sˆ0 − diag · (6) ∂ˆ s2 ¯sˆ0 =s ∂ˆ s ¯sˆ=s " # P P X X T −1 = I − M(s) · (Cm,n − I) W[m,n] (s) (Cm,n − I) s. m=−P n=−P
The above represents an operator that multiplies the image s. This operator applies a weighted sum of the input pixels in a stencil of (2P + 1)-by-(2P + 1) pixels to compute the output, and these weights are dependent on W[m,n] (s) and the local differences between the center pixel s[k, j] and its neighbors s[k − m, j − n]. Thus, this is a spatially adaptive FIR filter of some sort. In [10] it was shown that if the [m, n] weight at the pixel [k, j] is chosen as W[m,n] (k, j) =
ρ0 {s[k, j] − s[k − m, j − n]} · V [m, n], s[k, j] − s[k − m, j − n]
(7)
then we obtain the very filter that Tomasi and Manduchi proposed in [9]. For this equivalence we have to choose λ = 1, ρ(x) = 1−exp(−x2 /2σ 2 ), and V [m, n] being a Gaussian kernel. Still, we can consider many other robust functions and weights V [m, n] that give a filter very much in line with the spirit of the bilateral filter. Interestingly, this filter is a discrete version of the short-time effective kernel of the Beltrami flow as discussed in [14, 15]. This implies that this algorithm has deep roots in the geometric understanding of images as manifolds. While the above analysis is helpful in understanding the origins of the bilateral filter, it is hard to understand how it is applied in practice. As shown in [9], the effective filter computes every output pixel sˆ1 [k, j] by sˆ1 [k, j] =
P X
P X
a[m, n, k, j]s[k − m, j − n],
m=−P n=−P 1
The additional 0.5 comes to null the factor 2 in the gradient term.
(8)
Lecture Notes in Computer Science
where
a[m, n, k, j] =
³ 2 2 exp − m2σ+n − 2 s
(s[k,j]−s[k−m,j−n])2 2σr2
Z[k, j]
5
´ .
(9)
The term Z[k, j] normalizes these weights to sum to one. This filter assigns per every neighbor a weight inversely proportional to its Euclidean distance (m2 +n2 ) and inversely proportional to its distance in gray-value from the center pixel. The parameters σr and σs governs the behavior of the filters - more on those can be found in [9, 10].
3
Retinex by Two Bilateral Filters
In this section we present the new model for the retinex problem that uses the bilateral smoothness term. We use this model to develop the two bilateral filters that compose our novel retinex algorithm. Hallows in the retinex result could be avoided by allowing ` to be piece-wise 2 smooth. This could be easily accomplished by replacing the terms kDx `k2 + 2 kDy `k2 with kDx `k1 + kDy `k1 , TV [16], or any other robust statistics based penalty, and there are numerous options of the like. However, adopting such local terms implies a need for many iterations in the numerical solution. Thus, we consider instead the bilateral smoothness. For brevity of notations, we denote hereafter BW,P {x} =
P X
P X
T
(Cm,n x − x) W[m,n] (s) (Cm,n x − x) .
(10)
m=−P n=−P
Starting from the quadratic programming problem posed in (2), we propose the following alternative model for retinex n o n o 2 2 min λ` k` − sk2 + BW` ,P` {`} + α λr kr − s + `k2 + BWr ,Pr {r} .(11) `, r: `≥s
The first part handles the smoothness of the illumination ` and its proximity to s, while bounding it from above. The second part introduces the smoothness of the reflectance r, and requires it to be close to the residual image s − `. Thus, noise can be discarded by becoming the residual s − ` − r. Note also that our notations hints to the fact that we will consider different weights and parameters in the smoothness terms for ` and r. The formulation given in (11) leads to a decomposition that seeks both ` and r as unknown, and one does not imply the other as before. Instead of optimizing with respect to both ` and r in parallel (which is an option we have not explored in our work, but one that can certainly be addressed based on the model we have posed), we adopt a two stage process, first estimating `, based on the first part in (11), and then given `, we evaluate r. Starting with the quest for `, let us attempt to evaluate it such that it addresses only the first term in (11). Thus, we seek a solution to the problem min
`: `≥s
2
λ` · k` − sk2 + BW` ,P` {`} .
(12)
6
Michael Elad
Clearly, without the constraint ` ≥ s, the above is equivalent to the problem posed in (3), and as such, the bilateral filter is an excellent solver candidate. Thus, the natural question we should pose here is how the constraint should be accommodated, in a way that preserves the convenience of the bilateral filter. We propose to introduce a special choice of weights W` that handle the constraint implicitly. These new weights are based on Equation (7), but using a one-sided robust function ρ, ½ 1 − exp(−x2 /2σr2 ) x ≤ 0 ρ(x) = . (13) ∞ x>0 This alternative choice of weights introduces a simple modification to the bilateral filter, where, among the (2P + 1)-by-(2P + 1) neighbors per each pixel, we consider only those that satisfy s[k, j] ≤ s[k − m, j − n]. This way, the local averaging is done with non-negative normalized weights, while combining only pixels that have higher gray values than the center pixel, resulting with a final outcome that must satisfy `1 [k, j] ≥ s[k, j]. Thus, this new filter will necessarily achieve both a satisfaction of the constraint (by virtue of the weights), while reducing the newly defined penalty term that still considers smoothness as we desire. We refer hereafter to this filter as the envelope-bilateral filter. In practice, the above implies that the bilateral filter as presented in section 2 is slightly changed. Parallel to (8) and (9), in the envelope-bilateral filter every output pixel `1 [k, j] is evaluated by `1 [k, j] =
P X
P X
a[m, n, k, j]s[k − m, j − n],
(14)
m=−P n=−P
where a[m, n, k, j] =
³ 2 2 exp − m2σ+n − 2 s
(s[k,j]−s[k−m,j−n])2 2σr2
´ · µ{s[k − m, j − n] − s[k, j]}
Z[k, j]
The notation µ{x} stands for the step-function, being 1 for non-negative x and zero elsewhere. The term Z[k, j] normalizes these weights to sum to one, as before. Note that from the above description it is clear that if s[k, j] is the peak of its (2P + 1)2 neighborhood, then its filtering amounts to `1 [k, j] = s[k, j], since in this case all weights are zero and only a[0, 0, k, j] = 1. Assuming that the above stage has been completed, we have an estimate of ` and we now turn to evaluate r. We consider the second term in (11), solving 2
min λr · kr − (s − `)k2 + BWr ,Pr {r} . r
(15)
Since the image s − ` is given, this is the very bilateral filter formulation in (2). Thus, an application of the bilateral filter on the image s − ` should lead to the desired r. However, due to the transform to the log-domain, the noise that should be discarded from the reflectance image resides mostly in the regions where s
.
Lecture Notes in Computer Science
7
is low. Thus, we can better direct the above bilateral filter by using σr to be inversely proportional to s to reflect this matter. A choice of the form σr [k, j] = (C1 · s[k, j]p + C2 )−1 could be used to this effect2 . Nothing in the definition or the implementation of the bilateral filter prevents having such spatially adaptive parameter. This will ensure that r is hardly smoothed in regions where s is bright, while it is being smoothed in darker regions.
4
Speeding Up the Retinex Algorithm
In their paper, Durand and Dorsey proposed a wonderful speedup algorithm for the bilateral filter, and this algorithm can be applied directly to both our two bilateral filter versions. Here we outline the basic ideas of this speedup, starting from Equation (14), although everything said applies just as well to the second bilateral filter. Referring to s[k, j] in these equations as a constant c, we can re-write these equations as ¶ µ P P X X 1 m2 + n2 · · (16) exp − Z[k, j] 2σs2 m=−P n=−P · µ ¶ ¸ (c − s[k − m, j − n])2 · exp − · µ{s[k − m, j − n] − c}s[k − m, j − n] 2σr2 ¶ µ P P X X 1 m2 + n2 = · · g[k − m, j − n]. exp − Z[k, j] 2σs2
`1 [k, j] =
m=−P n=−P
This expression is a convolution between the image g[k, j], being · ½ ¾ ¸ (c − s[k, j])2 g[k, j] = exp − · µ{s[k, j] − c}s[k, j] , 2σr2
(17)
and the Gaussian blur. Thus, we could apply a sequence of such convolutions, scanning the values of s[k, j] in the range [0, loge 255], and then merging the results, choosing the proper values from each output, based on the s[k, j] values. Note that the value of Z[k, j] is given by a similar expression Z[k, j] =
P X
P X
m=−P n=−P
µ ¶ m 2 + n2 exp − u[k − m, j − n], 2σs2
(18)
where ½
(c − s[k, j])2 u[k, j] = exp − 2σr2 2
¾ · µ{s[k, j] − c}.
(19)
Recall that S[k, j] ∈ [1, 255] as we shift by 1 to avoid singularities, and we have also 0 ≤ s[k, j] ≤ 5.54.
8
Michael Elad
Thus, its computation can also be done using a sequence of similar convolutions. Durand and Dorsey proposed two ways to further speed-up the evaluation of `1 : (i) piece-wise linear approximation; and (ii) multi-scale implementation. The first idea is to scan the values of s[k, j] in the range [0, loge 255] with jumps, and interpolate in between. Practically speaking, using 30−50 equispaced jumps in the range [0, loge 255] are found to induce almost no change to the outcome. Since our weights include a step-function discontinuity, the interpolation should be done as a one-sided operation, always preferring to adopt the larger c to avoid a violation of the ` ≥ s constraint. This causes the interpolation to loose some of of its accuracy, but our experiments show that this lose is mild and unnoticeable. As to the multi-scale option, since images are convolved in the above expressions with wide-range Gaussian smoothers, a pre down-scale and post up-scale yield a substantial gain in run-time with almost no change in the outcome. The gain is especially noticed for wide supports (P À 1, and σs À 1). On top of these two ideas, note that the required convolutions required are all separable. Furthermore, when σs is large enough, the effective convolving kernel is the square step function. In such a case further speedup can be obtained using the computation of the integral image [17].
5
Results
An interesting idea reported in [6] is to return some of the illumination to the reflectance when presenting the final output image. Thus, the output image is computed as Out = R[k, j] · L[k, j]1/γ = S[k, j] · L[k, j]1/γ−1 . Reflectance images are typically unrealistic looking, and with a modest and reduced effect of illumination returned to it, the final image enjoys both the desired brightness and the natural appearance. We have made use of this idea in the following presented results. The illumination is returned to the original image by applying Gamma-correction on it using γ = 3, and multiplying it back by the estimated reflectance. Figures 1-2 present two pairs of original images3 and their retinex results. In these two cases, the use of spatially varying σr has very little effect because the images are of high quality, and thus we do not show it. Figure 3 presents results for a third image, where the dark region is noisy, and thus the two versions are shown side-by-side for comparison (with parameters p = 8, C1 = 5e − 3, C2 = 0.3). We should note that in processing color images we apply the retinex algorithm to the luminance (V) layer in the HSV color representation, and leave the chromatic layers unchanged. In all three cases we used the following setup parameters: the envelope bilateral filter parameters used are P` = 15, σr = 0.3, and σs = 100. The second bilateral filter used Pr = 4, σr = 0.3 or an adaptive method as described earlier, and σs = 100. The speedup algorithm was used with scale down factor of 2 : 1, and grey-value steps of 0.1. Figures 4 and 5 return to the first example, presenting several accompanying results. Figure 4 shows the obtained reflectance and illumination results (gray3
These images and the one in Figure 3 are from the NASA retinex web-page.
Lecture Notes in Computer Science
9
Fig. 1. Example 1 - An original image (left) and its retinex result (right).
Fig. 2. Example 2 - An original image (left) and its retinex result (right).
Fig. 3. Example 3 - An original image (left) and its retinex results using a regular bilateral filter for computing r (middle) and using the spatially adaptive σr (right).
10
Michael Elad
value images referring to the V-layer). As can be seen, what we call ‘reflectance’ is far from being satisfactory to describe the image, and indeed there is room to return of illumination. This Figure indicates that our separation is not perfect and there is a leakage between r and `. In fact, our ` are r stand for large scale intensity components and small scale corrections, respectively, both estimated with preservation of discontinuities. Still, the final outcome is satisfactory because of the illumination return. Figure 5 describes the Out-to-In correspondence of the overall retinex algorithm, showing that while the retinex process generally resembles a Gammacorrection effect, is has a different effect of varying Gamma. This idea is further expanded, showing in Figure 5 the image γeffective = log In/ log Out as a function of the location. This gives the effective Gamma correction that should be applied in every pixel to reproduce the obtained result.
Fig. 4. The reflectance and the illumination images in Example 1.
3
2.8
2.6
2.4
2.2
2
1.8
1.6
1.4
1.2
1
Fig. 5. Example 1 - Right: The Input-Output mapping, and overlayed on it is the Gamma-correction that corresponds to γ = 3; Left: The effective Gamma correction value per pixel.
The overall improvement in speed introduced by the speedup algorithm depends on many of the parameters that are mentioned above, and on implemen-
Lecture Notes in Computer Science
11
tation issues. We compared two efficient implementations of the bilateral filter - both implemented with Matlab. The first sweeps through the support of the filter, applying operations on complete images, and the other being the speedup algorithm mentioned above. For the parameters used here we obtained a factor of 5 − 10 shorter run time with the speedup algorithm. Figure 6 presents a comparison between the new algorithm and the one reported in [6] on a severely degraded image4 . For this comparison we changed the color space to YCbCr, and choose γ = 2.3 in the illumination return, both done to match with the alternative algorithm. The results show strong hallows in the previous method, while those are fully suppressed by our algorithm.
Fig. 6. An original image (top), the new retinex algorithm (bottom left), and the one reported in[6] (bottom right).
6
Conclusion
In this paper we have presented a new model for the retinex problem – removal of undesired illumination effects from an image. The new model enables a better handling of edges in the illumination that causes hallow effects, and it enables the 4
Curtesy of Eyal Gordon, The CS department - The Technion.
12
Michael Elad
suppression of noise in dark areas. An algorithm based on this model has been developed, leading to two specially tailored bilateral filters, the first evaluates the illumination and the second is used for the computation of the reflectance. Our work stands as a theoretic justification and refinement for the recently proposed heuristic use of the bilateral filter for retinex by Durand and Dorsey. We have used their way of speeding up the bilateral to propose a similar speedup methods for our filters.
References 1. Faugeras, O.D. (1979) “Digital image color processing within the framework of a human visual system”, IEEE Trans.on ASSP, 27: 380–393. 2. Land, E.H. (1983) “Recent advances in the Retinex theory and some implications for cortical computations: Color vision and the natural image”, Proc. Nat. Acad. Sci. USA, 80: 5163–5169. 3. Blake, A. (1985)“Boundary conditions of lightness computation in mondrian world”, Computer Vision Graphics and Image Processing, 32: 314–327. 4. Terzopoulos, D. (1986) “Image analysis using multigrid relaxation methods”, IEEE Trans. on PAMI, 8: 129–139. 5. Funt, B., Ciurea, F., and McCann, J. (2004) Retinex in Matlab, Journal of the Electronic Imaging, 13(1):48-57. Terzopoulos, D. (1986) “Image analysis using multigrid relaxation methods”, IEEE Trans. on PAMI, 8: 129–139. 6. Kimmel, R., Elad, M., Shaked, D., Keshet, R., and Sobel. I. (2003) “A variational framework for retinex”, International Journal of Computer Vision 52(1): 7–23. 7. M. Elad, M., Kimmel, R., Shaked, D., and Keshet, R. (2003) “Reduced complexity retinex algorithm via the variational approach”, Journal on Visual Communication and Image Representation, 14(4): 369–388. 8. Shaked, D. and Keshet, R. (2002) “Robust recursive envelope operators for fast retinex”, Hewlett-Packard Research Laboratories Technical Report, HPL-2002-74R1. 9. Tomasi, C. and Manduchi, R. (1998) “Bilateral filtering for gray and color images”, Proc. 6th Int. Conf. Computer Vision, New Delhi,India, 839–846. 10. Elad, M. (2002) “On the bilateral filter and ways to improve it”, IEEE Trans. On Image Processing, 11(10): 1141–1151. 11. Durand, F. and Dorsey, J. (2002) “Fast bilateral filtering for the display of highdynamic-range images”, SIGGRAPH 2002: 257–266. 12. Barash, D. (2001) “A fundamental relationship between bilateral filtering, adaptive smoothing, and the nonlinear diffusion equation”, IEEE Trans. Pattern Anal. Mach. Intell., 24(6): 844–847. 13. Mrazek, P., Weickert, J., and Bruhn, A. (2005) On robust estimation and smoothing with spatial and tonal kernels. In Geometric Properties from Incomplete Data. Kluwer, Dordrecht, 2005, to appear. 14. Sochen, N., Kimmel, R., and Malladi, R. (1998) “A geometrical framework for low level vision”, IEEE Trans. Image Processing, 7: 310–318. 15. Sochen, N., Kimmel, R. and Bruckstein, A.M. (2001) “Diffusions and confusions in signal and image processing”, J. Math. Imag. Vis., 14(3): 195–209. 16. Chan, T.F., Osher, S., and Shen, J. (2001) “The digital TV filter and nonlinear denoising”, IEEE Trans. Image Processing, 10: 231–241. 17. Viola, P. and Jones, M. (2001) “Rapid object detection using a boosted cascade of simple features”, Conference on Computer Vision and Pattern recognition.
