SpringerOpen Newsletter

Receive periodic news and updates relating to SpringerOpen.

This article is part of the series Advanced Statistical Tools for Enhanced Quality Digital Imaging with Realistic Capture Models.

Open Access Research

Robust flash denoising/deblurring by iterative guided filtering

Hae-Jong Seo1* and Peyman Milanfar2

Author Affiliations

1 Sharp Labs of America, Camas, WA 98683, USA

2 University of California-Santa Cruz, 1156 High Street, Santa Cruz, CA 95064, USA

For all author emails, please log on.

EURASIP Journal on Advances in Signal Processing 2012, 2012:3  doi:10.1186/1687-6180-2012-3


The electronic version of this article is the complete one and can be found online at: http://asp.eurasipjournals.com/content/2012/1/3


Received:23 June 2011
Accepted:6 January 2012
Published:6 January 2012

© 2012 Seo and Milanfar; licensee Springer.

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

A practical problem addressed recently in computational photography is that of producing a good picture of a poorly lit scene. The consensus approach for solving this problem involves capturing two images and merging them. In particular, using a flash produces one (typically high signal-to-noise ratio [SNR]) image and turning off the flash produces a second (typically low SNR) image. In this article, we present a novel approach for merging two such images. Our method is a generalization of the guided filter approach of He et al., significantly improving its performance. In particular, we analyze the spectral behavior of the guided filter kernel using a matrix formulation, and introduce a novel iterative application of the guided filter. These iterations consist of two parts: a nonlinear anisotropic diffusion of the noisier image, and a nonlinear reaction-diffusion (residual) iteration of the less noisy one. The results of these two processes are combined in an unsupervised manner. We demonstrate that the proposed approach outperforms state-of-the-art methods for both flash/no-flash denoising, and deblurring.

1 Introduction

Recently, several techniques [1-5] to enhance the quality of flash/no-flash image pairs have been proposed. While the flash image is better exposed, the lighting is not soft, and generally results in specularities and unnatural appearance. Meanwhile, the no-flash image tends to have a relatively low signal-to-noise ratio (SNR) while containing the natural ambient lighting of the scene. The key idea of flash/no-flash photography is to create a new image that is closest to the look of the real scene by having details of the flash image while maintaining the ambient illumination of the no-flash image. Eisemann and Durand [3] used bilateral filtering [6] to give the flash image the ambient tones from the no-flash image. On the other hand, Petschnigg et al. [2] focused on reducing noise in the no-flash image and transferring details from the flash image to the no-flash image by applying joint (or cross) bilateral filtering [3]. Agrawal et al. [4] removed flash artifacts, but did not test their method on no-flash images containing severe noise. As opposed to a visible flash used by [2-4], recently Krishnan and Fergus [7] used both near-infrared and near-ultraviolet illumination for low light image enhancement. Their so-called "dark flash" provides high-frequency detail in a less intrusive way than a visible flash does even though it results in incomplete color information. All these methods ignored any motion blur by either depending on a tripod setting or choosing sufficiently fast shutter speed. However, in practice, the captured images under low-light conditions using a hand-held camera often suffer from motion blur caused by camera shake.

More recently, Zhuo et al. [5] proposed a flash deblurring method that recovers a sharp image by combining a blurry image and a corresponding flash image. They integrated a so-called flash gradient into a maximum-a-posteriori framework and solved the optimization problem by alternating between blur kernel estimation and sharp image reconstruction. This method outperformed many states-of-the-art single image deblurring [8-10] and color transfer methods [11]. However, the final output of this method looks somewhat blurry because the model only deals with a spatially invariant motion blur.

Others have used multiple pictures of a scene taken at different exposures to generate high dynamic range images. This is called multi-exposure image fusion [12] which shares some similarity with our problem in that it seeks a new image that is of better quality than any of the input images. However, the flash/no-flash photography is generally more difficult due to the fact that there are only a pair of images. Enhancing a low SNR no-flash image with a spatially variant motion blur only with the help of a single flash image is still a challenging open problem.

2 Overview of the proposed approach

We address the problem of generating a high quality image from two captured images: a flash image (Z ) and a no-flash image (Y ; Figure 1). We treat these two images, Z and Y , as random variables. The task at hand is to generate a new image (X) that contains the ambient lighting of the no-flash image (Y ) and preserves the details of the flash-image (Z ). As in [2], the new image X can be decomposed into two layers: a base layer and a detail layer;

thumbnailFigure 1. Flash/no-flash pairs. No-flash image can be noisy or blurry.

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M1','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M1">View MathML</a>

(1)

Here, Y might be noisy or blurry (possibly both), and <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M2','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M2">View MathML</a> is an estimated version of Y, enhanced with the help of Z. Meanwhile, <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M3','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M3">View MathML</a> represents a nonlinear, (low-pass) filtered version of Z so that <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M4','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M4">View MathML</a> can provide details. Note that τ is a constant that strikes a balance between the two parts. In order to estimate <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M2','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M2">View MathML</a> and <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M3','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M3">View MathML</a>, we employ local linear minimum mean square error (LMMSE) predictorsa which explain, justify, and generalize the idea of guided filteringb as proposed in [1]. More specifically, we assumed that <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M2','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M2">View MathML</a> and <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M3','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M3">View MathML</a> are a liner (affine) function of Z in a window ωk centered at the pixel k:

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M5','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M5">View MathML</a>

(2)

where G(·) is the guided filtering (LMMSE) operator, <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M6','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M6">View MathML</a>, <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M7','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M7">View MathML</a>, zi are samples of <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M2','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M2">View MathML</a>, <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M3','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M3">View MathML</a>,Z respectively, at pixel i, and (a, b, c, d) are coefficients assumed to be constant in ωk (a square window of size p × p) and space-variant. Once we estimate a, b, c, d, Equation 1 can be rewritten as

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M8','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M8">View MathML</a>

(3)

In fact, <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M9','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M9">View MathML</a> is a linear function of zi. While it is not possible to estimate α and β directly from (equation (3); since they in turn depend on xi), the coefficients α, β can be expressed in terms of a, b, c, d which are optimally estimated from two different local linear models shown in Equation 2. Naturally, the simple linear model has its limitations in capturing complex behavior. Hence, we propose an iterative approach to boost its performance as follows:

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M10','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M10">View MathML</a>

(4)

where <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M11','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M11">View MathML</a> and αn, βn, and τn are functions of the iteration number n. A block-diagram of our approach is shown in Figure 2. The proposed method effectively removes noise and deals well with spatially variant motion blur without the need to estimate any blur kernel or to accurately register flash/no-flash image pairs when there is a modest displacement between them.

thumbnailFigure 2. System overview. Overview of our algorithm for flash/no-flash enhancement.

A preliminary version [13] of this article is appeared in the IEEE International Conference on Computer Vision (ICCV '11) workshop. This article is different from [13] in the following respects:

(1) We have provided a significantly expanded statistical derivation and description of the guided filter and its properties in Section 3 and Appendix.

(2) Figures 3 and 4 are provided to support the key idea of iterative guided filtering.

thumbnailFigure 3. Iteration improves accuracy. Accuracy of <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M67','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M67">View MathML</a> is improved upon <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M68','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M68">View MathML</a>.

thumbnailFigure 4. Effect of iteration. As opposed to <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M69','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M69">View MathML</a>, noise in Y was more effectively removed in <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M70','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M70">View MathML</a>.

(3) We provide many more experimental results for both flash/no-flash denoising and de- blurring in Section 5.

(4) We describe the key ideas of diffusion and residual iteration and their novel relevance to iterative guided filtering in the Appendix.

(5) We prove the convergence of the proposed iterative estimator in the Appendix.

(6) As supplemental material, we share our project websitec where flash/no-flash relighting examples are also presented.

In Section 3, we outline the guided filter and study its statistical properties. We describe how we actually estimate the linear model coefficients a, b, c, d and α, β, and we provide an interpretation of the proposed iterative framework in matrix form in Section 4. In Section 5, we demonstrate the performance of the system with some experimental results, and finally we conclude the article in Section 6.

3 The guided filter and its properties

In general, space-variant, nonparametric filters such as the bilateral filter [6], nonlocal means filter [14], and locally adaptive regression kernels filter [15] are estimated from the given corrupted input image to perform denoising. The guided filter can be distinguished from these in the sense that the filter kernel weights are computed from a (second) guide image which is presumably cleaner. In other words, the idea is to apply filter kernels Wij computed from the guide (e.g., flash) image Z to the more noisy (e.g., no-flash) image Y. Specifically, the filter output sample <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M74','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M74">View MathML</a> a pixel i is computed as a weighted averaged:

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M12','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M12">View MathML</a>

(5)

Note that the filter kernel Wij is a function of the guide image Z, but is independent of Y. The guided filter kernele can be explicitly written as

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M13','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M13">View MathML</a>

(6)

where |ω| is the total number of pixels (= p2) in ωk, ε is a global smoothing parameter, <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M14','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M14">View MathML</a>, and <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M15','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M15">View MathML</a>. Note that Wij are normalized weights, that is, ∑jWij(Z) = 1 Figure 5 shows examples of guided filter weights in four different patches. We can see that the guided filter kernel weights neatly capture underlying geometric structures as do other data-adaptive kernel weights [6,14,15]. It is worth noting that the use of the specific form of the guided filter here may not be critical in the sense that any other data-adaptive kernel weights such as non-local means kernels [16] and locally adaptive regression kernels [15] could be used.

thumbnailFigure 5. Examples of guided filter kernel weights in four different patches. The kernel weights well represent underlying structures.

Next, we study some fundamental properties of the guided filter kernel in matrix form.

We adopt a convenient vector form of Equation 5 as follows:

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M16','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M16">View MathML</a>

(7)

where y is a column vector of pixels in Y and <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M17','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M17">View MathML</a> is a vector of weights for each i. Note that N is the dimensionf of y. Writing the above at once for all i we have,

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M18','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M18">View MathML</a>

(8)

where z is a vector of pixels in Z and W is only a function of z. The filter output can be analyzed as the product of a matrix of weights W with the vector of the given the input image y.

The matrix W is symmetric as shown in Equation 8 and the sum of each row of W is equal to one (W1N = 1N ) by definition. However, as seen in Equation 6, the definition of the weights does not necessarily imply that the elements of the matrix W are positive in general. While this is not necessarily a problem in practice, we find it useful for our purposes to approximate this kernel with a proper admissible kernel [17]. That is, for the purposes of analysis, we approximate W as a positive valued, symmetric positive definite matrix with rows summing to one, as similarly done in [18]. For the details, we refer the reader to the Appendix A.

With this technical approximation in place, all eigenvalues λi (i = 1, ..., N) are real, and the largest eigenvalue of W is exactly one (λ1 = 1), with corresponding eigenvector <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M19','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M19">View MathML</a> as shown in Figure 6. Intuitively, this means that filtering by W will leave a constant signal (i.e., a "flat" image) unchanged. In fact, with the rest of its spectrum inside the unit disk, powers of W converge to a matrix of rank one, with identical rows, which (still) sum to one:

thumbnailFigure 6. Examples of W in one patch of size 25 × 25. All eigenvalues of the matrix W are nonnegative, thus the guided filter kernel matrix is positive definite matrix. The largest eigenvalue of W is one and the rank of W asymptotically becomes one. This figure is better viewed in color.

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M20','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M20">View MathML</a>

(9)

So u1 summarizes the asymptotic effect of applying the filter W many times. Figure 7 shows what a typical u1 looks like.

thumbnailFigure 7. Examples of the first left eigenvector u in three patches. The vector was reshaped into an image for illustration purpose.

Figure 8 shows examples of the (center) row vector (wT) from W's powers in one patch of size 25 × 25. The vector was reshaped into an image for illustration purposes. We can see that powers of W provide even better structure by generating larger (and more sophisticated) kernels. This insight reveals that applying W multiple times can improve the guided filtering performance, which leads us to the iterative use of the guided filter. This approach will produce the evolving coefficients αn, βn introduced in (4). In the following section, we describe how we actually compute these coefficients based on Bayesian mean square error (MSE) predictions.

thumbnailFigure 8. The guided filter kernel matrix. The guided filter kernel matrix W captures the underlying data structure, but powers of W provides even better structure by generating larger (but more sophisticated) kernel shapes. w is the (center) row vector of W. w was reshaped into an image for illustration purposes.

4 Iterative application of local LMMSE predictors

The coefficientsg ak, bk, ck, dk in (3) are chosen so that "on average" the estimated value <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M2','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M2">View MathML</a> is close to the observed value of Y (= yi) in ωk, and the estimated value <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M3','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M3">View MathML</a> is close to the observed value of Z (=zi) in ωk. More specifically, we adopt a stabilized MSE criterion in the window ωk as our measure of closenessh:

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M21','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M21">View MathML</a>

(10)

where ε1 and ε2 are small constants that prevent <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M22','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M22">View MathML</a>, <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M23','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M23">View MathML</a> from being too large. Note that ck and dk become simply 1 and 0 by setting ε2 = 0. By setting partial derivatives of MSE(ak, bk) with respect to ak, bk, and partial derivatives of MSE(ck, dk) with respect to ck, dk, respectively, to zero, the solutions to minimum MSE prediction in (10) are

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M24','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M24">View MathML</a>

(11)

where we compute <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M25','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M25">View MathML</a>.

Note that the use of different ωk results in different predictions of these coefficients. Hence, one must compute an aggregate estimate of these coefficients coming from all windows that contain the pixel of interest. As an illustration, consider a case where we predict <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M6','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M6">View MathML</a> using observed values of Y in ωk of size 3 × 3 as shown in Figure 9. There are nine possible windows that involve the pixel of interest i. Therefore, one takes into account all nine ak, bk's to predict <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M6','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M6">View MathML</a>. The simple strategy suggested by He et al. [1] is to average them as follows:

thumbnailFigure 9. Explanation of LMMSE. <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M71','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M71">View MathML</a> are estimated from nine different windows ωk and averaged coefficients <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M72','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M72">View MathML</a> are used to predict <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M73','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M73">View MathML</a>. This figure is better viewed in color.

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M26','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M26">View MathML</a>

(12)

As such, the resulting prediction of <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M2','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M2">View MathML</a> given the outcome Z = zi is

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M27','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M27">View MathML</a>

(13)

The idea of using these averaged coefficients <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M28','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M28">View MathML</a>, <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M29','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M29">View MathML</a> is analogous to the simplest form of aggregating multiple local estimates from overlapped patches in image denoising and super-resolution literature [19]. The aggregation helps the filter output look locally smooth and contain fewer artifacts.i Recall that <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M6','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M6">View MathML</a> and <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M30','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M30">View MathML</a> correspond to the base layer and the detail layer, respectively. The effect of the regularization parameters ε1 and ε2 is quite the opposite in each case in the sense that the higher ε2 is, the more detail through <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M30','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M30">View MathML</a> can be obtained; whereas the lower ε1 ensures that the image content in <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M2','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M2">View MathML</a> is not over-smoothed.

These local linear models work well when the window size p is small and the underlying data have a simple pattern. However, the linear models are too simple to deal effectively with more complicated structures, and thus there is a need to use larger window sizes. As we alluded to earlier, the estimation of these linear coefficients in an iterative fashion can deal well with more complex behavior of the image content. More specifically, by initializing <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M31','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M31">View MathML</a>, Equation 3 can be updated as follows

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M32','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M32">View MathML</a>

(14)

where n is the iteration number and τn >0 is set to be a monotonically decaying functionk of n such that <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M33','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M33">View MathML</a> converges. Figure 3 shows an example to illustrate that the resulting coefficients at the 20th iteration predict the underlying data better than α1, β1 do. Similarly, <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M34','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M34">View MathML</a> improves upon <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M35','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M35">View MathML</a> as shown in Figure 4. This iteration is closely related to diffusion and residual iteration which are two important methods [18] which we describe briefly below, and with more detail in Appendix.

Recall that Equation 14 can also be written in matrix form as done in Section 3:

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M36','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M36">View MathML</a>

(15)

where W and Wd are guided filter kernel matrices composed of the guided filter kernels W and Wd respectively.l Explicitly writing the iterations, we observe

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M37','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M37">View MathML</a>

(16)

where Pn is a polynomial function of W. The block-diagram in Figure 2 can be redrawn in terms of the matrix formulation as shown in Figure 10. The first term <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M38','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M38">View MathML</a> in Equation 16 is called the diffusion process that enhances SNR. The net effect of each application of W is essentially a step of anisotropic diffusion [20]. Note that this diffusion is applied to the no-flash image y which has a low SNR. On the other hand, the second term <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M39','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M39">View MathML</a> is connected with the idea of residual iteration [21]. The key idea behind this iteration is to filter the residual signalsm to extract detail. We refer the reader to Appendix B and [18] for more detail. By effectively combining the diffusion and residual iteration [as in (16)], we can achieve the goal of flash/no-flash pair enhancement which is to generate an image somewhere between the flash image z and the no-flash image y, but of better quality than both.n

thumbnailFigure 10. Diagram of the proposed iterative approach in matrix form. Note that the iteration can be divided into two parts: diffusion and residual iteration process.

5 Experimental results

In this section, we apply the proposed approach to flash/no-flash image pairs for denoising and deblurring. We convert images Z and Y from RGB color space to CIE Lab, and perform iterative guided filtering separately in each resulting channel. The final result is converted back to RGB space for display. We used the implementation of the guided filter [1] from the author's website.o All figures in this section are best viewed in color.p

5.1 Flash/no-flash denoising

5.1.1 Visible flash [2]

We show experimental results on a couple of flash/no-flash image pairs where no-flash images suffer from noiseq. We compare our results with the method based on the joint bilateral filter [2] in Figures 11, 12, and 13. Our proposed method effectively denoised the no-flash image while transferring the fine detail of the flash image and maintaining the ambient lighting of the no-flash image. We point out that the proposed iterative application of the guided filtering in terms of diffusion and residual iteration yielded much better results than one application of either the joint bilateral filtering [2] or the guided filter [1].

thumbnailFigure 11. Flash/no-flash denoising example compared to the state of the art method [2]. The iteration n for this example is 10.

thumbnailFigure 12. Flash/no-flash denoising example compared to the state of the art method [2]. The iteration n for this example is 10.

thumbnailFigure 13. Flash/no-flash denoising example compared to the state of the art method [2]. The iteration n for this example is 2.

5.1.2 Dark flash [7]

In this section, we use the dark flash method proposed in [7]. Let us call the dark flash image Z. Dark flash may introduce shadows and specularities in images, which affect the results of both the denoising and detail transfer. We detect those regions using the same methods proposed by [2]. Shadows are detected by finding the regions where |Z - Y| is small, and specularities are found by detecting saturated pixels in Z. After combining the shadow and specularities mask, we blur it using a Gaussian filter to feather the boundaries. By using the resulting mask, the output <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M40','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M40">View MathML</a> at each iteration is alpha-blended with a low-pass filter version of Y as similarly done in [2,7]. In order to realize ambient lighting conditions, we applied the same mapping function to the final output as in [7]. Figures 14, 15, 16, 17, 18, and 19 show that our results yield better detail with less color artifacts than the results of [7].

thumbnailFigure 14. Dark flash/no-flash (low noise) denoising example compared to the state of the art method [7]. The iteration n for this example is 10.

thumbnailFigure 15. Dark flash/no-flash (mid noise) denoising example compared to the state of the art method [7]. The iteration n for this example is 10.

thumbnailFigure 16. Dark flash/no-flash (high noise) denoising example compared to the state of the art method [7]. The iteration n for this example is 10.

thumbnailFigure 17. Dark flash/no-flash (low noise) denoising example compared to the state of the art method [7]. The iteration n for this example is 10.

thumbnailFigure 18. Dark flash/no-flash (mid noise) denoising example compared to the state of the art method [7]. The iteration n for this example is 10.

thumbnailFigure 19. Dark flash/no-flash (high noise) denoising example compared to the state-of-the-art method [7]. The iteration n for this example is 10.

5.2 Flash/no-flash deblurring

Motion blur due to camera shake is an annoying yet common problem in low-light photography. Our proposed method can also be applied to flash/no-flash deblurringr. Here, we show experimental results on a couple of flash/no-flash image pairs where no-flash images suffer from mild noise and strong motion blur. We compare our method with Zhuo et al. [5]. As shown in Figures 20, 21, 22, 23, and 24, our method outperforms the method by [5], obtaining much finer details with better color contrast even though our method does not estimate a blur kernel at all. The results by Zhuo et al. [5] tend to be somewhat blurry and distort the ambient lighting of the real scene. We point out that we only use a single blurred image in Figure 24 while Zhuo et al. [5] used two blurred images and one flash image.

thumbnailFigure 20. Flash/no-flash deblurring example compared to the state-of-the-art method [5]. The iteration n for this example is 20.

thumbnailFigure 21. Flash/no-flash deblurring example compared to the state-of-the-art method [5]. The iteration n for this example is 20.

thumbnailFigure 22. Flash/no-flash deblurring example compared to the state-of-the-art method [5]. The iteration n for this example is 20.

thumbnailFigure 23. Flash/no-flash deblurring example compared to the state-of-the-art method [5]. The iteration n for this example is 5.

thumbnailFigure 24. Flash/no-flash deblurring example compared to the state-of-the-art method [5]. The iteration n for this example is 20.

6 Summary and future work

The guided filter has proved to be more effective than the joint bilateral filter in several applications. Yet we have shown that it can be improved significantly more still. We analyzed the spectral behavior of the guided filter kernel using a matrix formulation and improved its performance by applying it iteratively. Iterations of the proposed method consist of a combination of diffusion and residual iteration. We demonstrated that the proposed approach yields outputs that not only preserve fine details of the flash image, but also the ambient lighting of the no-flash image. The proposed method outperforms state-of-the-art methods for flash/no-flash image denoising and deblurring. It would be interesting to see if the performance of other nonparametric filer kernels such as bilateral filters and locally adaptive regression kernels [15] can be further improved in our iterative framework. It is also worthwhile to explore several other applications such as joint upsampling [22], image matting [23], mesh smoothing [24,25], and specular highlight removal [26] where the proposed approach can be employed.

Appendix

Positive definite and symmetric row-stochastic approximation of W

In this section, we describe how we approximate W with a symmetric, positive definite, and row-stochastic matrix. First, as we mentioned earlier, the matrix W can contain negative values as shown in Figure 3. We employ the Taylor series approximation (exp (t) ≈ 1 + t) to ensure that W has both positive elements, and is positive-definite. To be more concrete, consider a simple example; namely, a local patch of size 5 × 5 as follows:

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M41','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M41">View MathML</a>

(17)

In this case, we can have up to 9 ωk of size |ω| = 3 × 3. Let <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M42','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M42">View MathML</a> in Equation 6. Then, W centered at the index 13 can be written and approximated as follows:

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M43','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M43">View MathML</a>

Next, we convert the matrix W (composed of strictly positive elements now) to a doubly-stochastic, symmetric, positive definite matrix as again done in [18]. The algorithm we use to effect this approximation is due to Sinkhorn [27,28], who proved that given a matrix with strictly positive elements, there exist diagonal matrices R = diag(r) and C = diag(c) such that

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M44','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M44">View MathML</a>

is doubly stochastic. That is,

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M45','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M45">View MathML</a>

(18)

Furthermore, the vectors r and c are unique to within a scalar (i.e., α r, c.) Sinkhorn's algorithm for obtaining r and c in effect involves repeated normalization of the rows and columns (see Algorithm 1 for details) so that they sum to one, and is provably convergent and optimal in the cross-entropy sense [29].

Algorithm 1 Algorithm for scaling a matrix A to a nearby doubly-stochastic matrix <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M46','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M46">View MathML</a>

   Given a matrix A, let (N, N) be size(A) and initialize r = ones(N, 1);

   for k = 1 : iter;

      c = 1./(AT r);

      r = 1./(A c);

   end

   C = diag(c); R = diag(r);

   <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M47','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M47">View MathML</a>

Diffusion and residual iteration

Diffusion

Here, we describe how multiple direct applications of the filter given by W is in effect equivalent to a nonlinear anisotropic diffusion process [20,30]. We define <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M48','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M48">View MathML</a>, and

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M49','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M49">View MathML</a>

(19)

From the iteration (19) we have

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M50','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M50">View MathML</a>

(20)

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M51','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M51">View MathML</a>

(21)

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M52','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M52">View MathML</a>

(22)

which we can rewrite as

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M53','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M53">View MathML</a>

(23)

where <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M54','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M54">View MathML</a> is a scaled version of by, and therefore the left-hand side of the above is a discretization of the derivative operator <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M55','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M55">View MathML</a>, and as detailed in [18], W - I is effectively the nonlinear Laplacian operator corresponding to the kernel in (6).

Residual iteration

An alternative to repeated applications of the filter W is to consider the residual signals, defined as the difference between the estimated signal and the measured signal. This results in a variation of the diffusion estimator which uses the residuals as an additional forcing term. The net result is a type of reaction-diffusion process [31]. In statistics, the use of the residuals in improving estimates has a rather long history, dating at least back to the study of Tukey [21] who termed the idea "twicing". More recently, the idea has been suggested in the applied mathematics community under the rubric of Bregman iterations [32], and in the machine learning and statistics literature [33] as L2-boosting.

Formally, the residuals are defined as the difference between the estimated signal and the measured signal: <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M56','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M56">View MathML</a>, where here we define the initializations <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M57','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M57">View MathML</a>. With this definition, we write the iterated estimates as

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M58','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M58">View MathML</a>

(24)

Explicitly writing the iterations, we observe:

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M59','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M59">View MathML</a>

(25)

where Fn is a polynomial function of W of order n + 1. The first iterate <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M60','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M60">View MathML</a> is precisely the "twicing" estimate of Tukey [21].

Convergence of the proposed iterative estimator

Recall the iterations:

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M61','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M61">View MathML</a>

(26)

where Pn is a polynomial function of W. The first (direct diffusion) term in the iteration is clearly stable and convergent as described in (9). Hence, we need to show the stability of the second part. To do this, we note that the spectrum of Pn can be written, and bounded in terms of the eigenvalues λi of W as follows:

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M62','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M62">View MathML</a>

(27)

where the inequality follows from the knowledge that 0 ≤ λN ≤ ... λ3 λ2 <λ1 = 1. Furthermore, in Section 4 we defined τn to be a monotonically decreasing sequence such that <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M63','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M63">View MathML</a>. Hence, all eigenvalues λi(Pn) are upper bounded by the constant c, independent of the number of iterations n, ensuring the stability of the iterative process.

End notes

aMore detail is provided in Section 4. bThe guided filter [1] reduces noise while preserving edges as bilateral filter [6] does. However, the guided filter outperforms the bilateral filter by avoiding the gradient reversal artifacts that may appear in such applications as detail enhancement, high dynamic range (HDR) compression, and flash/no-flash denoising. chttp://users.soe.ucsc.edu/~milanfar/research/rokaf/.html/IGF/ webcite. d <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M64','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M64">View MathML</a> in Equation 2 can be rewritten in terms of filtering and we refer the reader to a supplemental material http://personal.ie.cuhk.edu.hk/~hkm007/eccv10/eccv10supp.pdf webcite for derivation. eCross (or joint) bilateral filter [2,3] is defined in a similar way. f N is different from the window size p (N p). g Note that k is used to clarify that the coefficients are estimated for the window ωk. h Shan et al. [8] recently proposed a similar approach for high dynamic range compression. i It is worthwhile to note that we can benefit from more adaptive way of combining multiple estimates of coefficients, but this subject is not treated in this article. k We use <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M65','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M65">View MathML</a> throughout the all experiments. l Recall W in Equation 6. The difference between W and Wd lies in the parameter ε as follows (ε2 > ε1):

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M66','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M66">View MathML</a>

(28)

m This is generally defined as the difference between the estimated signal <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/3/mathml/M3','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/3/mathml/M3">View MathML</a> and the measured signal Z, but in our context refers to the detail signal. n We refer the reader to Appendix C for proof of convergence of the proposed iterative estimator. o http://personal.ie.cuhk.edu.hk/~hkm007/ webcite. p We refer the reader to the project Website http://users.soe.ucsc.edu/~milanfar/research/rokaf/.html/IGF/ webcite. q The window size p for Wd and W was set to 21 and 5 respectively for all the denoising examples. r The window size p for Wd and W was set to 41 and 81, respectively, for the deblurring examples to deal with displacement between the flash and no-flash image. s Note that due to the use of residuals, this is a different initialization than the one used in the diffusion iterations.

Competing interests

The authors declare that they have no competing interests.

Authors' contributions

HS carried out the design of iterative guided filtering and drafted the manuscript. PM participated in the design of iterative guided filtering and performed the statistical analysis. All authors read and approved the final manuscript.

Acknowledgements

We thank Dilip Krishnan for sharing the post-processing code [7] for dark-flash examples. This study was supported by the AFOSR Grant FA 9550-07-01-0365 and NSF Grant CCF-1016018. This study was done while the first author was at the University of California.

References

  1. K He, J Sun, X Tang, Guided image filtering. Proceedings of European Conference Computer Vision (ECCV) (2010)

  2. G Petschnigg, M Agrawala, H Hoppe, R Szeliski, M Cohen, K Toyama, Digital photography with flash and no-flash image pairs. ACM Trans Graph 21(3), 664–672 (2004)

  3. E Eisemann, F Durand, Flash photography enhancement via intrinsic relighting. ACM Trans Graph 21(3), 673–678 (2004)

  4. A Agrawal, R Raskar, S Nayar, Y Li, Removing photography artifacts using gradient projection and flash-exposure sampling. ACM Trans Graph 24, 828–835 (2005). Publisher Full Text OpenURL

  5. S Zhuo, D Guo, T Sim, Robust flash deblurring. IEEE Conference on Computer Vison and Pattern Recognition (2010)

  6. C Tomasi, R Manduchi, Bilateral Filtering for Gray and Color Images. Proceedings of the 1998 IEEE International Conference of Compute Vision Bombay, India, 836–846 (1998)

  7. D Krishnan, R Fergus, Dark flash photography. ACM Trans Graph 28(4), 594–611 (2009)

  8. Q Shan, J Jia, MS Brown, Globally optimized linear windowed tone-mapping. IEEE Trans Vis Comput Graph 16(4), 663–675 (2010). PubMed Abstract | Publisher Full Text OpenURL

  9. R Fergus, B Singh, A Hertsmann, ST Roweis, WT Freeman, Removing camera shake from a single image. ACM Trans Graph (SIGGRAPH) 25, 787–794 (2006). Publisher Full Text OpenURL

  10. L Yuan, J Sun, L Quan, HY Shum, Progressive inter-scale and intra-scale non-blind image deconvolution. ACM Trans Graph 27(3), 1–10 (2008)

  11. YW Tai, J Jia, CK Tang, Local color transfer via probabilistic segmentation by expectation-maximization. IEEE Conference on Computer Vison and Pattern Recognition (2005)

  12. W Hasinoff, Variable-aperture photography (PhD Thesis, Department of Computer Science, University of Toronto, 2008)

  13. H Seo, P Milanfar, Computational photography using a pair of flash/no-flash images by iterative guided filtering. IEEE International Conference on Computer Vision (ICCV) (2011) Submitted OpenURL

  14. A Buades, B Coll, JM Morel, A review of image denoising algorithms, with a new one. Multi-scale Model Simulat (SIAM) 4(2), 490–530 (2005)

  15. H Takeda, S Farsiu, P Milanfar, Kernel regression for image processing and reconstruction. IEEE Trans Image Process 16(2), 349–366 (2007). PubMed Abstract OpenURL

  16. A Buades, B Coll, JM Morel, Nonlocal image and movie denoising. Int J Comput Vis 76(2), 123–139 (2008). Publisher Full Text OpenURL

  17. T Hofmann, B Scholkopf, AJ Smola, Kernel methods in machine learning. Ann Stat 36(3), 1171–1220 (2008). Publisher Full Text OpenURL

  18. P Milanfar, A tour of modern image processing.. IEEE Signal Process Mag (2011) [http:/ / users.soe.ucsc.edu/ ~milanfar/ publications/ journal/ ModernTour_FinalSubmission.pdf] webcite OpenURL

  19. M Protter, M Elad, H Takeda, P Milanfar, Generalizing the non-local-means to super-resolution reconstruction. IEEE Trans Image Process 18, 36–51 (2009). PubMed Abstract | Publisher Full Text OpenURL

  20. P Perona, J Malik, Scale-space and edge detection using anistropic diffusion. IEEE Trans Pattern Anal Mach Intell 12(9), 629–639 (1990)

  21. JW Tukey, Exploratory Data Analysis (Addison Wesley, Reading, MA, 1977)

  22. J Kopf, MF Cohen, D Lischinski, M Uyttendaele, Joint bilateral upsampling. ACM Trans Graph 26(3), 96 (2007). Publisher Full Text OpenURL

  23. K He, J Sun, X Tang, Fast matting using large kernel matting Laplacian matrices. IEEE Conference on Computer Vision and Pattern Recognition (CVPR) (2010)

  24. S Fleishman, I Drori, D Cohen-Or, Bilateral mesh denoising. ACM Trans Graph 22(3), 950–953 (2003). Publisher Full Text OpenURL

  25. T Jones, F Durand, M Desbrun, Non-iterative feature preserving mesh smoothing. ACM Trans Graph 22(3), 943–949 (2003). Publisher Full Text OpenURL

  26. Q Yang, S Wang, N Ahuja, Real-time specular highlight removal using bilateral filtering. ECCV (2010)

  27. PA Knight, The Sinkhorn-Knopp algorithm: convergence and applications. SIAM J Matrix Anal Appl 30, 261–275 (2008). Publisher Full Text OpenURL

  28. R Sinkhorn, A relationship between arbitrary positive matrices and doubly stochastic matrices. Ann Math Stat 35(2), 876–879 (1964). Publisher Full Text OpenURL

  29. J Darroch, D Ratcliff, Generalized iterative scaling for log-linear models. Ann Math Stat 43, 1470–1480 (1972). Publisher Full Text OpenURL

  30. A Singer, Y Shkolinsky, B Nadler, Diffusion interpretation of nonlocal neighborhood filters for signal denoising. SIAM J Imaging Sci 2, 118–139 (2009). Publisher Full Text OpenURL

  31. G Cottet, L Germain, Image processing through reaction combined with nonlinear diffusion. Math Comput 61, 659–673 (1993). Publisher Full Text OpenURL

  32. S Osher, M Burger, D Goldfarb, J Xu, W Yin, An iterative regularization method for total variation-based image restoration. Multiscale Model Simulat 4(2), 460–489 (2005). Publisher Full Text OpenURL

  33. P Buhlmann, B Yu, Boosting with the L2 loss: regression and classification. J Am Stat Assoc 98(462), 324–339 (2003). Publisher Full Text OpenURL