Oliver Whyte
email: oliverw@microsoft.com

Josef Sivic
email: josef.sivic@ens.fr

Andrew Zisserman

Deblurring Shaken and Partially Saturated Images

Keywords: Non-Blind Deblurring, Saturation, Ringing, Outliers

We address the problem of deblurring images degraded by camera shake blur and saturated (over-exposed) pixels. Saturated pixels violate the common assumption that the image-formation process is linear, and often cause ringing in deblurred outputs. We provide an analysis of ringing in general, and show that in order to prevent ringing, it is insufficient to simply discard saturated pixels. We show that even when saturated pixels are removed, ringing is caused by attempting to estimate the values of latent pixels that are brighter than the sensor's maximum output. Estimating these latent pixels is likely to cause large errors, and these errors propagate across the rest of the image in the form of ringing. We propose a new deblurring algorithm that locates these error-prone bright pixels in the latent sharp image, and by decoupling them from the remainder of the latent image, greatly reduces ringing. In addition, we propose an approximate forward model for saturated images, which allows us to estimate these error-prone pixels separately without causing artefacts. Results are shown for non-blind deblurring of

Introduction

The task of deblurring "shaken" images has received considerable attention recently [START_REF] Fergus | Removing camera shake from a single photograph[END_REF][START_REF] Cho | Fast motion deblurring[END_REF][START_REF] Gupta | Single image deblurring using motion density functions[END_REF][START_REF] Joshi | Image deblurring using inertial measurement sensors[END_REF][START_REF] Whyte | Non-uniform deblurring for shaken images[END_REF][START_REF] Shan | High-quality motion deblurring from a single image[END_REF]. Significant progress has been made towards reliably estimating the point spread function (PSF) for a given blurred image, which describes how the image was blurred. Likewise, when the PSF for an image is known, many authors have proposed methods to invert the blur process and recover a high quality sharp image (referred to as "non-blind deblurring").

One problematic feature of blurred images, and in particular "shaken" images, which has received relatively little attention is the presence of saturated (over-exposed) pixels. These arise when the radiance of the scene exceeds the range of the camera's sensor, leaving bright highlights clipped at the maximum output value (e.g. 255 for an image with 8 bits per pixel). To anyone who has attempted to take handheld photographs at night, such pixels should be familiar as the conspicuous bright streaks left by electric lights, such as in Figure 1a. These bright pixels, with their clipped values, violate the assumption made by most deblurring algorithms that the image formation process is linear, and as a result can cause obtrusive artefacts in the deblurred images. This can be seen in the deblurred images in Figures 1b and1c.

In this paper we address the problem of deblurring images containing saturated pixels, offering an analysis of the artefacts caused by existing algorithms, and a new algorithm which avoids such artefacts by explicitly handling saturated pixels. Our method is applicable for all causes of blur, however in this work we focus on blur caused by camera shake (motion of the camera during the exposure).

The process of deblurring an image typically involves two steps. First, the PSF is estimated, either using a blind deblurring algorithm [START_REF] Fergus | Removing camera shake from a single photograph[END_REF][START_REF] Yuan | Image deblurring with blurred/noisy image pairs[END_REF][START_REF] Shan | High-quality motion deblurring from a single image[END_REF][START_REF] Cho | Fast motion deblurring[END_REF] to estimate the PSF from the blurred image itself, or by using additional hardware attached to the camera [START_REF] Joshi | Image deblurring using inertial measurement sensors[END_REF][START_REF] Tai | Image/video deblurring using a hybrid camera[END_REF]. Second, a non-blind deblurring algorithm is applied to estimate the sharp image, given the PSF. In this work we address the second of these two steps for the case of saturated images, and assume that the PSF is known or has been estimated already. Unless otherwise stated, all the results in this work use the algorithm of [START_REF] Whyte | Deblurring shaken and partially saturated images[END_REF] to estimate a spatially-variant PSF. The algorithm is based on the method of [START_REF] Cho | Fast motion deblurring[END_REF], and estimates the PSF directly from the blurred image itself. Figure 1d shows the output of the proposed algorithm, which contains far fewer artefacts than the two existing algorithms shown for comparison.

Related Work

Saturation has not received wide attention in the literature, although several authors have cited it as the cause of artefacts in deblurred images [START_REF] Fergus | Removing camera shake from a single photograph[END_REF][START_REF] Cho | Fast motion deblurring[END_REF][START_REF] Tai | Richardson-Lucy deblurring for scenes under a projective motion path[END_REF]. Harmeling et al. (2010b) address the issue in the setting of multi-frame blind deblurring by thresholding the blurred image to detect saturated pixels, and ignoring these in the deblurring process. When multiple blurred images of the same scene are available, these pixels can be safely discarded, since there will generally remain unsaturated pixels covering the same area in other images.

Recently, [START_REF] Cho | Handling outliers in non-blind image deconvolution[END_REF] have also considered saturated pixels, in the more general context of non-blind deblurring with outliers, and propose an expectation-maximisation algorithm to iteratively identify and exclude outliers in the blurred image. Saturated pixels are detected by blurring the current estimate of the sharp image and finding places where the result exceeds the range of the camera. Those blurred pixels detected as saturated are ignored in the subsequent iterations of the deblurring algorithm. In Section 4 we discuss why simply ignoring saturated pixels is, in general, not sufficient to prevent artefacts from appearing in single-image deblurring.

In an alternative line of work, several authors have proposed algorithms for non-blind deblurring that are robust against various types of modeling errors, without directly addressing the sources of those errors. [START_REF] Yuan | Progressive inter-scale and intra-scale non-blind image deconvolution[END_REF] propose a non-blind deblurring algorithm capable of suppressing "ringing" artefacts during deblurring, using multi-scale regularisation. [START_REF] Yang | An efficient TVL1 algorithm for deblurring multichannel images corrupted by impulsive noise[END_REF] and [START_REF] Xu | Two-phase kernel estimation for robust motion deblurring[END_REF] also consider non-blind deblurring with robust data-fidelity terms, to handle non-Gaussian impulse noise, however their formulations do not handle arbitrarily large deviations from the linear model, such as can be caused by saturation.

Many algorithms exist for non-blind deblurring in the linear (non-saturated) case, perhaps most famously the Wiener filter [START_REF] Wiener | Extrapolation, Interpolation, and Smoothing of Stationary Time Series[END_REF]) and the Richardson-Lucy algorithm [START_REF] Richardson | Bayesian-based iterative method of image restoration[END_REF][START_REF] Lucy | An iterative technique for the rectification of observed distributions[END_REF]. Recently, many authors have focused on the use of regularisation, derived from natural image statistics, to suppress noise in the output while encouraging sharp edges to appear [START_REF] Krishnan | Fast image deconvolution using hyper-Laplacian priors[END_REF][START_REF] Levin | Image and depth from a conventional camera with a coded aperture[END_REF][START_REF] Joshi | Image deblurring using inertial measurement sensors[END_REF][START_REF] Afonso | Fast image recovery using variable splitting and constrained optimization[END_REF][START_REF] Tai | Richardson-Lucy deblurring for scenes under a projective motion path[END_REF]Zoran and Weiss 2011).

For the problem of "blind" deblurring, where the PSF is unknown, single-image blind PSF estimation for camera shake has been widely studied using variational and maximum a posteriori (MAP) algorithms [START_REF] Fergus | Removing camera shake from a single photograph[END_REF][START_REF] Shan | High-quality motion deblurring from a single image[END_REF][START_REF] Cho | Fast motion deblurring[END_REF][START_REF] Cai | Blind motion deblurring from a single image using sparse approximation[END_REF][START_REF] Xu | Two-phase kernel estimation for robust motion deblurring[END_REF][START_REF] Levin | Efficient marginal likelihood optimization in blind deconvolution[END_REF][START_REF] Krishnan | Blind deconvolution using a normalized sparsity measure[END_REF]. [START_REF] Levin | Understanding and evaluating blind deconvolution algorithms[END_REF] review several approaches and provide a ground-truth dataset for comparison on spatially-invariant blur. While most work has focused on spatially-invariant blur, several approaches have also been proposed for spatially-varying blur [START_REF] Whyte | Non-uniform deblurring for shaken images[END_REF][START_REF] Gupta | Single image deblurring using motion density functions[END_REF]Harmeling et al. 2010a;[START_REF] Joshi | Image deblurring using inertial measurement sensors[END_REF][START_REF] Tai | Richardson-Lucy deblurring for scenes under a projective motion path[END_REF].

The remainder of this paper proceeds as follows: We begin in Section 2 by providing some background on non-blind deblurring and saturation in cameras. In Section 3 we analyse some of the properties and causes of "ringing" artefacts (which are common when deblurring saturated images), and discuss the implications of this analysis in Section 4. Based on this discussion, in Section 5 we describe our proposed approach. We present deblurring results and comparison to related work in Section 6.

Background

In most existing work on deblurring, the observed image produced by a camera is modelled as a linear blur operation applied to a sharp image, followed by a random noise process. Under this model, an observed blurred image g (written as an N × 1 vector, where N is the number of pixels in the image) can be written in terms of a (latent) sharp image f (also an N × 1 vector) as

g * = Af
(1)

g = g * + ε, (2)
where A is an N × N matrix representing the discrete PSF, g * represents the "noiseless" blurred image, and ε is some random noise affecting the image. Typically, the noise ε is modelled as following either a Poisson or Gaussian distribution, independent at each pixel.

For many causes of blur, the matrix A can be parameterised by a small set of weights w, often referred to as a blur kernel, such that

A = k w k T k , (3)
where each N ×N matrix T k applies some geometric transformation to the sharp image f . Classically, T k have been chosen to model translations of the sharp image, allowing Equation (1) to be written as a 2D convolution of f with w.

For blur caused by camera shake (motion of the camera during exposure), recent work [START_REF] Gupta | Single image deblurring using motion density functions[END_REF][START_REF] Joshi | Image deblurring using inertial measurement sensors[END_REF][START_REF] Whyte | Non-uniform deblurring for shaken images[END_REF][START_REF] Tai | Richardson-Lucy deblurring for scenes under a projective motion path[END_REF] has shown that using projective transformations for T k is more appropriate, and leads to more accurate modeling of the spatially-variant blur caused by camera shake. The remainder of this work is agnostic to the form of A, and thus can be applied equally to spatially-variant and spatially-invariant blur.

Non-blind deblurring (where A is known) is generally performed by attempting to solve a minimisation problem of the form

min f L(g, Af ) + αφ(f ), (4)
where the data-fidelity term L penalises sharp images that do not closely fit the observed data (i.e. L is a measure of "distance" between g and Af ), and the regularisation term φ penalises sharp images that do not adhere to some prior model of sharp images. The scalar weight α balances the contributions of the two terms in the optimisation.

In a probabilistic setting, where the random noise ε is assumed to follow a known distribution, the data-fidelity term can be derived from the negative log-likelihood:

L(g, Af ) = -log p(g|Af ), (5)
where p(g|Af ) denotes the probability of observing the blurry image g, given a sharp image f and PSF A (often referred to as the likelihood). If the noise follows pixelindependent Gaussian distributions with uniform variance, the appropriate data-fidelity term is

L G (g, Af ) = i g i -(Af ) i 2 , (6)
where (Af ) i indicates the i th element of the vector Af . With Gaussian noise, Equation ( 4) is typically solved using standard linear least-squares algorithms, such as conjugate gradient descent [START_REF] Levin | Image and depth from a conventional camera with a coded aperture[END_REF]). For the special case of spatially-invariant blur, and provided that the regularisation term φ can also be written as a quadratic function of f , Equation (4) has a closed-form solution in the frequency domain, which can be computed efficiently using the fast Fourier transform [START_REF] Wiener | Extrapolation, Interpolation, and Smoothing of Stationary Time Series[END_REF][START_REF] Gamelin | Complex Analysis[END_REF]).

If the noise follows a Poisson distribution, the appropriate data-fidelity term is

L P (g, Af ) = - i g i log(Af ) i -(Af ) i . (7)
The classic algorithm for deblurring images with Poisson noise is the Richardson-Lucy algorithm [START_REF] Richardson | Bayesian-based iterative method of image restoration[END_REF][START_REF] Lucy | An iterative technique for the rectification of observed distributions[END_REF]), an iterative algorithm described by a simple multiplicative update equation. The incorporation of regularisation terms into this algorithm has been addressed by [START_REF] Tai | Richardson-Lucy deblurring for scenes under a projective motion path[END_REF] and [START_REF] Welk | Robust variational approaches to positivity-constrained image deconvolution[END_REF]. We discuss this algorithm further in Section 5.

A third data-fidelity term that is more robust to outliers than the two mentioned above, and which has been applied for image deblurring with impulse noise is the ℓ 1 norm [START_REF] Bar | Image deblurring in the presence of impulsive noise[END_REF][START_REF] Yang | An efficient TVL1 algorithm for deblurring multichannel images corrupted by impulsive noise[END_REF]) (corresponding to noise with a Laplacian distribution):

L L (g, Af ) = i g i -(Af ) i . (8)
Although this data-fidelity term is more robust against noise values ε i with large magnitudes, compared to the Gaussian or Poisson data-fidelity terms, it still produces artefacts in the presence of saturation [START_REF] Cho | Handling outliers in non-blind image deconvolution[END_REF]. For clarity, in the remainder of the paper we denote the data-fidelity term L(g, Af ) simply as L(f ), since we consider the blurred image g and the PSF matrix A to be fixed.

Sensor Saturation

Sensor saturation occurs when the radiance of the scene within a pixel exceeds the camera sensor's range, at which point the sensor ceases to integrate the incident light, and produces an output that is clamped to the largest output value. This introduces a non-linearity into the image formation process that is not modelled by Equations ( 1) and (2). To correctly describe this effect, our model must include a nonlinear function R, which reflects the sensor's non-linear response to incident light. This function is applied to each pixel of the image before it is output by the sensor, i.e.

g i = R g * i + ε i , (9)
where ε i represents the random noise on pixel i.

For the purpose of describing saturation, we model the non-linear response function R as a truncated linear function, i.e. R(x) = min(x, 1), for intensities scaled to lie in the range [0, 1]. This model is supported by the data in Figure 2, which shows the relationship between pixel intensities in three different exposures of a bright light source. The pixel values in the short exposure (with no saturation) and the longer exposures (with saturation) clearly exhibit this clipped linear relationship. As the length of the exposure increases, more pixels saturate.

Due to the non-linearity in the image formation process, simply applying existing methods for non-blind deblurring (which assume a linear image formation model) to images affected by saturation, produces deblurred images which exhibit severe artefacts in the form of "ringing" -medium frequency ripples that spread across the image, e.g. in Figures 1b and 1c. Clearly, the saturation must be taken into The plots show the likelihood p(g i |g * i ) of observed intensity g i given the noiseless value g * i , under a Poisson noise model for the case with and without sensor saturation. The top row shows the 2D likelihood as an intensity-map, where black is zero and brighter is larger. The second row shows vertical slices for several values of g * i . Without saturation, the likelihood remains uni-modal around g * i even for bright pixels. With saturation (far-right plot), bright pixels have multi-modal likelihoods (one mode at g * i , and a second at 255). This makes the inversion of such a forward model particularly difficult. account during non-blind deblurring in order to avoid such artefacts.

Given the non-linear forward model in Equation ( 9), it is tempting to modify the data-fidelity term L to take into account saturation, and thus prevent artefacts from appearing in the deblurred image. The model in Equation ( 9), however, is not tractable to invert. Since the noise term ǫ i lies inside R, the likelihood p(g i |g * i ) (from which L is derived) is distorted and, in general, can no longer be written in closed-form. Figure 3 shows the difference between the likelihoods with and without saturation for Poisson noise. In the saturated likelihood, pixels near the top of the camera's range have distributions that are no longer smooth and uni-modal, but instead have a second sharp peak at 1. Furthermore, for some pixels this second mode at 1 is the maximum of the likelihood, i.e.

P (g i = 1|g * i ) > P (g i = g * i |g * i )
, which clearly contradicts the normal behaviour of Poisson or Gaussian noise, where the likelihood is smooth and has a single mode at g i = g * i .

Given the difficulty of inverting the true non-linear forward model, alternative approaches to handling saturated pixels are needed, in order to prevent ringing artefacts from appearing. Before discussing our approach to this, we provide some analysis of ringing in general.

Ringing

Ringing is common in deblurred images, and has been attributed to various causes, including the Gibbs phenomenon [START_REF] Yuan | Image deblurring with blurred/noisy image pairs[END_REF]) (arising from the inability of a finite Fourier basis to reconstruct perfect step edges), incorrectlymodelled image noise (i.e. using the incorrect likelihood for the data-fidelity term L in Equation ( 4)), or errors in the estimated PSF [START_REF] Shan | High-quality motion deblurring from a single image[END_REF]). In the following we show that the root cause of ringing is the fact that, in general, blur annihilates certain spatial frequencies in the image. These spatial frequencies are very poorly constrained by the observed data (the blurred image) and can become artificially amplified in the deblurred image. Incorrectly-modelled noise and PSF errors are two causes of such amplification.

What does ringing look like?

For a discrete PSF given by a matrix A, there is typically some set of vectors that are almost completely suppressed by A. These are vectors that lie in, or close to, the nullspace of A, i.e. {x : Ax ≃ 0}. When we attempt to invert the blur process by estimating the sharp image f that minimises L( f ), these directions in the solution space are very poorly constrained by the data, since A( f + λx) ≃ A f , where λ is a scalar. As such, they can become artificially amplified in the deblurred image. These amplified components x, lying close to the PSF's nullspace, cause visible artefacts in the deblurred image.

The reason that poorly-constrained directions appear as ringing, and not as some other kind of visual artefact, is that the nullspace of a PSF tends to be spanned by vectors having small support in the frequency domain, but large spatial support. This can be seen from a spectral decomposition of the PSF. In the simplest case, where A corresponds to a convolution with cyclic boundary conditions, its eigendecomposition can be obtained by the discrete Fourier transform. Thus, its eigenvectors are sinusoids, and its nullspace is spanned by a set of these sinusoids with eigenvalues close to zero. When the cyclic boundary conditions are removed, or the PSF is spatially-variant, e.g. due to camera shake, the exact correspondence with the Fourier transform and sinusoids no longer holds, however the nullspace is still spanned by vectors resembling spatial frequencies, as shown in Figure 4. These are vectors with a small support in the frequency domain but a large spatial support.

A concrete example is shown in Figure 5, which shows a blurred image with outliers added, deblurred using the Richardson-Lucy algorithm. Looking at the ringing artefacts in both the spatial domain and the frequency domain, it is clear that ringing occurs at frequencies that are poorly constrained by the blur kernel (frequencies at which the magni- tude of the discrete Fourier transform (DFT) of the kernel is small).

How does ringing emerge?

Assuming that our estimate of the PSF contains no gross errors, the amplification of poorly-constrained directions in the solution may occur in two major ways. First, as mentioned above, there are spatial frequencies in the sharp image that are almost completely annihilated by A. During nonblind deblurring, any such spatial frequency x in the solution will be essentially unconstrained, since L(f +λx) ≃ L(f ). If anything occurs during deblurring to cause λ to be non-zero, there is no force in the optimisation to reduce it again, and the deblurred result will be the superposition of the sharp image with some combination of sinusoids x. Figures 6g and 6h demonstrate this with a synthetic example, where the deblurred result depends highly on the initialisation. With some initialisations, the deblurred image contains ringing, while with other initialisations, no ringing appears. This indicates that the ringing components of the solution are indeed poorly constrained, but are not actively amplified by the data-fidelity cost. In Section 4 we will discuss the cause of the amplification in saturated images such as these.

The second, and harder to tackle, way that ringing can arise is if there are some unmodelled / outlier pixels in the blurred image g, e.g. due to saturation or impulse noise. In this case, there may exist some x such that L(f * + λx) < L(f * ), where f * denotes the true sharp image. This is possible because, by their very nature, the outliers cannot be well-explained by the true sharp image f * . In order for this to happen, the addition of λx must (a) decrease the datafidelity cost L for the outlier pixels, while (b) not significantly increasing it for the remainder of the image. Clearly, to satisfy (a), x cannot lie exactly in the nullspace of A, otherwise Ax = 0, however to satisfy (b), it also cannot lie too far from the nullspace, otherwise the data-fidelity cost would grow quickly with λ. The result is that x lies close to, but not in, the nullspace of A, and the optimisation is moreor-less free to change λ in order to explain the outliers more accurately.

Figure 6f shows an example of this, where for all initialisations, some ringing appears, and the deblurred image f has a lower data-fidelity cost than the true sharp image f * . Note that the ringing in Figure 6f is visually similar to that in Figures 6g and6h, underlining the fact that in both cases, ringing appears in the poorly-constrained frequencies of the result.

One additional way that ringing may emerge in a deblurred image, which we do not address in this work, is if the PSF is incorrectly estimated. In this case, the PSF used to deblur is different from the one which caused the blur, and the deblurred image will be incorrect due to this discrepancy. We do not go into detail here, but for an estimated PSF  that is related to the true PSF A by  = A + ∆A, the deblurred image will contain ringing in the spatial frequencies where  has a small frequency response (i.e. close to the null-space of Â) [START_REF] Whyte | Removing Camera Shake Blur and Unwanted Occluders from Photographs[END_REF]. This echoes the conclusions of the previous paragraphs, with the difference that the ringing frequencies are determined by the incorrect PSF Â, instead of the true PSF A. The problem of deblurring with an erroneous PSF has also been addressed recently by [START_REF] Ji | Robust image deblurring with an inaccurate blur kernel[END_REF], who introduce additional latent variables in order to estimate f in a manner that is invariant to a certain class of PSF errors.

Why doesn't regularisation suppress ringing?

Often, ringing appears in deblurred images despite the inclusion of some regularisation term φ in the optimisation problem in Equation (4). This seems counterintuitive, since the purpose of the regularisation is to provide some additional constraints on the solution, particularly in directions which are poorly constrained by the data. However, most popular forms of regularisation used in non-blind deblurring penalise only the magnitudes of the first or second spatial derivatives of the deblurred image, e.g.

φ(f ) = i ρ |(d x * f ) i | + ρ |(d y * f ) i | , ( 10
)
where the sum is taken over all pixels i in the image, the filters d x and d y compute derivates, and ρ is a non-decreasing scalar function [START_REF] Blake | Visual Reconstruction[END_REF][START_REF] Bouman | A generalized Gaussian image model for edge-preserving MAP estimation[END_REF][START_REF] Schultz | A bayesian approach to image expansion for improved definition[END_REF][START_REF] Levin | Image and depth from a conventional camera with a coded aperture[END_REF]). These derivatives are computed from the differences between neighbouring pixels, and are essentially high-pass filters. Unsurprisingly, such regularisation works best at suppressing high spatial frequencies. At medium-to-low frequencies, the power of the regularisation decreases. Figure 7 demonstrates this, showing how the power of the regularisation is greatest at high frequencies, and falls to zero at low frequencies. The blurred data on the other hand constrains the lowest frequencies (i.e. at scales larger than the blur), but not the high frequencies. There is a region of medium frequencies that are poorly constrained by both. It is at these intermediate frequencies that some ringing x can easily appear, since L(f * + λx) ≃ L(f * ), and φ(f * + λf ) ≃ φ(f * ). Although the regularisation weight α in Equation ( 4) can be increased in order to reduce ringing, this will also begin to over-smooth the deblurred image. [START_REF] Yuan | Progressive inter-scale and intra-scale non-blind image deconvolution[END_REF] observe this fact and propose a multiscale non-blind deblurring algorithm capable of preventing ringing caused by noise and numerical inaccuracies. However, although their method can suppress ringing at a wide range of frequencies, it is still generally unable to handle ringing caused by saturation, as shown in Figure 14.

Preventing Ringing when Deblurring Saturated Images

When applying existing non-blind deblurring algorithms to saturated images, which assume the linear model g = Af + ε, using Gaussian or Poisson data-fidelity terms L, ringing is almost certain to appear. As discussed in the previous section, saturated pixels can cause the phenomenon where a deblurred image f containing ringing actually has a lower data-fidelity cost than the true sharp image f * , due to the fact that the assumption of linearity is violated.

Although the saturation can be modelled with a nonlinear response R (as discussed in Section 2.1) in order to prevent the case where L(f * + λx) < L(f * ), the model is not tractable to invert. For this reason, some recent authors have instead chosen to treat saturated pixels simply as outliers, and model the rest of the image as linear (Harmeling et al. 2010b;[START_REF] Cho | Handling outliers in non-blind image deconvolution[END_REF]. The saturated pixels are discarded and treated as missing data. Not only is this much more tractable to optimise, it arguably does not sacrifice much -saturated pixels are largely uninformative, due to being clamped at the maximum output value.

In order to perform non-blind deblurring with missing data, we can define a mask m of binary weights, where m i = 0 if pixel i is missing, and m i = 1 otherwise. We then construct a weighted form of the data-fidelity term, using m as the weights, and denote this L m . For example, for Poisson noise, Equation (7) becomes:

L P m (f ) = - i m i g i log(Af ) i -(Af ) i . ( 11
)
By removing all the data that does not follow the linear model, it should no longer be possible for L m (f * + λx) < L m (f * ).

Harmeling et al. (2010b

) estimate the mask m directly from the blurred image, by finding all the pixels in g that are above some threshold ϕ close to 1, i.e.

m i = 1 if g i < ϕ 0 otherwise. ( 12
)
On the other hand, [START_REF] Cho | Handling outliers in non-blind image deconvolution[END_REF] estimate m by blurring the current estimate of the sharp image f :

m i = 1 if (A f ) i < 1 0 otherwise. ( 13
)
Note that although some of the blurred pixels are discarded, the entire latent image f is still estimated.

Although detecting and ignoring the (outlier) saturated pixels can reduce ringing in the deblurred images, it does not necessarily remove it entirely. By discarding saturated data, the data-fidelity term L m no longer actively encourages ringing. However, it may still be possible that L m (f * + λx) ≃ L m (f * ) for some ringing x. As discussed in the previous section, if anything in the data, the initialisation, or in the deblurring algorithm serves to increase λ, there is no force to reduce it again, and the ringing will appear in the final result. Figure 6g demonstrates this with a synthetic example. Even when the outliers are known exactly, and are completely removed from the deblurring process, ringing can still appear. On the other hand, when we initialise the deblurring algorithm with the true sharp image, no ringing appears, indicating that there is nothing about L m that causes ringing, there is simply nothing to prevent or suppress it.

From this discussion, we conclude that in order to deblur saturated images without introducing ringing, we must do more than simply remove the saturated pixels from the observed image g. We must also avoid amplifying the poorlyconstrained spatial frequencies in the latent image f . To achieve this, we could either introduce some form of regularisation to place additional constraints on f , such as in [START_REF] Yuan | Progressive inter-scale and intra-scale non-blind image deconvolution[END_REF], or we could directly modify the data-fidelity term. In this work we choose the latter approach, eschewing the use of regularisation to prevent ringing (although we revisit regularisation later in order to reduce other sources of noise in the deblurred results).

In this work, we posit that for the case of saturation, the main factor causing poorly-constrained spatial frequencies to become amplified is that there exist pixels in the sharp image that exceed the image's range, i.e. ∃ j : f j > 1. During non-blind deblurring, we must estimate these "bright" pixels' values, and it is the act of estimating these values that destabilises our estimates of other pixels. Note that removing the saturated pixels from g does not, in general, remove the influence of these "bright" latent pixels from the blurred image; a latent pixel with intensity greater than 1 can contribute to a blurred pixel whose intensity is less than 1.

Indeed, a blurred image with no saturation at all may still correspond to a latent image containing pixels brighter than 1.

When we attempt to estimate a "bright" pixel's value f j , we are likely to make a relatively inaccurate estimate. The main reason for this is that the set of observations (blurred pixels) concerning f j is likely to be incomplete. In the blurring process, each latent pixel in the sharp image is spread across multiple pixels in the blurred image g. Given that f j > 1, there is a good chance that some of these blurred pixels will be saturated, and thus uninformative as to the value of f j . With fewer observations from which to estimate f j , its value will be more susceptible to noise, and contain larger error than if we had a full set of unsaturated observations available. This is supported by Figure 6g, where the sharp image contains pixels with intensity greater than 1, and the appearance of ringing depends on the initialisation. When we initialise with the true sharp image, which already contains the bright (> 1) pixels, no ringing is caused. When we initialise with the blurred image, or random values in [0, 1], the algorithm is forced to attempt to estimate the bright values, and in doing so causes ringing.

Since this amplification is caused when we attempt to estimate the "bright" pixels in f , our approach to preventing ringing is based on the idea that we should avoid (or rather, delay) estimating such bright pixels. If we only estimate pixels that can be accurately estimated, we will not make gross errors, which in turn will avoid introducing ringing. For the synthetic example, Figure 6h shows the result of this approach. By refusing to estimate the bright (> 1) pixels in f , the ringing is almost removed for all initialisations. In the next section we discuss how we do this in practice.

Note that this idea is similar to the notion of "Highest Confidence First" [START_REF] Chou | The theory and practice of Bayesian image labeling[END_REF]. Chou and Brown point out that due to the coupling of pixels in a Markov Random Field image model, the estimate for a pixel with strong observations may be negatively impacted by a poor decision at a neighbouring pixel with weak observations. In such cases, the pixel with strong observations "can do better without the incorrect information of a neighbor."

Proposed Method

In this section we describe our proposed algorithm for deblurring images containing saturated pixels. We begin by describing our approach to estimating the latent image without introducing ringing, motivated by the discussion in the previous sections. In addition, we propose a method for estimating the bright pixels separately, without introducing artefacts, and finally describe how these elements are combined into the complete proposed method. As can be seen, in (f), ringing appears regardless of the initialisation, indicating that L is encouraging this to occur. In (g), ringing appears with some initialisations but not others, indicating that L does not encourage it, but does not suppress it either. In (h), ringing is almost gone, since we have removed the most destabilising data. Counterintuitively, although we have discarded more data, we end up with less ringing. Note that by removing all the blurred pixels shown in (d), we have no information about the bright pixels in f , and they simply retain their initial values. Note also that the deblurred images f may contain pixel values greater than 1, however they are clamped to 1 before writing the image to file

Preventing ringing by refusing to make bad estimates

As we have seen, it is clear that when we attempt to perform non-blind deblurring on an image containing saturation, we will incorrectly estimate the "bright" latent pixels (with values close to or greater than 1) that caused the saturation. As discussed in the previous section, the errors we make at these bright pixels will cause ringing. In this section we propose a method that iteratively attempts to classify these latent pixels, and remove them from the estimation process, thereby mitigating their effect on the rest of the image.

Assume for now that we already know which pixels in the latent image are bright, and will thus be poorly estimated. We denote this set of pixels S, and denote its set complement (containing the rest of the image, which can be estimated accurately) by U . We can write the latent image in terms of these two disjoint sets:

f = f U + f S .
Given that we are unable to estimate f S accurately, our aim is to prevent the propagation of errors from f S to f U . To achieve this, we decouple our estimate of f U from our estimate of f S . First, note that we can decompose the noiseless blurred image as:

g * = Af U + Af S . ( 14
)
We denote by V the set of blurred pixels that are not affected by f S , i.e. V = i (Af S ) i = 0 , and construct the corresponding binary mask v (where v i = 1 if i ∈ V, and v i = 0 otherwise). Then, given that v • Af S = 0, we can The blur kernel has minima in its power spectrum, which correspond to poorly-constrained frequencies in the deblurred solution. At the low-frequency minima (e.g. near 5 cycles/image), the gradient filter has also lost most of its power, and so ringing at these frequencies is unlikely to be suppressed by gradientbased regularisation write

v • g * = v • Af U (15)
where • • • represents the element-wise product between two vectors. From this, we can estimate f U independently of f S . Note that we can obtain the mask v simply by constructing the binary mask u that corresponds to the set U , and performing a binary erosion of u with the non-zeros of the PSF.

Furthermore, given that the set U does not contain any bright latent pixels, Af U should not cause any saturation. Thus, we can estimate f U without modelling the non-linearity caused by saturation, and write the observed blurred image, including noise, as

v • g = v • Af U + v • ε. ( 16
)
To estimate f U from v • g, we can adapt existing nonblind deblurring algorithms to handle missing data, as seen in Section 4. For Poisson noise, the data-fidelity term is simply weighted, as in Equation ( 11), using v as the weights to form L P v . The classical algorithm for non-blind deblurring with Poisson noise is the Richardson-Lucy algorithm, which is described by the multiplicative update equation:

f t+1 = f t • A ⊤ g A f t , (17)
where the division is performed element-wise.

We can adapt the Richardson-Lucy update equation from Equation ( 17) to minimise the binary-weighted data-fidelity term L P v , instead of the homogeneous L P , leading to the fol-lowing update equation for f U :

f t+1 U = f t U • A ⊤ g • v A f t U + 1 -v . (18)
Note that in order to avoid division by zero, we add a small positive constant to the denominator of the fraction. This update equation is applicable to any kind of missing data in g, the important distinction is in how we determine v, the mask of missing pixels. In Figure 6h, we show the result of deblurring with v determined as described above, by eroding u. For comparison, Figure 6g shows the result when v is determined simply by finding saturated pixels in g, i.e. v i = 1 if g i < 1. Despite removing all the saturated pixels, the results in Figure 6g still contain ringing, while significantly less ringing appears in Figure 6h using the proposed approach.

Preventing artefacts inside bright regions

While the update method in the previous section is effective at preventing ringing from propagating outwards from bright regions into other parts of the image, we still wish to estimate values for the pixels in S. Since these pixels are bright, and unlikely to be estimated accurately, we are not concerned with preventing propagation of information from U to S, and choose to update them using all the available data.

Using the standard Richardson-Lucy update for these bright regions however, can cause dark artefacts to appear, due to the fact that the linear model cannot explain the saturated data. Such artefacts are visible in Figure 10d. To prevent these artefacts, we propose a second modification of the Richardson-Lucy algorithm, to include a non-linear response function. Since the true forward model, discussed in Section 2.1 is not tractable to invert, we propose the following tractable but approximate alternative, which is sufficient for preventing artefacts in bright regions:

g i = R(g * i ) + ε i , (19)
where now the response R function comes before the noise.

Re-deriving the Richardson-Lucy algorithm with this model gives the following update equation for the pixels in S:

f t+1 S = f t S • A ⊤ g • R ′ (A f t ) R(A f t ) + 1 -R ′ (A f t ) , (20)
where R ′ (•) indicates the derivative of the response function R.

Since the ideal response function R(x) = min(x, 1) is not differentiable, we use a smooth, continuously differentiable approximation to this function (Chen and Mangasar-

0 0.5 1 1.5 2 0 0.5 1 x R(x) R ′ (x)
Fig. 8. Modelling the saturated sensor response. Smooth and differentiable approximation to the non-differentiable function min(x, 1) used to model sensor saturation, defined in Equation ( 21). The derivative is also smooth and defined everywhere.

ian 1996), where

R(x) = x - 1 a log 1 + exp a(x -1) (21) R ′ (x) = 1 1 + exp a(x -1) . ( 22
)
The parameter a controls the smoothness of the approximation, and we have found a = 50 to be a suitable compromise between smoothness and accuracy (we use this value in all results presented in this paper). Figure 8 shows the shape of these smooth versions of R and R ′ . Equation ( 20) can be roughly interpreted as weighting the blurry pixels according to the value of R ′ : in the linear (unsaturated) portion where x < 1, R(x) ≃ x and R ′ (x) ≃ 1, so that the term in parentheses is the same as for the standard RL algorithm. In the saturated portion where x > 1, R(x) ≃ 1 and R ′ (x) ≃ 0, so that the term in parentheses is equal to unity and has no influence on f . We note that this behaviour is very similar to the method used by [START_REF] Cho | Handling outliers in non-blind image deconvolution[END_REF] to handle saturation -given the current estimate f of the sharp image, they compute the blurred estimate A f , and ignore any blurry pixel i for which (A f ) i > 1.

Segmenting the latent image

So far in this section, we have assumed that S and U are known. Given a real blurred image however, we do not know a priori which parts of f belong in U and which in S. We thus treat U as another latent variable, in addition to the latent image f , and perform the segmentation at each iteration t. Given the discussion in previous sections, we adopt a simple segmentation process, which consists of thresholding the latent image at some level ϕ close to 1:

U = j f t j ≤ ϕ . (23)
We decompose f t according to

f t U = u • f t (24) f t S = f t -f t U , (25)
where u is the binary mask corresponding to U . Since our aim is to ensure that no large errors are introduced in f U , we set the threshold low enough that most potentially-bright pixels are assigned to S. Empirically, we choose ϕ = 0.9 for the results in this paper, although we have not found the results to be particularly sensitive to the exact value of ϕ.

Adding regularisation

Although the focus of this work is mitigating ringing by modifying the data-fidelity term, rather than by designing new forms of regularisation, it may nonetheless still be useful to apply some form of regularisation to reduce other noise throughout the deblurred image. As discussed by [START_REF] Tai | Richardson-Lucy deblurring for scenes under a projective motion path[END_REF] and [START_REF] Welk | Robust variational approaches to positivity-constrained image deconvolution[END_REF] it is possible to include a regularisation term in the Richardson-Lucy algorithm, and this remains true for our algorithm. In this work we include the ℓ 1 -norm of the gradients of the deblurred image as regularisation. Using the form of Equation ( 10), this is written

φ(f ) = i (d x * f ) i + (d y * f ) i , (26)
where the filters d x and d y compute derivates. We incorporate this in the manner described by [START_REF] Tai | Richardson-Lucy deblurring for scenes under a projective motion path[END_REF]. Denoting the unregularised update of f as f t+1 unreg (computed as in Equations ( 18) and (20) of the previous sections), we compute the update for the regularised problem as

f t+1 = f t+1 unreg 1 + α∇φ( f t ) , (27)
where ∇φ( f t ) is the vector of partial derivatives of the regularisation1 φ with respect to each pixel of the latent image f , evaluated at f t , i.e.

∇φ( f t ) j = ∂φ(f ) ∂f j f t . (28)
In Figure 9 we compare results obtained using our method with and without this regularisation. As can be seen, the ringing is largely suppressed by the steps described in the preceding sections. The addition of the regularisation term further improves the results by reducing noise elsewhere in the deblurred result. We summarise our complete proposed algorithm in Algorithm 1. Figure 10 shows the contributions of the two pro- posed modifications for a synthetic 1D example. As is visible, the use of the decoupled update prevents ringing from spreading across the deblurred result, while the use of the non-linear forward model prevents dark artefacts from appearing inside the bright regions. By combining these two methods, we produce our best estimate of the pixels in S and U , while preventing the errors we make in S from affecting the rest of the image.

Implementation

In this section we describe some of the implementation details of the proposed algorithm. When segmenting the current estimate of the latent image, we take additional steps to ensure that we make a conservative estimate of which pixels Compute regularised update f t+1 using ( 27)

9 end
can be estimated accurately. First, after thresholding the latent image in Equation ( 23), we perform a binary erosion on U , such that

U = j f j ≤ ϕ ⊖ M, (29)
where ⊖ denotes binary erosion, and the structuring element M used for erosion is a disk of radius 3 pixels. This ensures that all poorly-estimated pixels are correctly assigned to S (perhaps at the expense of wrongly including some wellestimated pixels too). Performing this step improves the deblurred results, since it is not only the bright pixels whose value is likely to be inaccurate due to saturation, but their neighbours too. Fewer artefacts arise from wrongly assigning a well-estimated pixel into S than the other way around. Second, in order to avoid introducing visible boundaries between the two regions, we blur the mask u slightly using a Gaussian filter with standard deviation 3 pixels when decomposing the current latent image f t into f t U and f t S in Equations ( 24) and (25).

Results

Figures 1 and 12 show results of non-blind deblurring using the proposed algorithm described in Section 5 on real handheld photographs. For comparison, in Figure 1 we show results produced with the standard Richardson-Lucy algorithm, and the algorithm of [START_REF] Krishnan | Fast image deconvolution using hyper-Laplacian priors[END_REF], neither of which is designed to take account of saturation. As such, both algorithms produce large amounts of ringing in the deblurred images. In Figure 12 we compare against a baseline approach, where saturated regions are simply masked out by detecting pixels in the blurred image that exceed a threshold value of ϕ = 0.9, as in Equation ( 12), and discarding those pixels. We found that dilating the masked regions using a In both Figures 1 and12, the (spatially-variant) PSFs for these images were estimated from the blurred images themselves using our MAP-type blind deblurring algorithm [START_REF] Whyte | Non-uniform deblurring for shaken images[END_REF], which is based on the algorithm of [START_REF] Cho | Fast motion deblurring[END_REF]. The only modification required to estimate PSFs for saturated images using this blind algorithm is to discard potentially-saturated regions of the blurred image using a threshold. Since, in this case, the aim is only to estimate the PSF (and not a complete deblurred image), we can safely discard all of these pixels, since the number of saturated pixels in an image is typically small compared to the total number of pixels. There will typically remain sufficient unsaturated pixels from which to estimate the PSF.

Note in Figure 1 that the standard Richardson-Lucy algorithm and the algorithm of [START_REF] Levin | Image and depth from a conventional camera with a coded aperture[END_REF] produce ringing around the saturated regions, while the proposed algorithm reduces this without sacrificing deblurring quality elsewhere. In all results in this paper we performed 50 iterations of all Richardson-Lucy variants.

Figures 13 and14 show the results of non-blind deblurring using our algorithm, alongside those of the recentlyproposed algorithm of [START_REF] Cho | Handling outliers in non-blind image deconvolution[END_REF]. The blurred images and their spatially-invariant PSFs are provided by Cho et al., along with their deblurred results2 . In most cases our results exhibit less ringing than those of [START_REF] Cho | Handling outliers in non-blind image deconvolution[END_REF], which is due to the fact that we explicitly decouple the estimates of bright pixels from other pixels, in addition to removing saturated blurred pixels.

In order to gauge the quantitative difference between deblurred images containing ringing and the results of our algorithm, we generated a set of synthetically blurred and saturated images, with varying degrees of saturation and blur. We then deblurred using the original Richardson-Lucy algorithm, which produces ringing, and our proposed algorithm. Figure 11 shows the results, indicating that our method pro-duces a measurable improvement in deblurred image quality, which increases as the amount of saturation increases.

Conclusion

In this paper we have developed an approach to the problem of non-blind deblurring of images blurred by camera shake and suffering from saturation. We have provided an analysis of the properties and causes of "ringing" artefacts in deblurred images, as they apply to saturated images. Based on this analysis, we have proposed a non-blind deblurring algorithm, derived from the Richardson-Lucy algorithm, which is able to deblur images containing saturated regions without introducing ringing, and without sacrificing detail in the result. The algorithm is underpinned by the principle that we should prevent gross errors in bright regions from propagating to other regions. We provide an implementation of our algorithm online at http://www.di.ens.fr/willow/ research/saturation/.

As future work, our algorithm could potentially be extended to handle other sources of ringing, such as moving objects, impulse noise, or post-capture non-linearities (such as JPEG compression). Whenever it is possible to identify poorly-estimated latent pixels, our approach has the potential to reduce artefacts by decoupling these pixels from the rest of the image. In addition, this underlying principle of decoupling poorly-estimated latent pixels could also be applied within other non-blind deblurring algorithms, that are faster or more suitable for different noise models than the Richardson-Lucy algorithm, e.g. [START_REF] Levin | Image and depth from a conventional camera with a coded aperture[END_REF][START_REF] Wang | A new alternating minimization algorithm for total variation image reconstruction[END_REF][START_REF] Krishnan | Fast image deconvolution using hyper-Laplacian priors[END_REF][START_REF] Afonso | Fast image recovery using variable splitting and constrained optimization[END_REF][START_REF] Almeida | Deconvolving images with unknown boundaries using the alternating direction method of multipliers[END_REF].

One alternative direction for future work is the investigation of new regularisers that are targeted at suppressing ringing. As discussed in Section 3.3, gradient-based regularisation is generally insufficient to suppress ringing, and regularisers with a larger bandwidth in the frequency domain are needed. The multi-scale method of [START_REF] Yuan | Progressive inter-scale and intra-scale non-blind image deconvolution[END_REF] has this property, and patch-based methods such as that of Zoran and Weiss (2011) might also prove effective. As the scale factor increases, more pixels in the synthetic blurred image become saturated. In the table, we report the structural similarity index (SSIM) [START_REF] Wang | Image quality assessment: From error visibility to structural similarity[END_REF]) between the true sharp image and the deblurred image, using the original Richardson-Lucy (RL) algorithm, which causes ringing, and our proposed algorithm. The SSIM of the deblurred images decreases both when the amount of blur increases, and when the amount of saturation increases. However, in almost all cases, our algorithm achieves the same or better results, and the improvement over the original Richardson-Lucy algorithm increases with the amount of saturation.

Fig. 1 .

 1 Fig.1. Deblurring in the presence of saturation. Existing deblurring methods, such as those in (b) and (c), do not take account of saturated pixels. This leads to large and unsightly artefacts in the results, such as the "ringing" around the bright lights in the zoomed section. Using the proposed method (d), the ringing is greatly reduced and the quality of the deblurred image is improved. The PSF for this 1024 × 768 pixel image causes a blur of about 35 pixels in width, and was estimated directly from the blurred image using the algorithm described by[START_REF] Whyte | Non-uniform deblurring for shaken images[END_REF].

Fig. 2 .

 2 Fig. 2. Saturated and unsaturated photos of the same scene. (a)-(c) 3 different exposure times for the same scene, with bright regions that saturate in the longer exposures. A small window has been extracted, which is unsaturated at the shortest exposure, and increasingly saturated in the longer two. (d) Scatter plot of the intensities in the small window in (b) against those in the window in (a), normalised by exposure time. (e) Scatter plot of the intensities in the window in (c) against the window in (a), normalised by exposure time. The scatter plots in (d) and (e) clearly show the clipped linear relationship expected.

Fig. 3 .

 3 Fig.3. Noise distributions with and without sensor saturation. The plots show the likelihood p(g i |g * i ) of observed intensity g i given the noiseless value g * i , under a Poisson noise model for the case with and without sensor saturation. The top row shows the 2D likelihood as an intensity-map, where black is zero and brighter is larger. The second row shows vertical slices for several values of g * i . Without saturation, the likelihood remains uni-modal around g * i even for bright pixels. With saturation (far-right plot), bright pixels have multi-modal likelihoods (one mode at g * i , and a second at 255). This makes the inversion of such a forward model particularly difficult.

Fig. 4 .

 4 Fig. 4. Poorly-constrained frequencies of a 1D blur kernel. (a) A blur kernel. (b) We construct the 50×50 matrix A that convolves a 50sample signal with (a) (with zero-padding boundary conditions), and plot the 10 singular vectors with smallest singular values. The singular vectors have large spatial support, and correspond closely to spatial frequencies and their harmonics.

  Fig.5. An example of ringing in the frequency domain. A sharp image f * containing several pixels with values greater than 1 is convolved with a kernel w (a), and clipped at 1 to produce a blurred, saturated image g (b). When we deconvolve g (c), ringing is produced (d). We can take the discrete Fourier transform (DFT) of the kernel F (w) (e), and the ringing artefacts F ( ff * ) (f), and compare their magnitudes. In these plots, large values are red, while small values are blue. It is visible from these plots that the ringing frequencies with largest magnitude are those frequencies where the kernel has the smallest magnitude (i.e. the red regions in (f) correspond to dark blue regions in (e)). By isolating the largest peaks in the ringing spectrum |F ( ff * )| (g), and visualizing the spatial frequencies that they represent (h)-(j), we can see that they do indeed correspond to the ringing that is visible in the deblurred image.

Fig. 6 .

 6 Fig. 6. Synthetic example of deblurring with saturated pixels. (a) A sharp image f * containing some bright pixels with intensity twice the maximum image value. (b) The sharp image f is convolved with w and clipped at 1 to simulate saturation. No noise is added in this example. (c) The mask of pixels in g that are not saturated, i.e. m i = 1 (white) if (f * * w) i < 1, and m i = 0 (black) otherwise. (d) The mask of pixels in g that are not influenced by the bright pixels in f * , i.e. m i = 1 if (f * ≥ 1) * w i = 0. See Section 5.1 for further explanation. The following rows show the results of deblurring g with 1000 iterations of the Richardson-Lucy algorithm, starting from three different initialisations. (1 st row) initialised with blurred image. (2 nd row) initialised with random values in [0, 1].(3 rd row) initialised with true sharp image. As can be seen, in (f), ringing appears regardless of the initialisation, indicating that L is encouraging this to occur. In (g), ringing appears with some initialisations but not others, indicating that L does not encourage it, but does not suppress it either. In (h), ringing is almost gone, since we have removed the most destabilising data. Counterintuitively, although we have discarded more data, we end up with less ringing. Note that by removing all the blurred pixels shown in (d), we have no information about the bright pixels in f , and they simply retain their initial values. Note also that the deblurred images f may contain pixel values greater than 1, however they are clamped to 1 before writing the image to file

Fig. 7 .

 7 Fig. 7. Power spectra of blur and gradient filters. The power spectrum of the 1D blur kernel in Figure4(solid blue line), and the power spectrum of the derivative filter [1, -1] (dashed red line), often used for regularisation in non-blind deblurring. The blur kernel has minima in its power spectrum, which correspond to poorly-constrained frequencies in the deblurred solution. At the low-frequency minima (e.g. near 5 cycles/image), the gradient filter has also lost most of its power, and so ringing at these frequencies is unlikely to be suppressed by gradientbased regularisation

  Fig. 9. Effect of regularisation. This figure shows a real shaken image, and the effect of our algorithm with and without regularisation. Note that the ringing around saturated regions, visible in (b) is removed by our method (c), even without regularisation. By adding regularisation (d), the remaining noise (visible in textureless regions) is ameliorated.

Algorithm 1 :

 1 Proposed algorithm Input: Blurred image g, blur descriptor w Output: Deblurred image f 1 f 0 ← g 2 for t = 0 to num_iterations do 3 Decompose f t into f t U and f t S using (23) to (25) 4 Compute set V in the blurred image by eroding U with

  Fig. 10. Synthetic 1D example of blur and saturation. Each row shows a sharp "top-hat" signal, blurred using the filter shown at the top. Gaussian noise is added and the blurred signal is clipped. The last four columns show the deblurred outputs for the Richardson-Lucy algorithm, and for our two proposed modifications (the split update step in Section 5.1 and the non-linear forward model in Equation (19) in Section 5.2), separately and together. (Top row) With no saturation, all algorithms produce similar results. (Middle and bottom rows) When some of the blurred signal is saturated (region B), the Richardson-Lucy algorithm (c) produces an output with large ringing artefacts. Although region A is not itself saturated, the ringing propagates outwards from B & C across the whole signal. (d) Performing the split update, decoupling poorly-estimated bright pixels from the rest of the image, reduces ringing, but dark artefacts remain inside the bright region. (e) Deblurring the whole image with the non-linear model prevents artefacts inside the bright regions, but ringing is caused. (f) By combining the split update step with the non-linear forward model, the deblurred signal contains the least ringing, and has no dark artefacts inside the bright region

  Fig. 11. Quantitative effect of our algorithm on reducing artefacts. Starting from the sharp image shown on the left, we synthesize a set of blurry images with varying degrees of blur and saturation by scaling the intensities by increasing amounts, and blurring with several sizes of horizontal linear blur. As the scale factor increases, more pixels in the synthetic blurred image become saturated. In the table, we report the structural similarity index (SSIM)[START_REF] Wang | Image quality assessment: From error visibility to structural similarity[END_REF]) between the true sharp image and the deblurred image, using the original Richardson-Lucy (RL) algorithm, which causes ringing, and our proposed algorithm. The SSIM of the deblurred images decreases both when the amount of blur increases, and when the amount of saturation increases. However, in almost all cases, our algorithm achieves the same or better results, and the improvement over the original Richardson-Lucy algorithm increases with the amount of saturation.

  Fig. 12. Deblurring saturated images. Note that the ringing around saturated regions, visible in (b) and (c) is removed by our method (d), without causing any loss in visual quality elsewhere.

  Fig.13. Comparison to the method of[START_REF] Cho | Handling outliers in non-blind image deconvolution[END_REF]. The blurred images, the spatially-invariant PSFs and the results of their method were provided by the authors. Note in the close-ups and elsewhere in the images, our results generally contain less ringing than those of[START_REF] Cho | Handling outliers in non-blind image deconvolution[END_REF], without sacrificing detail. For example, above the car in the left column, and along the top edge of the bikes in the middle column.

Fig. 14 .

 14 Fig. 14. Comparison to Cho et al. (2011) and Yuan et al. (2008). This figure compares non-blind deblurring results, on saturated images, for (c) the algorithm of Cho et al. (2011), (d) the Richardson-Lucy-based algorithm of Yuan et al. (2008), and (e) our proposed algorithm. While both our algorithm and the algorithm of Yuan et al. are based on the Richardson-Lucy algorithm, our results contain much less ringing, due to handling the saturated pixels explicitly. Compared to the method of Cho et al., our results contain similar or less ringing. Results in (c) and (d) provided by Cho et al..

We note that the regulariser φ in Equation (26) is not continuously differentiable, however this can be avoided by using the standard approximation |x| ≃ √ ǫ +

x 2 , for some small number ǫ.

http://cg.postech.ac.kr/research/deconv_ outliers/ (accessed November 12, 2011).

Acknowledgements This work was partly supported by the MSR-INRIA laboratory, the EIT ICT labs and ERC grant VisRec no. 228180.


A Cohen

Nathalie Fraysse

Jean Rajchenbach

Médéric Argentina

Yann Bouret

Christophe Raufaste

Capillary-hydraulic jump in liquid foams

come

Capillary-hydraulic jump in liquid foams

Liquid foams are dispersion of gas bubbles in a continuous liquid phase. The organization of the liquid-gas interfaces is of paramount importance to understand the mechanical stability of such media. While soap films result from the contact between two bubbles, the contacts between three bubbles form a network of micrometer-thick liquid channels called "Plateau borders". Due to their distinctive geometry, these unbounded channels (i.e. without rigid boundary conditions) are not prone to the Rayleigh-Plateau instability; this ensures the stability of the bubbles assembly.

The dynamics and stability of Plateau borders regarding a local perturbation are not well understood. Only few experiments were devoted to the Plateau border scale [START_REF] Panhuis | New variations on the soap film experiments of Plateau, I. experiments under forced drainage[END_REF][START_REF] Koehler | Foam drainage on the microscale, II. imaging flow through single Plateau border[END_REF][START_REF] Pitois | Hydrodynamic resistance of a single foam channel[END_REF][START_REF] Elias | Elasticity of a soap film junction[END_REF] and none of them dealt with transient regimes. We study a single, centimeter-long, horizontal, Plateau border (PB). We perturb it locally by dropping a droplet of liquid, which coalesces with the PB, and we follow the dynamics of the subsequent liquid imbibition. Because the PB is horizontal, gravity effects are negligible and the liquid flow is driven by capillarity only.

The transient flow is analyzed for several Plateau border radii and several bulk viscosities of the surfactant solution (in the high surface mobility regime). We find that an inertial regime precedes a viscous one, on a critical length that decreases as the bulk viscosity increases. The critical length can be as large as several centimeters for very thin Plateau borders and low bulk viscosities. This inertial regime is reported for the first time; it is characterized by the formation of a hydraulic jump, and by a fast (0.1-1 m/s) and constant imbibition velocity. The resulting Reynolds number is thus larger than one, even for micrometer-thick channels. We show that the liquid redistribution is triggered by an effective negative surface tension [START_REF] Géminard | Observation of negative line tensions from Plateau border regions in dry foam films[END_REF] and we propose a model based on the hydraulic jump geometry [START_REF] Cohen | Capillary hydraulic jump in a Plateau border and its applications for liquid foam stability[END_REF].

This study should have applications for liquid foam drainage [START_REF] Saint-Jalmes | Physical chemistry in foam drainage and coarsing[END_REF], in particular to account for the exact shape of the front zone in forced drainage experiments [START_REF] Grassia | Quasione-dimensional foam drainage[END_REF]. The inertial regime might also be relevant to understand the liquid reorganization following bubbles rearrangements [START_REF] Biance | How Topological Rearrangements and Liquid Fraction Control Liquid Foam Stability[END_REF].

Figures:

  Figures: As sketched (top left), the experiment consists in dropping a droplet of the same surfactant solution on a single horizontal PB. The set of snapshots on the right shows the relaxation of the local perturbation after coalescence of the droplet with the PB, in the capillary-inertial regime. An outline of the capillary hydraulic jump geometry is given on the bottom left.

Alexandre Cohen 1 , Nathalie Fraysse 1 , Jean Rajchenbach 1 , Médéric Argentina 2 , Yann Bouret 1 and Christophe Raufaste 1

	(1) Université Nice Sophia Antipolis, CNRS, LPMC, UMR7336, 06108 Nice cedex2, France,
	(2) Université Nice Sophia Antipolis, CNRS, INLN, UMR7335, 06560 Valbonne, France
	alexandre.cohen@unice.fr


Alexandre Cohen

Nathalie Fraysse

Jean Rajchenbach

Médéric Argentina

Yann Bouret

Christophe Raufaste

CAPILLARY HYDRAULIC JUMP IN A PLATEAU BORDER

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

CAPILLARY HYDRAULIC JUMP IN A PLATEAU BORDER

Applications to foam drainage

Alexandre Cohen (1), Nathalie Fraysse (1), Jean Rajchenbach (1), Médéric Argentina (2), Yann Bouret (1), and Christophe Raufaste (1) (1) Université Nice Sophia Antipolis, CNRS, LPMC, UMR7336, 06108 Nice cedex2, France,

(2) Université Nice Sophia Antipolis, CNRS, INLN, UMR7335, 06560 Valbonne, France alexandre.cohen@unice.fr

By means of experimental measurements, we study transient flows of liquid along a single Plateau border (PB). In opposition with a common belief, we show that a fast, inertial regime of flow, very little dependent on the viscosity, can be encountered in the high surface mobility limit.

This newly observed regime evolves into the regular viscous regime after a critical distance that can be as large as several centimeters for very dry foams. We show that the liquid redistribution is triggered by an effective negative surface tension [1] that results from the distinctive geometry of the PBs. The dynamics is characterized by the formation of a hydraulic jump, which travels at a constant and substantial velocity (0.1 to 1 m.s -1 ). The relevant Reynolds number is much larger than one, even for micrometer-thick PBs. Therefore, the now widely acknowledged "drainage equation" [2], which results from effective-medium models of flow through porous media and assumes low Reynolds numbers, fails to describe this capillary-inertial regime, even though it satisfactorily accounts for all previous experiments on foam drainage. However, most of these experiments were performed at the scale of the foam sample, considered steady regimes, and focused on the effects of gravity or of the interfacial rheology; only few experiments investigated the effect of the capillary suction [3,4].

To account for the observed motion of the hydraulic jump, we propose a theoretical description which accounts for our experimental data. The critical distance over which the inertial regime is observed is correctly predicted.

The impact of the present study on the foam drainage dynamics is also discussed. The inertial regime should be observed in macroscopic experiments as a transient feature of drainage under gravity, or in the local flow occuring along the Plateau borders during topological rearrangements of bubbles [5]. Moreover, this regime should be relevant to describe the exact shape of the imbibition front in forced drainage experiments [START_REF] Grassia | EUFOAM2014 Topic: E Preferred type of presentation: Oral[END_REF].

Figures: top

  Figures: top left: the experiment consists in dropping a droplet of the same surfactant solution on a single, horizontal PB. Right: The set of snapshots shows the relaxation of the local perturbation after coalescence of the droplet with the PB. Bottom left: outline of the capillary hydraulic jump geometry.